From 21dcd02d31d2fd2786b527c84ad4697543a4bd55 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 5 May 2021 15:45:28 +0000 Subject: [PATCH] Convert t_oriresdata to C++ --- admin/lsan-suppressions.txt | 1 - src/gromacs/listed_forces/listed_forces.cpp | 10 +- src/gromacs/listed_forces/orires.cpp | 328 +++++++++----------- src/gromacs/listed_forces/orires.h | 16 - src/gromacs/mdlib/energyoutput.cpp | 22 +- src/gromacs/mdrun/md.cpp | 5 +- src/gromacs/mdrun/runner.cpp | 11 +- src/gromacs/mdtypes/fcdata.h | 126 ++++++-- 8 files changed, 275 insertions(+), 244 deletions(-) diff --git a/admin/lsan-suppressions.txt b/admin/lsan-suppressions.txt index d064c4cbea..ce8a33ee2b 100644 --- a/admin/lsan-suppressions.txt +++ b/admin/lsan-suppressions.txt @@ -24,7 +24,6 @@ leak:init_df_history leak:init_disres leak:init_edsam leak:init_ekinstate -leak:init_orires leak:init_rot leak:init_swapcoords leak:make_bondeds_zone diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp index d8511e56b4..be7ba91f0f 100644 --- a/src/gromacs/listed_forces/listed_forces.cpp +++ b/src/gromacs/listed_forces/listed_forces.cpp @@ -492,7 +492,7 @@ real calc_one_bond(int thread, gmx::arrayRefFromArray(md->chargeA, md->nr), fcd, fcd->disres, - fcd->orires, + fcd->orires.get(), global_atom_index, flavor); } @@ -622,9 +622,9 @@ static void calcBondedForces(const InteractionDefinitions& idef, bool ListedForces::haveRestraints(const t_fcdata& fcdata) const { - GMX_ASSERT(fcdata.orires && fcdata.disres, "NMR restraints objects should be set up"); + GMX_ASSERT(fcdata.disres, "NMR distance restraint object should be set up"); - return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty() || fcdata.orires->nr > 0 + return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty() || fcdata.orires || fcdata.disres->nres > 0); } @@ -846,7 +846,7 @@ void ListedForces::calculate(struct gmx_wallcycle* wcycle, } /* Do pre force calculation stuff which might require communication */ - if (fcdata->orires->nr > 0) + if (fcdata->orires) { GMX_ASSERT(!xWholeMolecules.empty(), "Need whole molecules for orienation restraints"); enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, @@ -857,7 +857,7 @@ void ListedForces::calculate(struct gmx_wallcycle* wcycle, xWholeMolecules, x, fr->bMolPBC ? pbc : nullptr, - fcdata->orires, + fcdata->orires.get(), hist); } if (fcdata->disres->nres > 0) diff --git a/src/gromacs/listed_forces/orires.cpp b/src/gromacs/listed_forces/orires.cpp index 1ec9cef352..626d3996e7 100644 --- a/src/gromacs/listed_forces/orires.cpp +++ b/src/gromacs/listed_forces/orires.cpp @@ -59,7 +59,7 @@ #include "gromacs/topology/mtop_util.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/arrayref.h" -#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/exceptions.h" #include "gromacs/utility/pleasecite.h" #include "gromacs/utility/smalloc.h" @@ -69,87 +69,76 @@ using gmx::RVec; // TODO This implementation of ensemble orientation restraints is nasty because // a user can't just do multi-sim with single-sim orientation restraints. -void init_orires(FILE* fplog, - const gmx_mtop_t& mtop, - const t_inputrec* ir, - const t_commrec* cr, - const gmx_multisim_t* ms, - t_state* globalState, - t_oriresdata* od) +t_oriresdata::t_oriresdata(FILE* fplog, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const t_commrec* cr, + const gmx_multisim_t* ms, + t_state* globalState) : + numRestraints(gmx_mtop_ftype_count(mtop, F_ORIRES)) { - od->nr = gmx_mtop_ftype_count(mtop, F_ORIRES); - if (0 == od->nr) - { - /* Not doing orientation restraints */ - return; - } + GMX_RELEASE_ASSERT(numRestraints > 0, + "orires() should only be called with orientation restraints present"); const int numFitParams = 5; - if (od->nr <= numFitParams) + if (numRestraints <= numFitParams) { - gmx_fatal(FARGS, - "The system has %d orientation restraints, but at least %d are required, since " - "there are %d fitting parameters.", - od->nr, - numFitParams + 1, - numFitParams); + const std::string mesg = gmx::formatString( + "The system has %d orientation restraints, but at least %d are required, since " + "there are %d fitting parameters.", + numRestraints, + numFitParams + 1, + numFitParams); + GMX_THROW(gmx::InvalidInputError(mesg)); } - if (ir->bPeriodicMols) + if (ir.bPeriodicMols) { /* Since we apply fitting, we need to make molecules whole and this * can not be done when periodic molecules are present. */ - gmx_fatal(FARGS, - "Orientation restraints can not be applied when periodic molecules are present " - "in the system"); + GMX_THROW(gmx::InvalidInputError( + "Orientation restraints can not be applied when periodic molecules are present " + "in the system")); } - if (PAR(cr)) + if (cr && PAR(cr)) { - gmx_fatal(FARGS, - "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, " - "if possible."); + GMX_THROW(gmx::InvalidInputError( + "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, " + "if possible.")); } - GMX_RELEASE_ASSERT(globalState != nullptr, "We need a valid global state in init_orires"); + GMX_RELEASE_ASSERT(globalState != nullptr, "We need a valid global state in t_oriresdata()"); - od->fc = ir->orires_fc; - od->nex = 0; - od->S = nullptr; - od->M = nullptr; - od->eig = nullptr; - od->v = nullptr; + fc = ir.orires_fc; + numExperiments = 0; - int* nr_ex = nullptr; - int typeMin = INT_MAX; - int typeMax = 0; + std::vector nr_ex; + typeMin = INT_MAX; + int typeMax = 0; for (const auto il : IListRange(mtop)) { const int numOrires = il.list()[F_ORIRES].size(); if (il.nmol() > 1 && numOrires > 0) { - gmx_fatal(FARGS, - "Found %d copies of a molecule with orientation restrains while the current " - "code only supports a single copy. If you want to ensemble average, run " - "multiple copies of the system using the multi-sim feature of mdrun.", - il.nmol()); + const std::string mesg = gmx::formatString( + "Found %d copies of a molecule with orientation restrains while the current " + "code only supports a single copy. If you want to ensemble average, run " + "multiple copies of the system using the multi-sim feature of mdrun.", + il.nmol()); + GMX_THROW(gmx::InvalidInputError(mesg)); } for (int i = 0; i < numOrires; i += 3) { int type = il.list()[F_ORIRES].iatoms[i]; int ex = mtop.ffparams.iparams[type].orires.ex; - if (ex >= od->nex) + if (ex >= numExperiments) { - srenew(nr_ex, ex + 1); - for (int j = od->nex; j < ex + 1; j++) - { - nr_ex[j] = 0; - } - od->nex = ex + 1; + nr_ex.resize(ex + 1, 0); + numExperiments = ex + 1; } - GMX_ASSERT(nr_ex, "Check for allocated nr_ex to keep the static analyzer happy"); nr_ex[ex]++; typeMin = std::min(typeMin, type); @@ -158,77 +147,77 @@ void init_orires(FILE* fplog, } /* With domain decomposition we use the type index for indexing in global arrays */ GMX_RELEASE_ASSERT( - typeMax - typeMin + 1 == od->nr, + typeMax - typeMin + 1 == numRestraints, "All orientation restraint parameter entries in the topology should be consecutive"); - /* Store typeMin so we can index array with the type offset */ - od->typeMin = typeMin; - snew(od->S, od->nex); + snew(orderTensors, numExperiments); /* When not doing time averaging, the instaneous and time averaged data * are indentical and the pointers can point to the same memory. */ - snew(od->Dinsl, od->nr); + snew(DTensors, numRestraints); if (ms) { - snew(od->Dins, od->nr); + snew(DTensorsEnsembleAv, numRestraints); } else { - od->Dins = od->Dinsl; + DTensorsEnsembleAv = DTensors; } - if (ir->orires_tau == 0) + if (ir.orires_tau == 0) { - od->Dtav = od->Dins; - od->edt = 0.0; - od->edt_1 = 1.0; + DTensorsTimeAndEnsembleAv = DTensorsEnsembleAv; + edt = 0.0; + edt_1 = 1.0; } else { - snew(od->Dtav, od->nr); - od->edt = std::exp(-ir->delta_t / ir->orires_tau); - od->edt_1 = 1.0 - od->edt; + snew(DTensorsTimeAndEnsembleAv, numRestraints); + edt = std::exp(-ir.delta_t / ir.orires_tau); + edt_1 = 1.0 - edt; /* Extend the state with the orires history */ globalState->flags |= enumValueToBitMask(StateEntry::OrireInitF); globalState->hist.orire_initf = 1; globalState->flags |= enumValueToBitMask(StateEntry::OrireDtav); - globalState->hist.orire_Dtav.resize(od->nr * 5); + globalState->hist.orire_Dtav.resize(numRestraints * 5); } - snew(od->oinsl, od->nr); + orientations.resize(numRestraints); if (ms) { - snew(od->oins, od->nr); + orientationsEnsembleAvBuffer.resize(numRestraints); + orientationsEnsembleAv = orientationsEnsembleAvBuffer; } else { - od->oins = od->oinsl; + orientationsEnsembleAv = orientations; } - if (ir->orires_tau == 0) + if (ir.orires_tau == 0) { - od->otav = od->oins; + orientationsTimeAndEnsembleAv = orientationsEnsembleAv; } else { - snew(od->otav, od->nr); + orientationsTimeAndEnsembleAvBuffer.resize(numRestraints); + orientationsTimeAndEnsembleAv = orientationsTimeAndEnsembleAvBuffer; } - snew(od->tmpEq, od->nex); + tmpEq.resize(numExperiments); - od->nref = 0; + numReferenceAtoms = 0; for (int i = 0; i < mtop.natoms; i++) { if (getGroupType(mtop.groups, SimulationAtomGroupType::OrientationRestraintsFit, i) == 0) { - od->nref++; + numReferenceAtoms++; } } - snew(od->mref, od->nref); - snew(od->xref, od->nref); - snew(od->xtmp, od->nref); + mref.resize(numReferenceAtoms); + xref.resize(numReferenceAtoms); + xtmp.resize(numReferenceAtoms); - snew(od->eig, od->nex * 12); + eigenOutput.resize(numExperiments * c_numEigenRealsPerExperiment); /* Determine the reference structure on the master node. * Copy it to the other nodes after checking multi compatibility, @@ -246,76 +235,82 @@ void init_orires(FILE* fplog, || mtop.groups.groupNumbers[SimulationAtomGroupType::OrientationRestraintsFit][i] == 0) { /* Not correct for free-energy with changing masses */ - od->mref[j] = local.m; + mref[j] = local.m; // Note that only one rank per sim is supported. if (isMasterSim(ms)) { - copy_rvec(x[i], od->xref[j]); + copy_rvec(x[i], xref[j]); for (int d = 0; d < DIM; d++) { - com[d] += od->mref[j] * x[i][d]; + com[d] += mref[j] * x[i][d]; } } - mtot += od->mref[j]; + mtot += mref[j]; j++; } } svmul(1.0 / mtot, com, com); if (isMasterSim(ms)) { - for (int j = 0; j < od->nref; j++) + for (int j = 0; j < numReferenceAtoms; j++) { - rvec_dec(od->xref[j], com); + rvec_dec(xref[j], com); } } - fprintf(fplog, "Found %d orientation experiments\n", od->nex); - for (int i = 0; i < od->nex; i++) + if (fplog) { - fprintf(fplog, " experiment %d has %d restraints\n", i + 1, nr_ex[i]); - } - - sfree(nr_ex); + fprintf(fplog, "Found %d orientation experiments\n", numExperiments); + for (int i = 0; i < numExperiments; i++) + { + fprintf(fplog, " experiment %d has %d restraints\n", i + 1, nr_ex[i]); + } - fprintf(fplog, " the fit group consists of %d atoms and has total mass %g\n", od->nref, mtot); + fprintf(fplog, " the fit group consists of %d atoms and has total mass %g\n", numReferenceAtoms, mtot); + } if (ms) { - fprintf(fplog, " the orientation restraints are ensemble averaged over %d systems\n", ms->numSimulations_); + if (fplog) + { + fprintf(fplog, + " the orientation restraints are ensemble averaged over %d systems\n", + ms->numSimulations_); + } - check_multi_int(fplog, ms, od->nr, "the number of orientation restraints", FALSE); - check_multi_int(fplog, ms, od->nref, "the number of fit atoms for orientation restraining", FALSE); - check_multi_int(fplog, ms, ir->nsteps, "nsteps", FALSE); + check_multi_int(fplog, ms, numRestraints, "the number of orientation restraints", FALSE); + check_multi_int( + fplog, ms, numReferenceAtoms, "the number of fit atoms for orientation restraining", FALSE); + check_multi_int(fplog, ms, ir.nsteps, "nsteps", FALSE); /* Copy the reference coordinates from the master to the other nodes */ - gmx_sum_sim(DIM * od->nref, od->xref[0], ms); + gmx_sum_sim(DIM * numReferenceAtoms, xref[0], ms); } please_cite(fplog, "Hess2003"); } -void diagonalize_orires_tensors(t_oriresdata* od) +t_oriresdata::~t_oriresdata() { - if (od->M == nullptr) + sfree(orderTensors); + if (DTensorsTimeAndEnsembleAv != DTensorsEnsembleAv) { - snew(od->M, DIM); - for (int i = 0; i < DIM; i++) - { - snew(od->M[i], DIM); - } - snew(od->eig_diag, DIM); - snew(od->v, DIM); - for (int i = 0; i < DIM; i++) - { - snew(od->v[i], DIM); - } + sfree(DTensorsTimeAndEnsembleAv); } + if (DTensorsEnsembleAv != DTensors) + { + sfree(DTensorsEnsembleAv); + } + sfree(DTensors); +} - for (int ex = 0; ex < od->nex; ex++) +void diagonalize_orires_tensors(t_oriresdata* od) +{ + for (int ex = 0; ex < od->numExperiments; ex++) { /* Rotate the S tensor back to the reference frame */ matrix S, TMP; - mmul(od->R, od->S[ex], TMP); - mtmul(TMP, od->R, S); + mmul(od->rotationMatrix, od->orderTensors[ex], TMP); + mtmul(TMP, od->rotationMatrix, S); for (int i = 0; i < DIM; i++) { for (int j = 0; j < DIM; j++) @@ -324,8 +319,7 @@ void diagonalize_orires_tensors(t_oriresdata* od) } } - int nrot; - jacobi(od->M, DIM, od->eig_diag, od->v, &nrot); + jacobi(od->M, od->eig_diag, od->v); int ord[DIM]; for (int i = 0; i < DIM; i++) @@ -347,13 +341,14 @@ void diagonalize_orires_tensors(t_oriresdata* od) for (int i = 0; i < DIM; i++) { - od->eig[ex * 12 + i] = od->eig_diag[ord[i]]; + od->eigenOutput[ex * t_oriresdata::c_numEigenRealsPerExperiment + i] = od->eig_diag[ord[i]]; } for (int i = 0; i < DIM; i++) { for (int j = 0; j < DIM; j++) { - od->eig[ex * 12 + 3 + 3 * i + j] = od->v[j][ord[i]]; + od->eigenOutput[ex * t_oriresdata::c_numEigenRealsPerExperiment + 3 + 3 * i + j] = + od->v[j][ord[i]]; } } } @@ -361,13 +356,11 @@ void diagonalize_orires_tensors(t_oriresdata* od) void print_orires_log(FILE* log, t_oriresdata* od) { - real* eig; - diagonalize_orires_tensors(od); - for (int ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->numExperiments; ex++) { - eig = od->eig + ex * 12; + const real* eig = od->eigenOutput.data() + ex * t_oriresdata::c_numEigenRealsPerExperiment; fprintf(log, " Orientation experiment %d:\n", ex + 1); fprintf(log, " order parameter: %g\n", eig[0]); for (int i = 0; i < DIM; i++) @@ -394,30 +387,19 @@ real calc_orires_dev(const gmx_multisim_t* ms, t_oriresdata* od, const history_t* hist) { - int nref; - real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev; - OriresMatEq* matEq; - real* mref; - double mtot; - rvec * xref, *xtmp, com, r_unrot, r; - gmx_bool bTAV; - const real two_thr = 2.0 / 3.0; - - if (od->nr == 0) - { - /* This means that this is not the master node */ - gmx_fatal(FARGS, - "Orientation restraints are only supported on the master rank, use fewer ranks"); - } - - bTAV = (od->edt != 0); - edt = od->edt; - edt_1 = od->edt_1; - matEq = od->tmpEq; - nref = od->nref; - mref = od->mref; - xref = od->xref; - xtmp = od->xtmp; + real invn, pfac, r2, invr, corrfac, wsv2, sw, dev; + double mtot; + rvec com, r_unrot, r; + const real two_thr = 2.0 / 3.0; + + const bool bTAV = (od->edt != 0); + const real edt = od->edt; + const real edt_1 = od->edt_1; + gmx::ArrayRef matEq = od->tmpEq; + const int nref = od->numReferenceAtoms; + gmx::ArrayRef mref = od->mref; + gmx::ArrayRef xref = od->xref; + gmx::ArrayRef xtmp = od->xtmp; if (bTAV) { @@ -467,7 +449,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, rvec_dec(xtmp[j], com); } /* Calculate the rotation matrix to rotate x to the reference orientation */ - calc_fit_R(DIM, nref, mref, xref, xtmp, od->R); + calc_fit_R(DIM, nref, mref.data(), as_rvec_array(xref.data()), as_rvec_array(xtmp.data()), od->rotationMatrix); for (int fa = 0; fa < nfa; fa += 3) { @@ -481,7 +463,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, { rvec_sub(x[forceatoms[fa + 1]], x[forceatoms[fa + 2]], r_unrot); } - mvmul(od->R, r_unrot, r); + mvmul(od->rotationMatrix, r_unrot, r); r2 = norm2(r); invr = gmx::invsqrt(r2); /* Calculate the prefactor for the D tensor, this includes the factor 3! */ @@ -490,7 +472,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, { pfac *= invr; } - rvec5& Dinsl = od->Dinsl[restraintIndex]; + rvec5& Dinsl = od->DTensors[restraintIndex]; Dinsl[0] = pfac * (2 * r[0] * r[0] + r[1] * r[1] - r2); Dinsl[1] = pfac * (2 * r[0] * r[1]); Dinsl[2] = pfac * (2 * r[0] * r[2]); @@ -501,18 +483,18 @@ real calc_orires_dev(const gmx_multisim_t* ms, { for (int i = 0; i < 5; i++) { - od->Dins[restraintIndex][i] = Dinsl[i] * invn; + od->DTensorsEnsembleAv[restraintIndex][i] = Dinsl[i] * invn; } } } if (ms) { - gmx_sum_sim(5 * od->nr, od->Dins[0], ms); + gmx_sum_sim(5 * od->numRestraints, od->DTensorsEnsembleAv[0], ms); } /* Calculate the order tensor S for each experiment via optimization */ - for (int ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->numExperiments; ex++) { for (int i = 0; i < 5; i++) { @@ -528,17 +510,17 @@ real calc_orires_dev(const gmx_multisim_t* ms, { const int type = forceatoms[fa]; const int restraintIndex = type - od->typeMin; - rvec5& Dtav = od->Dtav[restraintIndex]; + rvec5& Dtav = od->DTensorsTimeAndEnsembleAv[restraintIndex]; if (bTAV) { - /* Here we update Dtav in t_fcdata using the data in history_t. + /* Here we update DTensorsTimeAndEnsembleAv in t_fcdata using the data in history_t. * Thus the results stay correct when this routine * is called multiple times. */ for (int i = 0; i < 5; i++) { Dtav[i] = edt * hist->orire_Dtav[restraintIndex * 5 + i] - + edt_1 * od->Dins[restraintIndex][i]; + + edt_1 * od->DTensorsEnsembleAv[restraintIndex][i]; } } @@ -555,7 +537,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, } } /* Now we have all the data we can calculate S */ - for (int ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->numExperiments; ex++) { OriresMatEq& eq = matEq[ex]; /* Correct corrfac and copy one half of T to the other half */ @@ -571,7 +553,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, } m_inv_gen(&eq.mat[0][0], 5, &eq.mat[0][0]); /* Calculate the orientation tensor S for this experiment */ - matrix& S = od->S[ex]; + matrix& S = od->orderTensors[ex]; S[0][0] = 0; S[0][1] = 0; S[0][2] = 0; @@ -591,7 +573,7 @@ real calc_orires_dev(const gmx_multisim_t* ms, S[2][2] = -S[0][0] - S[1][1]; } - const matrix* S = od->S; + const matrix* S = od->orderTensors; wsv2 = 0; sw = 0; @@ -602,15 +584,15 @@ real calc_orires_dev(const gmx_multisim_t* ms, const int restraintIndex = type - od->typeMin; const int ex = ip[type].orires.ex; - const rvec5& Dtav = od->Dtav[restraintIndex]; - od->otav[restraintIndex] = + const rvec5& Dtav = od->DTensorsTimeAndEnsembleAv[restraintIndex]; + od->orientationsTimeAndEnsembleAv[restraintIndex] = two_thr * corrfac * (S[ex][0][0] * Dtav[0] + S[ex][0][1] * Dtav[1] + S[ex][0][2] * Dtav[2] + S[ex][1][1] * Dtav[3] + S[ex][1][2] * Dtav[4]); if (bTAV) { - const rvec5& Dins = od->Dins[restraintIndex]; - od->oins[restraintIndex] = + const rvec5& Dins = od->DTensorsEnsembleAv[restraintIndex]; + od->orientationsEnsembleAv[restraintIndex] = two_thr * (S[ex][0][0] * Dins[0] + S[ex][0][1] * Dins[1] + S[ex][0][2] * Dins[2] + S[ex][1][1] * Dins[3] + S[ex][1][2] * Dins[4]); @@ -620,14 +602,14 @@ real calc_orires_dev(const gmx_multisim_t* ms, /* When ensemble averaging is used recalculate the local orientation * for output to the energy file. */ - const rvec5& Dinsl = od->Dinsl[restraintIndex]; - od->oinsl[restraintIndex] = + const rvec5& Dinsl = od->DTensors[restraintIndex]; + od->orientations[restraintIndex] = two_thr * (S[ex][0][0] * Dinsl[0] + S[ex][0][1] * Dinsl[1] + S[ex][0][2] * Dinsl[2] + S[ex][1][1] * Dinsl[3] + S[ex][1][2] * Dinsl[4]); } - dev = od->otav[restraintIndex] - ip[type].orires.obs; + dev = od->orientationsTimeAndEnsembleAv[restraintIndex] - ip[type].orires.obs; wsv2 += ip[type].orires.kfac * gmx::square(dev); sw += ip[type].orires.kfac; @@ -635,11 +617,11 @@ real calc_orires_dev(const gmx_multisim_t* ms, od->rmsdev = std::sqrt(wsv2 / sw); /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */ - for (int ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->numExperiments; ex++) { matrix RS; - tmmul(od->R, od->S[ex], RS); - mmul(RS, od->R, od->S[ex]); + tmmul(od->rotationMatrix, od->orderTensors[ex], RS); + mmul(RS, od->rotationMatrix, od->orderTensors[ex]); } return od->rmsdev; @@ -701,7 +683,7 @@ real orires(int nfa, ex = ip[type].orires.ex; power = ip[type].orires.power; fc = smooth_fc * ip[type].orires.kfac; - dev = oriresdata->otav[restraintIndex] - ip[type].orires.obs; + dev = oriresdata->orientationsTimeAndEnsembleAv[restraintIndex] - ip[type].orires.obs; /* NOTE: * there is no real potential when time averaging is applied @@ -711,7 +693,7 @@ real orires(int nfa, if (bTAV) { /* Calculate the force as the sqrt of tav times instantaneous */ - devins = oriresdata->oins[restraintIndex] - ip[type].orires.obs; + devins = oriresdata->orientationsEnsembleAv[restraintIndex] - ip[type].orires.obs; if (dev * devins <= 0) { dev = 0; @@ -731,7 +713,7 @@ real orires(int nfa, { pfac *= invr; } - mvmul(oriresdata->S[ex], r, Sr); + mvmul(oriresdata->orderTensors[ex], r, Sr); for (int i = 0; i < DIM; i++) { fij[i] = -pfac * dev * (4 * Sr[i] - 2 * (2 + power) * invr2 * iprod(Sr, r) * r[i]); @@ -763,11 +745,11 @@ void update_orires_history(const t_oriresdata& od, history_t* hist) * in calc_orires_dev. */ hist->orire_initf = od.exp_min_t_tau; - for (int pair = 0; pair < od.nr; pair++) + for (int pair = 0; pair < od.numRestraints; pair++) { for (int i = 0; i < 5; i++) { - hist->orire_Dtav[pair * 5 + i] = od.Dtav[pair][i]; + hist->orire_Dtav[pair * 5 + i] = od.DTensorsTimeAndEnsembleAv[pair][i]; } } } diff --git a/src/gromacs/listed_forces/orires.h b/src/gromacs/listed_forces/orires.h index 981cbed67e..3a298220d6 100644 --- a/src/gromacs/listed_forces/orires.h +++ b/src/gromacs/listed_forces/orires.h @@ -67,22 +67,6 @@ template class ArrayRef; } // namespace gmx -/*! \brief - * Decides whether orientation restraints can work, and initializes - * all the orientation restraint stuff in *od (and assumes *od is - * already allocated. - * If orientation restraint are used, globalState is read and modified - * on the master rank (which is the only rank, since orientation - * restraints can not run in parallel). - */ -void init_orires(FILE* fplog, - const gmx_mtop_t& mtop, - const t_inputrec* ir, - const t_commrec* cr, - const gmx_multisim_t* ms, - t_state* globalState, - t_oriresdata* od); - /*! \brief * Calculates the time averaged D matrices, the S matrix for each experiment. * diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index 374f12c13d..631e7b1dad 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -1191,18 +1191,20 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, nr[i] = 0; } - if (bOR && fcd->orires->nr > 0) + if (bOR && fcd->orires) { t_oriresdata& orires = *fcd->orires; diagonalize_orires_tensors(&orires); - nr[enxOR] = orires.nr; - block[enxOR] = orires.otav; - id[enxOR] = enxOR; - nr[enxORI] = (orires.oinsl != orires.otav) ? orires.nr : 0; - block[enxORI] = orires.oinsl; + nr[enxOR] = orires.numRestraints; + block[enxOR] = orires.orientationsTimeAndEnsembleAv.data(); + id[enxOR] = enxOR; + nr[enxORI] = (orires.orientations.data() != orires.orientationsTimeAndEnsembleAv.data()) + ? orires.numRestraints + : 0; + block[enxORI] = orires.orientations.data(); id[enxORI] = enxORI; - nr[enxORT] = orires.nex * 12; - block[enxORT] = orires.eig; + nr[enxORT] = ssize(orires.eigenOutput); + block[enxORT] = orires.eigenOutput.data(); id[enxORT] = enxORT; } @@ -1288,9 +1290,9 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene, free_enxframe(&fr); if (log) { - if (bOR && fcd->orires->nr > 0) + if (bOR && fcd->orires) { - print_orires_log(log, fcd->orires); + print_orires_log(log, fcd->orires.get()); } fprintf(log, " Energies (%s)\n", unit_energy); diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 789b892cba..186d257c48 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -283,8 +283,7 @@ void gmx::LegacySimulator::do_md() { // TODO This implementation of ensemble orientation restraints is nasty because // a user can't just do multi-sim with single-sim orientation restraints. - bool usingEnsembleRestraints = - (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0)); + bool usingEnsembleRestraints = (fcdata.disres->nsystems > 1) || ((ms != nullptr) && fcdata.orires); bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim() && (ms != nullptr)); // Replica exchange, ensemble restraints and AWH need all @@ -451,7 +450,7 @@ void gmx::LegacySimulator::do_md() "Essential dynamics is not supported with the GPU update.\n"); GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(*ir->pull), "Constraints pulling is not supported with the GPU update.\n"); - GMX_RELEASE_ASSERT(fcdata.orires->nr == 0, + GMX_RELEASE_ASSERT(fcdata.orires == nullptr, "Orientation restraints are not supported with the GPU update.\n"); GMX_RELEASE_ASSERT( ir->efep == FreeEnergyPerturbationType::No diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 40c538a60a..f309a376a2 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -1208,9 +1208,11 @@ int Mdrunner::mdrunner() globalState.get(), replExParams.exchangeInterval > 0); - t_oriresdata* oriresdata; - snew(oriresdata, 1); - init_orires(fplog, mtop, inputrec.get(), cr, ms, globalState.get(), oriresdata); + std::unique_ptr oriresData; + if (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0) + { + oriresData = std::make_unique(fplog, mtop, *inputrec, cr, ms, globalState.get()); + } auto deform = prepareBoxDeformation(globalState != nullptr ? globalState->box : box, MASTER(cr) ? DDRole::Master : DDRole::Agent, @@ -1661,7 +1663,7 @@ int Mdrunner::mdrunner() pforce); // Dirty hack, for fixing disres and orires should be made mdmodules fr->fcdata->disres = disresdata; - fr->fcdata->orires = oriresdata; + fr->fcdata->orires.swap(oriresData); // Save a handle to device stream manager to use elsewhere in the code // TODO: Forcerec is not a correct place to store it. @@ -2104,7 +2106,6 @@ int Mdrunner::mdrunner() fr.reset(nullptr); // destruct forcerec before gpu // TODO convert to C++ so we can get rid of these frees sfree(disresdata); - sfree(oriresdata); if (!hwinfo_->deviceInfoList.empty()) { diff --git a/src/gromacs/mdtypes/fcdata.h b/src/gromacs/mdtypes/fcdata.h index 2f36599bee..0da8e1a69e 100644 --- a/src/gromacs/mdtypes/fcdata.h +++ b/src/gromacs/mdtypes/fcdata.h @@ -38,13 +38,21 @@ #ifndef GMX_MDTYPES_FCDATA_H #define GMX_MDTYPES_FCDATA_H +#include #include #include "gromacs/math/vectypes.h" #include "gromacs/topology/idef.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/classhelpers.h" #include "gromacs/utility/real.h" enum class DistanceRestraintWeighting : int; +struct gmx_mtop_t; +struct gmx_multisim_t; +struct t_commrec; +struct t_inputrec; +class t_state; typedef real rvec5[5]; @@ -84,36 +92,92 @@ struct OriresMatEq }; /* Orientation restraining stuff */ -typedef struct t_oriresdata +struct t_oriresdata { - real fc; /* Force constant for the restraints */ - real edt; /* Multiplication factor for time averaging */ - real edt_1; /* 1 - edt */ - real exp_min_t_tau; /* Factor for slowly switching on the force */ - int nr; /* The number of orientation restraints */ - int nex; /* The number of experiments */ - int typeMin; /* The minimum iparam type index for restraints */ - int nref; /* The number of atoms for the fit */ - real* mref; /* The masses of the reference atoms */ - rvec* xref; /* The reference coordinates for the fit (nref) */ - rvec* xtmp; /* Temporary array for fitting (nref) */ - matrix R; /* Rotation matrix to rotate to the reference coor. */ - tensor* S; /* Array of order tensors for each experiment (nexp) */ - rvec5* Dinsl; /* The order matrix D for all restraints (nr x 5) */ - rvec5* Dins; /* The ensemble averaged D (nr x 5) */ - rvec5* Dtav; /* The time and ensemble averaged D (nr x 5) */ - real* oinsl; /* The calculated instantaneous orientations */ - real* oins; /* The calculated emsemble averaged orientations */ - real* otav; /* The calculated time and ensemble averaged orient. */ - real rmsdev; /* The weighted (using kfac) RMS deviation */ - OriresMatEq* tmpEq; /* An temporary array of matrix + rhs */ - real* eig; /* Eigenvalues/vectors, for output only (nex x 12) */ - - /* variables for diagonalization with diagonalize_orires_tensors()*/ - double** M; - double* eig_diag; - double** v; -} t_oriresdata; + /*! Constructor + * + * \param[in] fplog Log file, can be nullptr + * \param[in] mtop The global topology + * \param[in] ir The input record + * \param[in] cr The commrec, can be nullptr when not running in parallel + * \param[in] ms The multisim communicator, pass nullptr to avoid ensemble averaging + * \param[in,out] globalState The global state, orientation restraint entires are added + * + * \throws InvalidInputError when there is domain decomposition, fewer than 5 restraints, + * periodic molecules or more than 1 molecule for a moleculetype with restraints. + */ + t_oriresdata(FILE* fplog, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const t_commrec* cr, + const gmx_multisim_t* ms, + t_state* globalState); + + //! Destructor + ~t_oriresdata(); + + //! Force constant for the restraints + real fc; + //! Multiplication factor for time averaging + real edt; + //! 1 - edt + real edt_1; + //! Factor for slowly switching on the force + real exp_min_t_tau; + //! The number of orientation restraints + const int numRestraints; + //! The number of experiments + int numExperiments; + //! The minimum iparam type index for restraints + int typeMin; + //! The number of atoms for the fit + int numReferenceAtoms; + //! The masses of the reference atoms + std::vector mref; + //! The reference coordinates for the fit + std::vector xref; + //! Temporary array for fitting + std::vector xtmp; + //! Rotation matrix to rotate to the reference coordinates + matrix rotationMatrix; + //! Array of order tensors, one for each experiment + tensor* orderTensors = nullptr; + //! The order tensor D for all restraints + rvec5* DTensors = nullptr; + //! The ensemble averaged D for all restraints + rvec5* DTensorsEnsembleAv = nullptr; + //! The time and ensemble averaged D restraints + rvec5* DTensorsTimeAndEnsembleAv = nullptr; + //! The calculated instantaneous orientations + std::vector orientations; + //! The calculated emsemble averaged orientations + gmx::ArrayRef orientationsEnsembleAv; + //! Buffer for the calculated emsemble averaged orientations, only used with ensemble averaging + std::vector orientationsEnsembleAvBuffer; + //! The calculated time and ensemble averaged orientations + gmx::ArrayRef orientationsTimeAndEnsembleAv; + //! The weighted (using kfac) RMS deviation + std::vector orientationsTimeAndEnsembleAvBuffer; + //! Buffer for the weighted (using kfac) RMS deviation, only used with time averaging + real rmsdev; + //! An temporary array of matrix + rhs + std::vector tmpEq; + //! The number of eigenvalues + eigenvectors per experiment + static constexpr int c_numEigenRealsPerExperiment = 12; + //! Eigenvalues/vectors, for output only (numExperiments x 12) + std::vector eigenOutput; + + // variables for diagonalization with diagonalize_orires_tensors() + //! Tensor to diagonalize + std::array M; + //! Eigenvalues + std::array eig_diag; + //! Eigenvectors + std::array v; + + // Default copy and assign would be incorrect and manual versions are not yet implemented. + GMX_DISALLOW_COPY_AND_ASSIGN(t_oriresdata); +}; /* Cubic spline table for tabulated bonded interactions */ struct bondedtable_t @@ -137,8 +201,8 @@ struct t_fcdata std::vector dihtab; // TODO: Convert to C++ and unique_ptr (currently this data is not freed) - t_disresdata* disres = nullptr; - t_oriresdata* orires = nullptr; + t_disresdata* disres = nullptr; + std::unique_ptr orires; }; #endif -- 2.22.0