#* **************************************************************************
#*
#* 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 numpy as np
import CDTK.Tools.Conversion as cv
import CDTK.Tools.WaterModel as wm
[docs]
class MolecularMechanics(object):
"""
...
"""
def __init__(self, **opts):
"""
Initialize the MM region
...
"""
self.natoms_qm = None # Atoms in QM region
# Actual atom labels, positions and velocities
# MM region before point charge transformation
self.R = []
self.V = []
self.atomlist = []
self.atommasses = []
self.atommasses3 = []
# XMOLECULE input form
# Point charge values from water model
self.x_R = []
self.x_V = []
self.x_C = [] # point charges
self.x_M = [] # corresponding masses may be fractional
self.x_M3 = []
self.x_L = [] # point charge labels
# Gradients between regions
self.grad_qm_mm = [] # QM grad comes from QM part
self.grad_mm_mm = []
self.grad_mm_qm = []
# Coordinates of the QM region
self.qm_R = []
self.qm_V = []
self.qm_atomlist = []
self.qm_M = []
# Integration steps need also gradient from last time steps?
# MM values
self.dt = opts.get('dt', 0.1*cv.fs2au) # MM time step
self.mode = opts.get('mode', 'electron') # Water model or electron
self.t = opts.get('t', 0.0) # may have its own time if different dt
def _step(self):
"""
Move the whole MM region for one time step
This step may be different e.g. for the classical electron
"""
pass
[docs]
def update_grad_mm_mm(self):
"""
Gradient MM region -> MM region
"""
pass
[docs]
def update_grad_mm_qm(self):
"""
Gradient MM region -> QM region
"""
self.grad_mm_qm = np.zeros(len(self.qm_R))
if self.mode == 'SPC':
# Lennart-Jones potential
for i,l_qm in enumerate(self.qm_atomlist):
if l_qm=='O':
for j,l_mm in enumerate(self.atomlist):
if l_mm=='O':
r_OO = self.qm_R[3*i:3*i+3] - self.R[3*j:3*j+3]
self.grad_mm_qm[3*i:3*i+3] = wm.lennard_jones_gradient( r_OO, mode='SPC' )
else:
print('MM to QM gradient not yet implemented for ' + self.mode)
pass
[docs]
def mm_to_point_charge(self,qce='xmolecule',mode='electron'):
"""
Returns point charge positions and charge values from atomic input in MM region
Modes:
--------------------
electron - classical electron as point charge
SPC - single point charge water model
TIP4P/2005f TODO
"""
# TODO set mode here or in class
if qce=='xmolecule':
if mode=='electron':
if len(self.atomlist) == 1 and self.atomlist[0]=='e':
self.x_R = self.R
self.x_C = np.array([-1.0])
self.x_V = self.V
self.x_M = self.atommasses
self.x_M3 = self.atommasses3
self.x_L = ['X']
else:
raise ValueError('Can only have one classical electron')
elif mode=='SPC':
self.x_R = self.R
self.x_V = self.V
self.x_M = np.tile(1e3, len(self.atommasses))
self.x_M3 = np.tile(1e3, 3*len(self.atommasses))
self.x_L = np.tile('X', len(self.atommasses))
# point charges from SPC
tmp_C = []
for lab in self.atomlist:
if lab=='O':
tmp_C.append( wm.waterModel[mode]['q_O'] )
elif lab=='H':
tmp_C.append( wm.waterModel[mode]['q_H'] )
else:
raise ValueError('Unknown atom label for SPC: ' + lab)
self.x_C = np.asarray(tmp_C)
elif mode=='TIP4P/2005f':
raise ValueError( mode + ' not yet implemented')
else:
raise ValueError( mode + ' not valid')
else:
raise ValueError('mm_to_point charge only implemented for QCE XMolecule')
return self.x_R, self.x_C