* 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 <assert.h>
-#include "typedefs.h"
-#include "types/commrec.h"
-#include "gromacs/utility/smalloc.h"
-#include "update.h"
-#include "vec.h"
-#include "macros.h"
-#include "physics.h"
-#include "names.h"
-#include "gmx_fatal.h"
-#include "txtdump.h"
-#include "nrnb.h"
+#include <algorithm>
+
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/update.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
#include "gromacs/random/random.h"
-#include "update.h"
-#include "mdrun.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#define NTROTTERPARTS 3
{
/* general routine for both barostat and thermostat nose hoover chains */
- int i, j, mi, mj, jmax;
+ int i, j, mi, mj;
double Ekin, Efac, reft, kT, nd;
double dt;
t_grp_tcstat *tcstat;
{
iQinv = &(MassQ->QPinv[i*nh]);
nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
- reft = max(0.0, opts->ref_t[0]);
+ reft = std::max<real>(0, opts->ref_t[0]);
Ekin = sqr(*veta)/MassQ->Winv;
}
else
iQinv = &(MassQ->Qinv[i*nh]);
tcstat = &ekind->tcstat[i];
nd = opts->nrdf[i];
- reft = max(0.0, opts->ref_t[i]);
+ reft = std::max<real>(0, opts->ref_t[i]);
if (bEkinAveVel)
{
Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
real pscal;
double alpha;
- int i, j, d, n, nwall;
- real T, GW, vol;
- tensor Winvm, ekinmod, localpres;
+ int nwall;
+ real GW, vol;
+ tensor ekinmod, localpres;
/* The heat bath is coupled to a separate barostat, the last temperature group. In the
2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
* The pressure and compressibility always occur as a product,
* therefore the pressure unit drops out.
*/
- maxl = max(box[XX][XX], box[YY][YY]);
- maxl = max(maxl, box[ZZ][ZZ]);
+ maxl = std::max(box[XX][XX], box[YY][YY]);
+ maxl = std::max(maxl, box[ZZ][ZZ]);
for (d = 0; d < DIM; d++)
{
for (n = 0; n < DIM; n++)
{
int d, n;
real scalar_pressure, xy_pressure, p_corr_z;
- char *ptr, buf[STRLEN];
+ char buf[STRLEN];
/*
* Calculate the scaling matrix mu
t_nrnb *nrnb)
{
ivec *nFreeze = ir->opts.nFreeze;
- int n, d, g = 0;
+ int n, d;
+ int nthreads gmx_unused;
+
+#ifndef __clang_analyzer__
+ // cppcheck-suppress unreadVariable
+ nthreads = gmx_omp_nthreads_get(emntUpdate);
+#endif
/* Scale the positions */
+#pragma omp parallel for num_threads(nthreads) schedule(static)
for (n = start; n < start+nr_atoms; n++)
{
- if (cFREEZE)
+ int g;
+
+ if (cFREEZE == NULL)
+ {
+ g = 0;
+ }
+ else
{
g = cFREEZE[n];
}
if ((opts->tau_t[i] > 0) && (T > 0.0))
{
- reft = max(0.0, opts->ref_t[i]);
+ reft = std::max<real>(0, opts->ref_t[i]);
lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
- ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
+ ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
}
else
{
for (i = 0; (i < opts->ngtc); i++)
{
- reft = max(0.0, opts->ref_t[i]);
+ reft = std::max<real>(0, opts->ref_t[i]);
oldvxi = vxi[i];
vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
xi[i] += dt*(oldvxi + vxi[i])*0.5;
t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
{
- int n, i, j, d, ntgrp, ngtc, gc = 0;
+ int n, i, d, ngtc, gc = 0;
t_grp_tcstat *tcstat;
t_grpopts *opts;
gmx_int64_t step_eff;
- real ecorr, pcorr, dvdlcorr;
- real bmass, qmass, reft, kT, dt, nd;
- tensor dumpres, dumvir;
+ real dt;
double *scalefac, dtc;
int *trotter_seq;
- rvec sumv = {0, 0, 0}, consk;
+ rvec sumv = {0, 0, 0};
gmx_bool bCouple;
if (trotter_seqno <= ettTSEQ2)
extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
{
- int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
- t_grp_tcstat *tcstat;
+ int n, i, j, d, ngtc, nh;
t_grpopts *opts;
- real ecorr, pcorr, dvdlcorr;
- real bmass, qmass, reft, kT, dt, ndj, nd;
- tensor dumpres, dumvir;
+ real reft, kT, ndj, nd;
opts = &(ir->opts); /* just for ease of referencing */
ngtc = ir->opts.ngtc;
- nnhpres = state->nnhpres;
nh = state->nhchainlength;
if (ir->eI == eiMD)
{
if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
{
- reft = max(0.0, opts->ref_t[i]);
+ reft = std::max<real>(0, opts->ref_t[i]);
nd = opts->nrdf[i];
kT = BOLTZ*reft;
for (j = 0; j < nh; j++)
}
else
{
- reft = 0.0;
for (j = 0; j < nh; j++)
{
MassQ->Qinv[i*nh+j] = 0.0;
int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
{
- int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
- t_grp_tcstat *tcstat;
+ int i, j, nnhpres, nh;
t_grpopts *opts;
- real ecorr, pcorr, dvdlcorr;
- real bmass, qmass, reft, kT, dt, ndj, nd;
- tensor dumpres, dumvir;
+ real bmass, qmass, reft, kT;
int **trotter_seq;
opts = &(ir->opts); /* just for ease of referencing */
- ngtc = state->ngtc;
nnhpres = state->nnhpres;
nh = state->nhchainlength;
/* barostat temperature */
if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
{
- reft = max(0.0, opts->ref_t[0]);
+ reft = std::max<real>(0, opts->ref_t[0]);
kT = BOLTZ*reft;
for (i = 0; i < nnhpres; i++)
{
real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
{
- int i, j, nd, ndj, bmass, qmass, ngtcall;
- real ener_npt, reft, eta, kT, tau;
+ int i, j;
+ real nd, ndj;
+ real ener_npt, reft, kT;
double *ivxi, *ixi;
double *iQinv;
- real vol, dbaro, W, Q;
+ real vol;
int nh = state->nhchainlength;
ener_npt = 0;
ivxi = &state->nhpres_vxi[i*nh];
ixi = &state->nhpres_xi[i*nh];
iQinv = &(MassQ->QPinv[i*nh]);
- reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
+ reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
kT = BOLTZ * reft;
for (j = 0; j < nh; j++)
iQinv = &(MassQ->Qinv[i*nh]);
nd = ir->opts.nrdf[i];
- reft = max(ir->opts.ref_t[i], 0);
+ reft = std::max<real>(ir->opts.ref_t[i], 0);
kT = BOLTZ * reft;
- if (nd > 0)
+ if (nd > 0.0)
{
if (IR_NVT_TROTTER(ir))
{
}
else
{
- ndj = 1;
+ ndj = 1.0;
}
ener_npt += ndj*ixi[j]*kT;
}
case eannPERIODIC:
/* calculate time modulo the period */
pert = opts->anneal_time[i][npoints-1];
- n = t / pert;
+ n = static_cast<int>(t / pert);
thist = t - n*pert; /* modulo time */
/* Make sure rounding didn't get us outside the interval */
if (fabs(thist-pert) < GMX_REAL_EPS*100)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff