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,2016, 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.
41 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
42 #include "gromacs/legacyheaders/mdatoms.h"
43 #include "gromacs/legacyheaders/qmmm.h"
44 #include "gromacs/legacyheaders/typedefs.h"
45 #include "gromacs/topology/mtop_util.h"
46 #include "gromacs/utility/smalloc.h"
48 #define ALMOST_ZERO 1e-30
50 t_mdatoms
*init_mdatoms(FILE *fp
, gmx_mtop_t
*mtop
, gmx_bool bFreeEnergy
)
56 gmx_mtop_atomloop_all_t aloop
;
60 md
->nenergrp
= mtop
->groups
.grps
[egcENER
].nr
;
65 aloop
= gmx_mtop_atomloop_all_init(mtop
);
66 while (gmx_mtop_atomloop_all_next(aloop
, &a
, &atom
))
68 if (ggrpnr(&mtop
->groups
, egcVCM
, a
) > 0)
73 if (bFreeEnergy
&& PERTURBED(*atom
))
76 if (atom
->mB
!= atom
->m
)
80 if (atom
->qB
!= atom
->q
)
82 md
->nChargePerturbed
++;
84 if (atom
->typeB
!= atom
->type
)
97 if (bFreeEnergy
&& fp
)
100 "There are %d atoms and %d charges for free energy perturbation\n",
101 md
->nPerturbed
, md
->nChargePerturbed
);
104 md
->bOrires
= gmx_mtop_ftype_count(mtop
, F_ORIRES
);
109 void atoms2md(gmx_mtop_t
*mtop
, t_inputrec
*ir
,
110 int nindex
, int *index
,
115 gmx_mtop_atomlookup_t alook
;
118 gmx_groups_t
*groups
;
119 int nthreads gmx_unused
;
120 const real oneOverSix
= 1.0 / 6.0;
122 bLJPME
= EVDW_PME(ir
->vdwtype
);
126 groups
= &mtop
->groups
;
128 /* Index==NULL indicates no DD (unless we have a DD node with no
129 * atoms), so also check for homenr. This should be
130 * signaled properly with an extra parameter or nindex==-1.
132 if (index
== NULL
&& (homenr
> 0))
134 md
->nr
= mtop
->natoms
;
141 if (md
->nr
> md
->nalloc
)
143 md
->nalloc
= over_alloc_dd(md
->nr
);
145 if (md
->nMassPerturbed
)
147 srenew(md
->massA
, md
->nalloc
);
148 srenew(md
->massB
, md
->nalloc
);
150 srenew(md
->massT
, md
->nalloc
);
151 srenew(md
->invmass
, md
->nalloc
);
152 srenew(md
->chargeA
, md
->nalloc
);
153 srenew(md
->typeA
, md
->nalloc
);
156 srenew(md
->chargeB
, md
->nalloc
);
157 srenew(md
->typeB
, md
->nalloc
);
161 srenew(md
->sqrt_c6A
, md
->nalloc
);
162 srenew(md
->sigmaA
, md
->nalloc
);
163 srenew(md
->sigma3A
, md
->nalloc
);
166 srenew(md
->sqrt_c6B
, md
->nalloc
);
167 srenew(md
->sigmaB
, md
->nalloc
);
168 srenew(md
->sigma3B
, md
->nalloc
);
171 srenew(md
->ptype
, md
->nalloc
);
174 srenew(md
->cTC
, md
->nalloc
);
175 /* We always copy cTC with domain decomposition */
177 srenew(md
->cENER
, md
->nalloc
);
180 srenew(md
->cACC
, md
->nalloc
);
184 opts
->nFreeze
[0][XX
] || opts
->nFreeze
[0][YY
] || opts
->nFreeze
[0][ZZ
]))
186 srenew(md
->cFREEZE
, md
->nalloc
);
190 srenew(md
->cVCM
, md
->nalloc
);
194 srenew(md
->cORF
, md
->nalloc
);
198 srenew(md
->bPerturbed
, md
->nalloc
);
201 /* Note that these user t_mdatoms array pointers are NULL
202 * when there is only one group present.
203 * Therefore, when adding code, the user should use something like:
204 * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
206 if (mtop
->groups
.grpnr
[egcUser1
] != NULL
)
208 srenew(md
->cU1
, md
->nalloc
);
210 if (mtop
->groups
.grpnr
[egcUser2
] != NULL
)
212 srenew(md
->cU2
, md
->nalloc
);
217 srenew(md
->bQM
, md
->nalloc
);
221 srenew(md
->wf
, md
->nalloc
);
222 srenew(md
->tf_table_index
, md
->nalloc
);
226 alook
= gmx_mtop_atomlookup_init(mtop
);
228 // cppcheck-suppress unreadVariable
229 nthreads
= gmx_omp_nthreads_get(emntDefault
);
230 #pragma omp parallel for num_threads(nthreads) schedule(static)
231 for (i
= 0; i
< md
->nr
; i
++)
246 gmx_mtop_atomnr_to_atom(alook
, ag
, &atom
);
250 md
->cFREEZE
[i
] = ggrpnr(groups
, egcFREEZE
, ag
);
252 if (EI_ENERGY_MINIMIZATION(ir
->eI
))
254 /* Displacement is proportional to F, masses used for constraints */
258 else if (ir
->eI
== eiBD
)
260 /* With BD the physical masses are irrelevant.
261 * To keep the code simple we use most of the normal MD code path
262 * for BD. Thus for constraining the masses should be proportional
263 * to the friction coefficient. We set the absolute value such that
264 * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
265 * Then if we set the (meaningless) velocity to v=dx/dt, we get the
266 * correct kinetic energy and temperature using the usual code path.
267 * Thus with BD v*dt will give the displacement and the reported
268 * temperature can signal bad integration (too large time step).
272 mA
= 0.5*ir
->bd_fric
*ir
->delta_t
;
273 mB
= 0.5*ir
->bd_fric
*ir
->delta_t
;
277 /* The friction coefficient is mass/tau_t */
278 fac
= ir
->delta_t
/opts
->tau_t
[md
->cTC
? groups
->grpnr
[egcTC
][ag
] : 0];
279 mA
= 0.5*atom
->m
*fac
;
280 mB
= 0.5*atom
->mB
*fac
;
288 if (md
->nMassPerturbed
)
298 else if (md
->cFREEZE
)
301 if (opts
->nFreeze
[g
][XX
] && opts
->nFreeze
[g
][YY
] && opts
->nFreeze
[g
][ZZ
])
303 /* Set the mass of completely frozen particles to ALMOST_ZERO
304 * iso 0 to avoid div by zero in lincs or shake.
306 md
->invmass
[i
] = ALMOST_ZERO
;
310 /* Note: Partially frozen particles use the normal invmass.
311 * If such particles are constrained, the frozen dimensions
312 * should not be updated with the constrained coordinates.
314 md
->invmass
[i
] = 1.0/mA
;
319 md
->invmass
[i
] = 1.0/mA
;
321 md
->chargeA
[i
] = atom
->q
;
322 md
->typeA
[i
] = atom
->type
;
325 c6
= mtop
->ffparams
.iparams
[atom
->type
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
326 c12
= mtop
->ffparams
.iparams
[atom
->type
*(mtop
->ffparams
.atnr
+1)].lj
.c12
;
327 md
->sqrt_c6A
[i
] = sqrt(c6
);
328 if (c6
== 0.0 || c12
== 0)
334 md
->sigmaA
[i
] = pow(c12
/c6
, oneOverSix
);
336 md
->sigma3A
[i
] = 1/(md
->sigmaA
[i
]*md
->sigmaA
[i
]*md
->sigmaA
[i
]);
340 md
->bPerturbed
[i
] = PERTURBED(*atom
);
341 md
->chargeB
[i
] = atom
->qB
;
342 md
->typeB
[i
] = atom
->typeB
;
345 c6
= mtop
->ffparams
.iparams
[atom
->typeB
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
346 c12
= mtop
->ffparams
.iparams
[atom
->typeB
*(mtop
->ffparams
.atnr
+1)].lj
.c12
;
347 md
->sqrt_c6B
[i
] = sqrt(c6
);
348 if (c6
== 0.0 || c12
== 0)
354 md
->sigmaB
[i
] = pow(c12
/c6
, oneOverSix
);
356 md
->sigma3B
[i
] = 1/(md
->sigmaB
[i
]*md
->sigmaB
[i
]*md
->sigmaB
[i
]);
359 md
->ptype
[i
] = atom
->ptype
;
362 md
->cTC
[i
] = groups
->grpnr
[egcTC
][ag
];
365 (groups
->grpnr
[egcENER
] ? groups
->grpnr
[egcENER
][ag
] : 0);
368 md
->cACC
[i
] = groups
->grpnr
[egcACC
][ag
];
372 md
->cVCM
[i
] = groups
->grpnr
[egcVCM
][ag
];
376 md
->cORF
[i
] = groups
->grpnr
[egcORFIT
][ag
];
381 md
->cU1
[i
] = groups
->grpnr
[egcUser1
][ag
];
385 md
->cU2
[i
] = groups
->grpnr
[egcUser2
][ag
];
390 if (groups
->grpnr
[egcQMMM
] == 0 ||
391 groups
->grpnr
[egcQMMM
][ag
] < groups
->grps
[egcQMMM
].nr
-1)
400 /* Initialize AdResS weighting functions to adressw */
404 /* if no tf table groups specified, use default table */
405 md
->tf_table_index
[i
] = DEFAULT_TF_TABLE
;
406 if (ir
->adress
->n_tf_grps
> 0)
408 /* if tf table groups specified, tf is only applied to thoose energy groups*/
409 md
->tf_table_index
[i
] = NO_TF_TABLE
;
410 /* check wether atom is in one of the relevant energy groups and assign a table index */
411 for (g
= 0; g
< ir
->adress
->n_tf_grps
; g
++)
413 if (md
->cENER
[i
] == ir
->adress
->tf_table_index
[g
])
415 md
->tf_table_index
[i
] = g
;
422 gmx_mtop_atomlookup_destroy(alook
);
428 void update_mdatoms(t_mdatoms
*md
, real lambda
)
431 real L1
= 1.0-lambda
;
435 if (md
->nMassPerturbed
)
437 for (al
= 0; (al
< end
); al
++)
439 if (md
->bPerturbed
[al
])
441 md
->massT
[al
] = L1
*md
->massA
[al
]+ lambda
*md
->massB
[al
];
442 if (md
->invmass
[al
] > 1.1*ALMOST_ZERO
)
444 md
->invmass
[al
] = 1.0/md
->massT
[al
];
448 md
->tmass
= L1
*md
->tmassA
+ lambda
*md
->tmassB
;
452 md
->tmass
= md
->tmassA
;