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

[alberto santonocito]

Dear Richa Rashmi and MBX Team,

Thank you once again for always taking the time to respond,  I truly appreciate your support.

Your explanation regarding potential energy is clear to me, but I’m still unclear about the behavior of other quantities like kinetic energy. In particular, I don't observe the expected trend in internal energy (U) or enthalpy (H) as the number of beads increases.

To clarify my goal: I’m trying to determine the optimal number of beads needed to accurately capture quantum effects in path-integral molecular dynamics, specifically for calculating enthalpy. I expected that increasing the number of beads would cause H (and U) to decrease and eventually plateau,  indicating convergence. Identifying this point is crucial for accurately comparing the enthalpies of different ice phases.

Given that your semi-empirical potential is considered highly accurate for ice, I’m very interested in using it. However, I may be misinterpreting the results or having an error in the methodology.

I've attached an Excel file with my analysis.

Thank you again for your time and guidance.

Best regards,

Alberto Santonocito

[Richa Rashmi]

Dear Alberto,

Thank you for your email. The trend you observe in the internal energy with respect to the number of beads looks correct. As you increase the number of beads, you should see the values converge further.

It might be helpful to discuss this in more detail over a Zoom call so we can address your questions more effectively. Please let us know if that works for you, and we will be happy to coordinate a time that suits us both.

Thank you.
Best regards,
The MBX Team

[alberto santonocito]

Dear Richa Rashmi ,

Thank you once again for your feedback and prompt response. 

I think that now the trend of the data with beads is clear also to me.

Thank you again.

Best regards,

Alberto