energy is not conserved in PIMD

[Bence Balázs Mészáros]

Dear MBX team,

I would like to run PIMD water simulations in condensed phase using MBX and i-pi, but I have a problem with the conserved energy linearly decreasing (very significantly, ca. by 200 eV after 0.2 ps). I have attached all inputs and outputs in a compressed zip file. These are my main parameters:

256 water molecules (I got the input structure from the MBX examples)

32 beads

298.15 K

Pile_g thermostat, tau =100 fs

0.25 fs timestep

mbx.json:

{
   "Note" : "This is a cofiguration file",
   "MBX" : {
       "box" : [19.7295,0.0,0.0,0.0,19.7295,0.0,0.0,0.0,19.7295],
       "twobody_cutoff"   : 9.0,
       "threebody_cutoff" : 5.0,
       "max_n_eval_1b"    : 500,
       "max_n_eval_2b"    : 500,
       "max_n_eval_3b"    : 500,
       "dipole_tolerance" : 1E-12,
       "dipole_max_it"    : 100,
       "dipole_method"     : "aspc",
       "alpha_ewald_elec" : 0.6,
       "grid_density_elec" : 2.5,
       "spline_order_elec" : 6,
       "alpha_ewald_disp" : 0.6,
       "grid_density_disp" : 2.5,
       "spline_order_disp" : 6,
       "ignore_2b_poly" : [],
       "ignore_3b_poly" : [],
       "ttm_pairs" : []
   } ,
   "i-pi" : {
       "port" : 34567,
       "localhost" : "6351017"
   }
}

First, my problem appears to be with MBX here, as I've run a PIMD with stable conserved energy using qtip4pf in i-pi with the same settings.

I've also tried classical MD simulations, first NVE, then NVT. When using 0.5 fs timestep in the NVE simulation, I observed only a slight decrease in conserved energy (ca. 0.2 eV after 1 ps). This slight decrease disappeared when I used a timestep of 0.25 fs. When using a stochastic velocity rescaling thermostat in an NVT simulation (tau = 100 fs), with a time-step of 0.5 fs, results were similar (ca. 0.2 eV decrease after 1 ps).

These are the modifications in my PIMD input I've tried (separately) to resolve this problem, but they didn't help:

1. I have decreased the dipole tolerance to 1E-16

2. I have increased the three-body potential to 7.0 Å

3. I have used a slightly larger system of 278 waters, with a cell size of 21.0 Å

4. I have decreased the timestep to 0.1 fs (the decrease in conserved energy was smaller, but still large)

5. I have set the fixcom setting of i-pi to false

Do you see any error in my inputs that might cause this behavior?

Thank you in advance for your help,

Best Regards,

Bence Mészáros

[Philip Zhu]

One thing I noticed is that you’re using “aspc” in json file. I am guessing there might be a problem specific to how i-pi assign calculation jobs to MBX drivers. What happens if you change it to “cg” (conjugate gradient)?

[Bence Balázs Mészáros]

Thanks a lot for the quick reply. Now it appears to be working fine after a short simulation. Should I use "cg" in all of my simulations with i-pi?

Best regards,

Bence Mészáros

[Philip Zhu]

Yes, I would recommend using cg. In my understanding, cg always iteratively solve for the self-consistent induced dipoles for each timestep until the tolerance is met, whereas aspc is supposedly faster (but sometimes it’s tricky and not necessarily faster) by using information from previous timesteps to make a one-shot estimation to solve the induced dipoles and it doesn’t care much about the dipole tolerance. However, when i-pi sends multiple calculation tasks to mbx drivers, I’m suspecting that the mbx drivers was not able keep track of if the previous calculation is the “last timestep” from the same bead, or it’s just some other beads from the same step, so I think the result is probably incorrect. 

[Chongbin Wang]

Thanks a lot, I'm very happy to see your discussion! I have been confused by energy conservation in PIMD for one week. I will try to use "cg" in my simulation.