This only worked correctly with the group scheme and charge groups.
Although this could be implemented by using net charges of molecules,
having the interaction function depend on the system volume and
temperature complicates the code and correctness checking.
Change-Id: If29b28772cd83524c5984f64cfcd88cc4060ca36
atoms pairs, including self pairs, in which case the normal Coulomb term
in :eq:`eqns. %s <eqnvcrf>` and :eq:`%s <eqnfcrf>` is absent.
-Tironi *et al.* have introduced a generalized reaction field in which
-the dielectric continuum beyond the cut-off :math:`r_c` also has an
-ionic strength :math:`I` :ref:`71 <refTironi95>`. In this case we can
-rewrite the constants :math:`k_{rf}` and :math:`c_{rf}` using the
-inverse Debye screening length :math:`\kappa`:
-
-.. math:: \begin{aligned}
- \kappa^2 &=&
- \frac{2 I \,F^2}{\varepsilon_0 {\varepsilon_{rf}}RT}
- = \frac{F^2}{\varepsilon_0 {\varepsilon_{rf}}RT}\sum_{i=1}^{K} c_i z_i^2 \\
- k_{rf} &=& \frac{1}{r_c^3}\,
- \frac{({\varepsilon_{rf}}-{\varepsilon_r})(1 + \kappa r_c) + {\frac{1}{2}}{\varepsilon_{rf}}(\kappa r_c)^2}
- {(2{\varepsilon_{rf}}+ {\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2}
- \end{aligned}
- :label: eqnkgrf
-
-.. math:: \begin{aligned}
- c_{rf} &=& \frac{1}{r_c}\,
- \frac{3{\varepsilon_{rf}}(1 + \kappa r_c + {\frac{1}{2}}(\kappa r_c)^2)}
- {(2{\varepsilon_{rf}}+{\varepsilon_r})(1 + \kappa r_c) + {\varepsilon_{rf}}(\kappa r_c)^2}
- \end{aligned}
- :label: eqncgrf
-
-where :math:`F` is Faraday’s constant, :math:`R` is the ideal gas
-constant, :math:`T` the absolute temperature, :math:`c_i` the molar
-concentration for species :math:`i` and :math:`z_i` the charge number of
-species :math:`i` where we have :math:`K` different species. In the
-limit of zero ionic strength (:math:`\kappa=0`) :eq:`eqns. %s <eqnkgrf>` and
-:eq:`%s <eqncgrf>` reduce to the simple forms of :eq:`eqns. %s <eqnkrf>` and :eq:`%s <eqncrf>`
-respectively.
-
.. _modnbint:
Modified non-bonded interactions
for parallelization of particle simulations," *J. Chem. Phys.*, **124**
[18] 184109–184109 (2006).
-.. raw:: html
-
- </div>
-
-.. raw:: html
-
- <div id="ref-Tironi95">
-
-.. _refTironi95:
-
-:sup:`71` I.G. Tironi, R. Sperb, P.E. Smith, and W.F. van Gunsteren, "A
-generalized reaction field method for molecular dynamics simulations,"
-*J. Chem. Phys.*, **102** 5451–5459 (1995).
-
.. raw:: html
</div>
:issue:`1852`
+Generalized reaction-field
+""""""""""""""""""""""""""
+
+This only worked correctly with the group scheme. Note that generalized
+reaction-field simulations can still be performed using standard
+reaction field and computing the dielectric constant manually.
+
gmx anadock
"""""""""""
The gmx anadock tool was removed since it does not belong in gromacs
:mdp:`epsilon-rf`. The dielectric constant can be set to
infinity by setting :mdp:`epsilon-rf` =0.
- .. mdp-value:: Generalized-Reaction-Field
-
- Generalized reaction field with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rcoulomb`. The
- dielectric constant beyond the cut-off is
- :mdp:`epsilon-rf`. The ionic strength is computed from the
- number of charged (*i.e.* with non zero charge) charge
- groups. The temperature for the GRF potential is set with
- :mdp:`ref-t`.
-
.. mdp-value:: Reaction-Field-zero
In |Gromacs|, normal reaction-field electrostatics with
}
/* Generalized reaction field */
- if (ir->opts.ngtc == 0)
+ if (ir->coulombtype == eelGRF_NOTUSED)
{
- sprintf(err_buf, "No temperature coupling while using coulombtype %s",
- eel_names[eelGRF]);
- CHECK(ir->coulombtype == eelGRF);
- }
- else
- {
- sprintf(err_buf, "When using coulombtype = %s"
- " ref-t for temperature coupling should be > 0",
- eel_names[eelGRF]);
- CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
+ warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
+ "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
}
bAcc = FALSE;
std::abs(fr->qsum[1]) > 1e-4);
}
-void update_forcerec(t_forcerec *fr, matrix box)
-{
- if (fr->ic->eeltype == eelGRF)
- {
- calc_rffac(nullptr, fr->ic->eeltype, fr->ic->epsilon_r, fr->ic->epsilon_rf,
- fr->ic->rcoulomb, fr->temp, fr->zsquare, box,
- &fr->ic->k_rf, &fr->ic->c_rf);
- }
-}
-
static real calcBuckinghamBMax(FILE *fplog, const gmx_mtop_t *mtop)
{
const t_atoms *at1, *at2;
/* Reaction-field */
if (EEL_RF(ic->eeltype))
{
+ GMX_RELEASE_ASSERT(ic->eeltype != eelGRF_NOTUSED, "GRF is no longer supported");
ic->epsilon_rf = ir->epsilon_rf;
- /* Generalized reaction field parameters are updated every step */
- if (ic->eeltype != eelGRF)
- {
- calc_rffac(fp, ic->eeltype, ic->epsilon_r, ic->epsilon_rf,
- ic->rcoulomb, 0, 0, nullptr,
- &ic->k_rf, &ic->c_rf);
- }
+
+ calc_rffac(fp, ic->epsilon_r, ic->epsilon_rf,
+ ic->rcoulomb,
+ &ic->k_rf, &ic->c_rf);
}
else
{
fr->n_tpi = 0;
}
- if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
+ if (ir->coulombtype == eelRF_NEC_UNSUPPORTED ||
+ ir->coulombtype == eelGRF_NOTUSED)
{
gmx_fatal(FARGS, "%s electrostatics is no longer supported",
eel_names[ir->coulombtype]);
break;
case eelRF:
- case eelGRF:
fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_REACTIONFIELD;
break;
/* 1-4 interaction electrostatics */
fr->fudgeQQ = mtop->ffparams.fudgeQQ;
- /* Parameters for generalized RF */
- fr->zsquare = 0.0;
- fr->temp = 0.0;
-
- if (ic->eeltype == eelGRF)
- {
- init_generalized_rf(fp, mtop, ir, fr);
- }
-
fr->haveDirectVirialContributions =
(EEL_FULL(ic->eeltype) || EVDW_PME(ic->vdwtype) ||
fr->forceProviders->hasForceProvider() ||
void forcerec_set_excl_load(t_forcerec *fr,
const gmx_localtop_t *top);
-/*! \brief Update parameters dependent on box
- *
- * Updates parameters in the forcerec that are time dependent
- * \param[out] fr The force record
- * \param[in] box The simulation box
- */
-void update_forcerec(t_forcerec *fr, matrix box);
-
void free_gpu_resources(t_forcerec *fr,
const gmx::PhysicalNodeCommunicator &physicalNodeCommunicator,
const gmx_gpu_info_t &gpu_info);
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/pleasecite.h"
-void calc_rffac(FILE *fplog, int eel, real eps_r, real eps_rf, real Rc, real Temp,
- real zsq, matrix box,
+void calc_rffac(FILE *fplog, real eps_r, real eps_rf, real Rc,
real *krf, real *crf)
{
- /* Compute constants for Generalized reaction field */
- real kappa, k1, k2, I, rmin;
- real vol = 0;
-
- if (EEL_RF(eel))
- {
- if (eel == eelGRF)
- {
- /* Consistency check */
- if (Temp <= 0.0)
- {
- gmx_fatal(FARGS, "Temperature is %f while using"
- " Generalized Reaction Field\n", Temp);
- }
- /* Ionic strength (only needed for eelGRF */
- vol = det(box);
- I = 0.5*zsq/vol;
- kappa = std::sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
- }
- else
- {
- I = 0;
- kappa = 0;
- }
-
- /* eps == 0 signals infinite dielectric */
- if (eps_rf == 0)
- {
- *krf = 1/(2*Rc*Rc*Rc);
- }
- else
- {
- k1 = 1 + kappa*Rc;
- k2 = eps_rf*gmx::square(static_cast<real>(kappa*Rc));
-
- *krf = ((eps_rf - eps_r)*k1 + 0.5*k2)/((2*eps_rf + eps_r)*k1 + k2)/(Rc*Rc*Rc);
- }
- *crf = 1/Rc + *krf*Rc*Rc;
- // Make sure we don't lose resolution in pow() by casting real arg to double
- rmin = gmx::invcbrt(static_cast<double>(*krf*2.0));
-
- if (fplog)
- {
- if (eel == eelGRF)
- {
- please_cite(fplog, "Tironi95a");
- fprintf(fplog, "%s:\n"
- "epsRF = %10g, I = %10g, volume = %10g, kappa = %10g\n"
- "rc = %10g, krf = %10g, crf = %10g, epsfac = %10g\n",
- eel_names[eel], eps_rf, I, vol, kappa, Rc, *krf, *crf,
- ONE_4PI_EPS0/eps_r);
- }
- else
- {
- fprintf(fplog, "%s:\n"
- "epsRF = %g, rc = %g, krf = %g, crf = %g, epsfac = %g\n",
- eel_names[eel], eps_rf, Rc, *krf, *crf, ONE_4PI_EPS0/eps_r);
- }
- fprintf(fplog,
- "The electrostatics potential has its minimum at r = %g\n",
- rmin);
- }
- }
-}
-
-void init_generalized_rf(FILE *fplog,
- const gmx_mtop_t *mtop, const t_inputrec *ir,
- t_forcerec *fr)
-{
- int i, j;
- real q, zsq, nrdf, T;
- const gmx_moltype_t *molt;
- const t_block *cgs;
-
- if (ir->efep != efepNO && fplog)
+ /* eps == 0 signals infinite dielectric */
+ if (eps_rf == 0)
{
- fprintf(fplog, "\nWARNING: the generalized reaction field constants are determined from topology A only\n\n");
+ *krf = 1/(2*Rc*Rc*Rc);
}
- zsq = 0.0;
- for (const gmx_molblock_t &molb : mtop->molblock)
+ else
{
- molt = &mtop->moltype[molb.type];
- cgs = &molt->cgs;
- for (i = 0; (i < cgs->nr); i++)
- {
- q = 0;
- for (j = cgs->index[i]; (j < cgs->index[i+1]); j++)
- {
- q += molt->atoms.atom[j].q;
- }
- zsq += molb.nmol*q*q;
- }
- fr->zsquare = zsq;
+ *krf = (eps_rf - eps_r)/(2*eps_rf + eps_r)/(Rc*Rc*Rc);
}
+ *crf = 1/Rc + *krf*Rc*Rc;
- T = 0.0;
- nrdf = 0.0;
- for (i = 0; (i < ir->opts.ngtc); i++)
+ if (fplog)
{
- nrdf += ir->opts.nrdf[i];
- T += (ir->opts.nrdf[i] * ir->opts.ref_t[i]);
- }
- if (nrdf == 0)
- {
- gmx_fatal(FARGS, "No degrees of freedom!");
+ fprintf(fplog, "%s:\n"
+ "epsRF = %g, rc = %g, krf = %g, crf = %g, epsfac = %g\n",
+ eel_names[eelRF], eps_rf, Rc, *krf, *crf, ONE_4PI_EPS0/eps_r);
+ // Make sure we don't lose resolution in pow() by casting real arg to double
+ real rmin = gmx::invcbrt(static_cast<double>(*krf*2.0));
+ fprintf(fplog,
+ "The electrostatics potential has its minimum at r = %g\n",
+ rmin);
}
- fr->temp = T/nrdf;
}
struct t_forcerec;
struct t_inputrec;
-void calc_rffac(FILE *fplog, int eel, real eps_r, real eps_rf,
- real Rc, real Temp,
- real zsq, matrix box,
+void calc_rffac(FILE *fplog, real eps_r, real eps_rf,
+ real Rc,
real *krf, real *crf);
/* Determine the reaction-field constants */
-void init_generalized_rf(FILE *fplog,
- const gmx_mtop_t *mtop, const t_inputrec *ir,
- t_forcerec *fr);
-/* Initialize the generalized reaction field parameters */
-
-
#endif
if (bStateChanged)
{
- update_forcerec(fr, box);
-
if (inputrecNeedMutot(inputrec))
{
/* Calculate total (local) dipole moment in a temporary common array.
*/
real rlist = 0;
- /* Parameters for generalized reaction field */
- real zsquare = 0;
- real temp = 0;
-
/* Charge sum and dipole for topology A/B ([0]/[1]) for Ewald corrections */
double qsum[2] = { 0 };
double q2sum[2] = { 0 };
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
};
const char *eel_names[eelNR+1] = {
- "Cut-off", "Reaction-Field", "Generalized-Reaction-Field",
+ "Cut-off", "Reaction-Field", "Generalized-Reaction-Field (unused)",
"PME", "Ewald", "P3M-AD", "Poisson", "Switch", "Shift", "User",
"Generalized-Born (unused)", "Reaction-Field-nec", "Encad-shift",
"PME-User", "PME-Switch", "PME-User-Switch",
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/*! \brief Cut-off treatment for Coulomb */
enum {
- eelCUT, eelRF, eelGRF, eelPME, eelEWALD, eelP3M_AD,
+ eelCUT, eelRF, eelGRF_NOTUSED, eelPME, eelEWALD, eelP3M_AD,
eelPOISSON, eelSWITCH, eelSHIFT, eelUSER, eelGB_NOTUSED, eelRF_NEC_UNSUPPORTED, eelENCADSHIFT,
eelPMEUSER, eelPMESWITCH, eelPMEUSERSWITCH, eelRF_ZERO, eelNR
};
extern const char *eewg_names[eewgNR+1];
//! Macro telling us whether we use reaction field
-#define EEL_RF(e) ((e) == eelRF || (e) == eelGRF || (e) == eelRF_NEC_UNSUPPORTED || (e) == eelRF_ZERO )
+#define EEL_RF(e) ((e) == eelRF || (e) == eelGRF_NOTUSED || (e) == eelRF_NEC_UNSUPPORTED || (e) == eelRF_ZERO )
//! Macro telling us whether we use PME
#define EEL_PME(e) ((e) == eelPME || (e) == eelPMESWITCH || (e) == eelPMEUSER || (e) == eelPMEUSERSWITCH || (e) == eelP3M_AD)
tabsel[etiCOUL] = etabEwaldUserSwitch;
break;
case eelRF:
- case eelGRF:
case eelRF_ZERO:
tabsel[etiCOUL] = etabRF_ZERO;
break;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff