#* **************************************************************************
#*
#* 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
from CDTK.Tools import NumericalDerivatives as nd
[docs]
def velocityVerlet(x0,v0,g0,m,dt,grad,**opts):
"""
Perform a velocity-Verlet integration step
x(t) -> x(t+dt)
v(t) -> v(t+dt)
Input parameters:
x0 - positions x(t)
v0 - velocities v(t)
g0 - gradient g(t)
m - mass of each coordinate
dt - integration step
grad - function returning the gradient
Options:
constraints - List of constraints of the form
g(x) = 0
where g(x) are functions of all coordinates
...
The RATTLE algorithm is used to enforce the constraints:
H.C.Andersen, J.Comput.Phys. 52, 24-34, (1983)
A modified version is implemented here that uses numerical differentiation of the
constraints. This makes the procedure not restricted to particular types of
constraints (e.g. bonds or bond angles): any holonomic constraint on the system of
point particles can be defined.
tol - tolerance for the constraints enforcement in RATTLE
maxiter - maximum number of iterations in RATTLE step
Output:
x1 - x(t+dt)
v1 - v(t+dt)
g1 - g(t+dt)
"""
x1 = np.zeros(x0.size,float)
v1 = np.zeros(x0.size,float)
g1 = np.zeros(x0.size,float)
constraints = opts.get('constraints',None)
tol = opts.get('tol',1.e-6)
maxiter = opts.get('maxiter',500)
if constraints is None:
v1 = v0 - 0.5*dt*g0/m # v(t) -> v(t+dt/2)
x1 = x0 + v1*dt # x(t) -> x(t+dt)
g1 = grad(x1)
v1 = v1 - 0.5*dt*g1/m # v(t+dt/2) -> v(t+dt)
else: # RATTLE iterations
nc = len(constraints)
xlagr = np.zeros(nc,float)
vlagr = np.zeros(nc,float)
# Positions
x1 = x0 + v0*dt - 0.5*g0/m*dt**2 # x(t) -> x(t+dt)
converged = False
it = 0
gconstraints_x0 = _get_gradient_constraints(constraints,x0)
for i in range(nc):
gconstraints_x0[i] = 0.5*gconstraints_x0[i]/m
while not converged:
if it > maxiter:
raise ValueError('RATTLE(x) did not converge in %i iterations' % (maxiter))
converged = True
for i,cfunc in enumerate(constraints):
xold = x1.copy()
s = cfunc(xold)
xlagr[i] = s/np.dot(_ngrad(cfunc,xold),gconstraints_x0[i])
x1 = xold - xlagr[i] * gconstraints_x0[i]
if abs(s) > tol:
converged = False
it += 1
# Velocities
g1 = grad(x1)
v1 = v0 - 0.5*dt*(g0+g1)/m # v(t+dt/2) -> v(t+dt)
for i,gc in enumerate(gconstraints_x0):
v1 = v1 - xlagr[i]*gc*dt
converged = False
it = 0
gconstraints_x1 = _get_gradient_constraints(constraints,x1)
m_gconstraints_x1 = []
for i in range(nc):
m_gconstraints_x1.append(0.5*gconstraints_x1[i]/m)
u_gconstraints_x1 = map(_unitvector,gconstraints_x1)
while not converged:
if it > maxiter:
raise ValueError('RATTLE(v) did not converge in %i iterations' % (maxiter))
converged = True
for i,cfunc in enumerate(constraints):
vold = v1.copy()
vlagr[i] = np.dot(vold,gconstraints_x1[i])/np.dot(m_gconstraints_x1[i],gconstraints_x1[i])
s = np.dot(vold,gconstraints_x1[i])
v1 = vold - vlagr[i]*m_gconstraints_x1[i]
if abs(s) > tol:
converged = False
it += 1
return x1,v1,g1
def _get_gradient_constraints(constraints,x):
glist = []
for cfunc in constraints:
g = _ngrad(cfunc,x)
glist.append(g)
return glist
def _ngrad(f,x):
return nd.gradient(f,x,fast=True,step=0.001)
def _unitvector(v):
return v/np.linalg.norm(v)