<div>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.</div><div> </div><div>Does anyone know where I can find some details on this?</div><div> </div><div>My current setup is the following:</div><div>_________________________________________________________________________________ </div> <div>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}") </div><div> </div><div> </div><div> </div>