QEQ with COMB or ReaxFF: "net charge" feasible ?

LAMMPS LAMMPS Mailing List Mirror
1

Dear Ray,

I have read that thread, but still cannot get how it works.
So I would be grateful if you could enlighten me…

I am using COMB potential.
Hence I am using fix qeq/comb.
I have prepared a cell with all charges equal to zero except 1 atom with charge +1, just to check.
I run the simulation with 1 step of CG (with zero steps turns out that no charge equilibration is performed),
and I get
TotEng Step Fmax E_vdwl E_coul E_long q1 q2 q3 qtot
-960.2194 0 0.46281597 -968.92244 8.7030403 0 0.010416667 0 0 1
-1002.4015 1 0.45350848 -978.01715 -24.384306 0 0.15562117 -0.15626975 -0.094006096 -1.5376589e-14

where at the beginning the total charge is 1 and then turns into zero after the equilibration.

From this data I deduce that fix/comb is imposing a zero net charge at the beginning of the charge equilibration …

I am using the Lammps version of 16th of February 2016.

What command/settings do I have to use to impose a constant net global charge to the system in fix/comb ?

Or do I have to use another, more recent version of Lammps ?

Portions of my inputs:
pair_style comb3 polar_off
fix 2 all qeq/comb 1.0 1e-3 verlet file fq.out 1 (should I remove verlet ?)
min_style cg
minimize 1.0e-20 1.0e-20 1 1

With this thermo style:
group type1 type 1
compute charge1 type1 property/atom q
compute q1 type1 reduce ave c_charge1

group type2 type 2
compute charge2 type2 property/atom q
compute q2 type2 reduce ave c_charge2

group type3 type 3
compute charge3 type3 property/atom q
compute q3 type3 reduce ave c_charge3

variable qtot equal count(type1)*c_q1+count(type2)*c_q2+count(type3)*c_q3

thermo_style custom etotal step fmax evdwl ecoul elong c_q1 c_q2 c_q3 v_qtot

Many thanks in advance for your help

Bests,
Daniele

2

Daniele,

You must be using a customized fix qeq/comb since there is no keyword “verlet” in the LAMMPS distribution. Where did you get your source code from? I can only speak for the one shipped with LAMMPS, which behaved as expected. Below is an example I modified from /examples/comb/in.comb.HfO2, with either run or minimization the global net charge is preserved:

AMMPS (14 Apr 2016)

Reading data file …

triclinic box = (0 0 0) to (25.642 25.957 26.4845) with tilt (0 -4.46691 0)

1 by 1 by 1 MPI processor grid

reading atoms …

1500 atoms

500 atoms in group type1

1000 atoms in group type2

Setting atom values …

1500 settings made for charge

Setting atom values …

1 settings made for charge

Reading potential file ffield.comb with DATE: 2011-02-22

Pair COMB:

generating Coulomb integral lookup table …

element[1] = Hf, z = 0.679131

element[2] = O , z = 2.24307

will apply over-coordination correction …

WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:168)

Neighbor list info …

1 neighbor list requests

update every 1 steps, delay 0 steps, check yes

max neighbors/atom: 2000, page size: 100000

master list distance cutoff = 12.5

ghost atom cutoff = 12.5

binsize = 6.25, bins = 5 5 5

Setting up cg style minimization …

Unit style: metal

Memory usage per processor = 7.63908 Mbytes

Step q1 q2 qtot

0 0.0002 0 0.1

1 3.3578348 -1.6788174 0.1

Loop time of 5.7059 on 1 procs for 1 steps with 1500 atoms

99.9% CPU use with 1 MPI tasks x no OpenMP threads

I suggest you try again with the official LAMMPS version. My input script is attached for your reference.

in.comb.HfO2 (1.46 KB)

3

Dear Ray,

thank you for your reply.

Yes: indeed I am using a version optimized for the latest COMB3 potential.

But actually it turned out that this was not the problem.

The minimization step was the problem in my case.

In order to “fix the net charge” I have to put a “run 0” before the minimization (as in your script you put a “run 1”).

And then it works and keeps the net charge fixed during the minimize command.

Maybe this information can be useful to someone else.

Thanks to you and to Tao Lin for this information.

Bests,

Daniele

4

Sorry, Daniele, but it worked just as fine if I don’t have that run 1 or run 0 and go straight to minimization. There must be something that was changed in your customized fix qeq/comb.

In sum, fix qeq/comb in the official LAMMPS distribution conserves global net charge.