#* **************************************************************************
#*
#* 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
pid2 = math.pi/2.0
HUGE = 999999999999.9
TINY = 1.e-9
SMALL = 1.e-3
[docs]
class SurfaceHoppingLZ(object):
"""
Control a trajectory surface-hopping evolution
"""
def __init__(self,**opts):
"""
Init a SurfaceHoppingLZ object
optional arguments:
- initial_state: initially populated electronic state
default = 0
- nstates: number of electronic states
default = 1
- external_E : function returning array of adiabatic energies
default = None
- external_G : function returning array of adiabatic gradients
default = None
- gradient_method : method for adiabatic gradients
default = numerical
optional : analytical
- 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
"""
self.initial_state = opts.get('initial_state',0)
self.nstates = opts.get('nstates',1)
self.external_E = opts.get('external_E',None)
self.external_G = opts.get('external_G',None)
self.gradient_method = opts.get('gradient_method','numerical')
self.simbox = opts.get('simbox',None)
self.dxgrad = opts.get('dxgrad',0.001)
self.egap = opts.get('egap',0.005)
self.state = self.initial_state
self.lockstate = opts.get('lockstate',False)
self._logstep = {} # step info is stored here for logging, etc.
self.delta_e = np.zeros((self.nstates,self.nstates),float)
self.delta_e1 = np.zeros((self.nstates,self.nstates),float)
self.delta_e2 = np.zeros((self.nstates,self.nstates),float)
self.dEdtmax = np.zeros(self.nstates,float)
self.dEdt = np.zeros(self.nstates,float)
self.is_hopp = False
self.hop_state = 0
self._DEBUG = False
[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,**opts):
"""
Return the gradient function of current electronic state
to be used from a SimulationBox object
optional argument:
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):
nstates_g = opts.get('nstates_g',1)
state_g = opts.get('state_g',None)
state_type = opts.get('state_type',None)
e, g = self.gradient_LZ(x, state_type='all_electronic_states')
return g
return f
[docs]
def gradient_LZ(self,a_x,**opts):
"""
Return the energy and gradient on the current electronic state
If a surface hopp occurs the velocities will be rescaled accordingly
and the energy and gradient on the new surface provided.
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 argument:
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
default: all_electronic_states
options: current_electronic_state - for SimulationBox object
"""
nstates_g = opts.get('nstates_g',1)
state_g = opts.get('state_g',None)
state_type = opts.get('state_type','all_electronic_states')
self.is_hopp = False # reset hopp switch
# get list of adiabatic energies
e_adiabatic = self.external_E(a_x)
self.e_adiabatic = e_adiabatic.copy() # record the potential energy of each electronic state
# get adiabatic gradients
#if self.external_G is not None:
# g_adiabatic = self.external_G(a_x)
#else:
# g_adiabatic = self._gradients_num(a_x)
if self.gradient_method is 'analytical':
if self.external_G is not None:
if state_type is 'current_electronic_state':
state_g = [self.state]
g_adiabatic = self.external_G(a_x, state_g=state_g)
return g_adiabatic
# LZSH method requires adiabatic gradients for all electronic states
elif state_type is 'all_electronic_states':
nstates_g = self.nstates
state_g = range(nstates_g)
g_adiabatic = self.external_G(a_x, nstates_g=nstates_g, state_g=state_g)
# correct the possible error in Molcas. the error makes the gradients for all electronic states computed as the same
if np.linalg.norm(g_adiabatic[0] - g_adiabatic[1]) < SMALL:
raise ValueError('Analytical gradient calculation incorrect')
else:
raise ValueError('state_type not defined')
else:
raise ValueError('Analytical gradient calculation requires interface function external_G')
elif self.gradient_method is 'numerical':
if self.external_G is not None:
g_adiabatic = self.external_G(a_x)
else:
g_adiabatic = self._gradients_num(a_x)
else:
raise ValueError('Method for adiabatic gradient not defined')
# compute energy and gradient differences
self.delta_e2 = self.delta_e1
self.delta_e1 = abs(self.delta_e)
for i in range(self.nstates):
for j in range(self.nstates):
self.delta_e[i,j] = e_adiabatic[i] - e_adiabatic[j]
delta_g = g_adiabatic - g_adiabatic[self.state]
# compute dDE/dt and keep track of the maximum for every state
dEdt = abs(np.dot(delta_g,self.simbox.V))
dEdt[self.state] = 1.0 # avoid zerodivision, this entry is useless
for i,d in enumerate(dEdt):
if d > self.dEdtmax[i]: self.dEdtmax[i] = d
delta_e_s = self.delta_e[:,self.state]
delta_e1_s = self.delta_e1[:,self.state]
delta_e2_s = self.delta_e2[:,self.state]
# keep track of minimum (and maximum) energy differences
indxmin = np.zeros(self.nstates,int)
minfound = False
for i,e in enumerate(abs(delta_e_s)):
if delta_e1_s[i] < delta_e2_s[i] and delta_e1_s[i] < e: # loc. min
minfound = True
indxmin[i] = 1
if delta_e1_s[i] > delta_e2_s[i] and delta_e1_s[i] > e: # loc. max
self.dEdtmax[i] = 0.0
indxmin[self.state] = 0 # avoid selfjumps!
# if local minima in DE found, compute hoping prob.
if minfound and not self.lockstate:
minDE = HUGE
for i,indx in enumerate(indxmin):
if indx == 1:
if abs(delta_e_s[i]) < minDE:
minDE = abs(delta_e_s[i])
ipmax = i
if abs(delta_e_s[ipmax]) < self.egap:
# use LZ formula for hopping probability
ip_de = abs(delta_e_s[ipmax])
prob_LZ = np.exp(-pid2*ip_de*ip_de/self.dEdtmax[ipmax])
randnum = np.random.random()
self._logstep['prob_LZ'] = prob_LZ
self._logstep['randnum'] = randnum
self._logstep['ipmax'] = ipmax
self._logstep['dEdtmax'] = self.dEdtmax[ipmax]
self.dEdtmax[ipmax] = 0.0
if prob_LZ > randnum:
# save local information to be used from the trajectory object
self.is_hopp = True
self.hop_state = ipmax
self.de_s = delta_e_s[ipmax]
self.dg_s = delta_g[ipmax]
self.gradient_s = g_adiabatic[ipmax]
f = self.simbox.logfile
f.write('## Hopping probabilities\n')
f.write(str(self._logstep['prob_LZ'])+'\n')
f.write('## ipmax prob rand dEdtmax\n')
f.write(' %i %8.6f %8.6f %8.6f\n' % (
self._logstep['ipmax'],
self._logstep['prob_LZ'],
self._logstep['randnum'],
self._logstep['dEdtmax']))
self._logstep['delta_e'] = delta_e_s
self._logstep['delta_g'] = delta_g
self._logstep['dEdt'] = dEdt
if self._DEBUG:
sys.stdout.write('***DEBUG gradients ***\n')
for i,gradient_s in enumerate(g_adiabatic):
sys.stdout.write('state '+str(i)+'\n')
for g in gradient_s:
sys.stdout.write(str(g)+' ')
sys.stdout.write('\n')
sys.stdout.write('***DEBUG gradients ***\n')
return e_adiabatic[self.state],g_adiabatic[self.state]
# xnew,vnew,is_hopp = self.switch_condition(delta_e_s[ipmax],g_adiabatic[ipmax])
[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
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 _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')
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']))