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,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/compat/make_unique.h"
46 #include "gromacs/gpu_utils/hostallocator.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/mdlib/gmx_omp_nthreads.h"
49 #include "gromacs/mdlib/qmmm.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/topology/mtop_lookup.h"
53 #include "gromacs/topology/mtop_util.h"
54 #include "gromacs/topology/topology.h"
55 #include "gromacs/utility/exceptions.h"
56 #include "gromacs/utility/smalloc.h"
58 #define ALMOST_ZERO 1e-30
64 : mdatoms_(nullptr), chargeA_()
68 void MDAtoms::resize(int newSize
)
70 chargeA_
.resize(newSize
);
71 mdatoms_
->chargeA
= chargeA_
.data();
74 void MDAtoms::reserve(int newCapacity
)
76 chargeA_
.reserve(newCapacity
);
77 mdatoms_
->chargeA
= chargeA_
.data();
80 std::unique_ptr
<MDAtoms
>
81 makeMDAtoms(FILE *fp
, const gmx_mtop_t
&mtop
, const t_inputrec
&ir
,
84 auto mdAtoms
= compat::make_unique
<MDAtoms
>();
85 // GPU transfers want to use the pinning mode.
86 changePinningPolicy(&mdAtoms
->chargeA_
, useGpuForPme
? PinningPolicy::CanBePinned
: PinningPolicy::CannotBePinned
);
89 mdAtoms
->mdatoms_
.reset(md
);
91 md
->nenergrp
= mtop
.groups
.grps
[egcENER
].nr
;
92 md
->bVCMgrps
= (mtop
.groups
.grps
[egcVCM
].nr
> 1);
94 /* Determine the total system mass and perturbed atom counts */
95 double totalMassA
= 0.0;
96 double totalMassB
= 0.0;
98 gmx_mtop_atomloop_block_t aloop
= gmx_mtop_atomloop_block_init(&mtop
);
101 while (gmx_mtop_atomloop_block_next(aloop
, &atom
, &nmol
))
103 totalMassA
+= nmol
*atom
->m
;
104 totalMassB
+= nmol
*atom
->mB
;
106 if (ir
.efep
!= efepNO
&& PERTURBED(*atom
))
109 if (atom
->mB
!= atom
->m
)
111 md
->nMassPerturbed
+= nmol
;
113 if (atom
->qB
!= atom
->q
)
115 md
->nChargePerturbed
+= nmol
;
117 if (atom
->typeB
!= atom
->type
)
119 md
->nTypePerturbed
+= nmol
;
124 md
->tmassA
= totalMassA
;
125 md
->tmassB
= totalMassB
;
127 if (ir
.efep
!= efepNO
&& fp
)
130 "There are %d atoms and %d charges for free energy perturbation\n",
131 md
->nPerturbed
, md
->nChargePerturbed
);
134 md
->havePartiallyFrozenAtoms
= FALSE
;
135 for (int g
= 0; g
< ir
.opts
.ngfrz
; g
++)
137 for (int d
= YY
; d
< DIM
; d
++)
139 if (ir
.opts
.nFreeze
[d
] != ir
.opts
.nFreeze
[XX
])
141 md
->havePartiallyFrozenAtoms
= TRUE
;
146 md
->bOrires
= gmx_mtop_ftype_count(&mtop
, F_ORIRES
);
153 void atoms2md(const gmx_mtop_t
*mtop
, const t_inputrec
*ir
,
154 int nindex
, const int *index
,
156 gmx::MDAtoms
*mdAtoms
)
159 const t_grpopts
*opts
;
160 const gmx_groups_t
*groups
;
161 int nthreads gmx_unused
;
163 bLJPME
= EVDW_PME(ir
->vdwtype
);
167 groups
= &mtop
->groups
;
169 auto md
= mdAtoms
->mdatoms();
170 /* nindex>=0 indicates DD where we use an index */
177 md
->nr
= mtop
->natoms
;
180 if (md
->nr
> md
->nalloc
)
182 md
->nalloc
= over_alloc_dd(md
->nr
);
184 if (md
->nMassPerturbed
)
186 srenew(md
->massA
, md
->nalloc
);
187 srenew(md
->massB
, md
->nalloc
);
189 srenew(md
->massT
, md
->nalloc
);
190 /* The SIMD version of the integrator needs this aligned and padded.
191 * The padding needs to be with zeros, which we set later below.
193 gmx::AlignedAllocationPolicy::free(md
->invmass
);
194 md
->invmass
= new(gmx::AlignedAllocationPolicy::malloc((md
->nalloc
+ GMX_REAL_MAX_SIMD_WIDTH
)*sizeof(*md
->invmass
)))real
;
195 srenew(md
->invMassPerDim
, md
->nalloc
);
196 // TODO eventually we will have vectors and just resize
197 // everything, but for now the semantics of md->nalloc being
198 // the capacity are preserved by keeping vectors within
199 // mdAtoms having the same properties as the other arrays.
200 mdAtoms
->reserve(md
->nalloc
);
201 mdAtoms
->resize(md
->nr
);
202 srenew(md
->typeA
, md
->nalloc
);
205 srenew(md
->chargeB
, md
->nalloc
);
206 srenew(md
->typeB
, md
->nalloc
);
210 srenew(md
->sqrt_c6A
, md
->nalloc
);
211 srenew(md
->sigmaA
, md
->nalloc
);
212 srenew(md
->sigma3A
, md
->nalloc
);
215 srenew(md
->sqrt_c6B
, md
->nalloc
);
216 srenew(md
->sigmaB
, md
->nalloc
);
217 srenew(md
->sigma3B
, md
->nalloc
);
220 srenew(md
->ptype
, md
->nalloc
);
223 srenew(md
->cTC
, md
->nalloc
);
224 /* We always copy cTC with domain decomposition */
226 srenew(md
->cENER
, md
->nalloc
);
229 srenew(md
->cACC
, md
->nalloc
);
233 opts
->nFreeze
[0][XX
] || opts
->nFreeze
[0][YY
] || opts
->nFreeze
[0][ZZ
]))
235 srenew(md
->cFREEZE
, md
->nalloc
);
239 srenew(md
->cVCM
, md
->nalloc
);
243 srenew(md
->cORF
, md
->nalloc
);
247 srenew(md
->bPerturbed
, md
->nalloc
);
250 /* Note that these user t_mdatoms array pointers are NULL
251 * when there is only one group present.
252 * Therefore, when adding code, the user should use something like:
253 * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
255 if (mtop
->groups
.grpnr
[egcUser1
] != nullptr)
257 srenew(md
->cU1
, md
->nalloc
);
259 if (mtop
->groups
.grpnr
[egcUser2
] != nullptr)
261 srenew(md
->cU2
, md
->nalloc
);
266 srenew(md
->bQM
, md
->nalloc
);
272 // cppcheck-suppress unreadVariable
273 nthreads
= gmx_omp_nthreads_get(emntDefault
);
274 #pragma omp parallel for num_threads(nthreads) schedule(static) firstprivate(molb)
275 for (int i
= 0; i
< md
->nr
; i
++)
283 if (index
== nullptr)
291 const t_atom
&atom
= mtopGetAtomParameters(mtop
, ag
, &molb
);
295 md
->cFREEZE
[i
] = ggrpnr(groups
, egcFREEZE
, ag
);
297 if (EI_ENERGY_MINIMIZATION(ir
->eI
))
299 /* Displacement is proportional to F, masses used for constraints */
303 else if (ir
->eI
== eiBD
)
305 /* With BD the physical masses are irrelevant.
306 * To keep the code simple we use most of the normal MD code path
307 * for BD. Thus for constraining the masses should be proportional
308 * to the friction coefficient. We set the absolute value such that
309 * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
310 * Then if we set the (meaningless) velocity to v=dx/dt, we get the
311 * correct kinetic energy and temperature using the usual code path.
312 * Thus with BD v*dt will give the displacement and the reported
313 * temperature can signal bad integration (too large time step).
317 mA
= 0.5*ir
->bd_fric
*ir
->delta_t
;
318 mB
= 0.5*ir
->bd_fric
*ir
->delta_t
;
322 /* The friction coefficient is mass/tau_t */
323 fac
= ir
->delta_t
/opts
->tau_t
[md
->cTC
? groups
->grpnr
[egcTC
][ag
] : 0];
325 mB
= 0.5*atom
.mB
*fac
;
333 if (md
->nMassPerturbed
)
343 md
->invMassPerDim
[i
][XX
] = 0;
344 md
->invMassPerDim
[i
][YY
] = 0;
345 md
->invMassPerDim
[i
][ZZ
] = 0;
347 else if (md
->cFREEZE
)
350 if (opts
->nFreeze
[g
][XX
] && opts
->nFreeze
[g
][YY
] && opts
->nFreeze
[g
][ZZ
])
352 /* Set the mass of completely frozen particles to ALMOST_ZERO
353 * iso 0 to avoid div by zero in lincs or shake.
355 md
->invmass
[i
] = ALMOST_ZERO
;
359 /* Note: Partially frozen particles use the normal invmass.
360 * If such particles are constrained, the frozen dimensions
361 * should not be updated with the constrained coordinates.
363 md
->invmass
[i
] = 1.0/mA
;
365 for (int d
= 0; d
< DIM
; d
++)
367 md
->invMassPerDim
[i
][d
] = (opts
->nFreeze
[g
][d
] ? 0 : 1.0/mA
);
372 md
->invmass
[i
] = 1.0/mA
;
373 for (int d
= 0; d
< DIM
; d
++)
375 md
->invMassPerDim
[i
][d
] = 1.0/mA
;
379 md
->chargeA
[i
] = atom
.q
;
380 md
->typeA
[i
] = atom
.type
;
383 c6
= mtop
->ffparams
.iparams
[atom
.type
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
384 c12
= mtop
->ffparams
.iparams
[atom
.type
*(mtop
->ffparams
.atnr
+1)].lj
.c12
;
385 md
->sqrt_c6A
[i
] = sqrt(c6
);
386 if (c6
== 0.0 || c12
== 0)
392 md
->sigmaA
[i
] = gmx::sixthroot(c12
/c6
);
394 md
->sigma3A
[i
] = 1/(md
->sigmaA
[i
]*md
->sigmaA
[i
]*md
->sigmaA
[i
]);
398 md
->bPerturbed
[i
] = PERTURBED(atom
);
399 md
->chargeB
[i
] = atom
.qB
;
400 md
->typeB
[i
] = atom
.typeB
;
403 c6
= mtop
->ffparams
.iparams
[atom
.typeB
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
404 c12
= mtop
->ffparams
.iparams
[atom
.typeB
*(mtop
->ffparams
.atnr
+1)].lj
.c12
;
405 md
->sqrt_c6B
[i
] = sqrt(c6
);
406 if (c6
== 0.0 || c12
== 0)
412 md
->sigmaB
[i
] = gmx::sixthroot(c12
/c6
);
414 md
->sigma3B
[i
] = 1/(md
->sigmaB
[i
]*md
->sigmaB
[i
]*md
->sigmaB
[i
]);
417 md
->ptype
[i
] = atom
.ptype
;
420 md
->cTC
[i
] = groups
->grpnr
[egcTC
][ag
];
422 md
->cENER
[i
] = ggrpnr(groups
, egcENER
, ag
);
425 md
->cACC
[i
] = groups
->grpnr
[egcACC
][ag
];
429 md
->cVCM
[i
] = groups
->grpnr
[egcVCM
][ag
];
433 md
->cORF
[i
] = ggrpnr(groups
, egcORFIT
, ag
);
438 md
->cU1
[i
] = groups
->grpnr
[egcUser1
][ag
];
442 md
->cU2
[i
] = groups
->grpnr
[egcUser2
][ag
];
447 if (groups
->grpnr
[egcQMMM
] == nullptr ||
448 groups
->grpnr
[egcQMMM
][ag
] < groups
->grps
[egcQMMM
].nr
-1)
458 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
463 /* Pad invmass with 0 so a SIMD MD update does not change v and x */
464 for (int i
= md
->nr
; i
< md
->nr
+ GMX_REAL_MAX_SIMD_WIDTH
; i
++)
471 /* We set mass, invmass, invMassPerDim and tmass for lambda=0.
472 * For free-energy runs, these should be updated using update_mdatoms().
474 md
->tmass
= md
->tmassA
;
478 void update_mdatoms(t_mdatoms
*md
, real lambda
)
480 if (md
->nMassPerturbed
&& lambda
!= md
->lambda
)
482 real L1
= 1 - lambda
;
484 /* Update masses of perturbed atoms for the change in lambda */
485 // cppcheck-suppress unreadVariable
486 int gmx_unused nthreads
= gmx_omp_nthreads_get(emntDefault
);
487 #pragma omp parallel for num_threads(nthreads) schedule(static)
488 for (int i
= 0; i
< md
->nr
; i
++)
490 if (md
->bPerturbed
[i
])
492 md
->massT
[i
] = L1
*md
->massA
[i
] + lambda
*md
->massB
[i
];
493 /* Atoms with invmass 0 or ALMOST_ZERO are massless or frozen
494 * and their invmass does not depend on lambda.
496 if (md
->invmass
[i
] > 1.1*ALMOST_ZERO
)
498 md
->invmass
[i
] = 1.0/md
->massT
[i
];
499 for (int d
= 0; d
< DIM
; d
++)
501 if (md
->invMassPerDim
[i
][d
] > 1.1*ALMOST_ZERO
)
503 md
->invMassPerDim
[i
][d
] = md
->invmass
[i
];
510 /* Update the system mass for the change in lambda */
511 md
->tmass
= L1
*md
->tmassA
+ lambda
*md
->tmassB
;