Energy is not conserved

[Nagaprasad]

Dear MBX users and developers,

I have run an MD of liquid water with 64 water molecules, equilibration with NPT followed by a production run with NVT. I have run around 1 ns trajectory and analyzed the energy conservation. It seems the total energy is not conserved, and it is increasing with time. The same was observed with NVT data as well. I think I am doing things correct only. Here are the images of the NPT and NVT energy conservation.  Also attached are the xml files of NPT and NVT. For the NVT simulations, I just copied the checkpoint file of the NPT simulations to xml and changed NPT to NVT, hope that's enough.

That would be great if someone assist and explain the possible reasons for this odd behaviour.

Thanks in advance.

Naga

[Francesco Paesani]

Hello Naga,

Thank you for your interest in MBX.

This is not the expected behavior. The box size is too small for all the 2-body and 3-body of MB-pol to be calculated correctly. As discussed in the original papers (https://doi.org/10.1021/ct400863t and https://doi.org/10.1021/ct500079y), due to the 2-body and 3-body cutoffs, the box dimensions for simulations in periodic boundary conditions must be at least 18 A. 

Our suggestion is to use a box with 256 water molecules. Running NVE simulations would also be a better way to check the energy conservation. Should you need further clarification or run into any more issues, please do not hesitate to reach out. 

All the best,

The MBX Team

[Nagaprasad]

Dear Francesco,

Thank you so much for your quick response.

I thought that 64 water molecules were easier to play around, so I started with this. I will increase the cutoff and see.

Thank you for the suggestion of using 256 waters, I will try this as well.

Thank you

Naga

[Xuanyu Zhu]

Hi Naga,

Too add to what Francesco said, also note that the quantity you plotted is not the total energy. It was the "conserved quantity" associated with the thermostat. ( https://ipi-code.org/i-pi/tutorials.html#:~:text=In%20general%2C%20the%20time%20step%20should%20be%20made%20as%20large%20as%20possible%20without%20there%20being%20a%20drift%20in%20the%20conserved%20quantity. )

Any major drift like yours is not the expected behavior, but from this i-pi tutuorial, a "small shift" seems to be ok. ( https://ipi-code.org/i-pi/tutorials.html#:~:text=Finally%2C%20let%20us%20check,sufficiently%20short%20time%20step )

[Nagaprasad]

Thank you Xuanyu!

Dear Francesco,

I have tried as you suggested by changing the cutoff of 2-body and 3-body terms, and I always encountered the errors "** Error ** : ** ERROR ** C++ Exception in function BoxVecToBoxABCabc, in line 102 of file tools/math_tools.cpp. Y and Z components of first vector in box must be 0. Please double check your box definition." I have found a similar post and understand that the error is related to odd water formation during the simulations. I tried to correct the waters which looked odd and continued the simulations. The error repeats every 1000-5000 steps and I don't understand the solution. That would be great if someone could have any hints to avoid the terminations. Here I have copied the mbx.json file followed by the error message.

{

   "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"   : 19.0,
       "threebody_cutoff" : 19.0,
       "max_n_eval_1b"    : 500,
       "max_n_eval_2b"    : 500,
       "max_n_eval_3b"    : 500,
       "dipole_tolerance" : 1E-16,
       "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,
       "ttm_pairs" : [],
       "ignore_2b_poly" : [],
       "ignore_3b_poly" : []

   } ,
   "i-pi" : {
       "port" : 34543,
       "localhost" : "localhost"
   }
}

 # i-PI loaded input file:  config.xml
 @SOCKET:   Client asked for connection from . Now hand-shaking.
 @SOCKET:   Handshaking was successful. Added to the client list.
 # Average timings at MD step       0. t/step: 4.57780e-01
 # Average timings at MD step    1000. t/step: 1.62587e-01
 # Average timings at MD step    2000. t/step: 1.38150e-01
 # Average timings at MD step    3000. t/step: 1.35500e-01
 # Average timings at MD step    4000. t/step: 1.73780e-01
/workspace/naga/package/i-pi-master/ipi/engine/barostats.py:604: RuntimeWarning: overflow encountered in exp
  expq, expp = (np.exp(v * halfdt), np.exp(-v * halfdt))
/workspace/naga/package/i-pi-master/ipi/engine/barostats.py:609: RuntimeWarning: invalid value encountered in add
  self.nm.qnm[0, :] += ((expq - expp) / (2.0 * v)) * (
/workspace/naga/package/i-pi-master/ipi/engine/barostats.py:614: RuntimeWarning: invalid value encountered in multiply
  self.cell.h *= expq
/workspace/naga/package/i-pi-master/ipi/engine/barostats.py:609: RuntimeWarning: invalid value encountered in multiply
  self.nm.qnm[0, :] += ((expq - expp) / (2.0 * v)) * (
/workspace/naga/package/i-pi-master/ipi/utils/mathtools.py:259: RuntimeWarning: invalid value encountered in scalar multiply
  ih[0, 1] = -ih[0, 0] * h[0, 1] * ih[1, 1]
/workspace/naga/package/i-pi-master/ipi/utils/mathtools.py:260: RuntimeWarning: invalid value encountered in scalar multiply
  ih[1, 2] = -ih[1, 1] * h[1, 2] * ih[2, 2]
/workspace/naga/package/i-pi-master/ipi/utils/mathtools.py:261: RuntimeWarning: invalid value encountered in scalar multiply
  ih[0, 2] = -ih[1, 2] * h[0, 1] * ih[0, 0] - ih[0, 0] * h[0, 2] * ih[2, 2]
 ** Error ** : ** ERROR ** C++ Exception in function BoxVecToBoxABCabc, in line 102 of file tools/math_tools.cpp. Y and Z components of first vector in box must be 0. Please double check your box definition.
 !W!  @SOCKET:   Inconsistent client state in dispatch thread! (III)
 !W!  @SOCKET:   Client  died or got unresponsive(C). Removing from the list.

Thanks

Naga

[Francesco Paesani]

Hi Naga,

I believe that your input has some problems. The cutoff should be 9.0 A. Please refer to the examples we provide for i-PI.

As an additional note, you can also use "dipole_tolerance" : 1E-8, which makes the simulations faster without losing accuracy.

All the best,

The MBX Team

[Nagaprasad]

Dear Francesco,

I think I misunderstood the conversations. From the first conversation what I understood is that the box size and the cutoffs should be at least 18, so I modified cutoffs to 19 A. In the first conversation I was using the box size > 19 A, is relatively good. So I thought that the problem is with cutoffs, and changed them to 19 A. Here I am attaching all my first inputs. Please have a look if something is wrong.

Thank you so much for your time

Naga

[Francesco Paesani]

Hi Naga,

The box dimensions must be at least 18 Å, which implies that the cutoff can be, at most, as large as half of that value to prevent each molecule from interacting with itself in real space. Please refer to the examples, as they should guide you in setting up the inputs for MBX.

If you still encounter problems with your simulation, it is likely that there may be some issues with the initial configuration. We usually suggest using Packmol to generate the initial configuration and define in the input box dimensions that are ~1.5-2.0 Å larger than the values used to set up the box in Packmol. This helps to avoid clashes at the boundaries, as Packmol does not account for periodic boundary conditions and only packs molecules within a box without considering their interactions across the boundaries.

All the best,
The MBX Team

[Nagaprasad]

Dear Francesco,

I have tried with 256 water molecules, and I copied the inputs from the examples of MBX_HOME/plugins/i-pi/examples/molecular_dynamics/condensed_phase/npt/256h2o/

and here also the conserved quantity (kilocal/mol) is decreasing constantly. I don't understand whether I can proceed further with production runs.

That would be great if someone could give some clarification on this

Thanks

[Francesco Paesani]

Hi Naga,

Thank you for your interest in MBX.

I encourage you to refer to Xuanyu's post above regarding the conserved quantity in i-PI since MBX doesn't control it. We always recommend performing NVE simulations as a preliminary step to check for energy conservation and to ensure that the equations of motion are correctly integrated for a given timestep. Once this test is successful, which implies that MBX is calculating both energies and forces correctly, the same timestep can be confidently used in other ensembles. For NVT and NPT simulations, adjusting the timescales of the thermostat and barostat may be necessary to suit a specific system, as with any i-PI simulation (MBX doesn't control that).

For simulations other than path integral molecular dynamics (PIMD), we suggest using LAMMPS, which benefits from its own MPI parallelization capabilities to speed up all classical MD simulations.

If you have any further questions or need additional guidance, please feel free to reach out.

All the best,
The MBX Team

[Daniel Count]

Dear MBX-fellows,

maybe do just add a bit more to the initial discussion upstairs:

We do get drift-free energy conservation (within the numerical stability of the underlying velocity Verlet integrator) for a NVE run with a time step of 0.5 fs for a PBC system consisting of 64 H2O molecules with a box length of 12.4 Angstroem (and cutt-off of 6.2 Angstroem ) using MBX-1.0 as force driver in our in-house MD code. So no real need to move to a box size of 18 Angstroem for testing/troubleshooting.

Regards,

Daniel

[Francesco Paesani]

Hi Daniel,

Thank you, this is reassuring. Nevertheless, a cutoff smaller than 9.0 A implies that some 2-body and 3-body PIP contributions are not properly taken into account. Since the PIPs in MB-pol are only active in the real space (i.e., there is no Ewald analog for them) and have a range up to 9.0 A, some of their contributions are not accounted for if the box is too small.

All the best,

The MBX Team