--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2014, 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.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+# - Define macro to check if linking to zlib actually works,
+# because the find_package macro is content if one exists. This can
+# fail in cross-compilation environments, and we want to know about
+# zlib so the zlib TNG support is built only when it will work.
+#
+# GMX_TEST_ZLIB(VARIABLE)
+#
+# VARIABLE will be set to true if zlib support is present
+
+include(CheckLibraryExists)
+include(gmxOptionUtilities)
+function(GMX_TEST_ZLIB VARIABLE)
+ if(ZLIB_FOUND)
+ gmx_check_if_changed(_do_zlib_recompile ZLIB_INCLUDE_DIR ZLIB_LIBRARIES)
+ if(_do_zlib_recompile)
+ unset(ZLIB_LINKS_OK CACHE)
+ endif()
+ check_library_exists("${ZLIB_LIBRARIES}" "zlibVersion" "" ZLIB_LINKS_OK)
+ set(${VARIABLE} ${ZLIB_LINKS_OK} PARENT_SCOPE)
+ else()
+ set(${VARIABLE} OFF PARENT_SCOPE)
+ endif()
+endfunction()
+
+
+
#define TPX_TAG_RELEASE "release"
-/* This is the tag string which is stored in the tpx file.
- * Change this if you want to change the tpx format in a feature branch.
- * This ensures that there will not be different tpx formats around which
- * can not be distinguished.
+/*! \brief Tag string for the file format written to run input files
+ * written by this version of the code.
+ *
+ * Change this if you want to change the run input format in a feature
+ * branch. This ensures that there will not be different run input
+ * formats around which cannot be distinguished, while not causing
+ * problems rebasing the feature branch onto upstream changes. When
+ * merging with mainstream GROMACS, set this tag string back to
+ * TPX_TAG_RELEASE, and instead add an element to tpxv and set
+ * tpx_version to that.
*/
static const char *tpx_tag = TPX_TAG_RELEASE;
-
-/* The tpx_version number should be increased whenever the file format changes!
+/*! \brief Enum of values that describe the contents of a tpr file
+ * whose format matches a version number
*
- * The following comment section helps to keep track of which feature has been
- * added in which version.
+ * The enum helps the code be more self-documenting and ensure merges
+ * do not silently resolve when two patches make the same bump. When
+ * adding new functionality, add a new element to the end of this
+ * enumeration, change the definition of tpx_version, and write code
+ * below that does the right thing according to the value of
+ * file_version. */
+enum tpxv {
+ tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
+ tpxv_Use64BitRandomSeed /**< change ld_seed from int to gmx_int64_t */
+};
+
+/*! \brief Version number of the file format written to run input
+ * files by this version of the code.
*
- * version feature added
- * 96 support for ion/water position swaps (computational electrophysiology)
- * 97 switch ld_seed from int to gmx_int64_t
+ * The tpx_version number should be increased whenever the file format
+ * in the main development branch changes, generally to the highest
+ * value present in tpxv. Backward compatibility for reading old run
+ * input files is maintained by checking this version number against
+ * that of the file and then using the correct code path.
*
- * The following macros help the code be more self-documenting and
- * ensure merges do not silently resolve when two patches make the
- * same bump to the number. Unfortunately, compilers don't like
- * initializing a const int with a const int, so we have to be a bit
- * dirty and use a macro.
- */
-#define tpx_version_use_64_bit_random_seed 97
-static const int tpx_version = tpx_version_use_64_bit_random_seed;
+ * When developing a feature branch that needs to change the run input
+ * file format, change tpx_tag instead. */
+static const int tpx_version = tpxv_Use64BitRandomSeed;
/* This number should only be increased when you edit the TOPOLOGY section
gmx_fio_do_real(fio, bd_temp);
}
gmx_fio_do_real(fio, ir->bd_fric);
- if (file_version >= tpx_version_use_64_bit_random_seed)
+ if (file_version >= tpxv_Use64BitRandomSeed)
{
gmx_fio_do_int64(fio, ir->ld_seed);
}
}
/* Swap ions */
- if (file_version >= 96)
+ if (file_version >= tpxv_ComputationalElectrophysiology)
{
gmx_fio_do_int(fio, ir->eSwapCoords);
if (ir->eSwapCoords != eswapNO)
static real ener_drift(const verletbuf_atomtype_t *att, int natt,
const gmx_ffparams_t *ffp,
real kT_fac,
- real md_ljd, real dd_ljd,
- real md_ljr, real dd_ljr,
- real md_el, real dd_el,
+ real md1_ljd, real d2_ljd, real md3_ljd,
+ real md1_ljr, real d2_ljr, real md3_ljr,
+ real md1_el, real d2_el,
real r_buffer,
real rlist, real boxvol)
{
- double drift_tot, pot1, pot2, pot;
+ double drift_tot, pot1, pot2, pot3, pot;
int i, j;
real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
int ti, tj;
- real md, dd;
- real sc_fac, rsh;
+ real md1, d2, md3;
+ real sc_fac, rsh, rsh2;
double c_exp, c_erfc;
drift_tot = 0;
* pairs will partially cancel.
*/
/* -dV/dr at the cut-off for LJ + Coulomb */
- md =
- md_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
- md_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
- md_el*att[i].prop.q*att[j].prop.q;
-
- /* d2V/dr2 at the cut-off for Coulomb, we neglect LJ */
- dd =
- dd_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
- dd_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
- dd_el*att[i].prop.q*att[j].prop.q;
+ md1 =
+ md1_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ md1_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
+ md1_el*att[i].prop.q*att[j].prop.q;
+
+ /* d2V/dr2 at the cut-off for LJ + Coulomb */
+ d2 =
+ d2_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ d2_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
+ d2_el*att[i].prop.q*att[j].prop.q;
+
+ /* -d3V/dr3 at the cut-off for LJ, we neglect Coulomb */
+ md3 =
+ md3_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ md3_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12;
rsh = r_buffer;
sc_fac = 1.0;
c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
s = sqrt(s2);
+ rsh2 = rsh*rsh;
pot1 = sc_fac*
- md/2*((rsh*rsh + s2)*c_erfc - rsh*s*c_exp);
+ md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
pot2 = sc_fac*
- dd/6*(s*(rsh*rsh + 2*s2)*c_exp - rsh*(rsh*rsh + 3*s2)*c_erfc);
- pot = pot1 + pot2;
+ d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
+ pot3 =
+ md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
+ pot = pot1 + pot2 + pot3;
if (gmx_debug_at)
{
- fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d md %8.1e dd %8.1e pot1 %8.1e pot2 %8.1e pot %8.1e\n",
+ fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d -d1 %8.1e d2 %8.1e -d3 %8.1e pot1 %8.1e pot2 %8.1e pot3 %8.1e pot %8.1e\n",
att[i].n, att[j].n,
sqrt(s2i_2d), sqrt(s2i_3d),
sqrt(s2j_2d), sqrt(s2j_3d),
att[i].prop.bConstr+att[j].prop.bConstr,
- md, dd, pot1, pot2, pot);
+ md1, d2, md3,
+ pot1, pot2, pot3, pot);
}
/* Multiply by the number of atom pairs */
return area_rel/cluster_size;
}
+/* Returns the negative of the third derivative of a potential r^-p
+ * with a force-switch function, evaluated at the cut-off rc.
+ */
+static real md3_force_switch(real p, real rswitch, real rc)
+{
+ /* The switched force function is:
+ * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
+ */
+ real a, b;
+ real md3_pot, md3_sw;
+
+ a = -((p + 4)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 2));
+ b = ((p + 3)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 3));
+
+ md3_pot = (p + 2)*(p + 1)*p*pow(rc, p+3);
+ md3_sw = 2*a + 6*b*(rc - rswitch);
+
+ return md3_pot + md3_sw;
+}
+
void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
const t_inputrec *ir,
real reference_temperature,
verletbuf_atomtype_t *att = NULL;
int natt = -1, i;
double reppow;
- real md_ljd, dd_ljd, md_ljr, dd_ljr, md_el, dd_el;
+ real md1_ljd, d2_ljd, md3_ljd;
+ real md1_ljr, d2_ljr, md3_ljr;
+ real md1_el, d2_el;
real elfac;
real kT_fac, mass_min;
int ib0, ib1, ib;
fprintf(debug, "energy drift atom types: %d\n", natt);
}
- reppow = mtop->ffparams.reppow;
- md_ljd = 0;
- dd_ljd = 0;
- md_ljr = 0;
- dd_ljr = 0;
+ reppow = mtop->ffparams.reppow;
+ md1_ljd = 0;
+ d2_ljd = 0;
+ md3_ljd = 0;
+ md1_ljr = 0;
+ d2_ljr = 0;
+ md3_ljr = 0;
if (ir->vdwtype == evdwCUT)
{
- real sw_range, sw_range2;
+ real sw_range, md3_pswf;
switch (ir->vdw_modifier)
{
case eintmodNONE:
case eintmodPOTSHIFT:
/* -dV/dr of -r^-6 and r^-reppow */
- md_ljd = -6*pow(ir->rvdw, -7.0);
- md_ljr = reppow*pow(ir->rvdw, -(reppow+1));
- /* The contribution of the second derivative is negligible */
+ md1_ljd = -6*pow(ir->rvdw, -7.0);
+ md1_ljr = reppow*pow(ir->rvdw, -(reppow+1));
+ /* The contribution of the higher derivatives is negligible */
break;
case eintmodFORCESWITCH:
- /* At the cut-off: V=V'=V''=0.
- * We choose to approximate the potential over the switch
- * region using a linear force, thus quadratic potential.
- * This is a tight overestimate for too short switching
- * regions and not more than a factor 2 higher otherwise.
- */
- sw_range = ir->rvdw - ir->rvdw_switch;
- dd_ljd = -6*pow(ir->rvdw_switch, -7.0 )/sw_range;
- dd_ljr = reppow*pow(ir->rvdw_switch, -(reppow+1))/sw_range;
+ /* At the cut-off: V=V'=V''=0, so we use only V''' */
+ md3_ljd = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw);
+ md3_ljr = md3_force_switch(reppow, ir->rvdw_switch, ir->rvdw);
break;
case eintmodPOTSWITCH:
/* At the cut-off: V=V'=V''=0.
- * We choose to approximate the potential over the switch
- * region using a quadratic potential.
- * This is an overestimate close to the cut-off and can be
- * a slight underestimate close to rswitch.
+ * V''' is given by the original potential times
+ * the third derivative of the switch function.
*/
sw_range = ir->rvdw - ir->rvdw_switch;
- sw_range2 = sw_range*sw_range;
- dd_ljd = -12*pow(ir->rvdw_switch, -6.0 )/sw_range2;
- dd_ljr = 2*reppow*pow(ir->rvdw_switch, -reppow)/sw_range2;
+ md3_pswf = 60.0*pow(sw_range, -3.0);
+
+ md3_ljd = -pow(ir->rvdw, -6.0 )*md3_pswf;
+ md3_ljr = pow(ir->rvdw, -reppow)*md3_pswf;
break;
default:
gmx_incons("Unimplemented VdW modifier");
br4 = br2*br2;
br6 = br4*br2;
/* -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2), see LJ-PME equations in manual] and r^-reppow */
- md_ljd = -exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*pow(r, -7.0);
- md_ljr = reppow*pow(r, -(reppow+1));
- /* The contribution of the second derivative is negligible */
+ md1_ljd = -exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*pow(r, -7.0);
+ md1_ljr = reppow*pow(r, -(reppow+1));
+ /* The contribution of the higher derivatives is negligible */
}
else
{
elfac = ONE_4PI_EPS0/ir->epsilon_r;
/* Determine md=-dV/dr and dd=d^2V/dr^2 */
- md_el = 0;
- dd_el = 0;
+ md1_el = 0;
+ d2_el = 0;
if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
{
real eps_rf, k_rf;
if (eps_rf > 0)
{
- md_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
+ md1_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
}
- dd_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
+ d2_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
}
else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
{
real b, rc, br;
- b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
- rc = ir->rcoulomb;
- br = b*rc;
- md_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
- dd_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
+ b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
+ rc = ir->rcoulomb;
+ br = b*rc;
+ md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
+ d2_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
}
else
{
if (debug)
{
- fprintf(debug, "md_ljd %9.2e dd_ljd %9.2e\n", md_ljd, dd_ljd);
- fprintf(debug, "md_ljr %9.2e dd_ljr %9.2e\n", md_ljr, dd_ljr);
- fprintf(debug, "md_el %9.2e dd_el %9.2e\n", md_el, dd_el);
+ fprintf(debug, "md1_ljd %9.2e d2_ljd %9.2e md3_ljd %9.2e\n", md1_ljd, d2_ljd, md3_ljd);
+ fprintf(debug, "md1_ljr %9.2e d2_ljr %9.2e md3_ljr %9.2e\n", md1_ljr, d2_ljr, md3_ljr);
+ fprintf(debug, "md1_el %9.2e d2_el %9.2e\n", md1_el, d2_el);
fprintf(debug, "sqrt(kT_fac) %f\n", sqrt(kT_fac));
fprintf(debug, "mass_min %f\n", mass_min);
}
*/
drift = ener_drift(att, natt, &mtop->ffparams,
kT_fac,
- md_ljd, dd_ljd,
- md_ljr, dd_ljr,
- md_el, dd_el,
+ md1_ljd, d2_ljd, md3_ljd,
+ md1_ljr, d2_ljr, md3_ljr,
+ md1_el, d2_el,
rb,
rl, boxvol);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff