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

[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

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/20210609/b8a4c1e3/attachment.htm>