<div>I created another dataset, where the units were returned as:</div><div> </div><div>position [Angstrom]</div><div>velocity [Angstrom/fs]</div><div>Force [amu * Angstrom / fs^2] </div><div> </div><div>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.</div><div> </div><div>Here is a plot of the absolute position error and velocity error, clearly there is still drift/errors that are not just precision errors. </div><div> </div><div><img alt="position_velocity_error.png" data-iml="42994" width="1603px" height="545px" src="cid:c6ac9729-2531-456f-82f8-51670efe26f7"></div><div> </div><div> </div><div> </div><div> </div><div class="gmail_quote"><div dir="auto" class="gmail_attr">On Wednesday, 9 June 2021 at 17:38:23 UTC-7 Tue Boesen wrote: </div><blockquote class="gmail_quote" style="margin: 0 0 0 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div>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:</div><div><a href="https://paste.ofcode.org/gNtfRYgasX3bzzMpJRLaLB" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=en-GB&q=https://paste.ofcode.org/gNtfRYgasX3bzzMpJRLaLB&source=gmail&ust=1623429204241000&usg=AFQjCNGZGePmlcxhCXGijUZS0Oubf-INJg">https://paste.ofcode.org/gNtfRYgasX3bzzMpJRLaLB</a> </div> <div class="gmail_quote"><div dir="auto" class="gmail_attr">On Wednesday, 9 June 2021 at 17:35:17 UTC-7 Tue Boesen wrote: </div><blockquote class="gmail_quote" style="margin:0 0 0 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><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[<a href="http://z.to" rel="nofollow" target="_blank" data-saferedirecturl="https://www.google.com/url?hl=en-GB&q=http://z.to&source=gmail&ust=1623429204241000&usg=AFQjCNF_sH2tcz5kU_KPvQnnmTCR0OBv9w">z.to</a>(dtype=torch.int64)] return <a href="http://masses.to" rel="nofollow" target="_blank" data-saferedirecturl="https://www.google.com/url?hl=en-GB&q=http://masses.to&source=gmail&ust=1623429204241000&usg=AFQjCNGakX0tjVWkjqLo12gEDL1hSrlgVQ">masses.to</a>(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> </blockquote></div></blockquote></div>