# * **************************************************************************
# *
# * 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 numpy as np
from CDTK.Tools.Utils import getUniqueID
import CDTK.Tools.Utils as uti
import CDTK.Tools.EField as efi
import CDTK.Tools.Conversion as conv
import CDTK.Interfaces.ForceField as ff
import CDTK.Interfaces.WrapperInterface as win
TINY = 1.e-8
# - Default definitions
# !! alteration might lead to unittest failures !!
[docs]
class embedding(object):
"""
This defines a class to embedd a QM interface into a MM region.
"""
def __init__(self, atomlist, atomnums, atommass3, xcart, **opts):
"""
Initialize and set parameters.
TODO: Say here something about the possible input parameters
"""
# todo: Define here flag for electrostatic embedding calculation and its initialization
self.engine = 'embedding' # Used in the wrapper class
self._counter = 0 # Execution counter
self.keepdir = opts.get('keepdir', None)
# -------------------------------------------------
self.atomicSymbols = opts.get('atomicSymbols', [])
self.atomicNumbers = opts.get('atomicNumbers', [])
self.atomicMasses = opts.get('atomicMasses', [])
self.efields = opts.get('efields', [])
# -------------------------------------------------
self.quantum = opts.get('quantum', ['xmolecule'])[0]
self.qregion = int(opts.get('region', ['3'])[0])
self.noMregion = 'nomregion' in opts
self.forcefield = opts.get('forcefield', [None])[0]
self.electrostatic = bool(opts.get('electrostatic', None)) # Flag for electrostatic embedding calculation.
self.sigma_case = opts.get('sigma_case', None)
if self.sigma_case is not None and len(self.sigma_case) == 1: # ! there should be a better way to do this...
self.sigma_case = self.sigma_case[0]
self.gndm = bool(self.sigma_case) # Flag for Gaussian Nuclear Distribution Model.
self.qmmmRepulsionOH = opts.get('repoh', False) # extra repulsion interaction between QM hydrogens and MM oxygens
self.molecule2ffstr = None
if opts.get('molecule2ff', None) is not None:
self.molecule2ffstr = ''.join(opts.get('molecule2ff', None))
self.quantum_options = opts.get('quantum_options', None)
self.calcOverlaps = opts.get('calcOverlaps', False)
if 'topol' in opts:
topol = opts['topol']
elif os.path.exists('topol'):
topol = uti.readColumnFile('topol')
else:
raise ValueError('no file ./topol or "topol" section in input file found')
self.topol = np.array(topol)
# Define slices for the postions of QM region for FF or QCE.
self.qmXslice = slice(0, 3 * self.qregion)
self.mmXslice = slice(3 * self.qregion, None)
if self.electrostatic:
self.qceXslice = slice(0, None)
else:
self.qceXslice = slice(0, 3 * self.qregion)
# --------------------------------------------------------------------------
# Create interface to quantum chemistry engine (QCE)
# Currently on xmolecule is allowed - others need to be implemented
# --------------------------------------------------------------------------
is_keepinp = False
is_keeplog = False
# if I['system'].has_key('keepinp'): is_keepinp = True
# if I['system'].has_key('keeplog'): is_keeplog = True
print("QM region %d atoms" % self.topol[:self.qregion].shape[0])
print("MM only region %d atoms" % self.topol[self.qregion:].shape[0])
self.MMFull = ff.jointForceFields(
self.topol,
molecule2ffstr = self.molecule2ffstr,
waterForcefield=self.forcefield,
sigma_case=self.sigma_case,
coulomb = not self.electrostatic)
self.MMQM = ff.jointForceFields(
self.topol[:self.qregion],
waterForcefield=self.forcefield,
molecule2ffstr = self.molecule2ffstr,
coulomb = not self.electrostatic)
if self.electrostatic:
self.MMonly = ff.jointForceFields(
self.topol[self.qregion:],
waterForcefield=self.forcefield,
molecule2ffstr = self.molecule2ffstr,
bonded=False, lennard_jones=False)
if self.qmmmRepulsionOH:
self.qmmmRepPairs = []
for m1 in self.MMonly.molecules:
for m2 in self.MMQM.molecules:
for j1,t1 in enumerate(self.MMonly.topolMs[m1]):
for j2,t2 in enumerate(self.MMQM.topolMs[m2]):
if t1.split('_')[1] == 'OW' and t2.split('_')[1].startswith("HW"):
self.qmmmRepPairs.append((self.MMFull.atomsOfMs[m1][j1],
self.MMFull.atomsOfMs[m2][j2]))
print(f"adding extra repulsion between MM_{t1}={self.MMFull.atomsOfMs[m1][j1]} QM_{t2}={self.MMFull.atomsOfMs[m2][j2]}")
self.qmmmRepPairs = np.array(self.qmmmRepPairs)
if self.qmmmRepPairs.shape[0] > 0:
self._qmmmRepPairInd1 = [
np.where(self.qmmmRepPairs[:, 0] == i)[0]
for i in range(self.topol.shape[0])]
self._qmmmRepPairInd2 = [
np.where(self.qmmmRepPairs[:, 1] == i)[0]
for i in range(self.topol.shape[0])]
#if self.qregion == 0:
# return
# for XMolecule as a quantum chemistry tool
if self.quantum == 'xmolecule':
import CDTK.Interfaces.XMoleculeInterface as xmim
QCE = xmim.xmolecule()
QCE.is_storeinp = is_keepinp # Keep inp file of every calculation performed
QCE.is_storelog = is_keeplog # Keep log file of every calculation performed
QCE.atomicSymbols = uti.changeIsotopeSymbols(atomlist[:self.qregion])
QCE.atomicNumbers = atomnums[:self.qregion] # ? Is this doing anything
QCE.atomicMasses = np.array(atommass3[:self.qregion*3], float)*conv.am2au
QCE.calcOverlaps = self.calcOverlaps
if self.quantum_options:
if not self.electrostatic and not self.gndm:
x_pos = None
x_charges = None
x_sigma = None
elif self.electrostatic and not self.gndm:
self.quantum_options['qmmm'] = ['yes']
x_pos = xcart[self.qregion*3:]
x_charges = self.MMFull.charges[self.qregion:]
x_sigma = None
elif self.electrostatic and self.gndm:
self.quantum_options['qmmm'] = ['yes']
x_pos = xcart[self.qregion*3:]
x_charges = self.MMFull.charges[self.qregion:]
x_sigma = self.MMFull.sigma_nuc[self.qregion:]
elif not self.electrostatic and self.gndm:
raise ValueError("using nuclear spreading (sigma_case) makes only sense with electrostatic=yes")
mm_charges = {'x_pos': x_pos, 'x_charges': x_charges, 'x_sigma': x_sigma}
QCE.init_input_options(self.quantum_options,
atomlist[:self.qregion], xcart[:self.qregion*3],
**mm_charges)
else:
print("field $quantum is missing in the inputfile")
exit(-1)
else:
"QCE required not implemented in embedding: ", self.quantum
exit(-1)
self.W = win.wrapperInterface()
self.W.QCE = QCE
# here wrapper forces has to be set somehow
#self.W.forces = self.forces
#W.forces = forces[0]
#else:
# W.forces = 'analytical'
return
[docs]
def f_overlap(self, X1, X2, S):
self.W.QCE.calcOverlaps = self.calcOverlaps
return self.W.f_overlap(X1[self.qceXslice], X2[self.qceXslice], S)
[docs]
def hasConstraints(self):
"""
Checks wether there are constraints in the parts of the geometry modelled
by forcefields.
"""
if self.MMFull._numConstraints > 0:
return True
return False
[docs]
def enforceConstraints(self, x):
return self.MMFull.enforceConstraints(x)
[docs]
def getConstraintsCorrection(self, vec0, vec1, dt, case):
"""
Calculates the correction to the position or velocities
of atoms in case there are constraints. For now, only
RATTLE algorithm implemented.
Parameters
----------
vec0 : (3 * n_atoms) ndarray of float.
Array with position or velocities of atoms
vec1 : (3 * n_atoms) ndarray of float.
Array with position or velocities of atoms
dt : float
Timestep
case : str
Flag for different types of correction
Returns
-------
correction: (3 * n_atoms) ndarray of float.
Array with correction for position or velocities.
"""
correction = np.zeros(vec0.shape, dtype=float)
correction[3*self.qregion:] = self.MMFull.getRattleCorrection(vec0, vec1, dt, case)[3*self.qregion:]
return correction
[docs]
def getPartialCharges(self, chargetype='mulliken', state=0):
""" returns partial charges of all atoms
Args:
chargetype (str, optional): Type. Defaults to 'mulliken'.
state (int, optional): state index. Defaults to 0.
"""
charges = self.MMFull.charges[:].copy()
charges[:self.qregion] = self.W.QCE.getPartialCharges(chargetype=chargetype, state=state)[:self.qregion]
return charges
[docs]
def energy_gradient(self, X, V=None, S=0, Time=0.0, SH=False, **opts):
"""
Calculates the energy and gradient of the system using an ONIOM scheme:
:math: `E = E_{MM}(MM+QM) + E_{QM}(QM) - E_{MM}(QM)`
For the electrostatic embedding calculation, the Coulomb interaction in
the QM region and the QM-MM Coulomb interaction are calculated with QCE.
In this case, only the MM Coulomb energy is calculated with the forcefield.
Parameters
----------
X : (n_atoms, 3) ndarray of floats.
Atomic positions, in Bohr (au).
V : (n_atoms, 3) ndarray of floats., optional
Atomic velocities in au, by default None
S : int, optional
Current electronic state, by default 0
Time : float, optional
Current propagation time in au, by default 0.0
SH : bool, optional
Whether to run surface hopping or not for QM energy, by default False
Returns
-------
e: float
Energy in atomic units
g: (n_atoms, 3) ndarray of floats.
Gradient in atomic units.
"""
self.W.QCE.calcOverlaps = self.calcOverlaps
# ! REVISE: Virtual site exceptions for electrostatic embedding
# do not use virtual sites within quantum region
if self.noMregion == True:
listNoM = range(self.qregion)
else:
listNoM = []
# get MM energy for full system
e_mm_full, g_mm_full = self.MMFull.getEnGrad(X, listNoM)
if self.qregion == 0:
return [e_mm_full], [g_mm_full]
# get MM energy for quantum region
e_mm_region, g_mm_region = self.MMQM.getEnGrad(X[self.qmXslice], listNoM)
if self.electrostatic and self.MMFull.hasVirtualSites:
tmpX = X.copy()
X, vXderiv = self.MMFull.geomWithVirtualSites(X.reshape((-1, 3)), listNoM) # todo: better name for vXderiv
X = X.flatten()
# get qm energy for quantum region
if SH == False:
e_qm_region, g_qm_region = self.W.f_EGrad_GS(X[self.qceXslice])
else:
e_qm_region, g_qm_region = self.W.f_EGrad_NA(
X[self.qceXslice], V[self.qceXslice], Time, S, full=True)
# obtain full energy for every state
e = e_mm_full + e_qm_region - e_mm_region
# obtain number of states; correct for ground state
try:
e = np.array([e]).flatten()
n_states = len(e)
except TypeError:
n_states = 1
# translate gradient on virtual sites to gradient on real atoms.
# todo: Weird reshaping of g_qm_region, FIX THIS OR MAKE IT MORE CLEAR
if self.electrostatic and self.MMFull.hasVirtualSites:
X = tmpX
g_qm_region = np.reshape(g_qm_region, (n_states, -1)).copy() # for ground state calculations
for s in range(g_qm_region.shape[0]): # Change for n_states
g_qm_region[s, :] = self.MMFull.gradient_virt2orig(g_qm_region[s, :].reshape(-1, 3), vXderiv).flatten()
# repeat full MM gradient for every state
g = np.tile(g_mm_full, (n_states, 1)) # np.repeat(g_mm_full, n_states).reshape(-1, n_states).T
# ! REVISE: Slicing
# substract MM gradient within quantum region for every state
g[:, self.qmXslice] -= np.tile(g_mm_region, (n_states, 1))
# add gradients calculated with QCE to the relevant region (QM or QM+MM) for every state.
g[:, self.qceXslice] += g_qm_region
eCoulMM = 0
if self.electrostatic: # add Coulomb interaction in the MM region
eCoulMM, gCoulMM = self.MMonly.getEnGrad(X[self.mmXslice], listNoM)
e += eCoulMM
g[:, self.mmXslice] += np.tile(gCoulMM, (n_states, 1))
eQMMMRep = 0
if self.qmmmRepulsionOH and self.qmmmRepPairs.shape[0] > 0:
geom = X.reshape((-1, 3))
dist = geom[self.qmmmRepPairs[:, 0],:] - geom[self.qmmmRepPairs[:, 1],:]
d = np.linalg.norm(dist, axis=1)
distOverD = dist * (1/d)[:, np.newaxis]
eP, dEdD = self.qmmmRepEgrad(d)
eQMMMRep += eP
e += eP
dEdR = distOverD[:] * dEdD[:, np.newaxis]
for a1 in range(self.topol.shape[0]):
if self._qmmmRepPairInd1[a1].shape[0] > 0:
tmp = (dEdR[self._qmmmRepPairInd1[a1], :]).sum(axis=0)
g[:, 3 * a1: 3 * (a1 + 1)] += np.tile(tmp, (n_states, 1))
if self._qmmmRepPairInd2[a1].shape[0] > 0:
tmp = (dEdR[self._qmmmRepPairInd2[a1], :]).sum(axis=0)
g[:, 3 * a1: 3 * (a1 + 1)] += - np.tile(tmp, (n_states, 1))
return e, g
[docs]
def qmmmRepEgrad(self, d, alpha = 0.5*1/(0.4*conv.an2au)**2, epsilon = 500.0 * conv.kc2au):
"""computes Pauli exchance repulsion
E = epsilon * exp(-alpha d**2)
Args:
d (double array shape(n)): distances
Returns:
double: energy
np.array of doubles: derivatives of energy with respect to distances
"""
#alpha = 0.15 * 1./(conv.an2au)**2
#epsilon = 20.0 * conv.kc2au
d2 = d*d
eTerm = epsilon * np.exp(-alpha * d2)
dEdD = - (2. * alpha * d ) * eTerm
return eTerm.sum(), dEdD
[docs]
def nac(self, X, SI, SJ, V, t=0.0, **opts):
"""
Nonadiabatic coupling vector between states SI and SJ.
Parameters
----------
X : (n_atoms, 3) ndarray of floats.
Atomic positions, in Bohr (au).
SI : int
Current electronic state index.
SJ : int
Index of state to which the coupling is calculated.
V : (n_atoms, 3) ndarray of floats.
Atomic velocities in au.
t : float, optional
Current propagation time in au, by default 0.0
Returns
-------
coupling: (n_atoms) ndarray of complex.
Nonadiabatic coupling vector.
"""
self.W.QCE.calcOverlaps = self.calcOverlaps
if self.electrostatic and self.MMFull.hasVirtualSites:
# ! REVISE: Virtual site exceptions for electrostatic embedding
# do not use virtual sites within quantum region
if self.noMregion == True:
listNoM = range(self.qregion)
else:
listNoM = []
tmpX = X.copy()
X, vXderiv = self.MMFull.geomWithVirtualSites(X.reshape((-1, 3)), listNoM)
X = X.flatten()
# ! REVISE: Check slicing of the result
coupling = np.zeros(len(X))
coupling[self.qmXslice] = self.W.f_NAC(X[self.qceXslice], V[self.qceXslice], t, SI, SJ)[self.qmXslice]
if self.electrostatic and self.MMFull.hasVirtualSites:
coupling = self.MMFull.gradient_virt2orig(coupling.reshape(-1, 3), vXderiv).flatten()
return coupling
[docs]
def getW(self, X, V=None, t=0.0, **opts):
"""
Diagonal matrix with all adiabatic energies.
Parameters
----------
X : (n_atoms, 3) ndarray of floats.
Atomic positions, in Bohr (au).
V : (n_atoms, 3) ndarray of floats., optional
Atomic velocities in au, by default None
t : float, optional
Current propagation time in au, by default 0.0
Returns
-------
e: __type__
Diagonal matrix with adiabatic energies in atomic units.
"""
# do not use virtual sites within quantum region
if self.noMregion == True:
listNoM = range(self.qregion)
else:
listNoM = []
self.W.QCE.calcOverlaps = self.calcOverlaps
# get MM energy for full system
e_mm_full, _ = self.MMFull.getEnGrad(X, listNoM)
# get MM energy for quantum region
e_mm_region, _ = self.MMQM.getEnGrad(X[self.qmXslice], listNoM)
# calculate total ONIOM MM energy
e_mm = e_mm_full - e_mm_region
if self.electrostatic: # add Coulomb interaction in the MM region
eCoulMM, _ = self.MMonly.getEnGrad(X[self.mmXslice], listNoM)
e_mm += eCoulMM
if self.electrostatic and self.MMFull.hasVirtualSites:
X, vXderiv = self.MMFull.geomWithVirtualSites(X.reshape((-1, 3)), listNoM)
V = self.MMFull.velWithVirtualSites(V.reshape(-1, 3), vXderiv)
X, V = X.flatten(), V.flatten()
# get qm energies for quantum region
e_qm_region = self.W.f_W(X[self.qceXslice], V[self.qceXslice], t)
# create diagonal matrix with mm energies for all qm states
e = np.diag([e_mm]*len(e_qm_region))
# obtain full adiabatic energies
e += e_qm_region
return e
[docs]
def getD(self, X, V, t=0.0, **opts):
"""
Return matrix with kinetic couplings.
Parameters
----------
X : (n_atoms, 3) ndarray of floats.
Array of atomic positions, in Bohr (au).
V : (n_atoms, 3) ndarray of floats.
Array with atomic velocities in au.
t : float, optional
Current propagation time in au, by default 0.0
Returns
-------
_type_
Matrix with kinetic couplings in au.
"""
self.W.QCE.calcOverlaps = self.calcOverlaps
if self.electrostatic and self.MMFull.hasVirtualSites:
# ! REVISE: Virtual site exceptions for electrostatic embedding
# do not use virtual sites within quantum region
if self.noMregion == True:
listNoM = range(len(self.topol))
else:
listNoM = []
X, vXderiv = self.MMFull.geomWithVirtualSites(X.reshape((-1, 3)), listNoM)
V = self.MMFull.velWithVirtualSites(V.reshape(-1, 3), vXderiv)
X, V = X.flatten(), V.flatten()
return self.W.f_D(X[self.qceXslice], V[self.qceXslice], t)