2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
7 * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/ewald/pme.h"
47 #include "gromacs/gpu_utils/hostallocator.h"
48 #include "gromacs/math/functions.h"
49 #include "gromacs/mdlib/gmx_omp_nthreads.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/mdtypes/mdatom.h"
53 #include "gromacs/topology/mtop_lookup.h"
54 #include "gromacs/topology/mtop_util.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/exceptions.h"
57 #include "gromacs/utility/smalloc.h"
59 #define ALMOST_ZERO 1e-30
64 MDAtoms::MDAtoms() : mdatoms_(nullptr) {}
68 if (mdatoms_ == nullptr)
72 sfree(mdatoms_->massA);
73 sfree(mdatoms_->massB);
74 sfree(mdatoms_->massT);
75 gmx::AlignedAllocationPolicy::free(mdatoms_->invmass);
76 sfree(mdatoms_->invMassPerDim);
77 sfree(mdatoms_->typeA);
78 sfree(mdatoms_->typeB);
79 /* mdatoms->chargeA and mdatoms->chargeB point at chargeA_.data()
80 * and chargeB_.data() respectively. They get freed automatically. */
81 sfree(mdatoms_->sqrt_c6A);
82 sfree(mdatoms_->sigmaA);
83 sfree(mdatoms_->sigma3A);
84 sfree(mdatoms_->sqrt_c6B);
85 sfree(mdatoms_->sigmaB);
86 sfree(mdatoms_->sigma3B);
87 sfree(mdatoms_->ptype);
89 sfree(mdatoms_->cENER);
90 sfree(mdatoms_->cFREEZE);
91 sfree(mdatoms_->cVCM);
92 sfree(mdatoms_->cORF);
93 sfree(mdatoms_->bPerturbed);
98 void MDAtoms::resizeChargeA(const int newSize)
100 chargeA_.resizeWithPadding(newSize);
101 mdatoms_->chargeA = chargeA_.data();
104 void MDAtoms::resizeChargeB(const int newSize)
106 chargeB_.resizeWithPadding(newSize);
107 mdatoms_->chargeB = chargeB_.data();
110 void MDAtoms::reserveChargeA(const int newCapacity)
112 chargeA_.reserveWithPadding(newCapacity);
113 mdatoms_->chargeA = chargeA_.data();
116 void MDAtoms::reserveChargeB(const int newCapacity)
118 chargeB_.reserveWithPadding(newCapacity);
119 mdatoms_->chargeB = chargeB_.data();
122 std::unique_ptr<MDAtoms> makeMDAtoms(FILE* fp, const gmx_mtop_t& mtop, const t_inputrec& ir, const bool rankHasPmeGpuTask)
124 auto mdAtoms = std::make_unique<MDAtoms>();
125 // GPU transfers may want to use a suitable pinning mode.
126 if (rankHasPmeGpuTask)
128 changePinningPolicy(&mdAtoms->chargeA_, pme_get_pinning_policy());
129 changePinningPolicy(&mdAtoms->chargeB_, pme_get_pinning_policy());
133 mdAtoms->mdatoms_.reset(md);
135 md->nenergrp = mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size();
136 md->bVCMgrps = FALSE;
137 for (int i = 0; i < mtop.natoms; i++)
139 if (getGroupType(mtop.groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i) > 0)
145 /* Determine the total system mass and perturbed atom counts */
146 double totalMassA = 0.0;
147 double totalMassB = 0.0;
149 md->haveVsites = FALSE;
150 gmx_mtop_atomloop_block_t aloop = gmx_mtop_atomloop_block_init(mtop);
153 while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
155 totalMassA += nmol * atom->m;
156 totalMassB += nmol * atom->mB;
158 if (atom->ptype == ParticleType::VSite)
160 md->haveVsites = TRUE;
163 if (ir.efep != FreeEnergyPerturbationType::No && PERTURBED(*atom))
166 if (atom->mB != atom->m)
168 md->nMassPerturbed += nmol;
170 if (atom->qB != atom->q)
172 md->nChargePerturbed += nmol;
174 if (atom->typeB != atom->type)
176 md->nTypePerturbed += nmol;
181 md->tmassA = totalMassA;
182 md->tmassB = totalMassB;
184 if (ir.efep != FreeEnergyPerturbationType::No && fp)
187 "There are %d atoms and %d charges for free energy perturbation\n",
189 md->nChargePerturbed);
192 md->havePartiallyFrozenAtoms = FALSE;
193 for (int g = 0; g < ir.opts.ngfrz; g++)
195 for (int d = YY; d < DIM; d++)
197 if (ir.opts.nFreeze[g][d] != ir.opts.nFreeze[g][XX])
199 md->havePartiallyFrozenAtoms = TRUE;
204 md->bOrires = (gmx_mtop_ftype_count(mtop, F_ORIRES) != 0);
211 void atoms2md(const gmx_mtop_t& mtop,
212 const t_inputrec& inputrec,
214 gmx::ArrayRef<int> index,
216 gmx::MDAtoms* mdAtoms)
219 const t_grpopts* opts;
220 int nthreads gmx_unused;
222 bLJPME = EVDW_PME(inputrec.vdwtype);
224 opts = &inputrec.opts;
226 const SimulationGroups& groups = mtop.groups;
228 auto* md = mdAtoms->mdatoms();
229 /* nindex>=0 indicates DD where we use an index */
236 md->nr = mtop.natoms;
239 if (md->nr > md->nalloc)
241 md->nalloc = over_alloc_dd(md->nr);
243 if (md->nMassPerturbed)
245 srenew(md->massA, md->nalloc);
246 srenew(md->massB, md->nalloc);
248 srenew(md->massT, md->nalloc);
249 /* The SIMD version of the integrator needs this aligned and padded.
250 * The padding needs to be with zeros, which we set later below.
252 gmx::AlignedAllocationPolicy::free(md->invmass);
253 md->invmass = new (gmx::AlignedAllocationPolicy::malloc(
254 (md->nalloc + GMX_REAL_MAX_SIMD_WIDTH) * sizeof(*md->invmass))) real;
255 srenew(md->invMassPerDim, md->nalloc);
256 // TODO eventually we will have vectors and just resize
257 // everything, but for now the semantics of md->nalloc being
258 // the capacity are preserved by keeping vectors within
259 // mdAtoms having the same properties as the other arrays.
260 mdAtoms->reserveChargeA(md->nalloc);
261 mdAtoms->resizeChargeA(md->nr);
262 if (md->nPerturbed > 0)
264 mdAtoms->reserveChargeB(md->nalloc);
265 mdAtoms->resizeChargeB(md->nr);
267 srenew(md->typeA, md->nalloc);
270 srenew(md->typeB, md->nalloc);
274 srenew(md->sqrt_c6A, md->nalloc);
275 srenew(md->sigmaA, md->nalloc);
276 srenew(md->sigma3A, md->nalloc);
279 srenew(md->sqrt_c6B, md->nalloc);
280 srenew(md->sigmaB, md->nalloc);
281 srenew(md->sigma3B, md->nalloc);
284 srenew(md->ptype, md->nalloc);
287 srenew(md->cTC, md->nalloc);
288 /* We always copy cTC with domain decomposition */
290 srenew(md->cENER, md->nalloc);
291 if (inputrecFrozenAtoms(&inputrec))
293 srenew(md->cFREEZE, md->nalloc);
297 srenew(md->cVCM, md->nalloc);
301 srenew(md->cORF, md->nalloc);
305 srenew(md->bPerturbed, md->nalloc);
308 /* Note that these user t_mdatoms array pointers are NULL
309 * when there is only one group present.
310 * Therefore, when adding code, the user should use something like:
311 * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
313 if (!mtop.groups.groupNumbers[SimulationAtomGroupType::User1].empty())
315 srenew(md->cU1, md->nalloc);
317 if (!mtop.groups.groupNumbers[SimulationAtomGroupType::User2].empty())
319 srenew(md->cU2, md->nalloc);
325 nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Default);
326 #pragma omp parallel for num_threads(nthreads) schedule(static) firstprivate(molb)
327 for (int i = 0; i < md->nr; i++)
343 const t_atom& atom = mtopGetAtomParameters(mtop, ag, &molb);
347 md->cFREEZE[i] = getGroupType(groups, SimulationAtomGroupType::Freeze, ag);
349 if (EI_ENERGY_MINIMIZATION(inputrec.eI))
351 /* Displacement is proportional to F, masses used for constraints */
355 else if (inputrec.eI == IntegrationAlgorithm::BD)
357 /* With BD the physical masses are irrelevant.
358 * To keep the code simple we use most of the normal MD code path
359 * for BD. Thus for constraining the masses should be proportional
360 * to the friction coefficient. We set the absolute value such that
361 * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
362 * Then if we set the (meaningless) velocity to v=dx/dt, we get the
363 * correct kinetic energy and temperature using the usual code path.
364 * Thus with BD v*dt will give the displacement and the reported
365 * temperature can signal bad integration (too large time step).
367 if (inputrec.bd_fric > 0)
369 mA = 0.5 * inputrec.bd_fric * inputrec.delta_t;
370 mB = 0.5 * inputrec.bd_fric * inputrec.delta_t;
374 /* The friction coefficient is mass/tau_t */
375 fac = inputrec.delta_t
376 / opts->tau_t[md->cTC ? groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ag] : 0];
377 mA = 0.5 * atom.m * fac;
378 mB = 0.5 * atom.mB * fac;
386 if (md->nMassPerturbed)
396 md->invMassPerDim[i][XX] = 0;
397 md->invMassPerDim[i][YY] = 0;
398 md->invMassPerDim[i][ZZ] = 0;
400 else if (md->cFREEZE)
403 GMX_ASSERT(opts->nFreeze != nullptr, "Must have freeze groups to initialize masses");
404 if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
406 /* Set the mass of completely frozen particles to ALMOST_ZERO
407 * iso 0 to avoid div by zero in lincs or shake.
409 md->invmass[i] = ALMOST_ZERO;
413 /* Note: Partially frozen particles use the normal invmass.
414 * If such particles are constrained, the frozen dimensions
415 * should not be updated with the constrained coordinates.
417 md->invmass[i] = 1.0 / mA;
419 for (int d = 0; d < DIM; d++)
421 md->invMassPerDim[i][d] = (opts->nFreeze[g][d] ? 0 : 1.0 / mA);
426 md->invmass[i] = 1.0 / mA;
427 for (int d = 0; d < DIM; d++)
429 md->invMassPerDim[i][d] = 1.0 / mA;
433 md->chargeA[i] = atom.q;
434 md->typeA[i] = atom.type;
437 c6 = mtop.ffparams.iparams[atom.type * (mtop.ffparams.atnr + 1)].lj.c6;
438 c12 = mtop.ffparams.iparams[atom.type * (mtop.ffparams.atnr + 1)].lj.c12;
439 md->sqrt_c6A[i] = std::sqrt(c6);
440 if (c6 == 0.0 || c12 == 0)
446 md->sigmaA[i] = gmx::sixthroot(c12 / c6);
448 md->sigma3A[i] = 1 / (md->sigmaA[i] * md->sigmaA[i] * md->sigmaA[i]);
452 md->bPerturbed[i] = PERTURBED(atom);
453 md->chargeB[i] = atom.qB;
454 md->typeB[i] = atom.typeB;
457 c6 = mtop.ffparams.iparams[atom.typeB * (mtop.ffparams.atnr + 1)].lj.c6;
458 c12 = mtop.ffparams.iparams[atom.typeB * (mtop.ffparams.atnr + 1)].lj.c12;
459 md->sqrt_c6B[i] = std::sqrt(c6);
460 if (c6 == 0.0 || c12 == 0)
466 md->sigmaB[i] = gmx::sixthroot(c12 / c6);
468 md->sigma3B[i] = 1 / (md->sigmaB[i] * md->sigmaB[i] * md->sigmaB[i]);
471 md->ptype[i] = atom.ptype;
474 md->cTC[i] = groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ag];
476 md->cENER[i] = getGroupType(groups, SimulationAtomGroupType::EnergyOutput, ag);
479 md->cVCM[i] = groups.groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][ag];
483 md->cORF[i] = getGroupType(groups, SimulationAtomGroupType::OrientationRestraintsFit, ag);
488 md->cU1[i] = groups.groupNumbers[SimulationAtomGroupType::User1][ag];
492 md->cU2[i] = groups.groupNumbers[SimulationAtomGroupType::User2][ag];
495 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
500 /* Pad invmass with 0 so a SIMD MD update does not change v and x */
501 for (int i = md->nr; i < md->nr + GMX_REAL_MAX_SIMD_WIDTH; i++)
508 /* We set mass, invmass, invMassPerDim and tmass for lambda=0.
509 * For free-energy runs, these should be updated using update_mdatoms().
511 md->tmass = md->tmassA;
515 void update_mdatoms(t_mdatoms* md, real lambda)
517 if (md->nMassPerturbed && lambda != md->lambda)
519 real L1 = 1 - lambda;
521 /* Update masses of perturbed atoms for the change in lambda */
522 int gmx_unused nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Default);
523 #pragma omp parallel for num_threads(nthreads) schedule(static)
524 for (int i = 0; i < md->nr; i++)
526 if (md->bPerturbed[i])
528 md->massT[i] = L1 * md->massA[i] + lambda * md->massB[i];
529 /* Atoms with invmass 0 or ALMOST_ZERO are massless or frozen
530 * and their invmass does not depend on lambda.
532 if (md->invmass[i] > 1.1 * ALMOST_ZERO)
534 md->invmass[i] = 1.0 / md->massT[i];
535 for (int d = 0; d < DIM; d++)
537 if (md->invMassPerDim[i][d] > 1.1 * ALMOST_ZERO)
539 md->invMassPerDim[i][d] = md->invmass[i];
546 /* Update the system mass for the change in lambda */
547 md->tmass = L1 * md->tmassA + lambda * md->tmassB;