Equilibrate systems with fix pimd/nvt and pimd/langevin


I’m working on an NPT equilibration for cryogenic components using the PIMD method.

The issue can be seen in both the stable version (2 Aug 2023) and te latest version (7 Feb 2024 - Update 1).

LAMMPS documentation gives a pretty clear explanation of the two available methods: The older one, “fix pimd/nvt” applies the Nose-Hoover chain thermostat, and the newer one, “fix pimd/langevin” applies Langevin thermostat, which allows it to run in nve, nph, not and nve ensembles.

I tried both of the ensembles but found “fix pimd/nvt” could give physical movement to the particles (appears similar to the classical MD, probably because the thermostat applies in cartaisian coordinates).

fix 1 all pimd/nvt method nmpimd fmass 1.0 sp 1.0 temp 20.0 nhc 24 # 3 x N=1 x P=8 x Nc=1 ?
run 1000000

with the input command specifying an 8-bead simulation

mpiexec -np 8 lmp_mpi -in test.in -p 8x1

However, the position of the particles for the “fix pimd/langevin” remains fixed in space and vibrates (probably caused by the interaction with other beads). In npt ensemble the box changed in dimension very slightly and I’m pretty sure the particles are not moving physically like what we see in classical MD.

fix 1 all pimd/langevin ensemble npt temp 20 thermostat PILE_L 1234 tau 1.0 iso 1.97385 barostat MTTK taup 1.0
run 1000000

with the same input command.

mpiexec -np 8 lmp_mpi -in test.in -p 8x1

I have several questions:

Q1. Because we’re aiming to equilibrate the system in an npt ensemble, I assume fix pimd/langevin is the preferred option. But since it cannot show the system ergodicity, I tried the following hybrid fix to make the simulation seem correct:

fix 1 all npt temp 20 20 5.0 iso 1.97385 1.97385 500.0
fix 2 all pimd/langevin ensemble nvt temp 20 thermostat PILE_L 1234 tau 1.0


fix 1 all npt temp 20 20 5.0 iso 1.97385 1.97385 500.0
fix 2 all pimd/nvt method nmpimd fmass 1.0 sp 1.0 temp 20.0 nhc 24

Will it be physically accurate? If not, is there any recommended setup for npt equilibrations?

Q2. Our job runs on mpi condition, so every time we submit the PIMD simulation, multiple log files (In my case, 8) will be generated. If we want to evaluate the thermodynamic properties within the system (e.g. density, enthalpy), shall we average their values?

Q3. Additionally, I’m aware the partition functions etc, are calculated in a slightly different manner in PIMD. Are the energy terms like enthalpies recorded in the log.lammps file reflecting the real system properties? Or are they simply using the classical ways to compute them so we cannot trust those values?

Any advice is greatly appreciated.

No, you must not have two different time-integration fixes acting on the same atoms at the same time.

All thermodynamic output in LAMMPS uses the properties of the classical system in each replica. It knows nothing about this being a path-integral simulation.

I know too little about the rest, but I have notified the author of fix pimd/langevin to comment further on your post.

1 Like

Thanks a lot! I appreciate your help.

Hi Oliver,

I am the one who wrote fix pimd/langevin. I am sorry that I have just seen your post. I hope this reply is not too late.

Regarding your questions:

Q1. fix pimd/langevin is the only option to run NpT PIMD in LAMMPS. I have not encountered any erdodicity issue with it. Could you please tell me how you observed the particles are not moving in a physical way? This issue sounds a little bit serious and I need to fix it anyway. If you can let me test your system, you can write me an email at [email protected].

Q2. The density is the same among all the beads (since the barostat changed the volume of all the beads on the same footing), so you can use density in any log file. As for enthalpy, since the particle positions for all the beads differ, the enthalpy in each log file represents U(R)+P(R)V for each bead, where R represents the coordinate for each bead.

Q3. If you need PIMD-specific thermodynamic quantities, you can use the output of the fix style “pimd/langevin”. I provided a couple of thermodynamic quantities. The list is in the “Restart, fix_modify, output, run start/stop, minimize info” part of the document page.

If you need any help on using fix pimd/langevin, please let me know!