#include <stdio.h>
#include <string.h>
+#include <cstdint>
+
#include <array>
#include "gromacs/domdec/domdec.h"
#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/ewald/pme.h"
#include "gromacs/gmxlib/chargegroup.h"
-#include "gromacs/gmxlib/md_logging.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/timing/wallcyclereporting.h"
#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/sysinfo.h"
walltime_accounting_get_start_time_stamp(walltime_accounting));
}
-static void sum_forces(int start, int end, rvec f[], rvec flr[])
+static void sum_forces(rvec f[], const PaddedRVecVector *forceToAdd)
{
- int i;
+ /* TODO: remove this - 1 when padding is properly implemented */
+ int end = forceToAdd->size() - 1;
+ const rvec *fAdd = as_rvec_array(forceToAdd->data());
- if (gmx_debug_at)
- {
- pr_rvecs(debug, 0, "fsr", f+start, end-start);
- pr_rvecs(debug, 0, "flr", flr+start, end-start);
- }
// cppcheck-suppress unreadVariable
int gmx_unused nt = gmx_omp_nthreads_get(emntDefault);
#pragma omp parallel for num_threads(nt) schedule(static)
- for (i = start; i < end; i++)
+ for (int i = 0; i < end; i++)
{
- rvec_inc(f[i], flr[i]);
- }
-}
-
-/*
- * calc_f_el calculates forces due to an electric field.
- *
- * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
- *
- * Et[] contains the parameters for the time dependent
- * part of the field.
- * Ex[] contains the parameters for
- * the spatial dependent part of the field.
- * The function should return the energy due to the electric field
- * (if any) but for now returns 0.
- *
- * WARNING:
- * There can be problems with the virial.
- * Since the field is not self-consistent this is unavoidable.
- * For neutral molecules the virial is correct within this approximation.
- * For neutral systems with many charged molecules the error is small.
- * But for systems with a net charge or a few charged molecules
- * the error can be significant when the field is high.
- * Solution: implement a self-consistent electric field into PME.
- */
-static void calc_f_el(FILE *fp, int start, int homenr,
- real charge[], rvec f[],
- t_cosines Ex[], t_cosines Et[], double t)
-{
- rvec Ext;
- real t0;
- int i, m;
-
- for (m = 0; (m < DIM); m++)
- {
- if (Et[m].n > 0)
- {
- if (Et[m].n == 3)
- {
- t0 = Et[m].a[1];
- Ext[m] = std::cos(Et[m].a[0]*(t-t0))*std::exp(-gmx::square(t-t0)/(2.0*gmx::square(Et[m].a[2])));
- }
- else
- {
- Ext[m] = std::cos(Et[m].a[0]*t);
- }
- }
- else
- {
- Ext[m] = 1.0;
- }
- if (Ex[m].n > 0)
- {
- /* Convert the field strength from V/nm to MD-units */
- Ext[m] *= Ex[m].a[0]*FIELDFAC;
- for (i = start; (i < start+homenr); i++)
- {
- f[i][m] += charge[i]*Ext[m];
- }
- }
- else
- {
- Ext[m] = 0;
- }
- }
- if (fp != NULL)
- {
- fprintf(fp, "%10g %10g %10g %10g #FIELD\n", t,
- Ext[XX]/FIELDFAC, Ext[YY]/FIELDFAC, Ext[ZZ]/FIELDFAC);
+ rvec_inc(f[i], fAdd[i]);
}
}
wallcycle_start(wcycle, ewcPP_PMEWAITRECVF);
dvdl_q = 0;
dvdl_lj = 0;
- gmx_pme_receive_f(cr, fr->f_novirsum, fr->vir_el_recip, &e_q,
+ gmx_pme_receive_f(cr, as_rvec_array(fr->f_novirsum->data()), fr->vir_el_recip, &e_q,
fr->vir_lj_recip, &e_lj, &dvdl_q, &dvdl_lj,
&cycles_seppme);
enerd->term[F_COUL_RECIP] += e_q;
}
static void print_large_forces(FILE *fp, t_mdatoms *md, t_commrec *cr,
- gmx_int64_t step, real pforce, rvec *x, rvec *f)
+ gmx_int64_t step, real forceTolerance,
+ const rvec *x, const rvec *f)
{
- int i;
- real pf2, fn2;
- char buf[STEPSTRSIZE];
-
- pf2 = gmx::square(pforce);
- for (i = 0; i < md->homenr; i++)
+ real force2Tolerance = gmx::square(forceTolerance);
+ std::uintmax_t numNonFinite = 0;
+ for (int i = 0; i < md->homenr; i++)
{
- fn2 = norm2(f[i]);
- /* We also catch NAN, if the compiler does not optimize this away. */
- if (fn2 >= pf2 || fn2 != fn2)
+ real force2 = norm2(f[i]);
+ bool nonFinite = !std::isfinite(force2);
+ if (force2 >= force2Tolerance || nonFinite)
+ {
+ fprintf(fp, "step %" GMX_PRId64 " atom %6d x %8.3f %8.3f %8.3f force %12.5e\n",
+ step,
+ ddglatnr(cr->dd, i), x[i][XX], x[i][YY], x[i][ZZ], std::sqrt(force2));
+ }
+ if (nonFinite)
{
- fprintf(fp, "step %s atom %6d x %8.3f %8.3f %8.3f force %12.5e\n",
- gmx_step_str(step, buf),
- ddglatnr(cr->dd, i), x[i][XX], x[i][YY], x[i][ZZ], std::sqrt(fn2));
+ numNonFinite++;
}
}
+ if (numNonFinite > 0)
+ {
+ /* Note that with MPI this fatal call on one rank might interrupt
+ * the printing on other ranks. But we can only avoid that with
+ * an expensive MPI barrier that we would need at each step.
+ */
+ gmx_fatal(FARGS, "At step %" GMX_PRId64 " detected non-finite forces on %ju atoms", step, numNonFinite);
+ }
}
static void post_process_forces(t_commrec *cr,
* if the constructing atoms aren't local.
*/
wallcycle_start(wcycle, ewcVSITESPREAD);
- spread_vsite_f(vsite, x, fr->f_novirsum, NULL,
+ spread_vsite_f(vsite, x, as_rvec_array(fr->f_novirsum->data()), nullptr,
(flags & GMX_FORCE_VIRIAL), fr->vir_el_recip,
nrnb,
&top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
if (flags & GMX_FORCE_VIRIAL)
{
/* Now add the forces, this is local */
- if (fr->bDomDec)
- {
- sum_forces(0, fr->f_novirsum_n, f, fr->f_novirsum);
- }
- else
- {
- sum_forces(0, mdatoms->homenr,
- f, fr->f_novirsum);
- }
+ sum_forces(f, fr->f_novirsum);
+
if (EEL_FULL(fr->eeltype))
{
/* Add the mesh contribution to the virial */
gmx_bool use_GPU(const nonbonded_verlet_t *nbv)
{
- return nbv != NULL && nbv->bUseGPU;
+ return nbv != nullptr && nbv->bUseGPU;
}
static gmx_inline void clear_rvecs_omp(int n, rvec v[])
gmx_localtop_t *top,
gmx_groups_t gmx_unused *groups,
matrix box, rvec x[], history_t *hist,
- rvec f[],
+ PaddedRVecVector *force,
tensor vir_force,
t_mdatoms *mdatoms,
gmx_enerdata_t *enerd, t_fcdata *fcd,
real *lambda, t_graph *graph,
t_forcerec *fr, interaction_const_t *ic,
gmx_vsite_t *vsite, rvec mu_tot,
- double t, FILE *field, gmx_edsam_t ed,
+ double t, gmx_edsam_t ed,
gmx_bool bBornRadii,
int flags)
{
int cg1, i, j;
- int start, homenr;
double mu[2*DIM];
gmx_bool bStateChanged, bNS, bFillGrid, bCalcCGCM;
gmx_bool bDoForces, bUseGPU, bUseOrEmulGPU;
cycles_wait_gpu = 0;
nbv = fr->nbv;
- start = 0;
- homenr = mdatoms->homenr;
+ const int start = 0;
+ const int homenr = mdatoms->homenr;
clear_mat(vir_force);
* we do not need to worry about shifting.
*/
- int pme_flags = 0;
-
wallcycle_start(wcycle, ewcPP_PMESENDX);
bBS = (inputrec->nwall == 2);
svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
}
- if (EEL_PME(fr->eeltype))
- {
- pme_flags |= GMX_PME_DO_COULOMB;
- }
-
- if (EVDW_PME(fr->vdwtype))
- {
- pme_flags |= GMX_PME_DO_LJ;
- }
-
gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
- mdatoms->nChargePerturbed, mdatoms->nTypePerturbed, lambda[efptCOUL], lambda[efptVDW],
+ lambda[efptCOUL], lambda[efptVDW],
(flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
- pme_flags, step);
+ step);
wallcycle_stop(wcycle, ewcPP_PMESENDX);
}
nbnxn_put_on_grid(nbv->nbs, fr->ePBC, box,
0, vzero, box_diag,
0, mdatoms->homenr, -1, fr->cginfo, x,
- 0, NULL,
+ 0, nullptr,
nbv->grp[eintLocal].kernel_type,
nbv->grp[eintLocal].nbat);
wallcycle_sub_stop(wcycle, ewcsNBS_GRID_LOCAL);
wallcycle_stop(wcycle, ewcROT);
}
+ /* Temporary solution until all routines take PaddedRVecVector */
+ rvec *f = as_rvec_array(force->data());
+
/* Start the force cycle counter.
* This counter is stopped after do_force_lowlevel.
* No parallel communication should occur while this counter is running,
{
if (flags & GMX_FORCE_VIRIAL)
{
- fr->f_novirsum = fr->f_novirsum_alloc;
+ fr->f_novirsum = fr->forceBufferNoVirialSummation;
}
else
{
* a separate array for forces that do not contribute
* to the pressure.
*/
- fr->f_novirsum = f;
+ fr->f_novirsum = force;
}
}
{
if (flags & GMX_FORCE_VIRIAL)
{
- if (fr->bDomDec)
- {
- clear_rvecs_omp(fr->f_novirsum_n, fr->f_novirsum);
- }
- else
- {
- clear_rvecs_omp(homenr, fr->f_novirsum+start);
- }
+ /* TODO: remove this - 1 when padding is properly implemented */
+ clear_rvecs_omp(fr->f_novirsum->size() - 1,
+ as_rvec_array(fr->f_novirsum->data()));
}
}
/* Clear the short- and long-range forces */
if (bDoForces)
{
- if (inputrecElecField(inputrec))
+ /* Compute forces due to electric field */
+ if (fr->efield != nullptr)
{
- /* Compute forces due to electric field */
- calc_f_el(MASTER(cr) ? field : NULL,
- start, homenr, mdatoms->chargeA, fr->f_novirsum,
- inputrec->ex, inputrec->et, t);
+ fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
}
/* If we have NoVirSum forces, but we do not calculate the virial,
if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
{
wallcycle_start(wcycle, ewcVSITESPREAD);
- spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
+ spread_vsite_f(vsite, x, f, fr->fshift, FALSE, nullptr, nrnb,
&top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
wallcycle_stop(wcycle, ewcVSITESPREAD);
}
gmx_localtop_t *top,
gmx_groups_t *groups,
matrix box, rvec x[], history_t *hist,
- rvec f[],
+ PaddedRVecVector *force,
tensor vir_force,
t_mdatoms *mdatoms,
gmx_enerdata_t *enerd, t_fcdata *fcd,
real *lambda, t_graph *graph,
t_forcerec *fr, gmx_vsite_t *vsite, rvec mu_tot,
- double t, FILE *field, gmx_edsam_t ed,
+ double t, gmx_edsam_t ed,
gmx_bool bBornRadii,
int flags)
{
int cg0, cg1, i, j;
- int start, homenr;
double mu[2*DIM];
gmx_bool bStateChanged, bNS, bFillGrid, bCalcCGCM;
gmx_bool bDoForces;
float cycles_pme, cycles_force;
- start = 0;
- homenr = mdatoms->homenr;
+ const int start = 0;
+ const int homenr = mdatoms->homenr;
clear_mat(vir_force);
* we do not need to worry about shifting.
*/
- int pme_flags = 0;
-
wallcycle_start(wcycle, ewcPP_PMESENDX);
bBS = (inputrec->nwall == 2);
svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
}
- if (EEL_PME(fr->eeltype))
- {
- pme_flags |= GMX_PME_DO_COULOMB;
- }
-
- if (EVDW_PME(fr->vdwtype))
- {
- pme_flags |= GMX_PME_DO_LJ;
- }
-
gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
- mdatoms->nChargePerturbed, mdatoms->nTypePerturbed, lambda[efptCOUL], lambda[efptVDW],
+ lambda[efptCOUL], lambda[efptVDW],
(flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
- pme_flags, step);
+ step);
wallcycle_stop(wcycle, ewcPP_PMESENDX);
}
wallcycle_stop(wcycle, ewcROT);
}
+ /* Temporary solution until all routines take PaddedRVecVector */
+ rvec *f = as_rvec_array(force->data());
+
/* Start the force cycle counter.
* This counter is stopped after do_force_lowlevel.
* No parallel communication should occur while this counter is running,
{
if (flags & GMX_FORCE_VIRIAL)
{
- fr->f_novirsum = fr->f_novirsum_alloc;
- if (fr->bDomDec)
- {
- clear_rvecs(fr->f_novirsum_n, fr->f_novirsum);
- }
- else
- {
- clear_rvecs(homenr, fr->f_novirsum+start);
- }
+ fr->f_novirsum = fr->forceBufferNoVirialSummation;
+ /* TODO: remove this - 1 when padding is properly implemented */
+ clear_rvecs(fr->f_novirsum->size() - 1,
+ as_rvec_array(fr->f_novirsum->data()));
}
else
{
* a separate array for forces that do not contribute
* to the pressure.
*/
- fr->f_novirsum = f;
+ fr->f_novirsum = force;
}
}
if (bDoForces)
{
- if (inputrecElecField(inputrec))
+ /* Compute forces due to electric field */
+ if (fr->efield != nullptr)
{
- /* Compute forces due to electric field */
- calc_f_el(MASTER(cr) ? field : NULL,
- start, homenr, mdatoms->chargeA, fr->f_novirsum,
- inputrec->ex, inputrec->et, t);
+ fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
}
/* Communicate the forces */
if (EEL_FULL(fr->eeltype) && cr->dd->n_intercg_excl &&
(flags & GMX_FORCE_VIRIAL))
{
- dd_move_f(cr->dd, fr->f_novirsum, NULL);
+ dd_move_f(cr->dd, as_rvec_array(fr->f_novirsum->data()), nullptr);
}
wallcycle_stop(wcycle, ewcMOVEF);
}
if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
{
wallcycle_start(wcycle, ewcVSITESPREAD);
- spread_vsite_f(vsite, x, f, fr->fshift, FALSE, NULL, nrnb,
+ spread_vsite_f(vsite, x, f, fr->fshift, FALSE, nullptr, nrnb,
&top->idef, fr->ePBC, fr->bMolPBC, graph, box, cr);
wallcycle_stop(wcycle, ewcVSITESPREAD);
}
gmx_int64_t step, t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_localtop_t *top,
gmx_groups_t *groups,
- matrix box, rvec x[], history_t *hist,
- rvec f[],
+ matrix box, PaddedRVecVector *coordinates, history_t *hist,
+ PaddedRVecVector *force,
tensor vir_force,
t_mdatoms *mdatoms,
gmx_enerdata_t *enerd, t_fcdata *fcd,
- real *lambda, t_graph *graph,
+ gmx::ArrayRef<real> lambda, t_graph *graph,
t_forcerec *fr,
gmx_vsite_t *vsite, rvec mu_tot,
- double t, FILE *field, gmx_edsam_t ed,
+ double t, gmx_edsam_t ed,
gmx_bool bBornRadii,
int flags)
{
flags &= ~GMX_FORCE_NONBONDED;
}
+ GMX_ASSERT(coordinates->size() >= static_cast<unsigned int>(fr->natoms_force + 1), "We might need 1 element extra for SIMD");
+ GMX_ASSERT(force->size() >= static_cast<unsigned int>(fr->natoms_force + 1), "We might need 1 element extra for SIMD");
+
+ rvec *x = as_rvec_array(coordinates->data());
+
switch (inputrec->cutoff_scheme)
{
case ecutsVERLET:
top,
groups,
box, x, hist,
- f, vir_force,
+ force, vir_force,
mdatoms,
enerd, fcd,
- lambda, graph,
+ lambda.data(), graph,
fr, fr->ic,
vsite, mu_tot,
- t, field, ed,
+ t, ed,
bBornRadii,
flags);
break;
top,
groups,
box, x, hist,
- f, vir_force,
+ force, vir_force,
mdatoms,
enerd, fcd,
- lambda, graph,
+ lambda.data(), graph,
fr, vsite, mu_tot,
- t, field, ed,
+ t, ed,
bBornRadii,
flags);
break;
dvdl_dum = 0;
/* constrain the current position */
- constrain(NULL, TRUE, FALSE, constr, &(top->idef),
+ constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
ir, cr, step, 0, 1.0, md,
- state->x, state->x, NULL,
+ as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), nullptr,
fr->bMolPBC, state->box,
state->lambda[efptBONDED], &dvdl_dum,
- NULL, NULL, nrnb, econqCoord);
+ nullptr, nullptr, nrnb, econqCoord);
if (EI_VV(ir->eI))
{
/* constrain the inital velocity, and save it */
/* also may be useful if we need the ekin from the halfstep for velocity verlet */
- constrain(NULL, TRUE, FALSE, constr, &(top->idef),
+ constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
ir, cr, step, 0, 1.0, md,
- state->x, state->v, state->v,
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
fr->bMolPBC, state->box,
state->lambda[efptBONDED], &dvdl_dum,
- NULL, NULL, nrnb, econqVeloc);
+ nullptr, nullptr, nrnb, econqVeloc);
}
/* constrain the inital velocities at t-dt/2 */
if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
gmx_step_str(step, buf));
}
dvdl_dum = 0;
- constrain(NULL, TRUE, FALSE, constr, &(top->idef),
+ constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
ir, cr, step, -1, 1.0, md,
- state->x, savex, NULL,
+ as_rvec_array(state->x.data()), savex, nullptr,
fr->bMolPBC, state->box,
state->lambda[efptBONDED], &dvdl_dum,
- state->v, NULL, nrnb, econqCoord);
+ as_rvec_array(state->v.data()), nullptr, nrnb, econqCoord);
for (i = start; i < end; i++)
{
else
{
/* Pass NULL iso fplog to avoid graph prints for each molecule type */
- mk_graph_ilist(NULL, mtop->moltype[molb->type].ilist,
+ mk_graph_ilist(nullptr, mtop->moltype[molb->type].ilist,
0, molb->natoms_mol, FALSE, FALSE, graph);
for (mol = 0; mol < molb->nmol; mol++)
}
// TODO This can be cleaned up a lot, and move back to runner.cpp
-void finish_run(FILE *fplog, t_commrec *cr,
+void finish_run(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
t_inputrec *inputrec,
t_nrnb nrnb[], gmx_wallcycle_t wcycle,
gmx_walltime_accounting_t walltime_accounting,
nonbonded_verlet_t *nbv,
gmx_bool bWriteStat)
{
- t_nrnb *nrnb_tot = NULL;
+ t_nrnb *nrnb_tot = nullptr;
double delta_t = 0;
double nbfs = 0, mflop = 0;
double elapsed_time,
if (!walltime_accounting_get_valid_finish(walltime_accounting))
{
- md_print_warn(cr, fplog,
- "Simulation ended prematurely, no performance report will be written.");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("Simulation ended prematurely, no performance report will be written.");
return;
}
if (SIMMASTER(cr))
{
- struct gmx_wallclock_gpu_t* gputimes = use_GPU(nbv) ? nbnxn_gpu_get_timings(nbv->gpu_nbv) : NULL;
+ struct gmx_wallclock_gpu_t* gputimes = use_GPU(nbv) ? nbnxn_gpu_get_timings(nbv->gpu_nbv) : nullptr;
- wallcycle_print(fplog, cr->nnodes, cr->npmenodes, nthreads_pp, nthreads_pme,
+ wallcycle_print(fplog, mdlog, cr->nnodes, cr->npmenodes, nthreads_pp, nthreads_pme,
elapsed_time_over_all_ranks,
wcycle, cycle_sum, gputimes);
}
}
-extern void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, real *lambda, double *lam0)
+extern void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, gmx::ArrayRef<real> lambda, double *lam0)
{
/* this function works, but could probably use a logic rewrite to keep all the different
types of efep straight. */
- int i;
+ if ((ir->efep == efepNO) && (ir->bSimTemp == FALSE))
+ {
+ return;
+ }
+
t_lambda *fep = ir->fepvals;
+ *fep_state = fep->init_fep_state; /* this might overwrite the checkpoint
+ if checkpoint is set -- a kludge is in for now
+ to prevent this.*/
- if ((ir->efep == efepNO) && (ir->bSimTemp == FALSE))
+ for (int i = 0; i < efptNR; i++)
{
- for (i = 0; i < efptNR; i++)
+ /* overwrite lambda state with init_lambda for now for backwards compatibility */
+ if (fep->init_lambda >= 0) /* if it's -1, it was never initializd */
{
- lambda[i] = 0.0;
+ lambda[i] = fep->init_lambda;
if (lam0)
{
- lam0[i] = 0.0;
+ lam0[i] = lambda[i];
}
}
- return;
- }
- else
- {
- *fep_state = fep->init_fep_state; /* this might overwrite the checkpoint
- if checkpoint is set -- a kludge is in for now
- to prevent this.*/
- for (i = 0; i < efptNR; i++)
+ else
{
- /* overwrite lambda state with init_lambda for now for backwards compatibility */
- if (fep->init_lambda >= 0) /* if it's -1, it was never initializd */
- {
- lambda[i] = fep->init_lambda;
- if (lam0)
- {
- lam0[i] = lambda[i];
- }
- }
- else
+ lambda[i] = fep->all_lambda[i][*fep_state];
+ if (lam0)
{
- lambda[i] = fep->all_lambda[i][*fep_state];
- if (lam0)
- {
- lam0[i] = lambda[i];
- }
+ lam0[i] = lambda[i];
}
}
- if (ir->bSimTemp)
+ }
+ if (ir->bSimTemp)
+ {
+ /* need to rescale control temperatures to match current state */
+ for (int i = 0; i < ir->opts.ngtc; i++)
{
- /* need to rescale control temperatures to match current state */
- for (i = 0; i < ir->opts.ngtc; i++)
+ if (ir->opts.ref_t[i] > 0)
{
- if (ir->opts.ref_t[i] > 0)
- {
- ir->opts.ref_t[i] = ir->simtempvals->temperatures[*fep_state];
- }
+ ir->opts.ref_t[i] = ir->simtempvals->temperatures[*fep_state];
}
}
}
/* Send to the log the information on the current lambdas */
- if (fplog != NULL)
+ if (fplog != nullptr)
{
fprintf(fplog, "Initial vector of lambda components:[ ");
- for (i = 0; i < efptNR; i++)
+ for (const auto &l : lambda)
{
- fprintf(fplog, "%10.4f ", lambda[i]);
+ fprintf(fplog, "%10.4f ", l);
}
fprintf(fplog, "]\n");
}
void init_md(FILE *fplog,
t_commrec *cr, t_inputrec *ir, const gmx_output_env_t *oenv,
double *t, double *t0,
- real *lambda, int *fep_state, double *lam0,
+ gmx::ArrayRef<real> lambda, int *fep_state, double *lam0,
t_nrnb *nrnb, gmx_mtop_t *mtop,
gmx_update_t **upd,
int nfile, const t_filenm fnm[],
}
if (*bSimAnn)
{
- update_annealing_target_temp(ir, ir->init_t, upd ? *upd : NULL);
+ update_annealing_target_temp(ir, ir->init_t, upd ? *upd : nullptr);
}
- if (vcm != NULL)
+ if (vcm != nullptr)
{
*vcm = init_vcm(fplog, &mtop->groups, ir);
}
please_cite(fplog, "Goga2012");
}
}
- if ((ir->et[XX].n > 0) || (ir->et[YY].n > 0) || (ir->et[ZZ].n > 0))
- {
- please_cite(fplog, "Caleman2008a");
- }
init_nrnb(nrnb);
if (nfile != -1)
{
*outf = init_mdoutf(fplog, nfile, fnm, Flags, cr, ir, mtop, oenv, wcycle);
- *mdebin = init_mdebin((Flags & MD_APPENDFILES) ? NULL : mdoutf_get_fp_ene(*outf),
+ *mdebin = init_mdebin((Flags & MD_APPENDFILES) ? nullptr : mdoutf_get_fp_ene(*outf),
mtop, ir, mdoutf_get_fp_dhdl(*outf));
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff