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) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, 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.
38 /* This file is completely threadsafe - keep it that way! */
45 #include "gromacs/gmxlib/network.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/mdlib/gmx_omp_nthreads.h"
48 #include "gromacs/mdlib/rbin.h"
49 #include "gromacs/mdlib/update.h"
50 #include "gromacs/mdtypes/group.h"
51 #include "gromacs/mdtypes/inputrec.h"
52 #include "gromacs/mdtypes/mdatom.h"
53 #include "gromacs/topology/mtop_util.h"
54 #include "gromacs/topology/topology.h"
55 #include "gromacs/utility/exceptions.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 static void init_grpstat(const gmx_mtop_t
* mtop
, int ngacc
, t_grp_acc gstat
[])
64 const SimulationGroups
& groups
= mtop
->groups
;
65 for (const AtomProxy atomP
: AtomRange(*mtop
))
67 const t_atom
& local
= atomP
.atom();
68 int i
= atomP
.globalAtomNumber();
69 int grp
= getGroupType(groups
, SimulationAtomGroupType::Acceleration
, i
);
70 if ((grp
< 0) && (grp
>= ngacc
))
72 gmx_incons("Input for acceleration groups wrong");
75 /* This will not work for integrator BD */
76 gstat
[grp
].mA
+= local
.m
;
77 gstat
[grp
].mB
+= local
.mB
;
82 void init_ekindata(FILE gmx_unused
* log
,
83 const gmx_mtop_t
* mtop
,
84 const t_grpopts
* opts
,
85 gmx_ekindata_t
* ekind
,
90 /* bNEMD tells if we should remove remove the COM velocity
91 * from the velocities during velocity scaling in T-coupling.
92 * Turn this on when we have multiple acceleration groups
93 * or one accelerated group.
95 ekind
->bNEMD
= (opts
->ngacc
> 1 || norm2(opts
->acc
[0]) > 0);
97 ekind
->ngtc
= opts
->ngtc
;
98 ekind
->tcstat
.resize(opts
->ngtc
);
99 /* Set Berendsen tcoupl lambda's to 1,
100 * so runs without Berendsen coupling are not affected.
102 for (i
= 0; i
< opts
->ngtc
; i
++)
104 ekind
->tcstat
[i
].lambda
= 1.0;
105 ekind
->tcstat
[i
].vscale_nhc
= 1.0;
106 ekind
->tcstat
[i
].ekinscaleh_nhc
= 1.0;
107 ekind
->tcstat
[i
].ekinscalef_nhc
= 1.0;
110 int nthread
= gmx_omp_nthreads_get(emntUpdate
);
111 ekind
->nthreads
= nthread
;
112 snew(ekind
->ekin_work_alloc
, nthread
);
113 snew(ekind
->ekin_work
, nthread
);
114 snew(ekind
->dekindl_work
, nthread
);
115 #pragma omp parallel for num_threads(nthread) schedule(static)
116 for (int thread
= 0; thread
< nthread
; thread
++)
120 #define EKIN_WORK_BUFFER_SIZE 2
121 /* Allocate 2 extra elements on both sides, so in single
123 * EKIN_WORK_BUFFER_SIZE*DIM*DIM*sizeof(real) = 72/144 bytes
124 * buffer on both sides to avoid cache pollution.
126 snew(ekind
->ekin_work_alloc
[thread
], ekind
->ngtc
+ 2 * EKIN_WORK_BUFFER_SIZE
);
127 ekind
->ekin_work
[thread
] = ekind
->ekin_work_alloc
[thread
] + EKIN_WORK_BUFFER_SIZE
;
128 /* Nasty hack so we can have the per-thread accumulation
129 * variable for dekindl in the same thread-local cache lines
130 * as the per-thread accumulation tensors for ekin[fh],
131 * because they are accumulated in the same loop. */
132 ekind
->dekindl_work
[thread
] = &(ekind
->ekin_work
[thread
][ekind
->ngtc
][0][0]);
133 #undef EKIN_WORK_BUFFER_SIZE
135 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
138 ekind
->ngacc
= opts
->ngacc
;
139 ekind
->grpstat
.resize(opts
->ngacc
);
140 init_grpstat(mtop
, opts
->ngacc
, ekind
->grpstat
.data());
142 ekind
->cosacc
.cos_accel
= cos_accel
;
145 void accumulate_u(const t_commrec
* cr
, const t_grpopts
* opts
, gmx_ekindata_t
* ekind
)
147 /* This routine will only be called when it's necessary */
153 for (g
= 0; (g
< opts
->ngacc
); g
++)
155 add_binr(rb
, DIM
, ekind
->grpstat
[g
].u
);
159 for (g
= 0; (g
< opts
->ngacc
); g
++)
161 extract_binr(rb
, DIM
* g
, DIM
, ekind
->grpstat
[g
].u
);
166 void update_ekindata(int start
,
168 gmx_ekindata_t
* ekind
,
169 const t_grpopts
* opts
,
177 /* calculate mean velocities at whole timestep */
178 for (g
= 0; (g
< opts
->ngtc
); g
++)
180 ekind
->tcstat
[g
].T
= 0;
185 for (g
= 0; (g
< opts
->ngacc
); g
++)
187 clear_rvec(ekind
->grpstat
[g
].u
);
191 for (n
= start
; (n
< start
+ homenr
); n
++)
197 for (d
= 0; (d
< DIM
); d
++)
199 mv
= md
->massT
[n
] * v
[n
][d
];
200 ekind
->grpstat
[g
].u
[d
] += mv
;
204 for (g
= 0; (g
< opts
->ngacc
); g
++)
206 for (d
= 0; (d
< DIM
); d
++)
208 ekind
->grpstat
[g
].u
[d
] /=
209 (1 - lambda
) * ekind
->grpstat
[g
].mA
+ lambda
* ekind
->grpstat
[g
].mB
;
215 real
sum_ekin(const t_grpopts
* opts
, gmx_ekindata_t
* ekind
, real
* dekindlambda
, gmx_bool bEkinAveVel
, gmx_bool bScaleEkin
)
219 t_grp_tcstat
* tcstat
;
228 clear_mat(ekind
->ekin
);
230 for (i
= 0; (i
< ngtc
); i
++)
234 tcstat
= &ekind
->tcstat
[i
];
235 /* Sometimes a group does not have degrees of freedom, e.g.
236 * when it consists of shells and virtual sites, then we just
237 * set the temperatue to 0 and also neglect the kinetic
238 * energy, which should be zero anyway.
247 /* in this case, kinetic energy is from the current velocities already */
248 msmul(tcstat
->ekinf
, tcstat
->ekinscalef_nhc
, tcstat
->ekinf
);
253 /* Calculate the full step Ekin as the average of the half steps */
254 for (j
= 0; (j
< DIM
); j
++)
256 for (m
= 0; (m
< DIM
); m
++)
258 tcstat
->ekinf
[j
][m
] = 0.5
259 * (tcstat
->ekinh
[j
][m
] * tcstat
->ekinscaleh_nhc
260 + tcstat
->ekinh_old
[j
][m
]);
264 m_add(tcstat
->ekinf
, ekind
->ekin
, ekind
->ekin
);
266 tcstat
->Th
= calc_temp(trace(tcstat
->ekinh
), nd
);
267 tcstat
->T
= calc_temp(trace(tcstat
->ekinf
), nd
);
269 /* after the scaling factors have been multiplied in, we can remove them */
272 tcstat
->ekinscalef_nhc
= 1.0;
276 tcstat
->ekinscaleh_nhc
= 1.0;
295 *dekindlambda
= ekind
->dekindl
;
299 *dekindlambda
= 0.5 * (ekind
->dekindl
+ ekind
->dekindl_old
);