More cleanup to fix polarize*, access to the per-atom variables from the Atom class, and updates to doc pages
This commit is contained in:
Trung Nguyen
104
doc/src/fix_polarize_bem_gmres.rst
Normal file
104
doc/src/fix_polarize_bem_gmres.rst
Normal file
@ -0,0 +1,104 @@
|
||||
.. index:: fix polarize/bem/gmres
|
||||
|
||||
fix polarize/bem/gmres command
|
||||
===================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID polarize/bem/gmres nevery tolerance ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* polarize/bem/gmres = style name of this fix command
|
||||
* Nevery = this fixed is invoked every this many timesteps
|
||||
* tolerance = the tolerance for the iterative solver to stop
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 2 all polarize/bem/gmres 5 0.0001
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The fix polarize/bem/gmres computes the induced charges
|
||||
at the interface between two impermeable media with
|
||||
different dielectric constants.
|
||||
|
||||
There are some example scripts for using this fix
|
||||
with LAMMPS in the examples/USER/dielectric directory.
|
||||
|
||||
----------
|
||||
|
||||
The induced charges of the atoms in the specified group, which are
|
||||
the vertices on the interface, are computed using the equation:
|
||||
|
||||
..math::
|
||||
|
||||
\sigma_b(\mathbf{s}) = \dfrac{1 - \bar{\epsilon}}{\bar{\epsilon}}
|
||||
\sigma_f(\mathbf{s}) - \epsilon_0 \dfrac{\Delta \epsilon}{\bar{\epsilon}}
|
||||
\mathbf{E}(\mathbf{s}) \cdot \mathbf{n}(\mathbf{s})
|
||||
|
||||
* :math:`\sigma_b` is the induced charge density at the interface vertex
|
||||
:math:`\mathbf{s}`.
|
||||
* :math:`\bar{\epsilon}` is the mean dielectric constant at the interface vertex:
|
||||
:math:`\bar{\epsilon} = (\epsilon_1 + \epsilon_2)/2`.
|
||||
* :math:`\Delta \epsilon` is the dielectric constant difference at the interface vertex:
|
||||
:math:`\Delta \epsilon = \epsilon_1 - \epsilon_2`
|
||||
* :math:`\sigma_f` is the free charge density at the interface vertex
|
||||
* :math:`\mathbf{E}(\mathbf{s})` is the electrical field at the vertex
|
||||
* :math:`\mathbf{n}(\mathbf{s})` is the unit normal vector at the vertex
|
||||
pointing from medium with :math:`\epsilon_2` to that with :math:`\epsilon_1`
|
||||
|
||||
|
||||
The fix polarize/bem/gmres employs the Generalized Minimum Residual (GMRES)
|
||||
as described in :ref:`(Barros) <Barros>` to solve for :math:`\sigma_b`.
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
The :doc:`fix_modify <fix_modify>`
|
||||
|
||||
This fix computes a global scalar which can be accessed by various
|
||||
:doc:`output commands <Howto_output>`. The scalar is the Colvars
|
||||
energy mentioned above. The scalar value calculated by this fix is
|
||||
"extensive".
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix is part of the USER-DIELECTRIC package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` doc page for more info.
|
||||
|
||||
There can only be one colvars fix active at a time. Since the interface
|
||||
communicates only the minimum amount of information and colvars module
|
||||
itself can handle an arbitrary number of collective variables, this is
|
||||
not a limitation of functionality.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix polarize/bem/icc <fix_polarize_bem_icc>`, :doc:`fix polarize/functional <fix_polarize_functional>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
None.
|
||||
|
||||
----------
|
||||
|
||||
.. _Barros:
|
||||
|
||||
**(Barros)** Barros, Sinkovits, Luijten, J. Chem. Phys, 140, 064903 (2014)
|
||||
|
||||
|
||||
.. _NguyenTD:
|
||||
|
||||
**(NguyenTD)** Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-19 (2019)
|
||||
|
||||
103
doc/src/fix_polarize_bem_icc.rst
Normal file
103
doc/src/fix_polarize_bem_icc.rst
Normal file
@ -0,0 +1,103 @@
|
||||
.. index:: fix polarize/bem/icc
|
||||
|
||||
fix polarize/bem/icc command
|
||||
===================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID polarize nevery tolerance ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* polarize/bem/icc = style name of this fix command
|
||||
* Nevery = this fixed is invoked every this many timesteps
|
||||
* tolerance = the tolerance for the iterative solver to stop
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 interface polarize/bem/icc 1 0.0001
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The fix polarize/bem/icc computes the induced charges
|
||||
at the interface between two impermeable media with
|
||||
different dielectric constants.
|
||||
|
||||
There are some example scripts for using this fix
|
||||
with LAMMPS in the examples/USER/dielectric directory.
|
||||
|
||||
----------
|
||||
|
||||
The induced charges of the atoms in the specified group, which are
|
||||
the vertices on the interface, are computed using the equation:
|
||||
|
||||
..math::
|
||||
|
||||
\sigma_b(\mathbf{s}) = \dfrac{1 - \bar{\epsilon}}{\bar{\epsilon}}
|
||||
\sigma_f(\mathbf{s}) - \epsilon_0 \dfrac{\Delta \epsilon}{\bar{\epsilon}}
|
||||
\mathbf{E}(\mathbf{s}) \cdot \mathbf{n}(\mathbf{s})
|
||||
|
||||
* :math:`\sigma_b` is the induced charge density at the interface vertex
|
||||
:math:`\mathbf{s}`.
|
||||
* :math:`\bar{\epsilon}` is the mean dielectric constant at the interface vertex:
|
||||
:math:`\bar{\epsilon} = (\epsilon_1 + \epsilon_2)/2`.
|
||||
* :math:`\Delta \epsilon` is the dielectric constant difference at the interface vertex:
|
||||
:math:`\Delta \epsilon = \epsilon_1 - \epsilon_2`
|
||||
* :math:`\sigma_f` is the free charge density at the interface vertex
|
||||
* :math:`\mathbf{E}(\mathbf{s})` is the electrical field at the vertex
|
||||
* :math:`\mathbf{n}(\mathbf{s})` is the unit normal vector at the vertex
|
||||
pointing from medium with :math:`\epsilon_2` to that with :math:`\epsilon_1`
|
||||
|
||||
|
||||
The fix polarize/bem/icc employs the successive overrelaxation algorithm
|
||||
as described in :ref:`(Tyagi) <Tyagi>` to solve for :math:`\sigma_b`.
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
The :doc:`fix_modify <fix_modify>`
|
||||
|
||||
This fix computes a global scalar which can be accessed by various
|
||||
:doc:`output commands <Howto_output>`. The scalar is the Colvars
|
||||
energy mentioned above. The scalar value calculated by this fix is
|
||||
"extensive".
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix is part of the USER-DIELECTRIC package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` doc page for more info.
|
||||
|
||||
There can only be one colvars fix active at a time. Since the interface
|
||||
communicates only the minimum amount of information and colvars module
|
||||
itself can handle an arbitrary number of collective variables, this is
|
||||
not a limitation of functionality.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix polarize/bem/icc <fix_polarize_bem_gmres>`, :doc:`fix polarize/functional <fix_polarize_functional>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
None.
|
||||
|
||||
----------
|
||||
|
||||
.. _Tyagi:
|
||||
|
||||
**(Tyagi)** Tyagi, Suzen, Sega, Barbosa, Kantorovich, Holm, J Chem Phys, 132, 154112 (2010)
|
||||
|
||||
.. _NguyenTD:
|
||||
|
||||
**(NguyenTD)** Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-19 (2019)
|
||||
|
||||
@ -1,8 +1,6 @@
|
||||
.. index:: fix polarize/bem/icc
|
||||
.. index:: fix polarize/bem/gmres
|
||||
.. index:: fix polarize/functional
|
||||
|
||||
fix polarize command
|
||||
fix polarize/functional command
|
||||
===================
|
||||
|
||||
Syntax
|
||||
@ -13,7 +11,7 @@ Syntax
|
||||
fix ID group-ID polarize nevery tolerance ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* polarize/bem/icc, polarize/bem/gmres, or polarize/functional = style name of this fix command
|
||||
* polarize/functional = style name of this fix command
|
||||
* Nevery = this fixed is invoked every this many timesteps
|
||||
* tolerance = the tolerance for the iterative solver to stop
|
||||
|
||||
@ -23,9 +21,7 @@ Examples
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all polarize/bem/icc 1 0.0001
|
||||
fix 2 all polarize/bem/gmres 5 0.0001
|
||||
fix 3 all polarize/bem/functional 1 0.0001
|
||||
fix 3 all polarize/functional 1 0.001
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
@ -39,12 +35,6 @@ examples/USER/dielectric directory.
|
||||
|
||||
----------
|
||||
|
||||
The charges of the atoms in the specified group will be computed by the solver.
|
||||
fix polarize bem/icc computes the induced charges at the boundary elements
|
||||
(i.e. interface vertices) using the successive overrelaxation as described
|
||||
in (Tyagi). fix polarize bem/gmres computes the induced charges at
|
||||
the interface vertices using the successive overrelaxation
|
||||
as described in (Barros).
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
@ -71,7 +61,7 @@ not a limitation of functionality.
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix smd <fix_smd>`, :doc:`fix spring <fix_spring>`,
|
||||
:doc:`fix polarize/bem/icc <fix_polarize_bem_icc>`, :doc:`fix polarize/functional <fix_polarize_bem_gmres>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
@ -80,7 +70,11 @@ None.
|
||||
|
||||
----------
|
||||
|
||||
.. _Jadhao:
|
||||
|
||||
**(Jadhao)** Jadhao, Solis, Olvera de la Cruz, J Chem Phys, 138, 054119 (2013)
|
||||
|
||||
.. _NguyenTD:
|
||||
|
||||
**(NguyenTD)** Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Mol. Phys., DOI:10.1016/j.cpc.2019.03.006
|
||||
**(NguyenTD)** Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-19 (2019)
|
||||
|
||||
@ -63,7 +63,7 @@ elif "${method} == icc"&
|
||||
"fix 3 interface polarize/bem/icc 1 1.0e-4 itr_max 50" &
|
||||
"fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &
|
||||
elif "${method} == dof" &
|
||||
"fix 3 interface polarize/functional" &
|
||||
"fix 3 interface polarize/functional 1 0.001" &
|
||||
"fix_modify 3 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &
|
||||
else &
|
||||
"print 'Unsupported polarization solver' "
|
||||
|
||||
@ -35,9 +35,9 @@ fix 1 ions nve
|
||||
if "${method} == gmres" then &
|
||||
"fix 3 interface polarize/bem/gmres 1 1.0e-4" &
|
||||
elif "${method} == icc"&
|
||||
"fix 3 interface polarize/bem/icc 1 1.0e-4 itr_max 50" &
|
||||
"fix 3 interface polarize/bem/icc 1 1.0e-4" &
|
||||
elif "${method} == dof" &
|
||||
"fix 3 interface polarize/functional" &
|
||||
"fix 3 interface polarize/functional 1 1.0e-4" &
|
||||
else &
|
||||
"print 'Unsupported method for polarization' "
|
||||
|
||||
|
||||
@ -40,7 +40,7 @@ if "${method} == gmres" then &
|
||||
elif "${method} == icc"&
|
||||
"fix 3 interface polarize/bem/icc 1 1.0e-4 itr_max 50" &
|
||||
elif "${method} == dof" &
|
||||
"fix 3 interface polarize/functional" &
|
||||
"fix 3 interface polarize/functional 1 0.0001" &
|
||||
else &
|
||||
"print 'Unsupported method for polarization' "
|
||||
|
||||
|
||||
@ -155,3 +155,69 @@ void AtomVecDielectric::unpack_restart_init(int ilocal)
|
||||
nspecial[ilocal][2] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
assign an index to named atom property and return index
|
||||
return -1 if name is unknown to this atom style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int AtomVecDielectric::property_atom(char *name)
|
||||
{
|
||||
if (strcmp(name,"area") == 0) return 0;
|
||||
if (strcmp(name,"ed") == 0) return 1;
|
||||
if (strcmp(name,"em") == 0) return 2;
|
||||
if (strcmp(name,"epsilon") == 0) return 3;
|
||||
if (strcmp(name,"curvature") == 0) return 4;
|
||||
if (strcmp(name,"q_unscaled") == 0) return 5;
|
||||
return -1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack per-atom data into buf for ComputePropertyAtom
|
||||
index maps to data specific to this atom style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AtomVecDielectric::pack_property_atom(int index, double *buf,
|
||||
int nvalues, int groupbit)
|
||||
{
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int n = 0;
|
||||
|
||||
if (index == 0) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = area[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
} else if (index == 1) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = ed[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
} else if (index == 2) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = em[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
} else if (index == 3) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = epsilon[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
} else if (index == 4) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = curvature[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
} else if (index == 5) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) buf[n] = q_unscaled[i];
|
||||
else buf[n] = 0.0;
|
||||
n += nvalues;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -453,7 +453,8 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
|
||||
|
||||
// the inner loop k = 1..(n-1)
|
||||
// build up the k-th Krylov space,
|
||||
// actually build the Q_k matrix of size n by k, whose the columns are k vectors v(1)...v(k)
|
||||
// actually build the Q_k matrix of size n by k,
|
||||
// whose the columns are k vectors v(1)...v(k)
|
||||
// remember that v[0] is computed from the updated residual as above
|
||||
|
||||
for (k = 1; k <= mr; k++) {
|
||||
@ -503,7 +504,8 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
|
||||
}
|
||||
|
||||
// if k is not the first iteration,
|
||||
// find the vector y that minimizes the norm of the residual using the least square method
|
||||
// find the vector y that minimizes the norm of the residual
|
||||
// using the least square method
|
||||
|
||||
if (k > 1) {
|
||||
|
||||
@ -547,7 +549,8 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
|
||||
#ifdef _POLARIZE_DEBUG
|
||||
if (comm->me == 0) {
|
||||
char message[256];
|
||||
sprintf(message, "itr = %d: k = %d, norm(r) = %g norm(b) = %g", itr, k, rho, normb);
|
||||
sprintf(message, "itr = %d: k = %d, norm(r) = %g norm(b) = %g",
|
||||
itr, k, rho, normb);
|
||||
error->warning(FLERR, message);
|
||||
}
|
||||
#endif
|
||||
@ -586,7 +589,8 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
|
||||
#ifdef _POLARIZE_DEBUG
|
||||
if (comm->me == 0) {
|
||||
char message[256];
|
||||
sprintf(message, "itr = %d: norm(r) = %g norm(b) = %g", itr, rho, normb);
|
||||
sprintf(message, "itr = %d: norm(r) = %g norm(b) = %g",
|
||||
itr, rho, normb);
|
||||
error->warning(FLERR, message);
|
||||
}
|
||||
#endif
|
||||
@ -739,7 +743,8 @@ void FixPolarizeBEMGMRES::update_residual(double* w, double* r, int n)
|
||||
}
|
||||
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / epsilon[i];
|
||||
double sigma_f = q_real[i] / area[i];
|
||||
buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] - epsilon0 * ed[i] * dot / (4*MY_PI);
|
||||
buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] -
|
||||
epsilon0 * ed[i] * dot / (4*MY_PI);
|
||||
}
|
||||
|
||||
MPI_Allreduce(buffer,r,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
|
||||
@ -843,8 +848,10 @@ int FixPolarizeBEMGMRES::modify_param(int narg, char **arg)
|
||||
int set_charge=0;
|
||||
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
if (strcmp(arg[iarg+4],"NULL") != 0) areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
|
||||
if (strcmp(arg[iarg+3],"NULL") != 0)
|
||||
epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
if (strcmp(arg[iarg+4],"NULL") != 0)
|
||||
areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
|
||||
if (strcmp(arg[iarg+5],"NULL") != 0) {
|
||||
qreali = utils::numeric(FLERR,arg[iarg+5],false,lmp);
|
||||
set_charge = 1;
|
||||
|
||||
@ -391,8 +391,10 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
|
||||
int set_charge=0;
|
||||
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
if (strcmp(arg[iarg+4],"NULL") != 0) areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
|
||||
if (strcmp(arg[iarg+3],"NULL") != 0)
|
||||
epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
if (strcmp(arg[iarg+4],"NULL") != 0)
|
||||
areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
|
||||
if (strcmp(arg[iarg+5],"NULL") != 0) {
|
||||
qreali = utils::numeric(FLERR,arg[iarg+5],false,lmp);
|
||||
set_charge = 1;
|
||||
|
||||
@ -22,11 +22,7 @@
|
||||
(Rww + Rww^T) w = q Rwq
|
||||
at every time step, the vector (q Rwq) is computed, and so
|
||||
w = [Rww + Rww^T)^(-1)] (q Rwq)
|
||||
NOTE: Oct 7, 2019: switch from matrix inversion to a conjugate gradient solver
|
||||
|
||||
Contributing author: Trung Nguyen (Northwestern)
|
||||
based on the full-matrix implementation by Honghao Li (Northwestern University)
|
||||
Reference: Jadhao, Solis, Olvera de la Cruz, J. Chem. Phys. 138, 054119, 2013
|
||||
NOTE: Oct 7, 2019: switch to using a conjugate gradient solver
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_polarize_functional.h"
|
||||
@ -93,11 +89,17 @@ FixPolarizeFunctional::FixPolarizeFunctional(LAMMPS *lmp, int narg, char **arg)
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_user_dielectric_package);
|
||||
|
||||
if (narg != 3) error->all(FLERR,"Illegal fix polarize/functional command");
|
||||
if (narg < 4) error->all(FLERR,"Illegal fix polarize/functional command");
|
||||
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Fix polarize/functional requires atom style dielectric");
|
||||
|
||||
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
|
||||
if (nevery < 0) error->all(FLERR,"Illegal fix polarize/functional command");
|
||||
|
||||
tolerance = EPSILON;
|
||||
if (narg == 5) tolerance = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
|
||||
comm_forward = 1;
|
||||
nmax = 0;
|
||||
allocated = 0;
|
||||
@ -362,6 +364,9 @@ void FixPolarizeFunctional::setup_pre_force(int vflag)
|
||||
|
||||
void FixPolarizeFunctional::pre_force(int)
|
||||
{
|
||||
if (nevery == 0) return;
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// solve for the induced charges
|
||||
|
||||
update_induced_charges();
|
||||
@ -631,14 +636,14 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
// invoke full neighbor list
|
||||
|
||||
int inum,jnum,*ilist,*jlist,*numneigh,**firstneigh;
|
||||
inum = list->inum; // number of entries in the neighbor list
|
||||
ilist = list->ilist; // ilist[ii] gives the atom index of the entry ii in the list
|
||||
numneigh = list->numneigh; // numneigh[i] gives the number of neighbors of local atom i
|
||||
firstneigh = list->firstneigh; // firstneigh[i] gives the pointer to the neighbors of i
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// calculate G1ww, gradG1ww, ndotG1ww
|
||||
// fill up buffer1 with local G1ww and buffer2 with local ndotGww
|
||||
// seperate into two loops to let the compiler optimize / or later vectorization
|
||||
// seperate into two loops to let the compiler optimize/or later vectorization
|
||||
|
||||
for (int i = 0; i < num_induced_charges; i++)
|
||||
for (int j = 0; j < num_induced_charges; j++) buffer1[i][j] = 0;
|
||||
@ -660,6 +665,7 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
for (int kk = 0; kk < jnum; kk++) {
|
||||
int k = jlist[kk] & NEIGHMASK;
|
||||
if (mask[k] & groupbit) {
|
||||
|
||||
// interface particles: k can be ghost atoms
|
||||
double delx = xtmp - x[k][0];
|
||||
double dely = ytmp - x[k][1];
|
||||
@ -667,28 +673,31 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
int mk = tag2mat[tag[k]];
|
||||
|
||||
//G1ww[mi][mk] = calculate_greens_ewald(delx, dely, delz);
|
||||
// G1ww[mi][mk] = calculate_greens_ewald(delx, dely, delz);
|
||||
buffer1[mi][mk] = calculate_greens_ewald(delx, dely, delz);
|
||||
|
||||
// gradG1ww is vector, directly change it in the function
|
||||
double gradG1ww[3];
|
||||
//calculate_grad_greens_ewald(gradG1ww[mi][mk], delx, dely, delz);
|
||||
calculate_grad_greens_ewald(gradG1ww, delx, dely, delz);
|
||||
|
||||
// use mu to store the normal vector of interface vertex
|
||||
//ndotGww[mi][mk] = MathExtra::dot3(norm[i], gradG1ww[mi][mk]) / (4*MY_PI);
|
||||
buffer2[mi][mk] = MathExtra::dot3(norm[i], gradG1ww) / (4*MY_PI);
|
||||
}
|
||||
}
|
||||
|
||||
// special treatment for the diagonal terms, even though in the above loop there is mk == mi
|
||||
//G1ww[mi][mi] = calculate_greens_ewald_self_vertex(area[i]);
|
||||
// special treatment for the diagonal terms,
|
||||
// even though in the above loop there is mk == mi
|
||||
|
||||
buffer1[mi][mi] = calculate_greens_ewald_self_vertex(area[i]);
|
||||
//ndotGww[mi][mi] = calculate_ndotgreens_ewald_self_vertex(area[i], curvature[i]) / (4*MY_PI);
|
||||
buffer2[mi][mi] = calculate_ndotgreens_ewald_self_vertex(area[i], curvature[i]) / (4*MY_PI);
|
||||
buffer2[mi][mi] = calculate_ndotgreens_ewald_self_vertex(area[i], curvature[i]) /
|
||||
(4*MY_PI);
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(buffer1[0],G1ww[0],num_induced_charges*num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(buffer2[0],ndotGww[0],num_induced_charges*num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(buffer1[0],G1ww[0],num_induced_charges*num_induced_charges,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(buffer2[0],ndotGww[0],num_induced_charges*num_induced_charges,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// calculate G2ww
|
||||
// fill up buffer1 with local G2ww
|
||||
@ -738,10 +747,12 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(buffer1[0],G2ww[0],num_induced_charges*num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(buffer1[0],G2ww[0],num_induced_charges*num_induced_charges,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// calculate G3ww and Rww
|
||||
// G3ww is implemented as in _exact(), but can be optionally excluded due to its minor contribution
|
||||
// G3ww is implemented as in _exact(), but can be optionally excluded
|
||||
// due to its minor contribution
|
||||
// fill up buffer1 with local G3ww
|
||||
|
||||
for (int i = 0; i < num_induced_charges; i++)
|
||||
@ -759,13 +770,17 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
int k = jlist[kk] & NEIGHMASK;
|
||||
|
||||
if (mask[k] & groupbit) {
|
||||
|
||||
// interface particles: k can be ghost atoms
|
||||
|
||||
int mk = tag2mat[tag[k]];
|
||||
|
||||
double a1 = em[i] * (em[k] - 1.0);
|
||||
double a2 = 1.0 - em[i] - em[k];
|
||||
|
||||
// The first term (w/ G1ww) contributes the most to Rww
|
||||
// the second term (w/ G2ww) includes certain correction
|
||||
|
||||
//Rww[mi][mk] = a1 * G1ww[mi][mk] + a2 * G2ww[mi][mk];
|
||||
buffer1[mi][mk] = a1 * G1ww[mi][mk] + a2 * G2ww[mi][mk];
|
||||
|
||||
@ -804,9 +819,11 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
|
||||
// including the diagonal term
|
||||
double a1 = em[i] * (em[i] - 1.0);
|
||||
double a2 = 1.0 - em[i] - em[i];
|
||||
|
||||
// The first term (w/ G1ww) contributes the most to Rww
|
||||
// the second term (w/ G2ww) includes certain correction
|
||||
//Rww[mi][mi] = a1 * G1ww[mi][mi] + a2 * G2ww[mi][mi];
|
||||
|
||||
buffer1[mi][mi] = a1 * G1ww[mi][mi] + a2 * G2ww[mi][mi];
|
||||
}
|
||||
}
|
||||
@ -846,10 +863,10 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
|
||||
// invoke full neighbor list
|
||||
|
||||
int inum,*ilist,*jlist,*numneigh,**firstneigh;
|
||||
inum = list->inum; // number of entries in the neighbor list
|
||||
ilist = list->ilist; // ilist[ii] gives the atom index of the entry ii in the list
|
||||
numneigh = list->numneigh; // numneigh[i] gives the number of neighbors of local atom i
|
||||
firstneigh = list->firstneigh; // firstneigh[i] gives the pointer to the neighbors of i
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// calculate G1qw_real
|
||||
// fill up buffer1 with local G1qw_real
|
||||
@ -886,7 +903,8 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(&buffer1[0][0],&G1qw_real[0][0],num_ions*num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(&buffer1[0][0],&G1qw_real[0][0],num_ions*num_induced_charges,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// the following loop need the above results: gradG1wq_real
|
||||
// calculate sum1G1qw_epsilon and sum2ndotGwq_epsilon
|
||||
@ -917,7 +935,8 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
|
||||
delz = x[i][2] - ztmp;
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
|
||||
int mi = tag2mat_ions[tag[i]]; //ion_idx[i]; //get_matrix_index_from_local_index(i);
|
||||
int mi = tag2mat_ions[tag[i]]; //ion_idx[i];
|
||||
|
||||
//calculate_grad_greens_real(gradG1wq_real[mk][mi], delx, dely, delz);
|
||||
double gradG1wq[3];
|
||||
calculate_grad_greens_real(gradG1wq, delx, dely, delz);
|
||||
@ -929,8 +948,10 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
|
||||
temp_sum1 += G1qw_real[mi][mk] * q[i] / epsilon[i];
|
||||
}
|
||||
}
|
||||
|
||||
//sum1G1qw_epsilon[mk] = temp_sum1;// + ewaldDielectric->sum1G1qw_k_epsilon[mk];
|
||||
rhs1[mk] = temp_sum1;
|
||||
|
||||
//sum2ndotGwq_epsilon[mk] = MathExtra::dot3(norm[k], tempndotG);
|
||||
rhs2[mk] = MathExtra::dot3(norm[k], tempndotG);
|
||||
}
|
||||
@ -1007,14 +1028,16 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
|
||||
qiRwwVectorTemp1 += q[i] * (1.0 - em[k] / epsilon[i]) * G1qw_real[mi][mk];
|
||||
}
|
||||
}
|
||||
// qiRwwVectorTemp1 += ewaldDielectric->sum1G1qw_k[mk] - em[k] * ewaldDielectric->sum1G1qw_k_epsilon[mk];
|
||||
|
||||
// add the diagonal term
|
||||
|
||||
sum1G3qw += sum2ndotGwq_epsilon[mk] * G2ww[mk][mk] * area[k] * ed[k];
|
||||
|
||||
// qiRwwVectorTemp2 is a significant contribution, of which sum2G2wq is significant
|
||||
double qiRwwVectorTemp2 = (1.0 - 2.0 * em[k]) * sum2G2wq[mk] + sum1G2qw[mk] + 2.0 * sum1G3qw;
|
||||
// qiRqwVector[mk] = qiRwwVectorTemp1 + qiRwwVectorTemp2;
|
||||
double qiRwwVectorTemp2 = (1.0 - 2.0 * em[k]) * sum2G2wq[mk] +
|
||||
sum1G2qw[mk] + 2.0 * sum1G3qw;
|
||||
|
||||
// qiRqwVector[mk] = qiRwwVectorTemp1 + qiRwwVectorTemp2;
|
||||
rhs1[mk] = qiRwwVectorTemp1 + qiRwwVectorTemp2;
|
||||
}
|
||||
}
|
||||
@ -1071,14 +1094,16 @@ double FixPolarizeFunctional::greens_real(double r)
|
||||
double FixPolarizeFunctional::grad_greens_real_factor(double r)
|
||||
{
|
||||
double alpharij = g_ewald * r;
|
||||
double factor = erfc(alpharij) + 2.0 * alpharij / MY_PIS * exp(-(alpharij * alpharij));
|
||||
double factor = erfc(alpharij) + 2.0 * alpharij / MY_PIS *
|
||||
exp(-(alpharij * alpharij));
|
||||
double r3 = r*r*r;
|
||||
return (factor * (-1.0 / r3));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPolarizeFunctional::calculate_grad_greens_real(double *vec, double dx, double dy, double dz)
|
||||
void FixPolarizeFunctional::calculate_grad_greens_real(double *vec,
|
||||
double dx, double dy, double dz)
|
||||
{
|
||||
double r = sqrt(dx * dx + dy * dy + dz * dz);
|
||||
double real = grad_greens_real_factor(r);
|
||||
@ -1089,7 +1114,8 @@ void FixPolarizeFunctional::calculate_grad_greens_real(double *vec, double dx, d
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy, double dz)
|
||||
double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy,
|
||||
double dz)
|
||||
{
|
||||
// excluding the reciprocal term
|
||||
double r = sqrt(dx * dx + dy * dy + dz * dz);
|
||||
@ -1098,8 +1124,9 @@ double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy, doubl
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec, double dx, double dy, double dz)
|
||||
{
|
||||
void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec,
|
||||
double dx, double dy, double dz)
|
||||
{
|
||||
// real part of grad greens, excluding the reciprocal term
|
||||
calculate_grad_greens_real(vec, dx, dy, dz);
|
||||
}
|
||||
@ -1131,11 +1158,12 @@ double FixPolarizeFunctional::calculate_greens_ewald_self_vertex(double area)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixPolarizeFunctional::calculate_ndotgreens_ewald_self_vertex(double area, double curvature)
|
||||
double FixPolarizeFunctional::calculate_ndotgreens_ewald_self_vertex(double area,
|
||||
double curvature)
|
||||
{
|
||||
// this term is important, cannot be set to zero. see <Bugs in lammps implementation>.pptx
|
||||
// this term is important, cannot be set to zero
|
||||
// curvature = 1 / R, minus if norm is inverse of R to center.
|
||||
// Honghao Li's result, the same with Erik's paper, J Chem Phys 140 064903 (2014)
|
||||
|
||||
return curvature * MY_PIS / sqrt(area);
|
||||
}
|
||||
|
||||
@ -1180,8 +1208,9 @@ void FixPolarizeFunctional::cg_solver(double** A, double* b, double* x, int N)
|
||||
x[i] = x[i] + alpha*cg_p[i];
|
||||
cg_r[i] = cg_r[i] - alpha*cg_Ap[i];
|
||||
}
|
||||
|
||||
double rsq_new = inner_product(cg_r, cg_r, N);
|
||||
if (rsq_new < EPSILON) break;
|
||||
if (rsq_new < tolerance) break;
|
||||
|
||||
// beta = rsq_new / rsq
|
||||
double beta = rsq_new / rsq;
|
||||
@ -1189,5 +1218,4 @@ void FixPolarizeFunctional::cg_solver(double** A, double* b, double* x, int N)
|
||||
cg_p[i] = cg_r[i] + beta*cg_p[i];
|
||||
rsq = rsq_new;
|
||||
}
|
||||
|
||||
}
|
||||
@ -85,6 +85,7 @@ class FixPolarizeFunctional : public Fix {
|
||||
void calculate_qiRqw_cutoff();
|
||||
|
||||
// qw, qq ion-interface terms
|
||||
|
||||
double *qiRqwVector;
|
||||
double **G1qw_real;
|
||||
double *sum2G2wq;
|
||||
@ -93,10 +94,12 @@ class FixPolarizeFunctional : public Fix {
|
||||
double *sum2ndotGwq_epsilon;
|
||||
|
||||
// conjugate gradient solver
|
||||
|
||||
double *cg_r;
|
||||
double *cg_p;
|
||||
double *cg_Ap;
|
||||
double **cg_A;
|
||||
double tolerance;
|
||||
|
||||
void calculate_matrix_multiply_vector(double **, double *, double *, int);
|
||||
double inner_product(double*, double*, int);
|
||||
|
||||
16
src/atom.cpp
16
src/atom.cpp
@ -2661,6 +2661,14 @@ void *Atom::extract(const char *name)
|
||||
if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
|
||||
if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp;
|
||||
|
||||
// USER-DIELECTRIC
|
||||
if (strcmp(name,"area") == 0) return (void *) area;
|
||||
if (strcmp(name,"ed") == 0) return (void *) ed;
|
||||
if (strcmp(name,"em") == 0) return (void *) em;
|
||||
if (strcmp(name,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(name,"curvature") == 0) return (void *) curvature;
|
||||
if (strcmp(name,"q_unscaled") == 0) return (void *) q_unscaled;
|
||||
|
||||
// end of customization section
|
||||
// --------------------------------------------------------------------
|
||||
|
||||
@ -2744,6 +2752,14 @@ int Atom::extract_datatype(const char *name)
|
||||
if (strcmp(name,"dpdTheta") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"edpd_temp") == 0) return LAMMPS_DOUBLE;
|
||||
|
||||
// USER-DIELECTRIC
|
||||
if (strcmp(name,"area") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"ed") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"em") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"epsilon") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"curvature") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"q_unscaled") == 0) return LAMMPS_DOUBLE;
|
||||
|
||||
// end of customization section
|
||||
// --------------------------------------------------------------------
|
||||
|
||||
|
||||
Reference in New Issue
Block a user