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;
* 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()), NULL,
(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_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);
}
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,
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()), NULL);
}
wallcycle_stop(wcycle, ewcMOVEF);
}
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,
+ std::vector<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;
/* constrain the current position */
constrain(NULL, 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()), NULL,
fr->bMolPBC, state->box,
state->lambda[efptBONDED], &dvdl_dum,
NULL, NULL, nrnb, econqCoord);
/* also may be useful if we need the ekin from the halfstep for velocity verlet */
constrain(NULL, 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);
dvdl_dum = 0;
constrain(NULL, TRUE, FALSE, constr, &(top->idef),
ir, cr, step, -1, 1.0, md,
- state->x, savex, NULL,
+ as_rvec_array(state->x.data()), savex, NULL,
fr->bMolPBC, state->box,
state->lambda[efptBONDED], &dvdl_dum,
- state->v, NULL, nrnb, econqCoord);
+ as_rvec_array(state->v.data()), NULL, nrnb, econqCoord);
for (i = start; i < end; i++)
{
}
// 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,
{
struct gmx_wallclock_gpu_t* gputimes = use_GPU(nbv) ? nbnxn_gpu_get_timings(nbv->gpu_nbv) : NULL;
- 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, std::vector<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))
+ lambda->resize(efptNR);
+
+ 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];
}
}
}
if (fplog != NULL)
{
fprintf(fplog, "Initial vector of lambda components:[ ");
- for (i = 0; i < efptNR; i++)
+ for (int i = 0; i < efptNR; i++)
{
- fprintf(fplog, "%10.4f ", lambda[i]);
+ fprintf(fplog, "%10.4f ", (*lambda)[i]);
}
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,
+ std::vector<real> *lambda, int *fep_state, double *lam0,
t_nrnb *nrnb, gmx_mtop_t *mtop,
gmx_update_t **upd,
int nfile, const t_filenm fnm[],
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)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff