#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
#! /usr/bin/env python
import sys
import string
import os
import random
import math
import numpy as np
import CDTK.Tools.Conversion as conv
import CDTK.Dynamics.MDIntegrators as mdi
from scipy.integrate import ode
from CDTK.Tools.Mathematics import ovrMatrix2
pid2 = math.pi/2.0
HUGE = 999999999999.9
TINY = 1.e-12
SMALL = 1.e-3
[docs]
class SurfaceHoppingFS(object):
"""
Control a trajectory surface-hopping evolution
"""
def __init__(self,**opts):
"""
Init a SurfaceHoppingFS object
optional arguments:
- initial_state: initially populated electronic state
default = 0
- nstates: number of electronic states
default = 1
- N_q : quantum time steps inside a classical step
default = 10
- external_E : function returning array of adiabatic energies
default = None
- external_G : function returning array of adiabatic gradients
default = None
- external_D : function returning array of non-adiabatic spatial derivative couplings
default = None
- external_VD : function returning array of non-adiabatic couplings
default = None
- external_DPL: function returning array of electronic dipole vectors
default = None
- gradient_method : method for adiabatic gradients
default = numerical
optional : analytical
- cplVD_method : method for the non-adiabatic couplings
default = wavefunction_overlap
optional : explicit_d_matrix
- simbox: SimulationBox object
default = None
- dxgrad: step for fallback numerical differenciation if no gradients
default = 0.001, [au]
- lockstate: whether to disable electronic transitions
default = False
- _stepnum_intp_ref: classical step number reference for interpolation procedure
initialized from step zero
- ishop_dec: decision for a hop. decision made by switch_condition
default = False
- cplBlock_radius: radius [half width] of the block of electronic states
considered for wavefunction overlap calculation
default = 1
- is_hfkoopmans: whether applying the Hartree-Fock-Koopmans method
default = False
- is_DoubleHole: whether applying electronic structure computation in
double hole configuration space
default = False
- is_switch_field: whether turning on the electric field
default = False
- cplFIELD_method: method for the field coupling
default = None
options : [direct_coupling] [dressed_state_coupling]
- MuE_MatElt: the electric field dipole couplings
"""
self.initial_state = opts.get('initial_state',0)
self.nstates = opts.get('nstates',1)
self.N_q = opts.get('n_q',10)
self.external_E = opts.get('external_E',None)
self.external_G = opts.get('external_G',None)
self.external_D = opts.get('external_D',None)
self.external_VD = opts.get('external_VD',None)
self.external_DPL = opts.get('external_DPL',None)
self.gradient_method = opts.get('gradient_method','numerical')
self.cplVD_method = opts.get('cplVD_method','wavefunction_overlap')
self.cplBlock_radius = opts.get('cplBlock_radius',1)
self.simbox = opts.get('simbox',None)
self.dxgrad = opts.get('dxgrad',0.001)
self.state = self.initial_state
self.lockstate = opts.get('lockstate',False)
self._logstep = {} # step info is stored here for logging, etc.
self.ishop_dec = False
self.is_hfkoopmans = opts.get('HFKoopmans_method',False)
self.is_DoubleHole = opts.get('DoubleHole',False)
self.is_switch_field = opts.get('switch_field',False)
self.cplFIELD_method = opts.get('cplFIELD_method',None)
self.MuE_MatElt = np.zeros((self.nstates,self.nstates),float)
self.delta_e = np.zeros((self.nstates,self.nstates),float)
if self.gradient_method == 'numerical':
self.delta_g = np.zeros(self.nstates,float)
self.dEdtmax = np.zeros(self.nstates,float)
self.dEdt = np.zeros(self.nstates,float)
self._stepnum_intp_ref = 0
self.is_hopp = False
self.hop_state = 0
[docs]
def get_funcE(self,**opts):
"""
Return the energy function to be used from a SimulationBox object
"""
if self.external_E is None:
raise ValueError('external energy function not defined')
def f(x):
e = self.external_E(x)[self.state]
return e
return f
[docs]
def get_funcG(self):
"""
Return the gradient function of current electronic state
to be used from a SimulationBox object
optional arguments:
nstates_g -- number of electronic states requiring analytical gradient
state_g -- the certain designated electronic state(s) requiring analytical gradient
currently for ab initio analytical gradient method
state_type -- specify the type of electronic state that requires analytical gradient
options: current_electronic_state - for SimulationBox object
"""
def f(x,**opts):
nstates_g = opts.get('nstates_g',1)
state_g = opts.get('state_g',None)
state_type = opts.get('state_type',None)
state_g = [self.state]
e, g = self.gradient_FS(x, state_g=state_g, state_type='current_electronic_state')
return g
return f
[docs]
def get_funcVD(self,a_x,**opts):
"""
Return the coupling matrix to be used from a SurfaceHoppingFS object
optional arguments:
state_block - electronic state block for the Hartree-Fock-Koopmans method
"""
state_block = opts.get('state_block',[])
if self.is_hfkoopmans:
d = self.external_D(a_x, state_block)
else:
d = self.external_D(a_x)
ncoords = self.simbox.V.size
nvec_v = 1 # velocity vector
nvec_d = self.nstates * self.nstates # non-adiabatic spatial derivative coupling vectors
veclen = ncoords
vd = np.zeros((nvec_v, nvec_d),float)
vd = ovrMatrix2(self.simbox.V, d, nvec_v, nvec_d, veclen, vd)
vd.shape = (self.nstates, self.nstates)
return vd
[docs]
def get_funcDPL(self,a_x,**opts):
"""
Return the dipole function to be used from a SurfaceHoppingFS object
"""
if self.external_DPL is None:
raise ValueError('external dipole function not defined')
dipole = self.external_DPL(a_x)
return np.array(dipole)
[docs]
def gradient_FS(self,a_x,**opts):
"""
Return the energy and gradient on the current electronic state
A function for the fewest switching FS surface hopp method
Provided auxillary data for the switch_condition function
Input arguments:
a_x -- positions of the system
Returns:
e_adiabatic -- potential energy
g_adiabatic -- energy gradient
self.delta_e -- energy difference
self.delta_g -- gradient difference
optional arguments:
nstates_g -- number of electronic states requiring analytical gradient
state_g -- the certain designated electronic states requiring analytical gradient
currently for ab initio analytical gradient method
state_type -- specify the type of electronic state that requires analytical gradient
options: current_electronic_state - for SimulationBox object
hopping_electronic_state - for switch_condition function
"""
nstates_g = opts.get('nstates_g',1)
state_g = opts.get('state_g',None)
state_type = opts.get('state_type',None)
# get list of adiabatic energies
if self.is_DoubleHole:
if self.simbox._stepnum == 0:
self.e_adiabatic = self.external_E(a_x, state_g=state_g, nstates_g=len(state_g), init=True)
else:
self.e_adiabatic = self.external_E(a_x, state_g=state_g, nstates_g=len(state_g))
else:
self.e_adiabatic = self.external_E(a_x)
# get adiabatic gradients
if self.gradient_method == 'analytical':
if self.external_G is not None:
if state_type == 'current_electronic_state':
self.g_adiabatic = self.external_G(a_x, state_g=state_g)
# only calculate the gradient on the is-going-to-hop electronic state
# when the switch_condition adjustment is required
elif self.is_hopp and state_type == 'hopping_electronic_state':
self.g_hopping = self.external_G(a_x, state_g=state_g, clean=True)
else:
raise ValueError('state_type not defined')
else:
raise ValueError('Analytical gradient calculation requires interface function external_G')
elif self.gradient_method == 'numerical':
if self.external_G is not None:
self.g_adiabatic = self.external_G(a_x)
else:
self.g_adiabatic = self._gradients_num(a_x)
else:
raise ValueError('Method for adiabatic gradient not defined')
# compute energy differeces for switch_condition
for i in range(self.nstates):
for j in range(self.nstates):
self.delta_e[i,j] = self.e_adiabatic[i] - self.e_adiabatic[j]
# compute gradient differences for switch_condition
if self.gradient_method == 'analytical':
# for analytical gradient method
# return adiabatic gradient
# -- on current electronic state
# compute adiabatic gradient difference in the case is_hopp is True
# -- between current and is-going-to-hop electronic states
current_state = 0
if self.is_hopp and state_type == 'hopping_electronic_state':
hopping_state = 0
self.delta_g = self.g_hopping - self.g_adiabatic
return self.e_adiabatic[self.state],self.g_adiabatic[current_state]
elif self.gradient_method == 'numerical':
# for numerical gradient method
# return adiabatic gradient and gradient difference on all electronic states
self.delta_g = self.g_adiabatic - self.g_adiabatic[self.state]
return self.e_adiabatic[self.state],self.g_adiabatic[self.state]
def _is_Hopp(self,a_x0,**opts):
"""
Determinew the judgement of surface hopp. is_hopp
A function for the fewest switching FS surface hopp method
Input arguments:
a_x0 -- positions of the system
Returns
is_hopp -- judgement of hopp. before kinetic energy check - global
optional arguments
integrator -- ode integrator for quantum electronic SE
default: zvode
is_first -- whether to be the initialization step
default: false
is_restart -- whether to the initialization step for [restart] calculation
default: false
Implements Eqs. (15,19) in Tully ; JCP 93, 1061 (1990)
"""
self.is_hopp = False # reset hopp switch
integrator = opts.get('integrator', 'zvode')
is_first = opts.get('init', False)
is_restart = opts.get('restart',False)
if self.cplVD_method == 'wavefunction_overlap':
self._is_Hopp_wavefunction_overlap(a_x0,
init=is_first,
restart=is_restart)
elif self.cplVD_method == 'explicit_d_matrix':
self._is_Hopp_explicit_d_matrix (a_x0,
init=is_first,
restart=is_restart)
else:
raise ValueError('Method for coupling matrix calculation not defined')
def _is_Hopp_wavefunction_overlap(self,a_x0,**opts):
"""
Performs _is_Hopp procedure for wavefunction_overlap method of
non-adiabatic coupling matrix
"""
integrator = opts.get('integrator', 'zvode')
is_first = opts.get('init', False)
is_restart = opts.get('restart',False)
if is_first:
# electronic amplitude
self.c = np.zeros(self.nstates, dtype=complex)
# non-adiabatic couplings
self.Vd = np.zeros((self.nstates,self.nstates),dtype=float)
# parameters to start SE integration from TIME=0
self.c[self.initial_state] = 1.0
if self.external_VD is None:
raise ValueError('Wavefunction overlap calculation requires interface function external_VD')
self.a_x0_TIME = a_x0.copy() # X at TIME=0
self.energy = self.e_adiabatic.copy() # E at TIME=0
# global energy shift to avoid fast phase rotation in the integration
if not is_restart:
self.energy_reference = np.mean(self.energy)
elif is_restart:
self.energy_reference = self.simbox.energy_reference.copy()
# the quantum step
self.DT_q = self.simbox.DT / float(self.N_q)
# for restarting calculation
if is_restart:
self.state = self.simbox.state_save
self.c = self.simbox.c_wavefunc_save.copy()
# initialize the ode integrator for one quantum step of electronic SE
self.TIME_quantum = self.simbox.TIME
if integrator == 'zvode':
f = self._get_func()
self.q = ode(f)
self.q.set_integrator('zvode', nsteps=500, method='bdf')
self.q.set_initial_value(self.c, self.TIME_quantum)
else:
raise ValueError('ODE integrator not defined')
return
# end [if is_first]
# perform adiabatic dynamics calculation.
# hopping between electronic states is locked.
if self.lockstate:
self._logstep['delta_e'] = 0.0
self._logstep['dEdt'] = 0.0
return
# An integration step of classical nuclei motion for classical interval DT is done in SimulationBox
# x(t) -> x(t+dt)
# v(t) -> v(t+dt)
# Do integration of quantum electronic amplitude over DT with quantum interval DT_q
self.TIME_quantum = self.simbox.TIME
self.TIME_quantum_final = self.simbox.TIME + self.simbox.DT
# quantities to be interpolated for integration
# e -- potential energy
self.e_TIME = self.energy.copy() # E at TIME
self.energy = self.e_adiabatic.copy() # E at TIME+DT updated in external_G. record for t+dt->(t+dt)+dt
self.e_TIMEpDT = self.energy.copy() # E at TIME+DT
# Vd -- non-adiabatic coupling V.d
# The first classical step serves preparing VD at DT-DT/2 of quantum integration dt->(dt)+dt
# The quantum integration starts from the second classical step
#
# Construct the block of electronic states considered for non-adiabatic coupling
#
state_block = []
for state_s in range(self.nstates):
radius = abs(state_s - self.state)
if radius <= self.cplBlock_radius:
state_block.append(state_s)
#
if self.simbox._stepnum == 1:
self.a_x0_TIMEpDT = a_x0.copy() # X at 0+DT
if self.is_hfkoopmans or self.is_DoubleHole:
#MAX_RADIUS = 3 # Molcas parameter
#if self.cplBlock_radius > MAX_RADIUS:
# raise ValueError('Error: Hartree-Fock-Koopmans block radius too large')
self.Vd = self.external_VD(self.a_x0_TIME, self.a_x0_TIMEpDT, self.simbox.DT, self.simbox.phase_track, state_block, is_init=True, nstates=self.nstates, unit=self.simbox.unit) # VD at DT-DT/2 [wavefunction_overlap]
else:
self.Vd = self.external_VD(self.a_x0_TIME, self.a_x0_TIMEpDT, self.simbox.DT, self.simbox.phase_track, is_init=True, nstates=self.nstates, unit=self.simbox.unit) # VD at DT-DT/2 [wavefunction_overlap]
self.a_x0_TIME = self.a_x0_TIMEpDT.copy() # X updated at 0+DT for classical step (0+dt)->(0+dt)+dt
elif self.simbox._stepnum > 1:
self.Vd_TIME = self.Vd.copy() # V.D at TIME-DT/2 [wavefunction_overlap]
self.a_x0_TIMEpDT = a_x0.copy() # X at TIME+DT
if self.is_hfkoopmans or self.is_DoubleHole:
# In the case the block is shifted. recompute the wavefunction overlap for two molecular geometries
if self.ishop_dec:
self.Vd = self.external_VD(self.a_x0_TIME, self.a_x0_TIMEpDT, self.simbox.DT, self.simbox.phase_track, state_block, is_init=True, is_reinit=True, nstates=self.nstates, unit=self.simbox.unit) # V.D at (TIME+DT)-DT/2
else:
self.Vd = self.external_VD(self.a_x0_TIME, self.a_x0_TIMEpDT, self.simbox.DT, self.simbox.phase_track, state_block, nstates=self.nstates, unit=self.simbox.unit) # V.D at (TIME+DT)-DT/2
else:
self.Vd = self.external_VD(self.a_x0_TIME, self.a_x0_TIMEpDT, self.simbox.DT, self.simbox.phase_track, nstates=self.nstates, unit=self.simbox.unit) # V.D at (TIME+DT)-DT/2
self.a_x0_TIME = self.a_x0_TIMEpDT.copy() # X updated at TIME+DT for classical step (t+dt)->(t+dt)+dt
else:
raise ValueError('Error: incorrect step number')
self.Vd_TIMEpDT = self.Vd.copy() # V.D at TIME+DT
sys.stdout.write('\n')
# forced correction of wave function phase
# the absolute phase is unknown
if self.simbox._stepnum > 1:
for i in range(self.nstates):
for j in range(i+1,self.nstates):
r_ph = self.Vd_TIMEpDT[i,j] / self.Vd_TIME[i,j]
if (r_ph < 0.0) and (0.5 < abs(r_ph) < 2.0):
# assumed to be a matrix element with flipped sign
self.Vd_TIMEpDT[i,j] *= -1.0
self.Vd_TIMEpDT[j,i] *= -1.0
self.Vd[i,j] *= -1.0
self.Vd[j,i] *= -1.0
# electronic population A at TIME
A_TIME = np.zeros(self.nstates,float)
for i in range(self.nstates):
A_TIME[i] = abs(self.c[i].conjugate() * self.c[i])
# the surface hopp probability matrix g at the end of a classical step TIME+DT
# integrate within one classical step DT
prob_FS = np.zeros((self.nstates,self.nstates),float)
# re-initialize the integrator in the case the block is shifted due to hopping
# avoid stiffness to the integrator
if self.is_hfkoopmans:
if self.ishop_dec:
self.TIME_quantum = self.simbox.TIME
if integrator == 'zvode':
f = self._get_func()
self.q = ode(f)
self.q.set_integrator('zvode', nsteps=500, method='bdf')
self.q.set_initial_value(self.c, self.TIME_quantum)
while self.TIME_quantum < (self.TIME_quantum_final - SMALL):
# compute electronic amplitude c via integrating SE
self._step_quantum()
self.TIME_quantum = self.TIME_quantum + self.DT_q
# compute FS hopping probabilities via B from electronic amplitude and integrate over quantum steps within DT
time = self.TIME_quantum
prob_FS = prob_FS + self.get_prob_FS(time)
prob_FS = prob_FS * self.DT_q
prob_FS[self.state,:] = prob_FS[self.state,:] / A_TIME[self.state]
# adjust unphysical value due to num err
for j in range(self.nstates):
if prob_FS[self.state,j] > 1.0:
prob_FS[self.state,j] = 1.0
randnum = np.random.random()
norm_el = self._chk_norm_el(A_TIME)
self._logstep['prob_FS'] = prob_FS
self._logstep['randnum'] = randnum
self._logstep['norm_el'] = norm_el
# determine the pre-switching using FS formula before kinetic energy check in switch_condition
ipmax = -1
Sigma_prob_FS = 0.0
for i in range(self.nstates):
if i != self.state: # avoid selfjumps
Sigma_prob_FS = Sigma_prob_FS + prob_FS[self.state,i]
if Sigma_prob_FS > randnum:
self.is_hopp = True
ipmax = i
break
# compute dDE/dt and keep track of the maximum for every state
# only for numerical gradient
if self.gradient_method == 'numerical':
dEdt = abs(np.dot(self.delta_g,self.simbox.V))
dEdt[self.state] = 1.0 # avoid zerodivision
for i, dE in enumerate(dEdt):
if dE > self.dEdtmax[i]: self.dEdtmax[i] = dE
# save local information for velocity adjustment and K.E. check in switch_condition
# self.de_s
# self.dg_s
# is_hopp
delta_e_s = self.delta_e[:,self.state]
if self.is_hopp:
# save local information to be used from the trajectory object
self.hop_state = ipmax
self.de_s = delta_e_s[ipmax]
if self.gradient_method == 'analytical':
state_g = [ipmax]
self.gradient_FS(a_x0, state_g=state_g, state_type='hopping_electronic_state')
hopping_state = 0
self.dg_s = self.delta_g[hopping_state]
self.gradient_s = self.g_hopping[hopping_state]
elif self.gradient_method == 'numerical':
self.dg_s = self.delta_g[ipmax]
self.gradient_s = self.g_adiabatic[ipmax]
f = self.simbox.logfile
f.write('## Quantum electronic amplitude ')
f.write('%8.2f a.u. to %8.2f a.u. \n' % (
(self.TIME_quantum-self.simbox.DT),
self.TIME_quantum_final))
f.write('## Hopping probabilities \n')
for i, prob_FS_s in enumerate(self._logstep['prob_FS'][self.state]):
f.write(' %i %8.6e\n' % (
i, prob_FS_s))
f.write('## rand\n')
f.write(' %8.6f\n' % (
self._logstep['randnum']))
f.write('## Norm of electronic wavefunction \n')
f.write(str(self._logstep['norm_el'])+'\n')
self._logstep['delta_e'] = delta_e_s
if self.gradient_method == 'numerical':
self._logstep['delta_g'] = self.delta_g
self._logstep['dEdt'] = dEdt
def _is_Hopp_explicit_d_matrix(self,a_x0,**opts):
"""
Performs _is_Hopp procedure for explicit_d_matrix method of
non-adiabatic coupling matrix
"""
integrator = opts.get('integrator', 'zvode')
is_first = opts.get('init', False)
is_restart = opts.get('restart',False)
# coupling matrix element from electric field
if self.is_switch_field:
dipole = self.get_funcDPL(a_x0)
efield = self.simbox._fc.getElectricField(self.simbox.TIME)
# inner product. dipole and electric field
# The dipole elements [numpy array] are structured as
# | [0.0] [0.1] ... [0.n] |
# | ... [k.l] ... |
# | ... ... ... |
# each dipole element [k.l] is n-dimensional vector for the dipole of k- and l-th electronic state.
# The field coupling matrix element [numpy array] are identically structured.
# each field coupling matrix element [k.l] is scalar for the coupling of k- and l-th electronic state.
self.MuE_MatElt = np.dot(dipole,efield)
# Construct the block of electronic states considered for non-adiabatic coupling
if self.is_hfkoopmans:
state_block = []
for state_s in range(self.nstates):
radius = abs(state_s - self.state)
if radius <= self.cplBlock_radius:
state_block.append(state_s)
if is_first:
# electronic amplitude
self.c = np.zeros(self.nstates, dtype=complex)
# non-adiabatic couplings
self.Vd = np.zeros((self.nstates,self.nstates),dtype=float)
# parameters to start SE integration from TIME=0
self.c[self.initial_state] = 1.0
if self.external_D is None:
raise ValueError('Explicit d matrix calculation requires interface function external_D')
if self.is_hfkoopmans:
self.Vd = self.get_funcVD(a_x0,state_block=state_block) # V.D at TIME=0
else:
self.Vd = self.get_funcVD(a_x0) # V.D at TIME=0
self.energy = self.e_adiabatic.copy() # E at TIME=0
# global energy shift to avoid fast phase rotation in the integration
if not is_restart:
self.energy_reference = np.mean(self.energy)
elif is_restart:
self.energy_reference = self.simbox.energy_reference.copy()
# the quantum step
self.DT_q = self.simbox.DT / float(self.N_q)
# for restarting calculation
if is_restart:
self.state = self.simbox.state_save
self.c = self.simbox.c_wavefunc_save.copy()
# interpolated electric field couplings
if self.is_switch_field:
if self.cplFIELD_method == 'direct_coupling':
self.Vd = self.Vd - 1.0j * self.MuE_MatElt
elif self.cplFIELD_method == 'dressed_state_coupling':
# adiabatic energy plus field coupling
# in adiabatic representation
EpFmat_adiabatic = np.diag(self.energy) - self.MuE_MatElt
# in dressed representation
# dress energy and Vd non-adiabatic coupling matrix
self.EpFmat_dressed, self.EpF_Umat = np.linalg.eig(EpFmat_adiabatic)
self.energy = self.EpFmat_dressed.copy()
# calculate gradient D_Umat
EpF_D_Umat = self.simbox._fc.getEpF_D_Umat(a_x0, self.nstates, self.dxgrad, self.external_E, self.MuE_MatElt)
# dress the Vd matrix
self.Vd = self.simbox._fc.getDressedVd(self.Vd, self.simbox.V, self.nstates, self.EpF_Umat, EpF_D_Umat)
# transform electronic wavefunction from adiabatic to dressed representation
self.c = np.dot(self.EpF_Umat, self.c)
else:
raise ValueError('Method for electric field coupling not defined')
# initialize the ode integrator for one quantum step of electronic SE
self.TIME_quantum = self.simbox.TIME
if integrator == 'zvode':
f = self._get_func()
self.q = ode(f)
self.q.set_integrator('zvode', nsteps=500, method='bdf')
self.q.set_initial_value(self.c, self.TIME_quantum)
else:
raise ValueError('ODE integrator not defined')
return
# end [if is_first]
# perform adiabatic dynamics calculation.
# hopping between electronic states is locked.
if self.lockstate:
self._logstep['delta_e'] = 0.0
self._logstep['dEdt'] = 0.0
return
# An integration step of classical nuclei motion for classical interval DT is done in SimulationBox
# x(t) -> x(t+dt)
# v(t) -> v(t+dt)
# Do integration of quantum electronic amplitude over DT with quantum interval DT_q
self.TIME_quantum = self.simbox.TIME
self.TIME_quantum_final = self.simbox.TIME + self.simbox.DT
# quantities to be interpolated for integration
# e -- potential energy
self.e_TIME = self.energy.copy() # E at TIME
self.energy = self.e_adiabatic.copy() # E at TIME+DT updated in external_G. record for t+dt->(t+dt)+dt
self.e_TIMEpDT = self.energy.copy() # E at TIME+DT
# Vd -- non-adiabatic coupling V.d
self.Vd_TIME = self.Vd.copy() # V.D at TIME [explicit_d_matrix]
if self.is_hfkoopmans:
self.Vd = self.get_funcVD(a_x0,state_block=state_block) # V.D at TIME+DT. explicitly mult velocity V
else:
self.Vd = self.get_funcVD(a_x0) # V.D at TIME+DT. explicitly mult velocity V
# interpolated electric field coupling
if self.is_switch_field:
if self.cplFIELD_method == 'direct_coupling':
self.Vd = self.Vd - 1.0j * self.MuE_MatElt
elif self.cplFIELD_method == 'dressed_state_coupling':
# adiabatic energy plus field coupling
# in adiabatic representation
EpFmat_adiabatic = np.diag(self.energy) - self.MuE_MatElt
# in dressed representation
# dress energy and Vd non-adiabatic coupling matrix
EpF_Umat_TIME = self.EpF_Umat.copy() # direct dUdt at TIME
self.EpFmat_dressed, self.EpF_Umat = np.linalg.eig(EpFmat_adiabatic)
EpF_Umat_TIMEpDT = self.EpF_Umat.copy() # direct dUdt at TIME+DT
# adjust phase
if EpF_Umat_TIME[0,0] * EpF_Umat_TIMEpDT[0,0] < 0.0:
EpF_Umat_TIME = - EpF_Umat_TIME
EpF_dUdt = (EpF_Umat_TIMEpDT - EpF_Umat_TIME) / self.simbox.DT # direct dUdt
self.energy = self.EpFmat_dressed.copy()
# calculate gradient D_Umat
#EpF_D_Umat = self.simbox._fc.getEpF_D_Umat(a_x0, self.nstates, self.dxgrad, self.external_E, self.MuE_MatElt)
# dress the Vd matrix
#self.Vd = self.simbox._fc.getDressedVd(self.Vd, self.simbox.V, self.nstates, EpF_Umat, EpF_D_Umat)
self.Vd = self.simbox._fc.getDressedVd_dUdt(self.Vd, self.nstates, self.EpF_Umat, EpF_dUdt)
# transform electronic wavefunction from adiabatic to dressed representation
if self.simbox._stepnum > 1:
self.c = np.dot(self.EpF_Umat, self.c)
self.Vd_TIMEpDT = self.Vd.copy() # V.D at TIME+DT
# electronic population A at TIME
A_TIME = np.zeros(self.nstates,float)
for i in range(self.nstates):
A_TIME[i] = abs(self.c[i].conjugate() * self.c[i])
# the surface hopp probability matrix g at the end of a classical step TIME+DT
# integrate within one classical step DT
prob_FS = np.zeros((self.nstates,self.nstates),float)
# re-initialize the integrator in the case the block is shifted due to hopping
# avoid stiffness to the integrator
if self.is_hfkoopmans:
if self.ishop_dec:
self.TIME_quantum = self.simbox.TIME
if integrator == 'zvode':
f = self._get_func()
self.q = ode(f)
self.q.set_integrator('zvode', nsteps=500, method='bdf')
self.q.set_initial_value(self.c, self.TIME_quantum)
while self.TIME_quantum < (self.TIME_quantum_final - SMALL):
# compute electronic amplitude c via integrating SE
self._step_quantum()
self.TIME_quantum = self.TIME_quantum + self.DT_q
# compute FS hopping probabilities via B from electronic amplitude and integrate over quantum steps within DT
time = self.TIME_quantum
prob_FS = prob_FS + self.get_prob_FS(time)
prob_FS = prob_FS * self.DT_q
prob_FS[self.state,:] = prob_FS[self.state,:] / A_TIME[self.state]
# adjust unphysical value due to num err
for j in range(self.nstates):
if prob_FS[self.state,j] > 1.0:
prob_FS[self.state,j] = 1.0
randnum = np.random.random()
norm_el = self._chk_norm_el(A_TIME)
self._logstep['prob_FS'] = prob_FS
self._logstep['randnum'] = randnum
self._logstep['norm_el'] = norm_el
# determine the pre-switching using FS formula before kinetic energy check in switch_condition
ipmax = -1
Sigma_prob_FS = 0.0
for i in range(self.nstates):
if i != self.state: # avoid selfjumps
Sigma_prob_FS = Sigma_prob_FS + prob_FS[self.state,i]
if Sigma_prob_FS > randnum:
self.is_hopp = True
ipmax = i
break
# compute dDE/dt and keep track of the maximum for every state
# only for numerical gradient
if self.gradient_method == 'numerical':
dEdt = abs(np.dot(self.delta_g,self.simbox.V))
dEdt[self.state] = 1.0 # avoid zerodivision
for i, dE in enumerate(dEdt):
if dE > self.dEdtmax[i]: self.dEdtmax[i] = dE
# save local information for velocity adjustment and K.E. check in switch_condition
# self.de_s
# self.dg_s
# is_hopp
delta_e_s = self.delta_e[:,self.state]
if self.is_hopp:
# save local information to be used from the trajectory object
self.hop_state = ipmax
self.de_s = delta_e_s[ipmax]
if self.gradient_method == 'analytical':
state_g = [ipmax]
self.gradient_FS(a_x0, state_g=state_g, state_type='hopping_electronic_state')
hopping_state = 0
self.dg_s = self.delta_g[hopping_state]
self.gradient_s = self.g_hopping[hopping_state]
elif self.gradient_method == 'numerical':
self.dg_s = self.delta_g[ipmax]
self.gradient_s = self.g_adiabatic[ipmax]
f = self.simbox.logfile
f.write('## Quantum electronic amplitude ')
f.write('%8.2f a.u. to %8.2f a.u. \n' % (
(self.TIME_quantum-self.simbox.DT),
self.TIME_quantum_final))
f.write('## Hopping probabilities \n')
for i, prob_FS_s in enumerate(self._logstep['prob_FS'][self.state]):
f.write(' %i %8.6e\n' % (
i, prob_FS_s))
f.write('## rand\n')
f.write(' %8.6f\n' % (
self._logstep['randnum']))
f.write('## Norm of electronic wavefunction \n')
f.write(str(self._logstep['norm_el'])+'\n')
self._logstep['delta_e'] = delta_e_s
if self.gradient_method == 'numerical':
self._logstep['delta_g'] = self.delta_g
self._logstep['dEdt'] = dEdt
# electric field coupling
if self.is_switch_field:
if self.cplFIELD_method == 'dressed_state_coupling':
# transform electronic wavefunction from dressed to adiabatic representation
# for analyse functions
self.c = np.dot(np.transpose(self.EpF_Umat), self.c)
def _step_quantum(self,**opts):
"""
Updates the electronic amplitudes c from integrating SE
the adiabatic couplings d
The RHS functions of the ODE are computed from _get_func
Input arguments:
self.e_TIME -- potential energies of each state at TIME
self.e_TIMEpDT -- potential energies of each state at TIME+DT
self.Vd_TIME -- non-adiabatic couplings of each state at TIME
self.Vd_TIMEpDT -- non-adiabatic couplings of each state at TIME+DT
self.c -- electronic amplitudes within the integration procedure at TIME
Returns:
self.c -- updated electronic amplitude at each quantum steps within TIME to TIME+DT
"""
while self.q.successful() and self.q.t < (self.TIME_quantum + self.DT_q):
self.q.integrate(self.q.t + self.DT_q)
self.c = self.q.y
if not self.q.successful():
raise ValueError('The integration over quantum step fails')
# dump current status
def _funcf(self,c,time):
"""
This function is a thin wrapper over RHS of SE
"""
# interpolate E and V.D
e_intm, Vd_intm = self.linear_interpolate(time,retEVd='returnAll')
e_intm = e_intm - self.energy_reference
Hmat = -1.0j * np.diag(e_intm) - Vd_intm
HmatC = np.dot(Hmat,c)
return HmatC
def _get_func(self):
"""
Return the RHS function of the SE ODE
"""
def f(time,c):
return self._funcf(c,time)
return f
[docs]
def get_prob_FS(self,time):
"""
Return the FS hopp probability matrix g
"""
Vd_intp = self.linear_interpolate(time,retEVd='returnVd')
Gmat = np.zeros((self.nstates,self.nstates),float)
Gmat[self.state,self.state] = 0.0 # avoid selfjumps!
for j in range(self.nstates):
if j != self.state:
BMatElt = -2.0 * self.c[self.state].conjugate() * self.c[j] * Vd_intp[self.state,j]
Gmat[self.state,j] = BMatElt.real
if Gmat[self.state,j] < 0.0:
Gmat[self.state,j] = 0.0 # adjust unphysical value
return Gmat
[docs]
def linear_interpolate(self,time,**opts):
"""
Returns interpolated E and V.D
optional arguments:
- retEVd: variables to be interpolated
default: returnAll - interpolated E and V.D
retuenVd - interpolated V.D
"""
retEVd = opts.get('retEVd','returnAll')
if self.cplVD_method == 'explicit_d_matrix':
# non-adiabatic couplings from explicit d matrix method
# from self.simbox.TIME to self.simbox.TIME+self.simbox.DT
t0 = self.simbox.TIME
#t1 = self.simbox.TIME + self.simbox.DT
delta_t = time - t0
# linear interpolation
if retEVd == 'returnAll':
e_intm = self.e_TIME + delta_t * ( self.e_TIMEpDT - self.e_TIME ) / (self.simbox.DT)
Vd_intm = self.Vd_TIME + delta_t * ( self.Vd_TIMEpDT - self.Vd_TIME ) / (self.simbox.DT)
return e_intm, Vd_intm
elif retEVd == 'returnVd':
Vd_intm = self.Vd_TIME + delta_t * ( self.Vd_TIMEpDT - self.Vd_TIME ) / (self.simbox.DT)
return Vd_intm
elif self.cplVD_method == 'wavefunction_overlap':
# non-adiabatic couplings from wavefunction overlap method
# E from self.simbox.TIME to self.simbox.TIME+self.simbox.DT
# Vd from self.simbox.TIME-self.simbox.DT/2.0 to self.simbox.TIME+self.simbox.DT/2.0
t0_e = self.simbox.TIME
#t1_e = self.simbox.TIME + self.simbox.DT
delta_t_e = time - t0_e
# linear interpolation
if retEVd == 'returnAll':
e_intm = self.e_TIME + delta_t_e * ( self.e_TIMEpDT - self.e_TIME ) / (self.simbox.DT)
# the interpolation of non-adiabatic coupling is invoked
# when and after classical integration proceeds to the second step
if self.simbox._stepnum == 1:
# interpolate through [dt-dt/2, dt+dt/2]
Vd_intm = self.linear_interpolate_Vd(time,is_init=True)
else:
# interpolate through [t, t+dt/2]
Vd_intm = self.linear_interpolate_Vd(time)
return e_intm, Vd_intm
elif retEVd == 'returnVd':
if self.simbox._stepnum == 1:
Vd_intm = self.linear_interpolate_Vd(time,is_init=True)
else:
Vd_intm = self.linear_interpolate_Vd(time)
return Vd_intm
[docs]
def linear_interpolate_Vd(self,time,**opts):
"""
Continuously interpolate the non-adiabatic coupling Vd after the initializing classical step.
adapted to the wavefunction overlap method.
Inputs:
time: current time of quantum integration
uses:
self.Vd_intp_TIME: non-adiabatic coupling at t
self.Vd_intp_TIMEpDT: non-adiabatic coupling at t+dt
self.Vd_TIMEpDT: non-adiabatic coupling at t+dt/2
Returns:
Vd_intm: the interpolated non-adiabatic coupling
optional arguments:
is_init: whether to be the initialization classical step
default - False
"""
is_init = opts.get('is_init',False)
if is_init:
# for the first step the non-adiabatic coupling is set to be constant
Vd_intm = self.Vd.copy()
if self.simbox._stepnum != self._stepnum_intp_ref:
self.Vd_intp_TIMEpDT = self.Vd.copy()
self._stepnum_intp_ref = self.simbox._stepnum
else:
if self.simbox._stepnum != self._stepnum_intp_ref:
# classical step is advanced. update interpolation references
self.Vd_intp_TIME = self.Vd_intp_TIMEpDT.copy()
self.Vd_intp_TIMEpDT = self.Vd_TIMEpDT.copy()
self._stepnum_intp_ref = self.simbox._stepnum
delta_t_vd = time - self.simbox.TIME
if delta_t_vd <= (self.simbox.DT / 2.0):
Vd_intm = self.Vd_intp_TIME + delta_t_vd * ( self.Vd_TIMEpDT - self.Vd_intp_TIME ) / (self.simbox.DT / 2.0)
else:
Vd_intm = self.Vd_TIMEpDT.copy()
return Vd_intm
[docs]
def switch_condition(self,**opts):
"""
Compute initial conditions on the new PES after a hop.
uses:
self.de_s: energy difference (Vf - Vi) between the two adiabatic states
self.dg_s: non-adiabatic coupling gradient.
this gradient may be different in different flavors
of switching algorithms.
Common conventions are:
dg_s = Grad(Vf - Vi)
dg_s = <f| Grad |i> (|i>:electronic wfn)
Returns: (xnew,vnew,ishop)
xnew: new positions
vnew: new velocities
ishop: True if switch took place
Implements Eq.(3.44) in Miller and George; JCP 56, 5637 (1972)
Implements Eqs.(25,29) in Stine et al. ; JCP 65, 3975 (1976)
"""
# self.dEdtmax[:] = 0.0 !! dont forget
de = self.de_s
dg = self.dg_s
pold = self.simbox.M * self.simbox.V
inv_mass = 1.0 / self.simbox.M
gnorm = dg/np.linalg.norm(dg)
a = np.dot(gnorm,inv_mass*pold)
b = np.dot(gnorm,inv_mass*gnorm)
c = 1.0 - 2.0*de*b/(a**2.0)
# Ek_gnorm is the kinetic energy along gnorm before the hopping attempt
# Ek_gnorm = a*(np.dot(gnorm,pold))/2.0
if c > 0: # enough kinetic energy in gnorm direction for hop
# if the kinetic energy,Ek_gnorm, is less than
# the gap between the two potential energy surfaces, the hop is rejected
pnew = pold - gnorm*a*(1.0 - c**0.5)/b
xnew = self.simbox.X
vnew = inv_mass * pnew
ishop = True
else:
vnew = self.simbox.V
xnew = self.simbox.X
ishop = False
if self.is_switch_field and self.cplFIELD_method == 'direct_coupling' or self.cplFIELD_method == 'dressed_state_coupling':
if self.hop_state > self.state:
# field coupling
# hopping upwards without energy conservation check
# acquiring energy from photons
vnew = self.simbox.V
xnew = self.simbox.X
ishop = True
self.ishop_dec = ishop # decision for a hop
return xnew,vnew,ishop
def _gradients_num(self,x,**opts):
"""
Fallback to numerical gradients
"""
dx = self.dxgrad # step for numerical differentiation
grads = np.zeros((self.nstates,x.size),float)
for i in range(x.size):
x1 = x.copy()
x1[i] = x1[i] + 0.5*dx
e1 = self.external_E(x1)
x0 = x.copy()
x0[i] = x0[i] - 0.5*dx
e0 = self.external_E(x0)
grads[:,i] = (e1-e0)/dx
return grads
def _chk_norm_el(self,pop_el):
"""
Returns norm of electronic wavefunction
"""
norm_el = 0.0
for pop in pop_el:
norm_el = norm_el + pop
return norm_el
def _update_logfile(self):
f = self.simbox.logfile
f.write('## E_i - E_state\n')
f.write(str(self._logstep['delta_e'])+'\n')
f.write('## d DeltaE/dt\n')
if self.gradient_method == 'numerical':
f.write(str(self._logstep['dEdt'])+'\n')
# f.write('## Hopping probabilities\n')
# f.write(str(self._logstep['prob_LZ'])+'\n')
# f.write('## ipmax prob rand\n')
# f.write(' %i %8.6f %8.6f\n' % (
# self._logstep['ipmax'],
# self._logstep['prob_max'],
# self._logstep['randnum']))
[docs]
def gradient_simple(self,a_x,**opts):
"""
Return the energy and gradient on the current electronic state
A function for the SE integration
Merely compute the energies and gradients without preparing data for switch_condition
Input arguments:
a_x -- positions of the system
Returns:
e_adiabatic -- potential energy
g_adiabatic -- energy gradient
"""
# get list of adiabatic energies
self.e_adiabatic = self.external_E(a_x)
# get adiabatic gradients
if self.external_G is not None:
self.g_adiabatic = self.external_G(a_x)
else:
self.g_adiabatic = self._gradients_num(a_x)
return self.e_adiabatic[self.state],self.g_adiabatic[self.state]