#* **************************************************************************
#*
#* 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 sys
import numpy as np
TINY = 1.e-8
[docs]
class wrapperInterface(object):
def __init__(self,**opts):
"""
This class wraps around specific quantum chemistry interfaces and
exposes functions with exactly the argument structure
iexpected by dynamics tools.
It is also responsible for achieving an optimal use of the quantum
chemistry programs in terms of the number of calls. This may involve
caching of the energies, gradients, etc. of the latest geometry or even
further back in time if necessary.
In the future interpolation optimizations might be incorporated at this
level.
"""
self.QCE = None # Interface to a quantum chemistry program
self._oldX = None
self._oldX2 = None
self._oldS = -1
self.forces = 'analytical'
self._oldDnac = None
self._oldDmu = None
self._stateOverlap = None
#self.calcOverlaps = opts.get('calcOverlaps',False)
[docs]
def f_E(self, x):
"""
Return energies using a ground state method.
X -- 3N array with atomic coordinates
"""
return self.f_EGrad_GS(x)[0]
[docs]
def f_EGrad_GS(self,X,S=0,t=0):
"""
Return energies and gradients using a ground state method
X -- 3N array with atomic coordinates
S -- electronic state. This is disregarded here
t -- time
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
return self._oldE,self._oldG
else:
# if self.forces == 'numerical': raise ValueError('some message')
e,g = self.QCE.energy_gradient_HF(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
return e,g
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X):
return self._oldE,self._oldG
else:
e, g = self.QCE.energy_gradient_HF(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
return e, g
elif self.QCE.engine == 'embedding':
if self._is_old_geometry(X):
return self._oldE,self._oldG
else:
e, g = self.QCE.energy_gradient(X,Time=t)
self._oldX = X.copy()
self._oldE = e[0]
self._oldG = g[0]
return e[0], g[0]
[docs]
def f_EGrad_ES(self,X,S,t=0.0):
"""
Return energies and gradients for excited states
Used in single state (Born-Oppenheimer) calculations
X -- 3N array with atomic coordinates
S -- electronic state
t -- time
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
return self._oldE[S],self._oldG[S]
else:
e,g = self.QCE.mcscf(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
return e[S], g[S]
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X) and S == self._oldS:
return self._oldE[S],self._oldG[S]
else:
if self.QCE.cis:
if self.QCE.calcOverlaps:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X,Time=t,state=S, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac = self.QCE.CIS_properties(X,Time=t,state=S)
stateOverlap = None
elif self.QCE.part:
occ, e, g, nac = self.QCE.part_properties(X,Time=t)
self.QCE.set_occ(occ[S])
elif self.QCE.hole:
occ, e, g, nac = self.QCE.hole_properties(X,Time=t)
self.QCE.set_occ(occ[S])
else:
raise ValueError("Please specify either gs_occ = yes or gs_occ_part = yes in xmolecule section")
if self._oldX is not None:
self._oldX2 = X2.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldS = S
return e[S], g[S]
[docs]
def f_overlap(self, X1, X2, S):
"""
Returns the overlap between electronic states at geometry X1 and X2
"""
if self.QCE.engine == 'XMolecule':
if self.QCE.cis:
assert(self.QCE.calcOverlaps)
if self._is_old_geometry(X1, oldX=self._oldX2) and self._is_old_geometry(X2, oldX=self._oldX):
return self._stateOverlap
elif self._is_old_geometry(X1, oldX=self._oldX):
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X2, S=S, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X1, S=S, overlap=self.QCE.calcOverlaps)
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X2, S=S, overlap=self.QCE.calcOverlaps)
self._stateOverlap = stateOverlap
self._oldX2 = X1.copy()
self._oldX = X2.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldS = S
return stateOverlap
elif self.QCE.engine == 'embedding':
return self.QCE.f_overlap(X1, X2, S)
print("f_overlap not implemented")
return(np.eye(self.nstates))
[docs]
def f_EGrad_NA(self,X,V,t,S, full=False):
"""
Return energies and gradients for excited states
Used in non-adiabtic calculations
X -- 3N array with atomic coordinates
V -- 3N array with atomic velocities
t -- time
S -- electronic state
full -- optional; bool. Return all energies and gradients; only implemented for XMolecule.
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
return self._oldE[S],self._oldG[S]
else:
e,g,nac,mu = self.QCE.nacme(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldMU = mu
return e[S], g[S]
elif self.QCE.engine == 'molcas':
if self._is_old_geometry(X) and S == self._oldS:
return self._oldE[S], self._oldG
else:
if self.forces == 'analytical':
e,gs,D = self.QCE.rasci_D(X,V,Time=t,root_g=S+1)
elif self.forces == 'numerical':
e,gs,D = self.QCE.rasci_D_num(X,V,Time=t,root_g=S+1)
elif self.forces == 'compare':
ea,gsa,Da = self.QCE.rasci_D(X,V,Time=t,root_g=S+1)
en,gsn,Dn = self.QCE.rasci_D_num(X,V,Time=t,root_g=S+1)
sys.stdout.write('g_analytical\n')
# sys.stdout.write(gsa)
print(gsa)
sys.stdout.write('\n')
sys.stdout.write('g_numerical\n')
# sys.stdout.write(gsn)
print(gsn)
sys.stdout.write('\n')
sys.exit(0)
self._oldX = X.copy()
self._oldE = e
self._oldG = gs
self._oldD = D
return e[S],gs
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X) and S == self._oldS:
return self._oldE[S], self._oldG[S]
else:
if self.QCE.cis:
if self.QCE.calcOverlaps:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X,Time=t,state=S, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac = self.QCE.CIS_properties(X,Time=t,state=S)
stateOverlap = None
elif self.QCE.part:
stateOverlap = None
occ, e, g, nac = self.QCE.part_properties(X,Time=t)
self.QCE.set_occ(occ[S, :].tolist())
elif self.QCE.hole:
stateOverlap = None
occ, e, g, nac = self.QCE.hole_properties(X,Time=t)
self.QCE.set_occ(occ[S, :].tolist())
else:
raise ValueError("Please specify either gs_occ = yes or gs_occ_part = yes in xmolecule section")
if self._oldX is not None:
self._oldX2 = self._oldX.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldS = S
if full == True:
return e, g
return e[S], g[S]
elif self.QCE.engine == 'embedding':
if self._is_old_geometry(X):
return self._oldE[S], self._oldG[S]
else:
e, g = self.QCE.energy_gradient(X, V=V, S=S, Time=t, SH=True)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
return e[S], g[S]
[docs]
def f_NAC(self,X,V,t,SI,SJ):
"""
Return NAC vector between two states
X -- 3N array with atomic coordinates
V -- 3N array with atomic velocities
t -- time
SI -- Ith electronic state
SJ -- Jth electronic state
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
return self._oldNAC[(SI,SJ)]
else:
e,g,nac,mu = self.QCE.nacme(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldMU = mu
return nac[(SI,SJ)]
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X):
return self._oldNAC[SI, SJ]
else:
if self.QCE.cis:
if self.QCE.calcOverlaps:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X,Time=t, S=SI, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac = self.QCE.CIS_properties(X,Time=t, S=SI)
stateOverlap = None
elif self.QCE.part:
stateOverlap = None
occ, e, g, nac = self.QCE.part_properties(X,Time=t)
elif self.QCE.hole:
stateOverlap = None
occ, e, g, nac = self.QCE.hole_properties(X,Time=t)
else:
raise ValueError("Please specify either gs_occ = yes or gs_occ_part = yes in xmolecule section")
if self._oldX is not None:
self._oldX2 = self._oldX.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldS = SI
return nac[SI, SJ]
elif self.QCE.engine == 'embedding':
# if self._is_old_geometry(X):
# return self._oldNAC
# else:
nac = self.QCE.nac(X, SI, SJ, V=V, t=t)
self._oldX = X.copy()
self._oldNAC = nac
return nac
[docs]
def f_W(self,X,V,t):
"""
Return diagonal matrix with adiabatic energies
X -- 3N array with atomic coordinates
t -- time
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
return np.diag(self._oldE)
else:
e,g,nac,mu = self.QCE.nacme(X,Time=t)
if self._oldX is not None:
self._oldX2 = self._oldX.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldNAC = nac
self._oldMU = mu
self._oldS = SI # fix this?
return np.diag(e)
elif self.QCE.engine == 'molcas':
if self._is_old_geometry(X):
return np.diag(self._oldE)
else:
raise ValueError('f_W')
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X):
return np.diag(self._oldE)
else:
if self.QCE.cis:
if self.QCE.calcOverlaps:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X, Time=t, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac = self.QCE.CIS_properties(X, Time=t)
stateOverlap = None
self._oldS = 0 # if S is omitted in CIS_properties, we return gradients for the ground state
elif self.QCE.part:
stateOverlap = None
occ, e, g, nac = self.QCE.part_properties(X,Time=t)
elif self.QCE.hole:
stateOverlap = None
occ, e, g, nac = self.QCE.hole_properties(X,Time=t)
else:
raise ValueError("Please specify either gs_occ = yes or gs_occ_part = yes in xmolecule section")
if self._oldX is not None:
self._oldX2 = self._oldX.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldNAC = nac
self._oldG = g
return np.diag(e)
elif self.QCE.engine == 'embedding':
return self.QCE.getW(X, V=V, t=t)
[docs]
def f_D(self,X,V,t):
"""
Return matrix with kinetic couplings
X -- 3N array with atomic coordinates
V -- 3N array with atomic velocities
t -- time
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
D = self._kineticCouplingsMatrix(self._oldNAC,V,t,len(self._oldE))
return D
else:
e,g,nac,mu = self.QCE.nacme(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldMU = mu
D = self._kineticCouplingsMatrix(nac,V,t,len(e))
self._oldDnac = D
return D
elif self.QCE.engine == 'molcas':
if self._is_old_geometry(X):
return -1.0j*self._oldD
else:
raise ValueError('f_D')
elif self.QCE.engine == 'XMolecule':
if self._is_old_geometry(X):
D = self._kineticCouplingsMatrix(self._oldNAC,V,t,len(self._oldE))
return D
else:
if self.QCE.cis:
if self.QCE.calcOverlaps:
l, e, g, nac, stateOverlap = self.QCE.CIS_properties(X,Time=t, overlap=self.QCE.calcOverlaps)
else:
l, e, g, nac = self.QCE.CIS_properties(X,Time=t)
stateOverlap = None
self._oldS = 0 # if S is omitted in CIS_properties, we return gradients for the ground state
elif self.QCE.part:
stateOverlap = None
occ, e, g, nac = self.QCE.part_properties(X,Time=t)
elif self.QCE.hole:
stateOverlap = None
occ, e, g, nac = self.QCE.hole_properties(X,Time=t)
else:
raise ValueError("Please specify either gs_occ = yes or gs_occ_part = yes in xmolecule section")
if self._oldX is not None:
self._oldX2 = self._oldX.copy()
self._oldX = X.copy()
self._stateOverlap = stateOverlap
self._oldE = e
self._oldNAC = nac
self._oldG = g
D = self._kineticCouplingsMatrix(nac,V,t,len(e))
self._oldDnac = D
return D
elif self.QCE.engine == 'embedding':
D = self.QCE.getD(X, V, t=t)
return D
[docs]
def f_Dmu(self,X,V,t):
"""
Return matrix with kinetic couplings generated by an external E field
X -- 3N array with atomic coordinates
V -- 3N array with atomic velocities
t -- time
"""
if self.QCE.engine == 'GamessUS':
if self._is_old_geometry(X):
D = self._kineticCouplingsMatrixMu(self._oldMU,V,t,len(self._oldE))
return D
else:
e,g,nac,mu = self.QCE.nacme(X,Time=t)
self._oldX = X.copy()
self._oldE = e
self._oldG = g
self._oldNAC = nac
self._oldMU = mu
D = self._kineticCouplingsMatrixMu(mu,V,t,len(e))
self._oldDmu = D
return D
elif self.QCE.engine == 'molcas':
raise ValueError('f_Dmu not implemented for the molcas interface')
[docs]
def f_MOE(self):
"""
Return array with molecular orbital energies of occupied orbitals
"""
if self.QCE.engine == 'XMolecule':
occ = np.asarray(self.QCE.get_occ())
moe = np.asarray(self.QCE.getOrbitalEnergies())
return moe[occ>0]
else:
raise ValueError('f_MOE not implemented for this Interface.')
def _is_old_geometry(self, X, oldX=None):
"""
This function checks if a new geometry is the same as the already stored geometry from the
previous step. Please note that similarity is meassured up to a norm difference of
TINY = 1.e-8.
Input:
X -- 3N array with geometry to be checked.
Output:
Returns False if given geometry is different or no old geometry is saved.
Returns True if given geometry is the same (within 1e-8) as old geometry.
"""
if oldX is None:
oldX = self._oldX
if oldX is None:
return False
else:
n = np.linalg.norm( X - oldX )
if n > TINY:
return False
else:
return True
def _kineticCouplingsMatrix(self,nac,v,t,ns):
"""
This function returns the coupling matrix due to nuclear motion. This matrix is used
in the fewest switches algorithm to propagate the nuclear Schroedinger equation and
to calculate the hopping probabilities. Please note that this matrix includes a factor
of -i. Also the phases of the NAC is tracked and corrected.
Input:
nac -- Non-adiabatic couplings obtained from electronic structure calculations
v -- 3N vector containing the current velocity
t -- Current time
ns -- Number of states
Output:
D -- -i * (NAC * V); Matrix containing the couplings due to nuclear motion.
"""
D = np.zeros((ns,ns),complex)
for i in range(nac.shape[0]):
for l in range(nac.shape[1]):
D[i,l] = -1.0j * np.dot(nac[i, l, :],v )
for i in range(ns):
for l in range(i):
if( abs(D[i,l] - D[l,i].conjugate()) > 1e-6):
raise ValueError('D is not hermitian: {:s} !'.format(D))
if self._oldDnac == None:
return D
else:
p = D*self._oldDnac
mask = np.piecewise(p, [p.imag < 0, p.imag >= 0], [-1, 1])
return D*mask
def _kineticCouplingsMatrixMu(self,mu,v,t,ns):
"""
This function returns the coupling matrix due to an electric field explicitly included in the
electronic structure calculations. This matrix is used in the fewest switches algorithm to
propagate the nuclear Schroedinger equation and to calculate the hopping probabilities.
Please note that this matrix includes a factor of -i. Also track the phases of the coupling.
Input:
mu -- Transition-dipole matrix elements
v -- 3N vector containing the current velocity
t -- Current time
ns -- Number of states
Output:
D -- i * (MU * (dt eps)/Delta E); Matrix containing the couplings due to the electric field explicitly included in the
electronic structure calculations
"""
D = np.zeros((ns,ns),complex)
for i in range(mu.shape[0]):
for j in range(mu.shape[1]):
if i == j: continue
eji = self._oldE[j] - self._oldE[i]
D[i,j] = -1.0j * np.dot(mu[i,j],self.QCE.dt_efield(t) )/(eji + TINY)
if self._oldDmu == None:
return D
else:
p = D*self._oldDmu
mask = np.piecewise(p, [p.imag < 0, p.imag >= 0], [-1, 1])
return D*mask