#include "gromacs/gmxlib/network.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
* started, but even when called, the prediction was always
* over-written by a subsequent call in the MD loop, so has been
* removed. */
-static void predict_shells(FILE* fplog,
- ArrayRef<RVec> x,
- ArrayRef<RVec> v,
- real dt,
- ArrayRef<const t_shell> shells,
- const real mass[],
- gmx_mtop_t* mtop,
- gmx_bool bInit)
+static void predict_shells(FILE* fplog,
+ ArrayRef<RVec> x,
+ ArrayRef<RVec> v,
+ real dt,
+ ArrayRef<const t_shell> shells,
+ gmx::ArrayRef<const real> mass,
+ gmx_bool bInit)
{
int m, n1, n2, n3;
real dt_1, fudge, tm, m1, m2, m3;
- GMX_RELEASE_ASSERT(mass || mtop, "Must have masses or a way to look them up");
-
/* We introduce a fudge factor for performance reasons: with this choice
* the initial force on the shells is about a factor of two lower than
* without
dt_1 = fudge * dt;
}
- int molb = 0;
for (const t_shell& shell : shells)
{
const int s1 = shell.shellIndex;
case 2:
n1 = shell.nucl1;
n2 = shell.nucl2;
- if (mass)
- {
- m1 = mass[n1];
- m2 = mass[n2];
- }
- else
- {
- /* Not the correct masses with FE, but it is just a prediction... */
- m1 = mtopGetAtomMass(mtop, n1, &molb);
- m2 = mtopGetAtomMass(mtop, n2, &molb);
- }
+ m1 = mass[n1];
+ m2 = mass[n2];
tm = dt_1 / (m1 + m2);
for (m = 0; (m < DIM); m++)
{
n1 = shell.nucl1;
n2 = shell.nucl2;
n3 = shell.nucl3;
- if (mass)
- {
- m1 = mass[n1];
- m2 = mass[n2];
- m3 = mass[n3];
- }
- else
- {
- /* Not the correct masses with FE, but it is just a prediction... */
- m1 = mtopGetAtomMass(mtop, n1, &molb);
- m2 = mtopGetAtomMass(mtop, n2, &molb);
- m3 = mtopGetAtomMass(mtop, n3, &molb);
- }
+ m1 = mass[n1];
+ m2 = mass[n2];
+ m3 = mass[n3];
tm = dt_1 / (m1 + m2 + m3);
for (m = 0; (m < DIM); m++)
{
}
gmx_shellfc_t* init_shell_flexcon(FILE* fplog,
- const gmx_mtop_t* mtop,
+ const gmx_mtop_t& mtop,
int nflexcon,
int nstcalcenergy,
bool usingDomainDecomposition,
#define NBT asize(bondtypes)
const gmx_ffparams_t* ffparams;
- const std::array<int, eptNR> numParticles = gmx_mtop_particletype_count(*mtop);
+ const gmx::EnumerationArray<ParticleType, int> numParticles = gmx_mtop_particletype_count(mtop);
if (fplog)
{
/* Print the number of each particle type */
- int pType = 0;
- for (const auto& n : numParticles)
+ for (const auto entry : gmx::keysOf(numParticles))
{
- if (n != 0)
+ const int number = numParticles[entry];
+ if (number != 0)
{
- fprintf(fplog, "There are: %d %ss\n", n, ptype_str[pType]);
+ fprintf(fplog, "There are: %d %ss\n", number, enumValueToString(entry));
}
- pType++;
}
}
- nshell = numParticles[eptShell];
+ nshell = numParticles[ParticleType::Shell];
if (nshell == 0 && nflexcon == 0)
{
"not supported in combination with shell particles.\nPlease make a new tpr file.",
nstcalcenergy);
}
- if (usingDomainDecomposition)
+ if (nshell > 0 && usingDomainDecomposition)
{
gmx_fatal(
FARGS,
/* We have shells: fill the shell data structure */
/* Global system sized array, this should be avoided */
- std::vector<int> shell_index(mtop->natoms);
+ std::vector<int> shell_index(mtop.natoms);
nshell = 0;
- for (const AtomProxy atomP : AtomRange(*mtop))
+ for (const AtomProxy atomP : AtomRange(mtop))
{
const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
- if (local.ptype == eptShell)
+ if (local.ptype == ParticleType::Shell)
{
shell_index[i] = nshell++;
}
std::vector<t_shell> shell(nshell);
- ffparams = &mtop->ffparams;
+ ffparams = &mtop.ffparams;
/* Now fill the structures */
/* TODO: See if we can use update groups that cover shell constructions */
shfc->bInterCG = FALSE;
ns = 0;
a_offset = 0;
- for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
+ for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
{
- const gmx_molblock_t* molb = &mtop->molblock[mb];
- const gmx_moltype_t* molt = &mtop->moltype[molb->type];
+ const gmx_molblock_t& molb = mtop.molblock[mb];
+ const gmx_moltype_t& molt = mtop.moltype[molb.type];
- const t_atom* atom = molt->atoms.atom;
- for (mol = 0; mol < molb->nmol; mol++)
+ const t_atom* atom = molt.atoms.atom;
+ for (mol = 0; mol < molb.nmol; mol++)
{
for (j = 0; (j < NBT); j++)
{
- const int* ia = molt->ilist[bondtypes[j]].iatoms.data();
- for (i = 0; (i < molt->ilist[bondtypes[j]].size());)
+ const int* ia = molt.ilist[bondtypes[j]].iatoms.data();
+ for (i = 0; (i < molt.ilist[bondtypes[j]].size());)
{
type = ia[0];
ftype = ffparams->functype[type];
case F_CUBICBONDS:
case F_POLARIZATION:
case F_ANHARM_POL:
- if (atom[ia[1]].ptype == eptShell)
+ if (atom[ia[1]].ptype == ParticleType::Shell)
{
aS = ia[1];
aN = ia[2];
}
- else if (atom[ia[2]].ptype == eptShell)
+ else if (atom[ia[2]].ptype == ParticleType::Shell)
{
aS = ia[2];
aN = ia[1];
aS + 1,
mb + 1);
}
- shell[nsi].k += gmx::square(qS) * ONE_4PI_EPS0
+ shell[nsi].k += gmx::square(qS) * gmx::c_one4PiEps0
/ ffparams->iparams[type].polarize.alpha;
break;
case F_WATER_POL:
+ ffparams->iparams[type].wpol.al_y
+ ffparams->iparams[type].wpol.al_z)
/ 3.0;
- shell[nsi].k += gmx::square(qS) * ONE_4PI_EPS0 / alpha;
+ shell[nsi].k += gmx::square(qS) * gmx::c_one4PiEps0 / alpha;
break;
default: gmx_fatal(FARGS, "Death Horror: %s, %d", __FILE__, __LINE__);
}
i += nra + 1;
}
}
- a_offset += molt->atoms.nr;
+ a_offset += molt.atoms.nr;
}
/* Done with this molecule type */
}
return shfc;
}
-void gmx::make_local_shells(const t_commrec* cr, const t_mdatoms* md, gmx_shellfc_t* shfc)
+void gmx::make_local_shells(const t_commrec* cr, const t_mdatoms& md, gmx_shellfc_t* shfc)
{
int a0, a1;
gmx_domdec_t* dd = nullptr;
- if (DOMAINDECOMP(cr))
+ if (haveDDAtomOrdering(*cr))
{
dd = cr->dd;
a0 = 0;
std::vector<t_shell>& shells = shfc->shells;
shells.clear();
+ auto* ptype = md.ptype;
for (int i = a0; i < a1; i++)
{
- if (md->ptype[i] == eptShell)
+ if (ptype[i] == ParticleType::Shell)
{
if (dd)
{
const t_commrec* cr,
int dd_ac1,
int64_t step,
- const t_mdatoms* md,
+ const t_mdatoms& md,
int end,
ArrayRefWithPadding<RVec> xOld,
ArrayRef<RVec> x_init,
ArrayRef<const real> lambda,
real* dvdlambda)
{
- double dt, w_dt;
- int n, d;
- unsigned short* ptype;
+ double dt, w_dt;
+ int n, d;
- if (DOMAINDECOMP(cr))
+ if (haveDDAtomOrdering(*cr))
{
n = dd_ac1;
}
rvec* x_old = as_rvec_array(xOld.paddedArrayRef().data());
rvec* x = as_rvec_array(xCurrent.paddedArrayRef().data());
- ptype = md->ptype;
-
- dt = ir->delta_t;
+ auto* ptype = md.ptype;
+ auto invmass = gmx::arrayRefFromArray(md.invmass, md.nr);
+ dt = ir->delta_t;
- /* Does NOT work with freeze groups (yet) */
+ /* Does NOT work with freeze or acceleration groups (yet) */
for (n = 0; n < end; n++)
{
- w_dt = md->invmass[n] * dt;
+ w_dt = invmass[n] * dt;
for (d = 0; d < DIM; d++)
{
- if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
+ if ((ptype[n] != ParticleType::VSite) && (ptype[n] != ParticleType::Shell))
{
xnold[n][d] = x[n][d] - (x_init[n][d] - x_old[n][d]);
xnew[n][d] = 2 * x[n][d] - x_old[n][d] + f[n][d] * w_dt * dt;
shfc->adir_xnold.arrayRefWithPadding(),
{},
box,
- lambda[efptBONDED],
- &(dvdlambda[efptBONDED]),
+ lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)],
+ &(dvdlambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)]),
{},
computeVirial,
nullptr,
shfc->adir_xnew.arrayRefWithPadding(),
{},
box,
- lambda[efptBONDED],
- &(dvdlambda[efptBONDED]),
+ lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)],
+ &(dvdlambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)]),
{},
computeVirial,
nullptr,
for (d = 0; d < DIM; d++)
{
xnew[n][d] = -(2 * x[n][d] - xnold[n][d] - xnew[n][d]) / gmx::square(dt)
- - f[n][d] * md->invmass[n];
+ - f[n][d] * invmass[n];
}
clear_rvec(acc_dir[n]);
}
shfc->adir_xnew.arrayRefWithPadding(),
acc_dir,
box,
- lambda[efptBONDED],
- &(dvdlambda[efptBONDED]),
+ lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)],
+ &(dvdlambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)]),
{},
computeVirial,
nullptr,
const history_t* hist,
gmx::ForceBuffersView* f,
tensor force_vir,
- const t_mdatoms* md,
+ const t_mdatoms& md,
+ CpuPpLongRangeNonbondeds* longRangeNonbondeds,
t_nrnb* nrnb,
- gmx_wallcycle_t wcycle,
+ gmx_wallcycle* wcycle,
gmx_shellfc_t* shfc,
t_forcerec* fr,
gmx::MdrunScheduleWorkload* runScheduleWork,
real dum = 0;
char sbuf[22];
int nat, dd_ac0, dd_ac1 = 0, i;
- int homenr = md->homenr, end = homenr;
+ int homenr = md.homenr, end = homenr;
int d, Min = 0, count = 0;
#define Try (1 - Min) /* At start Try = 1 */
ArrayRef<t_shell> shells = shfc->shells;
const int nflexcon = shfc->nflexcon;
- if (DOMAINDECOMP(cr))
+ if (haveDDAtomOrdering(*cr))
{
- nat = dd_natoms_vsite(cr->dd);
+ nat = dd_natoms_vsite(*cr->dd);
if (nflexcon > 0)
{
- dd_get_constraint_range(cr->dd, &dd_ac0, &dd_ac1);
+ dd_get_constraint_range(*cr->dd, &dd_ac0, &dd_ac1);
nat = std::max(nat, dd_ac1);
}
}
ArrayRef<RVec> x = xPadded.unpaddedArrayRef();
ArrayRef<RVec> v = vPadded.unpaddedArrayRef();
- if (bDoNS && inputrec->pbcType != PbcType::No && !DOMAINDECOMP(cr))
+ if (bDoNS && inputrec->pbcType != PbcType::No && !haveDDAtomOrdering(*cr))
{
/* This is the only time where the coordinates are used
* before do_force is called, which normally puts all
* charge groups in the box.
*/
put_atoms_in_box_omp(
- fr->pbcType, box, x.subArray(0, md->homenr), gmx_omp_nthreads_get(emntDefault));
+ fr->pbcType, box, x.subArray(0, md.homenr), gmx_omp_nthreads_get(ModuleMultiThread::Default));
}
if (nflexcon)
}
}
+ auto massT = gmx::arrayRefFromArray(md.massT, md.nr);
/* Do a prediction of the shell positions, when appropriate.
* Without velocities (EM, NM, BD) we only do initial prediction.
*/
if (shfc->predictShells && !bCont && (EI_STATE_VELOCITY(inputrec->eI) || bInit))
{
- predict_shells(fplog, x, v, inputrec->delta_t, shells, md->massT, nullptr, bInit);
+ predict_shells(fplog, x, v, inputrec->delta_t, shells, massT, bInit);
}
/* Calculate the forces first time around */
hist,
&forceViewInit,
force_vir,
- md,
+ &md,
enerd,
lambda,
fr,
mu_tot,
t,
nullptr,
+ longRangeNonbondeds,
(bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags,
ddBalanceRegionHandler);
for (i = 0; i < end; i++)
{
- sf_dir += md->massT[i] * norm2(shfc->acc_dir[i]);
+ sf_dir += massT[i] * norm2(shfc->acc_dir[i]);
}
}
- accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals);
+ accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals.get());
Epot[Min] = enerd->term[F_EPOT];
df[Min] = rms_force(cr, forceWithPadding[Min].paddedArrayRef(), shells, nflexcon, &sf_dir, &Epot[Min]);
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "force0", as_rvec_array(force[Min].data()), md->nr);
+ pr_rvecs(debug, 0, "force0", as_rvec_array(force[Min].data()), md.nr);
}
if (!shells.empty() || nflexcon > 0)
{
if (vsite)
{
- vsite->construct(pos[Min], inputrec->delta_t, v, box);
+ vsite->construct(pos[Min], v, box, gmx::VSiteOperation::PositionsAndVelocities);
}
if (nflexcon)
hist,
&forceViewTry,
force_vir,
- md,
+ &md,
enerd,
lambda,
fr,
mu_tot,
t,
nullptr,
+ longRangeNonbondeds,
shellfc_flags,
ddBalanceRegionHandler);
- accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals);
+ accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals.get());
if (gmx_debug_at)
{
pr_rvecs(debug, 0, "RELAX: force[Min]", as_rvec_array(force[Min].data()), homenr);
ArrayRef<const RVec> acc_dir = shfc->acc_dir;
for (i = 0; i < end; i++)
{
- sf_dir += md->massT[i] * norm2(acc_dir[i]);
+ sf_dir += massT[i] * norm2(acc_dir[i]);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff