* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <ctype.h>
+#include <stdlib.h>
#include <string.h>
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
-#include "gromacs/fileio/confio.h"
#include "gromacs/commandline/pargs.h"
-#include "pbc.h"
-#include "force.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
-#include "gromacs/math/utilities.h"
-#include "macros.h"
-#include "vec.h"
+#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/tpxio.h"
-#include "mdrun.h"
-#include "main.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/random/random.h"
-#include "index.h"
-#include "mtop_util.h"
-#include "gmx_ana.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
static void insert_ion(int nsa, int *nwater,
gmx_bool bSet[], int repl[], atom_id index[],
char *str;
int i;
- str = strdup(mname);
+ str = gmx_strdup(mname);
i = strlen(str)-1;
while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
{
if (np)
{
snew(pptr, 1);
- pptr[0] = strdup(p_name);
+ pptr[0] = gmx_strdup(p_name);
snew(paptr, 1);
paptr[0] = aname(p_name);
}
if (nn)
{
snew(nptr, 1);
- nptr[0] = strdup(n_name);
+ nptr[0] = gmx_strdup(n_name);
snew(naptr, 1);
naptr[0] = aname(n_name);
}
static void update_topol(const char *topinout, int p_num, int n_num,
const char *p_name, const char *n_name, char *grpname)
{
- #define TEMP_FILENM "temp.top"
FILE *fpin, *fpout;
char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
int line, i, nsol, nmol_line, sol_line, nsol_last;
gmx_bool bMolecules;
+ char temporary_filename[STRLEN];
printf("\nProcessing topology\n");
fpin = gmx_ffopen(topinout, "r");
- fpout = gmx_ffopen(TEMP_FILENM, "w");
+ strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
+ gmx_tmpnam(temporary_filename);
+ fpout = gmx_ffopen(temporary_filename, "w");
line = 0;
bMolecules = FALSE;
}
/* Store this molecules section line */
srenew(mol_line, nmol_line+1);
- mol_line[nmol_line] = strdup(buf);
+ mol_line[nmol_line] = gmx_strdup(buf);
nmol_line++;
}
}
/* use gmx_ffopen to generate backup of topinout */
fpout = gmx_ffopen(topinout, "w");
gmx_ffclose(fpout);
- rename(TEMP_FILENM, topinout);
- #undef TEMP_FILENM
+ rename(temporary_filename, topinout);
}
int gmx_genion(int argc, char *argv[])
{ "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
{ "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
{ "-conc", FALSE, etREAL, {&conc},
- "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [TT].tpr[tt] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
+ "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [REF].tpr[ref] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
{ "-neutral", FALSE, etBOOL, {&bNeutral}, "This option will add enough ions to neutralize the system. These ions are added on top of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. "}
};
t_topology top;
output_env_t oenv;
gmx_rng_t rng;
t_filenm fnm[] = {
- { efTPX, NULL, NULL, ffREAD },
+ { efTPR, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
{ efSTO, "-o", NULL, ffWRITE },
{ efTOP, "-p", "topol", ffOPTRW }
};
#define NFILE asize(fnm)
- if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
+ if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
asize(desc), desc, asize(bugs), bugs, &oenv))
{
return 0;
}
/* Read atom positions and charges */
- read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
+ read_tps_conf(ftp2fn(efTPR, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
atoms = top.atoms;
/* Compute total charge */
* the research papers on the package. Check out http://www.gromacs.org.
*/
/* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
-#include "typedefs.h"
-#include "main.h"
-#include "network.h"
-#include "rbin.h"
+#include "gromacs/legacyheaders/rbin.h"
+
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/typedefs.h"
#include "gromacs/utility/smalloc.h"
t_bin *mk_bin(void)
/* Copy pointer */
rbuf = b->rbuf+b->nreal;
- #if (__ICC == 1500 || __ICL == 1500) && defined __MIC__
+ #if (__ICC >= 1500 || __ICL >= 1500) && defined __MIC__
#pragma novector /* Work-around for incorrect vectorization */
#endif
for (i = 0; (i < nr); i++)
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, 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.
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#include "solvate.h"
+#include "gmxpre.h"
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "solvate.h"
#include <string.h>
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/fileio/confio.h"
-#include "macros.h"
-#include "gromacs/fileio/futil.h"
-#include "atomprop.h"
-#include "names.h"
-#include "vec.h"
-#include "gmx_fatal.h"
+#include <algorithm>
+
#include "gromacs/commandline/pargs.h"
-#include "gromacs/gmxlib/conformation-utilities.h"
-#include "addconf.h"
-#include "read-conformation.h"
+#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/pdbio.h"
-#include "pbc.h"
-
-#ifdef DEBUG
-static void print_stat(rvec *x, int natoms, matrix box)
-{
- int i, m;
- rvec xmin, xmax;
- for (m = 0; (m < DIM); m++)
- {
- xmin[m] = x[0][m];
- xmax[m] = x[0][m];
- }
- for (i = 0; (i < natoms); i++)
- {
- for (m = 0; m < DIM; m++)
- {
- xmin[m] = min(xmin[m], x[i][m]);
- xmax[m] = max(xmax[m], x[i][m]);
- }
- }
- for (m = 0; (m < DIM); m++)
- {
- fprintf(stderr, "DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
- m, xmin[m], xmax[m], box[m][m]);
- }
-}
-#endif
+#include "gromacs/gmxlib/conformation-utilities.h"
+#include "gromacs/gmxpreprocess/read-conformation.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/selection/nbsearch.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/atomsbuilder.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
char *name;
if ( (i == 0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) )
{
/* see if this was a molecule type we haven't had yet: */
- moltp = NOTSET;
- for (j = 0; (j < nrmoltypes) && (moltp == NOTSET); j++)
+ moltp = -1;
+ for (j = 0; (j < nrmoltypes) && (moltp == -1); j++)
{
- if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
- moltypes[j].name) == 0)
+ /* cppcheck-suppress nullPointer
+ * moltypes is guaranteed to be allocated because otherwise
+ * nrmoltypes is 0. */
+ if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name), moltypes[j].name) == 0)
{
moltp = j;
}
}
- if (moltp == NOTSET)
+ if (moltp == -1)
{
moltp = nrmoltypes;
nrmoltypes++;
}
}
-/* Make a new configuration by adding boxes*/
-static void make_new_conformation(t_atoms *atoms, rvec *x, rvec *v, real *r, matrix box, ivec n_box)
+/*! \brief
+ * Generates a solvent configuration of desired size by stacking solvent boxes.
+ *
+ * \param[in,out] atoms Solvent atoms.
+ * \param[in,out] x Solvent positions.
+ * \param[in,out] v Solvent velocities (`*v` can be NULL).
+ * \param[in,out] r Solvent exclusion radii.
+ * \param[in] box Initial solvent box.
+ * \param[in] boxTarget Target box size.
+ *
+ * The solvent box of desired size is created by stacking the initial box in
+ * the smallest k*l*m array that covers the box, and then removing any residue
+ * where all atoms are outside the target box (with a small margin).
+ * This function does not remove overlap between solvent atoms across the
+ * edges.
+ *
+ * Note that the input configuration should be in the rectangular unit cell and
+ * have whole residues.
+ */
+static void replicateSolventBox(t_atoms *atoms, rvec **x, rvec **v, real **r,
+ const matrix box, const matrix boxTarget)
{
- int i, ix, iy, iz, m, j, imol, offset;
- rvec delta;
- int nmol;
+ // Calculate the box multiplication factors.
+ ivec n_box;
+ int nmol = 1;
+ for (int i = 0; i < DIM; ++i)
+ {
+ n_box[i] = 1;
+ while (n_box[i] * box[i][i] < boxTarget[i][i])
+ {
+ n_box[i]++;
+ }
+ nmol *= n_box[i];
+ }
+ fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
+ n_box[XX], n_box[YY], n_box[ZZ]);
+
+ // Create arrays for storing the generated system (cannot be done in-place
+ // in case the target box is smaller than the original in one dimension,
+ // but not in all).
+ t_atoms newAtoms;
+ init_t_atoms(&newAtoms, 0, FALSE);
+ gmx::AtomsBuilder builder(&newAtoms);
+ builder.reserve(atoms->nr * nmol, atoms->nres * nmol);
+ rvec *newX;
+ rvec *newV = NULL;
+ real *newR;
+ snew(newX, atoms->nr * nmol);
+ snew(newR, atoms->nr * nmol);
+ if (*v)
+ {
+ snew(newV, atoms->nr * nmol);
+ }
- nmol = n_box[XX]*n_box[YY]*n_box[ZZ];
+ const real maxRadius = *std::max_element(*r, *r + atoms->nr);
+ rvec boxWithMargin;
+ for (int i = 0; i < DIM; ++i)
+ {
+ // The code below is only interested about the box diagonal.
+ boxWithMargin[i] = boxTarget[i][i] + 3*maxRadius;
+ }
- /*print message*/
- fprintf(stderr, "Generating configuration\n");
- imol = 0;
- for (ix = 0; (ix < n_box[XX]); ix++)
+ for (int ix = 0; ix < n_box[XX]; ++ix)
{
+ rvec delta;
delta[XX] = ix*box[XX][XX];
- for (iy = 0; (iy < n_box[YY]); iy++)
+ for (int iy = 0; iy < n_box[YY]; ++iy)
{
delta[YY] = iy*box[YY][YY];
- for (iz = 0; (iz < n_box[ZZ]); iz++)
+ for (int iz = 0; iz < n_box[ZZ]; ++iz)
{
delta[ZZ] = iz*box[ZZ][ZZ];
- offset = imol*atoms->nr;
- for (i = 0; (i < atoms->nr); i++)
+ bool bKeepResidue = false;
+ for (int i = 0; i < atoms->nr; ++i)
{
- for (m = 0; (m < DIM); m++)
+ const int newIndex = builder.currentAtomCount();
+ bool bKeepAtom = true;
+ for (int m = 0; m < DIM; ++m)
{
- x[offset+i][m] = delta[m]+x[i][m];
+ const real newCoord = delta[m] + (*x)[i][m];
+ bKeepAtom = bKeepAtom && (newCoord < boxWithMargin[m]);
+ newX[newIndex][m] = newCoord;
}
- if (v)
+ bKeepResidue = bKeepResidue || bKeepAtom;
+ if (newV)
{
- for (m = 0; (m < DIM); m++)
+ copy_rvec((*v)[i], newV[newIndex]);
+ }
+ newR[newIndex] = (*r)[i];
+ builder.addAtom(*atoms, i);
+ if (i == atoms->nr - 1
+ || atoms->atom[i+1].resind != atoms->atom[i].resind)
+ {
+ if (bKeepResidue)
+ {
+ builder.finishResidue(atoms->resinfo[atoms->atom[i].resind]);
+ }
+ else
{
- v[offset+i][m] = v[i][m];
+ builder.discardCurrentResidue();
}
+ // Reset state for the next residue.
+ bKeepResidue = false;
}
- r[offset+i] = r[i];
}
- imol++;
}
}
}
- for (i = 1; (i < nmol); i++)
+ sfree(atoms->atom);
+ sfree(atoms->atomname);
+ sfree(atoms->resinfo);
+ atoms->nr = newAtoms.nr;
+ atoms->nres = newAtoms.nres;
+ atoms->atom = newAtoms.atom;
+ atoms->atomname = newAtoms.atomname;
+ atoms->resinfo = newAtoms.resinfo;
+ sfree(*x);
+ sfree(*v);
+ sfree(*r);
+ *x = newX;
+ *v = newV;
+ *r = newR;
+ fprintf(stderr, "Solvent box contains %d atoms in %d residues\n",
+ atoms->nr, atoms->nres);
+}
+
+/*! \brief
+ * Removes overlap of solvent atoms across the edges.
+ *
+ * \param[in,out] atoms Solvent atoms.
+ * \param[in,out] x Solvent positions.
+ * \param[in,out] v Solvent velocities (can be NULL).
+ * \param[in,out] r Solvent exclusion radii.
+ * \param[in] pbc PBC information.
+ *
+ * Solvent residues that lay on the edges that do not touch the origin are
+ * removed if they overlap with other solvent atoms across the PBC.
+ * This is done in this way as the assumption is that the input solvent
+ * configuration is already equilibrated, and so does not contain any
+ * undesirable overlap. The only overlap that should be removed is caused by
+ * cutting the box in half in replicateSolventBox() and leaving a margin of
+ * solvent outside those box edges; these atoms can then overlap with those on
+ * the opposite box edge in a way that is not part of the pre-equilibrated
+ * configuration.
+ */
+static void removeSolventBoxOverlap(t_atoms *atoms, rvec *x, rvec *v, real *r,
+ const t_pbc &pbc)
+{
+ gmx::AtomsRemover remover(*atoms);
+
+ // TODO: We could limit the amount of pairs searched significantly,
+ // since we are only interested in pairs where the positions are on
+ // opposite edges.
+ const real maxRadius = *std::max_element(r, r + atoms->nr);
+ gmx::AnalysisNeighborhood nb;
+ nb.setCutoff(2*maxRadius);
+ gmx::AnalysisNeighborhoodPositions pos(x, atoms->nr);
+ gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, pos);
+ gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
+ gmx::AnalysisNeighborhoodPair pair;
+ while (pairSearch.findNextPair(&pair))
{
- int offs = i*atoms->nr;
- int offsres = i*atoms->nres;
- for (j = 0; (j < atoms->nr); j++)
+ const int i1 = pair.refIndex();
+ const int i2 = pair.testIndex();
+ if (remover.isMarked(i2))
+ {
+ pairSearch.skipRemainingPairsForTestPosition();
+ continue;
+ }
+ if (remover.isMarked(i1) || atoms->atom[i1].resind == atoms->atom[i2].resind)
{
- atoms->atomname[offs+j] = atoms->atomname[j];
- atoms->atom[offs+j].resind = atoms->atom[j].resind + offsres;
- atoms->resinfo[atoms->atom[offs+j].resind] =
- atoms->resinfo[atoms->atom[j].resind];
- atoms->resinfo[atoms->atom[offs+j].resind].nr += offsres;
+ continue;
}
+ if (pair.distance2() < sqr(r[i1] + r[i2]))
+ {
+ rvec dx;
+ rvec_sub(x[i2], x[i1], dx);
+ bool bCandidate1 = false, bCandidate2 = false;
+ // To satisfy Clang static analyzer.
+ GMX_ASSERT(pbc.ndim_ePBC <= DIM, "Too many periodic dimensions");
+ for (int d = 0; d < pbc.ndim_ePBC; ++d)
+ {
+ // If the distance in some dimension is larger than the
+ // cutoff, then it means that the distance has been computed
+ // over the PBC. Mark the position with a larger coordinate
+ // for potential removal.
+ if (dx[d] > maxRadius)
+ {
+ bCandidate2 = true;
+ }
+ else if (dx[d] < -maxRadius)
+ {
+ bCandidate1 = true;
+ }
+ }
+ // Only mark one of the positions for removal if both were
+ // candidates.
+ if (bCandidate2 && (!bCandidate1 || i2 > i1))
+ {
+ remover.markResidue(*atoms, i2, true);
+ pairSearch.skipRemainingPairsForTestPosition();
+ }
+ else if (bCandidate1)
+ {
+ remover.markResidue(*atoms, i1, true);
+ }
+ }
+ }
+
+ remover.removeMarkedVectors(x);
+ if (v != NULL)
+ {
+ remover.removeMarkedVectors(v);
+ }
+ remover.removeMarkedValues(r);
+ const int originalAtomCount = atoms->nr;
+ remover.removeMarkedAtoms(atoms);
+ fprintf(stderr, "Removed %d solvent atoms due to solvent-solvent overlap\n",
+ originalAtomCount - atoms->nr);
+}
+
+/*! \brief
+ * Removes solvent molecules that overlap with the solute, and optionally also
+ * those that are outside a given shell radius from the solute.
+ *
+ * \param[in,out] atoms Solvent atoms.
+ * \param[in,out] x Solvent positions.
+ * \param[in,out] v Solvent velocities (can be NULL).
+ * \param[in,out] r Solvent exclusion radii.
+ * \param[in] pbc PBC information.
+ * \param[in] soluteAtomCount Number of solute atoms.
+ * \param[in] x_solute Solute positions.
+ * \param[in] r_solute Solute exclusion radii.
+ * \param[in] rshell If >0, only keep solvent atoms within a shell of
+ * this size from the solute.
+ */
+static void removeSoluteOverlap(t_atoms *atoms, rvec *x, rvec *v, real *r,
+ const t_pbc &pbc, int soluteAtomCount,
+ const rvec *x_solute, const real *r_solute,
+ real rshell)
+{
+ const real maxRadius1
+ = *std::max_element(r, r + atoms->nr);
+ const real maxRadius2
+ = *std::max_element(r_solute, r_solute + soluteAtomCount);
+
+ gmx::AtomsRemover remover(*atoms);
+ // If rshell is >0, the neighborhood search looks at all pairs
+ // within rshell, and unmarks those that are within the cutoff.
+ // This line marks everything, so that solvent outside rshell remains
+ // marked after the loop.
+ // Without rshell, the neighborhood search only marks the overlapping
+ // solvent atoms, and all others are left alone.
+ if (rshell > 0.0)
+ {
+ remover.markAll();
+ }
+
+ gmx::AnalysisNeighborhood nb;
+ nb.setCutoff(std::max(maxRadius1 + maxRadius2, rshell));
+ gmx::AnalysisNeighborhoodPositions posSolute(x_solute, soluteAtomCount);
+ gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, posSolute);
+ gmx::AnalysisNeighborhoodPositions pos(x, atoms->nr);
+ gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
+ gmx::AnalysisNeighborhoodPair pair;
+ while (pairSearch.findNextPair(&pair))
+ {
+ if (remover.isMarked(pair.testIndex()))
+ {
+ pairSearch.skipRemainingPairsForTestPosition();
+ continue;
+ }
+ const real r1 = r_solute[pair.refIndex()];
+ const real r2 = r[pair.testIndex()];
+ const bool bRemove = (pair.distance2() < sqr(r1 + r2));
+ remover.markResidue(*atoms, pair.testIndex(), bRemove);
}
- atoms->nr *= nmol;
- atoms->nres *= nmol;
- for (i = 0; i < DIM; i++)
+
+ remover.removeMarkedVectors(x);
+ if (v != NULL)
+ {
+ remover.removeMarkedVectors(v);
+ }
+ remover.removeMarkedValues(r);
+ const int originalAtomCount = atoms->nr;
+ remover.removeMarkedAtoms(atoms);
+ fprintf(stderr, "Removed %d solvent atoms due to solute-solvent overlap\n",
+ originalAtomCount - atoms->nr);
+}
+
+/*! \brief
+ * Removes a given number of solvent residues.
+ *
+ * \param[in,out] atoms Solvent atoms.
+ * \param[in,out] x Solvent positions.
+ * \param[in,out] v Solvent velocities (can be NULL).
+ * \param[in] numberToRemove Number of residues to remove.
+ *
+ * This function is called last in the process of creating the solvent box,
+ * so it does not operate on the exclusion radii, as no code after this needs
+ * them.
+ */
+static void removeExtraSolventMolecules(t_atoms *atoms, rvec *x, rvec *v,
+ int numberToRemove)
+{
+ gmx::AtomsRemover remover(*atoms);
+ // TODO: It might be nicer to remove a random set of residues, but
+ // in practice this should give a roughly uniform spatial distribution.
+ const int stride = atoms->nr / numberToRemove;
+ for (int i = 0; i < numberToRemove; ++i)
{
- for (j = 0; j < DIM; j++)
+ int atomIndex = (i+1)*stride - 1;
+ while (remover.isMarked(atomIndex))
{
- box[j][i] *= n_box[j];
+ ++atomIndex;
+ if (atomIndex == atoms->nr)
+ {
+ atomIndex = 0;
+ }
}
+ remover.markResidue(*atoms, atomIndex, true);
+ }
+ remover.removeMarkedVectors(x);
+ if (v != NULL)
+ {
+ remover.removeMarkedVectors(v);
}
+ remover.removeMarkedAtoms(atoms);
}
-static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **exclusionDistances,
- int ePBC, matrix box,
- gmx_atomprop_t aps,
- real defaultDistance, real scaleFactor, int *atoms_added,
- int *residues_added, real rshell, int max_sol,
- const output_env_t oenv)
+static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v,
+ int ePBC, matrix box, gmx_atomprop_t aps,
+ real defaultDistance, real scaleFactor,
+ real rshell, int max_sol)
{
- int i, nmol;
- ivec n_box;
- char filename[STRLEN];
- char title_solvt[STRLEN];
t_atoms *atoms_solvt;
rvec *x_solvt, *v_solvt = NULL;
- real *exclusionDistances_solvt;
int ePBC_solvt;
matrix box_solvt;
- int onr, onres;
- char *lfn;
- lfn = gmxlibfn(fn);
- strncpy(filename, lfn, STRLEN);
- sfree(lfn);
+ char *filename = gmxlibfn(fn);
+ snew(atoms_solvt, 1);
+ char *title_solvt
+ = readConformation(filename, atoms_solvt, &x_solvt, &v_solvt,
+ &ePBC_solvt, box_solvt, "solvent");
+ sfree(title_solvt);
+ if (0 == atoms_solvt->nr)
{
- int natoms;
- get_stx_coordnum(filename, &natoms);
- if (0 == natoms)
- {
- gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
- }
- snew(atoms_solvt, 1);
- init_t_atoms(atoms_solvt, natoms, FALSE);
- }
- snew(x_solvt, atoms_solvt->nr);
- if (v)
- {
- snew(v_solvt, atoms_solvt->nr);
- }
- snew(exclusionDistances_solvt, atoms_solvt->nr);
- snew(atoms_solvt->resinfo, atoms_solvt->nr);
- snew(atoms_solvt->atomname, atoms_solvt->nr);
- snew(atoms_solvt->atom, atoms_solvt->nr);
- atoms_solvt->pdbinfo = NULL;
- fprintf(stderr, "Reading solvent configuration%s\n",
- v_solvt ? " and velocities" : "");
- read_stx_conf(filename, title_solvt, atoms_solvt, x_solvt, v_solvt,
- &ePBC_solvt, box_solvt);
- fprintf(stderr, "\"%s\"\n", title_solvt);
- fprintf(stderr, "solvent configuration contains %d atoms in %d residues\n",
- atoms_solvt->nr, atoms_solvt->nres);
+ gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
+ }
+ sfree(filename);
fprintf(stderr, "\n");
/* apply pbc for solvent configuration for whole molecules */
rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
/* initialise distance arrays for solvent configuration */
- exclusionDistances_solvt = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);
+ fprintf(stderr, "Initialising inter-atomic distances...\n");
+ real *exclusionDistances
+ = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
+ real *exclusionDistances_solvt
+ = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);
- /* calculate the box multiplication factors n_box[0...DIM] */
- nmol = 1;
- for (i = 0; (i < DIM); i++)
+ /* generate a new solvent configuration */
+ fprintf(stderr, "Generating solvent configuration\n");
+ t_pbc pbc;
+ set_pbc(&pbc, ePBC, box);
+ replicateSolventBox(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt,
+ box_solvt, box);
+ if (ePBC != epbcNONE)
{
- n_box[i] = 1;
- while (n_box[i]*box_solvt[i][i] < box[i][i])
- {
- n_box[i]++;
- }
- nmol *= n_box[i];
+ removeSolventBoxOverlap(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, pbc);
}
- fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
- n_box[XX], n_box[YY], n_box[ZZ]);
-
- /* realloc atoms_solvt for the new solvent configuration */
- srenew(atoms_solvt->resinfo, atoms_solvt->nres*nmol);
- srenew(atoms_solvt->atomname, atoms_solvt->nr*nmol);
- srenew(atoms_solvt->atom, atoms_solvt->nr*nmol);
- srenew(x_solvt, atoms_solvt->nr*nmol);
- if (v_solvt)
+ if (atoms->nr > 0)
{
- srenew(v_solvt, atoms_solvt->nr*nmol);
+ removeSoluteOverlap(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, pbc,
+ atoms->nr, *x, exclusionDistances, rshell);
}
- srenew(exclusionDistances_solvt, atoms_solvt->nr*nmol);
-
- /* generate a new solvent configuration */
- make_new_conformation(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, box_solvt, n_box);
-#ifdef DEBUG
- print_stat(x_solvt, atoms_solvt->nr, box_solvt);
-#endif
+ if (max_sol > 0 && atoms_solvt->nres > max_sol)
+ {
+ const int numberToRemove = atoms_solvt->nres - max_sol;
+ removeExtraSolventMolecules(atoms_solvt, x_solvt, v_solvt, numberToRemove);
+ }
-#ifdef DEBUG
- print_stat(x_solvt, atoms_solvt->nr, box_solvt);
-#endif
/* Sort the solvent mixture, not the protein... */
sort_molecule(&atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt);
- /* add the two configurations */
- onr = atoms->nr;
- onres = atoms->nres;
- add_conf(atoms, x, v, exclusionDistances, TRUE, ePBC, box, FALSE,
- atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, TRUE, rshell, max_sol, oenv);
- *atoms_added = atoms->nr-onr;
- *residues_added = atoms->nres-onres;
+ // Merge the two configurations.
+ srenew(*x, atoms->nr + atoms_solvt->nr);
+ if (v != NULL)
+ {
+ srenew(*v, atoms->nr + atoms_solvt->nr);
+ }
+ for (int i = 0; i < atoms_solvt->nr; ++i)
+ {
+ const int index = atoms->nr + i;
+ copy_rvec(x_solvt[i], (*x)[index]);
+ if (v != NULL)
+ {
+ copy_rvec(v_solvt[i], (*v)[index]);
+ }
+ }
+ {
+ gmx::AtomsBuilder builder(atoms);
+ builder.mergeAtoms(*atoms_solvt);
+ }
+ fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
+ atoms_solvt->nr, atoms_solvt->nres);
sfree(x_solvt);
+ sfree(v_solvt);
+ sfree(exclusionDistances);
sfree(exclusionDistances_solvt);
done_atom(atoms_solvt);
sfree(atoms_solvt);
-
- fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
- *atoms_added, *residues_added);
}
static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
gmx_atomprop_t aps)
{
- #define TEMP_FILENM "temp.top"
FILE *fpin, *fpout;
char buf[STRLEN], buf2[STRLEN], *temp;
const char *topinout;
topinout = ftp2fn(efTOP, NFILE, fnm);
if (ftp2bSet(efTOP, NFILE, fnm) )
{
+ char temporary_filename[STRLEN];
+ strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
+
fprintf(stderr, "Processing topology\n");
fpin = gmx_ffopen(topinout, "r");
- fpout = gmx_ffopen(TEMP_FILENM, "w");
+ gmx_tmpnam(temporary_filename);
+ fpout = gmx_ffopen(temporary_filename, "w");
line = 0;
bSystem = bMolecules = FALSE;
while (fgets(buf, STRLEN, fpin))
/* use gmx_ffopen to generate backup of topinout */
fpout = gmx_ffopen(topinout, "w");
gmx_ffclose(fpout);
- rename(TEMP_FILENM, topinout);
+ rename(temporary_filename, topinout);
}
- #undef TEMP_FILENM
}
int gmx_solvate(int argc, char *argv[])
"into the box. This can create a void that can cause problems later.",
"Choose your volume wisely.[PAR]",
- "The program can optionally rotate the solute molecule to align the",
- "longest molecule axis along a box edge. This way the amount of solvent",
- "molecules necessary is reduced.",
- "It should be kept in mind that this only works for",
- "short simulations, as e.g. an alpha-helical peptide in solution can ",
- "rotate over 90 degrees, within 500 ps. In general it is therefore ",
- "better to make a more or less cubic box.[PAR]",
-
"Setting [TT]-shell[tt] larger than zero will place a layer of water of",
"the specified thickness (nm) around the solute. Hint: it is a good",
"idea to put the protein in the center of a box first (using [gmx-editconf]).",
/* parameter data */
gmx_bool bProt, bBox;
const char *conf_prot, *confout;
- real *exclusionDistances = NULL;
gmx_atomprop_t aps;
/* protein configuration data */
int ePBC = -1;
matrix box;
- /* other data types */
- int atoms_added, residues_added;
-
t_filenm fnm[] = {
{ efSTX, "-cp", "protein", ffOPTRD },
{ efSTX, "-cs", "spc216", ffLIBRD},
"Keep velocities from input solute and solvent" },
};
- if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
+ if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
asize(desc), desc, asize(bugs), bugs, &oenv))
{
return 0;
/* Generate a solute configuration */
conf_prot = opt2fn("-cp", NFILE, fnm);
title = readConformation(conf_prot, atoms, &x,
- bReadV ? &v : NULL, &ePBC, box);
- exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
-
+ bReadV ? &v : NULL, &ePBC, box, "solute");
if (bReadV && !v)
{
fprintf(stderr, "Note: no velocities found\n");
"or give explicit -box command line option");
}
- add_solv(solventFileName, atoms, &x, v ? &v : NULL, &exclusionDistances, ePBC, box,
- aps, defaultDistance, scaleFactor, &atoms_added, &residues_added, r_shell, max_sol,
- oenv);
+ add_solv(solventFileName, atoms, &x, v ? &v : NULL, ePBC, box,
+ aps, defaultDistance, scaleFactor, r_shell, max_sol);
/* write new configuration 1 to file confout */
confout = ftp2fn(efSTO, NFILE, fnm);
if (bProt)
{
write_sto_conf(confout, title, atoms, x, v, ePBC, box);
- /* print box sizes and box type to stderr */
- fprintf(stderr, "%s\n", title);
}
else
{
update_top(atoms, box, NFILE, fnm, aps);
gmx_atomprop_destroy(aps);
- sfree(exclusionDistances);
sfree(x);
sfree(v);
done_atom(atoms);
sfree(atoms);
+ sfree(title);
+ output_env_done(oenv);
+ done_filenms(NFILE, fnm);
return 0;
}
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "membed.h"
#include <signal.h>
#include <stdlib.h>
-#include "typedefs.h"
-#include "types/commrec.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
-#include "vec.h"
-#include "macros.h"
-#include "main.h"
-#include "gromacs/fileio/futil.h"
+
#include "gromacs/essentialdynamics/edsam.h"
-#include "index.h"
-#include "physics.h"
-#include "names.h"
-#include "mtop_util.h"
#include "gromacs/fileio/tpxio.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/readinp.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/cstringutil.h"
-#include "membed.h"
-#include "pbc.h"
-#include "readinp.h"
-#include "gromacs/gmxpreprocess/readir.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
/* information about scaling center */
typedef struct {
/* Write a topology where the number of molecules is correct for the system after embedding */
static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
{
- #define TEMP_FILENM "temp.top"
- int bMolecules = 0;
+ int bMolecules = 0;
FILE *fpin, *fpout;
char buf[STRLEN], buf2[STRLEN], *temp;
int i, *nmol_rm, nmol, line;
+ char temporary_filename[STRLEN];
fpin = gmx_ffopen(topfile, "r");
- fpout = gmx_ffopen(TEMP_FILENM, "w");
+ strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
+ gmx_tmpnam(temporary_filename);
+ fpout = gmx_ffopen(temporary_filename, "w");
snew(nmol_rm, mtop->nmoltype);
for (i = 0; i < rm_p->nr; i++)
/* use gmx_ffopen to generate backup of topinout */
fpout = gmx_ffopen(topfile, "w");
gmx_ffclose(fpout);
- rename(TEMP_FILENM, topfile);
- #undef TEMP_FILENM
+ rename(temporary_filename, topfile);
}
void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x)
resize(membed->r_ins, x, membed->pos_ins, membed->fac);
}
+/* We would like gn to be const as well, but C doesn't allow this */
+/* TODO this is utility functionality (search for the index of a
+ string in a collection), so should be refactored and located more
+ centrally. Also, it nearly duplicates the same string in readir.c */
+static int search_string(const char *s, int ng, char *gn[])
+{
+ int i;
+
+ for (i = 0; (i < ng); i++)
+ {
+ if (gmx_strcasecmp(s, gn[i]) == 0)
+ {
+ return i;
+ }
+ }
+
+ gmx_fatal(FARGS,
+ "Group %s selected for embedding was not found in the index file.\n"
+ "Group names must match either [moleculetype] names or custom index group\n"
+ "names, in which case you must supply an index file to the '-n' option\n"
+ "of grompp.",
+ s);
+
+ return -1;
+}
+
gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
{
get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
&maxwarn, &pieces, &bALLOW_ASYMMETRY);
- tpr_version = get_tpr_version(ftp2fn(efTPX, nfile, fnm));
+ tpr_version = get_tpr_version(ftp2fn(efTPR, nfile, fnm));
if (tpr_version < membed_version)
{
gmx_fatal(FARGS, "Version of *.tpr file to old (%d). "
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+
+#include "gmxpre.h"
+
+#include "config.h"
+
+#include <assert.h>
#include <signal.h>
#include <stdlib.h>
-#ifdef HAVE_UNISTD_H
-#include <unistd.h>
-#endif
#include <string.h>
-#include <assert.h>
-#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
-#include "copyrite.h"
-#include "force.h"
-#include "mdrun.h"
-#include "md_logging.h"
-#include "md_support.h"
-#include "network.h"
-#include "names.h"
-#include "disre.h"
-#include "orires.h"
-#include "pme.h"
-#include "mdatoms.h"
-#include "repl_ex.h"
-#include "deform.h"
-#include "qmmm.h"
-#include "domdec.h"
-#include "coulomb.h"
-#include "constr.h"
-#include "mvdata.h"
-#include "checkpoint.h"
-#include "mtop_util.h"
-#include "sighandler.h"
-#include "txtdump.h"
-#include "gmx_detect_hardware.h"
-#include "gmx_omp_nthreads.h"
-#include "gromacs/gmxpreprocess/calc_verletbuf.h"
-#include "gmx_fatal_collective.h"
-#include "membed.h"
-#include "macros.h"
-#include "gmx_thread_affinity.h"
-#include "inputrec.h"
+#include <algorithm>
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/essentialdynamics/edsam.h"
+#include "gromacs/ewald/pme.h"
#include "gromacs/fileio/tpxio.h"
-#include "gromacs/mdlib/nbnxn_search.h"
+#include "gromacs/gmxlib/gpu_utils/gpu_utils.h"
+#include "gromacs/legacyheaders/checkpoint.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/disre.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/gmx_detect_hardware.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/gmx_thread_affinity.h"
+#include "gromacs/legacyheaders/inputrec.h"
+#include "gromacs/legacyheaders/main.h"
+#include "gromacs/legacyheaders/md_logging.h"
+#include "gromacs/legacyheaders/md_support.h"
+#include "gromacs/legacyheaders/mdatoms.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/oenv.h"
+#include "gromacs/legacyheaders/orires.h"
+#include "gromacs/legacyheaders/qmmm.h"
+#include "gromacs/legacyheaders/sighandler.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/calc_verletbuf.h"
#include "gromacs/mdlib/nbnxn_consts.h"
-#include "gromacs/timing/wallcycle.h"
-#include "gromacs/utility/gmxmpi.h"
-#include "gromacs/utility/gmxomp.h"
-#include "gromacs/swap/swapcoords.h"
-#include "gromacs/essentialdynamics/edsam.h"
+#include "gromacs/mdlib/nbnxn_search.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/pulling/pull_rotation.h"
+#include "gromacs/swap/swapcoords.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "deform.h"
+#include "membed.h"
+#include "repl_ex.h"
#ifdef GMX_FAHCORE
#include "corewrap.h"
#endif
-#include "gpu_utils.h"
-#include "nbnxn_cuda_data_mgmt.h"
-
typedef struct {
gmx_integrator_t *func;
} gmx_intp_t;
real pforce;
real cpt_period;
real max_hours;
- const char *deviceOptions;
int imdport;
unsigned long Flags;
};
mc.nbpu_opt, mc.nstlist_cmdline,
mc.nsteps_cmdline, mc.nstepout, mc.resetstep,
mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_nex, mc.repl_ex_seed, mc.pforce,
- mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.imdport, mc.Flags);
+ mc.cpt_period, mc.max_hours, mc.imdport, mc.Flags);
}
/* called by mdrunner() to start a specific number of threads (including
int nstepout, int resetstep,
int nmultisim, int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
real pforce, real cpt_period, real max_hours,
- const char *deviceOptions, unsigned long Flags)
+ unsigned long Flags)
{
int ret;
struct mdrunner_arglist *mda;
mda->pforce = pforce;
mda->cpt_period = cpt_period;
mda->max_hours = max_hours;
- mda->deviceOptions = deviceOptions;
mda->Flags = Flags;
/* now spawn new threads that start mdrunner_start_fn(), while
else if (hw_opt->nthreads_omp > 0)
{
/* Here we could oversubscribe, when we do, we issue a warning later */
- nthreads_tmpi = max(1, nthreads_tot/hw_opt->nthreads_omp);
+ nthreads_tmpi = std::max(1, nthreads_tot/hw_opt->nthreads_omp);
}
else
{
{
int nthreads_hw, nthreads_tot_max, nthreads_tmpi, nthreads_new, ngpu;
int min_atoms_per_mpi_thread;
- char *env;
- char sbuf[STRLEN];
gmx_bool bCanUseGPU;
if (hw_opt->nthreads_tmpi > 0)
}
bCanUseGPU = (inputrec->cutoff_scheme == ecutsVERLET &&
- hwinfo->gpu_info.ncuda_dev_compatible > 0);
+ hwinfo->gpu_info.n_dev_compatible > 0);
+
if (bCanUseGPU)
{
- ngpu = hwinfo->gpu_info.ncuda_dev_compatible;
+ ngpu = hwinfo->gpu_info.n_dev_compatible;
}
else
{
{
/* the thread number was chosen automatically, but there are too many
threads (too few atoms per thread) */
- nthreads_new = max(1, mtop->natoms/min_atoms_per_mpi_thread);
+ nthreads_new = std::max(1, mtop->natoms/min_atoms_per_mpi_thread);
/* Avoid partial use of Hyper-Threading */
if (gmx_cpuid_x86_smt(hwinfo->cpuid_info) == GMX_CPUID_X86_SMT_ENABLED &&
/* We determine the extra cost of the non-bonded kernels compared to
* a reference nstlist value of 10 (which is the default in grompp).
*/
-static const int nbnxn_reference_nstlist = 10;
+static const int nbnxnReferenceNstlist = 10;
/* The values to try when switching */
const int nstlist_try[] = { 20, 25, 40 };
#define NNSTL sizeof(nstlist_try)/sizeof(nstlist_try[0])
float listfac_ok, listfac_max;
int nstlist_orig, nstlist_prev;
verletbuf_list_setup_t ls;
- real rlist_nstlist10, rlist_inc, rlist_ok, rlist_max;
+ real rlistWithReferenceNstlist, rlist_inc, rlist_ok, rlist_max;
real rlist_new, rlist_prev;
- int nstlist_ind = 0;
+ size_t nstlist_ind = 0;
t_state state_tmp;
gmx_bool bBox, bDD, bCont;
const char *nstl_gpu = "\nFor optimal performance with a GPU nstlist (now %d) should be larger.\nThe optimum depends on your CPU and GPU resources.\nYou might want to try several nstlist values.\n";
const char *box_err = "Can not increase nstlist because the box is too small";
const char *dd_err = "Can not increase nstlist because of domain decomposition limitations";
char buf[STRLEN];
+ const float oneThird = 1.0f / 3.0f;
if (nstlist_cmdline <= 0)
{
verletbuf_get_list_setup(bGPU, &ls);
/* Allow rlist to make the list a given factor larger than the list
- * would be with nstlist=10.
+ * would be with the reference value for nstlist (10).
*/
nstlist_prev = ir->nstlist;
- ir->nstlist = 10;
+ ir->nstlist = nbnxnReferenceNstlist;
calc_verlet_buffer_size(mtop, det(box), ir, -1, &ls, NULL,
- &rlist_nstlist10);
+ &rlistWithReferenceNstlist);
ir->nstlist = nstlist_prev;
/* Determine the pair list size increase due to zero interactions */
rlist_inc = nbnxn_get_rlist_effective_inc(ls.cluster_size_j,
mtop->natoms/det(box));
- rlist_ok = (rlist_nstlist10 + rlist_inc)*pow(listfac_ok, 1.0/3.0) - rlist_inc;
- rlist_max = (rlist_nstlist10 + rlist_inc)*pow(listfac_max, 1.0/3.0) - rlist_inc;
+ rlist_ok = (rlistWithReferenceNstlist + rlist_inc)*pow(listfac_ok, oneThird) - rlist_inc;
+ rlist_max = (rlistWithReferenceNstlist + rlist_inc)*pow(listfac_max, oneThird) - rlist_inc;
if (debug)
{
fprintf(debug, "nstlist tuning: rlist_inc %.3f rlist_ok %.3f rlist_max %.3f\n",
}
}
-static void convert_to_verlet_scheme(FILE *fplog,
- t_inputrec *ir,
- gmx_mtop_t *mtop, real box_vol)
-{
- char *conv_mesg = "Converting input file with group cut-off scheme to the Verlet cut-off scheme";
-
- md_print_warn(NULL, fplog, "%s\n", conv_mesg);
-
- ir->cutoff_scheme = ecutsVERLET;
- ir->verletbuf_tol = 0.005;
-
- if (ir->rcoulomb != ir->rvdw)
- {
- gmx_fatal(FARGS, "The VdW and Coulomb cut-offs are different, whereas the Verlet scheme only supports equal cut-offs");
- }
-
- if (ir->vdwtype == evdwUSER || EEL_USER(ir->coulombtype))
- {
- gmx_fatal(FARGS, "User non-bonded potentials are not (yet) supported with the Verlet scheme");
- }
- else if (ir_vdw_switched(ir) || ir_coulomb_switched(ir))
- {
- if (ir_vdw_switched(ir) && ir->vdw_modifier == eintmodNONE)
- {
- ir->vdwtype = evdwCUT;
-
- switch (ir->vdwtype)
- {
- case evdwSHIFT: ir->vdw_modifier = eintmodFORCESWITCH; break;
- case evdwSWITCH: ir->vdw_modifier = eintmodPOTSWITCH; break;
- default: gmx_fatal(FARGS, "The Verlet scheme does not support Van der Waals interactions of type '%s'", evdw_names[ir->vdwtype]);
- }
- }
- if (ir_coulomb_switched(ir) && ir->coulomb_modifier == eintmodNONE)
- {
- if (EEL_FULL(ir->coulombtype))
- {
- /* With full electrostatic only PME can be switched */
- ir->coulombtype = eelPME;
- ir->coulomb_modifier = eintmodPOTSHIFT;
- }
- else
- {
- md_print_warn(NULL, fplog, "NOTE: Replacing %s electrostatics with reaction-field with epsilon-rf=inf\n", eel_names[ir->coulombtype]);
- ir->coulombtype = eelRF;
- ir->epsilon_rf = 0.0;
- ir->coulomb_modifier = eintmodPOTSHIFT;
- }
- }
-
- /* We set the pair energy error tolerance to a small number.
- * Note that this is only for testing. For production the user
- * should think about this and set the mdp options.
- */
- ir->verletbuf_tol = 1e-4;
- }
-
- if (inputrec2nboundeddim(ir) != 3)
- {
- gmx_fatal(FARGS, "Can only convert old tpr files to the Verlet cut-off scheme with 3D pbc");
- }
-
- if (ir->efep != efepNO || ir->implicit_solvent != eisNO)
- {
- gmx_fatal(FARGS, "Will not convert old tpr files to the Verlet cut-off scheme with free-energy calculations or implicit solvent");
- }
-
- if (EI_DYNAMICS(ir->eI) && !(EI_MD(ir->eI) && ir->etc == etcNO))
- {
- verletbuf_list_setup_t ls;
-
- verletbuf_get_list_setup(FALSE, &ls);
- calc_verlet_buffer_size(mtop, box_vol, ir, -1, &ls, NULL, &ir->rlist);
- }
- else
- {
- real rlist_fac;
-
- if (EI_MD(ir->eI))
- {
- rlist_fac = 1 + verlet_buffer_ratio_NVE_T0;
- }
- else
- {
- rlist_fac = 1 + verlet_buffer_ratio_nodynamics;
- }
- ir->verletbuf_tol = -1;
- ir->rlist = rlist_fac*max(ir->rvdw, ir->rcoulomb);
- }
-
- gmx_mtop_remove_chargegroups(mtop);
-}
-
static void print_hw_opt(FILE *fp, const gmx_hw_opt_t *hw_opt)
{
fprintf(fp, "hw_opt: nt %d ntmpi %d ntomp %d ntomp_pme %d gpu_id '%s'\n",
#ifndef GMX_THREAD_MPI
if (hw_opt->nthreads_tot > 0)
{
- gmx_fatal(FARGS, "Setting the total number of threads is only supported with thread-MPI and Gromacs was compiled without thread-MPI");
+ gmx_fatal(FARGS, "Setting the total number of threads is only supported with thread-MPI and GROMACS was compiled without thread-MPI");
}
if (hw_opt->nthreads_tmpi > 0)
{
- gmx_fatal(FARGS, "Setting the number of thread-MPI threads is only supported with thread-MPI and Gromacs was compiled without thread-MPI");
+ gmx_fatal(FARGS, "Setting the number of thread-MPI threads is only supported with thread-MPI and GROMACS was compiled without thread-MPI");
}
#endif
#ifndef GMX_OPENMP
if (hw_opt->nthreads_omp > 1)
{
- gmx_fatal(FARGS, "More than 1 OpenMP thread requested, but Gromacs was compiled without OpenMP support");
+ gmx_fatal(FARGS, "More than 1 OpenMP thread requested, but GROMACS was compiled without OpenMP support");
}
hw_opt->nthreads_omp = 1;
#endif
#ifndef GMX_OPENMP
if (hw_opt->nthreads_omp > 1)
{
- gmx_fatal(FARGS, "OpenMP threads are requested, but Gromacs was compiled without OpenMP support");
+ gmx_fatal(FARGS, "OpenMP threads are requested, but GROMACS was compiled without OpenMP support");
}
#endif
gmx_parse_gpu_ids(&hw_opt->gpu_opt);
#ifdef GMX_THREAD_MPI
- if (hw_opt->gpu_opt.ncuda_dev_use > 0 && hw_opt->nthreads_tmpi == 0)
+ if (hw_opt->gpu_opt.n_dev_use > 0
+ &&
+ hw_opt->nthreads_tmpi == 0)
{
/* Set the number of MPI threads equal to the number of GPUs */
- hw_opt->nthreads_tmpi = hw_opt->gpu_opt.ncuda_dev_use;
+ hw_opt->nthreads_tmpi = hw_opt->gpu_opt.n_dev_use;
if (hw_opt->nthreads_tot > 0 &&
hw_opt->nthreads_tmpi > hw_opt->nthreads_tot)
t_inputrec *ir,
const t_commrec *cr)
{
- char sbuf[STEPSTRSIZE];
-
assert(ir);
assert(cr);
/* override with anything else than the default -2 */
if (nsteps_cmdline > -2)
{
- char stmp[STRLEN];
+ char sbuf_steps[STEPSTRSIZE];
+ char sbuf_msg[STRLEN];
ir->nsteps = nsteps_cmdline;
if (EI_DYNAMICS(ir->eI) && nsteps_cmdline != -1)
{
- sprintf(stmp, "Overriding nsteps with value passed on the command line: %s steps, %.3g ps",
- gmx_step_str(nsteps_cmdline, sbuf),
+ sprintf(sbuf_msg, "Overriding nsteps with value passed on the command line: %s steps, %.3g ps",
+ gmx_step_str(nsteps_cmdline, sbuf_steps),
fabs(nsteps_cmdline*ir->delta_t));
}
else
{
- sprintf(stmp, "Overriding nsteps with value passed on the command line: %s steps",
- gmx_step_str(nsteps_cmdline, sbuf));
+ sprintf(sbuf_msg, "Overriding nsteps with value passed on the command line: %s steps",
+ gmx_step_str(nsteps_cmdline, sbuf_steps));
}
- md_print_warn(cr, fplog, "%s\n", stmp);
+ md_print_warn(cr, fplog, "%s\n", sbuf_msg);
}
-}
-
-/* Frees GPU memory and destroys the CUDA context.
- *
- * Note that this function needs to be called even if GPUs are not used
- * in this run because the PME ranks have no knowledge of whether GPUs
- * are used or not, but all ranks need to enter the barrier below.
- */
-static void free_gpu_resources(const t_forcerec *fr,
- const t_commrec *cr)
-{
- gmx_bool bIsPPrankUsingGPU;
- char gpu_err_str[STRLEN];
-
- bIsPPrankUsingGPU = (cr->duty & DUTY_PP) && fr->nbv != NULL && fr->nbv->bUseGPU;
-
- if (bIsPPrankUsingGPU)
+ else if (nsteps_cmdline < -2)
{
- /* free nbnxn data in GPU memory */
- nbnxn_cuda_free(fr->nbv->cu_nbv);
-
- /* With tMPI we need to wait for all ranks to finish deallocation before
- * destroying the context in free_gpu() as some ranks may be sharing
- * GPU and context.
- * Note: as only PP ranks need to free GPU resources, so it is safe to
- * not call the barrier on PME ranks.
- */
-#ifdef GMX_THREAD_MPI
- if (PAR(cr))
- {
- gmx_barrier(cr);
- }
-#endif /* GMX_THREAD_MPI */
-
- /* uninitialize GPU (by destroying the context) */
- if (!free_gpu(gpu_err_str))
- {
- gmx_warning("On rank %d failed to free GPU #%d: %s",
- cr->nodeid, get_current_gpu_device_id(), gpu_err_str);
- }
+ gmx_fatal(FARGS, "Invalid nsteps value passed on the command line: %d",
+ nsteps_cmdline);
}
+ /* Do nothing if nsteps_cmdline == -2 */
}
int mdrunner(gmx_hw_opt_t *hw_opt,
gmx_int64_t nsteps_cmdline, int nstepout, int resetstep,
int gmx_unused nmultisim, int repl_ex_nst, int repl_ex_nex,
int repl_ex_seed, real pforce, real cpt_period, real max_hours,
- const char *deviceOptions, int imdport, unsigned long Flags)
+ int imdport, unsigned long Flags)
{
gmx_bool bForceUseGPU, bTryUseGPU, bRerunMD, bCantUseGPU;
- double nodetime = 0, realtime;
t_inputrec *inputrec;
t_state *state = NULL;
matrix box;
gmx_ddbox_t ddbox = {0};
int npme_major, npme_minor;
- real tmpr1, tmpr2;
t_nrnb *nrnb;
gmx_mtop_t *mtop = NULL;
t_mdatoms *mdatoms = NULL;
t_fcdata *fcd = NULL;
real ewaldcoeff_q = 0;
real ewaldcoeff_lj = 0;
- gmx_pme_t *pmedata = NULL;
+ struct gmx_pme_t **pmedata = NULL;
gmx_vsite_t *vsite = NULL;
gmx_constr_t constr;
- int i, m, nChargePerturbed = -1, nTypePerturbed = 0, status, nalloc;
- char *gro;
+ int nChargePerturbed = -1, nTypePerturbed = 0, status;
gmx_wallcycle_t wcycle;
gmx_bool bReadEkin;
- int list;
gmx_walltime_accounting_t walltime_accounting = NULL;
int rc;
gmx_int64_t reset_counters;
gmx_edsam_t ed = NULL;
- t_commrec *cr_old = cr;
int nthreads_pme = 1;
int nthreads_pp = 1;
gmx_membed_t membed = NULL;
if (SIMMASTER(cr))
{
/* Read (nearly) all data required for the simulation */
- read_tpx_state(ftp2fn(efTPX, nfile, fnm), inputrec, state, NULL, mtop);
-
- if (inputrec->cutoff_scheme != ecutsVERLET &&
- ((Flags & MD_TESTVERLET) || getenv("GMX_VERLET_SCHEME") != NULL))
- {
- convert_to_verlet_scheme(fplog, inputrec, mtop, det(state->box));
- }
+ read_tpx_state(ftp2fn(efTPR, nfile, fnm), inputrec, state, NULL, mtop);
if (inputrec->cutoff_scheme == ecutsVERLET)
{
/* Here the master rank decides if all ranks will use GPUs */
- bUseGPU = (hwinfo->gpu_info.ncuda_dev_compatible > 0 ||
+ bUseGPU = (hwinfo->gpu_info.n_dev_compatible > 0 ||
getenv("GMX_EMULATE_GPU") != NULL);
/* TODO add GPU kernels for this and replace this check by:
gmx_fatal(FARGS, "Can not set nstlist with the group cut-off scheme");
}
- if (hwinfo->gpu_info.ncuda_dev_compatible > 0)
+ if (hwinfo->gpu_info.n_dev_compatible > 0)
{
md_print_warn(cr, fplog,
"NOTE: GPU(s) found, but the current simulation can not use GPUs\n"
- " To use a GPU, set the mdp option: cutoff-scheme = Verlet\n"
- " (for quick performance testing you can use the -testverlet option)\n");
+ " To use a GPU, set the mdp option: cutoff-scheme = Verlet\n");
}
if (bForceUseGPU)
}
}
- /* Check for externally set OpenMP affinity and turn off internal
- * pinning if any is found. We need to do this check early to tell
- * thread-MPI whether it should do pinning when spawning threads.
- * TODO: the above no longer holds, we should move these checks down
- */
- gmx_omp_check_thread_affinity(fplog, cr, hw_opt);
-
/* Check and update the hardware options for internal consistency */
check_and_update_hw_opt_1(hw_opt, SIMMASTER(cr));
+ /* Early check for externally set process affinity. */
+ gmx_check_thread_affinity_set(fplog, cr,
+ hw_opt, hwinfo->nthreads_hw_avail, FALSE);
if (SIMMASTER(cr))
{
-#ifdef GMX_THREAD_MPI
- /* Early check for externally set process affinity.
- * With thread-MPI this is needed as pinning might get turned off,
- * which needs to be known before starting thread-MPI.
- * With thread-MPI hw_opt is processed here on the master rank
- * and passed to the other ranks later, so we only do this on master.
- */
- gmx_check_thread_affinity_set(fplog,
- NULL,
- hw_opt, hwinfo->nthreads_hw_avail, FALSE);
-#endif
#ifdef GMX_THREAD_MPI
if (cr->npmenodes > 0 && hw_opt->nthreads_tmpi <= 0)
if (hw_opt->nthreads_tmpi > 1)
{
+ t_commrec *cr_old = cr;
/* now start the threads. */
cr = mdrunner_start_threads(hw_opt, fplog, cr_old, nfile, fnm,
oenv, bVerbose, bCompact, nstglobalcomm,
nbpu_opt, nstlist_cmdline,
nsteps_cmdline, nstepout, resetstep, nmultisim,
repl_ex_nst, repl_ex_nex, repl_ex_seed, pforce,
- cpt_period, max_hours, deviceOptions,
+ cpt_period, max_hours,
Flags);
/* the main thread continues here with a new cr. We don't deallocate
the old cr because other threads may still be reading it. */
{
/* now broadcast everything to the non-master nodes/threads: */
init_parallel(cr, inputrec, mtop);
+
+ /* The master rank decided on the use of GPUs,
+ * broadcast this information to all ranks.
+ */
+ gmx_bcast_sim(sizeof(bUseGPU), &bUseGPU, cr);
}
+
if (fplog != NULL)
{
pr_inputrec(fplog, 0, "Input Parameters", inputrec, FALSE);
"but %s was not started through mpirun/mpiexec or only one rank was requested through mpirun/mpiexec"
#endif
#endif
- , ShortProgram()
+ , output_env_get_program_display_name(oenv)
);
}
cr->npmenodes = 0;
}
+ if (bUseGPU && cr->npmenodes < 0)
+ {
+ /* With GPUs we don't automatically use PME-only ranks. PME ranks can
+ * improve performance with many threads per GPU, since our OpenMP
+ * scaling is bad, but it's difficult to automate the setup.
+ */
+ cr->npmenodes = 0;
+ }
+
#ifdef GMX_FAHCORE
if (MASTER(cr))
{
}
}
- if (((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
-#ifdef GMX_THREAD_MPI
- /* With thread MPI only the master node/thread exists in mdrun.c,
- * therefore non-master nodes need to open the "seppot" log file here.
- */
- || (!MASTER(cr) && (Flags & MD_SEPPOT))
-#endif
- )
+ if (MASTER(cr) && (Flags & MD_APPENDFILES))
{
- gmx_log_open(ftp2fn(efLOG, nfile, fnm), cr, !(Flags & MD_SEPPOT),
+ gmx_log_open(ftp2fn(efLOG, nfile, fnm), cr,
Flags, &fplog);
}
if (MULTISIM(cr))
{
md_print_info(cr, fplog,
- "This is simulation %d out of %d running as a composite Gromacs\n"
+ "This is simulation %d out of %d running as a composite GROMACS\n"
"multi-simulation job. Setup for this simulation:\n\n",
cr->ms->sim, cr->ms->nsim);
}
(cr->duty & DUTY_PP) == 0,
inputrec->cutoff_scheme == ecutsVERLET);
- if (PAR(cr))
- {
- /* The master rank decided on the use of GPUs,
- * broadcast this information to all ranks.
- */
- gmx_bcast_sim(sizeof(bUseGPU), &bUseGPU, cr);
- }
+#ifndef NDEBUG
+ if (integrator[inputrec->eI].func != do_tpi &&
+ inputrec->cutoff_scheme == ecutsVERLET)
+ {
+ gmx_feenableexcept();
+ }
+#endif
+
if (bUseGPU)
{
- if (cr->npmenodes == -1)
- {
- /* Don't automatically use PME-only nodes with GPUs */
- cr->npmenodes = 0;
- }
-
/* Select GPU id's to use */
gmx_select_gpu_ids(fplog, cr, &hwinfo->gpu_info, bForceUseGPU,
&hw_opt->gpu_opt);
else
{
/* Ignore (potentially) manually selected GPUs */
- hw_opt->gpu_opt.ncuda_dev_use = 0;
+ hw_opt->gpu_opt.n_dev_use = 0;
}
/* check consistency across ranks of things like SIMD
/*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
"nofile","nofile","nofile","nofile",FALSE,pforce);
*/
- fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
/* Initialize QM-MM */
if (fr->bQMMM)
/* Assumes uniform use of the number of OpenMP threads */
walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntDefault));
- if (inputrec->ePull != epullNO)
+ if (inputrec->bPull)
{
/* Initialize pull code */
init_pull(fplog, inputrec, nfile, fnm, mtop, cr, oenv, inputrec->fepvals->init_lambda,
if (DOMAINDECOMP(cr))
{
+ GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP");
dd_init_bondeds(fplog, cr->dd, mtop, vsite, inputrec,
Flags & MD_DDBONDCHECK, fr->cginfo_mb);
repl_ex_nst, repl_ex_nex, repl_ex_seed,
membed,
cpt_period, max_hours,
- deviceOptions,
imdport,
Flags,
walltime_accounting);
- if (inputrec->ePull != epullNO)
+ if (inputrec->bPull)
{
finish_pull(inputrec->pull);
}
}
else
{
+ GMX_RELEASE_ASSERT(pmedata, "pmedata was NULL while cr->duty was not DUTY_PP");
/* do PME only */
walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntPME));
gmx_pmeonly(*pmedata, cr, nrnb, wcycle, walltime_accounting, ewaldcoeff_q, ewaldcoeff_lj, inputrec);
*/
finish_run(fplog, cr,
inputrec, nrnb, wcycle, walltime_accounting,
- fr != NULL && fr->nbv != NULL && fr->nbv->bUseGPU ?
- nbnxn_cuda_get_timings(fr->nbv->cu_nbv) : NULL,
+ fr ? fr->nbv : NULL,
EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
/* Free GPU memory and context */
- free_gpu_resources(fr, cr);
+ free_gpu_resources(fr, cr, &hwinfo->gpu_info, fr ? fr->gpu_opt : NULL);
if (opt2bSet("-membed", nfile, fnm))
{
rc = (int)gmx_get_stop_condition();
+ done_ed(&ed);
+
#ifdef GMX_THREAD_MPI
/* we need to join all threads. The sub-threads join when they
exit this function, but the master thread needs to be told to
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff