Question regarding Lost atoms and setforce commands

Hey everyone!

I am trying to do free energy calculations of ice for which I am following the Einstein molecule method. For this I need to attach each atom of all the water molecules to springs, which I am doing using spring/self command. Furthermore, I need to fix the oxygen atom of one of the molecules in my system. To do so, I am using the ‘setforce’ fix and setting the initial velocity to zero for the atom chosen to be fixed (i.e. freeze).

But, I get an error saying atoms are lost. I am unable to understand what is going wrong in my system. Although I should mention that each molecule does not interact with each other, i.e. these are ideal gas like molecules, which experience only spring force.

The timestep I am using is 1 fs. The error reads:
ERROR: Lost atoms: original 408 current 3 (…/thermo.cpp:439)
{I don’t understand if the current atom is the number of atoms lost or remaining}

The simulation runs without any error if I reduce the timestep, but I do not want to do that as I have to run a really long simulation. Also, I do not understand why should atoms be lost if there is a spring attached with a substantial spring constant of 2700 kcal/A^2.

Please explain why it happens and what is the solution apart from reducing the timestep.
I am attaching my input file screenshots.


Uploading screenshots is a waste of bandwidth and is making it a much bigger effort to check out your input compared to attaching the file
in.melt (736 Bytes)
or inserting the text as quoted text (using triple backquotes (```)) as shown below for the “melt example”. One cannot easily cut-n-paste text from images.

If you want the best help, you need to make it easy to help you. The more effort it is to read and test what you are doing, the less likely is it that people will respond (and that is an exponential or worse relation).

# 3d Lennard-Jones melt

units           lj
atom_style      atomic

lattice         fcc 0.8442
region          box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 box
mass            1 1.0

velocity        all create 3.0 87287 loop geom

pair_style      lj/cut 2.5
pair_coeff      1 1 1.0 1.0 2.5

neighbor        0.3 bin
neigh_modify    every 20 delay 0 check no

fix             1 all nve

#dump           id all atom 50 dump.melt

#dump           2 all image 25 image.*.jpg type type &
#               axes yes 0.8 0.02 view 60 -30
#dump_modify    2 pad 3

#dump           3 all movie 25 movie.mpg type type &
#               axes yes 0.8 0.02 view 60 -30
#dump_modify    3 pad 3

thermo          50
run             250

A high spring constant means atoms vibrate with a high fundamental frequency (\omega \sim \sqrt{k/m}), which means you need a shorter integration timestep to accurately integrate their motion. Typically ten timesteps per oscillation is a good starting point (you still need to double check).

A quick order-of-magnitude estimate suggests that your fundamental oscillation period would be on the order of 1fs, which suggests that your timestep would need to be on the order of 0.1fs.

1 Like

Sorry for the inconvenience. I am attaching the unit file as text.

boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic

pair_style lj/cut 8.4
pair_modify tail yes

read_data data.lmp

pair_coeff    1    1   0.000000     0.000000  # O O
pair_coeff    2    2   0.000000     0.000000  # H H
pair_coeff    1    2   0.000000     0.000000  # O H

minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

fix SHAKE all shake 0.0001 20 0 b 1 a 1

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

timestep 0.8

group one id 1:3
group rest id 4:408

variable TK equal 271.1
variable PBAR equal 25e5

velocity all create ${TK} 12345
velocity one set 0.0 0.0 0.0

fix Aspring all spring/self 2693.75
fix_modify Aspring energy yes

fix PSEUDONVT all nve

fix TSTAT all langevin 271.1 271.1 200 22547

fix fixmol one setforce 0.0 0.0 0.0

thermo 1000
thermo_style custom step cpu etotal ke pe f_fixmol[1] f_fixmol[2] f_fixmol[3] f_Aspring evdwl ecoul elong temp press vol density

dump TRAJ all custom 100 dump.lammpstrj id mol type element q xu yu zu
dump_modify TRAJ element O H

restart 100000 restart*.delta_A1

run 1250000

write_data data_delta_A1.eq.lmp

Thank you for pointing it out. It is basic physics.

So does this mean there is no alternative way but to reduce the timestep?

Or reduce the spring constant. I do not know your method (Einstein crystals), so I can’t advise you on what works and doesn’t for your method. Good luck!

1 Like

You could speed up your computation substantially by using pair style zero and picking a much shorter cutoff. Possibly, even pair style none could work. In that case, you have to set the communication cutoff appropriately via comm_modify.
It just needs to be long enough that the neighbor list code finds all bond and angle atoms for ghost atoms from the neighboring subdomains.
Perhaps 3 angstroms may suffice. LAMMPS will print a warning if it is too short.

1 Like


Thank you very much.


Fair enough. Thanks.

ERROR: Lost atoms: original 408 current 3 (…/thermo.cpp:439)

Please tell me what is 3 in this statement? 3 atoms are lost or 3 are remaining or is it the atom index of the lost atom?

Only one of those three possible explanations makes sense, if you really think about it.