Hi Lammps Users,
I investigated some time ago how to calculate the cohesion energy of a system using lammps and I found a workable method in the literature “https://www.sciencedirect.com/science/article/pii/S0032386117311576#sec3”.
cohesion energy = system potential energy - Σ intermolecular interactions.
Where the intermolecular interactions are written in a loop implementation as follows,
# -------------------- Run --------------------
timestep 1
variable M equal mass(all)
fix print all print 1 "${mID} ${Elocal} $(pe) $(vol)" append output.txt title "# molID Elocal PE Vol"
variable i loop 100
label loopa
group m1 molecule $i
compute $i m1 group/group m1 kspace yes
variable mID equal $i
variable Elocal equal c_$i
fix 1 all nvt temp 300 300 100
thermo_style custom step ke pe v_mID v_Elocal vol v_M
#thermo_modify line multi flush yes
thermo 1
run 0
group m1 clear
next i
jump SELF loopa
The cohesion energy density was obtained by dividing the cohesion energy by the volume of the system.
However, the cohesive energy density calculated by lammps is not the same as the one calculated by materials studio. （The cohesive energy density calculated by lammps is about 5kJ/m3, whereas the cohesive energy calculated by MS is 8kJ/m3, and I think that MS may be more accurate as it is a packaged module.)
I think my current method of calculating cohesion energy via lammps is still flawed, any improvements please.
Thanks,
Cao
The sum of the compute group/group in the loop should be the self-energy of all the molecules in a periodic lattice. This energy also includes the interaction of every molecule with its periodic image, which is not how you define the intra-molecular energy.
As stated in the paper you cited, the cohesive energy is computed by subtracting “the average potential energy of the isolated chains (molecules)” from the “average potential energy of the entire simulation system in the condensed phase”.
Remember that the energy computed by compute group/group:
This compute does not calculate any bond or angle or dihedral or improper interactions between atoms in the two groups.
Therefore, the potential energy of the condensed phase must include only the pairwise terms.
Could you run again your script and print the cohesion energy by computing it in the following way?
compute intra all group/group all pair yes kspace yes molecule intra
variable enb equal evdwl+ecoul+elong
variable ced equal (v_enb-c_intra)/$(vol)
run 0
print "Cohesion energy: $(ced)"
Thank you for your suggestion. After I try, I will report on the progress
I still haven’t obtained a result that matches MS, and the van der Waals component calculated by LAMMPS is only consistent with MS.
In MS, cohesive energy=van der Waals component + electrostatic component.
So I think the problem lies in the calculation of electrostatic components using LAMMPS.
Yes, I believe the problem is the correction energy term due to the long-range solver. If you want to help debug this issue, could you provide a small test sample (e.g. a small box of water) for which you can compute the cohesion energy with MS and LAMMPS? That would be very useful.
Sure, There are in/data files for computing the correction energy
test.in (702 Bytes)
reax.data (2.5 MB)
hesion energy.
You should provide the reference value obtained for the same system with Material Studio.
With the input provided, the cohesive energy density obtained with the compute intra command (maintained locally) is 2.32E6 KJ/m3
The total of cohesive energy density is 7.619 1e8 J/m^3.
the part of van der Waals is 2.5591e8 J/m^3.
the part of electrostatic is 4.9211e8 J/m^3.
the other is 0.13961e8 J/m^3.