Applied edits/optimizations suggested by Axel. Further simplifified/fixed MC acceptance equations, few clarifications to documentation.
This commit is contained in:
tc387
@ -1,22 +1,18 @@
|
||||
|
||||
.. Yuan documentation master file, created by
|
||||
sphinx-quickstart on Sat Jan 30 14:06:22 2021.
|
||||
You can adapt this file completely to your liking, but it should at least
|
||||
contain the root `toctree` directive.
|
||||
tc387: Multiple text additions/changes, Feb 2 2021
|
||||
.. index:: fix fix_charge_regulation
|
||||
.. index:: fix charge/regulation
|
||||
|
||||
fix_charge_regulation command
|
||||
fix charge/regulation command
|
||||
=============================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID charge_regulation cation_type anion_type keyword value(s)
|
||||
fix ID group-ID charge/regulation cation_type anion_type keyword value(s)
|
||||
|
||||
* ID, group-ID are documented in fix command
|
||||
* charge_regulation = style name of this fix command
|
||||
* charge/regulation = style name of this fix command
|
||||
* cation_type = atom type of free cations
|
||||
* anion_type = atom type of free anions
|
||||
|
||||
@ -51,9 +47,9 @@ Examples
|
||||
""""""""
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix chareg all charge_regulation 1 2 acid_type 3 base_type 4 pKa 5 pKb 7 lb 1.0 nevery 200 nexchange 200 seed 123 tempfixid fT
|
||||
fix chareg all charge/regulation 1 2 acid_type 3 base_type 4 pKa 5 pKb 7 lb 1.0 nevery 200 nexchange 200 seed 123 tempfixid fT
|
||||
|
||||
fix chareg all charge_regulation 1 2 pIp 3 pIm 3 tempfixid fT tag yes onlysalt yes 2 -1
|
||||
fix chareg all charge/regulation 1 2 pIp 3 pIm 3 onlysalt yes 2 -1 seed 123 tag yes temp 1.0
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
@ -75,7 +71,7 @@ An acid ionization reaction (:math:`\mathrm{A} \rightleftharpoons \mathrm{A}^-+\
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix acid_reaction all charge_regulation 2 3 acid_type 1 pH 7.0 pKa 5.0 pIp 7.0 pIm 7.0
|
||||
fix acid_reaction all charge/regulation 2 3 acid_type 1 pH 7.0 pKa 5.0 pIp 7.0 pIm 7.0
|
||||
|
||||
where the fix attempts to charge :math:`\mathrm{A}` (discharge :math:`\mathrm{A}^-`) to :math:`\mathrm{A}^-` (:math:`\mathrm{A}`) and insert (delete) a proton (atom type 2). Besides, the fix implements self-ionization reaction of water :math:`\emptyset \rightleftharpoons \mathrm{H}^++\mathrm{OH}^-`.
|
||||
However, this approach is highly inefficient at :math:`\mathrm{pH} \approx 7` when the concentration of both protons and hydroxyl ions is low, resulting in a relatively low acceptance rate of MC moves.
|
||||
@ -85,7 +81,7 @@ participate in ionization reactions, which can be easily achieved via
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix acid_reaction all charge_regulation 4 5 acid_type 1 pH 7.0 pKa 5.0 pIp 2.0 pIm 2.0
|
||||
fix acid_reaction all charge/regulation 4 5 acid_type 1 pH 7.0 pKa 5.0 pIp 2.0 pIm 2.0
|
||||
|
||||
where particles of atom type 4 and 5 are the salt cations and anions, both at chemical potential pI=2.0, see :ref:`(Curk1) <Curk1>` and :ref:`(Landsgesell) <Landsgesell>` for more details.
|
||||
|
||||
@ -94,14 +90,14 @@ where particles of atom type 4 and 5 are the salt cations and anions, both at ch
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix base_reaction all charge_regulation 2 3 base_type 6 pH 7.0 pKb 6.0 pIp 7.0 pIm 7.0
|
||||
fix base_reaction all charge/regulation 2 3 base_type 6 pH 7.0 pKb 6.0 pIp 7.0 pIm 7.0
|
||||
|
||||
where the fix will attempt to charge :math:`\mathrm{B}` (discharge :math:`\mathrm{B}^+`) to :math:`\mathrm{B}^+` (:math:`\mathrm{B}`) and insert (delete) a hydroxyl ion :math:`\mathrm{OH}^-` of atom type 3.
|
||||
If neither the acid or the base type is specified, for example,
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix salt_reaction all charge_regulation 4 5 pIp 2.0 pIm 2.0
|
||||
fix salt_reaction all charge/regulation 4 5 pIp 2.0 pIm 2.0
|
||||
|
||||
the fix simply inserts or deletes an ion pair of a free cation (atom type 4) and a free anion (atom type 5) as done in a conventional grand-canonical MC simulation.
|
||||
|
||||
@ -119,7 +115,7 @@ Note that LAMMPS implicitly assumes a constant number of particles (degrees of f
|
||||
|
||||
The chemical potential units (e.g. pH) are in the standard log10 representation assuming reference concentration :math:`\rho_0 = \mathrm{mol}/\mathrm{l}`.
|
||||
Therefore, to perform the internal unit conversion, the length (in nanometers) of the LAMMPS unit length
|
||||
must be specified via *lunit_nm* (default is set to the Bjerrum length in water at room temprature *lunit_nm* = 0.72nm). For example, in the dilute ideal solution limit, the concentration of free ions
|
||||
must be specified via *lunit_nm* (default is set to the Bjerrum length in water at room temprature *lunit_nm* = 0.71nm). For example, in the dilute ideal solution limit, the concentration of free ions
|
||||
will be :math:`c_\mathrm{I} = 10^{-\mathrm{pIp}}\mathrm{mol}/\mathrm{l}`.
|
||||
|
||||
The temperature used in MC acceptance probability is set by *temp*. This temperature should be the same as the temperature set by the molecular dynamics thermostat. For most purposes, it is probably best to use *tempfixid* keyword which dynamically sets the temperature equal to the chosen MD thermostat temperature, in the example above we assumed the thermostat fix-ID is *fT*. The inserted particles attain a random velocity corresponding to the specified temperature. Using *tempfixid* overrides any fixed temperature set by *temp*.
|
||||
@ -138,14 +134,14 @@ The *group* keyword can be used to add inserted particles to a specific group-ID
|
||||
|
||||
Output
|
||||
""""""
|
||||
This fix computes a global vector of length 8, which can be accessed by various output commands. The vector values are the following global cumulative quantities:
|
||||
This fix computes a global vector of length 8, which can be accessed by various output commands. The vector values are the following global quantities:
|
||||
|
||||
* 1 = cumulative MC attempts
|
||||
* 2 = cumulative MC successes
|
||||
* 3 = current # of neutral acid atoms
|
||||
* 4 = current # of -1 charged acid atoms
|
||||
* 5 = current # of neutral base atoms
|
||||
* 6 = current # of +1 charged acid atoms
|
||||
* 6 = current # of +1 charged base atoms
|
||||
* 7 = current # of free cations
|
||||
* 8 = current # of free anions
|
||||
|
||||
@ -157,6 +153,8 @@ See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
The :doc:`atom_style <atom_style>`, used must contain the charge property, for example, the style could be *charge* or *full*. Only usable for 3D simulations. Atoms specified as free ions cannot be part of rigid bodies or molecules and cannot have bonding interactions. The scheme is limited to integer charges, any atoms with non-integer charges will not be considered by the fix.
|
||||
|
||||
All interaction potentials used must be continuous, otherwise the MD integration and the particle exchange MC moves do not correspond to the same equilibrium ensemble. For example, if an lj/cut pair style is used, the LJ potential must be shifted so that it vanishes at the cutoff. This can be easily achieved using the :doc:`pair_modify <pair_modify>` command, i.e., by using: *pair_modify shift yes*.
|
||||
|
||||
Note: Region restrictions are not yet implemented.
|
||||
|
||||
Related commands
|
||||
@ -167,7 +165,7 @@ Related commands
|
||||
|
||||
Default
|
||||
"""""""
|
||||
pH = 7.0; pKa = 100.0; pKb = 100.0; pIp = 5.0; pIm = 5.0; pKs=14.0; acid_type = -1; base_type = -1; lunit_nm = 0.72; temp = 1.0; nevery = 100; nmc = 100; xrd = 0; seed = 2345; tag = no; onlysalt = no, pmcmoves = 0.33 0.33 0.33, group-ID = all
|
||||
pH = 7.0; pKa = 100.0; pKb = 100.0; pIp = 5.0; pIm = 5.0; pKs = 14.0; acid_type = -1; base_type = -1; lunit_nm = 0.71; temp = 1.0; nevery = 100; nmc = 100; xrd = 0; seed = 0; tag = no; onlysalt = no, pmcmoves = [1/3, 1/3, 1/3], group-ID = all
|
||||
|
||||
----------
|
||||
|
||||
@ -181,4 +179,4 @@ pH = 7.0; pKa = 100.0; pKb = 100.0; pIp = 5.0; pIm = 5.0; pKs=14.0; acid_type =
|
||||
|
||||
.. _Landsgesell:
|
||||
|
||||
**(Landsgesell)** J. Landsgesell, P. Hebbeker, O. Rud, R. Lunkad, P. Kosovan, and C. Holm, “Grand-reaction method for simulations of ionization equilibria coupled to ion partitioning,” Macromolecules 53, 3007–3020 (2020).
|
||||
**(Landsgesell)** J. Landsgesell, P. Hebbeker, O. Rud, R. Lunkad, P. Kosovan, and C. Holm, "Grand-reaction method for simulations of ionization equilibria coupled to ion partitioning", Macromolecules 53, 3007–3020 (2020).
|
||||
|
||||
@ -2,6 +2,7 @@ This directory has two input scripts that illustrates how to use fix charge_regu
|
||||
|
||||
The script `in_acid.chreg` sets up a simple weak acid electrolyte (pH=7,pKa=6,pI=3). Four different types of MC moves are implemented: acid protonation & de-protonation, and monovalent ion pair insertion and deletion. Note here we have grouped all free monovalent ions into a single type, a physically natural choice on the level of coarse-grained primitive electrolyte models, which increases the calculation performance but has no effects on thermodynamic observables. The variables such as pH, pKa, pI, and lb at the top of the input script can be adjusted to play with various simulation parameters. The cumulative MC attempted moves and cumulative number of accepted moves, as well as, current number of neutral and charged acid particles, neutral and charged base particles (in this example always 0), and the current number of free cations and anions in the system are printed in the `log_acid.lammps`.
|
||||
|
||||
The script `in_polymer.chreg` sets up a weak polyelectrolyte chain of N=80 beads. Each bead is a weak acid with pKa=5 and solution has pH=7 and salt pI=3. In this example, we choose to treat salt ions, protons, and hydroxyl ions separately, which results in 5 types of MC moves: acid [type 1] protonation & de-protonation (with protons [type 4] insertion & deletion), acid [type 1] protonation & de-protonation (with salt cation [type 2] insertion & deletion), water self-ionization (insertion and deletion of proton [type4] and hydroxyl ion [type 5] pair), insertion and deletion of monovalent salt pair [type 2 and type 3] , insertion and deletion
|
||||
Of a proton [type4] and salt anion [type 3]. The current number of neutral and charged acid particles, the current number of free salt cations and anions, and the current number of protons and hydroxyl ions are printed in the `log_polymer.lammps`.
|
||||
The script `in_polymer.chreg` sets up a weak polyelectrolyte chain of N=80 beads. Each bead is a weak acid with pKa=5 and solution has pH=7 and monovalent salt chemical potential pI=3. In this example, we choose to treat salt ions, protons, and hydroxyl ions separately, which results in 5 types of MC moves: acid [type 1] protonation & de-protonation (with protons [type 4] insertion & deletion), acid [type 1] protonation & de-protonation (with salt cation [type 2] insertion & deletion), water self-ionization (insertion and deletion of proton [type4] and hydroxyl ion [type 5] pair), insertion and deletion of monovalent salt pair [type 2 and type 3] , insertion and deletion
|
||||
Of a proton [type4] and salt anion [type 3].
|
||||
The current number of neutral and charged acid particles, the current number of free salt cations and anions, and the current number of protons and hydroxyl ions are printed in the `log_polymer.lammps`.
|
||||
|
||||
|
||||
@ -24,7 +24,7 @@ pair_style lj/cut/coul/long ${cut_lj} ${cut_long}
|
||||
pair_coeff * * 1.0 1.0 ${cut_lj} # charges
|
||||
kspace_style ewald 1.0e-3
|
||||
dielectric 1.0
|
||||
pair_modify shift no
|
||||
pair_modify shift yes
|
||||
|
||||
######### VERLET INTEGRATION WITH LANGEVIN THERMOSTAT ###########
|
||||
fix fnve all nve
|
||||
@ -33,8 +33,9 @@ compute_modify dtemp dynamic yes
|
||||
fix fT all langevin 1.0 1.0 1.0 123
|
||||
fix_modify fT temp dtemp
|
||||
|
||||
fix chareg all charge_regulation 2 3 acid_type 1 pH ${pH} pKa ${pKa} pIp ${pIp} pIm ${pIm} lunit_nm ${lunit_nm} nevery ${nevery} nmc ${nmc} seed 2345 tempfixid fT
|
||||
fix chareg all charge/regulation 2 3 acid_type 1 pH ${pH} pKa ${pKa} pIp ${pIp} pIm ${pIm} lunit_nm ${lunit_nm} nevery ${nevery} nmc ${nmc} seed 2345 tempfixid fT
|
||||
thermo 100
|
||||
thermo_style custom step pe c_dtemp f_chareg[1] f_chareg[2] f_chareg[3] f_chareg[4] f_chareg[5] f_chareg[6] f_chareg[7] f_chareg[8]
|
||||
log log_acid.lammps
|
||||
timestep 0.005
|
||||
run 10000
|
||||
run 20000
|
||||
|
||||
@ -14,7 +14,7 @@ velocity all create 1.0 8008 loop geom
|
||||
pair_style lj/cut/coul/long 1.122462 20
|
||||
pair_coeff * * 1.0 1.0 1.122462 # charges
|
||||
kspace_style pppm 1.0e-3
|
||||
pair_modify shift no
|
||||
pair_modify shift yes
|
||||
dielectric 1.0
|
||||
|
||||
######### VERLET INTEGRATION WITH LANGEVIN THERMOSTAT ###########
|
||||
@ -24,12 +24,14 @@ compute_modify dtemp dynamic yes
|
||||
fix fT all langevin 1.0 1.0 1.0 123
|
||||
fix_modify fT temp dtemp
|
||||
|
||||
fix chareg1 all charge_regulation 2 3 acid_type 1 pH 7.0 pKa 6.5 pIp 3.0 pIm 3.0 temp 1.0 nmc 40
|
||||
fix chareg2 all charge_regulation 4 5 acid_type 1 pH 7.0 pKa 6.5 pIp 7.0 pIm 7.0 temp 1.0 nmc 40
|
||||
fix chareg3 all charge_regulation 4 3 pIp 7.0 pIm 3.0 temp 1.0 nmc 20
|
||||
fix chareg1 all charge/regulation 2 3 acid_type 1 pH 7.0 pKa 6.5 pIp 3.0 pIm 3.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg2 all charge/regulation 4 5 acid_type 1 pH 7.0 pKa 6.5 pIp 7.0 pIm 7.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg3 all charge/regulation 4 3 pIp 7.0 pIm 3.0 temp 1.0 nmc 20 seed 2345
|
||||
|
||||
thermo 100
|
||||
# print: step, potential energy, temperature, neutral acids, charged acids, salt cations, salt anions, H+ ions, OH- ions
|
||||
thermo_style custom step pe c_dtemp f_chareg1[3] f_chareg1[4] f_chareg1[7] f_chareg1[8] f_chareg2[7] f_chareg2[8]
|
||||
log log_polymer.lammps
|
||||
|
||||
timestep 0.005
|
||||
run 10000
|
||||
run 20000
|
||||
|
||||
@ -15,50 +15,56 @@
|
||||
Contributing author: Tine Curk (tcurk5@gmail.com) and Jiaxing Yuan (yuanjiaxing123@hotmail.com)
|
||||
------------------------------------------------------------------------- */
|
||||
#include "fix_charge_regulation.h"
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
#include "angle.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "molecule.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "bond.h"
|
||||
#include "comm.h"
|
||||
#include "compute.h"
|
||||
#include "group.h"
|
||||
#include "domain.h"
|
||||
#include "region.h"
|
||||
#include "random_park.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "bond.h"
|
||||
#include "angle.h"
|
||||
#include "dihedral.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "improper.h"
|
||||
#include "kspace.h"
|
||||
#include "math_extra.h"
|
||||
#include "math_const.h"
|
||||
#include "math_extra.h"
|
||||
#include "math_special.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "modify.h"
|
||||
#include "molecule.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "random_park.h"
|
||||
#include "region.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
using namespace MathSpecial;
|
||||
|
||||
// large energy value used to signal overlap
|
||||
#define MAXENERGYSIGNAL 1.0e100
|
||||
#define MAXENERGYTEST 1.0e50
|
||||
#define small 0.0000001
|
||||
#define PI 3.1415926
|
||||
#define SMALL 0.0000001
|
||||
#define NA_RHO0 0.602214 // Avogadro's constant times reference concentration (N_A * mol / liter) [nm^-3]
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
Fix_charge_regulation::Fix_charge_regulation(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
ngroups(0), groupstrings(NULL),
|
||||
random_equal(NULL), random_unequal(NULL),
|
||||
idftemp(NULL), ptype_ID(NULL) {
|
||||
|
||||
FixChargeRegulation::FixChargeRegulation(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
ngroups(0), groupstrings(NULL),
|
||||
random_equal(NULL), random_unequal(NULL),
|
||||
idftemp(NULL), ptype_ID(NULL)
|
||||
{
|
||||
|
||||
// Region restrictions not yet implemented ..
|
||||
|
||||
@ -79,28 +85,28 @@ Fix_charge_regulation::Fix_charge_regulation(LAMMPS *lmp, int narg, char **arg)
|
||||
// set defaults and read optional arguments
|
||||
options(narg - 5, &arg[5]);
|
||||
|
||||
if (nevery <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (nmc < 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (llength_unit_in_nm < 0.0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (*target_temperature_tcp < 0.0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (seed <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (cation_type <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (anion_type <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (reaction_distance < 0.0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (salt_charge[0] <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (salt_charge[1] >= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (nevery <= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (nmc < 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (llength_unit_in_nm < 0.0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (*target_temperature_tcp < 0.0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (seed <= 0) error->all(FLERR, "Illegal fix charge/regulation command: Seed value (positive integer) must be provided ");
|
||||
if (cation_type <= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (anion_type <= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (reaction_distance < 0.0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (salt_charge[0] <= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (salt_charge[1] >= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if ((salt_charge[1] % salt_charge[0] != 0) && (salt_charge[0] % salt_charge[1] != 0))
|
||||
error->all(FLERR,
|
||||
"Illegal fix charge_regulation command, multivalent cation/anion charges are allowed, "
|
||||
"Illegal fix charge/regulation command, multivalent cation/anion charges are allowed, "
|
||||
"but must be divisible, e.g. (3,-1) is fine, but (3,-2) is not implemented");
|
||||
|
||||
if (pmcmoves[0] < 0 || pmcmoves[1] < 0 || pmcmoves[2] < 0)
|
||||
error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
if (acid_type < 0) pmcmoves[0] = 0;
|
||||
if (base_type < 0) pmcmoves[1] = 0;
|
||||
// normalize
|
||||
double psum = pmcmoves[0] + pmcmoves[1] + pmcmoves[2];
|
||||
if (psum <= 0) error->all(FLERR, "Illegal fix charge_regulation command");
|
||||
if (psum <= 0) error->all(FLERR, "Illegal fix charge/regulation command");
|
||||
pmcmoves[0] /= psum;
|
||||
pmcmoves[1] /= psum;
|
||||
pmcmoves[2] /= psum;
|
||||
@ -117,7 +123,7 @@ Fix_charge_regulation::Fix_charge_regulation(LAMMPS *lmp, int narg, char **arg)
|
||||
nsalt_successes = 0;
|
||||
}
|
||||
|
||||
Fix_charge_regulation::~Fix_charge_regulation() {
|
||||
FixChargeRegulation::~FixChargeRegulation() {
|
||||
|
||||
memory->destroy(ptype_ID);
|
||||
|
||||
@ -130,18 +136,16 @@ Fix_charge_regulation::~Fix_charge_regulation() {
|
||||
}
|
||||
}
|
||||
|
||||
int Fix_charge_regulation::setmask() {
|
||||
int FixChargeRegulation::setmask() {
|
||||
int mask = 0;
|
||||
mask |= PRE_EXCHANGE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::init() {
|
||||
void FixChargeRegulation::init() {
|
||||
|
||||
triclinic = domain->triclinic;
|
||||
|
||||
char *id_pe = (char *) "thermo_pe";
|
||||
int ipe = modify->find_compute(id_pe);
|
||||
triclinic = domain->triclinic;
|
||||
int ipe = modify->find_compute("thermo_pe");
|
||||
c_pe = modify->compute[ipe];
|
||||
|
||||
if (atom->molecule_flag) {
|
||||
@ -165,34 +169,27 @@ void Fix_charge_regulation::init() {
|
||||
// skip if already exists from previous init()
|
||||
|
||||
if (!exclusion_group_bit) {
|
||||
char **group_arg = new char *[4];
|
||||
|
||||
// create unique group name for atoms to be excluded
|
||||
int len = strlen(id) + 30;
|
||||
group_arg[0] = new char[len];
|
||||
sprintf(group_arg[0], "Fix_CR:exclusion_group:%s", id);
|
||||
group_arg[1] = (char *) "subtract";
|
||||
group_arg[2] = (char *) "all";
|
||||
group_arg[3] = (char *) "all";
|
||||
group->assign(4, group_arg);
|
||||
exclusion_group = group->find(group_arg[0]);
|
||||
|
||||
auto group_id = std::string("FixChargeRegulation:CR_exclusion_group:") + id;
|
||||
group->assign(group_id + " subtract all all");
|
||||
exclusion_group = group->find(group_id);
|
||||
if (exclusion_group == -1)
|
||||
error->all(FLERR, "Could not find fix CR exclusion group ID");
|
||||
error->all(FLERR,"Could not find fix charge/regulation exclusion group ID");
|
||||
exclusion_group_bit = group->bitmask[exclusion_group];
|
||||
|
||||
// neighbor list exclusion setup
|
||||
// turn off interactions between group all and the exclusion group
|
||||
|
||||
int narg = 4;
|
||||
char **arg = new char *[narg];;
|
||||
char **arg = new char*[narg];;
|
||||
arg[0] = (char *) "exclude";
|
||||
arg[1] = (char *) "group";
|
||||
arg[2] = group_arg[0];
|
||||
arg[2] = (char *) group_id.c_str();
|
||||
arg[3] = (char *) "all";
|
||||
neighbor->modify_params(narg, arg);
|
||||
delete[] group_arg[0];
|
||||
delete[] group_arg;
|
||||
delete[] arg;
|
||||
neighbor->modify_params(narg,arg);
|
||||
delete [] arg;
|
||||
}
|
||||
|
||||
// check that no deletable atoms are in atom->firstgroup
|
||||
@ -226,7 +223,7 @@ void Fix_charge_regulation::init() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::pre_exchange() {
|
||||
void FixChargeRegulation::pre_exchange() {
|
||||
|
||||
if (next_reneighbor != update->ntimestep) return;
|
||||
xlo = domain->boxlo[0];
|
||||
@ -263,14 +260,22 @@ void Fix_charge_regulation::pre_exchange() {
|
||||
reaction_distance = 0;
|
||||
}
|
||||
// volume in units of (N_A * mol / liter)
|
||||
volume_rx = (xhi - xlo) * (yhi - ylo) * (zhi - zlo) * pow(llength_unit_in_nm, 3) * 0.602214;
|
||||
if (reaction_distance < small) {
|
||||
volume_rx = (xhi - xlo) * (yhi - ylo) * (zhi - zlo) * cube(llength_unit_in_nm) * NA_RHO0;
|
||||
if (reaction_distance < SMALL) {
|
||||
vlocal_xrd = volume_rx;
|
||||
} else {
|
||||
vlocal_xrd = 4.0 * PI * pow(reaction_distance, 3) / 3.0 * pow(llength_unit_in_nm, 3) * 0.602214;
|
||||
vlocal_xrd = 4.0 * MY_PI * cube(reaction_distance) / 3.0 * cube(llength_unit_in_nm) * NA_RHO0;
|
||||
}
|
||||
beta = 1.0 / (force->boltz * *target_temperature_tcp);
|
||||
|
||||
// pre-compute powers
|
||||
c10pH = pow(10.0,-pH); // dissociated ion (H+) activity
|
||||
c10pKa = pow(10.0,-pKa); // acid dissociation constant
|
||||
c10pKb = pow(10.0,-pKb); // base dissociation constant
|
||||
c10pOH = pow(10.0,-pKs + pH); // dissociated anion (OH-) activity
|
||||
c10pI_plus = pow(10.0,-pI_plus); // free cation activity
|
||||
c10pI_minus = pow(10.0,-pI_minus); // free anion activity
|
||||
|
||||
// reinitialize counters
|
||||
nacid_neutral = particle_number(acid_type, 0);
|
||||
nacid_charged = particle_number(acid_type, -1);
|
||||
@ -337,7 +342,7 @@ void Fix_charge_regulation::pre_exchange() {
|
||||
next_reneighbor = update->ntimestep + nevery;
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::forward_acid() {
|
||||
void FixChargeRegulation::forward_acid() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
@ -360,7 +365,7 @@ void Fix_charge_regulation::forward_acid() {
|
||||
pos[2] = atom->x[m1][2];
|
||||
}
|
||||
npart_xrd2 = ncation;
|
||||
if (reaction_distance >= small) {
|
||||
if (reaction_distance >= SMALL) {
|
||||
pos_all[0] = pos[0];
|
||||
pos_all[1] = pos[1];
|
||||
pos_all[2] = pos[2];
|
||||
@ -368,8 +373,8 @@ void Fix_charge_regulation::forward_acid() {
|
||||
npart_xrd2 = particle_number_xrd(cation_type, 1, reaction_distance, pos_all);
|
||||
}
|
||||
m2 = insert_particle(cation_type, 1, reaction_distance, pos_all);
|
||||
factor = nacid_neutral * vlocal_xrd * pow(10, -pKa)
|
||||
* (1 + pow(10, pH - pI_plus)) / ((1 + nacid_charged) * (1 + npart_xrd2));
|
||||
factor = nacid_neutral * vlocal_xrd * c10pKa * c10pI_plus /
|
||||
(c10pH * (1 + nacid_charged) * (1 + npart_xrd2));
|
||||
|
||||
double energy_after = energy_full();
|
||||
|
||||
@ -395,7 +400,7 @@ void Fix_charge_regulation::forward_acid() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::backward_acid() {
|
||||
void FixChargeRegulation::backward_acid() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
@ -418,7 +423,7 @@ void Fix_charge_regulation::backward_acid() {
|
||||
pos[1] = atom->x[m1][1];
|
||||
pos[2] = atom->x[m1][2];
|
||||
}
|
||||
if (reaction_distance >= small) {
|
||||
if (reaction_distance >= SMALL) {
|
||||
pos_all[0] = pos[0];
|
||||
pos_all[1] = pos[1];
|
||||
pos_all[2] = pos[2];
|
||||
@ -433,8 +438,8 @@ void Fix_charge_regulation::backward_acid() {
|
||||
mask_tmp = atom->mask[m2]; // remember group bits.
|
||||
atom->mask[m2] = exclusion_group_bit;
|
||||
}
|
||||
factor = (1 + nacid_neutral) * vlocal_xrd * pow(10, -pKa)
|
||||
* (1 + pow(10, pH - pI_plus)) / (nacid_charged * npart_xrd);
|
||||
factor = (1 + nacid_neutral) * vlocal_xrd * c10pKa * c10pI_plus /
|
||||
(c10pH * nacid_charged * npart_xrd);
|
||||
|
||||
double energy_after = energy_full();
|
||||
|
||||
@ -468,7 +473,7 @@ void Fix_charge_regulation::backward_acid() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::forward_base() {
|
||||
void FixChargeRegulation::forward_base() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
@ -491,16 +496,15 @@ void Fix_charge_regulation::forward_base() {
|
||||
pos[2] = atom->x[m1][2];
|
||||
}
|
||||
npart_xrd2 = nanion;
|
||||
if (reaction_distance >= small) {
|
||||
if (reaction_distance >= SMALL) {
|
||||
pos_all[0] = pos[0];
|
||||
pos_all[1] = pos[1];
|
||||
pos_all[2] = pos[2];
|
||||
MPI_Allreduce(pos, pos_all, 3, MPI_DOUBLE, MPI_SUM, world);
|
||||
npart_xrd2 = particle_number_xrd(anion_type, -1, reaction_distance, pos_all);
|
||||
}
|
||||
factor = nbase_neutral * vlocal_xrd * pow(10, -pKb)
|
||||
* (1 + pow(10, pKs - pH - pI_minus)) /
|
||||
((1 + nbase_charged) * (1 + npart_xrd2));
|
||||
factor = nbase_neutral * vlocal_xrd * c10pKb * c10pI_minus /
|
||||
(c10pOH * (1 + nbase_charged) * (1 + npart_xrd2));
|
||||
m2 = insert_particle(anion_type, -1, reaction_distance, pos_all);
|
||||
|
||||
double energy_after = energy_full();
|
||||
@ -526,7 +530,7 @@ void Fix_charge_regulation::forward_base() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::backward_base() {
|
||||
void FixChargeRegulation::backward_base() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
@ -549,7 +553,7 @@ void Fix_charge_regulation::backward_base() {
|
||||
pos[1] = atom->x[m1][1];
|
||||
pos[2] = atom->x[m1][2];
|
||||
}
|
||||
if (reaction_distance >= small) {
|
||||
if (reaction_distance >= SMALL) {
|
||||
pos_all[0] = pos[0];
|
||||
pos_all[1] = pos[1];
|
||||
pos_all[2] = pos[2];
|
||||
@ -563,8 +567,8 @@ void Fix_charge_regulation::backward_base() {
|
||||
mask_tmp = atom->mask[m2]; // remember group bits.
|
||||
atom->mask[m2] = exclusion_group_bit;
|
||||
}
|
||||
factor = (1 + nbase_neutral) * vlocal_xrd * pow(10, -pKb)
|
||||
* (1 + pow(10, pKs - pH - pI_minus)) / (nbase_charged * npart_xrd);
|
||||
factor = (1 + nbase_neutral) * vlocal_xrd * c10pKb * c10pI_minus /
|
||||
(c10pOH * nbase_charged * npart_xrd);
|
||||
|
||||
double energy_after = energy_full();
|
||||
|
||||
@ -598,20 +602,20 @@ void Fix_charge_regulation::backward_base() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::forward_ions() {
|
||||
void FixChargeRegulation::forward_ions() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
double *dummyp;
|
||||
int m1 = -1, m2 = -1;
|
||||
factor = volume_rx * volume_rx * (pow(10, -pH) + pow(10, -pI_plus))
|
||||
* (pow(10, -pKs + pH) + pow(10, -pI_minus)) /
|
||||
factor = volume_rx * volume_rx * c10pI_plus * c10pI_minus /
|
||||
((1 + ncation) * (1 + nanion));
|
||||
|
||||
m1 = insert_particle(cation_type, +1, 0, dummyp);
|
||||
m2 = insert_particle(anion_type, -1, 0, dummyp);
|
||||
double energy_after = energy_full();
|
||||
if (energy_after < MAXENERGYTEST && random_equal->uniform() < factor * exp(beta * (energy_before - energy_after))) {
|
||||
if (energy_after < MAXENERGYTEST &&
|
||||
random_equal->uniform() < factor * exp(beta * (energy_before - energy_after))) {
|
||||
energy_stored = energy_after;
|
||||
nsalt_successes += 1;
|
||||
ncation++;
|
||||
@ -632,7 +636,7 @@ void Fix_charge_regulation::forward_ions() {
|
||||
}
|
||||
|
||||
|
||||
void Fix_charge_regulation::backward_ions() {
|
||||
void FixChargeRegulation::backward_ions() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
@ -659,8 +663,7 @@ void Fix_charge_regulation::backward_ions() {
|
||||
mask2_tmp = atom->mask[m2];
|
||||
atom->mask[m2] = exclusion_group_bit;
|
||||
}
|
||||
factor = (volume_rx * volume_rx * (pow(10, -pH) + pow(10, -pI_plus)) *
|
||||
(pow(10, -pKs + pH) + pow(10, -pI_minus))) / (ncation * nanion);
|
||||
factor = volume_rx * volume_rx * c10pI_plus * c10pI_minus / (ncation * nanion);
|
||||
|
||||
double energy_after = energy_full();
|
||||
if (energy_after < MAXENERGYTEST &&
|
||||
@ -717,7 +720,7 @@ void Fix_charge_regulation::backward_ions() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::forward_ions_multival() {
|
||||
void FixChargeRegulation::forward_ions_multival() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor = 1;
|
||||
@ -728,19 +731,19 @@ void Fix_charge_regulation::forward_ions_multival() {
|
||||
// insert one anion and (salt_charge_ratio) cations
|
||||
|
||||
mm[0] = insert_particle(anion_type, salt_charge[1], 0, dummyp);
|
||||
factor *= volume_rx * pow(10, -pI_minus) / (1 + nanion);
|
||||
factor *= volume_rx * c10pI_minus / (1 + nanion);
|
||||
for (int i = 0; i < salt_charge_ratio; i++) {
|
||||
mm[i + 1] = insert_particle(cation_type, salt_charge[0], 0, dummyp);
|
||||
factor *= volume_rx * pow(10, -pI_plus) / (1 + ncation + i);
|
||||
factor *= volume_rx *c10pI_plus / (1 + ncation + i);
|
||||
}
|
||||
} else {
|
||||
// insert one cation and (salt_charge_ratio) anions
|
||||
|
||||
mm[0] = insert_particle(cation_type, salt_charge[0], 0, dummyp);
|
||||
factor *= volume_rx * pow(10, -pI_plus) / (1 + ncation);
|
||||
factor *= volume_rx * c10pI_plus / (1 + ncation);
|
||||
for (int i = 0; i < salt_charge_ratio; i++) {
|
||||
mm[i + 1] = insert_particle(anion_type, salt_charge[1], 0, dummyp);
|
||||
factor *= volume_rx * pow(10, -pI_minus) / (1 + nanion + i);
|
||||
factor *= volume_rx * c10pI_minus / (1 + nanion + i);
|
||||
}
|
||||
}
|
||||
|
||||
@ -771,7 +774,7 @@ void Fix_charge_regulation::forward_ions_multival() {
|
||||
}
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::backward_ions_multival() {
|
||||
void FixChargeRegulation::backward_ions_multival() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor = 1;
|
||||
@ -786,7 +789,7 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
|
||||
mm[0] = get_random_particle(anion_type, salt_charge[1], 0, dummyp);
|
||||
if (npart_xrd != nanion) error->all(FLERR, "Fix charge regulation salt count inconsistent");
|
||||
factor *= volume_rx * pow(10, -pI_minus) / (nanion);
|
||||
factor *= volume_rx * c10pI_minus / (nanion);
|
||||
if (mm[0] >= 0) {
|
||||
qq[0] = atom->q[mm[0]];
|
||||
atom->q[mm[0]] = 0;
|
||||
@ -796,7 +799,7 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
for (int i = 0; i < salt_charge_ratio; i++) {
|
||||
mm[i + 1] = get_random_particle(cation_type, salt_charge[0], 0, dummyp);
|
||||
if (npart_xrd != ncation - i) error->all(FLERR, "Fix charge regulation salt count inconsistent");
|
||||
factor *= volume_rx * pow(10, -pI_plus) / (ncation - i);
|
||||
factor *= volume_rx * c10pI_plus / (ncation - i);
|
||||
if (mm[i + 1] >= 0) {
|
||||
qq[i + 1] = atom->q[mm[i + 1]];
|
||||
atom->q[mm[i + 1]] = 0;
|
||||
@ -810,7 +813,7 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
if (nanion < salt_charge_ratio || ncation < 1) return;
|
||||
mm[0] = get_random_particle(cation_type, salt_charge[0], 0, dummyp);
|
||||
if (npart_xrd != ncation) error->all(FLERR, "Fix charge regulation salt count inconsistent");
|
||||
factor *= volume_rx * pow(10, -pI_plus) / (ncation);
|
||||
factor *= volume_rx * c10pI_plus / (ncation);
|
||||
if (mm[0] >= 0) {
|
||||
qq[0] = atom->q[mm[0]];
|
||||
atom->q[mm[0]] = 0;
|
||||
@ -826,7 +829,7 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
mask_tmp[i + 1] = atom->mask[mm[i + 1]];
|
||||
atom->mask[mm[i + 1]] = exclusion_group_bit;
|
||||
}
|
||||
factor *= volume_rx * pow(10, -pI_minus) / (nanion - i);
|
||||
factor *= volume_rx * c10pI_minus / (nanion - i);
|
||||
}
|
||||
}
|
||||
|
||||
@ -840,7 +843,7 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
atom->natoms -= 1 + salt_charge_ratio;
|
||||
// ions must be deleted in order, otherwise index m could change upon the first deletion
|
||||
for (int i = 0; i < salt_charge_ratio + 1; i++) {
|
||||
// get max mm value, poor N^2 scaling, but charge ratio is a small number (2 or 3).
|
||||
// get max mm value, poor N^2 scaling, but charge ratio is a SMALL number (2 or 3).
|
||||
int maxmm = -1, jmaxm = -1;
|
||||
for (int j = 0; j < salt_charge_ratio + 1; j++) {
|
||||
if (mm[j] > maxmm) {
|
||||
@ -883,23 +886,42 @@ void Fix_charge_regulation::backward_ions_multival() {
|
||||
}
|
||||
}
|
||||
|
||||
int Fix_charge_regulation::insert_particle(int ptype, double charge, double rd, double *target) {
|
||||
int FixChargeRegulation::insert_particle(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
// insert a particle of type (ptype) with charge (charge) within distance (rd) of (target)
|
||||
|
||||
double coord[3];
|
||||
int m = -1;
|
||||
if (rd < small) {
|
||||
if (rd < SMALL) {
|
||||
coord[0] = xlo + random_equal->uniform() * (xhi - xlo);
|
||||
coord[1] = ylo + random_equal->uniform() * (yhi - ylo);
|
||||
coord[2] = zlo + random_equal->uniform() * (zhi - zlo);
|
||||
} else {
|
||||
double radius = reaction_distance * random_equal->uniform();
|
||||
double theta = random_equal->uniform() * PI;
|
||||
double phi = random_equal->uniform() * 2 * PI;
|
||||
coord[0] = target[0] + radius * sin(theta) * cos(phi);
|
||||
coord[1] = target[1] + radius * sin(theta) * sin(phi);
|
||||
coord[2] = target[2] + radius * cos(theta);
|
||||
// get a random point inside a sphere with radius rd
|
||||
// simple rejection sampling, probably the fastest method
|
||||
double dxx=1,dyy=1,dzz=1;
|
||||
while (dxx * dxx + dyy * dyy + dzz * dzz > 1.0) {
|
||||
dxx = 2 * random_equal->uniform() - 1.0;
|
||||
dyy = 2 * random_equal->uniform() - 1.0;
|
||||
dzz = 2 * random_equal->uniform() - 1.0;
|
||||
}
|
||||
coord[0] = target[0] + rd * dxx;
|
||||
coord[1] = target[1] + rd * dyy;
|
||||
coord[2] = target[2] + rd * dzz;
|
||||
|
||||
// Alternative way, but likely somewhat less efficient
|
||||
/*
|
||||
double radius = rd * pow(random_equal->uniform(), THIRD);
|
||||
double theta = acos(2 * random_equal->uniform() - 1);
|
||||
double phi = random_equal->uniform() * 2 * MY_PI;
|
||||
double sinphi = sin(phi);
|
||||
double cosphi = cos(phi);
|
||||
double sintheta = sin(theta);
|
||||
double costheta = cos(theta);
|
||||
coord[0] = target[0] + radius * sintheta * cosphi;
|
||||
coord[1] = target[1] + radius * sintheta * sinphi;
|
||||
coord[2] = target[2] + radius * costheta;
|
||||
*/
|
||||
coord[0] = coord[0] - floor(1.0 * (coord[0] - xlo) / (xhi - xlo)) * (xhi - xlo);
|
||||
coord[1] = coord[1] - floor(1.0 * (coord[1] - ylo) / (yhi - ylo)) * (yhi - ylo);
|
||||
coord[2] = coord[2] - floor(1.0 * (coord[2] - zlo) / (zhi - zlo)) * (zhi - zlo);
|
||||
@ -926,7 +948,7 @@ int Fix_charge_regulation::insert_particle(int ptype, double charge, double rd,
|
||||
return m;
|
||||
}
|
||||
|
||||
int Fix_charge_regulation::get_random_particle(int ptype, double charge, double rd, double *target) {
|
||||
int FixChargeRegulation::get_random_particle(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
// returns a randomly chosen particle of type (ptype) with charge (charge)
|
||||
// chosen among particles within distance (rd) of (target)
|
||||
@ -947,9 +969,9 @@ int Fix_charge_regulation::get_random_particle(int ptype, double charge, double
|
||||
count_global = 0;
|
||||
count_before = 0;
|
||||
|
||||
if (rd < small) { //reaction_distance < small: No geometry constraint on random particle choice
|
||||
if (rd < SMALL) { //reaction_distance < SMALL: No constraint on random particle choice
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < small &&
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < SMALL &&
|
||||
atom->mask[i] != exclusion_group_bit) {
|
||||
ptype_ID[count_local] = i;
|
||||
count_local++;
|
||||
@ -966,7 +988,7 @@ int Fix_charge_regulation::get_random_particle(int ptype, double charge, double
|
||||
dz -= static_cast<int>(1.0 * dz / (zhi - zlo) + 0.5) * (zhi - zlo);
|
||||
distance_check = dx * dx + dy * dy + dz * dz;
|
||||
if ((distance_check < rd * rd) && atom->type[i] == ptype &&
|
||||
fabs(atom->q[i] - charge) < small && atom->mask[i] != exclusion_group_bit) {
|
||||
fabs(atom->q[i] - charge) < SMALL && atom->mask[i] != exclusion_group_bit) {
|
||||
ptype_ID[count_local] = i;
|
||||
count_local++;
|
||||
}
|
||||
@ -990,7 +1012,7 @@ int Fix_charge_regulation::get_random_particle(int ptype, double charge, double
|
||||
return -1;
|
||||
}
|
||||
|
||||
double Fix_charge_regulation::energy_full() {
|
||||
double FixChargeRegulation::energy_full() {
|
||||
int imolecule;
|
||||
if (triclinic) domain->x2lamda(atom->nlocal);
|
||||
domain->pbc();
|
||||
@ -1050,12 +1072,12 @@ double Fix_charge_regulation::energy_full() {
|
||||
return total_energy;
|
||||
}
|
||||
|
||||
int Fix_charge_regulation::particle_number_xrd(int ptype, double charge, double rd, double *target) {
|
||||
int FixChargeRegulation::particle_number_xrd(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
int count = 0;
|
||||
if (rd < small) {
|
||||
if (rd < SMALL) {
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < small && atom->mask[i] != exclusion_group_bit)
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < SMALL && atom->mask[i] != exclusion_group_bit)
|
||||
count++;
|
||||
}
|
||||
} else {
|
||||
@ -1069,7 +1091,7 @@ int Fix_charge_regulation::particle_number_xrd(int ptype, double charge, double
|
||||
dz -= static_cast<int>(1.0 * dz / (zhi - zlo) + 0.5) * (zhi - zlo);
|
||||
distance_check = dx * dx + dy * dy + dz * dz;
|
||||
if ((distance_check < rd * rd) && atom->type[i] == ptype &&
|
||||
fabs(atom->q[i] - charge) < small && atom->mask[i] != exclusion_group_bit) {
|
||||
fabs(atom->q[i] - charge) < SMALL && atom->mask[i] != exclusion_group_bit) {
|
||||
count++;
|
||||
}
|
||||
}
|
||||
@ -1079,11 +1101,11 @@ int Fix_charge_regulation::particle_number_xrd(int ptype, double charge, double
|
||||
return count_sum;
|
||||
}
|
||||
|
||||
int Fix_charge_regulation::particle_number(int ptype, double charge) {
|
||||
int FixChargeRegulation::particle_number(int ptype, double charge) {
|
||||
|
||||
int count = 0;
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < small && atom->mask[i] != exclusion_group_bit)
|
||||
if (atom->type[i] == ptype && fabs(atom->q[i] - charge) < SMALL && atom->mask[i] != exclusion_group_bit)
|
||||
count = count + 1;
|
||||
}
|
||||
int count_sum = count;
|
||||
@ -1091,7 +1113,7 @@ int Fix_charge_regulation::particle_number(int ptype, double charge) {
|
||||
return count_sum;
|
||||
}
|
||||
|
||||
double Fix_charge_regulation::compute_vector(int n) {
|
||||
double FixChargeRegulation::compute_vector(int n) {
|
||||
double count_temp = 0;
|
||||
if (n == 0) {
|
||||
return nacid_attempts + nbase_attempts + nsalt_attempts;
|
||||
@ -1113,7 +1135,7 @@ double Fix_charge_regulation::compute_vector(int n) {
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::setThermoTemperaturePointer() {
|
||||
void FixChargeRegulation::setThermoTemperaturePointer() {
|
||||
int ifix = -1;
|
||||
ifix = modify->find_fix(idftemp);
|
||||
if (ifix == -1) {
|
||||
@ -1126,7 +1148,7 @@ void Fix_charge_regulation::setThermoTemperaturePointer() {
|
||||
|
||||
}
|
||||
|
||||
void Fix_charge_regulation::assign_tags() {
|
||||
void FixChargeRegulation::assign_tags() {
|
||||
// Assign tags to ions with zero tags
|
||||
if (atom->tag_enable) {
|
||||
tagint *tag = atom->tag;
|
||||
@ -1167,14 +1189,14 @@ void Fix_charge_regulation::assign_tags() {
|
||||
parse input options
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Fix_charge_regulation::options(int narg, char **arg) {
|
||||
void FixChargeRegulation::options(int narg, char **arg) {
|
||||
if (narg < 0) error->all(FLERR, "Illegal fix charge regulation command");
|
||||
|
||||
// defaults
|
||||
|
||||
pH = 7.0;
|
||||
pI_plus = 100;
|
||||
pI_minus = 100;
|
||||
pI_plus = 5;
|
||||
pI_minus = 5;
|
||||
acid_type = -1;
|
||||
base_type = -1;
|
||||
pKa = 100;
|
||||
@ -1182,12 +1204,12 @@ void Fix_charge_regulation::options(int narg, char **arg) {
|
||||
pKs = 14.0;
|
||||
nevery = 100;
|
||||
nmc = 100;
|
||||
pmcmoves[0] = pmcmoves[1] = pmcmoves[2] = 0.33;
|
||||
llength_unit_in_nm= 0.72;
|
||||
pmcmoves[0] = pmcmoves[1] = pmcmoves[2] = THIRD;
|
||||
llength_unit_in_nm= 0.71; // Default set to Bjerrum length in water at 20 degrees C [nm]
|
||||
|
||||
reservoir_temperature = 1.0;
|
||||
reaction_distance = 0;
|
||||
seed = 12345;
|
||||
seed = 0;
|
||||
target_temperature_tcp = &reservoir_temperature;
|
||||
add_tags_flag = false;
|
||||
only_salt_flag = false;
|
||||
@ -1297,9 +1319,9 @@ void Fix_charge_regulation::options(int narg, char **arg) {
|
||||
if (iarg + 4 > narg) error->all(FLERR, "Illegal fix charge regulation command");
|
||||
salt_charge[0] = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
|
||||
salt_charge[1] = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
|
||||
if (fabs(salt_charge[0] - utils::inumeric(FLERR, arg[iarg + 2], false, lmp)) > small)
|
||||
if (fabs(salt_charge[0] - utils::inumeric(FLERR, arg[iarg + 2], false, lmp)) > SMALL)
|
||||
error->all(FLERR, "Illegal fix charge regulation command, cation charge must be an integer");
|
||||
if (fabs(salt_charge[1] - utils::inumeric(FLERR, arg[iarg + 3], false, lmp)) > small)
|
||||
if (fabs(salt_charge[1] - utils::inumeric(FLERR, arg[iarg + 3], false, lmp)) > SMALL)
|
||||
error->all(FLERR, "Illegal fix charge regulation command, anion charge must be an integer");
|
||||
iarg += 4;
|
||||
} else if (strcmp(arg[iarg + 1], "no") == 0) {
|
||||
@ -1329,7 +1351,7 @@ void Fix_charge_regulation::options(int narg, char **arg) {
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double Fix_charge_regulation::memory_usage() {
|
||||
double FixChargeRegulation::memory_usage() {
|
||||
double bytes = cr_nmax * sizeof(int);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
@ -17,21 +17,21 @@
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(charge_regulation,Fix_charge_regulation)
|
||||
FixStyle(charge/regulation,FixChargeRegulation)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_charge_regulation_H
|
||||
#define LMP_FIX_charge_regulation_H
|
||||
#ifndef LMP_FIX_CHARGE_REGULATION_H
|
||||
#define LMP_FIX_CHARGE_REGULATION_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fix_charge_regulation : public Fix {
|
||||
class FixChargeRegulation : public Fix {
|
||||
public:
|
||||
Fix_charge_regulation(class LAMMPS *, int, char **);
|
||||
~Fix_charge_regulation();
|
||||
FixChargeRegulation(class LAMMPS *, int, char **);
|
||||
~FixChargeRegulation();
|
||||
int setmask();
|
||||
void init();
|
||||
void pre_exchange();
|
||||
@ -59,7 +59,8 @@ namespace LAMMPS_NS {
|
||||
int nevery, seed; // begin MC cycle every nevery MD timesteps, random seed
|
||||
int nmc; // MC move attempts per cycle
|
||||
double llength_unit_in_nm ; // LAMMPS unit of length in nm, needed since chemical potentials are in units of mol/l
|
||||
double pH, pKa, pKb, pKs, pI_plus, pI_minus; // chemical potentials
|
||||
double pH, pKa, pKb, pKs, pI_plus, pI_minus; // chemical potentials and equilibrium constant in log10 base
|
||||
double c10pH, c10pKa, c10pKb, c10pOH, c10pI_plus, c10pI_minus; // 10 raised to chemical potential value, in units of concentration [mol/liter]
|
||||
double pmcmoves[3]; // mc move attempt probability: acid, base, ion pair exchange
|
||||
double pmcc; // mc move cumulative attempt probability
|
||||
int npart_xrd; // # of particles (ions) within xrd
|
||||
@ -78,7 +79,7 @@ namespace LAMMPS_NS {
|
||||
double reservoir_temperature;
|
||||
double beta, sigma, volume, volume_rx; // inverse temperature, speed, total volume, reacting volume
|
||||
int salt_charge[2]; // charge of salt ions: [0] - cation, [1] - anion
|
||||
int salt_charge_ratio; // charge ration when using multivalent ion exchange
|
||||
int salt_charge_ratio; // charge ratio when using multivalent ion exchange
|
||||
double xlo, xhi, ylo, yhi, zlo, zhi; // box size
|
||||
double energy_stored; // full energy of old/current configuration
|
||||
int triclinic; // 0 = orthog box, 1 = triclinic
|
||||
|
||||
Reference in New Issue
Block a user