#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it 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.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
import os
import sys
import datetime
import uuid
import math as mt
import numpy as np
from numpy import linalg as LA
from scipy.integrate import ode
import CDTK.Tools.Conversion as conv
import CDTK.Tools.Mathematics as tmath
IDLEN = 8
TINY = 1.0e-8
[docs]
class Trajectory_SH(object):
"""
Surface Hopping trajectory
self.NDof -- int, Number of degrees of freedom
self.NStates -- int, Number of electronic states
self.R -- np.array, positions
self.V -- np.array, velocities
self.M -- np.array float, masses
self.dt -- float, time-step for nuclear displacements
self.rseed -- int, default None, random seed to be set before 1st step
self.C -- np.array complex, quantum coefficients
self.NAME -- Identification string of this trajectory. It is uniquely
generated at the initialization
self.f_EGrad -- function, provides energy and gradient on PES of state S
e,g = f(R,V,S,t)
self.f_NAC -- function, provides NAC vector between states sI,sJ
dIJ = f(R,sI,sJ,t)
self.f_W -- function, provides potential matrix for quantum evolution
Returns a diagonal matrix with adiabatic energies if purely
adiabatic representation
self.f_D -- function, provides the matrix of kinetic couplings
-i dIJ * V
This matrix would be zero in a purely diabatic representation
self.f_DD -- function, provides the matrix of kinetic couplings
-i dIJ * V
This matrix provides couplings to an external system,
for example an electric field, and therefore the
transitions induced by do not conserve total energy. In
practice, no velocity adjustments are performed when a
transition takes place.
self.ES -- electronic state for large hop electron dynamics
Notes on the supplied external functions:
f_EGrad : In the usual Tully SH prescription the gradients of adibatic states
are used, for example from an electronic structure program. This is
the usual prescription. The gradient is used exclusively for the
evolution of the classical DOF.
f_NAC : Non-adiabatic coupling vector between a pair of electronic states.
It is used in the adjustment of the velocities after a quantum
transition.
f_NAC is _not_ used in the quantum evolution
f_W : The potential part of the electronic Hamiltonian (depends only on
nuclear positions). In the usual SH prescription this is a diagonal
matrix with the adiabatic PES. The matrix provided will be used as it
is (e.g. if a diabatic matrix is provided) unless differently
specified.
f_D : The kinetic part of the electronic Hamiltonian (depends on positions
and velocities of the classical DOF). If a diabatic representation is
provided this function may remain undefined.
"""
def __init__(self,NDof=1,NStates=1,**opts):
self.NDof = NDof #### Degrees of freedom
self.NStates = NStates #### No. of electronic states
self.R = np.zeros(self.NDof,float) #### Position
self.V = np.zeros(self.NDof,float) #### Velocity
self.M = np.ones(self.NDof,float) #### Mass
self.dt = opts.get('dt',10.0) #### Time-step (au)
self.dt_max = opts.get('dt_max',10.0)
self.noreturn = opts.get('noreturn',None)
self.rseed = opts.get('rseed',None) #### Random seed
self.C = np.zeros(self.NStates,complex) #### Array of coefficients
self.C[0] = 1.0 + 0.0j
self.t = 0.0 #### Current time
self.S = 0 #### Electronic state
self.f_EGrad = opts.get('f_EGrad',None)
self.f_W = opts.get('f_W',None)
self.f_D = opts.get('f_D',None)
self.f_DD = opts.get('f_DD',None)
self.f_NAC = opts.get('f_NAC',None)
self.f_overlap = opts.get('f_overlap',None)
self._ID = uuid.uuid4().hex[0:IDLEN] #### unique ID
self._CTIME = datetime.datetime.today() #### creation date/time
self._BNAME = 'trj_%4i%02i%02i_' % (self._CTIME.year,self._CTIME.month,self._CTIME.day)
self.NAME = self._BNAME + self._ID #### Used to name output files, etc.
self.DIR = self.NAME #### Name directory
self._stepnum = 0 #### step counter
self.initiator = 0 #### Iniates the dynamics
self.rescaling = 'nac' #### (nac|grd) method to adjust velocities
self.is_sh = True #### whether SH is performed
self.first_step = 0.01 #### Initial step of quantum integrator
self.quantum_step_integrator = 'ipicture' #### (ipicture|unitary)
self.counter1 = 0 #### For adaptive step size
self.p_thresh = opts.get('p_thresh', 1.0) #### Parameter for the adaptive step size
self.field_implicit = opts.get('field_explicit','None') #### For the calculation of probability "b"
self.no_hops = opts.get('no_hops',False) #### whether hops between surfaces are allowed at all. If True, trajectories can be compared to adiabatic quantum dynamics
# whether overall translation should be removed
self.rcom = False
# whether overall rotation should be removed
self.rrot = False
# treatement of trivialCrossings
self.trivialCrossing = opts.get('trivialCrossing',None)
self.dE_thresh = None
self.ES = None # electronic state object for Monte-Carlo electron dynamics
self.constraints = False
self.getPartial = None # partial charge calculation
self.ekin = 0
[docs]
def set_constraints(self, f_constraints):
self.constraints = True
self.getConstraintsCorrection = f_constraints
[docs]
def restart(self):
"""
Function to restart trajectory. Reads in all important values from log-files.
"""
if type(1) == type(self.rseed):
np.random.seed(self.rseed)
self.logR = [] # List of positions along trajectory
self.logV = [] # List of velocities along trajectory
self.logC = [] # List of coefficient vectors
self.logS = [] # List of coefficient vectors
self.logt = [] # List of times along the simulation
self.logKE = [] # List of Total Energy along trajectory
self.logPE = [] # List of Total Energy along trajectory
self.logH = [] # diagonals of the Hamiltonian (ad PES)
self.logWarning = [] # List of warning messages
self.logSwitch = []
if self.getPartial is not None:
self.logPartial = [] # Log partial charges
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
self.logCrossings = []
# get R, V, t, S, stepnum, NAC
with open(self.DIR + '/R.log','r') as posfile:
for line in posfile:
self._stepnum += 1
vals = line.split()
self.logt.append( float(vals[0]) )
self.logR.append(np.array( list(map(float, vals[1:])) ))
with open(self.DIR + '/V.log','r') as velfile:
for line in velfile:
vals = line.split()
self.logV.append(np.array( list(map(float, vals[1:])) ))
with open(self.DIR + '/C.log','r') as coefffile:
for line in coefffile:
vals = line.split()
cc = np.zeros((len(vals)-1)//2, dtype=complex)
for i in range(1,len(vals),2):
cc[(i-1)//2] = complex(float(vals[i]), float(vals[i+1]))
self.logC.append(cc)
with open(self.DIR + '/S.log','r') as statefile:
for line in statefile:
vals = line.split()
self.logS.append(int(vals[1]))
with open(self.DIR + '/NAC.log','r') as nacfile:
hh = np.zeros((self.NStates, self.NStates), dtype=complex)
counter = 0
for line in nacfile:
if counter % (self.NStates*self.NStates) == 0:
hh = np.zeros((self.NStates, self.NStates), dtype=complex)
vals = line.split()
hh[int(vals[0]), int(vals[1])] = complex(float(vals[3]), float(vals[4]))
counter += 1
if counter % (self.NStates*self.NStates) == 0:
self.logH.append(hh)
with open(self.DIR + '/E.log','r') as efile:
for line in efile:
vals = line.split()
self.logKE.append(float(vals[1]))
self.logPE.append(float(vals[2]))
self._EREF = np.diag(np.ones(self.NStates,complex)) * self.logPE[0]
with open(self.DIR + '/Switch.log','r') as switchfile:
for line in switchfile:
self.logSwitch.append(line)
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
with open(self.DIR + '/trivialCrossings.log','r') as f:
i = 0
thisStepCrossing = []
for line in f:
if line.startswith('#'):
continue
line = line.rstrip().split()
tc = float(line[0])
a = int(line[1])
b = int(line[2])
while self.logt[i] < tc:
self.logCrossings.append(thisStepCrossing)
i = i + 1
thisStepCrossing = []
thisStepCrossing.append([a,b])
self.logCrossings.append(thisStepCrossing)
i = i + 1
thisStepCrossing = []
while i < len(self.logt):
self.logCrossings.append(thisStepCrossing)
i = i + 1
if self.getPartial is not None:
with open(self.DIR + '/partial.log','r') as partialfile:
for line in partialfile:
vals = line.split()
self.logPartial.append(np.array( list(map(float, vals[1:])) ))
self.t = self.logt[-1]
self.R = self.logR[-1]
self.V = self.logV[-1]
self.C = self.logC[-1]
self.S = self.logS[-1]
self.ekin = self.kinetic_energy()
if self.getPartial is not None:
self.partialCharges = self.logPartial[-1]
self._pot,self._grd = self.f_EGrad(self.R,self.V,self.t,self.S) # get initial POT, GRAD
# Get Hamiltonian and remove large diagonal
self._E0 = self.f_W( self.R,self.V,self.t ) - self._EREF
self._D0 = self.f_D( self.R,self.V,self.t )
self._H0 = self._E0 + self._D0
if self.f_DD:
self._D0 = self._D0 - self.f_DD( self.R,self.V,self.t )
self._H0 = self._E0 + self._D0
self._H0ad = self._H0.copy()
[docs]
def integrate(self, tfinal = None, **opts):
"""
Integrate trajectory until tfinal
tfinal -- final time to be reached by the integrator in steps of dt.
**opts:
dt -- delta time used in this segment.
It is not a good idea to change this when continuing from a
previous integration
"""
if self.initiator == 0: #### Initial step of the calculation
self.dt = opts.get( 'dt', self.dt )
if tfinal is None:
tfinal = self.t + self.dt
if self.f_EGrad is None:
raise ValueError('self.f_EGrad is not defined, cannot integrate')
if self._stepnum == 0:
if self.rseed is not None:
np.random.seed(self.rseed)
# Initialization
self._pot,self._grd = self.f_EGrad(self.R,self.V,self.t,self.S) #### get initial POT, GRAD
self.ekin = self.kinetic_energy()
# Diagonal matrix with large diagonal
self._EREF = np.diag(np.ones(self.NStates,complex)) * self._pot
# Get Hamiltonian and remove large diagonal
W = self.f_W( self.R,self.V,self.t )
if(W.shape[0] != self.NStates):
raise ValueError(f"The number of states requested for the propagation is {self.NStates}.\n" +
f"The number of states returned from the quantum chemistry enging is {W.shape[0]}.\n" +
f"Please check your input! [You may want to add nstates={self.NStates} to the xmolecule section]")
self._E0 = self.f_W( self.R,self.V,self.t ) - self._EREF
self._D0 = self.f_D( self.R,self.V,self.t )
if self.f_DD:
self._D0 = self._D0 - self.f_DD( self.R,self.V,self.t )
self._H0 = self._E0 + self._D0
self._H0ad = self._H0.copy()
# Start logging
if self.initiator == 0: #### Initiates the logging
self.init_log()
self.initiator += 1
while self.t < tfinal:
self._step()
if self.counter1 == 0: #### Steps
self._stepnum += 1
self.t = self.t + self.dt
self.log()
# Check for energy conservation
if self.dE_thresh is not None:
if abs((self.logKE[-1] + self.logPE[-1] - self.logKE[0] - self.logPE[0]) / \
( self.logKE[0] + self.logPE[0] ) ) > self.dE_thresh:
raise RuntimeError('Energy is not conserved. I will write the logfiles and exit')
self.dt = self.dt1
if self.counter1 == 1:
tfinal = self.t + self.dt1
self.dt = self.dt1
break
def _step(self,**opts):
"""
Do a Tully's surface hopping integration step
"""
self._classical_step()
if self.noreturn != None and LA.norm(self.R[0:3] - self.R[3:6]) > float(self.noreturn):
print("exceeding")
exit(-1)
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
self._trivial_crossing()
self._quantum_step()
#if len(self.crossings)==0:
self._surface_hopping()
if self.getPartial is not None:
self.partialCharges = self.getPartial(state=self.S)
def _check_trivial_crossing_step(self,x0,x1,**opts):
"""
checks for trivial crossings between geometry step x0 and x1
by maximizing the trace of the overlap matrix square
see JCP 137, 014512 (2012)
"""
# state overlap matrix
S=self.f_overlap(x0, x1, S=self.S)
# state overlap matrix squared
S2=S*S
converged=False
#print "State overlap^2 Matrix:"
#for i in range(0,self.NStates):
# print S2[i,0:self.NStates]
# rearange the columns in S2 such that Trace{S2} is maximal
pivot = np.arange(self.NStates)
while not converged:
converged = True
for i in range(0,len(S2)):
for j in range(0,len(S2)):
#print(i,j,S2[i,pivot[j]] + S2[j,pivot[i]], S2[i,pivot[i]] + S2[j,pivot[j]])
if( S2[i,pivot[i]] + S2[j,pivot[j]] < S2[i,pivot[j]] + S2[j,pivot[i]] ):
pivot[i],pivot[j]=pivot[j].copy(),pivot[i].copy()
converged=False
# find the smallest diagonal
min_diag = min([ S2[i,pivot[i]] for i in range(len(pivot)) ])
return pivot, S2[:, :]
def _trivial_crossing(self,**opts):
"""
Performs trivial crossing hops according to overlap criterion
see JCP 137, 014512 (2012)
"""
x0=self.x0
x1=self.R
self.crossings=[]
threshold=0.6
smallest_step = 0.5
N=1
d = 0.
step = 1.
pivot = np.arange(self.NStates)
while d < 1.:
if d + step > 1:
step = 1. - d
# print("taking a step from %f to %f " % (d, d + step))
sx0 = x0 + d * (x1-x0)
d = d + step
sx1 = x0 + d * (x1-x0)
step_pivot, S2 = self._check_trivial_crossing_step(sx0,sx1)
min_diag = min([ S2[i,step_pivot[i]] for i in range(len(step_pivot)) ])
idx = np.argmin(np.array([ S2[i,step_pivot[i]] for i in range(len(step_pivot)) ]))
# print("# min_diagonal_overlap=%f state=%d" % (min_diag, idx))
sys.stdout.flush()
if min_diag < threshold and step > smallest_step:
# go back to the previous step
d = d - step
# decrease step
step = step/2.
print("# decreasing step size to %f" % step)
continue
if min_diag < threshold:
print("loss of identity for state %d" % idx)
if min_diag > threshold and step < 0.5:
print("# increasing step size")
# increase step
step = step * 2
print("# increasing step size to %f" % step)
pivot[:] = pivot[step_pivot].copy()
sys.stdout.flush()
# pivot maps state at x0 onto states at x1
# the reverse mapping is pivot_reverse: states at x1 onto states at x0
pivot_reverse=np.arange(len(pivot))
for i in range(0,len(pivot)):
pivot_reverse[pivot[i]] = i
# Give some output and append states that are involved in self.crossings
for i in range(0,len(pivot)):
if i != pivot[i]:
print("# Detected a trivial crossing: state %d -> %d " % (i,pivot[i]))
self.crossings.append((i,pivot[i]))
# do a jump in case of our active state is changed
hopped=False
if pivot[self.S]!=self.S:
print("# state crossing from {:d} to {:d}".format(self.S,pivot[self.S]))
self.S=pivot[self.S]
hopped=True
# readjust all quantities to the new state indices at x1
H0old=self._H0.copy()
D0old=self._D0.copy()
E0old=self._E0.copy()
Cold=self.C.copy()
for i in range(self._H0.shape[0]):
self.C[pivot[i]] = Cold[i]
for j in range(self._H0.shape[1]):
self._H0[pivot[i],pivot[j]] = H0old[i,j]
self._D0[pivot[i],pivot[j]] = D0old[i,j]
self._E0[pivot[i],pivot[j]] = E0old[i,j]
# return whether we have hopped
return hopped
def _classical_step(self,**opts):
"""
Advance classical coordinates with a velocity Verlet step t0 -> t1.
"""
m = self.M
self.x0 = self.R
self.v0 = self.V
self.g0 = self._grd
self._potdum = self._pot
dt = self.dt
# Classical step (Vel. Verlet)
v1 = self.v0 - ((0.5*dt)*(self.g0/m)) #### v(t) -> v(t+dt/2)
if self.constraints:
x1 = self.x0 + v1*dt
v1 += self.getConstraintsCorrection(self.x0, x1, dt, case="pos")
v1 = self._remove_com_motion(self.x0, v1, m, self.rcom, self.rrot)
# include ionization dynamics steps here
if self.ES and self.ES.state != 'doNotRun':
hopped = self.ES.step(dt)
x1 = self.x0 + v1*dt #### x(t) -> x(t+dt)
e1,g1 = self.f_EGrad( x1,v1,self.t+dt,self.S )
v1 = v1 - (0.5*dt*g1)/m #### v(t+dt/2) -> v(t+dt)
# the kinetic energy is calculated before applying velocity constraints.
# this makes a smoother transition based on gromacs data
self.V = v1
self.ekin = self.kinetic_energy()
if self.constraints:
v1 += self.getConstraintsCorrection(x1, v1, dt, case="vel")
v1 = self._remove_com_motion(x1, v1, m, self.rcom, self.rrot)
self._grd = g1 #### save gradient for next step
self._pot = e1
self.R = x1
self.V = v1
def _remove_com_motion(self, r, v, m, rcom, rrot):
"""
TODO: The inertia tensor is currently regularized with 10^-10*Id-Matrix!
This function removes the center of mass motion by setting the center of mass velocity to zero after every
velocity step in velocity verlet.
Input:
r -- 3N array of positions in bohr
v -- 3N array of velocities in au
m -- 3N array of masses in au
rcom -- Whether center of mass motion should be removed
rrot -- Whether center of mass rotation should be removed
Output:
v -- 3N array of new velocities without center of mass motion/rotation in au
"""
if rcom == False:
return v
# n = number of atoms
n = len(m) // 3
# center of mass velocity
v_com = np.sum(v.reshape(n,3) * m.reshape(n,3), axis=0) / np.sum(m) * 3.0
v_new = np.zeros(n*3)
# remove center of mass velocity for each atom
for i in range(n):
v_new[i*3:(i+1)*3] = v[i*3:(i+1)*3] - v_com
if rrot == False:
return v_new
### calculate overall rotation with respect to center of mass and substract
# L is overall angular momentum
L = np.zeros(3)
# Theta is inertia tensoer
Theta = np.zeros((3,3))
# calculate center of mass
com = np.sum(r.reshape(n,3) * m.reshape(n,3), axis=0) / np.sum(m) * 3.0
for i in range(n):
# calculate particle position relative to center of mass
r_com = r[i*3:(i+1)*3] - com
# calculate overall angular momentum
L += (m[i*3] * np.cross(r_com, v_new[i*3:(i+1)*3]))
# calculate inertia tensor
Theta += m[i*3] * (np.dot(r_com, r_com) * np.identity(3) - np.outer(r_com, r_com))
# regularize Theta!
Theta += 1e-10 * np.identity(3)
for i in range(n):
r_com = r[i*3:(i+1)*3] - com
# v_p = v_p - Th^-1 * L (X) r_p
try:
v_new[i*3:(i+1)*3] -= np.cross(np.dot(np.linalg.inv(Theta), L), r_com)
except np.linalg.linalg.LinAlgError:
return v_new
return v_new
def _quantum_step(self,**opts):
# self._quantum_step_def()
if self.quantum_step_integrator == 'ipicture':
self._quantum_step_IP()
elif self.quantum_step_integrator == 'unitary':
self._quantum_step_unitary()
else:
self._quantum_step_def()
def _quantum_step_def(self,**opts):
"""
Advance expansion coefficients t0 -> t1
"""
self._H0dum = self._H0 #### To reuse the Hamitonian from previous time step after the adaptive step size
self._E0dum = self._E0 #### To reuse the Potential energy (diagonal terms) previous time step after the adaptive step size
self._D0dum = self._D0 #### To reuse the coupling (off diagonal) terms previous time step after the adaptive step size
self.Cdum = self.C #### To reuse the co-efficients after the adaptive step size
dt = self.dt
# Quantum step (Schroedinger) t0 -> t1
# Tully 1990, Eq. (8).
self._E1 = self.f_W( self.R,self.V,self.t+dt ) - self._EREF #### H(t1)
self._D1 = self.f_D( self.R,self.V,self.t+dt )
if self.f_DD:
self._D1 = self._D1 - self.f_DD( self.R,self.V,self.t+dt )
self._H1 = self._E1 + self._D1
self._Hh = 0.5*( self._H0 + self._H1 ) #### ( H(t0) + H(t1) )/2
def f(t,c):
r = (t-self.t)/dt
r = 0.5
H = (1.0-r)*self._H0 + r*self._H1
return -1.0j * np.dot(H,c)
self.q = ode(f)
self.q.set_integrator('zvode',
first_step=self.first_step)#, method='bdf')
self.q.set_initial_value(self.C, self.t)
while self.q.successful() and self.q.t < (self.t + self.dt):
self.q.integrate(self.q.t + self.dt)
self.C = self.q.y
self._H0 = self._H1 #### Store electr. Hamiltonian for reuse
self._H0ad = self._H0.copy()
self._E0 = self._E1 #### Store electr. Hamiltonian for reuse (diagn)
self._D0 = self._D1 #### Store electr. Hamiltonian for reuse (coupl)
def _quantum_step_IP(self,**opts):
"""
Advance expansion coefficients t0 -> t1
This function uses an interaction picture
"""
self._H0dum = self._H0 #### To reuse the Hamitonian from previous time step after the adaptive step size
self._E0dum = self._E0 #### To reuse the Potential energy (diagonal terms) previous time step after the adaptive step size
self._D0dum = self._D0 #### To reuse the coupling (off diagonal) terms previous time step after the adaptive step size
self.Cdum = self.C #### To reuse the co-efficients after the adaptive step size
dt = self.dt
# Quantum step (Schroedinger) t0 -> t1
# Tully 1990, Eq. (8).
self._E1 = self.f_W( self.R,self.V,self.t+dt ) - self._EREF #### H(t1)
self._D1 = self.f_D( self.R,self.V,self.t+dt ) #### It is negative beacause of the Tully 1980
if self.f_DD:
self._D1 = self._D1 - self.f_DD( self.R,self.V,self.t+dt ) #### It is positive after the Helman Feynman theorem
self._H1 = self._E1 + self._D1
# Schroedinger picture V
E1mE0 = self._E1 - self._E0
D1mD0 = self._D1 - self._D0
self._Hh = 0.5*( self._H0 + self._H1 ) #### ( H(t0) + H(t1) )/2
def f(t,c):
r = (t-self.t)/dt
VS = E1mE0*r + D1mD0*r + self._D0
VI = tmath.schroedinger_to_interaction(VS,self._E0,t-self.t)
return -1.0j*np.dot(VI,c)
self.q = ode(f)
self.q.set_integrator('zvode',
first_step=self.first_step) #, method='bdf')
self.q.set_initial_value(self.C, self.t)
while self.q.successful() and self.q.t < (self.t + self.dt):
self.q.integrate(self.q.t + (self.dt))
self.C = self.q.y
# Change back to Schroedinger picture
self.C = np.exp( -1j*dt*np.diag(self._H0) ) * self.C
self._H0 = self._H1 #### Store electr. Hamiltonian for reuse (total)
self._H0ad = self._H0.copy()
self._E0 = self._E1 #### Store electr. Hamiltonian for reuse (diagn)
self._D0 = self._D1 #### Store electr. Hamiltonian for reuse (coupl)
def _quantum_step_unitary(self,**opts):
"""
Advance expansion coefficients t0 -> t1
This function uses a unitary propagator based on the interpolated
Hamiltonian at t=(t1-t0)/2
TODO: There seems to be some bug here!!
"""
self._H0dum = self._H0 #### To reuse the Hamitonian from previous time step after the adaptive step size
self._E0dum = self._E0 #### To reuse the Potential energy (diagonal terms) previous time step after the adaptive step size
self._D0dum = self._D0 #### To reuse the coupling (off diagonal) terms previous time step after the adaptive step size
self.Cdum = self.C #### To reuse the co-efficients after the adaptive step size
# Quantum step (Schroedinger) t0 -> t1
# Tully 1990, Eq. (8).
dt = self.dt
self._E1 = self.f_W( self.R,self.V,self.t+dt ) - self._EREF #### H(t1)
self._D1 = self.f_D( self.R,self.V,self.t+dt )
if self.f_DD:
self._D1 = self._D1 - self.f_DD( self.R,self.V,self.t+dt )
self._H1 = self._E1 + self._D1
self._Hh = 0.5*( self._H0 + self._H1 ) #### ( H(t0) + H(t1) )/2
ener,trafo = np.linalg.eigh(self._Hh)
diag = np.diag( np.exp(-1.0j*ener*dt) )
trafoH = np.conjugate(trafo.T)
propagator = np.dot(trafo, np.dot(diag, trafoH))
self.C = np.dot(propagator,self.C)
self._H0 = self._H1 #### Store electr. Hamiltonian for reuse (total)
self._H0ad = self._H0.copy()
self._E0 = self._E1 #### Store electr. Hamiltonian for reuse (diagn)
self._D0 = self._D1 #### Store electr. Hamiltonian for reuse (coupl)
def _surface_hopping(self,**opts):
"""
Compute hopping probabilities out of state self.S, hopp if necessary
Tully 1990, Eq. (14).
"""
# save old state
s0 = self.S
# Here the probability for hopping is calculated
dt = self.dt
b = 2.0 * ( self.f_D( self.R,self.V,self.t+dt )[self.S,:]) * ( np.conj(self.C[self.S]) * self.C)
if self.f_DD:
if self.field_explicit:
b = b - 2.0 * ( self.f_DD( self.R,self.V,self.t+dt )[self.S,:]) * ( np.conj(self.C[self.S]) * self.C)
else:
b = b + 2.0 * ( self.f_DD( self.R,self.V,self.t+dt )[self.S,:]) * ( np.conj(self.C[self.S]) * self.C)
b[self.S] = 0.0 #### where "b" gives the coupling between the current state self.S and all other states
a = np.conj(self.C[self.S]) * self.C[self.S] #### population of state self.S
# RW: th sign in p has to be -; checked with Tully models!! TODO: write good unit tests and update docu
p = -self.dt*b.imag/(a.real) #### "p" gives the probability of hopping from the state self.S to other states ---> Tully 1990, Eq. (19)
p[p<0.0] = 0.0 #### Set to zero all negative probabilities
# Normalizing the probability to 1
psum = p.sum()
if psum > 1.0:
p = p/( psum + TINY )
self._issue_warning('Switch prob. larger than 1. / rescaled / smaller DT might be necessary')
# Adaptive step size calculation
self.counter1 = 0 #### This counter indicates the change in step size
p_trans = max(p) #### maximum probibility to hop from state self.S to any other state
if p_trans == 0:
self.dt1 = dt
else:
dt1 = (float(self.p_thresh)/p_trans)*dt
self.dt1 = dt1
if dt1 < dt:
self.dt1 = dt1
self.counter1 = 1
self.R = self.x0
self.V = self.v0
self.C = self.Cdum
self._H0 = self._H0dum
self._E0 = self._E0dum
self._D0 = self._D0dum
self._grd = self.g0
self._pot = self._potdum
return
if self.dt_max > dt1 > dt:
self.dt1 = dt1
if dt1 > self.dt_max:
self.dt1 = self.dt_max
# Finding the final state to hop using a random number
r = np.random.random() #### Here we check for total probability being smaller than 1 and warning issued for probability larger than 1.
self._finalS = _switch_state(p,r) #### self._finalS contains the final state number
# Attempt a switch from state self.S to state self._finalS
hopping_no_field = False #### Indicates the hopping due to no electric field
hopping_with_field = False #### Indicates the hopping due to electric field
steplog = '-' #### (-, no switch|* switch attempted|+ switch succeed)
if self._finalS > -1: #### state switch will take place if enough KE available
steplog = '*'
es,gs = self.f_EGrad( self.R, self.V, self.t+dt, self._finalS ) #### POT, GRAD for final state
if self.f_DD:
b_mol = (self.f_D( self.R,self.V,self.t+dt )[self.S,self._finalS ]) * (self.C[self._finalS ] * np.conj(self.C[self.S]))
b_ext = (self.f_DD( self.R,self.V,self.t+dt )[self.S,self._finalS ]) * (self.C[self._finalS ] * np.conj(self.C[self.S]))
x = abs(b_ext.imag)/abs(b_mol.imag)
pp = 0.5*(1+mt.tanh(x-(1/x)))
r1 = np.random.normal(0.5,0.2,1)
dec = pp+r1
if dec <= 1:
self._grd_finalS = gs
dKE = self._pot - es #### Kinetic E increase of _nuclei_ (can be negative if transition to upper state!)
is_hopping, V_NEW = self._adjust_velocity(dKE)
if is_hopping: #### Now we really changed state!
steplog = '+'
hopping_no_field = True
if self.is_sh: #### global switch for SH
self.S = self._finalS
self.V = V_NEW
self._grd = gs
self._pot = es
else:
steplog = '+'
hopping_with_field = True
self.S = self._finalS
self._grd = gs
self._pot = es
else:
# print(self.S, self._finalS)
self._grd_finalS = gs
dKE = self._pot - es #### Kinetic E increase of _nuclei_ (can be negative if transition to upper state!)
# print("DKE = ", dKE)
is_hopping, V_NEW = self._adjust_velocity(dKE)
# print("is hopping: ", is_hopping)
# print("New velcoties should be: ", V_NEW)
if is_hopping: #### Now we really changed state!
steplog = '+'
hopping_no_field = True
if self.is_sh: #### global switch for SH
self.S = self._finalS
self.V = V_NEW
self._grd = gs
self._pot = es
# print("New velcoties are: ", self.V)
logstr = '{0} {1:>13.3f} {2: >6.4f} {3: >4d} {4: >4d} {5} {6}'.format(steplog,self.t*conv.au2fs,r,s0,self.S,hopping_no_field,hopping_with_field)
for pp in p:
logstr = logstr + ' {0:5.3f}'.format(pp.real)
logstr = logstr + '\n'
self.logSwitch.append(logstr)
def _adjust_velocity(self,dKE,**opts):
"""
Adjust velocity to conserve total energy if hopping is possible
Return True/False, V
dKE -- float, kinetic _nuclear_ energy increase
"""
# Skip if hops are not allowed (adiabatic case)
if self.no_hops:
return False, self.V
KE = self.kinetic_energy()
if KE + dKE > 0: #### Enough KE for hopping is available
# This condition will not always guarantee a hopping event, e.g. if the
# kinetic energy needs to be available in a particular direction.
# Therefore the final decision is made below according to the
# direction in which the velocity vector needs to be adjusted.
if self.rescaling == 'nac':
# NAC direction based decision (standard in TullyFS)
d = self.f_NAC(self.R,self.V,self.t+self.dt,self.S,self._finalS) #### NAC vector
elif self.rescaling == 'grd':
# Gradient difference direction based decision
d = self._grd - self._grd_finalS
else:
raise ValueError('self.rescaling not valid: '+self.rescaling)
#! If other possibilities are implemented we'll add "if" statements or specific subroutines
u = d / mt.sqrt( np.dot(d,d) ) #### direction of NAC vector
a = 0.5*np.dot( u * self.M, u)
b = np.dot( self.V * self.M, u)
# We are solving for a*k^2 + b*k - dKE = 0
# k is the needed factor such that v_NEW = v_OLD + k*u
tmp = b*b + 4.0*a*dKE
if tmp < 0.0: #### Not enough KE in the NAC direction for hopping
return False, self.V
else:
k1 = 0.5*( -b + np.sqrt(tmp) )/a
k2 = 0.5*( -b - np.sqrt(tmp) )/a
# Each solution leads to the same final kinetic energy. We pick the
# solution that results in the smallest momentum change
if abs(k1) < abs(k2):
k = k1
else:
k = k2
V_NEW = self.V + k*u
return True, V_NEW
else: #### NOT enough KE for hopping available
return False, self.V
[docs]
def kinetic_energy(self):
return 0.5*np.dot( self.V * self.M, self.V )
[docs]
def init_log(self):
self.logR = [self.R,] #### List of positions along trajectory
self.logV = [self.V,] #### List of velocities along trajectory
self.logC = [self.C,] #### List of coefficient vectors
self.logS = [self.S,] #### List of coefficient vectors
self.logt = [self.t,] #### List of times along the simulation
self.logKE = [self.ekin,] #### List of Total Energy along trajectory
self.logPE = [self._pot] #### List of Total Energy along trajectory
self.logH = [self._H0ad,] #### diagonals of the Hamiltonian (ad PES)
self.logWarning = [] #### List of warning messages
self.logSwitch = []
self.logPartial = []
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
self.logCrossings = [[]]
[docs]
def log(self):
self.logR.append( self.R )
self.logV.append( self.V )
self.logC.append( self.C )
self.logS.append( self.S )
self.logt.append( self.t )
self.logKE.append( self.ekin )
self.logPE.append( self._pot )
self.logH.append( self._H0ad )
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
self.logCrossings.append(self.crossings)
if self.getPartial is not None:
self.logPartial.append(self.partialCharges)
[docs]
def log_to_file(self,**opts):
"""
Write all log information to a set of files
"""
if not os.path.exists(self.DIR):
os.mkdir(self.DIR)
self._R_to_file() #### positions
self._V_to_file() #### velocities
self._C_to_file() #### quantum coefficients
self._P_to_file() #### quantum populations
self._S_to_file() #### quantum state
self._E_to_file() #### kinetic, potential and total energies
self._Warning_to_file() #### Warning messages
self._AdiabPot_to_file()
self._Switch_to_file()
self._NAC_to_file()
if self.trivialCrossing == 'overlap' and self.f_overlap is not None:
self._Crossings_to_file()
self._crossPot_to_file()
if self.getPartial is not None:
self._partial_to_file()
def _issue_warning(self,msg):
"""
Write warning message to warning log-list
"""
strng = '%8i %7.4f %s\n' % ( self._stepnum,self.t,msg.rstrip() )
self.logWarning.append( strng )
def _Switch_to_file(self):
with open(self.DIR + '/Switch.log','w') as f:
f.write('# stepTaken time(fs) r s0 s hoppingNoField hoppingWiField pp\'s ... \n')
for L in self.logSwitch:
f.write(L)
f.flush()
def _Crossings_to_file(self):
with open(self.DIR + '/trivialCrossings.log','w') as f:
f.write('# step from state -> to state\n')
for i, cs in enumerate(self.logCrossings):
for c in cs:
f.write( "%f %d %d\n" % (self.logt[i], c[0], c[1]))
f.flush()
def _NAC_to_file(self):
with open(self.DIR + '/NAC.log','w') as f:
for i,H in enumerate(self.logH):
for k in range(self.NStates):
for l in range(self.NStates):
f.write( '%4i %4i %11.4f' % (k,l,(self.logt[i])) )
f.write(' %13.8f %13.8f' % ( H[k,l].real, H[k,l].imag) )
f.write('\n')
f.flush()
def _AdiabPot_to_file(self):
with open(self.DIR + '/V_ad.log','w') as f:
for i,H in enumerate(self.logH):
Vad = np.diag( H + self._EREF )
f.write( '%11.4f ' % ((self.logt[i]), ) )
for v in Vad:
f.write(' %13.8f' % (v.real,) )
f.write('\n')
f.flush()
def _crossPot_to_file(self):
with open(self.DIR + '/V_cross.log','w') as f:
for i,H in enumerate(self.logH):
idx = np.arange(self.NStates)
for j,cs in enumerate(self.logCrossings[i:0:-1]):
idxOld = idx.copy()
for c in cs:
idx[c[1]] = idxOld[c[0]]
Vad = np.diag( H + self._EREF )[idx]
f.write( '%11.4f ' % ((self.logt[i]), ) )
for v in Vad:
f.write(' %13.8f' % (v.real,) )
f.write('\n')
f.flush()
def _R_to_file(self):
f=open(self.DIR + '/R.log','w')
for i,R in enumerate(self.logR):
f.write( '%9.6f' % ((self.logt[i]),) )
for r in R:
f.write(' %12.6f' % (r,) )
f.write('\n')
f.flush()
f.close()
def _V_to_file(self):
f=open(self.DIR + '/V.log','w')
for i,V in enumerate(self.logV):
f.write( '%9.6f' % ((self.logt[i]),) )
for v in V:
f.write(' %12.6f' % (v,) )
f.write('\n')
f.flush()
f.close()
def _C_to_file(self):
f=open(self.DIR + '/C.log','w')
for i,C in enumerate(self.logC):
f.write( '%9.6f' % ((self.logt[i]),) )
for c in C:
f.write(' %8.5f %8.5f' % (c.real, c.imag) )
f.write('\n')
f.flush()
f.close()
def _P_to_file(self):
f=open(self.DIR + '/P.log','w')
for i,C in enumerate(self.logC):
f.write( '%10.6f' % ((self.logt[i]),) )
P = np.abs(C*C.conj())
N = P.sum()
for p in P:
f.write(' %8.5f' % (p,) )
f.write(' %8.5f' % (N,) )
f.write('\n')
f.flush()
f.close()
def _S_to_file(self):
f=open(self.DIR + '/S.log','w')
for i,S in enumerate(self.logS):
f.write( '%9.6f %4i\n' % ((self.logt[i]), S) )
f.flush()
f.close()
def _E_to_file(self):
f=open(self.DIR + '/E.log','w')
for i in range(len(self.logt)):
f.write( '%9.6f %12.6f %12.6f %12.6f\n' % ((self.logt[i]),
self.logKE[i],
self.logPE[i],
self.logKE[i] + self.logPE[i] )
)
f.flush()
f.close()
def _Warning_to_file(self):
f=open(self.DIR + '/Warnings.log','w')
for msg in self.logWarning:
f.write( msg.rstrip() + '\n' )
f.flush()
f.close()
def _partial_to_file(self):
with open(self.DIR + '/partial.log','w') as f:
for i, ps in enumerate(self.logPartial):
f.write( '%11.4f ' % ((self.logt[i]), ) )
for p in ps:
f.write(' %13.8f' % (p) )
f.write('\n')
f.flush()
f.close()
def _switch_state(p,r):
"""
Determine the final state to jump to from an array of probabilities and a rand num.
Return -1 in case no jump takes place, otherwise return final state number
p -- array with probabilities, e.g. [0.1,0.15,0,0,0.3]
r -- random number [0,1]
"""
if np.sum(p) > 1.0:
raise ValueError('Total switching prob. > 1')
sp = 0.0
for i,pi in enumerate(p):
sp = sp + pi
if r < sp: return i
return -1
[docs]
def dot(K,L):
if len(K)==len(L) and len(K)!=0:
return sum([K[n]*L[n] for n in range(len(K))])
else:
return 0
#######################################################################################
#### Tully's model A (single avoided crossing) ####
#### JCP 93, 1061 (1990) ####
#######################################################################################
# Diabatic model
#@np.vectorize
[docs]
def W11_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
if x >= 0:
return A*(1.0-np.exp(-B*x))
if x < 0:
return -A*(1.0-np.exp(B*x))
#@np.vectorize
[docs]
def W22_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return -W11_TA(x,A,B,C,D)
#@np.vectorize
[docs]
def W12_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return C*np.exp(-D*x**2)
#@np.vectorize
[docs]
def W_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
W = np.empty((2,2),float)
W[0,0] = W11_TA(x,A,B,C,D)
W[1,1] = -W[0,0]
W[0,1] = W12_TA(x,A,B,C,D)
W[1,0] = W[0,1]
return W
#@np.vectorize
[docs]
def dq_W11_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
# Derivative of diabatic matrix element wrt all coordinates (only one in this model)
# All derivatives are returned by convention in a list
if x >= 0:
return np.array( (A*B*np.exp(-B*x) ),float )
if x < 0:
return np.array( (A*B*np.exp(B*x) ),float )
#@np.vectorize
[docs]
def dq_W22_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return -1.0*dq_W11_TA(x,A,B,C,D)
#@np.vectorize
[docs]
def dq_W12_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return np.array( ( -C*np.exp(-D*x**2)*2.0*D*x ),float )
#######################################################################################
#######################################################################################
## Run a set of trajectories in Tully's model A at fixed incoming momentum (KE) ##
#######################################################################################
[docs]
def return_TullyA_traj(p=15.0):
"""
Return a trajectory object ready to be run for Tully's model A
p -- initial momentum in the positive direction
Trajectories start at R0 = -5
"""
from CDTK.Models.DiabaticToAdiabatic2States import DiabaticToAdiabatic2States as di2ad
mod=di2ad(W11_TA,W22_TA,W12_TA,dq_W11_TA,dq_W22_TA,dq_W12_TA)
v1 = mod.ret_V1()
v2 = mod.ret_V2()
dqv1 = mod.ret_dqV1()
dqv2 = mod.ret_dqV2()
d12 = mod.ret_NAC()
def egrad(x,v,t,s):
if s==0:
return v1(x[0]),dqv1(x[0])
elif s==1:
return v2(x[0]),dqv2(x[0])
else:
raise ValueError
def nac(x,v,t,i,j):
return d12(x[0])
def w(x,v,t):
W = np.zeros((2,2),float)
W[0,0] = v1(x[0])
W[1,1] = v2(x[0])
return W
def d(x,v,t):
D = np.zeros((2,2),complex)
D[0,1] = -1.0j*np.dot( d12(x[0]), v[0] )
D[1,0] = -D[0,1]
return D
M = 2000.0 #### particle's mass
R0 = -5.0
V = p/M
traj = Trajectory_SH(1,2)
traj.M[0] = M
traj.R[0] = R0
traj.V[0] = V
traj.f_EGrad = egrad
traj.f_NAC = nac
traj.f_W = w
traj.f_D = d
return traj
[docs]
def run_TullyA(p,ntraj=10,rseed=None,do_save=False,dt=10.0):
"""
Run a set of trajectories in Tully's model A
p -- incoming momentum (a.u.)
ntraj -- number of averaged trajectories
rseed -- int; seed to start random number generator
"""
from CDTK.Models.DiabaticToAdiabatic2States import DiabaticToAdiabatic2States as di2ad
mod=di2ad(W11_TA,W22_TA,W12_TA,dq_W11_TA,dq_W22_TA,dq_W12_TA)
if type(1) == type(rseed):
np.random.seed(rseed)
v1 = mod.ret_V1()
v2 = mod.ret_V2()
dqv1 = mod.ret_dqV1()
dqv2 = mod.ret_dqV2()
d12 = mod.ret_NAC()
def egrad(x,v,t,s):
if s==0:
return v1(x[0]),dqv1(x[0])
elif s==1:
return v2(x[0]),dqv2(x[0])
else:
raise ValueError
def nac(x,v,t,i,j):
return d12(x)
def w(x,v,t):
W = np.zeros((2,2),float)
W[0,0] = v1(x[0])
W[1,1] = v2(x[0])
return W
def d(x,v,t):
D = np.zeros((2,2),complex)
D[0,1] = -1.0j*np.dot( d12(x[0]), v[0] )
D[1,0] = -D[0,1]
return D
M = 2000.0 #### particle's mass
R0 = -5.0
V = p/M
pop_0_L = 0
pop_0_R = 0
pop_1_L = 0
pop_1_R = 0
for nt in range(ntraj):
traj = Trajectory_SH(1,2,dt=dt)
traj.M[0] = M
traj.R[0] = R0
traj.V[0] = V
traj.f_EGrad = egrad
traj.f_NAC = nac
traj.f_W = w
traj.f_D = d
traj.integrate(tfinal=6000.0)
if do_save: traj.log_to_file()
if traj.S == 0 and traj.R[0] < 0.0: pop_0_L += 1
if traj.S == 0 and traj.R[0] > 0.0: pop_0_R += 1
if traj.S == 1 and traj.R[0] < 0.0: pop_1_L += 1
if traj.S == 1 and traj.R[0] > 0.0: pop_1_R += 1
trans = np.array([pop_0_L,pop_0_R,pop_1_L,pop_1_R],float)/ntraj
return trans
[docs]
def run_TullyA_Fig4(pmin=1.0,pmax=30.0,dp=1.0,ntraj=100,rseed=None,dt=10.0):
"""
Reproduce Fig. 4 in Tully's 1990 paper
Note: 2000 trajectories per momentum were used there
"""
print("Do Tully Fig 4, ntraj=", ntraj)
if type(1) == type(rseed):
np.random.seed(rseed)
pvec = np.arange(pmin,pmax,dp)
for p in pvec:
r = run_TullyA(p,ntraj=ntraj,dt=dt)
out = ' %5.2f %6.3f %6.3f %6.3f %6.3f\n' % (p,r[0],r[1],r[2],r[3])
sys.stdout.write(out)
sys.stdout.flush()
#######################################################################################
#### Tully's model B (dual avoided crossing) ####
#### JCP 93, 1061 (1990) ####
#######################################################################################
# Diabatic model
#@np.vectorize
[docs]
def W11_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
return 0.0
#@np.vectorize
[docs]
def W22_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
return -A*np.exp(-B*x**2) + E0
#@np.vectorize
[docs]
def W12_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
return C*np.exp(-D*x**2)
#@np.vectorize
[docs]
def W_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
W = np.empty((2,2),float)
W[0,0] = W11_TB(x,A,B,E0,C,D)
W[1,1] = W22_TB(x,A,B,E0,C,D)
W[0,1] = W12_TB(x,A,B,E0,C,D)
W[1,0] = W[0,1]
return W
#@np.vectorize
[docs]
def dq_W11_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
# Derivative of diabatic matrix element wrt all coordinates (only one in this model)
# All derivatives are returned by convention in a list
return 0.0
#@np.vectorize
[docs]
def dq_W22_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
return np.array( ( A*np.exp(-B*x**2)*2.0*B*x ),float )
#@np.vectorize
[docs]
def dq_W12_TB(x,A=0.1,B=0.28,E0=0.05,C=0.015,D=0.06):
return np.array( ( -C*np.exp(-D*x**2)*2.0*D*x ),float )
#######################################################################################
#######################################################################################
## Run a set of trajectories in Tully's model B at fixed incoming momentum (KE) ##
#######################################################################################
[docs]
def return_TullyB_traj(p=15.0):
"""
Return a trajectory object ready to be run for Tully's model B
p -- initial momentum in the positive direction
Trajectories start at R0 = -5
"""
from CDTK.Models.DiabaticToAdiabatic2States import DiabaticToAdiabatic2States as di2ad
mod=di2ad(W11_TB,W22_TB,W12_TB,dq_W11_TB,dq_W22_TB,dq_W12_TB)
v1 = mod.ret_V1()
v2 = mod.ret_V2()
dqv1 = mod.ret_dqV1()
dqv2 = mod.ret_dqV2()
d12 = mod.ret_NAC()
def egrad(x,v,t,s):
if s==0:
return v1(x[0]),dqv1(x[0])
elif s==1:
return v2(x[0]),dqv2(x[0])
else:
raise ValueError
def nac(x,v,t,i,j):
return d12(x[0])
def w(x,v,t):
W = np.zeros((2,2),float)
W[0,0] = v1(x[0])
W[1,1] = v2(x[0])
return W
def d(x,v,t):
D = np.zeros((2,2),complex)
D[0,1] = -1.0j*np.dot( d12(x[0]), v[0] )
D[1,0] = -D[0,1]
return D
M = 2000.0 #### particle's mass
R0 = -5.0
V = p/M
traj = Trajectory_SH(1,2)
traj.M[0] = M
traj.R[0] = R0
traj.V[0] = V
traj.f_EGrad = egrad
traj.f_NAC = nac
traj.f_W = w
traj.f_D = d
return traj
[docs]
def run_TullyB(E,ntraj=10,rseed=None,do_save=False,dt=10.0):
"""
Run a set of trajectories in Tully's model B
E -- incoming energy (a.u.)
ntraj -- number of averaged trajectories
rseed -- int; seed to start random number generator
"""
from CDTK.Models.DiabaticToAdiabatic2States import DiabaticToAdiabatic2States as di2ad
mod=di2ad(W11_TB,W22_TB,W12_TB,dq_W11_TB,dq_W22_TB,dq_W12_TB)
if type(1) == type(rseed):
np.random.seed(rseed)
v1 = mod.ret_V1()
v2 = mod.ret_V2()
dqv1 = mod.ret_dqV1()
dqv2 = mod.ret_dqV2()
d12 = mod.ret_NAC()
def egrad(x,v,t,s):
if s==0:
return v1(x[0]),dqv1(x[0])
elif s==1:
return v2(x[0]),dqv2(x[0])
else:
raise ValueError
def nac(x,v,t,i,j):
return d12(x[0])
def w(x,v,t):
W = np.zeros((2,2),float)
W[0,0] = v1(x[0])
W[1,1] = v2(x[0])
return W
def d(x,v,t):
D = np.zeros((2,2),complex)
D[0,1] = -1.0j*np.dot( d12(x[0]), v[0] )
D[1,0] = -D[0,1]
return D
M = 2000.0 #### particle's mass
p = np.sqrt(2*M*E)
R0 = -8.0
V = p/M
pop_0_L = 0
pop_0_R = 0
pop_1_L = 0
pop_1_R = 0
for nt in range(ntraj):
traj = Trajectory_SH(1,2,dt=dt)
traj.M[0] = M
traj.R[0] = R0
traj.V[0] = V
traj.f_EGrad = egrad
traj.f_NAC = nac
traj.f_W = w
traj.f_D = d
traj.integrate(tfinal=6000.0)
if do_save: traj.log_to_file()
if traj.S == 0 and traj.R[0] < 0.0: pop_0_L += 1
if traj.S == 0 and traj.R[0] > 0.0: pop_0_R += 1
if traj.S == 1 and traj.R[0] < 0.0: pop_1_L += 1
if traj.S == 1 and traj.R[0] > 0.0: pop_1_R += 1
trans = np.array([pop_0_L,pop_0_R,pop_1_L,pop_1_R],float)/ntraj
return trans
[docs]
def run_TullyB_Fig5(lemin=-4.0,lemax=1.0,de=0.2,ntraj=2000,rseed=None,dt=10.0):
"""
Reproduce Fig. 5 in Tully's 1990 paper
Note: 2000 trajectories per momentum were used there
"""
print("Do Tully Fig 5, ntraj=", ntraj)
if type(1) == type(rseed):
np.random.seed(rseed)
evec = np.arange(lemin,lemax,de)
for le in evec:
E = np.exp(le)
r = run_TullyB(E,ntraj=ntraj,dt=dt)
out = ' %5.2f %8.6f %6.3f %6.3f %6.3f %6.3f\n' % (le,E,r[0],r[1],r[2],r[3])
sys.stdout.write(out)
sys.stdout.flush()
######################################################################################
#### Test of a single trajectory ####
######################################################################################
if __name__ == "__main__":
# run_TullyA(8.5,ntraj=1,rseed=1,do_save=True,dt=10.0)
n = 2000
model = 'A'
pmin = 1.0
pmax = 29.0
if len(sys.argv) < 2:
print("Using default: Model ", model, ", ", n, " trajectories")
elif len(sys.argv) < 3:
model = sys.argv[1]
print("Using : Model ", model ,", ", n, " trajectories")
elif len(sys.argv) < 4:
model = sys.argv[1]
n = int(sys.argv[2])
print("Using : Model ", model ,", ", n, " trajectories")
elif len(sys.argv) < 5:
model = sys.argv[1]
n = int(sys.argv[2])
pmin = float(sys.argv[3])
print("Using : Model ", model ,", ", n, " trajectories")
print("Starting at pmin = ", pmin)
else:
model = sys.argv[1]
n = int(sys.argv[2])
pmin = float(sys.argv[3])
pmax = float(sys.argv[4])
print("Using : Model ", model ,", ", n, " trajectories")
print("Starting at pmin = ", pmin)
print("Ending before pmax = ", pmax)
if model == 'A':
run_TullyA_Fig4(ntraj=n, pmin=pmin, pmax=pmax)
elif model == 'B':
run_TullyB_Fig5(ntraj=n, lemin=pmin, lemax=pmax)
else:
print("Incorrect model chosen")
# r = run_TullyA(17,ntraj=50,rseed=1)
# out = '%6.3f %6.3f %6.3f %6.3f\n' % (r[0],r[1],r[2],r[3])
# sys.stdout.write(out)