#*********************************************************************************
#* Copyright (C) 2018-2020 Alexey V. Akimov
#*
#* This file is distributed under the terms of the GNU General Public License
#* as published by the Free Software Foundation, either version 3 of
#* the License, or (at your option) any later version.
#* See the file LICENSE in the root directory of this distribution
#* or <http://www.gnu.org/licenses/>.
#***********************************************************************************
"""
.. module:: models_Tully
:platform: Unix, Windows
:synopsis: This module implements 3 Tully models for testing NA-MD dynamics
.. moduleauthor:: Alexey V. Akimov
"""
import os
import sys
import math
import copy
if sys.platform=="cygwin":
from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
from liblibra_core import *
import util.libutil as comn
import libra_py.units as units
class tmp:
pass
[docs]def Tully1_py(q, params):
"""
Pure Python implementation of the Tully model I = Simple Avoided Crossing (SAC):
H_00 = A*(1.0-exp(-B*x)) x>0,
= A*(exp(B*x)-1.0 ) x<0
H_11 = -H_00
H_01 = C*exp(-D*x^2)
Args:
q ( MATRIX(1,1) ): coordinates of the particle, ndof = 1
params ( dictionary ): model parameters
* **params["A"]** ( double ): [ default: 0.010, units: Ha]
* **params["B"]** ( double ): [ default: 1.600, units: Bohr^-1]
* **params["C"]** ( double ): [ default: 0.005, units: Ha]
* **params["D"]** ( double ): [ default: 1.000, units: Bohr^-2]
Returns:
PyObject: obj, with the members:
* obj.ham_dia ( CMATRIX(2,2) ): diabatic Hamiltonian
* obj.ovlp_dia ( CMATRIX(2,2) ): overlap of the basis (diabatic) states [ identity ]
* obj.d1ham_dia ( list of 1 CMATRIX(2,2) objects ):
derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate
* obj.dc1_dia ( list of 1 CMATRIX(2,2) objects ): derivative coupling in the diabatic basis [ zero ]
"""
critical_params = [ ]
default_params = {"A":0.010, "B":1.600, "C":0.005, "D":1.000 }
comn.check_input(params, default_params, critical_params)
A = params["A"]
B = params["B"]
C = params["C"]
D = params["D"]
Hdia = CMATRIX(2,2)
Sdia = CMATRIX(2,2)
d1ham_dia = CMATRIXList(); d1ham_dia.append( CMATRIX(2,2) )
dc1_dia = CMATRIXList(); dc1_dia.append( CMATRIX(2,2) )
x = q.get(0)
Sdia.set(0,0, 1.0+0.0j); Sdia.set(0,1, 0.0+0.0j);
Sdia.set(1,0, 0.0+0.0j); Sdia.set(1,1, 1.0+0.0j);
V11,dV11 = 0.0, 0.0
if x>0:
V11 = A*(1.0 - math.exp(-B*x))
dV11 = A*B*math.exp(-B*x)
else:
V11 = -A*(1.0 - math.exp(B*x))
dV11 = A*B*math.exp(B*x)
V = C * math.exp(-D*x*x)
dV = -2.0*x*C*D*math.exp(-D*x*x)
Hdia.set(0,0, V11*(1.0+0.0j) ); Hdia.set(0,1, V*(1.0+0.0j));
Hdia.set(1,0, V*(1.0+0.0j)); Hdia.set(1,1, V11*(-1.0+0.0j));
for i in [0]:
# d Hdia / dR_0
d1ham_dia[i].set(0,0, dV11*(1.0+0.0j) ); d1ham_dia[i].set(0,1, dV*(1.0+0.0j));
d1ham_dia[i].set(1,0, dV*(1.0+0.0j)); d1ham_dia[i].set(1,1, dV11*(-1.0+0.0j));
# <dia| d/dR_0| dia >
dc1_dia[i].set(0,0, 0.0+0.0j); dc1_dia[i].set(0,1, 0.0+0.0j);
dc1_dia[i].set(1,0, 0.0+0.0j); dc1_dia[i].set(1,1, 0.0+0.0j);
obj = tmp()
obj.ham_dia = Hdia
obj.ovlp_dia = Sdia
obj.d1ham_dia = d1ham_dia
obj.dc1_dia = dc1_dia
return obj
[docs]def Tully1(q, params):
"""
The implementation that calls the C++ implementation of Tully model I = Simple Avoided Crossing (SAC):
H_00 = A*(1.0-exp(-B*x)) x>0,
= A*(exp(B*x)-1.0 ) x<0
H_11 = -H_00
H_01 = C*exp(-D*x^2)
Args:
q ( MATRIX(1,1) ): coordinates of the particle, ndof = 1
params ( dictionary ): model parameters
* **params["A"]** ( double ): [ default: 0.010, units: Ha]
* **params["B"]** ( double ): [ default: 1.600, units: Bohr^-1]
* **params["C"]** ( double ): [ default: 0.005, units: Ha]
* **params["D"]** ( double ): [ default: 1.000, units: Bohr^-2]
Returns:
PyObject: obj, with the members:
* obj.ham_dia ( CMATRIX(2,2) ): diabatic Hamiltonian
* obj.ovlp_dia ( CMATRIX(2,2) ): overlap of the basis (diabatic) states [ identity ]
* obj.d1ham_dia ( list of 1 CMATRIX(2,2) objects ):
derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate
* obj.dc1_dia ( list of 1 CMATRIX(2,2) objects ): derivative coupling in the diabatic basis [ zero ]
"""
critical_params = [ ]
default_params = {"A":0.010, "B":1.600, "C":0.005, "D":1.000 }
comn.check_input(params, default_params, critical_params)
A = params["A"]
B = params["B"]
C = params["C"]
D = params["D"]
obj = tmp()
obj.ham_dia = CMATRIX(2,2)
obj.ovlp_dia = CMATRIX(2,2)
obj.d1ham_dia = CMATRIXList(); obj.d1ham_dia.append( CMATRIX(2,2) )
obj.dc1_dia = CMATRIXList(); obj.dc1_dia.append( CMATRIX(2,2) )
# Convert MATRIX to doubleList()
qq = doubleList(); qq.append(q.get(0))
prms = doubleList() # will use the default values
prms.append(A); prms.append(B); prms.append(C); prms.append(D);
model_SAC(obj.ham_dia, obj.ovlp_dia, obj.d1ham_dia, obj.dc1_dia, qq, prms);
return obj
[docs]def Tully2(q, params):
"""
The implementation that calls the C++ implementation of Tully model II = Double Avoided Crossing (DAC):
H_00 = 0.0
H_11 = E - A*exp(-B*x^2)
H_01 = C*exp(-D*x^2)
Args:
q ( MATRIX(1,1) ): coordinates of the particle, ndof = 1
params ( dictionary ): model parameters
* **params["A"]** ( double ): [ default: 0.100, units: Ha]
* **params["B"]** ( double ): [ default: 0.028, units: Bohr^-2]
* **params["C"]** ( double ): [ default: 0.015, units: Ha]
* **params["D"]** ( double ): [ default: 0.060, units: Bohr^-2]
* **params["E"]** ( double ): [ default: 0.050, units: Ha]
Returns:
PyObject: obj, with the members:
* obj.ham_dia ( CMATRIX(2,2) ): diabatic Hamiltonian
* obj.ovlp_dia ( CMATRIX(2,2) ): overlap of the basis (diabatic) states [ identity ]
* obj.d1ham_dia ( list of 1 CMATRIX(2,2) objects ):
derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate
* obj.dc1_dia ( list of 1 CMATRIX(2,2) objects ): derivative coupling in the diabatic basis [ zero ]
"""
critical_params = [ ]
default_params = {"A":0.10, "B":0.28, "C":0.015, "D":0.060, "E":0.050 }
comn.check_input(params, default_params, critical_params)
A = params["A"]
B = params["B"]
C = params["C"]
D = params["D"]
E = params["E"]
obj = tmp()
obj.ham_dia = CMATRIX(2,2)
obj.ovlp_dia = CMATRIX(2,2)
obj.d1ham_dia = CMATRIXList(); obj.d1ham_dia.append( CMATRIX(2,2) )
obj.dc1_dia = CMATRIXList(); obj.dc1_dia.append( CMATRIX(2,2) )
# Convert MATRIX to doubleList()
qq = doubleList(); qq.append(q.get(0))
prms = doubleList() # will use the default values
prms.append(A); prms.append(B); prms.append(C); prms.append(D); prms.append(E);
model_DAC(obj.ham_dia, obj.ovlp_dia, obj.d1ham_dia, obj.dc1_dia, qq, prms);
return obj
[docs]def Tully3(q, params):
"""
The implementation that calls the C++ implementation of Tully model III =
Extended Coupling With Reflection (ECWR):
H_00 = A
H_11 = -H_00
H_01 = B*exp(C*x); x <= 0
B*(2.0 - exp(-C*x)); x > 0
Args:
q ( MATRIX(1,1) ): coordinates of the particle, ndof = 1
params ( dictionary ): model parameters
* **params["A"]** ( double ): [ default: 0.0006, units: Ha ]
* **params["B"]** ( double ): [ default: 0.1000, units: Ha ]
* **params["C"]** ( double ): [ default: 0.9000, units: Bohr^-1 ]
Returns:
PyObject: obj, with the members:
* obj.ham_dia ( CMATRIX(2,2) ): diabatic Hamiltonian
* obj.ovlp_dia ( CMATRIX(2,2) ): overlap of the basis (diabatic) states [ identity ]
* obj.d1ham_dia ( list of 1 CMATRIX(2,2) objects ):
derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate
* obj.dc1_dia ( list of 1 CMATRIX(2,2) objects ): derivative coupling in the diabatic basis [ zero ]
"""
critical_params = [ ]
default_params = {"A":0.0006, "B":0.1000, "C":0.9000 }
comn.check_input(params, default_params, critical_params)
A = params["A"]
B = params["B"]
C = params["C"]
obj = tmp()
obj.ham_dia = CMATRIX(2,2)
obj.ovlp_dia = CMATRIX(2,2)
obj.d1ham_dia = CMATRIXList(); obj.d1ham_dia.append( CMATRIX(2,2) )
obj.dc1_dia = CMATRIXList(); obj.dc1_dia.append( CMATRIX(2,2) )
# Convert MATRIX to doubleList()
qq = doubleList(); qq.append(q.get(0))
prms = doubleList() # will use the default values
prms.append(A); prms.append(B); prms.append(C)
model_ECWR(obj.ham_dia, obj.ovlp_dia, obj.d1ham_dia, obj.dc1_dia, qq, prms);
return obj
[docs]def chain_potential(q, params, full_id):
"""
A 1D linear chain potential.
This is the Hamiltonian to mimic the one used by Tully and Parandekar to study
thermal equilibrium in quantum-classical systems
Although the indended interpretation is that the first DOFs is quantum, this potential
is totally generic, so all DOFs could be treated as quantum
Args:
q ( MATRIX(ndof,1) ): coordinates of the nuclear DOFs
params ( dictionary ): model parameters
* **params["A"]** ( double ): [ units: None ]
* **params["a"]** ( double ): [ units: bohr-1 ]
* **params["V0]** ( double ): [ units: Ha ]
* **params["E_n"]** ( list of doubles ): [ default: [0.0, 0.001, 0.001, 0.001], units: Ha]
* **params["x_n"]** ( list of doubles ): [ default: [0.0, 1.0, 1.0, 1.0], units: Bohr]
* **params["k_n"]** ( list of doubles ): [ default: [0.001, 0.001, 0.001, 0.001], units: Ha/Bohr^2]
* **params["V"]** ( list of lists of double ): [ default: [[0.001]*4]*4, units: Ha]
Returns:
PyObject: obj, with the members:
* obj.ham_dia ( CMATRIX(n,n) ): diabatic Hamiltonian
* obj.ovlp_dia ( CMATRIX(n,n) ): overlap of the basis (diabatic) states [ identity ]
* obj.d1ham_dia ( list of 1 CMATRIX(n, n) objects ):
derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate
* obj.dc1_dia ( list of 1 CMATRIX(n, n) objects ): derivative coupling in the diabatic basis [ zero ]
"""
Id = Cpp2Py(full_id)
indx = Id[-1]
X = q.col(indx) # coordinates of all particles for this trajectory
ndof = q.num_of_rows # the total number of DOFs, the first one is the quantum DOF
critical_params = ["E_n", "x_n", "k_n" ]
default_params = { "V": [ [0.001, 0.001, 0.001, 0.001], [0.001, 0.001, 0.001, 0.001],
[0.001, 0.001, 0.001, 0.001], [0.001, 0.001, 0.001, 0.001] ],
"a":0.25, "U0":0.06665
}
comn.check_input(params, default_params, critical_params)
# Parameters describing the diabatic PES for quantum DOF
E_n = params["E_n"] # energies of the minima of the diabatic PESs for quantum particle
x_n = params["x_n"] # positions of diabatic PESs
k_n = params["k_n"] # curvatures of diabatic PESs
V = params["V"] # diabatic couplings
nstates = len(E_n) # the number of electronic states
# Parameters describing classical particles and their coupling to the quantum particle
a = params["a"]
U0 = params["U0"]
a2 = U0*a**2
a3 = -U0*a**3
a4 = 0.58*U0*a**4
Hdia = CMATRIX(nstates,nstates)
Sdia = CMATRIX(nstates,nstates)
Sdia.identity()
basis_transform = CMATRIX(nstates,nstates)
basis_transform.identity()
d1ham_dia = CMATRIXList();
dc1_dia = CMATRIXList();
for k in range(ndof):
d1ham_dia.append( CMATRIX(nstates,nstates) )
dc1_dia.append( CMATRIX(nstates,nstates) )
#========= Classical (bath) contributions =========
H_bath = 0.0
dH_bath = []
for k in range(ndof-1):
x1 = X.get(k)
x2 = X.get(k+1)
d = x1 - x2
d2 = d*d
d3 = d2*d
d4 = d3*d
H_bath += (a2*d2 + a3*d3 + a4*d4)
dH_bath.append( 2.0*a2*d + 3.0*a3*d2 + 4.0*a4*d3 )
for i in range(nstates):
Hdia.add(i,i, H_bath*(1.0+0.0j))
k = 0
d1ham_dia[k].add(i, i, ( dH_bath[k] ) * (1.0+0.0j) )
for k in range(1, ndof-1):
d1ham_dia[k].add(i, i, ( dH_bath[k] - dH_bath[k-1] ) * (1.0+0.0j) )
k = ndof - 1
d1ham_dia[k].add(i, i, ( -dH_bath[k-1] ) * (1.0+0.0j) )
"""
for k in range(1, ndof):
d = (X.get(0) - X.get(k) )
d2 = d*d
H_bath = 0.5 * U0*d2
dH_bath = U0 * d
for i in range(nstates):
Hdia.add(i,i, H_bath*(1.0+0.0j) )
d1ham_dia[k].add(i, i, ( -dH_bath ) * (1.0+0.0j) )
d1ham_dia[0].add(i, i, ( dH_bath ) * (1.0+0.0j) )
"""
#========== Quantum system ==================
x = X.get(0)
for i in range(nstates):
# Energies
Hdia.add(i,i, (E_n[i] + 0.5*k_n[i]*(x - x_n[i])**2) * (1.0+0.0j) )
for j in range(nstates):
if i!=j:
Hdia.add(i,j, V[i][j] * (1.0+0.0j) )
# Derivatives
k = 0 # quantum DOF
d1ham_dia[k].add(i,i, (k_n[i] * (x - x_n[i]))*(1.0+0.0j) )
obj = tmp()
obj.ham_dia = Hdia
obj.ovlp_dia = Sdia
obj.d1ham_dia = d1ham_dia
obj.dc1_dia = dc1_dia
#obj.basis_transform = basis_transform
return obj