Should the quantum kinetic energy be time-step dependent?

[Xu Heseri]

Dear i-pi developers,

    I performed machine learning potential accelerated (32 beads) PIMD simulations for liquid water with time-step values of 0.5 fs, 0.25 fs, and 0.1 fs, respectively, based on the enclosed input file, using the pile_g thermostat with a time constant of 100 fs. The total simulation times were 100 ps, 100 ps, and 50 ps, respectively.

    I found that the average quantum kinetic energy values derived from these three simulations differed considerably. For H, the values from the simulations with time-step values of 0.5 fs, 0.25 fs, and 0.1 fs are 232032.0 K (1929.2 kJ/mol), 228334.7 K (1898.5 kJ/mol), and 227490.5 K (1891.5 kJ/mol), respectively. This discrepancy leads to different <exp(-beta*sc)> values (calculated using the direct scaled-coordinate estimator) for H, i.e., 9.51, 8.91, and 8.80, respectively.

    I am not sure whether the quantum kinetic energy should depend on the time-step, or if there might be some issues with my input file. Could anyone kindly provide some guidance? Thank you in advance!

Sincerely,

Xu

[Michele Ceriotti]

Hi Xu,

     it's not entirely unreasonable, 0.5fs is a bit at the edge of what is stable for water and 32 beads. One thing that MAY be making things worse is that you're using very aggressive thermostatting of the ring polymer internal degrees of freedom. You can try to run including <pile_lambda> 0.2 </pile_lambda> inside the <thermostat> block - this reduces a bit the Langevin strength and might improve the accuracy of the integrator. Let us know if this helps, otherwise you could look into the Cayley integrators, but they haven't been tested much so you're on your own there.

Best

Michele

[Mariana Rossi]

Hi Xu, in addition I did not parse your outputs, but what is the trend for the total kinetic_cv and the one for O? Do they follow the same trend or is the kinetic_cv (total) more stable wrt. time step?

[Xu Heseri]

Dear Prof. Ceriotti,

Thank you very much for your early reply! Following your guidance, I have submitted the jobs using the pile_lambda with time-step values of 0.25 fs and 0.1 fs. I will update the results as soon as they are finished.

In addition, could you please explain the “aggressive” thermostat in a bit more detail? The input file I used seems quite ordinary for PIMD simulations, and a time constant of 100 fs is less aggressive than what has been employed in other studies, i.e., about 100 times the time step. Thank you very much!

Sincerely,

Xu

[Xu Heseri]

Dear Dr. Rossi,

Thank you very much for your reply. Yes, the average total kinetic_cv values show the same trend. However, as expected, for O the differences among the average quantum kinetic energy values and <exp(-beta*sc)> values are negligible.

Specifically:

(1) The average total kinetic_cv values from the simulations with time-step values of 0.5 fs, 0.25 fs, and 0.1 fs are 275681.7 K (2292.2 kJ/mol), 271752.8 K (2259.5 kJ/mol), and 270834.3 K (2251.9 kJ/mol), respectively.

(2) For O, the average quantum kinetic energy values from the simulations with time-step values of 0.5 fs, 0.25 fs, and 0.1 fs are 43649.8 K (362.9 kJ/mol), 43418.2 K (361.0 kJ/mol), and 43343.8 K (360.4 kJ/mol), respectively; and the <exp(-beta*sc)> values are all equal to 1.06 within the error.

Best regards,

Xu

[Michele Ceriotti]

PILE is designed to act in ring-polymer normal modes basis. The tau defines the time constant of the centroid (and indeed 100fs global is not aggressive at all) while the internal modes are termostatted based on the optimal coupling time of the free ring polymer modes. However, given that for many of these this is smaller than the time step, I had found empirically it made the stability of the split integrator with time step a bit worse. pile_lambda scales down the coupling strength and can maybe (guess you can tell now?) improve a bit integrator stability.

Best

Michele

[Xu Heseri]

Dear Prof. Ceriotti,

Thank you very much for such a detailed explanation! Because of the use of the direct scaled-coordinate estimator, my simulations are slower than expected. The average values of total kinetic_cv, H kinetic_cv, O kinetic_cv, and the <exp(-beta*sc)> values for H and O are listed below. It seems that pile_lambda has little effect on these values. I am wondering if you think the differences between the results obtained from simulations with time-step values of 0.25 fs and 0.1 fs are acceptable. Thank you for your guidance!

Sincerely,

Xu

[Michele Ceriotti]

Thanks for testing. I am surprised it is not yet fully converged by 0.25fs, but the fractionation ratio depends on such tiny energy differences that maybe this is not too surprising. Personally, I'd consider a difference of 0.1 on the ratio to be negligible. The only other thing that I could suggest is these Cayley integrators, but the implementation has never been truly tested AFAIK. 

Best,

Michele