[CP2K-user] Details on cp2k's velocity verlet integrator?

[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.

[image: position_velocity_error.png]




On Wednesday, 9 June 2021 at 17:38:23 UTC-7 Tue Boesen wrote:

> 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
>
> On Wednesday, 9 June 2021 at 17:35:17 UTC-7 Tue Boesen wrote:
>
>> 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}")
>>
>>
>>
>>
>>
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20210610/b155e3c6/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: position_velocity_error.png
Type: image/png
Size: 49084 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20210610/b155e3c6/attachment.png>