Lost atoms during i-pi/lammps simulation

[Minho Kim]

Hi all,

I'm running the PIMD simulation of bulk water with 216 water molecules via lammps.

However, after 100,000 - 400,000 timestep, the lammps shows the error as "ERROR: Lost atoms: original 648 current 647 (../thermo.cpp:481)"

As I understand, this error usually occurs in lammps where the system is too unstable, but I couldn't find any system instability from trajectory or energy with time. It seems that one hydrogen atom suddenly disappears and the simulation is crashed. Shorter timestep, low number of CPUs, different neighbor atom settings, and different bead numbers show the same problem.

Do you have any idea about this error?

Here is the XML file for i-pi simulation.

<simulation verbosity="medium">
  <!-- specification of output files. property(options){units}, format, stride, etc. -->
  <output prefix="base_p-1">
          <properties stride="100" filename="out"> [ step, time{femtosecond}, conserved{kilocal/mol}, temperature{kelvin}, kinetic_cv{kilocal/mol}, potential{kilocal/mol}] </properties>
    <trajectory filename="pos" stride="1000" cell_units="angstrom" format="pdb"> positions{angstrom} </trajectory>
    <checkpoint filename="chk" stride="100000" overwrite="true"/>
  </output>
  <total_steps> 2000000 </total_steps>
  <prng>
     <seed>32345</seed>
  </prng>
  <!-- "forcefield" section. will use an external driver communicating over a UNIX socket -->
  <ffsocket name="lmpserial" mode="unix" matching="lock">
          <address> water566 </address> <latency> 0.001 </latency>
  </ffsocket>
  <!-- system specifications. this defines the physical system and the evolution "rules" (i.e. the integrator) -->
  <system>
    <!-- initialize from an external file, the nbeads attribute determines the size of the ring polymer -->
    <initialize nbeads="4">
      <file mode="pdb"> ./water-216.pdb </file>
      <velocities mode="thermal" units="kelvin"> 600 </velocities>
    </initialize>
    <!-- the <force> section specifies which ff should be used to describe interactions. as we will see,
         here one can combine multiple force providers, if desired -->
    <forces>
      <force forcefield="lmpserial"> </force>
    </forces>
    <!-- how will we evolve the system? with MD, of course! -->
    <motion mode='dynamics'>
      <dynamics mode="nvt">  <!-- constant temperature integrator - will use a simple PILE-G for the moment -->
        <thermostat mode="pile_g">
          <tau units="femtosecond"> 50 </tau>
          <pile_lambda> 0.5 </pile_lambda>
        </thermostat>
        <timestep units="femtosecond"> 0.25</timestep> <!-- this is a conservative choice for MD, but it is safe also with PIMD 300K -->
      </dynamics>
    </motion>
    <ensemble> <!-- ensemble specification. here we basically just state what is the temperature -->
      <temperature units="kelvin"> 300 </temperature>
    </ensemble>
  </system>
</simulation>

Thank you,

Minho

[Michele Ceriotti]

Hmmm. This might be due to water molecules drifting too far away. In i-PI 3.0 we switched 

from folding atoms back in the cell by default to not folding. This is because folding triggers

unnecessary neighbor list re-calculations. you should be able to fix this by using the dev version

of LAMMPS (from github) or by forcing folding

 <ffsocket name="lmpserial" mode="unix" matching="lock" pbc="true">

          <address> water566 </address> <latency> 0.001 </latency>
  </ffsocket>

Let me know if this fixes it. 

Michele

[Minho Kim]

With that setup, I started the simulation, which runs well for ~2,000,000 timesteps. (previously, stopped at 100,000 - 400,000 timesteps)

The problem seems to be fixed.

Thank you for your comment!

Minho

2025년 1월 6일 월요일 오전 6시 44분 33초 UTC+9에 michele....@gmail.com님이 작성:

[Michele Ceriotti]

Good to know. If you have the opportunity, try also the dev/ version of LAMMPS with `pbc="false"`. 

It should work, and be significantly faster as it avoids triggering a neighbor-list recalculation at (nearly) 

each step, which is why we made the change in the first place.

Michele

[Minho Kim]

Thank you for the suggestion! I'll try with dev/ version with 'pbc="false" to boost up the simulation.

Minho

2025년 1월 8일 수요일 오후 4시 38분 46초 UTC+9에 michele....@gmail.com님이 작성: