Large Enthalpy Change from 1 to 8 Beads in PIMD (MBX + i-PI)

[alberto santonocito]

 Dear MBX users,  

I'm running PIMD simulations of ice Ih (D₂O) using MBX + i-PI, varying the number of beads from 1 to 64 (in steps of 8), and calculating the enthalpy as:

H=potential+kinetic_cv+P⋅V

where all quantities (potential, kinetic_cv, pressure_cv, and volume) are averaged from the out file over the production run.

However, I observe a large variation in enthalpy (and in all related quantities) between 1 and 8 beads,  much larger than expected.

I'm using the NVT ensemble with a Langevin thermostat, a timestep of 0.25 fs, and 250,000 total steps (with 50,000 steps discarded for equilibration and 200,000 used for production). These choices were made due tocomputational resources constraints.

To vary the number of beads, I only modify the line <initialize nbeads='x'> (with x = 1, 8, 16, 32, 64), keeping the rest of the input unchanged.

I’ve attached the input and output files for the 1- and 8-bead cases (note: the 8-bead run is slightly shorter than 250,000 steps, but the differences in all  quantities related to enthalpy  is evident).

Is this setup appropriate for computing enthalpy  using PIMD in MBX + i-PI?

Apologies if this is a basic question , but this is my first time attempting to compute thermodynamic quantities using PIMD with MBX and i-PI.

Thank you very much in advance for your suggestions.

Best regards,

Alberto Santonocito

[Richa Rashmi]

Dear Alberto, Thank you for your interest in using MBX. MBX uses a 2-body cutoff of 9 Å. This is the distance at which 2-body interactions are truncated in real space. Therefore, in periodic boundary conditions, the simulation box length in each direction should be greater than twice the 2-body cutoff (> 18 Å). We noticed that your box is about 15 Å along the y-axis. Please try using a larger ice box, you could generate one with GenIce (https://github.com/vitroid/GenIce). Please let us know if this resolves the large enthalpy difference you are observing between 1 bead and 8 beads. If the issue persists, we will be happy to investigate further. Thank you. Best regards, The MBX Team

[alberto santonocito]

Dear MBX Team,

Thank you very much for your prompt and helpful response.

I will follow your suggestion and generate a larger ice box using GenIce to ensure that the box dimensions meet the required cutoff conditions. I’ll re-run the simulations accordingly and will update you as soon as I have new results.

I appreciate your support.

Best regards,
Alberto

[alberto santonocito]

Dear MBX Team,

Following your suggestion, I am now using, for testing purposes, a cubic ice box with a side length of 19.125942 Å, which is larger than twice the 2-body cutoff distance (> 18 Å).

However, I am still observing a significant difference  in potential energy and kinetic_cv energy (thus also in enthalpy = potential + kinetic_cv +pV) between the 1-bead and 8-bead simulations.

To help clarify the issue, I have attached both the input files and the corresponding outputs for the 1-bead and 8-beads cases.

Thank you in advance for your time and support.

Best regards,

Alberto Santonocito

[Richa Rashmi]

Dear Alberto,

We believe the difference you are observing in potential energy between nbead = 1 and nbead = 8 is expected.
As a rough estimate:

• The ZPE for a harmonic OD oscillator is about 1250 cm⁻¹ ≈ 3.6 kcal/mol.
  
• With 216 water molecules in the box, the total ZPE is approximately 3.6 × 216 ≈ 770 kcal/mol.
  

Your simulations show a potential energy difference of about 700 kcal/mol between nbead = 1 and nbead = 8, which is consistent with this estimate.
Please let us know if there are any specific reasons you believe the values you obtained might be higher than expected.

Thank you.
Best regards,
The MBX Team