Problems using hybrid potential for GAFF and reaxFF

Hongxia_Hao · July 20, 2020, 6:29pm

Dear LAMMPS users,

I have a system where I want to treat the organic molecules using GAFF force field and water molecules with reaxff force field, and the interactions between organic molecules and water with LJ potential. It’s working now but I have difficulty of controlling the temperature. Here is my input and data file.
I have 12 atom types. atom type 1 to 8 is for GAFF, and atom type 9 to 12 is for reaxFF.

Masses

1 16.0000 # oh

2 12.0100 # ca

3 1.0080 # ho

4 16.0000 # os

5 1.0080 # h1

6 1.0080 # hc

7 12.0100 # c3

8 1.0080 # ha

9 1.0080 # HW

10 15.9994 # OW

11 0.01 # El

12 1.0 # X

units real

atom_style full

boundary p p p

dielectric 1

special_bonds lj/coul 0.0 0.0 0.8

pair_style hybrid reax/c lmp_control checkqeq no lgvdw yes lj/charmm/coul/charmm 9.0 10.00000 lj/cut 9.0

bond_style harmonic

angle_style harmonic

dihedral_style harmonic

improper_style cvff

#kspace_style pppm 0.0001

read_data data.lammps

#read_restart file.restart

pair_coeff * * reax/c ffield.reax.eReaxFF_lg_3 NULL NULL NULL NULL NULL NULL NULL NULL H O El

pair_coeff 1 1 lj/charmm/coul/charmm 0.21040000 3.06647339 # oh oh

pair_coeff 2 2 lj/charmm/coul/charmm 0.08600000 3.39966951 # ca ca

pair_coeff 3 3 lj/charmm/coul/charmm 0.00000000 0.00000000 # ho ho

pair_coeff 4 4 lj/charmm/coul/charmm 0.17000000 3.00001234 # os os

pair_coeff 5 5 lj/charmm/coul/charmm 0.01570000 2.47135304 # h1 h1

pair_coeff 6 6 lj/charmm/coul/charmm 0.01570000 2.64953279 # hc hc

pair_coeff 7 7 lj/charmm/coul/charmm 0.10940000 3.39966951 # c3 c3

pair_coeff 8 8 lj/charmm/coul/charmm 0.01500000 2.59964246 # ha ha

pair_modify mix geometric

pair_coeff 1 9 lj/cut 0.000 0.000 #oh HW

pair_coeff 1 10 lj/cut 0.000 0.000 #oh OW

pair_coeff 1 11 lj/cut 0.000 0.000 #oh El

pair_coeff 2 9 lj/cut 0.000 0.000 #ca HW

pair_coeff 2 10 lj/cut 0.000 0.000 #ca OW

pair_coeff 2 11 lj/cut 0.000 0.000 #ca El

pair_coeff 3 9 lj/cut 0.000 0.000 #ho HW

pair_coeff 3 10 lj/cut 0.000 0.000 #ho OW

pair_coeff 3 11 lj/cut 0.000 0.000 #ho El

pair_coeff 4 9 lj/cut 0.000 0.000 #os HW

pair_coeff 4 10 lj/cut 0.000 0.000 #os OW

pair_coeff 4 11 lj/cut 0.000 0.000 #os El

pair_coeff 5 9 lj/cut 0.000 0.000 #h1 HW

pair_coeff 5 10 lj/cut 0.000 0.000 #h1 OW

pair_coeff 5 11 lj/cut 0.000 0.000 #h1 El

pair_coeff 6 9 lj/cut 0.000 0.000 #hc HW

pair_coeff 6 10 lj/cut 0.000 0.000 #hc OW

pair_coeff 6 11 lj/cut 0.000 0.000 #hc El

pair_coeff 7 9 lj/cut 0.000 0.000 #c3 HW

pair_coeff 7 10 lj/cut 0.000 0.000 #c3 OW

pair_coeff 7 11 lj/cut 0.000 0.000 #c3 El

pair_coeff 8 9 lj/cut 0.000 0.000 #ha HW

pair_coeff 8 10 lj/cut 0.000 0.000 #ha OW

pair_coeff 8 11 lj/cut 0.000 0.000 #ha El

With the pair_coeff between 1-8 and 9-11 as 0.0, I assume there is no interaction between the gaff organic molecules and reaxff water. But I found that the coulomb energy and vdWs energy of the whole system is not the sum of the organic molecules and the waters. The coulomb energy and vdWs energy are changed compared to the energy when I only calculate for water or for the organic molecules.

I’m wondering why that happens? Did I understand it correctly?

Thanks in advance!

Sincerely
Hongxia

akohlmey · July 20, 2020, 7:15pm

this is an incomplete input. important information is missing.

what energies exactly are you comparing and how did you determine those?

also, having no interactions between the two subsystem makes no sense. that way atoms can move “through” each other without “seeing” each other.
if that is what you want, you would be better off using pair style zero.

axel.

Hongxia_Hao · July 20, 2020, 11:00pm

Dear Axel,

Here is the remaining of the input.

neighbor 2 bin #distance

neigh_modify every 1 delay 0 check yes page 120000 one 12000 #frequency to update nblist, here every step

group cores type 1 2 3 4 5 6 7 8 9 10

group shells type 11

group solvent type 1 2 3 4 5 6 7 8

group reax type 9 10 11

compute 1 reax pair reax/c

variable eb equal c_1[1]

variable ea equal c_1[2]

variable elp equal c_1[3]

variable emol equal c_1[4]

variable ev equal c_1[5]

variable epen equal c_1[6]

variable ecoa equal c_1[7]

variable ehb equal c_1[8]

variable et equal c_1[9]

variable eco equal c_1[10]

variable ew equal c_1[11]

variable ep equal c_1[12] #this is the cgem energy

variable efi equal c_1[13]

variable eqeq equal c_1[14]

variable ereax equal v_eb+v_ea+v_elp+v_emol+v_ev+v_epen+v_ecoa+v_ehb+v_et+v_eco+v_ew+v_ep+v_efi+v_eqeq

dump 1 all custom 50 dumpmin.lammpstrj type element mass q xu yu zu

dump_modify 1 sort id

thermo 10

#thermo_style multi

thermo_style custom step etotal v_ereax pe ke density temp press ecoul v_ep

#thermo_style custom step v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq

#first min of the shells

fix 5 cores setforce 0.0 0.0 0.0 #fix the cores, all the forces are 0.

min_style cg #min of the shells with congugate gradient

minimize 1.0e-12 0.0 10000 10000 #min for shell

unfix 5 #unfix

fix 6 solvent setforce 0.0 0.0 0.0 #fix the cores, all the forces are 0.

min_style cg #min of the shells with congugate gradient

minimize 1.0e-12 0.0 10000 10000 #min for the whole system

unfix 6

minimize 1.0e-5 0.0 10000 10000 #min for the whole system

undump 1

write_data Smin.dat

write_restart Smin.rest

dump 1 all custom 50 dumpheat.lammpstrj type element mass q xu yu zu

dump_modify 1 sort id

fix 11 all momentum 1 linear 1 1 1 angular

compute 11 cores temp #print temp for cores

compute 12 shells temp

timestep 0.2

thermo 10

#thermo_style multi

thermo_style custom step etotal pe ke density temp press c_11 c_12

#fix 8 greax setforce 0.0 0.0 0.0

velocity cores create 0.0 12345

velocity shells create 0.0 12345

fix 13 cores nvt temp 1.0 300.0 100.0

fix 14 shells nvt temp 1.0 1.0 100.0

run 100000

unfix 13

unfix 14

undump 1

write_data Sheat.dat

write_restart Sheat.rest

dump 1 all custom 1000 dumpnvt.lammpstrj type element mass q xu yu zu

dump_modify 1 sort id

fix 13 cores nvt temp 300.0 300.0 100.0

fix 14 shells nvt temp 1.0 1.0 100.0

run 200000

I first defined v_reax, and do a single point calculation for just water, and found that E_tot from thermo multi is equal to v_reax, E_coul from thermo multi is equal to v_eq, and E_vdws from thermo multi is equal to the sum of c_1[1] to c_1[14] without c_1[12].

Then I added the organic molecule on top of the current water geometry, and do another single point calculation, and found that v_eq and v_ew for water changed. But I think it shouldn’t because there is no interaction between these the organic molecules and the water.

Similarly, I did a single point calculation for only the organic molecules, and then compare the E_vdws and E_coul with (E_vdws for whole system - E_vdws for water), and (E_coul for whole system - E_coul for water), these two numbers also changed.

Here I did a minimization first, and then try to heat up the system with group core temp to be 300K and group shell temp to be 1K. The shell temp becomes 2000K in 1 step.

I intend to have no interactions just for a debug purpose. They should have interactions.

Please let me know if I didn’t express it clearly.

Sincerely
Hongxia

akohlmey · July 21, 2020, 12:48am

[…]

I intend to have no interactions just for a debug purpose. They should have interactions.

Please let me know if I didn’t express it clearly.

sorry, this is all to convoluted and difficult to follow.

if you want me to look into this, i would need the following:

the data file for the “optimized” structure
the ReaxFF parameter file
3 input files using the same data file, one for each energy computation. they should only contain the minimum commands to compute and output the energies. nothing else.

some additional observations

you are not using charge equilibration. that is unusual for reaxff.
you are using an unusual reaxff model. are you sure this can work with the LAMMPS code and does not need an update to the implementation.
you don’t set the timestep. reaxff usually needs a much shorter timestep than the default for real units

axel.