Difference in single point energy between GULP and LAMMPS

GULP
1

Dear Julian and users,

I have computed the single point energy of my periodic MOF system in a cubic box with unit cell length of 26.3 Å.

The 1-4 electrostatic and van der Waals interactions in my force field are both scaled by 0.5.

I computed the single point energy in both GULP and lammps, the Interatomic potential (LJ+bond term) is well matched, but the coulombic term differs by 873.68 KJ/mol between the two programs. I would like to understand the cause of this discrepancy.

In the following, I provided more details.
the output is:

image
image1393×304 32.1 KB

in the GULP input file, the key lines to scale the interaction are:
lennard geometric 12 6 o14 kcal all 0.5
6.0
coul o14 0.5
X X 6.0

the same settings in lammps are:
pair_style lj/cut/coul/long 6.0 6.0
kspace_style ewald 1e-12
pair_modify mix geometric
special_bonds lj/coul 0.0 0.0 0.5

Is anyone able to spot any mistake in these two parts of the input file?

Thanks for in advance.

Best regards,
wang

2

Dear Wang,
It’s hard to know without seeing all of both inputs & it may need a LAMMPS expert, since the GULP part looks OK. The thing we can’t see is what keyword you’ve used to activate the molecular Coulomb subtraction. This needs to be “molmec” & so if you’ve put “mole” this would be a source of difference. The other possibility would be if the charges weren’t the same, obviously.
One other thing to note, is that there will usually be a small difference in the Coulomb terms (much smaller than you have) due to differences in the precise conversion factor from 1/Angstroms to eV in the codes. This can be removed by changing the value of “inverse_angstroms_to_ev” in modules.F90 before compilation (for recent versions of GULP) to whatever LAMMPS has.
Regards,
Julian

3

Dear Julian,

Thanks for your reply.

I have posted this question in the LAMMPS forum. https://matsci.org/t/coulombic-energy-comparsion-between-same-model-with-different-bond-terms/62050

The point is that the different bond topologies in lammps and gulp. therefore, the ‘o14 0.5’ excludes different non-bonded interactions, and the monopole-monopole (total) term in gulp gives a different coulombic energy.

the problem in my case is how to exclude the Cu-Cu bond from bond topology in gulp.

I have tried to output the connectivity for both cases with or without Cu-Cu bond term, but the results are the same.

here are my gin files for both cases with or without Cu-Cu bond term.
single_no-Cu-Cu.gin (10.7 KB)
single.gin (10.6 KB)
no_Cu-Cu.got (23.5 KB)
Cu-Cu.got (23.8 KB)

for both got files, you can see the monopole-monopole (total) term is the same, but actually, the bond topology is different, hence the excluded non-bonded interaction should be also different.

How can I make gulp distinguish that the bond topology is different in two gin files.

Thank you again.

wang