[lammps-users] potential energy

Thanks Axel,

I was using lammps to calculate the energy of a series of configurations produced by other software. So my input code is :

------- loops ----------------

variable a loop 1001

------- 3d melt of dimethoxyethane --------------

units real
dimension 3
atom_style full

neighbor 2.0 bin
neigh_modify delay 0 every 5

boundary p p p

read_data VolNrm$a.dme

pair_style lj/cut/coul/long 9.0

-------- LJ for 1-CH3 2-O 3-CH2 and Coulombic ---------

pair_coeff 1 1 0.194746 3.75
pair_coeff 2 2 0.109296 2.80
pair_coeff 3 3 0.0914113 3.95
pair_modify mix arithmetic # LJ mixing rules
pair_modify tail yes # enable long-range tail correction to energy and pressure of LJ portion
special_bonds lj/coul 0.0 0.0 0.5 dihedral no # 1-2,1-3,1-4 interactions
kspace_style pppm 0.0001 # long-range Coulombic force

-------- bond definition and fixed length, 1-CHxCHy, 2-CHxO --------

bond_style harmonic
bond_coeff 1 0.0 1.54 # Non bond energy calculated because of SHAKE
bond_coeff 2 0.0 1.41

#fix 1 sub shake 0.0001 20 50000 b 1 2 # record SHAKE statistics every 100ps

--------- angle potentials ----------

angle_style harmonic
angle_coeff 1 60.0136 112 #394026.226 112
angle_coeff 2 49.9783 112 #328137.843 112

---------- dihedral potentials -----------

dihedral_style multi/harmonic
dihedral_coeff 1 0.38456 -5.06931 -1.86879 5.11472 1.53537
dihedral_coeff 2 0.530115 -4.21128 1.99618 8.03059 1.20913

--------- system ------------------

velocity all create 303.15 2349852

compute MyPress all pressure thermo_temp

fix 2 all nvt 303.15 303.15 100

-------------- calculation of energy ---------------

thermo_style custom step pe evdwl etail ecoul elong ebond eangle edihed press
thermo_modify press MyPress

thermo 5
log log.$a.VolNrm
run 1
clear
next a
jump in.volnrm

You are adding in the tail energy twice, as it is already included in epair:

      epair = pairwise energy (evdwl + ecoul + elong + etail)

-575.19968 - (-1179.6621 + -140.00964 + 8303.8697 + -8092.3217 + 0 + 282.51382 + 110.40054 )
140.0097

2010/7/6 hong bingbing <cynthiahong1983@…215…>

Thanks, Aidan, for your calculation.

But I did not use ‘epair’ at all in my calculation. Neither did I output ‘epair’ in my log file. And the values of ‘evdwl’ and ‘etail’ and ‘ecoul’ … are generated by LAMMPS automatically. Does it mean ‘evdwl’ from LAMMPS have already taken into account ‘etail’? This is very strange. I haven’t done any modification on evdwl to ask it include etail.

Bingster

Good news, you are wrong. etail is included in evdwl, epair, epot and etotal. Just leave it out of thermo_style command and all will be well. And instead of telling us that you think this is strange, check out this doc page.

http://lammps.sandia.gov/doc/thermo_style.html

“A long-range tail correction etail for the VanderWaal pairwise energy will be non-zero only if the pair_modify tail option is turned on. The etail contribution is included in evdwl, pe, and etotal, and the corresponding tail correction to the pressure is included in press and pxx, pyy, etc.”

2010/7/7 hong bingbing <cynthiahong1983@…215…>