Total force on an atom is different using two different methods

LAMMPS LAMMPS General Discussion
1

I calculated the total force (in the z-direction) on an atom using two methods:
i) by dumping the fz forces:
dump 6 theatom custom 100 rerun.lammpstrj id type z fz
ii) by using the group/group computation:
compute interaction theatom group/group all kspace yes
variable fz equal c_interaction[3]
fix group_group_fz all print 100 “${fz}” file fz_group_group.txt

Interestingly, I get different values for the force (fz) using the two methods above. Does anyone know the reason.
Thank you

2

You are comparing apples and oranges. It is difficult to say more because there is no information about your system and the force field and more.

3

Thanks for your comment. I have the script below. In my simulation, I’m studying the interactions between an ion and other atoms/molecules in the system. The entire system is charge neutral.

units real
dimension 3
boundary p p p
atom_style full
neigh_modify delay 0 every 1 check yes exclude type 2 2 exclude type 6 6 #exclude type 1 2
processors * * *

variable timecoef equal 1.00
variable T equal 300
variable thermo_itv equal 50
variable dump_itv equal 50
variable dumpvel_itv equal 50
variable equi_steps equal 100000
variable flow_stepsA equal 40000
variable seed equal 2012
#variable fid equal 30604

atom definition:

read_data ./…/…/Neutral_1M_1nm_chloride_drag.data #change
###############################################

group definition:

group cnt type 2 6 # silica
group oxygen type 3
group hydrogens type 1
group sod type 4
group anion type 5 # CLA
group ions type 4 5
group water type 1 3
group fluidmols type 1 3 4 5
group theion id 78137 #chloride #change
group otherions id 58639:59104 77894:78137 #other ions #change
#################################################

force fields

variable rcF equal 12.0
pair_style lj/cut/coul/long {rcF} {rcF}
pair_coeff 1 1 0.0 0.0 # H
pair_coeff 2 2 0.054 3.47 # O of SiO2
pair_coeff 3 3 0.15539 3.1656 # O of water
pair_coeff 4 4 0.3526 2.159 # Na
pair_coeff 5 5 0.0128 4.83 # anion CL
pair_coeff 6 6 0.093 4.15 # Si
pair_modify mix arithmetic
############################################
kspace_style pppm 1e-05
dielectric 1.0
bond_style harmonic
bond_coeff 1 450.0 0.9572
angle_style harmonic
angle_coeff 1 100.0 104.52

velocity cnt set 0.0 0.0 0.0 units box
fix 5 cnt setforce 0.0 0.0 0.0

minimization

#min_style sd
#minimize 1.0e-4 1.0e-5 1000 10000
#print “Min”
###############################################3
reset_timestep 0
timestep {timecoef} compute fluidtemp fluidmols temp #velocity fluidmols create {T} ${seed} dist gaussian temp fluidtemp mom yes rot yes units box
######################################################################################
fix 1 water shake 0.0001 20 0 b 1 a 1
fix 3 fluidmols nvt temp 300 300 100.0
######################################################################################################

external bias

#fix 4 all efield 0.0 0.0 0.000494 #L=101.2A 1lammps=V/A*101.2A=101.2V 0.7084
####################################################################################################
#compute energy all pe/atom pair kspace
#compute coord all coord/atom 3.3 2 # oxygen rdf distance
#dump 5 solutions custom 1000 pe.txt mol id type c_energy c_coord
#dump 55 oxygen custom 1000 coord.txt mol id type c_coord
dump 6 theion custom 100 rerun.lammpstrj id type z fz #100 follows the original simulation dumping frequency
#######################################################################################################
##group/group

compute ion_water_short_range theion group/group water kspace no
compute ion_water_long_range theion group/group water kspace yes
variable ion_water_short_range_fz equal c_ion_water_short_range[3]
variable ion_water_long_range_fz equal c_ion_water_long_range[3]

compute ion_membr_short_range theion group/group cnt kspace no
compute ion_membr_long_range theion group/group cnt kspace yes
variable ion_membr_short_range_fz equal c_ion_membr_short_range[3]
variable ion_membr_long_range_fz equal c_ion_membr_long_range[3]

compute ion_ion_short_range theion group/group otherions kspace no
compute ion_ion_long_range theion group/group otherions kspace yes
variable ion_ion_short_range_fz equal c_ion_ion_short_range[3]
variable ion_ion_long_range_fz equal c_ion_ion_long_range[3]

compute ion_all_short_range theion group/group all kspace no
compute ion_all_long_range theion group/group all kspace yes
variable ion_all_short_range_fz equal c_ion_all_short_range[3]
variable ion_all_long_range_fz equal c_ion_all_long_range[3]

fix group_group_short_range_fz all print 100 “{ion_water_short_range_fz} {ion_membr_short_range_fz} {ion_ion_short_range_fz} {ion_all_short_range_fz}” file fz_group_group_short_range.txt #100 follows
the original simulation dumping frequenc

fix group_group_long_range_fz all print 100 “{ion_water_long_range_fz} {ion_membr_long_range_fz} {ion_ion_long_range_fz} {ion_all_long_range_fz}” file fz_group_group_long_range.txt
########################################################################################################i
#reset_timestep 0
#thermo ${thermo_itv}
#thermo_style multi
rerun ./…/all.lammpstrj dump x y z #is the name the same

4

Sorry, but you have to either attach the input or enclose it in triple back quotes “```” to stop the forum from trying to typeset it. This way it is not readable.

5

My guess is that the force printout includes thermostatting forces while compute group/group doesn’t. But it is very difficult to read your input without using backquotes to turn on verbatim mode.

6

I apologize. I was not allowed to attach the file so I enclose it in back quotes:

dimension	3
boundary 	p p p
atom_style 	full
neigh_modify delay 0 every 1 check yes exclude type 2 2 exclude type 6 6 
processors 	* * *

variable timecoef equal 1.00
variable T equal 300
variable thermo_itv equal 50
variable dump_itv equal 50
variable dumpvel_itv equal 50
variable equi_steps equal 100000
variable flow_stepsA equal 40000
variable seed equal 2012
#variable fid equal 30604
# atom definition:
read_data       ./../../Neutral_1M_1nm_chloride_drag.data    
###############################################
# group definition:
group           cnt       type 2 6 # silica
group           oxygen    type 3
group           hydrogens type 1
group           sod       type 4
group           anion     type 5 # CLA
group           ions      type 4 5
group           water     type 1 3
group           fluidmols type 1 3 4 5
group           theion    id  78137                        
group           otherions id  58639:59104 77894:78137      
#################################################
# force fields
variable        rcF equal 12.0
pair_style      lj/cut/coul/long  ${rcF} ${rcF}
pair_coeff      1 1     0.0     0.0                       # H
pair_coeff      2 2     0.054   3.47                      # O of SiO2
pair_coeff      3 3     0.15539 3.1656                    # O of water
pair_coeff      4 4     0.3526  2.159                     # Na
pair_coeff      5 5     0.0128  4.83                      # anion CL
pair_coeff      6 6     0.093   4.15                      # Si
pair_modify     mix arithmetic
############################################
kspace_style    pppm 1e-05
dielectric      1.0
bond_style harmonic
bond_coeff 1 450.0 0.9572
angle_style harmonic
angle_coeff 1 100.0 104.52 

velocity  cnt set 0.0 0.0 0.0 units box
fix 5 cnt setforce 0.0 0.0 0.0

# minimization
#min_style       sd
#minimize        1.0e-4 1.0e-5 1000 10000
#print "Min"
###############################################3
reset_timestep 0
timestep ${timecoef}
compute fluidtemp fluidmols temp
#velocity  fluidmols create ${T} ${seed} dist gaussian temp fluidtemp mom yes rot yes units box
######################################################################################
fix 1 water shake 0.0001 20 0 b 1 a 1 
fix 3 fluidmols  nvt temp 300 300 100.0    
######################################################################################################
dump 6 theion custom 100 rerun.lammpstrj id type z fz 
########################################################
##group/group


compute ion_water_short_range theion group/group water kspace no
compute ion_water_long_range theion group/group water kspace yes
variable ion_water_short_range_fz equal c_ion_water_short_range[3]
variable ion_water_long_range_fz equal c_ion_water_long_range[3]


compute ion_membr_short_range theion group/group cnt kspace no
compute ion_membr_long_range theion group/group cnt kspace yes
variable ion_membr_short_range_fz equal c_ion_membr_short_range[3]
variable ion_membr_long_range_fz equal c_ion_membr_long_range[3]


compute ion_ion_short_range theion group/group otherions kspace no
compute ion_ion_long_range theion group/group otherions kspace yes
variable ion_ion_short_range_fz equal c_ion_ion_short_range[3]
variable ion_ion_long_range_fz equal c_ion_ion_long_range[3]


compute ion_all_short_range theion group/group all kspace no
compute ion_all_long_range theion group/group all kspace yes
variable ion_all_short_range_fz equal c_ion_all_short_range[3]
variable ion_all_long_range_fz equal c_ion_all_long_range[3]


variable z_position equal xcm(theion,z)

fix group_group_short_range_fz all print 100 "${z_position} ${ion_water_short_range_fz} ${ion_membr_short_range_fz} ${ion_ion_short_range_fz} ${ion_all_short_range_fz}" file fz_group_group_short_range.txt   

fix group_group_long_range_fz all print 100 "${z_position} ${ion_water_long_range_fz} ${ion_membr_long_range_fz} ${ion_ion_long_range_fz} ${ion_all_long_range_fz}" file fz_group_group_long_range.txt 
########################################################################################################i
#reset_timestep 0
#thermo ${thermo_itv}
#thermo_style    multi 
rerun ./../all.lammpstrj dump x y z
7

That is a very good point. I tried removing the thermostat from the ions but it was still different.

Changing fix 3 fluidmols nvt temp 300 300 100.0 to fix 3 water nvt temp 300 300 100.0

fluidmols include water molecules and ions and water includes only water molecules. I calculated the force by all the atoms on one of the ions using the group/group computation and it was still different from the dumped fz forces.

8

I’ve made some tests with a simpler system to see for myself.
What was missing in the discussion was some assessment as to how different the numbers were.

The conclusion is mostly:

  • you get very similar forces from looking at the net force and what you get from compute group/group so the result is not completely wrong and follows the right trends
  • however, there is a much larger error when using kspace. and I think this is due to having a non-neutral total system when doing the compute group/group calculation since you are computing the interaction of a single charged atom with a group of atoms of neutral molecules. in kspace calculations that automatically has the effect as if there is an additional compensating charge added to the system, so the forces will be impacted by that
  • when switching to a cutoff pair style and dropping kspace the agreement between compute group/group and the net force is much better.
  • you will never get a perfect agreement because the computation and accumulation of the forces is slightly different and floating point math deficiencies will prevent that.

i am seeing a relative difference of around 2-5% with kspace and 0.1-0.3% without.
since this is an analysis run done via re-run, I would suggest to skip using kspace and switching to a cutoff coulomb for the analysis. due to being the same trajectory there will be no error from the lack of kspace on positions, but the analysis run will be faster and with a smaller error.

9

since this is a rerun, thermostatting should not make a difference until one looks at velocities.