Unable to run Hybrid forcefield with multi-cores system using LAMMPS

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS users,

I am trying to perform molecular dynamic simulation using Hybrid forcefield in the LAMMPS. My HybridFF MD run successfully with single core, however, when I increase number of cores, I get below warning and my simulation couldn’t start.

WARNING: Fix qeq/reax CG convergence failed after 200 iterations at 0 step (…/fix_qeq_reax_omp.cpp:567)
WARNING: Fix qeq/reax CG convergence failed after 200 iterations at 0 step (…/fix_qeq_reax_omp.cpp:567)

Anyone please help me how to resolve this issue?

Looking forward.

Regards
Khizar Hayat

2

Are you trying to use multiple Reaxff instances in a hybrid setting?
That would be a very bad idea.
The fact that the calculations completes with a single core is no indication that it is correct.

But to say anything more specific, you need to provide more information like the LAMMPS version, the platform you are running on, and your input.

3

Dear akohlmey,

Thank you for your response.

The LAMMPS version is “LAMMPS 64-bit 10Mar2021”. The workstation that I am using is based on two processors with;

Processor Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz, 2394 Mhz, 20 Core(s), 40 Logical Processor(s)
Processor Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz, 2394 Mhz, 20 Core(s), 40 Logical Processor(s)

Below is the input script for my HybridFF MD simulation.

REAXff - MFI type ZTC

-------------------- Init Section ----------------------

units real
dimension 3
boundary p p p
atom_style full
pair_style hybrid/overlay reax/c NULL checkqeq no lj/cut/coul/long 12
kspace_style ewald 1.0e-5

-------------------- Atom Definition Section ------------

read_data Cathode-min.data

#change_box all x final 0 20 y final 0 20 z final 0 20

Li-air battery electrolyte

group dmso molecule 1 # density is 1.0956 g/cm3
group Li molecule 2
group salt molecule 3

Cathode structure

group C molecule 4 # Zeolite templated carbon structure

Reactants for electrochemical reaction

group oxygen molecule 5

#--------------------- Pair Coefficients ------------------------

pair_coeff * * reax/c ffield.reax_mf.LIS NULL NULL NULL NULL Li NULL NULL C O NULL

Lennard jonnes parameters

pair_coeff 1 1 lj/cut/coul/long 0.1200 3.4000 # O(dmso)-O(dmso) LJ
pair_coeff 1 2 lj/cut/coul/long 0.2049 3.7000 # O(dmso)-S(dmso) LJ
pair_coeff 1 3 lj/cut/coul/long 0.0967 3.7400 # O(dmso)-C(dmso) LJ
pair_coeff 1 4 lj/cut/coul/long 0.0537 3.0400 # H(dmso)-O(dmso) LJ
pair_coeff 1 5 lj/cut/coul/long 0.1112 2.4212 # O(dmso)-Li+ LJ
pair_coeff 1 6 lj/cut/coul/long 0.1257 3.5475 # O(dmso)-P(pf6) LJ
pair_coeff 1 7 lj/cut/coul/long 0.0587 3.1674 # O(dmso)-F(pf6) LJ
pair_coeff 1 8 lj/cut/coul/long 0.0817 3.4000 # O(dmso)-C(ztc) LJ
pair_coeff 1 9 lj/cut/coul/long 0.1121 3.3485 # O(dmso)-O2 LJ
pair_coeff 1 10 lj/cut/coul/long 0.0000 1.7000 # O(dmso)-H(ztc) LJ

pair_coeff 2 2 lj/cut/coul/long 0.3500 4.0000 # S(dmso)-S(dmso) LJ
pair_coeff 2 3 lj/cut/coul/long 0.1652 4.0400 # S(dmso)-C(dmso) LJ
pair_coeff 2 4 lj/cut/coul/long 0.0916 3.3400 # S(dmso)-H(dmso) LJ
pair_coeff 2 5 lj/cut/coul/long 0.1899 2.7212 # S(dmso)-Li+ LJ
pair_coeff 2 6 lj/cut/coul/long 0.2147 3.8475 # S(dmso)-P(pf6) LJ
pair_coeff 2 7 lj/cut/coul/long 0.1002 3.4674 # S(dmso)-F(pf6) LJ
pair_coeff 2 8 lj/cut/coul/long 0.1395 3.7000 # S(dmso)-C (ztc) LJ
pair_coeff 2 9 lj/cut/coul/long 0.1914 3.6485 # S(dmso)-O2 LJ
pair_coeff 2 10 lj/cut/coul/long 0.0000 2.0000 # S(dmso)-H(ztc) LJ

pair_coeff 3 3 lj/cut/coul/long 0.0780 4.0800 # C(dmso)-C(dmso) LJ
pair_coeff 3 4 lj/cut/coul/long 0.0433 3.3800 # C(dmso)-H(dmso) LJ
pair_coeff 3 5 lj/cut/coul/long 0.0897 2.7612 # C(dmso)-Li+ LJ
pair_coeff 3 6 lj/cut/coul/long 0.1013 3.8875 # C(dmso)-P(pf6) LJ
pair_coeff 3 7 lj/cut/coul/long 0.0473 3.5074 # C(dmso)-F(pf6) LJ
pair_coeff 3 8 lj/cut/coul/long 0.0658 3.7400 # C(dmso)-C(ztc) LJ
pair_coeff 3 9 lj/cut/coul/long 0.0904 3.6885 # C(dmso)-O2 LJ
pair_coeff 3 10 lj/cut/coul/long 0.0000 2.0400 # C(dmso)-H(ztc) LJ

pair_coeff 4 4 lj/cut/coul/long 0.0240 2.6800 # H(dmso)-H(dmso) LJ
pair_coeff 4 5 lj/cut/coul/long 0.0498 2.0612 # H(dmso)-Li+ LJ
pair_coeff 4 6 lj/cut/coul/long 0.0562 3.1875 # H(dmso)-P(pf6) LJ
pair_coeff 4 7 lj/cut/coul/long 0.0262 2.8074 # H(dmso)-F(pf6) LJ
pair_coeff 4 8 lj/cut/coul/long 0.0365 3.0400 # H(dmso)-C(ztc) LJ
pair_coeff 4 9 lj/cut/coul/long 0.0501 2.9885 # H(dmso)-O2 LJ
pair_coeff 4 10 lj/cut/coul/long 0.0000 1.3400 # H(dmso)-H(ztc) LJ

#pair_coeff 5 5 lj/cut/coul/long 0.0000 1.4424 # Li±Li+ LJ
pair_coeff 5 6 lj/cut/coul/long 0.1165 2.5687 # Li±P(pf6) LJ
pair_coeff 5 7 lj/cut/coul/long 0.0544 2.1885 # Li±F(pf6) LJ
#pair_coeff 5 8 lj/cut/coul/long 0.0000 2.4212 # Li±C(ztc) LJ
#pair_coeff 5 9 lj/cut/coul/long 0.0000 2.3697 # Li±O2 LJ
pair_coeff 5 10 lj/cut/coul/long 0.0000 0.7212 # Li±H(ztc) LJ

pair_coeff 6 6 lj/cut/coul/long 0.1317 3.6950 # P(pf6)-P(pf6) LJ
pair_coeff 6 7 lj/cut/coul/long 0.0287 3.3148 # P(pf6)-F(pf6) LJ
pair_coeff 6 8 lj/cut/coul/long 0.0856 3.5475 # P(pf6)-C(ztc) LJ
pair_coeff 6 9 lj/cut/coul/long 0.1174 3.4960 # P(pf6)-O2 LJ
pair_coeff 6 10 lj/cut/coul/long 0.0000 1.8475 # P(pf6)-H(ztc) LJ

pair_coeff 7 7 lj/cut/coul/long 0.0287 2.9347 # F(pf6)-F(pf6) LJ
pair_coeff 7 8 lj/cut/coul/long 0.0399 3.1674 # F(pf6)-C(ztc) LJ
pair_coeff 7 9 lj/cut/coul/long 0.0548 3.1158 # F(pf6)-O2 LJ
pair_coeff 7 10 lj/cut/coul/long 0.0000 1.4674 # F(pf6)-H(ztc) LJ

#pair_coeff 8 8 lj/cut/coul/long 0.0000 3.4000 # C(ztc)-C(ztc) LJ
#pair_coeff 8 9 lj/cut/coul/long 0.0000 3.3485 # C(ztc)-O2 LJ
pair_coeff 8 10 lj/cut/coul/long 0.0000 1.7000 # C(ztc)-H (ztc) LJ

#pair_coeff 9 9 lj/cut/coul/long 0.0000 3.2970 # O2-O2 LJ
pair_coeff 9 10 lj/cut/coul/long 0.0000 1.6485 # O2-H(ztc) LJ

pair_coeff 10 10 lj/cut/coul/long 0.0000 0.0000 # H(ztc)-H(ztc) LJ

-------------------- Settings Section -------------------

neighbor 1 bin
neigh_modify every 1 delay 0 check no

timestep 0.01

fix 2 all qeq/reax 1 0.0 10.0 1e-4 reax/c
fix rffbonds all reax/c/bonds 100 Cathode.bonds

velocity dmso create 300 4928459 dist gaussian
velocity salt create 300 4928459 dist gaussian
velocity Li create 300 4928459 dist gaussian
velocity oxygen create 300 4928459 dist gaussian

-------------------- Run Section (NVT: density)------------------

fix 1 C rigid group 1 C reinit yes
fix C C setforce 0.0 0.0 0.0

fix 3 dmso nvt temp 300.0 300.0 10
fix 4 salt nvt temp 300.0 300.0 10
fix 5 Li nvt temp 300.0 300.0 10
fix 6 oxygen nvt temp 300.0 300.0 10

-------------------- MSD Section -------------------

compute msdelec dmso msd com yes
compute msdpf6 salt msd com yes
compute msdli Li msd com yes
compute msdo2 oxygen msd com yes

fix msdelec dmso ave/time 100 10 1000 c_msdelec[1] c_msdelec[2] c_msdelec[3] c_msdelec[4] file elec.msd
fix msdpf6 salt ave/time 100 10 1000 c_msdpf6[1] c_msdpf6[2] c_msdpf6[3] c_msdpf6[4] file pf6.msd
fix msdli Li ave/time 100 10 1000 c_msdli[1] c_msdli[2] c_msdli[3] c_msdli[4] file Li.msd
fix msdo2 oxygen ave/time 100 10 1000 c_msdo2[1] c_msdo2[2] c_msdo2[3] c_msdo2[4] file O2.msd

-------------------- Concentration Section -------------------

compute elec dmso chunk/atom bin/1d z lower 2.0
compute pf6 salt chunk/atom bin/1d z lower 2.0
compute lithium Li chunk/atom bin/1d z lower 2.0
compute o2 oxygen chunk/atom bin/1d z lower 2.0

fix 8 dmso ave/chunk 100 10 1000 elec density/mass density/number c_elec ave running file density_elec.out
fix 9 salt ave/chunk 100 10 1000 pf6 density/mass density/number c_pf6 ave running file density_pf6.out
fix 10 Li ave/chunk 100 10 1000 lithium density/mass density/number c_lithium ave running file density_lithium.out
fix 11 oxygen ave/chunk 100 10 1000 o2 density/mass density/number c_o2 ave running file density_o2.out

-------------------- Simulation data settings ------------------------

dump 1 all atom 100 Cathode.dump

dump_modify 1 scale no image yes

dump 2 all custom 100 Cathode.veldump vx vy vz
dump 3 all atom 1000 MD4-NVT*.lammpstrj

compute reax all pair reax/c
variable rffeb equal c_reax[1]
variable rffea equal c_reax[2]
variable rffelp equal c_reax[3]
variable rffemol equal c_reax[4]
variable rffev equal c_reax[5]
variable rffepen equal c_reax[6]
variable rffecoa equal c_reax[7]
variable rffehb equal c_reax[8]
variable rffet equal c_reax[9]
variable rffeco equal c_reax[10]
variable rffew equal c_reax[11]
variable rffep equal c_reax[12]
variable rffefi equal c_reax[13]
variable rffeqeq equal c_reax[14]

thermo 100

thermo_style custom etotal ke pe v_rffeb v_rffea v_rffelp v_rffev v_rffepen v_rffecoa v_rffet v_rffew v_rffep elong v_rffeqeq v_rffehb temp press

thermo_modify line multi

thermo_modify flush yes

restart 30000 MD.restart

run 40000

write_data Cathode-MD1.data

Regards
Khizar

4

This is bogus due to interaction double counting and incorrect application of QeQ to charges that are supposed to be static.

When using ReaxFF you usually need to have a parameter set for the entire system. There are a few special cases where you can expect reasonable results, but they require that you have separate objects and only the ReaxFF interactions may use charges.

This is not satisfied for your system:

  • you apply the ReaxFF charge equilibration to all atoms and for most of them you don’t have suitable parameters. So LAMMPS will use garbage values. Under normal circumstances, this should result in a crash.
  • if you had QeQ parameters, it would be bogus to do charge equilibration on the CHARMM interactions, since that is not a polarizable force field and the charges must be unchanged
  • if, however, you exclude those charges from charge equilibration you will also exclude them from ReaxFF and have it adapt to that and that is incorrect as well.
  • your CHARMM force field interactions require long-range coulomb interactions. you are using the ewald kspace style for that. However, that is counting long-range contributions double, since in ReaxFF those are already included

Bottom line, what you are trying to do cannot be made to work correctly due to fundamental incompatibility of the two force fields you are trying to combine.

5

Dear akohlmey,

Thank you for your valuable feedback. Sure, I try to fix this in light of your comments.

Regards

6

There is nothing to fix. You need to start over with a consistent choice of force field.