Details on cp2k's velocity verlet integrator?

[Tue Boesen]

I want to make sure I understand exactly how the equations of motions are propagated forward in cp2k, so I decided to code up my own propagator in python and propagate the positions/velocities forward in time using the forces predicted during the simulation. However either I'm doing something wrong or cp2k does not use a simple velocity verlet integrator to propagate a molecular system forward in time.

Does anyone know where I can find some details on this?

My current setup is the following:

_________________________________________________________________________________

import numpy as np
import torch

def atomic_masses(z):
atomic_masses = torch.tensor([0,1.008, 4.0026, 6.94, 9.0122, 10.81, 12.011, 14.007, 15.999, 18.998, 20.180, 22.990, 24.305, 26.982, 28.085])
masses = atomic_masses[z.to(dtype=torch.int64)]
return masses.to(device=z.device)


class VelocityVerletIntegrator(torch.nn.Module):
def __init__(self,dt,force_predictor):
super().__init__()
self.dt = dt
self.fp = force_predictor
return

def r_step(self,r,v,a):
dt = self.dt
r = r + v * dt + 0.5 * a * dt**2
return r

def v_step(self,v,a,a_new):
dt = self.dt
v = v + 0.5 * (a + a_new) * dt
return v

def forward(self,r,v,m,nsteps):
fp = self.fp

F = fp(r)
a = F / m

for i in range(nsteps):
r = self.r_step(r,v,a)
F_new = fp(r)
a_new = F_new / m
v = self.v_step(v,a,a_new)
a = a_new
return r,v


class VelocityVerletIntegrator(torch.nn.Module):
def __init__(self,dt):
super().__init__()
self.dt = dt
return

def r_step(self,r,v,a):
dt = self.dt
r = r + v * dt + 0.5 * a * dt**2
return r

def v_step(self,v,a,a_new):
dt = self.dt
v = v + 0.5 * (a + a_new) * dt
return v

def forward(self,r,v,m,a=None,F=None,F_new=None):
if a is None: # This is used on the first step
a = F / m
r_new = self.r_step(r,v,a)
a_new = F_new / m
v_new = self.v_step(v,a,a_new)
return r_new,v_new,a_new


def MDdataloader(file, position_converter=1,velocity_converter=1,force_converter=1):
data = np.load(file)
R = torch.from_numpy(data['R'])*position_converter
V = torch.from_numpy(data['V'])*velocity_converter
F = torch.from_numpy(data['F'])*force_converter
z = torch.from_numpy(data['z'])
m = atomic_masses(z)
return R,V,F,z,m

if __name__ == '__main__':
torch.set_default_dtype(torch.float64)
bohr = 5.29177208590000E-11 #Bohr radius [m]
hbar = 1.05457162825177E-34 #Planck constant/(2pi) [J*s]
Eh = 4.3597447222071E-18 #Hatree energy [J]
am = 1.66053906660E-27 #atomic mass [kg]
au_T = hbar / Eh #[s^-1]

#velocties are given in bohr / au_t , we want it in Angstrom/fs
vel_converter = bohr / au_T * 1e-5
# Next we have forces which are given in Hartree energy / bohr, we want it in [au*Angstrom/fs^2]
force_converter = Eh / bohr / (am * 1e20)

mddata_file = '/home/tue/data/MD/ethanol/ethanol.npz'

dt = 0.5 # [fs]
nsteps = 10000

Rt,Vt,Ft,z,m = MDdataloader(mddata_file,velocity_converter=vel_converter,force_converter=force_converter)
r = Rt[0]
v = Vt[0]
a = None

VVI = VelocityVerletIntegrator(dt)
for i in range(nsteps):
r,v,a = VVI(r,v,m[:,None],a=a,F=Ft[i],F_new=Ft[i+1])

dr = torch.mean(torch.norm(Rt[i+1] - r,p=2,dim=1))
dv = torch.mean(torch.norm(Vt[i+1] - v,p=2,dim=1))
print(f"{i}, dr={dr:0.2e}, dv={dv:0.2e}")

[Tue Boesen]

Sorry, just noticed that all the tabs disappeared from the code, in case anyone want to see or use it. I have uploaded it here:

https://paste.ofcode.org/gNtfRYgasX3bzzMpJRLaLB

[Tue Boesen]

I created another dataset, where the units were returned as:

position [Angstrom]

velocity [Angstrom/fs]

Force [amu * Angstrom / fs^2]

and while the error when down significantly using this dataset where I didn't have to do any unit conversions, I'm still seeing very clear drift.

Here is a plot of the absolute position error and velocity error, clearly there is still drift/errors that are not just precision errors.