proper usage of hybrid usage for reaxFF and GAFF force field

LAMMPS LAMMPS Mailing List Mirror
1

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. I tried but failed to make it work. Here is my input and data file. May you have a look and help me figure out how to make the pair_coeff command correct?
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

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

bond_style harmonic

angle_style harmonic

dihedral_style harmonic

improper_style cvff

#kspace_style pppm 0.0001

read_data sys.lammps

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

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 # 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 geometri

pair_coeff 1 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #oh HW

pair_coeff 1 10 lj/charmm/coul/charmm 0.146495 3.1266467 #oh OW

pair_coeff 1 11 lj/charmm/coul/charmm 1.000000 0.0000000 #oh El

pair_coeff 2 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #ca HW

pair_coeff 2 10 lj/charmm/coul/charmm 0.093658956 3.29213447 #ca OW

pair_coeff 2 11 lj/charmm/coul/charmm 1.000000 0.0000000 #ca El

pair_coeff 3 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #ho HW

pair_coeff 3 10 lj/charmm/coul/charmm 0.0000000 0.0000 #ho OW

pair_coeff 3 11 lj/charmm/coul/charmm 0.0000 0.0000 #ho El

pair_coeff 4 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #os HW

pair_coeff 4 10 lj/charmm/coul/charmm 0.13168 3.0925785 #os OW

pair_coeff 4 11 lj/charmm/coul/charmm 1.0000 0.0000 #os El

pair_coeff 5 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #h1 HW

pair_coeff 5 10 lj/charmm/coul/charmm 0.0400175 2.80689787 #h1 OW

pair_coeff 5 11 lj/charmm/coul/charmm 1.0000 0.0000 #h1 El

pair_coeff 6 9 lj/charmm/coul/charmm 0.0000 1.324766395 #hc HW

pair_coeff 6 10 lj/charmm/coul/charmm 0.0400174 2.918766835 #hc OW

pair_coeff 6 11 lj/charmm/coul/charmm 1.0000 0.0000 #hc El

pair_coeff 7 9 lj/charmm/coul/charmm 0.0000 1.699834755 #c3 HW

pair_coeff 7 10 lj/charmm/coul/charmm 0.0458 3.293835195 #c3 OW

pair_coeff 7 11 lj/charmm/coul/charmm 1.0000 0.0000 #c3 El

pair_coeff 8 9 lj/charmm/coul/charmm 0.0000000 0.0000000 #ha HW

pair_coeff 8 10 lj/charmm/coul/charmm 0.039115 2.878830 #ha OW

pair_coeff 8 11 lj/charmm/coul/charmm 1.0000 0.0000 #ha El

group greax type 9 10 11 12

compute 1 greax pair reax/c

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

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

dump_modify 1 sort id

thermo 10

thermo_style multi

#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 sd #min of the shells with congugate gradient

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

unfix 5 #unfix

With command pair_coeff * * reax/c ffield.reax.eReaxFF_lg_3 NULL NULL NULL NULL NULL NULL NULL NULL H O El X, it returned me error, All pair coeffs are not set.
With command pair_coeff * * reax/c ffield.reax.eReaxFF_lg_3 O C H O H H C H H O El X, it returned me nan error for energy term E_vdWs.
With command ,it returned me error pair_coeff * * reax/c ffield.reax.eReaxFF_lg_3 H O El X. Incorrect args for pair coefficients.

How may I do this correctly?

Sincerely
Hongxia

2

you have 12 atom types, but you define no interactions between the first 8 atom types and the 12th atom type. hence the error message.

axel.

3

Dear Axel,

Thank you very much! That works now. Regarding the input file above, may I ask: could I open the kspace energy calculation in this hybrid system by changing the lj potential from lj/charmm/coul/charmm to lj/charmm/coul/long? Would that physically correct?

Sincerely
Hongxia

4

Dear Axel,

Thank you very much! That works now. Regarding the input file above, may I ask: could I open the kspace energy calculation in this hybrid system by changing the lj potential from lj/charmm/coul/charmm to lj/charmm/coul/long? Would that physically correct?

No. ReaxFF has its own long-range handling based on Wolf summation. Adding a kspace style would double count the long-range contributions, since that cannot be limited to specific groups of atoms. Technically, you might switch from using lj/charmm/coul/charmm to lj/charmm plus coul/wolf, but that is only a minor change and I don’t think it is worth the complication, since you are facing other problems from your setup. The bigger issue is that you have charge equilibration for ReaxFF and you must not have it for the CHARMM atoms. That will lead to a systematic over-polarization of ReaxFF atoms close to the CHARMM atoms, since those include an averaged effective polarization that is not impacted by the charges in the ReaxFF part. That will also lead to energy leakage. Basically, you are facing a lot of the issues that also apply to QM/MM simulations, and those can have a substantial impact on the error margins. It is almost always preferable to do all-ReaxFF calculations over mixed systems.

Axel.

5

Dear Lammps users,

We now use Lammps to simulate a nickel microforming problem. lattice fcc 3.52. In order to reduce the calculation scale, I plan to simplify it to a two-dimensional problem. dimension 2, then lattcie sq2 3.52/2 (is this okay)? At the same time, if it is handled as a 2D problem, can Ovito or other post-processing software perform dislocation analysis (DXA) because there is no FCC, BCC, etc… Can the grain boundaries be displayed?
Any suggestions are appreciated very much!

Best,
Wugui
July 11, 2020

6

it really is up to you to decide whether a flat sheet of atoms will provide the same quality of answers that you are looking for than a full 3d system. if you are unsure, repeat some calculations where the results are known and see whether you can achieve satisfactory results.

my personal opinion on your approach is that it usually is not a good idea to simplify a model only because of the computational cost associated with doing it properly and without good physical justification. you just take a large risk on the credibility of your work, and these days it is much easier to get access to additional computational resources in local or external HPC centers, than arguing away flaws in a model after having done all the work.

axel.

7

Dear axel,

Thanks very much for your helpful suggestion.

Best,
WG

8

Another issue to consider is what interatomic potential you are using
in your model. If you are using EAM (embedded atom method) that
it is parameterized as a 3d model. I.e. the densities it expects to compute
come from a 3d neighborhood. You can use it to model a surface
or defects in 3d. But using it to model a 2d system with the
same lattice constant as a 3d system, will likely not give physical
results.

Steve