tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2802 from tc387/charge-regulation-update

Browse Source 
Added variable pH support and fixed Ewald self-energy updating
This commit is contained in:
Axel Kohlmeyer
2021-06-17 23:02:14 -04:00
committed by GitHub
parent ce79622897 6b56177e48
commit 6a07942428
3 changed files with 56 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/fix_charge_regulation.rst
View File

@ -21,7 +21,7 @@ Syntax
.. parsed-literal::
keyword = *pH*, *pKa*, *pKb*, *pIp*, *pIm*, *pKs*, *acid_type*, *base_type*, *lunit_nm*, *temp*, *tempfixid*, *nevery*, *nmc*, *xrd*, *seed*, *tag*, *group*, *onlysalt*, *pmcmoves*
*pH* value = pH of the solution
*pH* value = pH of the solution (can be specified as an equal-style variable)
*pKa* value = acid dissociation constant
*pKb* value = base dissociation constant
*pIp* value = chemical potential of free cations

52
src/MC/fix_charge_regulation.cpp
View File

@ -31,6 +31,7 @@
#include "force.h"
#include "group.h"
#include "improper.h"
#include "input.h"
#include "kspace.h"
#include "math_const.h"
#include "math_special.h"
@ -40,6 +41,7 @@
#include "pair.h"
#include "random_park.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
@ -50,6 +52,8 @@ using namespace FixConst;
using namespace MathConst;
using namespace MathSpecial;
enum{CONSTANT,EQUAL}; // parsing input variables
// large energy value used to signal overlap
#define MAXENERGYSIGNAL 1.0e100
#define MAXENERGYTEST 1.0e50
@ -62,7 +66,7 @@ FixChargeRegulation::FixChargeRegulation(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
ngroups(0), groupstrings(nullptr), ptype_ID(nullptr),
random_equal(nullptr), random_unequal(nullptr),
idftemp(nullptr)
idftemp(nullptr), pHstr(nullptr)
{
// Region restrictions not yet implemented ..
@ -130,6 +134,7 @@ FixChargeRegulation::~FixChargeRegulation() {
delete random_equal;
delete random_unequal;
delete[] pHstr;
delete[] idftemp;
if (group) {
@ -150,6 +155,14 @@ void FixChargeRegulation::init() {
int ipe = modify->find_compute("thermo_pe");
c_pe = modify->compute[ipe];
if (pHstr) {
pHvar = input->variable->find(pHstr);
if (pHvar < 0)
error->all(FLERR,"Variable name for fix charge/regulation does not exist");
if (input->variable->equalstyle(pHvar)) pHstyle = EQUAL;
else error->all(FLERR,"Variable for fix charge/regulation is invalid style");
}
if (atom->molecule_flag) {
int flag = 0;
@ -270,6 +283,9 @@ void FixChargeRegulation::pre_exchange() {
}
beta = 1.0 / (force->boltz * *target_temperature_tcp);
if (pHstyle == EQUAL)
pH = input->variable->compute_equal(pHvar);
// pre-compute powers
c10pH = pow(10.0,-pH); // dissociated ion (H+) activity
c10pKa = pow(10.0,-pKa); // acid dissociation constant
@ -378,6 +394,8 @@ void FixChargeRegulation::forward_acid() {
factor = nacid_neutral * vlocal_xrd * c10pKa * c10pI_plus /
(c10pH * (1 + nacid_charged) * (1 + npart_xrd2));
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
@ -444,6 +462,8 @@ void FixChargeRegulation::backward_acid() {
factor = (1 + nacid_neutral) * vlocal_xrd * c10pKa * c10pI_plus /
(c10pH * nacid_charged * npart_xrd);
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
@ -467,6 +487,8 @@ void FixChargeRegulation::backward_acid() {
atom->q[m2] = 1;
atom->mask[m2] = mask_tmp;
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
}
} else {
if (m1 >= 0) {
@ -511,6 +533,8 @@ void FixChargeRegulation::forward_base() {
(c10pOH * (1 + nbase_charged) * (1 + npart_xrd2));
m2 = insert_particle(anion_type, -1, reaction_distance, pos_all);
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
random_equal->uniform() < factor * exp(beta * (energy_before - energy_after))) {
@ -575,6 +599,8 @@ void FixChargeRegulation::backward_base() {
factor = (1 + nbase_neutral) * vlocal_xrd * c10pKb * c10pI_minus /
(c10pOH * nbase_charged * npart_xrd);
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
@ -598,6 +624,8 @@ void FixChargeRegulation::backward_base() {
atom->q[m2] = -1;
atom->mask[m2] = mask_tmp;
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
}
} else {
if (m1 >= 0) {
@ -672,6 +700,8 @@ void FixChargeRegulation::backward_ions() {
}
factor = volume_rx * volume_rx * c10pI_plus * c10pI_minus / (ncation * nanion);
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
random_equal->uniform() < (1.0 / factor) * exp(beta * (energy_before - energy_after))) {
@ -701,9 +731,6 @@ void FixChargeRegulation::backward_ions() {
atom->nlocal--;
}
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
} else {
energy_stored = energy_before;
@ -716,6 +743,8 @@ void FixChargeRegulation::backward_ions() {
atom->q[m2] = -1;
atom->mask[m2] = mask2_tmp;
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
}
}
}
@ -750,6 +779,8 @@ void FixChargeRegulation::forward_ions_multival() {
}
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST && random_equal->uniform() < factor * exp(beta * (energy_before - energy_after))) {
energy_stored = energy_after;
@ -845,6 +876,8 @@ void FixChargeRegulation::backward_ions_multival() {
// attempt deletion
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
double energy_after = energy_full();
if (energy_after < MAXENERGYTEST &&
random_equal->uniform() < (1.0 / factor) * exp(beta * (energy_before - energy_after))) {
@ -880,8 +913,6 @@ void FixChargeRegulation::backward_ions_multival() {
nanion -= salt_charge_ratio;
ncation--;
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
} else {
energy_stored = energy_before;
@ -893,6 +924,8 @@ void FixChargeRegulation::backward_ions_multival() {
atom->mask[mm[i]] = mask_tmp[i];
}
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
}
}
@ -1251,7 +1284,12 @@ void FixChargeRegulation::options(int narg, char **arg) {
} else if (strcmp(arg[iarg], "pH") == 0) {
if (iarg + 2 > narg)
error->all(FLERR, "Illegal fix charge/regulation command");
pH = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (strstr(arg[iarg + 1],"v_") == arg[iarg + 1]) {
pHstr = utils::strdup(&arg[iarg + 1][2]);
} else {
pH = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
pHstyle = CONSTANT;
}
iarg += 2;
} else if (strcmp(arg[iarg], "pIp") == 0) {
if (iarg + 2 > narg)

18
src/MC/fix_charge_regulation.h
View File

@ -61,9 +61,9 @@ class FixChargeRegulation : public Fix {
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 and equilibrium constant in log10 base
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]
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
@ -71,17 +71,17 @@ class FixChargeRegulation : public Fix {
double vlocal_xrd; // # local volume within xrd
bool
only_salt_flag; // true if performing only salt insertion/deletion, no acid/base dissociation.
bool add_tags_flag; // true if each inserted atom gets its unique atom tag
int groupbitall; // group bitmask for inserted atoms
int ngroups; // number of group-ids for inserted atoms
char **groupstrings; // list of group-ids for inserted atoms
bool add_tags_flag; // true if each inserted atom gets its unique atom tag
int groupbitall; // group bitmask for inserted atoms
int ngroups; // number of group-ids for inserted atoms
char **groupstrings; // list of group-ids for inserted atoms
// counters
unsigned long int nacid_attempts, nacid_successes, nbase_attempts, nbase_successes,
nsalt_attempts, nsalt_successes;
int nacid_neutral, nacid_charged, nbase_neutral, nbase_charged, ncation,
nanion; // particle type counts
int cr_nmax; // max number of local particles
nanion; // particle type counts
int cr_nmax; // max number of local particles
double reservoir_temperature;
double beta, sigma, volume,
volume_rx; // inverse temperature, speed, total volume, reacting volume
@ -97,6 +97,8 @@ class FixChargeRegulation : public Fix {
int acid_type, cation_type, base_type, anion_type; // reacting atom types
int reaction_distance_flag; // radial reaction restriction flag
double reaction_distance; // max radial distance from acid/base for ion insertion
int pHvar, pHstyle; // variable pH style
char *pHstr; // variable pH input parsing
class Pair *pair;
class Compute *c_pe; // energy compute pointer