fix qeq/comb command

Syntax:

fix ID group-ID qeq/comb dtq precision keyword value ... 
ID, group-ID are documented in fix command
qeq/comb = style name of this fix command
dtq = time step of this fix command (attosecond unit)
precision = the convergence and threshold criterion of the maximum charge force
zero or more keyword/value pairs may be appended
keyword = nevery or optimize or verlet or maxloop or appq or file 
  nevery value = Nevery
    Nevery = perform this fix every Nevery time steps
  optimize args = none
  verlet args = none
  maxloop value = Maxloop
    Maxloop = maximum loop of the charge dynamics
  appq value = Appq
    Appq = initial net charge that is applied to the system (e unit)
  file value = filename and qthermo
    filename = name of file to write charge dynamics info to
    qthermo = write the charge dynamics info to file every qthermo time steps
Examples:

fix 1 surface qeq/comb 2.0 0.1 file fq.out 1000
fix cathode cathode qeq/comb 2.0 0.1 appq 2.0
fix anode anode qeq/comb 2.0 0.1 appq -2.0
fix electrolyte electrolyte qeq/comb 2.0 0.1 appq 1.0

Description:

Perform the charge dynamics, in which the charge equations of motion uses the second
Newton's law in a similar manner of the molecular dynamics. The charge force
(UNIT: V) on an atom is the difference between mean value of the electronegativity of
the system and the electronegativity of the atom, where the electronegativity
of an atom is defined as the derivative of the total energy over the charge of
the atom. The fictious charge mass charge in this fix is 0.006 (UNIT: eV*fs^2/e^2). Though it can be applied to any variable atomic charge models, this fix is currently limited to the COMB (Charge-Optimized Many-Body) potential as described in (COMB_1), (COMB_2) and (COMB_3).

Only charges on the atoms in the specified group are involved this charge dynamics.

{dtq} is the time step of the charge dynamics in this fix. Since the charge
mass is smaller than ions, the {dtq} should be small. Typically, the {dtq}
should be smaller than 5 attoseconds.

The {precision} argument is the threshold value of the maximum charge force of
the system. During charge optimization or equilibration (or QEq), the charge
is equilibrated until the maximum charge force is smaller than {precision}.
During velocity verlet integration, the charge optimization is enforced if the
maximum charge force is larger than {precision}. The charge force of a
system scales with the temperature of the system. The typical charge force of
a system at 300 K is about 0.2 V. In this case, it is a waste of time to
equilibrate the charges such that the maximum charge force is way smaller than
0.2 V. The recommended {precision} value is 1.0e-3 for below 50 K, 0.1 for
room temperature and 0.5 for 1000 K or above. 

During a run, this fix is performed every {Nevery} time steps. The {Nevery} is defaulted to 1.
 
The keywords, {optimize} and {verlet} are designed to different styles of the charge dynamics.

For the {optimize} keyword, the charge optimization is performed in an
iterative fashion until the charge forces on all atoms is smaller than
{precision}.  See "Rappe_and_Goddard"_#Rappe_and_Goddard and
"Rick_and_Stuart"_#Rick_and_Stuart for details.

For the {verlet} keyword, the charge is determined by the velocity
verlet integration over the charge equations of motion. The charge
optimization is enforced under the following two conditions:
1) if it is  the first step of each run;
2) if the maximum charge force is larger than {precision} or the charge temperature
is larger than 5 K.

The {maxloop} argument gives the maximum number of iterations that are
allowed for the charge optimization. After the first step of each run, the
maximum number of iterations allowed is the smaller value between {maxloop}
and 20/{dtq}.

The {appq} argument gives the initial net charge that is applied to the system to simulate the charged system. The last three examples give a sample settings for an electrochemical system with applied voltage (V) between cathode and anode. The resulting applied charges on electrodes in this fix command is calculated as Appq = C*V, where C is the capacitance of electrolyte and V is the applied voltage. In this case, the electrolyte is charged with 1 e, which can reflect the unbalanced charge of cations and anions in the electrolyte.

This fix checks the charge neutrality if the net charge drifts
MIN(1.0e-4,{precision}*0.01). At the same
time, the charge velocity and force on each atom is reset to zero. The charge dynamics implemented in this fix can not gurantee the charge neutrality of the system. It is vital important to maintain the
neutrality to prevent the total energy from drifting.  

WARNING: Due to the thermal coupling between the nuclear and electronic degrees of freedom, the total energy of a system with "fix_qeq/comb" is not strictly conserved. For a "NVE" ensemble, the total energy of a system will decrease until the temperature of the system decreases to 0 K. For a "NVT" ensemble, the system normally reaches a steady state after 20 ps. No energy drift is observed.  

It is also advised to check the momentum conservation of the system every 200 ps. See "fix_momentum"_fix_momentum.html for details or use keywords in "fix_langevin"_fix_langevin.html thermostatting to check the momentum  conservation over time.

If the {file} keyword is used, then information about the charge dynamics is
written to the specifed file every {qthermo} time step. The information of the charge dynamics includes
"time step", "mean value of the electronegativity", "maximum force on
charges", "temperature of charges", "total charge", "kinetic energy from
charges" and "electrostatical energy shift". The last term is defined as the negative
value of the product of the "net charge" and the "mean value of the
electronegativity", where the "net charge" is the 
difference of the "total charge" between this time step and the first step. 
Please note that the "mean value of the electronegativity" is equal to the 
negative value of the electrostatical chemical potential. Thus, the last term
can be regarded as the energy needed to add the "net charge" to the system.

This fix computes a global scalar  which can be accessed by various
"output commands"_Section_howto.html#howto_15. The scalar consists of the kinetic
energy from charges and the electrostatical energy shift. The scalar value calculated by this fix is never
normalized. Note that the real total energy of a system is the total energy
from the ions + the compute scalar value of this fix.

Default:

fix 1 surface qeq/comb 2.0 0.1 nevery 1 optimize maxloop 200 appq 0.0
no file output