Why is spring energy non zero for the initial configuration at zeroth time-step after the fix self/spring is issued?

LAMMPS LAMMPS General Discussion
1

Hi !

I am using self/spring fix to attach springs to the atoms in my system to calculate the spring energy. This fix describes that the atoms are tethered to their initial positions. The initial positions are the positions when the fix is issued. Does this not mean that the initial position itself is the equilibrium positions of the atoms and hence the spring energy for the time-step when the fix is issued should be zero, since the spring would neither be compressed or stretched?

Thanks and regards,
Chaitanya

2

You have to provide more specific and an example that reproduces this.

If I add the line

fix  0 all spring/self 10.0 xyz

as the first fix to the example/peptide input deck and also use a custom thermo_style:

thermo_style custom step temp press epair emol f_0

Then there is no problem:

Setting up Verlet run ...
  Unit style    : real
  Current step  : 0
  Time step     : 2
SHAKE stats (type/ave/delta/count) on step 0
Bond:    4   1.111     1.44264e-05        9
Bond:    6   0.996998  7.26967e-06        6
Bond:    8   1.08      1.32536e-05        7
Bond:   10   1.111     1.22749e-05        8
Bond:   12   1.08      1.11767e-05        9
Bond:   14   0.96      0                  1
Bond:   18   0.957206  4.37979e-05     1280
Angle:  31   104.519   0.00396029       640
Per MPI rank memory allocation (min/avg/max) = 19.4 | 19.4 | 19.4 Mbytes
   Step          Temp          Press          E_pair         E_mol           f_0      
         0   282.10052     -837.01119     -6442.768       70.391395      0            
        50   259.04253     -2570.3897     -6826.1722      61.265128      659.09508

Of course, if I (incorrectly) place it after fix shake, then fix shake will have already modified the positions to fulfil the constraints and thus there will be some (small) restoring force.

   Step          Temp          Press          E_pair         E_mol           f_0      
         0   282.10052     -837.01119     -6442.768       70.391395      8.6124166e-07
        50   258.79098     -2539.4712     -6814.8979      61.261388      656.76934
3

Okay. Thanks for the explanation. But I am getting a rather large value of the spring energy.
This is my input script followed by output of the thermo command. Can you help me with where I am going wrong?


units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic

#special_bonds lj/coul 0.0 0.0 0.5

# remove hybrid if not necessary
pair_style lj/cut 8.4
pair_modify tail yes

read_data Ice.NPT.eq.lmp
# read_restart restart1.lmp

set atom * charge 0.0
# remove pair style if not using hybrid
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.5

group one id 1
group rest id 2:864

variable TK equal 200.0
variable PBAR equal 24.67

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

fix Aspring all spring/self 3975
fix_modify Aspring energy yes

compute 1 all temp

fix avgT all ave/time 100 1000 100000 c_1 mode scalar file avgT.dat

fix PSEUDONVT rest nve


### Output

f_spring

 1877.269
1088.6063
1136.3039
1141.9676
1098.9832
1125.4292
1146.8491
1135.0299
1117.8366
1108.4801
1137.2973
1126.7825
1128.3454
1118.4373
1109.6121
1127.3284
1111.6018
1134.1626
 1147.506
1138.0841
 1106.736
1070.8623
1152.6479
 1139.049
1133.6159
 1107.534
1108.4305
 1105.962
1124.8217
1138.9614
1117.3752
1138.0984
1127.7198
1124.4399
1131.9615
4

Please read my previous message again and more carefully.
Fix shake should be defined after any fixes that add forces to the system.
So your input is not correct.

The unexpected spring force will be larger when the constraints are not well fulfilled with the initial geometry, and thus fix shake will change the positions more than in the example that I have given. My example is started from the geometry of a previous run and thus the difference of the initial geometry to the constrained geometry is only due to the loss of precision from storing the geometry in a data file.

5

I understand that I have to apply the fix SHAKE after all the fixes that apply forces on the molecules, but I am reading a data file written at the end of another simulation in which a same SHAKE command was already used. My point being, the bond and angle constraints to be applied by SHAKE were already satisfied while generating the input configuration. Hence, wouldn’t it be correct to think that the SHAKE will not move atoms by much to fulfill the constraints? That is why I am a little surprised looking at the spring energy at 0th timestep.

6

I don’t know what is in your data file, I don’t know what you have run before, and don’t know how exactly you are running things now. It is very strange for example that you are using a LJ pair potential with energies set to zero. Did you do that before?

What I do know is that if you are doing the restarting correctly and are running the fixes in the correct order, there would be no spring force in the first step.

Bottom line, I have complete trust in the LAMMPS features you are using working correctly, so any discrepancy has to come from what you are doing, but since I have only limited knowledge of that, I cannot point out where the culprit is.

7

I tried running the simulation with fix SHAKE after all other force setting fixes as you have explained and indeed got zero energy at the zeroth time-step. However, my point was just that it is surprising to to get energies in the order of 10^3 kcal/mol just by the alleged little shifting of atoms due to SHAKE.

Thanks for your help!