Changed x, v and cg_p pointers in t_state to PaddedRVecVector.
PaddedRVecVector is currently std::vector<RVec> and is always
sized with one element extra for 4-wide SIMD (gather) loads.
But this will be replaced by a proper class.
t_state is now freed by the default destructor with some issues:
* ekinstate_t is fully converted to std::vector in a child change
* state struct of special algorithms still need conversion
Changes to t_state propagate over large parts of the core code.
This change extracts rvec pointers at intermediate levels to
make this single change manageable.
The state in EM also needed to be updated, which led to a lot
of changes in minimize.cpp, especially to do_lbfgs() which was not
yet converted to use em_state_t.
Change-Id: Ib6ead1f13741ead8dbeeddd118ecdd4bf1cc6584
gmx_incons("The state does not the domain decomposition state");
}
- state->ncg_gl = dd->ncg_home;
- if (state->ncg_gl > state->cg_gl_nalloc)
- {
- state->cg_gl_nalloc = over_alloc_dd(state->ncg_gl);
- srenew(state->cg_gl, state->cg_gl_nalloc);
- }
- for (i = 0; i < state->ncg_gl; i++)
+ state->cg_gl.resize(dd->ncg_home);
+ for (i = 0; i < dd->ncg_home; i++)
{
state->cg_gl[i] = dd->index_gl[i];
}
cgs_gl = &dd->comm->cgs_gl;
- ncg_home = state_local->ncg_gl;
- cg = state_local->cg_gl;
+ ncg_home = state_local->cg_gl.size();
+ cg = state_local->cg_gl.data();
nat_home = 0;
for (i = 0; i < ncg_home; i++)
{
}
static void dd_collect_vec_sendrecv(gmx_domdec_t *dd,
- rvec *lv, rvec *v)
+ const rvec *lv, rvec *v)
{
gmx_domdec_master_t *ma;
int n, i, c, a, nalloc = 0;
if (!DDMASTER(dd))
{
#if GMX_MPI
- MPI_Send(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, DDMASTERRANK(dd),
- dd->rank, dd->mpi_comm_all);
+ MPI_Send(const_cast<void *>(static_cast<const void *>(lv)), dd->nat_home*sizeof(rvec), MPI_BYTE,
+ DDMASTERRANK(dd), dd->rank, dd->mpi_comm_all);
#endif
}
else
}
static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
- rvec *lv, rvec *v)
+ const rvec *lv, rvec *v)
{
gmx_domdec_master_t *ma;
int *rcounts = NULL, *disps = NULL;
}
}
-void dd_collect_vec(gmx_domdec_t *dd,
- t_state *state_local, rvec *lv, rvec *v)
+void dd_collect_vec(gmx_domdec_t *dd,
+ t_state *state_local,
+ const PaddedRVecVector *localVector,
+ rvec *v)
{
dd_collect_cg(dd, state_local);
+ const rvec *lv = as_rvec_array(localVector->data());
+
if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
{
dd_collect_vec_sendrecv(dd, lv, v);
}
}
+void dd_collect_vec(gmx_domdec_t *dd,
+ t_state *state_local,
+ const PaddedRVecVector *localVector,
+ PaddedRVecVector *vector)
+{
+ dd_collect_vec(dd, state_local, localVector, as_rvec_array(vector->data()));
+}
+
void dd_collect_state(gmx_domdec_t *dd,
t_state *state_local, t_state *state)
switch (est)
{
case estX:
- dd_collect_vec(dd, state_local, state_local->x, state->x);
+ dd_collect_vec(dd, state_local, &state_local->x, &state->x);
break;
case estV:
- dd_collect_vec(dd, state_local, state_local->v, state->v);
+ dd_collect_vec(dd, state_local, &state_local->v, &state->v);
break;
case est_SDX_NOTSUPPORTED:
break;
case estCGP:
- dd_collect_vec(dd, state_local, state_local->cg_p, state->cg_p);
+ dd_collect_vec(dd, state_local, &state_local->cg_p, &state->cg_p);
break;
case estDISRE_INITF:
case estDISRE_RM3TAV:
}
}
-static void dd_realloc_state(t_state *state, rvec **f, int nalloc)
+static void dd_resize_state(t_state *state, PaddedRVecVector *f, int natoms)
{
int est;
if (debug)
{
- fprintf(debug, "Reallocating state: currently %d, required %d, allocating %d\n", state->nalloc, nalloc, over_alloc_dd(nalloc));
+ fprintf(debug, "Resizing state: currently %d, required %d\n", state->natoms, natoms);
}
- state->nalloc = over_alloc_dd(nalloc);
-
for (est = 0; est < estNR; est++)
{
if (EST_DISTR(est) && (state->flags & (1<<est)))
switch (est)
{
case estX:
- srenew(state->x, state->nalloc + 1);
+ state->x.resize(natoms + 1);
break;
case estV:
- srenew(state->v, state->nalloc + 1);
+ state->v.resize(natoms + 1);
break;
case est_SDX_NOTSUPPORTED:
break;
case estCGP:
- srenew(state->cg_p, state->nalloc + 1);
+ state->cg_p.resize(natoms + 1);
break;
case estDISRE_INITF:
case estDISRE_RM3TAV:
/* No reallocation required */
break;
default:
- gmx_incons("Unknown state entry encountered in dd_realloc_state");
+ gmx_incons("Unknown state entry encountered in dd_resize_state");
}
}
}
if (f != NULL)
{
- srenew(*f, state->nalloc);
+ (*f).resize(natoms + 1);
}
}
-static void dd_check_alloc_ncg(t_forcerec *fr, t_state *state, rvec **f,
- int nalloc)
+static void dd_check_alloc_ncg(t_forcerec *fr,
+ t_state *state,
+ PaddedRVecVector *f,
+ int numChargeGroups)
{
- if (nalloc > fr->cg_nalloc)
+ if (numChargeGroups > fr->cg_nalloc)
{
if (debug)
{
- fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, nalloc, over_alloc_dd(nalloc));
+ fprintf(debug, "Reallocating forcerec: currently %d, required %d, allocating %d\n", fr->cg_nalloc, numChargeGroups, over_alloc_dd(numChargeGroups));
}
- fr->cg_nalloc = over_alloc_dd(nalloc);
+ fr->cg_nalloc = over_alloc_dd(numChargeGroups);
srenew(fr->cginfo, fr->cg_nalloc);
if (fr->cutoff_scheme == ecutsGROUP)
{
srenew(fr->cg_cm, fr->cg_nalloc);
}
}
- if (fr->cutoff_scheme == ecutsVERLET && nalloc > state->nalloc)
+ if (fr->cutoff_scheme == ecutsVERLET)
{
/* We don't use charge groups, we use x in state to set up
* the atom communication.
*/
- dd_realloc_state(state, f, nalloc);
+ dd_resize_state(state, f, numChargeGroups);
}
}
static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
t_state *state, t_state *state_local,
- rvec **f)
+ PaddedRVecVector *f)
{
int i, j, nh;
}
}
}
- dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda);
+ dd_bcast(dd, ((efptNR)*sizeof(real)), state_local->lambda.data());
dd_bcast(dd, sizeof(int), &state_local->fep_state);
dd_bcast(dd, sizeof(real), &state_local->veta);
dd_bcast(dd, sizeof(real), &state_local->vol0);
dd_bcast(dd, sizeof(state_local->boxv), state_local->boxv);
dd_bcast(dd, sizeof(state_local->svir_prev), state_local->svir_prev);
dd_bcast(dd, sizeof(state_local->fvir_prev), state_local->fvir_prev);
- dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi);
- dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi);
- dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral);
- dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
- dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
+ dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_xi.data());
+ dd_bcast(dd, ((state_local->ngtc*nh)*sizeof(double)), state_local->nosehoover_vxi.data());
+ dd_bcast(dd, state_local->ngtc*sizeof(double), state_local->therm_integral.data());
+ dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi.data());
+ dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi.data());
/* communicate df_history -- required for restarting from checkpoint */
dd_distribute_dfhist(dd, state_local->dfhist);
- if (dd->nat_home > state_local->nalloc)
- {
- dd_realloc_state(state_local, f, dd->nat_home);
- }
+ dd_resize_state(state_local, f, dd->nat_home);
+
for (i = 0; i < estNR; i++)
{
if (EST_DISTR(i) && (state_local->flags & (1<<i)))
switch (i)
{
case estX:
- dd_distribute_vec(dd, cgs, state->x, state_local->x);
+ dd_distribute_vec(dd, cgs, as_rvec_array(state->x.data()), as_rvec_array(state_local->x.data()));
break;
case estV:
- dd_distribute_vec(dd, cgs, state->v, state_local->v);
+ dd_distribute_vec(dd, cgs, as_rvec_array(state->v.data()), as_rvec_array(state_local->v.data()));
break;
case est_SDX_NOTSUPPORTED:
break;
case estCGP:
- dd_distribute_vec(dd, cgs, state->cg_p, state_local->cg_p);
+ dd_distribute_vec(dd, cgs, as_rvec_array(state->cg_p.data()), as_rvec_array(state_local->cg_p.data()));
break;
case estDISRE_INITF:
case estDISRE_RM3TAV:
}
static void rebuild_cgindex(gmx_domdec_t *dd,
- const int *gcgs_index, t_state *state)
+ const int *gcgs_index, const t_state *state)
{
- int nat, i, *ind, *dd_cg_gl, *cgindex, cg_gl;
+ int * gmx_restrict dd_cg_gl = dd->index_gl;
+ int * gmx_restrict cgindex = dd->cgindex;
+ int nat = 0;
- ind = state->cg_gl;
- dd_cg_gl = dd->index_gl;
- cgindex = dd->cgindex;
- nat = 0;
+ /* Copy back the global charge group indices from state
+ * and rebuild the local charge group to atom index.
+ */
cgindex[0] = nat;
- for (i = 0; i < state->ncg_gl; i++)
+ for (unsigned int i = 0; i < state->cg_gl.size(); i++)
{
cgindex[i] = nat;
- cg_gl = ind[i];
+ int cg_gl = state->cg_gl[i];
dd_cg_gl[i] = cg_gl;
nat += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
}
- cgindex[i] = nat;
+ cgindex[state->cg_gl.size()] = nat;
- dd->ncg_home = state->ncg_gl;
+ dd->ncg_home = state->cg_gl.size();
dd->nat_home = nat;
set_zones_ncg_home(dd);
if (pos_d >= limit1[d])
{
cg_move_error(fplog, dd, step, cg, d, 1,
- cg_cm != state->x, limitd[d],
+ cg_cm != as_rvec_array(state->x.data()), limitd[d],
cg_cm[cg], cm_new, pos_d);
}
dev[d] = 1;
if (pos_d < limit0[d])
{
cg_move_error(fplog, dd, step, cg, d, -1,
- cg_cm != state->x, limitd[d],
+ cg_cm != as_rvec_array(state->x.data()), limitd[d],
cg_cm[cg], cm_new, pos_d);
}
dev[d] = -1;
static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
gmx_domdec_t *dd, ivec tric_dir,
- t_state *state, rvec **f,
+ t_state *state, PaddedRVecVector *f,
t_forcerec *fr,
gmx_bool bCompact,
t_nrnb *nrnb,
cgindex,
( thread *dd->ncg_home)/nthread,
((thread+1)*dd->ncg_home)/nthread,
- fr->cutoff_scheme == ecutsGROUP ? cg_cm : state->x,
+ fr->cutoff_scheme == ecutsGROUP ? cg_cm : as_rvec_array(state->x.data()),
move);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
*/
home_pos_cg =
compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
- nvec, state->x, comm, FALSE);
+ nvec, as_rvec_array(state->x.data()), comm, FALSE);
if (bCompact)
{
home_pos_cg -= *ncg_moved;
vec = 0;
home_pos_at =
compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
- nvec, vec++, state->x, comm, bCompact);
+ nvec, vec++, as_rvec_array(state->x.data()),
+ comm, bCompact);
if (bV)
{
compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
- nvec, vec++, state->v, comm, bCompact);
+ nvec, vec++, as_rvec_array(state->v.data()),
+ comm, bCompact);
}
if (bCGP)
{
compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
- nvec, vec++, state->cg_p, comm, bCompact);
+ nvec, vec++, as_rvec_array(state->cg_p.data()),
+ comm, bCompact);
}
if (bCompact)
nvr += i;
}
+ dd_check_alloc_ncg(fr, state, f, home_pos_cg + ncg_recv);
+ if (fr->cutoff_scheme == ecutsGROUP)
+ {
+ /* Here we resize to more than necessary and shrink later */
+ dd_resize_state(state, f, home_pos_at + ncg_recv*MAX_CGCGSIZE);
+ }
+
/* Process the received charge groups */
buf_pos = 0;
for (cg = 0; cg < ncg_recv; cg++)
dd->index_gl[home_pos_cg] = comm->buf_int[cg*DD_CGIBS];
dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
/* Copy the state from the buffer */
- dd_check_alloc_ncg(fr, state, f, home_pos_cg+1);
if (fr->cutoff_scheme == ecutsGROUP)
{
cg_cm = fr->cg_cm;
comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
}
- if (home_pos_at+nrcg > state->nalloc)
- {
- dd_realloc_state(state, f, home_pos_at+nrcg);
- }
for (i = 0; i < nrcg; i++)
{
copy_rvec(comm->vbuf.v[buf_pos++],
dd->ncg_home = home_pos_cg;
dd->nat_home = home_pos_at;
+ if (fr->cutoff_scheme == ecutsGROUP && !bCompact)
+ {
+ /* We overallocated before, we need to set the right size here */
+ dd_resize_state(state, f, dd->nat_home);
+ }
+
if (debug)
{
fprintf(debug,
dd = cr->dd;
set_ddbox(dd, FALSE, cr, ir, state->box,
- TRUE, &dd->comm->cgs_gl, state->x, &ddbox);
+ TRUE, &dd->comm->cgs_gl, as_rvec_array(state->x.data()), &ddbox);
LocallyLimited = 0;
static void setup_dd_communication(gmx_domdec_t *dd,
matrix box, gmx_ddbox_t *ddbox,
- t_forcerec *fr, t_state *state, rvec **f)
+ t_forcerec *fr,
+ t_state *state, PaddedRVecVector *f)
{
int dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
int nzone, nzone_send, zone, zonei, cg0, cg1;
cg_cm = fr->cg_cm;
break;
case ecutsVERLET:
- cg_cm = state->x;
+ cg_cm = as_rvec_array(state->x.data());
break;
default:
gmx_incons("unimplemented");
}
else
{
- cg_cm = state->x;
+ cg_cm = as_rvec_array(state->x.data());
}
/* Communicate cg_cm */
if (cd->bInPlace)
switch (i)
{
case estX:
- order_vec_atom(dd->ncg_home, cgindex, cgsort, state->x, vbuf);
+ order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->x.data()), vbuf);
break;
case estV:
- order_vec_atom(dd->ncg_home, cgindex, cgsort, state->v, vbuf);
+ order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->v.data()), vbuf);
break;
case est_SDX_NOTSUPPORTED:
break;
case estCGP:
- order_vec_atom(dd->ncg_home, cgindex, cgsort, state->cg_p, vbuf);
+ order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->cg_p.data()), vbuf);
break;
case estLD_RNG:
case estLD_RNGI:
const gmx_mtop_t *top_global,
const t_inputrec *ir,
t_state *state_local,
- rvec **f,
+ PaddedRVecVector *f,
t_mdatoms *mdatoms,
gmx_localtop_t *top_local,
t_forcerec *fr,
ncgindex_set = 0;
set_ddbox(dd, bMasterState, cr, ir, state_global->box,
- TRUE, cgs_gl, state_global->x, &ddbox);
+ TRUE, cgs_gl, as_rvec_array(state_global->x.data()), &ddbox);
get_cg_distribution(fplog, dd, cgs_gl,
- state_global->box, &ddbox, state_global->x);
+ state_global->box, &ddbox, as_rvec_array(state_global->x.data()));
dd_distribute_state(dd, cgs_gl,
state_global, state_local, f);
if (fr->cutoff_scheme == ecutsGROUP)
{
calc_cgcm(fplog, 0, dd->ncg_home,
- &top_local->cgs, state_local->x, fr->cg_cm);
+ &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
}
inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
{
/* Redetermine the cg COMs */
calc_cgcm(fplog, 0, dd->ncg_home,
- &top_local->cgs, state_local->x, fr->cg_cm);
+ &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
}
inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
set_ddbox(dd, bMasterState, cr, ir, state_local->box,
- TRUE, &top_local->cgs, state_local->x, &ddbox);
+ TRUE, &top_local->cgs, as_rvec_array(state_local->x.data()), &ddbox);
bRedist = dlbIsOn(comm);
}
copy_rvec(comm->box_size, ddbox.box_size);
}
set_ddbox(dd, bMasterState, cr, ir, state_local->box,
- bNStGlobalComm, &top_local->cgs, state_local->x, &ddbox);
+ bNStGlobalComm, &top_local->cgs, as_rvec_array(state_local->x.data()), &ddbox);
bBoxChanged = TRUE;
bRedist = TRUE;
0, dd->ncg_home,
comm->zones.dens_zone0,
fr->cginfo,
- state_local->x,
+ as_rvec_array(state_local->x.data()),
ncg_moved, bRedist ? comm->moved : NULL,
fr->nbv->grp[eintLocal].kernel_type,
fr->nbv->grp[eintLocal].nbat);
}
dd_sort_state(dd, fr->cg_cm, fr, state_local,
bResortAll ? -1 : ncg_home_old);
+
+ /* After sorting and compacting we set the correct size */
+ dd_resize_state(state_local, f, dd->nat_home);
+
/* Rebuild all the indices */
ga2la_clear(dd->ga2la);
ncgindex_set = 0;
/*
write_dd_pdb("dd_home",step,"dump",top_global,cr,
- -1,state_local->x,state_local->box);
+ -1,as_rvec_array(state_local->x.data()),state_local->box);
*/
wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
comm->cellsize_min, np,
fr,
- fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : state_local->x,
+ fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : as_rvec_array(state_local->x.data()),
vsite, top_global, top_local);
wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
* or constraint communication.
*/
state_local->natoms = comm->nat[ddnatNR-1];
- if (state_local->natoms > state_local->nalloc)
- {
- dd_realloc_state(state_local, f, state_local->natoms);
- }
+
+ dd_resize_state(state_local, f, state_local->natoms);
if (fr->bF_NoVirSum)
{
* the last vsite construction, we need to communicate the constructing
* atom coordinates again (for spreading the forces this MD step).
*/
- dd_move_x_vsites(dd, state_local->box, state_local->x);
+ dd_move_x_vsites(dd, state_local->box, as_rvec_array(state_local->x.data()));
wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
{
- dd_move_x(dd, state_local->box, state_local->x);
+ dd_move_x(dd, state_local->box, as_rvec_array(state_local->x.data()));
write_dd_pdb("dd_dump", step, "dump", top_global, cr,
- -1, state_local->x, state_local->box);
+ -1, as_rvec_array(state_local->x.data()), state_local->box);
}
/* Store the partitioning step */
struct t_commrec;
struct t_inputrec;
-#ifdef __cplusplus
-extern "C" {
-#endif
-
/*! \brief Returns the global topology atom number belonging to local atom index i.
*
* This function is intended for writing ASCII output
/*! \brief Collects local rvec arrays \p lv to \p v on the master rank */
void dd_collect_vec(struct gmx_domdec_t *dd,
- t_state *state_local, rvec *lv, rvec *v);
+ t_state *state_local, const PaddedRVecVector *lv, rvec *v);
+
+/*! \brief Collects local rvec arrays \p lv to \p v on the master rank */
+void dd_collect_vec(struct gmx_domdec_t *dd,
+ t_state *state_local, const PaddedRVecVector *lv, PaddedRVecVector *v);
/*! \brief Collects the local state \p state_local to \p state on the master rank */
void dd_collect_state(struct gmx_domdec_t *dd,
const gmx_mtop_t *top_global,
const t_inputrec *ir,
t_state *state_local,
- rvec **f,
+ PaddedRVecVector *f,
t_mdatoms *mdatoms,
gmx_localtop_t *top_local,
t_forcerec *fr,
const t_block *cgs, const rvec *x,
gmx_ddbox_t *ddbox);
-#ifdef __cplusplus
-}
-#endif
-
#endif
#endif
}
-void dd_scatter(gmx_domdec_t gmx_unused *dd, int gmx_unused nbytes, void gmx_unused *src, void *dest)
+void dd_scatter(gmx_domdec_t gmx_unused *dd, int gmx_unused nbytes, const void gmx_unused *src, void *dest)
{
#if GMX_MPI
if (dd->nnodes > 1)
{
- MPI_Scatter(src, nbytes, MPI_BYTE,
+ /* Some MPI implementions don't specify const */
+ MPI_Scatter(const_cast<void *>(src), nbytes, MPI_BYTE,
dest, nbytes, MPI_BYTE,
DDMASTERRANK(dd), dd->mpi_comm_all);
}
}
}
-void dd_gather(gmx_domdec_t gmx_unused *dd, int gmx_unused nbytes, void gmx_unused *src, void gmx_unused *dest)
+void dd_gather(gmx_domdec_t gmx_unused *dd, int gmx_unused nbytes, const void gmx_unused *src, void gmx_unused *dest)
{
#if GMX_MPI
- MPI_Gather(src, nbytes, MPI_BYTE,
+ /* Some MPI implementions don't specify const */
+ MPI_Gather(const_cast<void *>(src), nbytes, MPI_BYTE,
dest, nbytes, MPI_BYTE,
DDMASTERRANK(dd), dd->mpi_comm_all);
#endif
}
void dd_scatterv(gmx_domdec_t gmx_unused *dd,
- int gmx_unused *scounts, int gmx_unused *disps, void *sbuf,
+ int gmx_unused *scounts, int gmx_unused *disps, const void *sbuf,
int rcount, void *rbuf)
{
#if GMX_MPI
/* MPI does not allow NULL pointers */
rbuf = &dum;
}
- MPI_Scatterv(sbuf, scounts, disps, MPI_BYTE,
+ /* Some MPI implementions don't specify const */
+ MPI_Scatterv(const_cast<void *>(sbuf), scounts, disps, MPI_BYTE,
rbuf, rcount, MPI_BYTE,
DDMASTERRANK(dd), dd->mpi_comm_all);
}
}
void dd_gatherv(gmx_domdec_t gmx_unused *dd,
- int gmx_unused scount, void gmx_unused *sbuf,
+ int gmx_unused scount, const void gmx_unused *sbuf,
int gmx_unused *rcounts, int gmx_unused *disps, void gmx_unused *rbuf)
{
#if GMX_MPI
/* MPI does not allow NULL pointers */
sbuf = &dum;
}
- MPI_Gatherv(sbuf, scount, MPI_BYTE,
+ /* Some MPI implementions don't specify const */
+ MPI_Gatherv(const_cast<void *>(sbuf), scount, MPI_BYTE,
rbuf, rcounts, disps, MPI_BYTE,
DDMASTERRANK(dd), dd->mpi_comm_all);
#endif
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2008,2009,2010,2012,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2008,2009,2010,2012,2014,2015,2016, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/*! \brief Scatters \p nbytes from \p src on \p DDMASTERRANK to all PP ranks, received in \p dest */
void
-dd_scatter(struct gmx_domdec_t *dd, int nbytes, void *src, void *dest);
+dd_scatter(struct gmx_domdec_t *dd, int nbytes, const void *src, void *dest);
/*! \brief Gathers \p nbytes from \p src on all PP ranks, received in \p dest on \p DDMASTERRANK */
void
-dd_gather(struct gmx_domdec_t *dd, int nbytes, void *src, void *dest);
+dd_gather(struct gmx_domdec_t *dd, int nbytes, const void *src, void *dest);
/*! \brief Scatters \p scounts bytes from \p src on \p DDMASTERRANK to all PP ranks, receiving \p rcount bytes in \p dest.
*
* If rcount==0, rbuf is allowed to be NULL */
void
dd_scatterv(struct gmx_domdec_t *dd,
- int *scounts, int *disps, void *sbuf,
+ int *scounts, int *disps, const void *sbuf,
int rcount, void *rbuf);
/*! \brief Gathers \p rcount bytes from \p src on all PP ranks, received in \p scounts bytes in \p dest on \p DDMASTERRANK.
* If scount==0, sbuf is allowed to be NULL */
void
dd_gatherv(struct gmx_domdec_t *dd,
- int scount, void *sbuf,
+ int scount, const void *sbuf,
int *rcounts, int *disps, void *rbuf);
#endif
print_missing_interactions_atoms(fplog, cr, top_global, &top_local->idef);
write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
- -1, state_local->x, state_local->box);
+ -1, as_rvec_array(state_local->x.data()), state_local->box);
std::string errorMessage;
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
+#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/df_history.h"
#include "gromacs/mdtypes/energyhistory.h"
return do_cpte_reals_low(xd, cptp, ecpt, sflags, n, NULL, v, list, ecprREAL);
}
+/* This function stores n along with the reals for reading,
+ * but on reading it assumes that n matches the value in the checkpoint file,
+ * a fatal error is generated when this is not the case.
+ */
+static int do_cpte_reals(XDR *xd, int cptp, int ecpt, int sflags,
+ int n, std::vector<real> *v, FILE *list)
+{
+ real *v_real;
+ if (list == NULL && (sflags & (1 << ecpt)))
+ {
+ /* Resizes on read, on write the size should already be n */
+ v->resize(n);
+ v_real = v->data();
+ }
+ else
+ {
+ v_real = NULL;
+ }
+
+ return do_cpte_reals_low(xd, cptp, ecpt, sflags, n, NULL, &v_real, list, ecprREAL);
+}
+
/* This function does the same as do_cpte_reals,
* except that on reading it ignores the passed value of *n
* and stored the value read from the checkpoint file in *n.
return 0;
}
+static int do_cpte_doubles(XDR *xd, int cptp, int ecpt, int sflags,
+ int n, std::vector<double> *v, FILE *list)
+{
+ double *v_double;
+ if (list == NULL && (sflags & (1 << ecpt)))
+ {
+ /* Resizes on read, on write the size should already be n */
+ v->resize(n);
+ v_double = v->data();
+ }
+ else
+ {
+ v_double = NULL;
+ }
+
+ return do_cpte_doubles(xd, cptp, ecpt, sflags, n, &v_double, list);
+}
+
static int do_cpte_double(XDR *xd, int cptp, int ecpt, int sflags,
double *r, FILE *list)
{
static int do_cpte_rvecs(XDR *xd, int cptp, int ecpt, int sflags,
- int n, rvec **v, FILE *list)
+ int n, std::vector<gmx::RVec> *v, FILE *list)
{
+ rvec *v_rvec;
+
+ if (list == NULL && (sflags & (1 << ecpt)))
+ {
+ /* We resize the vector here to avoid pointer reallocation in
+ * do_cpte_reals_low. Note the we allocate 1 element extra for SIMD.
+ */
+ v->resize(n + 1);
+ v_rvec = as_rvec_array(v->data());
+ }
+ else
+ {
+ v_rvec = NULL;
+ }
+
return do_cpte_reals_low(xd, cptp, ecpt, sflags,
- n*DIM, NULL, (real **)v, list, ecprRVEC);
+ n*DIM, NULL, (real **)(&v_rvec), list, ecprRVEC);
}
static int do_cpte_matrix(XDR *xd, int cptp, int ecpt, int sflags,
fr->bX = (state.flags & (1<<estX));
if (fr->bX)
{
- fr->x = state.x;
- state.x = NULL;
+ fr->x = getRvecArrayFromPaddedRVecVector(&state.x, state.natoms);
}
fr->bV = (state.flags & (1<<estV));
if (fr->bV)
{
- fr->v = state.v;
- state.v = NULL;
+ fr->v = getRvecArrayFromPaddedRVecVector(&state.v, state.natoms);
}
fr->bF = FALSE;
fr->bBox = (state.flags & (1<<estBOX));
{
copy_mat(state.box, fr->box);
}
- done_state(&state);
}
void list_checkpoint(const char *fn, FILE *out)
{
gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
}
-
- done_state(&state);
}
/* This routine cannot print tons of data, since it is called before the log file is opened. */
{
gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
}
- done_state(&state);
}
}
static int do_tpx(t_fileio *fio, gmx_bool bRead,
- t_inputrec *ir, t_state *state, gmx_mtop_t *mtop,
- gmx_bool bXVallocated)
+ t_inputrec *ir, t_state *state, rvec *x, rvec *v,
+ gmx_mtop_t *mtop)
{
t_tpxheader tpx;
gmx_mtop_t dum_top;
gmx_bool TopOnlyOK;
- rvec *xptr, *vptr;
int ePBC;
gmx_bool bPeriodicMols;
if (!bRead)
{
+ GMX_RELEASE_ASSERT(x == NULL && v == NULL, "Passing separate x and v pointers to do_tpx() is not supported when writing");
+
tpx.natoms = state->natoms;
tpx.ngtc = state->ngtc;
tpx.fep_state = state->fep_state;
tpx.lambda = state->lambda[efptFEP];
tpx.bIr = (ir != NULL);
tpx.bTop = (mtop != NULL);
- tpx.bX = (state->x != NULL);
- tpx.bV = (state->v != NULL);
+ tpx.bX = (state->flags & (1 << estX));
+ tpx.bV = (state->flags & (1 << estV));
tpx.bF = FALSE;
tpx.bBox = TRUE;
}
+ else
+ {
+ GMX_RELEASE_ASSERT(!(x == NULL && v != NULL), "Passing x==NULL and v!=NULL is not supported");
+ }
TopOnlyOK = (ir == NULL);
if (bRead)
{
state->flags = 0;
- if (bXVallocated)
+ if (x != NULL)
{
- xptr = state->x;
- vptr = state->v;
init_state(state, 0, tpx.ngtc, 0, 0, 0);
- state->natoms = tpx.natoms;
- state->nalloc = tpx.natoms;
- state->x = xptr;
- state->v = vptr;
}
else
{
}
}
+ if (x == NULL)
+ {
+ x = as_rvec_array(state->x.data());
+ v = as_rvec_array(state->v.data());
+ }
+
#define do_test(fio, b, p) if (bRead && (p != NULL) && !b) gmx_fatal(FARGS, "No %s in %s",#p, gmx_fio_getname(fio))
do_test(fio, tpx.bBox, state->box);
done_mtop(&dum_top, TRUE);
}
}
- do_test(fio, tpx.bX, state->x);
+ do_test(fio, tpx.bX, x);
if (tpx.bX)
{
if (bRead)
{
state->flags |= (1<<estX);
}
- gmx_fio_ndo_rvec(fio, state->x, state->natoms);
+ gmx_fio_ndo_rvec(fio, x, tpx.natoms);
}
- do_test(fio, tpx.bV, state->v);
+ do_test(fio, tpx.bV, v);
if (tpx.bV)
{
if (bRead)
{
state->flags |= (1<<estV);
}
- gmx_fio_ndo_rvec(fio, state->v, state->natoms);
+ gmx_fio_ndo_rvec(fio, v, tpx.natoms);
}
// No need to run do_test when the last argument is NULL
{
rvec *dummyForces;
snew(dummyForces, state->natoms);
- gmx_fio_ndo_rvec(fio, dummyForces, state->natoms);
+ gmx_fio_ndo_rvec(fio, dummyForces, tpx.natoms);
sfree(dummyForces);
}
fio = open_tpx(fn, "w");
do_tpx(fio, FALSE,
const_cast<t_inputrec *>(ir),
- const_cast<t_state *>(state),
- const_cast<gmx_mtop_t *>(mtop), FALSE);
+ const_cast<t_state *>(state), NULL, NULL,
+ const_cast<gmx_mtop_t *>(mtop));
close_tpx(fio);
}
t_fileio *fio;
fio = open_tpx(fn, "r");
- do_tpx(fio, TRUE, ir, state, mtop, FALSE);
+ do_tpx(fio, TRUE, ir, state, NULL, NULL, mtop);
close_tpx(fio);
}
rvec *x, rvec *v, gmx_mtop_t *mtop)
{
t_fileio *fio;
- t_state state;
+ t_state state {};
int ePBC;
- state.x = x;
- state.v = v;
fio = open_tpx(fn, "r");
- ePBC = do_tpx(fio, TRUE, ir, &state, mtop, TRUE);
+ ePBC = do_tpx(fio, TRUE, ir, &state, x, v, mtop);
close_tpx(fio);
- *natoms = state.natoms;
+ *natoms = mtop->natoms;
if (box)
{
copy_mat(state.box, box);
}
- state.x = NULL;
- state.v = NULL;
- done_state(&state);
return ePBC;
}
}
/* Prepare an x and q array with only the charged atoms */
- ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
+ ncharges = prepare_x_q(&q, &x, mtop, as_rvec_array(state->x.data()), cr);
if (MASTER(cr))
{
calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
}
t_topology *conftop;
+ rvec *x = NULL;
+ rvec *v = NULL;
snew(conftop, 1);
init_state(state, 0, 0, 0, 0, 0);
- read_tps_conf(confin, conftop, NULL, &state->x, &state->v, state->box, FALSE);
- state->natoms = state->nalloc = conftop->atoms.nr;
+ read_tps_conf(confin, conftop, NULL, &x, &v, state->box, FALSE);
+ state->natoms = conftop->atoms.nr;
if (state->natoms != sys->natoms)
{
gmx_fatal(FARGS, "number of coordinates in coordinate file (%s, %d)\n"
" does not match topology (%s, %d)",
confin, state->natoms, topfile, sys->natoms);
}
+ /* It would be nice to get rid of the copies below, but we don't know
+ * a priori if the number of atoms in confin matches what we expect.
+ */
+ state->flags |= (1 << estX);
+ state->x.resize(state->natoms);
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(x[i], state->x[i]);
+ }
+ sfree(x);
+ if (v != NULL)
+ {
+ state->flags |= (1 << estV);
+ state->v.resize(state->natoms);
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(v[i], state->v[i]);
+ }
+ sfree(v);
+ }
/* This call fixes the box shape for runs with pressure scaling */
set_box_rel(ir, state);
opts->seed = static_cast<int>(gmx::makeRandomSeed());
fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
}
- maxwell_speed(opts->tempi, opts->seed, sys, state->v);
+ state->flags |= (1 << estV);
+ maxwell_speed(opts->tempi, opts->seed, sys, as_rvec_array(state->v.data()));
- stop_cm(stdout, state->natoms, mass, state->x, state->v);
+ stop_cm(stdout, state->natoms, mass, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()));
sfree(mass);
}
t_inputrec *ir;
int nvsite, comb, mt;
t_params *plist;
- t_state *state;
matrix box;
real fudgeQQ;
double reppow;
{
gmx_fatal(FARGS, "%s does not exist", fn);
}
- snew(state, 1);
+
+ t_state state {};
new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
- opts, ir, bZero, bGenVel, bVerbose, state,
+ opts, ir, bZero, bGenVel, bVerbose, &state,
atype, sys, &nmi, &mi, &intermolecular_interactions,
plist, &comb, &reppow, &fudgeQQ,
opts->bMorse,
/* Check velocity for virtual sites and shells */
if (bGenVel)
{
- check_vel(sys, state->v);
+ check_vel(sys, as_rvec_array(state.v.data()));
}
/* check for shells and inpurecs */
}
else
{
- buffer_temp = calc_temp(sys, ir, state->v);
+ buffer_temp = calc_temp(sys, ir, as_rvec_array(state.v.data()));
}
if (buffer_temp > 0)
{
}
}
- set_verlet_buffer(sys, ir, buffer_temp, state->box, wi);
+ set_verlet_buffer(sys, ir, buffer_temp, state.box, wi);
}
}
}
/* Init the temperature coupling state */
- init_gtc_state(state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
+ init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
if (bVerbose)
{
fprintf(stderr, "getting data from old trajectory ...\n");
}
cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
- bNeedVel, bGenVel, fr_time, ir, state, sys, oenv);
+ bNeedVel, bGenVel, fr_time, ir, &state, sys, oenv);
}
if (ir->ePBC == epbcXY && ir->nwall != 2)
{
- clear_rvec(state->box[ZZ]);
+ clear_rvec(state.box[ZZ]);
}
if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
{
set_warning_line(wi, mdparin, -1);
- check_chargegroup_radii(sys, ir, state->x, wi);
+ check_chargegroup_radii(sys, ir, as_rvec_array(state.x.data()), wi);
}
if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
{
/* Calculate the optimal grid dimensions */
- copy_mat(state->box, box);
+ copy_mat(state.box, box);
if (ir->ePBC == epbcXY && ir->nwall == 2)
{
svmul(ir->wall_ewald_zfac, box[ZZ], box[ZZ]);
potentially conflict if not handled correctly. */
if (ir->efep != efepNO)
{
- state->fep_state = ir->fepvals->init_fep_state;
+ state.fep_state = ir->fepvals->init_fep_state;
for (i = 0; i < efptNR; i++)
{
/* init_lambda trumps state definitions*/
if (ir->fepvals->init_lambda >= 0)
{
- state->lambda[i] = ir->fepvals->init_lambda;
+ state.lambda[i] = ir->fepvals->init_lambda;
}
else
{
}
else
{
- state->lambda[i] = ir->fepvals->all_lambda[i][state->fep_state];
+ state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
}
}
}
if (ir->bPull)
{
- pull = set_pull_init(ir, sys, state->x, state->box, state->lambda[efptMASS], oenv);
+ pull = set_pull_init(ir, sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], oenv);
}
/* Modules that supply external potential for pull coordinates
if (ir->bRot)
{
- set_reference_positions(ir->rot, state->x, state->box,
+ set_reference_positions(ir->rot, as_rvec_array(state.x.data()), state.box,
opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
wi);
}
if (EEL_PME(ir->coulombtype))
{
- float ratio = pme_load_estimate(sys, ir, state->box);
+ float ratio = pme_load_estimate(sys, ir, state.box);
fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
/* With free energy we might need to do PME both for the A and B state
* charges. This will double the cost, but the optimal performance will
}
done_warning(wi, FARGS);
- write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, state, sys);
+ write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
/* Output IMD group, if bIMD is TRUE */
- write_IMDgroup_to_file(ir->bIMD, ir, state, sys, NFILE, fnm);
+ write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
- done_state(state);
- sfree(state);
done_atomtype(atype);
done_mtop(sys, TRUE);
done_inputrec_strings();
{
IMDatoms = gmx_mtop_global_atoms(sys);
write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
- state->x, state->v, ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
}
}
#define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
#define nblock_bc(cr, nr, d) { gmx_bcast((nr)*sizeof((d)[0]), (d), (cr)); }
#define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
-/* Dirty macro with bAlloc not as an argument */
-#define nblock_abc(cr, nr, d) { if (bAlloc) {snew((d), (nr)); } nblock_bc(cr, (nr), (d)); }
+
+#if !GMX_DOUBLE
+static void nblock_abc(const t_commrec *cr, int numElements, real **v)
+{
+ if (!MASTER(cr))
+ {
+ snew(*v, numElements);
+ }
+ nblock_bc(cr, numElements, *v);
+}
+#endif
+
+static void nblock_abc(const t_commrec *cr, int numElements, double **v)
+{
+ if (!MASTER(cr))
+ {
+ snew(*v, numElements);
+ }
+ nblock_bc(cr, numElements, *v);
+}
+
+static void nblock_abc(const t_commrec *cr, int numElements, std::vector<double> *v)
+{
+ if (!MASTER(cr))
+ {
+ v->resize(numElements);
+ }
+ gmx_bcast(numElements*sizeof(double), v->data(), cr);
+}
static void bc_cstring(const t_commrec *cr, char **s)
{
}
}
+static void bcastPaddedRVecVector(const t_commrec *cr, PaddedRVecVector *v, unsigned int n)
+{
+ (*v).resize(n + 1);
+ nblock_bc(cr, n, as_rvec_array(v->data()));
+}
+
void bcast_state(const t_commrec *cr, t_state *state)
{
int i, nnht, nnhtp;
- gmx_bool bAlloc;
if (!PAR(cr) || (cr->nnodes - cr->npmenodes <= 1))
{
block_bc(cr, state->nnhpres);
block_bc(cr, state->nhchainlength);
block_bc(cr, state->flags);
- if (state->lambda == NULL)
- {
- snew_bc(cr, state->lambda, efptNR)
- }
+ state->lambda.resize(efptNR);
if (cr->dd)
{
nnht = (state->ngtc)*(state->nhchainlength);
nnhtp = (state->nnhpres)*(state->nhchainlength);
- /* We still need to allocate the arrays in state for non-master
- * ranks, which is done (implicitly via bAlloc) in the dirty,
- * dirty nblock_abc macro. */
- bAlloc = !MASTER(cr);
- if (bAlloc)
- {
- state->nalloc = state->natoms;
- }
for (i = 0; i < estNR; i++)
{
if (state->flags & (1<<i))
{
switch (i)
{
- case estLAMBDA: nblock_bc(cr, efptNR, state->lambda); break;
+ case estLAMBDA: nblock_bc(cr, efptNR, state->lambda.data()); break;
case estFEPSTATE: block_bc(cr, state->fep_state); break;
case estBOX: block_bc(cr, state->box); break;
case estBOX_REL: block_bc(cr, state->box_rel); break;
case estPRES_PREV: block_bc(cr, state->pres_prev); break;
case estSVIR_PREV: block_bc(cr, state->svir_prev); break;
case estFVIR_PREV: block_bc(cr, state->fvir_prev); break;
- case estNH_XI: nblock_abc(cr, nnht, state->nosehoover_xi); break;
- case estNH_VXI: nblock_abc(cr, nnht, state->nosehoover_vxi); break;
- case estNHPRES_XI: nblock_abc(cr, nnhtp, state->nhpres_xi); break;
- case estNHPRES_VXI: nblock_abc(cr, nnhtp, state->nhpres_vxi); break;
- case estTC_INT: nblock_abc(cr, state->ngtc, state->therm_integral); break;
+ case estNH_XI: nblock_abc(cr, nnht, &state->nosehoover_xi); break;
+ case estNH_VXI: nblock_abc(cr, nnht, &state->nosehoover_vxi); break;
+ case estNHPRES_XI: nblock_abc(cr, nnhtp, &state->nhpres_xi); break;
+ case estNHPRES_VXI: nblock_abc(cr, nnhtp, &state->nhpres_vxi); break;
+ case estTC_INT: nblock_abc(cr, state->ngtc, &state->therm_integral); break;
case estVETA: block_bc(cr, state->veta); break;
case estVOL0: block_bc(cr, state->vol0); break;
- case estX: nblock_abc(cr, state->natoms, state->x); break;
- case estV: nblock_abc(cr, state->natoms, state->v); break;
- case estCGP: nblock_abc(cr, state->natoms, state->cg_p); break;
+ case estX: bcastPaddedRVecVector(cr, &state->x, state->natoms);
+ case estV: bcastPaddedRVecVector(cr, &state->v, state->natoms);
+ case estCGP: bcastPaddedRVecVector(cr, &state->cg_p, state->natoms);
case estDISRE_INITF: block_bc(cr, state->hist.disre_initf); break;
case estDISRE_RM3TAV:
block_bc(cr, state->hist.ndisrepairs);
- nblock_abc(cr, state->hist.ndisrepairs, state->hist.disre_rm3tav);
+ nblock_abc(cr, state->hist.ndisrepairs, &state->hist.disre_rm3tav);
break;
case estORIRE_INITF: block_bc(cr, state->hist.orire_initf); break;
case estORIRE_DTAV:
block_bc(cr, state->hist.norire_Dtav);
- nblock_abc(cr, state->hist.norire_Dtav, state->hist.orire_Dtav);
+ nblock_abc(cr, state->hist.norire_Dtav, &state->hist.orire_Dtav);
break;
default:
gmx_fatal(FARGS,
constr->ed = ed;
if (ed != NULL || state->edsamstate != NULL)
{
- init_edsam(mtop, ir, cr, ed, state->x, state->box, state->edsamstate);
+ init_edsam(mtop, ir, cr, ed, as_rvec_array(state->x.data()), state->box, state->edsamstate);
}
constr->warn_mtop = mtop;
}
}
-t_state *init_bufstate(const t_state *template_state)
-{
- t_state *state;
- int nc = template_state->nhchainlength;
- snew(state, 1);
- snew(state->nosehoover_xi, nc*template_state->ngtc);
- snew(state->nosehoover_vxi, nc*template_state->ngtc);
- snew(state->therm_integral, template_state->ngtc);
- snew(state->nhpres_xi, nc*template_state->nnhpres);
- snew(state->nhpres_vxi, nc*template_state->nnhpres);
-
- return state;
-}
-
-void destroy_bufstate(t_state *state)
-{
- sfree(state->x);
- sfree(state->v);
- sfree(state->nosehoover_xi);
- sfree(state->nosehoover_vxi);
- sfree(state->therm_integral);
- sfree(state->nhpres_xi);
- sfree(state->nhpres_vxi);
- sfree(state);
-}
-
void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
gmx_enerdata_t *enerd, t_state *state,
tensor vir, t_mdatoms *md,
break;
case etrtBARONHC:
case etrtBARONHC2:
- NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
- state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
+ NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
+ state->nhpres_vxi.data(), NULL, &(state->veta), MassQ, FALSE);
break;
case etrtNHC:
case etrtNHC2:
- NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
- state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
+ NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
+ state->nosehoover_vxi.data(), scalefac, NULL, MassQ, (ir->eI == eiVV));
/* need to rescale the kinetic energies and velocities here. Could
scale the velocities later, but we need them scaled in order to
produce the correct outputs, so we'll scale them here. */
}
}
-void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
+void sum_dhdl(gmx_enerdata_t *enerd, const std::vector<real> *lambda, t_lambda *fepvals)
{
int i, j, index;
double dlam;
for (j = 0; j < efptNR; j++)
{
/* Note that this loop is over all dhdl components, not just the separated ones */
- dlam = (fepvals->all_lambda[j][i]-lambda[j]);
+ dlam = (fepvals->all_lambda[j][i] - (*lambda)[j]);
enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
if (debug)
{
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdtypes/fcdata.h"
#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/timing/wallcycle.h"
struct gmx_edsam;
void sum_epot(gmx_grppairener_t *grpp, real *epot);
/* Locally sum the non-bonded potential energy terms */
-void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals);
+void sum_dhdl(gmx_enerdata_t *enerd, const std::vector<real> *lambda, t_lambda *fepvals);
/* Sum the free energy contributions */
/* Compute the average C6 and C12 params for LJ corrections */
gmx_int64_t step, struct 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, struct t_graph *graph,
+ std::vector<real> *lambda, struct t_graph *graph,
t_forcerec *fr,
gmx_vsite_t *vsite, rvec mu_tot,
double t, FILE *field, struct gmx_edsam *ed,
return nmin;
}
-void copy_coupling_state(t_state *statea, t_state *stateb,
- gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts)
-{
-
- /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
-
- int i, j, nc;
-
- /* Make sure we have enough space for x and v */
- if (statea->nalloc > stateb->nalloc)
- {
- stateb->nalloc = statea->nalloc;
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- srenew(stateb->x, stateb->nalloc + 1);
- srenew(stateb->v, stateb->nalloc + 1);
- }
-
- stateb->natoms = statea->natoms;
- stateb->ngtc = statea->ngtc;
- stateb->nnhpres = statea->nnhpres;
- stateb->veta = statea->veta;
- if (ekinda)
- {
- copy_mat(ekinda->ekin, ekindb->ekin);
- for (i = 0; i < stateb->ngtc; i++)
- {
- ekindb->tcstat[i].T = ekinda->tcstat[i].T;
- ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
- copy_mat(ekinda->tcstat[i].ekinh, ekindb->tcstat[i].ekinh);
- copy_mat(ekinda->tcstat[i].ekinf, ekindb->tcstat[i].ekinf);
- ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
- ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
- ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
- }
- }
- copy_rvecn(statea->x, stateb->x, 0, stateb->natoms);
- copy_rvecn(statea->v, stateb->v, 0, stateb->natoms);
- copy_mat(statea->box, stateb->box);
- copy_mat(statea->box_rel, stateb->box_rel);
- copy_mat(statea->boxv, stateb->boxv);
-
- for (i = 0; i < stateb->ngtc; i++)
- {
- nc = i*opts->nhchainlength;
- for (j = 0; j < opts->nhchainlength; j++)
- {
- stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
- stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
- }
- }
- if (stateb->nhpres_xi != NULL)
- {
- for (i = 0; i < stateb->nnhpres; i++)
- {
- nc = i*opts->nhchainlength;
- for (j = 0; j < opts->nhchainlength; j++)
- {
- stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
- stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
- }
- }
- }
-}
-
real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
{
real quantity = 0;
quantity = NPT_energy(ir, state, MassQ);
break;
case etcVRESCALE:
- quantity = vrescale_energy(&(ir->opts), state->therm_integral);
+ quantity = vrescale_energy(&(ir->opts), state->therm_integral.data());
break;
default:
break;
if (bStopCM)
{
calc_vcm_grp(0, mdatoms->homenr, mdatoms,
- state->x, state->v, vcm);
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), vcm);
}
if (bTemp || bStopCM || bPres || bEner || bConstrain)
{
correct_ekin(debug,
0, mdatoms->homenr,
- state->v, vcm->group_p[0],
+ as_rvec_array(state->v.data()), vcm->group_p[0],
mdatoms->massT, mdatoms->tmass, ekind->ekin);
}
{
check_cm_grp(fplog, vcm, ir, 1);
do_stopcm_grp(0, mdatoms->homenr, mdatoms->cVCM,
- state->x, state->v, vcm);
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), vcm);
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
state->flags |= (1<<estFEPSTATE);
}
state->flags |= (1<<estX);
- if (state->lambda == NULL)
- {
- snew(state->lambda, efptNR);
- }
- if (state->x == NULL)
- {
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- snew(state->x, state->nalloc + 1);
- }
+ state->lambda.resize(efptNR);
+ GMX_RELEASE_ASSERT(state->x.size() >= static_cast<unsigned int>(state->natoms), "We should start a run with an initialized state->x");
if (EI_DYNAMICS(ir->eI))
{
state->flags |= (1<<estV);
- if (state->v == NULL)
- {
- snew(state->v, state->nalloc + 1);
- }
+ state->v.resize(state->natoms + 1);
}
if (ir->eI == eiCG)
{
state->flags |= (1<<estCGP);
- if (state->cg_p == NULL)
- {
- /* cg_p is not stored in the tpx file, so we need to allocate it */
- snew(state->cg_p, state->nalloc + 1);
- }
+ /* cg_p is not stored in the tpx file, so we need to allocate it */
+ state->cg_p.resize(state->natoms + 1);
}
state->nnhpres = 0;
int multisim_min(const gmx_multisim_t *ms, int nmin, int n);
/* Set an appropriate value for n across the whole multi-simulation */
-void copy_coupling_state(t_state *statea, t_state *stateb,
- gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts);
-/* Copy stuff from state A to state B */
-
void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
t_forcerec *fr, gmx_ekindata_t *ekind,
t_state *state, t_mdatoms *mdatoms,
gmx_mtop_t *top_global,
gmx_int64_t step, double t,
t_state *state_local, t_state *state_global,
- rvec *f_local)
+ PaddedRVecVector *f_local)
{
rvec *f_global;
{
if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))
{
- dd_collect_vec(cr->dd, state_local, state_local->x,
- state_global->x);
+ dd_collect_vec(cr->dd, state_local, &state_local->x,
+ &state_global->x);
}
if (mdof_flags & MDOF_V)
{
- dd_collect_vec(cr->dd, state_local, state_local->v,
- state_global->v);
+ dd_collect_vec(cr->dd, state_local, &state_local->v,
+ &state_global->v);
}
}
f_global = of->f_global;
/* We have the whole state locally: copy the local state pointer */
state_global = state_local;
- f_global = f_local;
+ f_global = as_rvec_array(f_local->data());
}
if (MASTER(cr))
if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
{
+ const rvec *x = (mdof_flags & MDOF_X) ? as_rvec_array(state_global->x.data()) : NULL;
+ const rvec *v = (mdof_flags & MDOF_V) ? as_rvec_array(state_global->v.data()) : NULL;
+ const rvec *f = (mdof_flags & MDOF_F) ? f_global : NULL;
+
if (of->fp_trn)
{
gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],
state_local->box, top_global->natoms,
- (mdof_flags & MDOF_X) ? state_global->x : NULL,
- (mdof_flags & MDOF_V) ? state_global->v : NULL,
- (mdof_flags & MDOF_F) ? f_global : NULL);
+ x, v, f);
if (gmx_fio_flush(of->fp_trn) != 0)
{
gmx_file("Cannot write trajectory; maybe you are out of disk space?");
gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],
state_local->box,
top_global->natoms,
- (mdof_flags & MDOF_X) ? state_global->x : NULL,
- (mdof_flags & MDOF_V) ? state_global->v : NULL,
- (mdof_flags & MDOF_F) ? f_global : NULL);
+ x, v, f);
}
/* If only a TNG file is open for compressed coordinate output (no uncompressed
coordinate output) also write forces and velocities to it. */
gmx_fwrite_tng(of->tng_low_prec, FALSE, step, t, state_local->lambda[efptFEP],
state_local->box,
top_global->natoms,
- (mdof_flags & MDOF_X) ? state_global->x : NULL,
- (mdof_flags & MDOF_V) ? state_global->v : NULL,
- (mdof_flags & MDOF_F) ? f_global : NULL);
+ x, v, f);
}
}
if (mdof_flags & MDOF_X_COMPRESSED)
{
/* We are writing the positions of all of the atoms to
the compressed output */
- xxtc = state_global->x;
+ xxtc = as_rvec_array(state_global->x.data());
}
else
{
gmx_mtop_t *top_global,
gmx_int64_t step, double t,
t_state *state_local, t_state *state_global,
- rvec *f_local);
+ PaddedRVecVector *f_local);
#define MDOF_X (1<<0)
#define MDOF_V (1<<1)
//! Utility structure for manipulating states during EM
typedef struct {
//! Copy of the global state
- t_state s;
+ t_state s;
//! Force array
- rvec *f;
+ PaddedRVecVector f;
//! Potential energy
- real epot;
+ real epot;
//! Norm of the force
- real fnorm;
+ real fnorm;
//! Maximum force
- real fmax;
+ real fmax;
//! Direction
- int a_fmax;
+ int a_fmax;
} em_state_t;
-//! Initiate em_state_t structure and return pointer to it
-static em_state_t *init_em_state()
-{
- em_state_t *ems;
-
- snew(ems, 1);
-
- /* does this need to be here? Should the array be declared differently (staticaly)in the state definition? */
- snew(ems->s.lambda, efptNR);
-
- return ems;
-}
-
//! Print the EM starting conditions
static void print_em_start(FILE *fplog,
t_commrec *cr,
//! Print message about convergence of the EM
static void print_converged(FILE *fp, const char *alg, real ftol,
gmx_int64_t count, gmx_bool bDone, gmx_int64_t nsteps,
- real epot, real fmax, int nfmax, real fnorm)
+ const em_state_t *ems, double sqrtNumAtoms)
{
char buf[STEPSTRSIZE];
}
#if GMX_DOUBLE
- fprintf(fp, "Potential Energy = %21.14e\n", epot);
- fprintf(fp, "Maximum force = %21.14e on atom %d\n", fmax, nfmax+1);
- fprintf(fp, "Norm of force = %21.14e\n", fnorm);
+ fprintf(fp, "Potential Energy = %21.14e\n", ems->epot);
+ fprintf(fp, "Maximum force = %21.14e on atom %d\n", ems->fmax, ems->a_fmax + 1);
+ fprintf(fp, "Norm of force = %21.14e\n", ems->fnorm/sqrtNumAtoms);
#else
- fprintf(fp, "Potential Energy = %14.7e\n", epot);
- fprintf(fp, "Maximum force = %14.7e on atom %d\n", fmax, nfmax+1);
- fprintf(fp, "Norm of force = %14.7e\n", fnorm);
+ fprintf(fp, "Potential Energy = %14.7e\n", ems->epot);
+ fprintf(fp, "Maximum force = %14.7e on atom %d\n", ems->fmax, ems->a_fmax + 1);
+ fprintf(fp, "Norm of force = %14.7e\n", ems->fnorm/sqrtNumAtoms);
#endif
}
//! Compute the norm and max of the force array in parallel
static void get_f_norm_max(t_commrec *cr,
- t_grpopts *opts, t_mdatoms *mdatoms, rvec *f,
+ t_grpopts *opts, t_mdatoms *mdatoms, const rvec *f,
real *fnorm, real *fmax, int *a_fmax)
{
double fnorm2, *sum;
t_grpopts *opts, t_mdatoms *mdatoms,
em_state_t *ems)
{
- get_f_norm_max(cr, opts, mdatoms, ems->f, &ems->fnorm, &ems->fmax, &ems->a_fmax);
+ get_f_norm_max(cr, opts, mdatoms, as_rvec_array(ems->f.data()),
+ &ems->fnorm, &ems->fmax, &ems->a_fmax);
}
//! Initialize the energy minimization
t_commrec *cr, t_inputrec *ir,
t_state *state_global, gmx_mtop_t *top_global,
em_state_t *ems, gmx_localtop_t **top,
- rvec **f,
t_nrnb *nrnb, rvec mu_tot,
t_forcerec *fr, gmx_enerdata_t **enerd,
t_graph **graph, t_mdatoms *mdatoms, gmx_global_stat_t *gstat,
int imdport, unsigned long gmx_unused Flags,
gmx_wallcycle_t wcycle)
{
- int i;
real dvdl_constr;
if (fplog)
state_global->ngtc = 0;
/* Initialize lambda variables */
- initialize_lambdas(fplog, ir, &(state_global->fep_state), state_global->lambda, NULL);
+ initialize_lambdas(fplog, ir, &(state_global->fep_state), &state_global->lambda, NULL);
init_nrnb(nrnb);
/* Interactive molecular dynamics */
- init_IMD(ir, cr, top_global, fplog, 1, state_global->x,
+ init_IMD(ir, cr, top_global, fplog, 1, as_rvec_array(state_global->x.data()),
nfile, fnm, NULL, imdport, Flags);
if (ir->eI == eiNM)
dd_init_local_state(cr->dd, state_global, &ems->s);
- *f = NULL;
-
/* Distribute the charge groups over the nodes from the master node */
dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
state_global, top_global, ir,
}
else
{
- snew(*f, top_global->natoms);
-
/* Just copy the state */
ems->s = *state_global;
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- snew(ems->s.x, ems->s.nalloc + 1);
- snew(ems->f, ems->s.nalloc+1);
- snew(ems->s.v, ems->s.nalloc+1);
- for (i = 0; i < state_global->natoms; i++)
- {
- copy_rvec(state_global->x[i], ems->s.x[i]);
- }
- copy_mat(state_global->box, ems->s.box);
-
+ ems->s.x.resize(ems->s.natoms + 1);
+ ems->f.resize(ems->s.natoms + 1);
snew(*top, 1);
mdAlgorithmsSetupAtomData(cr, ir, top_global, *top, fr,
dvdl_constr = 0;
constrain(PAR(cr) ? NULL : fplog, TRUE, TRUE, constr, &(*top)->idef,
ir, cr, -1, 0, 1.0, mdatoms,
- ems->s.x, ems->s.x, NULL, fr->bMolPBC, ems->s.box,
+ as_rvec_array(ems->s.x.data()),
+ as_rvec_array(ems->s.x.data()),
+ NULL,
+ fr->bMolPBC, ems->s.box,
ems->s.lambda[efptFEP], &dvdl_constr,
NULL, NULL, nrnb, econqCoord);
}
}
//! Swap two different EM states during minimization
-static void swap_em_state(em_state_t *ems1, em_state_t *ems2)
+static void swap_em_state(em_state_t **ems1, em_state_t **ems2)
{
- em_state_t tmp;
+ em_state_t *tmp;
tmp = *ems1;
*ems1 = *ems2;
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
top_global, step, (double)step,
- &state->s, state_global, state->f);
+ &state->s, state_global,
+ &state->f);
if (confout != NULL && MASTER(cr))
{
{
/* Make molecules whole only for confout writing */
do_pbc_mtop(fplog, ir->ePBC, state_global->box, top_global,
- state_global->x);
+ as_rvec_array(state_global->x.data()));
}
write_sto_conf_mtop(confout,
*top_global->name, top_global,
- state_global->x, NULL, ir->ePBC, state_global->box);
+ as_rvec_array(state_global->x.data()), NULL, ir->ePBC, state_global->box);
}
}
// \returns true when the step succeeded, false when a constraint error occurred
static bool do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
gmx_bool bMolPBC,
- em_state_t *ems1, real a, rvec *f, em_state_t *ems2,
+ em_state_t *ems1, real a, const PaddedRVecVector *force,
+ em_state_t *ems2,
gmx_constr_t constr, gmx_localtop_t *top,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_int64_t count)
s2->flags = s1->flags;
- if (s2->nalloc != s1->nalloc)
+ if (s2->natoms != s1->natoms)
{
- s2->nalloc = s1->nalloc;
+ s2->natoms = s1->natoms;
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- srenew(s2->x, s1->nalloc + 1);
- srenew(ems2->f, s1->nalloc);
+ s2->x.resize(s1->natoms + 1);
+ ems2->f.resize(s1->natoms + 1);
if (s2->flags & (1<<estCGP))
{
- srenew(s2->cg_p, s1->nalloc + 1);
+ s2->cg_p.resize(s1->natoms + 1);
}
}
+ if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size())
+ {
+ s2->cg_gl.resize(s1->cg_gl.size());
+ }
s2->natoms = s1->natoms;
copy_mat(s1->box, s2->box);
/* Copy free energy state */
- for (int i = 0; i < efptNR; i++)
- {
- s2->lambda[i] = s1->lambda[i];
- }
+ s2->lambda = s1->lambda;
copy_mat(s1->box, s2->box);
start = 0;
nthreads = gmx_omp_nthreads_get(emntUpdate);
#pragma omp parallel num_threads(nthreads)
{
- rvec *x1 = s1->x;
- rvec *x2 = s2->x;
+ const rvec *x1 = as_rvec_array(s1->x.data());
+ rvec *x2 = as_rvec_array(s2->x.data());
+ const rvec *f = as_rvec_array(force->data());
- int gf = 0;
+ int gf = 0;
#pragma omp for schedule(static) nowait
for (int i = start; i < end; i++)
{
if (s2->flags & (1<<estCGP))
{
/* Copy the CG p vector */
- rvec *p1 = s1->cg_p;
- rvec *p2 = s2->cg_p;
+ const rvec *p1 = as_rvec_array(s1->cg_p.data());
+ rvec *p2 = as_rvec_array(s2->cg_p.data());
#pragma omp for schedule(static) nowait
for (int i = start; i < end; i++)
{
if (DOMAINDECOMP(cr))
{
s2->ddp_count = s1->ddp_count;
- if (s2->cg_gl_nalloc < s1->cg_gl_nalloc)
- {
-#pragma omp barrier
- s2->cg_gl_nalloc = s1->cg_gl_nalloc;
- try
- {
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- srenew(s2->cg_gl, s2->cg_gl_nalloc + 1);
- }
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
-#pragma omp barrier
- }
- s2->ncg_gl = s1->ncg_gl;
+
+ /* OpenMP does not supported unsigned loop variables */
#pragma omp for schedule(static) nowait
- for (int i = 0; i < s2->ncg_gl; i++)
+ for (int i = 0; i < static_cast<int>(s2->cg_gl.size()); i++)
{
s2->cg_gl[i] = s1->cg_gl[i];
}
validStep =
constrain(NULL, TRUE, TRUE, constr, &top->idef,
ir, cr, count, 0, 1.0, md,
- s1->x, s2->x, NULL, bMolPBC, s2->box,
+ as_rvec_array(s1->x.data()), as_rvec_array(s2->x.data()),
+ NULL, bMolPBC, s2->box,
s2->lambda[efptBONDED], &dvdl_constr,
NULL, NULL, nrnb, econqCoord);
wallcycle_stop(wcycle, ewcCONSTR);
if (vsite)
{
- construct_vsites(vsite, ems->s.x, 1, NULL,
+ construct_vsites(vsite, as_rvec_array(ems->s.x.data()), 1, NULL,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, ems->s.box);
}
*/
do_force(fplog, cr, inputrec,
count, nrnb, wcycle, top, &top_global->groups,
- ems->s.box, ems->s.x, &ems->s.hist,
- ems->f, force_vir, mdatoms, enerd, fcd,
- ems->s.lambda, graph, fr, vsite, mu_tot, t, NULL, NULL, TRUE,
+ ems->s.box, &ems->s.x, &ems->s.hist,
+ &ems->f, force_vir, mdatoms, enerd, fcd,
+ &ems->s.lambda, graph, fr, vsite, mu_tot, t, NULL, NULL, TRUE,
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
(bNS ? GMX_FORCE_NS : 0));
/* Project out the constraint components of the force */
wallcycle_start(wcycle, ewcCONSTR);
dvdl_constr = 0;
+ rvec *f_rvec = as_rvec_array(ems->f.data());
constrain(NULL, FALSE, FALSE, constr, &top->idef,
inputrec, cr, count, 0, 1.0, mdatoms,
- ems->s.x, ems->f, ems->f, fr->bMolPBC, ems->s.box,
+ as_rvec_array(ems->s.x.data()), f_rvec, f_rvec,
+ fr->bMolPBC, ems->s.box,
ems->s.lambda[efptBONDED], &dvdl_constr,
NULL, &shake_vir, nrnb, econqForceDispl);
enerd->term[F_DVDL_CONSTR] += dvdl_constr;
enerd->term[F_PRES] =
calc_pres(fr->ePBC, inputrec->nwall, ems->s.box, ekin, vir, pres);
- sum_dhdl(enerd, ems->s.lambda, inputrec->fepvals);
+ sum_dhdl(enerd, &ems->s.lambda, inputrec->fepvals);
if (EI_ENERGY_MINIMIZATION(inputrec->eI))
{
gmx_mtop_t *top_global,
em_state_t *s_min, em_state_t *s_b)
{
- rvec *fm, *fb, *fmg;
t_block *cgs_gl;
int ncg, *cg_gl, *index, c, cg, i, a0, a1, a, gf, m;
double partsum;
fprintf(debug, "Doing reorder_partsum\n");
}
- fm = s_min->f;
- fb = s_b->f;
+ const rvec *fm = as_rvec_array(s_min->f.data());
+ const rvec *fb = as_rvec_array(s_b->f.data());
cgs_gl = dd_charge_groups_global(cr->dd);
index = cgs_gl->index;
* This conflicts with the spirit of domain decomposition,
* but to fully optimize this a much more complicated algorithm is required.
*/
+ rvec *fmg;
snew(fmg, top_global->natoms);
- ncg = s_min->s.ncg_gl;
- cg_gl = s_min->s.cg_gl;
+ ncg = s_min->s.cg_gl.size();
+ cg_gl = s_min->s.cg_gl.data();
i = 0;
for (c = 0; c < ncg; c++)
{
gmx_sum(top_global->natoms*3, fmg[0], cr);
/* Now we will determine the part of the sum for the cgs in state s_b */
- ncg = s_b->s.ncg_gl;
- cg_gl = s_b->s.cg_gl;
+ ncg = s_b->s.cg_gl.size();
+ cg_gl = s_b->s.cg_gl.data();
partsum = 0;
i = 0;
gf = 0;
gmx_mtop_t *top_global,
em_state_t *s_min, em_state_t *s_b)
{
- rvec *fm, *fb;
double sum;
- int gf, i, m;
/* This is just the classical Polak-Ribiere calculation of beta;
* it looks a bit complicated since we take freeze groups into account,
(s_min->s.ddp_count == cr->dd->ddp_count &&
s_b->s.ddp_count == cr->dd->ddp_count))
{
- fm = s_min->f;
- fb = s_b->f;
- sum = 0;
- gf = 0;
+ const rvec *fm = as_rvec_array(s_min->f.data());
+ const rvec *fb = as_rvec_array(s_b->f.data());
+ sum = 0;
+ int gf = 0;
/* This part of code can be incorrect with DD,
* since the atom ordering in s_b and s_min might differ.
*/
- for (i = 0; i < mdatoms->homenr; i++)
+ for (int i = 0; i < mdatoms->homenr; i++)
{
if (mdatoms->cFREEZE)
{
gf = mdatoms->cFREEZE[i];
}
- for (m = 0; m < DIM; m++)
+ for (int m = 0; m < DIM; m++)
{
if (!opts->nFreeze[gf][m])
{
{
const char *CG = "Polak-Ribiere Conjugate Gradients";
- em_state_t *s_min, *s_a, *s_b, *s_c;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
- rvec *p, *sf;
- double gpa, gpb, gpc, tmp, minstep;
- real fnormn;
+ double tmp, minstep;
real stepsize;
real a, b, c, beta = 0.0;
real epot_repl = 0;
tensor vir, pres;
int number_steps, neval = 0, nstcg = inputrec->nstcgsteep;
gmx_mdoutf_t outf;
- int i, m, gf, step, nminstep;
+ int m, step, nminstep;
step = 0;
- s_min = init_em_state();
- s_a = init_em_state();
- s_b = init_em_state();
- s_c = init_em_state();
+ /* Create 4 states on the stack and extract pointers that we will swap */
+ em_state_t s0 {}, s1 {}, s2 {}, s3 {};
+ em_state_t *s_min = &s0;
+ em_state_t *s_a = &s1;
+ em_state_t *s_b = &s2;
+ em_state_t *s_c = &s3;
/* Init em and store the local state in s_min */
init_em(fplog, CG, cr, inputrec,
- state_global, top_global, s_min, &top, &f,
+ state_global, top_global, s_min, &top,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
vsite, constr, NULL,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
*/
/* Calculate the new direction in p, and the gradient in this direction, gpa */
- p = s_min->s.cg_p;
- sf = s_min->f;
- gpa = 0;
- gf = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ rvec *pm = as_rvec_array(s_min->s.cg_p.data());
+ const rvec *sfm = as_rvec_array(s_min->f.data());
+ double gpa = 0;
+ int gf = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
if (mdatoms->cFREEZE)
{
{
if (!inputrec->opts.nFreeze[gf][m])
{
- p[i][m] = sf[i][m] + beta*p[i][m];
- gpa -= p[i][m]*sf[i][m];
+ pm[i][m] = sfm[i][m] + beta*pm[i][m];
+ gpa -= pm[i][m]*sfm[i][m];
/* f is negative gradient, thus the sign */
}
else
{
- p[i][m] = 0;
+ pm[i][m] = 0;
}
}
}
}
/* Calculate the norm of the search vector */
- get_f_norm_max(cr, &(inputrec->opts), mdatoms, p, &pnorm, NULL, NULL);
+ get_f_norm_max(cr, &(inputrec->opts), mdatoms, pm, &pnorm, NULL, NULL);
/* Just in case stepsize reaches zero due to numerical precision... */
if (stepsize <= 0)
* relative change in coordinate is smaller than precision
*/
minstep = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
{
tmp = 1.0;
}
- tmp = p[i][m]/tmp;
+ tmp = pm[i][m]/tmp;
minstep += tmp*tmp;
}
}
}
/* Take a trial step (new coords in s_c) */
- do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, c, s_min->s.cg_p, s_c,
+ do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, c, &s_min->s.cg_p, s_c,
constr, top, nrnb, wcycle, -1);
neval++;
mu_tot, enerd, vir, pres, -1, FALSE);
/* Calc derivative along line */
- p = s_c->s.cg_p;
- sf = s_c->f;
- gpc = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ const rvec *pc = as_rvec_array(s_c->s.cg_p.data());
+ const rvec *sfc = as_rvec_array(s_c->f.data());
+ double gpc = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
- gpc -= p[i][m]*sf[i][m]; /* f is negative gradient, thus the sign */
+ gpc -= pc[i][m]*sfc[i][m]; /* f is negative gradient, thus the sign */
}
}
/* Sum the gradient along the line across CPUs */
*
* If we already found a lower value we just skip this step and continue to the update.
*/
+ double gpb;
if (!foundlower)
{
nminstep = 0;
}
/* Take a trial step to this new point - new coords in s_b */
- do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, b, s_min->s.cg_p, s_b,
+ do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, b, &s_min->s.cg_p, s_b,
constr, top, nrnb, wcycle, -1);
neval++;
/* p does not change within a step, but since the domain decomposition
* might change, we have to use cg_p of s_b here.
*/
- p = s_b->s.cg_p;
- sf = s_b->f;
- gpb = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ const rvec *pb = as_rvec_array(s_b->s.cg_p.data());
+ const rvec *sfb = as_rvec_array(s_b->f.data());
+ gpb = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
- gpb -= p[i][m]*sf[i][m]; /* f is negative gradient, thus the sign */
+ gpb -= pb[i][m]*sfb[i][m]; /* f is negative gradient, thus the sign */
}
}
/* Sum the gradient along the line across CPUs */
if (gpb > 0)
{
/* Replace c endpoint with b */
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
c = b;
gpc = gpb;
}
else
{
/* Replace a endpoint with b */
- swap_em_state(s_b, s_a);
+ swap_em_state(&s_b, &s_a);
a = b;
gpa = gpb;
}
fprintf(debug, "CGE: C (%f) is lower than A (%f), moving C to B\n",
s_c->epot, s_a->epot);
}
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
gpb = gpc;
}
else
fprintf(debug, "CGE: A (%f) is lower than C (%f), moving A to B\n",
s_a->epot, s_c->epot);
}
- swap_em_state(s_b, s_a);
+ swap_em_state(&s_b, &s_a);
gpb = gpa;
}
fprintf(debug, "CGE: Found a lower energy %f, moving C to B\n",
s_c->epot);
}
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
gpb = gpc;
}
/* update positions */
- swap_em_state(s_min, s_b);
+ swap_em_state(&s_min, &s_b);
gpa = gpb;
/* Print it if necessary */
}
/* Send energies and positions to the IMD client if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
if (MASTER(cr))
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
- fnormn = s_min->fnorm/sqrtNumAtoms;
print_converged(stderr, CG, inputrec->em_tol, step, converged, number_steps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
print_converged(fplog, CG, inputrec->em_tol, step, converged, number_steps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
}
em_state_t ems;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
int ncorr, nmaxcorr, point, cp, neval, nminstep;
double stepsize, step_taken, gpa, gpb, gpc, tmp, minstep;
- real *rho, *alpha, *ff, *xx, *p, *s, *lastx, *lastf, **dx, **dg;
- real *xa, *xb, *xc, *fa, *fb, *fc, *xtmp, *ftmp;
+ real *rho, *alpha, *p, *s, **dx, **dg;
real a, b, c, maxdelta, delta;
- real diag, Epot0, Epot, EpotA, EpotB, EpotC;
+ real diag, Epot0;
real dgdx, dgdg, sq, yr, beta;
t_mdebin *mdebin;
gmx_bool converged;
rvec mu_tot;
- real fnorm, fmax;
gmx_bool do_log, do_ene, do_x, do_f, foundlower, *frozen;
tensor vir, pres;
int start, end, number_steps;
gmx_mdoutf_t outf;
- int i, k, m, n, nfmax, gf, step;
+ int i, k, m, n, gf, step;
int mdof_flags;
if (PAR(cr))
n = 3*state_global->natoms;
nmaxcorr = inputrec->nbfgscorr;
- /* Allocate memory */
- /* Use pointers to real so we dont have to loop over both atoms and
- * dimensions all the time...
- * x/f are allocated as rvec *, so make new x0/f0 pointers-to-real
- * that point to the same memory.
- */
- snew(xa, n);
- snew(xb, n);
- snew(xc, n);
- snew(fa, n);
- snew(fb, n);
- snew(fc, n);
snew(frozen, n);
snew(p, n);
- snew(lastx, n);
- snew(lastf, n);
snew(rho, nmaxcorr);
snew(alpha, nmaxcorr);
/* Init em */
init_em(fplog, LBFGS, cr, inputrec,
- state_global, top_global, &ems, &top, &f,
+ state_global, top_global, &ems, &top,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
vsite, constr, NULL,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
- /* Do_lbfgs is not completely updated like do_steep and do_cg,
- * so we free some memory again.
- */
- sfree(ems.s.x);
- sfree(ems.f);
-
- xx = (real *)state_global->x;
- ff = (real *)f;
start = 0;
end = mdatoms->homenr;
+ /* We need 4 working states */
+ em_state_t s0 {}, s1 {}, s2 {}, s3 {};
+ em_state_t *sa = &s0;
+ em_state_t *sb = &s1;
+ em_state_t *sc = &s2;
+ em_state_t *last = &s3;
+ /* Initialize by copying the state from ems (we could skip x and f here) */
+ *sa = ems;
+ *sb = ems;
+ *sc = ems;
+
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, LBFGS);
if (vsite)
{
- construct_vsites(vsite, state_global->x, 1, NULL,
+ construct_vsites(vsite, as_rvec_array(state_global->x.data()), 1, NULL,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state_global->box);
}
* We do not unshift, so molecules are always whole
*/
neval++;
- ems.s.x = state_global->x;
- ems.f = f;
evaluate_energy(fplog, cr,
top_global, &ems, top,
inputrec, nrnb, wcycle, gstat,
}
where();
- /* This is the starting energy */
- Epot = enerd->term[F_EPOT];
-
- fnorm = ems.fnorm;
- fmax = ems.fmax;
- nfmax = ems.a_fmax;
-
/* Set the initial step.
* since it will be multiplied by the non-normalized search direction
* vector (force vector the first time), we scale it by the
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
fprintf(stderr, "Using %d BFGS correction steps.\n\n", nmaxcorr);
- fprintf(stderr, " F-max = %12.5e on atom %d\n", fmax, nfmax+1);
- fprintf(stderr, " F-Norm = %12.5e\n", fnorm/sqrtNumAtoms);
+ fprintf(stderr, " F-max = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
+ fprintf(stderr, " F-Norm = %12.5e\n", ems.fnorm/sqrtNumAtoms);
fprintf(stderr, "\n");
/* and copy to the log file too... */
fprintf(fplog, "Using %d BFGS correction steps.\n\n", nmaxcorr);
- fprintf(fplog, " F-max = %12.5e on atom %d\n", fmax, nfmax+1);
- fprintf(fplog, " F-Norm = %12.5e\n", fnorm/sqrtNumAtoms);
+ fprintf(fplog, " F-max = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
+ fprintf(fplog, " F-Norm = %12.5e\n", ems.fnorm/sqrtNumAtoms);
fprintf(fplog, "\n");
}
point = 0;
// Set initial search direction to the force (-gradient), or 0 for frozen particles.
+ real *fInit = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
for (i = 0; i < n; i++)
{
if (!frozen[i])
{
- dx[point][i] = ff[i]; /* Initial search direction */
+ dx[point][i] = fInit[i]; /* Initial search direction */
}
else
{
// (the main efficiency in the algorithm comes from changing directions), but
// we still need an initial value, so estimate it as the inverse of the norm
// so we take small steps where the potential fluctuates a lot.
- stepsize = 1.0/fnorm;
+ stepsize = 1.0/ems.fnorm;
/* Start the loop over BFGS steps.
* Each successful step is counted, and we continue until
}
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
- top_global, step, (real)step, state_global, state_global, f);
+ top_global, step, (real)step, &ems.s, state_global, &ems.f);
/* Do the linesearching in the direction dx[point][0..(n-1)] */
/* make s a pointer to current search direction - point=0 first time we get here */
s = dx[point];
+ real *xx = static_cast<real *>(as_rvec_array(ems.s.x.data())[0]);
+ real *ff = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
+
// calculate line gradient in position A
for (gpa = 0, i = 0; i < n; i++)
{
}
// Before taking any steps along the line, store the old position
- for (i = 0; i < n; i++)
- {
- lastx[i] = xx[i];
- lastf[i] = ff[i];
- }
- Epot0 = Epot;
+ *last = ems;
+ real *lastx = static_cast<real *>(as_rvec_array(last->s.x.data())[0]);
+ real *lastf = static_cast<real *>(as_rvec_array(last->f.data())[0]);
+ Epot0 = ems.epot;
- for (i = 0; i < n; i++)
- {
- xa[i] = xx[i];
- }
+ *sa = ems;
/* Take a step downhill.
* In theory, we should find the actual minimum of the function in this
// State "A" is the first position along the line.
// reference position along line is initially zero
- EpotA = Epot0;
a = 0.0;
// Check stepsize first. We do not allow displacements
while (maxdelta > inputrec->em_stepsize);
// Take a trial step and move the coordinate array xc[] to position C
+ real *xc = static_cast<real *>(as_rvec_array(sc->s.x.data())[0]);
for (i = 0; i < n; i++)
{
xc[i] = lastx[i] + c*s[i];
neval++;
// Calculate energy for the trial step in position C
- ems.s.x = (rvec *)xc;
- ems.f = (rvec *)fc;
evaluate_energy(fplog, cr,
- top_global, &ems, top,
+ top_global, sc, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, step, FALSE);
- EpotC = ems.epot;
// Calc line gradient in position C
+ real *fc = static_cast<real *>(as_rvec_array(sc->f.data())[0]);
for (gpc = 0, i = 0; i < n; i++)
{
gpc -= s[i]*fc[i]; /* f is negative gradient, thus the sign */
// This is the max amount of increase in energy we tolerate.
// By allowing VERY small changes (close to numerical precision) we
// frequently find even better (lower) final energies.
- tmp = sqrt(GMX_REAL_EPS)*fabs(EpotA);
+ tmp = sqrt(GMX_REAL_EPS)*fabs(sa->epot);
// Accept the step if the energy is lower in the new position C (compared to A),
// or if it is not significantly higher and the line derivative is still negative.
- if (EpotC < EpotA || (gpc < 0 && EpotC < (EpotA+tmp)))
+ if (sc->epot < sa->epot || (gpc < 0 && sc->epot < (sa->epot + tmp)))
{
// Great, we found a better energy. We no longer try to alter the
// stepsize, but simply accept this new better position. The we select a new
// to take smaller steps along a line. Set fnorm based on the new C position,
// which will be used to update the stepsize to 1/fnorm further down.
foundlower = TRUE;
- fnorm = ems.fnorm;
}
else
{
// I also have a safeguard for potentially really pathological functions so we never
// take more than 20 steps before we give up.
// If we already found a lower value we just skip this step and continue to the update.
- nminstep = 0;
+ real fnorm = 0;
+ nminstep = 0;
do
{
// Select a new trial point B in the interval [A,C].
}
// Take a trial step to point B
+ real *xb = static_cast<real *>(as_rvec_array(sb->s.x.data())[0]);
for (i = 0; i < n; i++)
{
xb[i] = lastx[i] + b*s[i];
neval++;
// Calculate energy for the trial step in point B
- ems.s.x = (rvec *)xb;
- ems.f = (rvec *)fb;
evaluate_energy(fplog, cr,
- top_global, &ems, top,
+ top_global, sb, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, step, FALSE);
- EpotB = ems.epot;
- fnorm = ems.fnorm;
+ fnorm = sb->fnorm;
// Calculate gradient in point B
+ real *fb = static_cast<real *>(as_rvec_array(sb->f.data())[0]);
for (gpb = 0, i = 0; i < n; i++)
{
gpb -= s[i]*fb[i]; /* f is negative gradient, thus the sign */
if (gpb > 0)
{
/* Replace c endpoint with b */
- EpotC = EpotB;
- c = b;
- gpc = gpb;
- /* swap coord pointers b/c */
- xtmp = xb;
- ftmp = fb;
- xb = xc;
- fb = fc;
- xc = xtmp;
- fc = ftmp;
+ c = b;
+ /* swap states b and c */
+ swap_em_state(&sb, &sc);
}
else
{
/* Replace a endpoint with b */
- EpotA = EpotB;
- a = b;
- gpa = gpb;
- /* swap coord pointers a/b */
- xtmp = xb;
- ftmp = fb;
- xb = xa;
- fb = fa;
- xa = xtmp;
- fa = ftmp;
+ a = b;
+ /* swap states a and b */
+ swap_em_state(&sa, &sb);
}
/*
*/
nminstep++;
}
- while ((EpotB > EpotA || EpotB > EpotC) && (nminstep < 20));
+ while ((sb->epot > sa->epot || sb->epot > sc->epot) && (nminstep < 20));
- if (fabs(EpotB-Epot0) < GMX_REAL_EPS || nminstep >= 20)
+ if (fabs(sb->epot - Epot0) < GMX_REAL_EPS || nminstep >= 20)
{
/* OK. We couldn't find a significantly lower energy.
* If ncorr==0 this was steepest descent, and then we give up.
/* Select min energy state of A & C, put the best in xx/ff/Epot
*/
- if (EpotC < EpotA)
+ if (sc->epot < sa->epot)
{
- Epot = EpotC;
/* Use state C */
- for (i = 0; i < n; i++)
- {
- xx[i] = xc[i];
- ff[i] = fc[i];
- }
+ ems = *sc;
step_taken = c;
}
else
{
- Epot = EpotA;
/* Use state A */
- for (i = 0; i < n; i++)
- {
- xx[i] = xa[i];
- ff[i] = fa[i];
- }
+ ems = *sa;
step_taken = a;
}
else
{
/* found lower */
- Epot = EpotC;
/* Use state C */
- for (i = 0; i < n; i++)
- {
- xx[i] = xc[i];
- ff[i] = fc[i];
- }
+ ems = *sc;
step_taken = c;
}
}
}
- /* Test whether the convergence criterion is met */
- get_f_norm_max(cr, &(inputrec->opts), mdatoms, f, &fnorm, &fmax, &nfmax);
-
/* Print it if necessary */
if (MASTER(cr))
{
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
fprintf(stderr, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
- step, Epot, fnorm/sqrtNumAtoms, fmax, nfmax+1);
+ step, ems.epot, ems.fnorm/sqrtNumAtoms, ems.fmax, ems.a_fmax + 1);
fflush(stderr);
}
/* Store the new (lower) energies */
}
/* Send x and E to IMD client, if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
// Reset stepsize in we are doing more iterations
- stepsize = 1.0/fnorm;
+ stepsize = 1.0/ems.fnorm;
/* Stop when the maximum force lies below tolerance.
* If we have reached machine precision, converged is already set to true.
*/
- converged = converged || (fmax < inputrec->em_tol);
+ converged = converged || (ems.fmax < inputrec->em_tol);
} /* End of the loop */
step--; /* we never took that last step in this case */
}
- if (fmax > inputrec->em_tol)
+ if (ems.fmax > inputrec->em_tol)
{
if (MASTER(cr))
{
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
print_converged(stderr, LBFGS, inputrec->em_tol, step, converged,
- number_steps, Epot, fmax, nfmax, fnorm/sqrtNumAtoms);
+ number_steps, &ems, sqrtNumAtoms);
print_converged(fplog, LBFGS, inputrec->em_tol, step, converged,
- number_steps, Epot, fmax, nfmax, fnorm/sqrtNumAtoms);
+ number_steps, &ems, sqrtNumAtoms);
fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
}
gmx_walltime_accounting_t walltime_accounting)
{
const char *SD = "Steepest Descents";
- em_state_t *s_min, *s_try;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
real stepsize;
- real ustep, fnormn;
+ real ustep;
gmx_mdoutf_t outf;
t_mdebin *mdebin;
gmx_bool bDone, bAbort, do_x, do_f;
int count = 0;
int steps_accepted = 0;
- s_min = init_em_state();
- s_try = init_em_state();
+ /* Create 2 states on the stack and extract pointers that we will swap */
+ em_state_t s0 {}, s1 {};
+ em_state_t *s_min = &s0;
+ em_state_t *s_try = &s1;
/* Init em and store the local state in s_try */
init_em(fplog, SD, cr, inputrec,
- state_global, top_global, s_try, &top, &f,
+ state_global, top_global, s_try, &top,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
vsite, constr, NULL,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
{
validStep =
do_em_step(cr, inputrec, mdatoms, fr->bMolPBC,
- s_min, stepsize, s_min->f, s_try,
+ s_min, stepsize, &s_min->f, s_try,
constr, top, nrnb, wcycle, count);
}
/* Copy the arrays for force, positions and energy */
/* The 'Min' array always holds the coords and forces of the minimal
sampled energy */
- swap_em_state(s_min, s_try);
+ swap_em_state(&s_min, &s_try);
if (count > 0)
{
ustep *= 1.2;
}
/* Send IMD energies and positions, if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
if (MASTER(cr))
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
- fnormn = s_min->fnorm/sqrtNumAtoms;
print_converged(stderr, SD, inputrec->em_tol, count, bDone, nsteps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
print_converged(fplog, SD, inputrec->em_tol, count, bDone, nsteps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
}
finish_em(cr, outf, walltime_accounting, wcycle);
int nnodes, node;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
tensor vir, pres;
size_t sz;
gmx_sparsematrix_t * sparse_matrix = NULL;
real * full_matrix = NULL;
- em_state_t * state_work;
/* added with respect to mdrun */
int row, col;
gmx_fatal(FARGS, "Constraints present with Normal Mode Analysis, this combination is not supported");
}
- state_work = init_em_state();
-
gmx_shellfc_t *shellfc;
+ em_state_t state_work {};
+
/* Init em and store the local state in state_minimum */
init_em(fplog, NM, cr, inputrec,
- state_global, top_global, state_work, &top,
- &f,
+ state_global, top_global, &state_work, &top,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
vsite, constr, &shellfc,
nfile, fnm, &outf, NULL, imdport, Flags, wcycle);
/* Make evaluate_energy do a single node force calculation */
cr->nnodes = 1;
evaluate_energy(fplog, cr,
- top_global, state_work, top,
+ top_global, &state_work, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, -1, TRUE);
cr->nnodes = nnodes;
/* if forces are not small, warn user */
- get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, state_work);
+ get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, &state_work);
- GMX_LOG(mdlog.warning).appendTextFormatted("Maximum force:%12.5e", state_work->fmax);
- if (state_work->fmax > 1.0e-3)
+ GMX_LOG(mdlog.warning).appendTextFormatted("Maximum force:%12.5e", state_work.fmax);
+ if (state_work.fmax > 1.0e-3)
{
GMX_LOG(mdlog.warning).appendText(
"The force is probably not small enough to "
int force_flags = GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES;
double t = 0;
- x_min = state_work->s.x[atom][d];
+ x_min = state_work.s.x[atom][d];
for (unsigned int dx = 0; (dx < 2); dx++)
{
if (dx == 0)
{
- state_work->s.x[atom][d] = x_min - der_range;
+ state_work.s.x[atom][d] = x_min - der_range;
}
else
{
- state_work->s.x[atom][d] = x_min + der_range;
+ state_work.s.x[atom][d] = x_min + der_range;
}
/* Make evaluate_energy do a single node force calculation */
inputrec, bNS, force_flags,
top,
constr, enerd, fcd,
- &state_work->s, state_work->f, vir, mdatoms,
+ &state_work.s, &state_work.f, vir, mdatoms,
nrnb, wcycle, graph, &top_global->groups,
shellfc, fr, bBornRadii, t, mu_tot,
vsite, NULL);
else
{
evaluate_energy(fplog, cr,
- top_global, state_work, top,
+ top_global, &state_work, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, atom*2+dx, FALSE);
{
for (size_t i = 0; i < atom_index.size(); i++)
{
- copy_rvec(state_work->f[atom_index[i]], fneg[i]);
+ copy_rvec(state_work.f[atom_index[i]], fneg[i]);
}
}
}
/* x is restored to original */
- state_work->s.x[atom][d] = x_min;
+ state_work.s.x[atom][d] = x_min;
for (size_t j = 0; j < atom_index.size(); j++)
{
for (size_t k = 0; (k < DIM); k++)
{
dfdx[j][k] =
- -(state_work->f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
+ -(state_work.f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
}
}
} t_shell;
struct gmx_shellfc_t {
+ /* Shell counts, indices, parameters and working data */
int nshell_gl; /* The number of shells in the system */
t_shell *shell_gl; /* All the shells (for DD only) */
int *shell_index_gl; /* Global shell index (for DD only) */
gmx_bool bPredict; /* Predict shell positions */
gmx_bool bRequireInit; /* Require initialization of shell positions */
int nflexcon; /* The number of flexible constraints */
- rvec *x[2]; /* Array for iterative minimization */
- rvec *f[2]; /* Array for iterative minimization */
- int x_nalloc; /* The allocation size of x and f */
+
+ /* Temporary arrays, should be fixed size 2 when fully converted to C++ */
+ PaddedRVecVector *x; /* Array for iterative minimization */
+ PaddedRVecVector *f; /* Array for iterative minimization */
+
+ /* Flexible constraint working data */
rvec *acc_dir; /* Acceleration direction for flexcon */
rvec *x_old; /* Old coordinates for flexcon */
int flex_nalloc; /* The allocation size of acc_dir and x_old */
}
snew(shfc, 1);
+ shfc->x = new PaddedRVecVector[2] {};
+ shfc->f = new PaddedRVecVector[2] {};
shfc->nflexcon = nflexcon;
if (nshell == 0)
shfc->shell = shell;
}
-static void do_1pos(rvec xnew, rvec xold, rvec f, real step)
+static void do_1pos(rvec xnew, const rvec xold, const rvec f, real step)
{
real xo, yo, zo;
real dx, dy, dz;
xnew[ZZ] = zo+dz;
}
-static void do_1pos3(rvec xnew, rvec xold, rvec f, rvec step)
+static void do_1pos3(rvec xnew, const rvec xold, const rvec f, const rvec step)
{
real xo, yo, zo;
real dx, dy, dz;
xnew[ZZ] = zo+dz;
}
-static void directional_sd(rvec xold[], rvec xnew[], rvec acc_dir[],
- int start, int homenr, real step)
+static void directional_sd(const PaddedRVecVector *xold, PaddedRVecVector *xnew, const rvec acc_dir[],
+ int homenr, real step)
{
- int i;
+ const rvec *xo = as_rvec_array(xold->data());
+ rvec *xn = as_rvec_array(xnew->data());
- for (i = start; i < homenr; i++)
+ for (int i = 0; i < homenr; i++)
{
- do_1pos(xnew[i], xold[i], acc_dir[i], step);
+ do_1pos(xn[i], xo[i], acc_dir[i], step);
}
}
-static void shell_pos_sd(rvec xcur[], rvec xnew[], rvec f[],
+static void shell_pos_sd(const PaddedRVecVector * gmx_restrict xcur,
+ PaddedRVecVector * gmx_restrict xnew,
+ const PaddedRVecVector *f,
int ns, t_shell s[], int count)
{
const real step_scale_min = 0.8,
{
for (d = 0; d < DIM; d++)
{
- dx = xcur[shell][d] - s[i].xold[d];
- df = f[shell][d] - s[i].fold[d];
+ dx = (*xcur)[shell][d] - s[i].xold[d];
+ df = (*f)[shell][d] - s[i].fold[d];
/* -dx/df gets used to generate an interpolated value, but would
* cause a NaN if df were binary-equal to zero. Values close to
* zero won't cause problems (because of the min() and max()), so
#endif
}
}
- copy_rvec(xcur[shell], s[i].xold);
- copy_rvec(f[shell], s[i].fold);
+ copy_rvec((*xcur)[shell], s[i].xold);
+ copy_rvec((*f)[shell], s[i].fold);
- do_1pos3(xnew[shell], xcur[shell], f[shell], s[i].step);
+ do_1pos3((*xnew)[shell], (*xcur)[shell], (*f)[shell], s[i].step);
if (gmx_debug_at)
{
fprintf(debug, "shell[%d] = %d\n", i, shell);
- pr_rvec(debug, 0, "fshell", f[shell], DIM, TRUE);
- pr_rvec(debug, 0, "xold", xcur[shell], DIM, TRUE);
+ pr_rvec(debug, 0, "fshell", (*f)[shell], DIM, TRUE);
+ pr_rvec(debug, 0, "xold", (*xcur)[shell], DIM, TRUE);
pr_rvec(debug, 0, "step", s[i].step, DIM, TRUE);
- pr_rvec(debug, 0, "xnew", xnew[shell], DIM, TRUE);
+ pr_rvec(debug, 0, "xnew", (*xnew)[shell], DIM, TRUE);
}
}
#ifdef PRINT_STEP
}
-static real rms_force(t_commrec *cr, rvec f[], int ns, t_shell s[],
+static real rms_force(t_commrec *cr, const PaddedRVecVector *force, int ns, t_shell s[],
int ndir, real *sf_dir, real *Epot)
{
- int i, shell, ntot;
- double buf[4];
+ double buf[4];
+ const rvec *f = as_rvec_array(force->data());
buf[0] = *sf_dir;
- for (i = 0; i < ns; i++)
+ for (int i = 0; i < ns; i++)
{
- shell = s[i].shell;
- buf[0] += norm2(f[shell]);
+ int shell = s[i].shell;
+ buf[0] += norm2(f[shell]);
}
- ntot = ns;
+ int ntot = ns;
if (PAR(cr))
{
return (ntot ? std::sqrt(buf[0]/ntot) : 0);
}
-static void check_pbc(FILE *fp, rvec x[], int shell)
+static void check_pbc(FILE *fp, PaddedRVecVector x, int shell)
{
int m, now;
{
if (fabs(x[shell][m]-x[now][m]) > 0.3)
{
- pr_rvecs(fp, 0, "SHELL-X", x+now, 5);
+ pr_rvecs(fp, 0, "SHELL-X", as_rvec_array(x.data())+now, 5);
break;
}
}
}
-static void dump_shells(FILE *fp, rvec x[], rvec f[], real ftol, int ns, t_shell s[])
+static void dump_shells(FILE *fp, PaddedRVecVector x, PaddedRVecVector f, real ftol, int ns, t_shell s[])
{
int i, shell;
real ft2, ff2;
static void init_adir(FILE *log, gmx_shellfc_t *shfc,
gmx_constr_t constr, t_idef *idef, t_inputrec *ir,
t_commrec *cr, int dd_ac1,
- gmx_int64_t step, t_mdatoms *md, int start, int end,
+ gmx_int64_t step, t_mdatoms *md, int end,
rvec *x_old, rvec *x_init, rvec *x,
rvec *f, rvec *acc_dir,
gmx_bool bMolPBC, matrix box,
- real *lambda, real *dvdlambda, t_nrnb *nrnb)
+ const std::vector<real> *lambda, real *dvdlambda,
+ t_nrnb *nrnb)
{
rvec *xnold, *xnew;
double dt, w_dt;
}
else
{
- n = end - start;
+ n = end;
}
if (n > shfc->adir_nalloc)
{
dt = ir->delta_t;
/* Does NOT work with freeze or acceleration groups (yet) */
- for (n = start; n < end; n++)
+ for (n = 0; n < end; n++)
{
w_dt = md->invmass[n]*dt;
{
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
{
- xnold[n-start][d] = x[n][d] - (x_init[n][d] - x_old[n][d]);
- xnew[n-start][d] = 2*x[n][d] - x_old[n][d] + f[n][d]*w_dt*dt;
+ 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;
}
else
{
- xnold[n-start][d] = x[n][d];
- xnew[n-start][d] = x[n][d];
+ xnold[n][d] = x[n][d];
+ xnew[n][d] = x[n][d];
}
}
}
constrain(log, FALSE, FALSE, constr, idef, ir, cr, step, 0, 1.0, md,
- x, xnold-start, NULL, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ x, xnold, NULL, bMolPBC, box,
+ (*lambda)[efptBONDED], &(dvdlambda[efptBONDED]),
NULL, NULL, nrnb, econqCoord);
constrain(log, FALSE, FALSE, constr, idef, ir, cr, step, 0, 1.0, md,
- x, xnew-start, NULL, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ x, xnew, NULL, bMolPBC, box,
+ (*lambda)[efptBONDED], &(dvdlambda[efptBONDED]),
NULL, NULL, nrnb, econqCoord);
- for (n = start; n < end; n++)
+ for (n = 0; n < end; n++)
{
for (d = 0; d < DIM; d++)
{
- xnew[n-start][d] =
- -(2*x[n][d]-xnold[n-start][d]-xnew[n-start][d])/gmx::square(dt)
+ xnew[n][d] =
+ -(2*x[n][d]-xnold[n][d]-xnew[n][d])/gmx::square(dt)
- f[n][d]*md->invmass[n];
}
clear_rvec(acc_dir[n]);
/* Project the acceleration on the old bond directions */
constrain(log, FALSE, FALSE, constr, idef, ir, cr, step, 0, 1.0, md,
- x_old, xnew-start, acc_dir, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ x_old, xnew, acc_dir, bMolPBC, box,
+ (*lambda)[efptBONDED], &(dvdlambda[efptBONDED]),
NULL, NULL, nrnb, econqDeriv_FlexCon);
}
gmx_localtop_t *top,
gmx_constr_t constr,
gmx_enerdata_t *enerd, t_fcdata *fcd,
- t_state *state, rvec f[],
+ t_state *state, PaddedRVecVector *f,
tensor force_vir,
t_mdatoms *md,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
int nshell;
t_shell *shell;
t_idef *idef;
- rvec *pos[2], *force[2], *acc_dir = NULL, *x_old = NULL;
+ rvec *acc_dir = NULL, *x_old = NULL;
real Epot[2], df[2];
real sf_dir, invdt;
real ftol, dum = 0;
char sbuf[22];
gmx_bool bCont, bInit, bConverged;
int nat, dd_ac0, dd_ac1 = 0, i;
- int start = 0, homenr = md->homenr, end = start+homenr, cg0, cg1;
+ int homenr = md->homenr, end = homenr, cg0, cg1;
int nflexcon, number_steps, d, Min = 0, count = 0;
#define Try (1-Min) /* At start Try = 1 */
nat = state->natoms;
}
- if (nat > shfc->x_nalloc)
+ for (i = 0; (i < 2); i++)
{
- /* Allocate local arrays */
- shfc->x_nalloc = over_alloc_dd(nat);
- for (i = 0; (i < 2); i++)
- {
- srenew(shfc->x[i], shfc->x_nalloc);
- srenew(shfc->f[i], shfc->x_nalloc);
- }
+ shfc->x[i].resize(nat + 1);
+ shfc->f[i].resize(nat + 1);
}
+
+ /* Create pointer that we can swap */
+ PaddedRVecVector *pos[2];
+ PaddedRVecVector *force[2];
for (i = 0; (i < 2); i++)
{
- pos[i] = shfc->x[i];
- force[i] = shfc->f[i];
+ pos[i] = &shfc->x[i];
+ force[i] = &shfc->f[i];
}
if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr))
*/
if (inputrec->cutoff_scheme == ecutsVERLET)
{
- put_atoms_in_box_omp(fr->ePBC, state->box, md->homenr, state->x);
+ put_atoms_in_box_omp(fr->ePBC, state->box, md->homenr, as_rvec_array(state->x.data()));
}
else
{
cg0 = 0;
cg1 = top->cgs.nr;
put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
- &(top->cgs), state->x, fr->cg_cm);
+ &(top->cgs), as_rvec_array(state->x.data()), fr->cg_cm);
}
if (graph)
{
- mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
+ mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
}
}
/* After this all coordinate arrays will contain whole charge groups */
if (graph)
{
- shift_self(graph, state->box, state->x);
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
if (nflexcon)
for (d = 0; d < DIM; d++)
{
shfc->x_old[i][d] =
- state->x[start+i][d] - state->v[start+i][d]*inputrec->delta_t;
+ state->x[i][d] - state->v[i][d]*inputrec->delta_t;
}
}
}
/* Do a prediction of the shell positions */
if (shfc->bPredict && !bCont)
{
- predict_shells(fplog, state->x, state->v, inputrec->delta_t, nshell, shell,
+ predict_shells(fplog, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), inputrec->delta_t, nshell, shell,
md->massT, NULL, bInit);
}
/* do_force expected the charge groups to be in the box */
if (graph)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
/* Calculate the forces first time around */
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "x b4 do_force", state->x + start, homenr);
+ pr_rvecs(debug, 0, "x b4 do_force", as_rvec_array(state->x.data()), homenr);
}
do_force(fplog, cr, inputrec, mdstep, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
+ state->box, &state->x, &state->hist,
force[Min], force_vir, md, enerd, fcd,
- state->lambda, graph,
+ &state->lambda, graph,
fr, vsite, mu_tot, t, fp_field, NULL, bBornRadii,
(bDoNS ? GMX_FORCE_NS : 0) | force_flags);
if (nflexcon)
{
init_adir(fplog, shfc,
- constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end,
- shfc->x_old-start, state->x, state->x, force[Min],
- shfc->acc_dir-start,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ constr, idef, inputrec, cr, dd_ac1, mdstep, md, end,
+ shfc->x_old, as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), as_rvec_array(force[Min]->data()),
+ shfc->acc_dir,
+ fr->bMolPBC, state->box, &state->lambda, &dum, nrnb);
- for (i = start; i < end; i++)
+ for (i = 0; i < end; i++)
{
- sf_dir += md->massT[i]*norm2(shfc->acc_dir[i-start]);
+ sf_dir += md->massT[i]*norm2(shfc->acc_dir[i]);
}
}
Epot[Min] = enerd->term[F_EPOT];
- df[Min] = rms_force(cr, shfc->f[Min], nshell, shell, nflexcon, &sf_dir, &Epot[Min]);
+ df[Min] = rms_force(cr, &shfc->f[Min], nshell, shell, nflexcon, &sf_dir, &Epot[Min]);
df[Try] = 0;
if (debug)
{
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "force0", force[Min], md->nr);
+ pr_rvecs(debug, 0, "force0", as_rvec_array(force[Min]->data()), md->nr);
}
if (nshell+nflexcon > 0)
* shell positions are updated, therefore the other particles must
* be set here.
*/
- memcpy(pos[Min], state->x, nat*sizeof(state->x[0]));
- memcpy(pos[Try], state->x, nat*sizeof(state->x[0]));
+ *pos[Min] = state->x;
+ *pos[Try] = state->x;
}
if (bVerbose && MASTER(cr))
{
if (vsite)
{
- construct_vsites(vsite, pos[Min], inputrec->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(pos[Min]->data()),
+ inputrec->delta_t, as_rvec_array(state->v.data()),
idef->iparams, idef->il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
if (nflexcon)
{
init_adir(fplog, shfc,
- constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end,
- x_old-start, state->x, pos[Min], force[Min], acc_dir-start,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ constr, idef, inputrec, cr, dd_ac1, mdstep, md, end,
+ x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Min]->data()), as_rvec_array(force[Min]->data()), acc_dir,
+ fr->bMolPBC, state->box, &state->lambda, &dum, nrnb);
- directional_sd(pos[Min], pos[Try], acc_dir-start, start, end,
- fr->fc_stepsize);
+ directional_sd(pos[Min], pos[Try], acc_dir, end, fr->fc_stepsize);
}
/* New positions, Steepest descent */
/* do_force expected the charge groups to be in the box */
if (graph)
{
- unshift_self(graph, state->box, pos[Try]);
+ unshift_self(graph, state->box, as_rvec_array(pos[Try]->data()));
}
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "RELAX: pos[Min] ", pos[Min] + start, homenr);
- pr_rvecs(debug, 0, "RELAX: pos[Try] ", pos[Try] + start, homenr);
+ pr_rvecs(debug, 0, "RELAX: pos[Min] ", as_rvec_array(pos[Min]->data()), homenr);
+ pr_rvecs(debug, 0, "RELAX: pos[Try] ", as_rvec_array(pos[Try]->data()), homenr);
}
/* Try the new positions */
do_force(fplog, cr, inputrec, 1, nrnb, wcycle,
top, groups, state->box, pos[Try], &state->hist,
force[Try], force_vir,
- md, enerd, fcd, state->lambda, graph,
+ md, enerd, fcd, &state->lambda, graph,
fr, vsite, mu_tot, t, fp_field, NULL, bBornRadii,
force_flags);
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "RELAX: force[Min]", force[Min] + start, homenr);
- pr_rvecs(debug, 0, "RELAX: force[Try]", force[Try] + start, homenr);
+ pr_rvecs(debug, 0, "RELAX: force[Min]", as_rvec_array(force[Min]->data()), homenr);
+ pr_rvecs(debug, 0, "RELAX: force[Try]", as_rvec_array(force[Try]->data()), homenr);
}
sf_dir = 0;
if (nflexcon)
{
init_adir(fplog, shfc,
- constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end,
- x_old-start, state->x, pos[Try], force[Try], acc_dir-start,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ constr, idef, inputrec, cr, dd_ac1, mdstep, md, end,
+ x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Try]->data()), as_rvec_array(force[Try]->data()), acc_dir,
+ fr->bMolPBC, state->box, &state->lambda, &dum, nrnb);
- for (i = start; i < end; i++)
+ for (i = 0; i < end; i++)
{
- sf_dir += md->massT[i]*norm2(acc_dir[i-start]);
+ sf_dir += md->massT[i]*norm2(acc_dir[i]);
}
}
{
if (gmx_debug_at)
{
- pr_rvecs(debug, 0, "F na do_force", force[Try] + start, homenr);
+ pr_rvecs(debug, 0, "F na do_force", as_rvec_array(force[Try]->data()), homenr);
}
if (gmx_debug_at)
{
fprintf(debug, "SHELL ITER %d\n", count);
- dump_shells(debug, pos[Try], force[Try], ftol, nshell, shell);
+ dump_shells(debug, *pos[Try], *force[Try], ftol, nshell, shell);
}
}
{
/* Correct the velocities for the flexible constraints */
invdt = 1/inputrec->delta_t;
- for (i = start; i < end; i++)
+ for (i = 0; i < end; i++)
{
for (d = 0; d < DIM; d++)
{
}
/* Copy back the coordinates and the forces */
- memcpy(state->x, pos[Min], nat*sizeof(state->x[0]));
- memcpy(f, force[Min], nat*sizeof(f[0]));
+ state->x = *pos[Min];
+ *f = *force[Min];
}
void done_shellfc(FILE *fplog, gmx_shellfc_t *shfc, gmx_int64_t numSteps)
#include <cstdio>
#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/timing/wallcycle.h"
struct gmx_constr;
gmx_localtop_t *top,
gmx_constr *constr,
gmx_enerdata_t *enerd, t_fcdata *fcd,
- t_state *state, rvec f[],
+ t_state *state, PaddedRVecVector *f,
tensor force_vir,
t_mdatoms *md,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
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,
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());
+ rvec *f = as_rvec_array(force->data());
+
switch (inputrec->cutoff_scheme)
{
case ecutsVERLET:
f, vir_force,
mdatoms,
enerd, fcd,
- lambda, graph,
+ lambda->data(), graph,
fr, fr->ic,
vsite, mu_tot,
t, field, ed,
f, vir_force,
mdatoms,
enerd, fcd,
- lambda, graph,
+ lambda->data(), graph,
fr, vsite, mu_tot,
t, field, ed,
bBornRadii,
/* 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++)
{
}
}
-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->all_lambda[i][*fep_state];
+ if (lam0)
{
- lambda[i] = fep->init_lambda;
- if (lam0)
- {
- lam0[i] = lambda[i];
- }
- }
- else
- {
- 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[],
matrix box, real lambda, tensor pres, tensor virial,
real *prescorr, real *enercorr, real *dvdlcorr);
-void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, real *lambda, double *lam0);
+void initialize_lambdas(FILE *fplog, t_inputrec *ir, int *fep_state, std::vector<real> *lambda, double *lam0);
void do_constrain_first(FILE *log, gmx_constr *constr,
t_inputrec *inputrec, t_mdatoms *md,
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[],
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
- gmx_localtop_t *top;
- gmx_groups_t *groups;
- gmx_enerdata_t *enerd;
- rvec *f;
- real lambda, t, temp, beta, drmax, epot;
- double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
- t_trxstatus *status;
- t_trxframe rerun_fr;
- gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
- tensor force_vir, shake_vir, vir, pres;
- int cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
- rvec *x_mol;
- rvec mu_tot, x_init, dx, x_tp;
- int nnodes, frame;
- gmx_int64_t frame_step_prev, frame_step;
- gmx_int64_t nsteps, stepblocksize = 0, step;
- gmx_int64_t seed;
- int i;
- FILE *fp_tpi = NULL;
- char *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
- double dbl, dump_ener;
- gmx_bool bCavity;
- int nat_cavity = 0, d;
- real *mass_cavity = NULL, mass_tot;
- int nbin;
- double invbinw, *bin, refvolshift, logV, bUlogV;
- real prescorr, enercorr, dvdlcorr;
- gmx_bool bEnergyOutOfBounds;
- const char *tpid_leg[2] = {"direct", "reweighted"};
+ gmx_localtop_t *top;
+ gmx_groups_t *groups;
+ gmx_enerdata_t *enerd;
+ PaddedRVecVector f {};
+ real lambda, t, temp, beta, drmax, epot;
+ double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
+ t_trxstatus *status;
+ t_trxframe rerun_fr;
+ gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
+ tensor force_vir, shake_vir, vir, pres;
+ int cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
+ rvec *x_mol;
+ rvec mu_tot, x_init, dx, x_tp;
+ int nnodes, frame;
+ gmx_int64_t frame_step_prev, frame_step;
+ gmx_int64_t nsteps, stepblocksize = 0, step;
+ gmx_int64_t seed;
+ int i;
+ FILE *fp_tpi = NULL;
+ char *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
+ double dbl, dump_ener;
+ gmx_bool bCavity;
+ int nat_cavity = 0, d;
+ real *mass_cavity = NULL, mass_tot;
+ int nbin;
+ double invbinw, *bin, refvolshift, logV, bUlogV;
+ real prescorr, enercorr, dvdlcorr;
+ gmx_bool bEnergyOutOfBounds;
+ const char *tpid_leg[2] = {"direct", "reweighted"};
/* Since there is no upper limit to the insertion energies,
* we need to set an upper limit for the distribution output.
snew(enerd, 1);
init_enerdata(groups->grps[egcENER].nr, inputrec->fepvals->n_lambda, enerd);
- snew(f, top_global->natoms);
+ f.resize(top_global->natoms + 1);
/* Print to log file */
walltime_accounting_start(walltime_accounting);
}
bRFExcl = (bCharge && EEL_RF(fr->eeltype));
- calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), state_global->x, fr->cg_cm);
+ calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), as_rvec_array(state_global->x.data()), fr->cg_cm);
if (bCavity)
{
if (norm(fr->cg_cm[cg_tp]) > 0.5*inputrec->rlist && fplog)
copy_rvec(x_mol[i-a_tp0], state_global->x[i]);
}
/* Rotate the molecule randomly */
- rotate_conf(a_tp1-a_tp0, state_global->x+a_tp0, NULL,
+ rotate_conf(a_tp1-a_tp0, as_rvec_array(state_global->x.data())+a_tp0, NULL,
2*M_PI*dist(rng),
2*M_PI*dist(rng),
2*M_PI*dist(rng));
cr->nnodes = 1;
do_force(fplog, cr, inputrec,
step, nrnb, wcycle, top, &top_global->groups,
- state_global->box, state_global->x, &state_global->hist,
- f, force_vir, mdatoms, enerd, fcd,
- state_global->lambda,
+ state_global->box, &state_global->x, &state_global->hist,
+ &f, force_vir, mdatoms, enerd, fcd,
+ &state_global->lambda,
NULL, fr, NULL, mu_tot, t, NULL, NULL, FALSE,
GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
(bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS : 0) |
{
sprintf(str, "t%g_step%d.pdb", t, (int)step);
sprintf(str2, "t: %f step %d ener: %f", t, (int)step, epot);
- write_sto_conf_mtop(str, str2, top_global, state_global->x, state_global->v,
+ write_sto_conf_mtop(str, str2, top_global, as_rvec_array(state_global->x.data()), as_rvec_array(state_global->v.data()),
inputrec->ePBC, state_global->box);
}
#include "gromacs/utility/smalloc.h"
void
-do_md_trajectory_writing(FILE *fplog,
- t_commrec *cr,
- int nfile,
- const t_filenm fnm[],
- gmx_int64_t step,
- gmx_int64_t step_rel,
- double t,
- t_inputrec *ir,
- t_state *state,
- t_state *state_global,
- gmx_mtop_t *top_global,
- t_forcerec *fr,
- gmx_mdoutf_t outf,
- t_mdebin *mdebin,
- gmx_ekindata_t *ekind,
- rvec *f,
- int *nchkpt,
- gmx_bool bCPT,
- gmx_bool bRerunMD,
- gmx_bool bLastStep,
- gmx_bool bDoConfOut,
- gmx_bool bSumEkinhOld
+do_md_trajectory_writing(FILE *fplog,
+ t_commrec *cr,
+ int nfile,
+ const t_filenm fnm[],
+ gmx_int64_t step,
+ gmx_int64_t step_rel,
+ double t,
+ t_inputrec *ir,
+ t_state *state,
+ t_state *state_global,
+ gmx_mtop_t *top_global,
+ t_forcerec *fr,
+ gmx_mdoutf_t outf,
+ t_mdebin *mdebin,
+ gmx_ekindata_t *ekind,
+ PaddedRVecVector *f,
+ int *nchkpt,
+ gmx_bool bCPT,
+ gmx_bool bRerunMD,
+ gmx_bool bLastStep,
+ gmx_bool bDoConfOut,
+ gmx_bool bSumEkinhOld
)
{
int mdof_flags;
bDoConfOut && MASTER(cr) &&
!bRerunMD)
{
- if (fr->bMolPBC && state->x == state_global->x)
+ if (fr->bMolPBC && state == state_global)
{
/* This (single-rank) run needs to allocate a
temporary array of size natoms so that any
identical, and makes .edr restarts binary
identical. */
snew(x_for_confout, state_global->natoms);
- copy_rvecn(state_global->x, x_for_confout, 0, state_global->natoms);
+ copy_rvecn(as_rvec_array(state_global->x.data()), x_for_confout, 0, state_global->natoms);
}
else
{
/* With DD, or no bMolPBC, it doesn't matter if
- we change state_global->x */
- x_for_confout = state_global->x;
+ we change as_rvec_array(state_global->x.data()) */
+ x_for_confout = as_rvec_array(state_global->x.data());
}
/* x and v have been collected in mdoutf_write_to_trajectory_files,
}
write_sto_conf_mtop(ftp2fn(efSTO, nfile, fnm),
*top_global->name, top_global,
- x_for_confout, state_global->v,
+ x_for_confout, as_rvec_array(state_global->v.data()),
ir->ePBC, state->box);
- if (fr->bMolPBC && state->x == state_global->x)
+ if (fr->bMolPBC && state == state_global)
{
sfree(x_for_confout);
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
gmx_mdoutf_t outf,
t_mdebin *mdebin,
struct gmx_ekindata_t *ekind,
- rvec *f,
+ PaddedRVecVector *f,
int *nchkpt,
gmx_bool bCPT,
gmx_bool bRerunMD,
struct gmx_update_t
{
- gmx_stochd_t *sd;
+ gmx_stochd_t *sd;
/* xprime for constraint algorithms */
- rvec *xp;
- int xp_nalloc;
+ PaddedRVecVector xp;
/* Variables for the deform algorithm */
- gmx_int64_t deformref_step;
- matrix deformref_box;
+ gmx_int64_t deformref_step;
+ matrix deformref_box;
};
real invmass[],
unsigned short ptype[], unsigned short cFREEZE[],
unsigned short cACC[], unsigned short cTC[],
- rvec x[], rvec xprime[], rvec v[],
- rvec f[], matrix M,
+ const rvec x[], rvec xprime[], rvec v[],
+ const rvec f[], matrix M,
gmx_bool bNH, gmx_bool bPR)
{
double imass, w_dt;
static void do_update_vv_vel(int start, int nrend, double dt,
rvec accel[], ivec nFreeze[], real invmass[],
unsigned short ptype[], unsigned short cFREEZE[],
- unsigned short cACC[], rvec v[], rvec f[],
+ unsigned short cACC[], rvec v[], const rvec f[],
gmx_bool bExtended, real veta, real alpha)
{
double w_dt;
static void do_update_vv_pos(int start, int nrend, double dt,
ivec nFreeze[],
unsigned short ptype[], unsigned short cFREEZE[],
- rvec x[], rvec xprime[], rvec v[],
+ const rvec x[], rvec xprime[], rvec v[],
gmx_bool bExtended, real veta)
{
int gf = 0;
double nh_vxi[],
real invmass[],
unsigned short ptype[], unsigned short cTC[],
- rvec x[], rvec xprime[], rvec v[],
- rvec f[], matrix M, matrix box, real
+ const rvec x[], rvec xprime[], rvec v[],
+ const rvec f[], matrix M, matrix box, real
cos_accel, real vcos,
gmx_bool bNH, gmx_bool bPR)
{
gmx_update_t *init_update(const t_inputrec *ir)
{
- gmx_update_t *upd;
-
- snew(upd, 1);
+ gmx_update_t *upd = new(gmx_update_t);
if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
{
update_temperature_constants(upd, ir);
- upd->xp = NULL;
- upd->xp_nalloc = 0;
+ upd->xp.resize(0);
return upd;
}
-void update_realloc(gmx_update_t *upd, int state_nalloc)
+void update_realloc(gmx_update_t *upd, int natoms)
{
GMX_ASSERT(upd, "upd must be allocated before its fields can be reallocated");
- if (state_nalloc > upd->xp_nalloc)
- {
- upd->xp_nalloc = state_nalloc;
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries. */
- srenew(upd->xp, upd->xp_nalloc + 1);
- }
+
+ /* We need to allocate one element extra, since we might use
+ * (unaligned) 4-wide SIMD loads to access rvec entries. */
+ upd->xp.resize(natoms + 1);
}
static void do_update_sd1(gmx_stochd_t *sd,
real invmass[], unsigned short ptype[],
unsigned short cFREEZE[], unsigned short cACC[],
unsigned short cTC[],
- rvec x[], rvec xprime[], rvec v[], rvec f[],
+ const rvec x[], rvec xprime[], rvec v[], const rvec f[],
gmx_bool bDoConstr,
gmx_bool bFirstHalfConstr,
gmx_int64_t step, int seed, int* gatindex)
ivec nFreeze[],
real invmass[], unsigned short ptype[],
unsigned short cFREEZE[], unsigned short cTC[],
- rvec x[], rvec xprime[], rvec v[],
- rvec f[], real friction_coefficient,
+ const rvec x[], rvec xprime[], rvec v[],
+ const rvec f[], real friction_coefficient,
real *rf, gmx_int64_t step, int seed,
int* gatindex)
{
}
static void dump_it_all(FILE gmx_unused *fp, const char gmx_unused *title,
- int gmx_unused natoms, rvec gmx_unused x[], rvec gmx_unused xp[],
- rvec gmx_unused v[], rvec gmx_unused f[])
+ int gmx_unused natoms,
+ const gmx_unused PaddedRVecVector *x,
+ const gmx_unused PaddedRVecVector *xp,
+ const gmx_unused PaddedRVecVector *v,
+ const gmx_unused PaddedRVecVector *f)
{
#ifdef DEBUG
if (fp)
{
fprintf(fp, "%s\n", title);
- pr_rvecs(fp, 0, "x", x, natoms);
- pr_rvecs(fp, 0, "xp", xp, natoms);
- pr_rvecs(fp, 0, "v", v, natoms);
- pr_rvecs(fp, 0, "f", f, natoms);
+ pr_rvecs(fp, 0, "x", as_rvec_array(x->data()), natoms);
+ pr_rvecs(fp, 0, "xp", as_rvec_array(xp->data()), natoms);
+ pr_rvecs(fp, 0, "v", as_rvec_array(v->data()), natoms);
+ if (f != NULL)
+ {
+ pr_rvecs(fp, 0, "f", as_rvec_array(f->data()), natoms);
+ }
}
#endif
}
{
if (ekind->cosacc.cos_accel == 0)
{
- calc_ke_part_normal(state->v, opts, md, ekind, nrnb, bEkinAveVel);
+ calc_ke_part_normal(as_rvec_array(state->v.data()), opts, md, ekind, nrnb, bEkinAveVel);
}
else
{
- calc_ke_part_visc(state->box, state->x, state->v, opts, md, ekind, nrnb, bEkinAveVel);
+ calc_ke_part_visc(state->box, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), opts, md, ekind, nrnb, bEkinAveVel);
}
}
snew(ekinstate->ekinh, ekinstate->ekin_n);
snew(ekinstate->ekinf, ekinstate->ekin_n);
snew(ekinstate->ekinh_old, ekinstate->ekin_n);
- snew(ekinstate->ekinscalef_nhc, ekinstate->ekin_n);
- snew(ekinstate->ekinscaleh_nhc, ekinstate->ekin_n);
- snew(ekinstate->vscale_nhc, ekinstate->ekin_n);
+ ekinstate->ekinscalef_nhc.resize(ekinstate->ekin_n);
+ ekinstate->ekinscaleh_nhc.resize(ekinstate->ekin_n);
+ ekinstate->vscale_nhc.resize(ekinstate->ekin_n);
ekinstate->dekindl = 0;
ekinstate->mvcos = 0;
}
break;
case etcNOSEHOOVER:
nosehoover_tcoupl(&(inputrec->opts), ekind, dttc,
- state->nosehoover_xi, state->nosehoover_vxi, MassQ);
+ state->nosehoover_xi.data(), state->nosehoover_vxi.data(), MassQ);
break;
case etcVRESCALE:
vrescale_tcoupl(inputrec, step, ekind, dttc,
- state->therm_integral);
+ state->therm_integral.data());
break;
}
/* rescale in place here */
if (EI_VV(inputrec->eI))
{
- rescale_velocities(ekind, md, 0, md->homenr, state->v);
+ rescale_velocities(ekind, md, 0, md->homenr, as_rvec_array(state->v.data()));
}
}
else
t_state *state,
gmx_bool bMolPBC,
t_graph *graph,
- rvec force[], /* forces on home particles */
+ PaddedRVecVector *force, /* forces on home particles */
t_idef *idef,
tensor vir_part,
t_commrec *cr,
{
constrain(NULL, bLog, bEner, constr, idef,
inputrec, cr, step, 1, 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()),
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
NULL, bCalcVir ? &vir_con : NULL, nrnb, econqVeloc);
{
constrain(NULL, bLog, bEner, constr, idef,
inputrec, cr, step, 1, 1.0, md,
- state->x, upd->xp, NULL,
+ as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), NULL,
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
- state->v, bCalcVir ? &vir_con : NULL, nrnb, econqCoord);
+ as_rvec_array(state->v.data()), bCalcVir ? &vir_con : NULL, nrnb, econqCoord);
}
wallcycle_stop(wcycle, ewcCONSTR);
where();
dump_it_all(fplog, "After Shake",
- state->natoms, state->x, upd->xp, state->v, force);
+ state->natoms, &state->x, &upd->xp, &state->v, force);
if (bCalcVir)
{
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, upd->xp, state->v, force,
+ as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), as_rvec_array(state->v.data()), as_rvec_array(force->data()),
bDoConstr, FALSE,
step, inputrec->ld_seed,
DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
constrain(NULL, bLog, bEner, constr, idef,
inputrec, cr, step, 1, 0.5, md,
- state->x, upd->xp, NULL,
+ as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), NULL,
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
- state->v, NULL, nrnb, econqCoord);
+ as_rvec_array(state->v.data()), NULL, nrnb, econqCoord);
wallcycle_stop(wcycle, ewcCONSTR);
}
if (graph && (graph->nnodes > 0))
{
- unshift_x(graph, state->box, state->x, upd->xp);
+ unshift_x(graph, state->box, as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()));
if (TRICLINIC(state->box))
{
inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
}
else
{
+ /* The copy is performance sensitive, so use a bare pointer */
+ rvec *xp = as_rvec_array(upd->xp.data());
#ifndef __clang_analyzer__
// cppcheck-suppress unreadVariable
nth = gmx_omp_nthreads_get(emntUpdate);
for (i = start; i < nrend; i++)
{
// Trivial statement, does not throw
- copy_rvec(upd->xp[i], state->x[i]);
+ copy_rvec(xp[i], state->x[i]);
}
}
wallcycle_stop(wcycle, ewcUPDATE);
dump_it_all(fplog, "After unshift",
- state->natoms, state->x, upd->xp, state->v, force);
+ state->natoms, &state->x, &upd->xp, &state->v, force);
}
/* ############# END the update of velocities and positions ######### */
}
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- rvec force[], /* forces on home particles */
matrix pcoupl_mu,
t_nrnb *nrnb,
gmx_update_t *upd)
if (inputrec->nstpcouple == 1 || (step % inputrec->nstpcouple == 1))
{
berendsen_pscale(inputrec, pcoupl_mu, state->box, state->box_rel,
- start, homenr, state->x, md->cFREEZE, nrnb);
+ start, homenr, as_rvec_array(state->x.data()), md->cFREEZE, nrnb);
}
break;
case (epcPARRINELLORAHMAN):
if (inputrecDeform(inputrec))
{
- deform(upd, start, homenr, state->x, state->box, inputrec, step);
+ deform(upd, start, homenr, as_rvec_array(state->x.data()), state->box, inputrec, step);
}
where();
dump_it_all(fplog, "After update",
- state->natoms, state->x, upd->xp, state->v, force);
+ state->natoms, &state->x, &upd->xp, &state->v, NULL);
}
void update_coords(FILE *fplog,
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- rvec *f, /* forces on home particles */
+ PaddedRVecVector *f, /* forces on home particles */
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
/* ############# START The update of velocities and positions ######### */
where();
dump_it_all(fplog, "Before update",
- state->natoms, state->x, upd->xp, state->v, f);
+ state->natoms, &state->x, &upd->xp, &state->v, f);
nth = gmx_omp_nthreads_get(emntUpdate);
start_th = start + ((nrend-start)* th )/nth;
end_th = start + ((nrend-start)*(th+1))/nth;
+ const rvec *x_rvec = as_rvec_array(state->x.data());
+ rvec *xp_rvec = as_rvec_array(upd->xp.data());
+ rvec *v_rvec = as_rvec_array(state->v.data());
+ const rvec *f_rvec = as_rvec_array(f->data());
+
switch (inputrec->eI)
{
case (eiMD):
{
do_update_md(start_th, end_th,
dt, inputrec->nstpcouple,
- ekind->tcstat, state->nosehoover_vxi,
+ ekind->tcstat, state->nosehoover_vxi.data(),
ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, upd->xp, state->v, f, M,
+ x_rvec, xp_rvec, v_rvec, f_rvec, M,
bNH, bPR);
}
else
{
do_update_visc(start_th, end_th,
dt, inputrec->nstpcouple,
- ekind->tcstat, state->nosehoover_vxi,
+ ekind->tcstat, state->nosehoover_vxi.data(),
md->invmass, md->ptype,
- md->cTC, state->x, upd->xp, state->v, f, M,
+ md->cTC, x_rvec, xp_rvec, v_rvec, f_rvec, M,
state->box,
ekind->cosacc.cos_accel,
ekind->cosacc.vcos,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
- state->x, upd->xp, state->v, f,
+ x_rvec, xp_rvec, v_rvec, f_rvec,
bDoConstr, TRUE,
step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
break;
do_update_bd(start_th, end_th, dt,
inputrec->opts.nFreeze, md->invmass, md->ptype,
md->cFREEZE, md->cTC,
- state->x, upd->xp, state->v, f,
+ x_rvec, xp_rvec, v_rvec, f_rvec,
inputrec->bd_fric,
upd->sd->bd_rf,
step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC,
- state->v, f,
+ v_rvec, f_rvec,
(bNH || bPR), state->veta, alpha);
break;
case etrtPOSITION:
do_update_vv_pos(start_th, end_th, dt,
inputrec->opts.nFreeze,
md->ptype, md->cFREEZE,
- state->x, upd->xp, state->v,
+ x_rvec, xp_rvec, v_rvec,
(bNH || bPR), state->veta);
break;
}
#define GMX_MDLIB_UPDATE_H
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
/* Update the size of per-atom arrays (e.g. after DD re-partitioning,
which might increase the number of home atoms). */
-void update_realloc(gmx_update_t *upd, int state_nalloc);
+void update_realloc(gmx_update_t *upd, int natoms);
/* Store the box at step step
* as a reference state for simulations with box deformation.
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- rvec *f, /* forces on home particles */
+ PaddedRVecVector *f, /* forces on home particles */
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
t_state *state,
gmx_bool bMolPBC,
t_graph *graph,
- rvec force[], /* forces on home particles */
+ PaddedRVecVector *force, /* forces on home particles */
t_idef *idef,
tensor vir_part,
t_commrec *cr,
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- rvec force[], /* forces on home particles */
matrix pcoupl_mu,
t_nrnb *nrnb,
gmx_update_t *upd);
void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
double xi[], double vxi[], t_extmass *MassQ);
-t_state *init_bufstate(const t_state *template_state);
-
-void destroy_bufstate(t_state *state);
-
void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
gmx_enerdata_t *enerd, t_state *state, tensor vir, t_mdatoms *md,
t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno);
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/utility/compare.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
/* The source code in this file should be thread-safe.
eks->ekinh = NULL;
eks->ekinf = NULL;
eks->ekinh_old = NULL;
- eks->ekinscalef_nhc = NULL;
- eks->ekinscaleh_nhc = NULL;
- eks->vscale_nhc = NULL;
+ eks->ekinscalef_nhc.resize(0);
+ eks->ekinscaleh_nhc.resize(0);
+ eks->vscale_nhc.resize(0);
eks->dekindl = 0;
eks->mvcos = 0;
}
void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
{
- int i, j;
-
state->ngtc = ngtc;
state->nnhpres = nnhpres;
state->nhchainlength = nhchainlength;
- if (state->ngtc > 0)
- {
- snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
- snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
- snew(state->therm_integral, state->ngtc);
- for (i = 0; i < state->ngtc; i++)
- {
- for (j = 0; j < state->nhchainlength; j++)
- {
- state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
- state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
- }
- }
- for (i = 0; i < state->ngtc; i++)
- {
- state->therm_integral[i] = 0.0;
- }
- }
- else
- {
- state->nosehoover_xi = NULL;
- state->nosehoover_vxi = NULL;
- state->therm_integral = NULL;
- }
-
- if (state->nnhpres > 0)
- {
- snew(state->nhpres_xi, state->nhchainlength*nnhpres);
- snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
- for (i = 0; i < nnhpres; i++)
- {
- for (j = 0; j < state->nhchainlength; j++)
- {
- state->nhpres_xi[i*nhchainlength + j] = 0.0;
- state->nhpres_vxi[i*nhchainlength + j] = 0.0;
- }
- }
- }
- else
- {
- state->nhpres_xi = NULL;
- state->nhpres_vxi = NULL;
- }
+ state->nosehoover_xi.resize(state->nhchainlength*state->ngtc, 0);
+ state->nosehoover_vxi.resize(state->nhchainlength*state->ngtc, 0);
+ state->therm_integral.resize(state->ngtc, 0);
+ state->nhpres_xi.resize(state->nhchainlength*nnhpres, 0);
+ state->nhpres_vxi.resize(state->nhchainlength*nnhpres, 0);
}
void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int dfhistNumLambda)
{
- int i;
-
state->natoms = natoms;
state->flags = 0;
state->fep_state = 0;
- state->lambda = 0;
- snew(state->lambda, efptNR);
- for (i = 0; i < efptNR; i++)
- {
- state->lambda[i] = 0;
- }
+ state->lambda.resize(efptNR, 0);
state->veta = 0;
clear_mat(state->box);
clear_mat(state->box_rel);
clear_mat(state->svir_prev);
clear_mat(state->fvir_prev);
init_gtc_state(state, ngtc, nnhpres, nhchainlength);
- state->nalloc = state->natoms;
- if (state->nalloc > 0)
+ if (state->natoms > 0)
{
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- snew(state->x, state->nalloc + 1);
- snew(state->v, state->nalloc + 1);
+ state->x.resize(state->natoms + 1);
+ state->v.resize(state->natoms + 1);
}
else
{
- state->x = NULL;
- state->v = NULL;
+ state->x.resize(0);
+ state->v.resize(0);
}
- state->cg_p = NULL;
+ state->cg_p.resize(0);
zero_history(&state->hist);
zero_ekinstate(&state->ekinstate);
snew(state->enerhist, 1);
state->edsamstate = NULL;
state->ddp_count = 0;
state->ddp_count_cg_gl = 0;
- state->cg_gl = NULL;
- state->cg_gl_nalloc = 0;
-}
-
-void done_state(t_state *state)
-{
- sfree(state->x);
- sfree(state->v);
- sfree(state->cg_p);
- sfree(state->enerhist);
- state->nalloc = 0;
- sfree(state->cg_gl);
- state->cg_gl_nalloc = 0;
- sfree(state->lambda);
- if (state->ngtc > 0)
- {
- sfree(state->nosehoover_xi);
- sfree(state->nosehoover_vxi);
- sfree(state->therm_integral);
- }
+ state->cg_gl.resize(0);
}
void comp_state(const t_state *st1, const t_state *st2,
if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX)))
{
fprintf(stdout, "comparing x\n");
- cmp_rvecs(stdout, "x", st1->natoms, st1->x, st2->x, bRMSD, ftol, abstol);
+ cmp_rvecs(stdout, "x", st1->natoms, as_rvec_array(st1->x.data()), as_rvec_array(st2->x.data()), bRMSD, ftol, abstol);
}
if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV)))
{
fprintf(stdout, "comparing v\n");
- cmp_rvecs(stdout, "v", st1->natoms, st1->v, st2->v, bRMSD, ftol, abstol);
+ cmp_rvecs(stdout, "v", st1->natoms, as_rvec_array(st1->v.data()), as_rvec_array(st2->v.data()), bRMSD, ftol, abstol);
}
}
}
+
+rvec *getRvecArrayFromPaddedRVecVector(const PaddedRVecVector *v,
+ unsigned int n)
+{
+ GMX_ASSERT(v->size() >= n, "We can't copy more elements than the vector size");
+
+ rvec *dest;
+
+ snew(dest, n);
+
+ const rvec *vPtr = as_rvec_array(v->data());
+ for (unsigned int i = 0; i < n; i++)
+ {
+ copy_rvec(vPtr[i], dest[i]);
+ }
+
+ return dest;
+}
#ifndef GMX_MDTYPES_STATE_H
#define GMX_MDTYPES_STATE_H
+#include <vector>
+
#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/utility/basedefinitions.h"
/* The names of the state entries, defined in src/gmxlib/checkpoint.c */
extern const char *est_names[estNR];
+/* This vector is not padded yet, padding will be added soon */
+typedef std::vector<gmx::RVec> PaddedRVecVector;
+
typedef struct history_t
{
real disre_initf; /* The scaling factor for initializing the time av. */
*/
typedef struct ekinstate_t
{
- gmx_bool bUpToDate;
- int ekin_n;
- tensor *ekinh;
- tensor *ekinf;
- tensor *ekinh_old;
- tensor ekin_total;
- double *ekinscalef_nhc;
- double *ekinscaleh_nhc;
- double *vscale_nhc;
- real dekindl;
- real mvcos;
+ gmx_bool bUpToDate;
+ int ekin_n;
+ tensor *ekinh;
+ tensor *ekinf;
+ tensor *ekinh_old;
+ tensor ekin_total;
+ std::vector<double> ekinscalef_nhc;
+ std::vector<double> ekinscaleh_nhc;
+ std::vector<double> vscale_nhc;
+ real dekindl;
+ real mvcos;
} ekinstate_t;
typedef struct df_history_t
int nhchainlength; /* number of nose-hoover chains */
int flags; /* Flags telling which entries are present */
int fep_state; /* indicates which of the alchemical states we are in */
- real *lambda; /* lambda vector */
+ std::vector<real> lambda; /* lambda vector */
matrix box; /* box vector coordinates */
matrix box_rel; /* Relitaive box vectors to preserve shape */
matrix boxv; /* box velocitites for Parrinello-Rahman pcoupl */
matrix pres_prev; /* Pressure of the previous step for pcoupl */
matrix svir_prev; /* Shake virial for previous step for pcoupl */
matrix fvir_prev; /* Force virial of the previous step for pcoupl */
- double *nosehoover_xi; /* for Nose-Hoover tcoupl (ngtc) */
- double *nosehoover_vxi; /* for N-H tcoupl (ngtc) */
- double *nhpres_xi; /* for Nose-Hoover pcoupl for barostat */
- double *nhpres_vxi; /* for Nose-Hoover pcoupl for barostat */
- double *therm_integral; /* for N-H/V-rescale tcoupl (ngtc) */
+ std::vector<double> nosehoover_xi; /* for Nose-Hoover tcoupl (ngtc) */
+ std::vector<double> nosehoover_vxi; /* for N-H tcoupl (ngtc) */
+ std::vector<double> nhpres_xi; /* for Nose-Hoover pcoupl for barostat */
+ std::vector<double> nhpres_vxi; /* for Nose-Hoover pcoupl for barostat */
+ std::vector<double> therm_integral; /* for N-H/V-rescale tcoupl (ngtc) */
real veta; /* trotter based isotropic P-coupling */
real vol0; /* initial volume,required for computing NPT conserverd quantity */
- int nalloc; /* Allocation size for x and v when !=NULL*/
- rvec *x; /* the coordinates (natoms) */
- rvec *v; /* the velocities (natoms) */
- rvec *cg_p; /* p vector for conjugate gradient minimization */
+ PaddedRVecVector x; /* the coordinates (natoms) */
+ PaddedRVecVector v; /* the velocities (natoms) */
+ PaddedRVecVector cg_p; /* p vector for conjugate gradient minimization */
history_t hist; /* Time history for restraints */
int ddp_count; /* The DD partitioning count for this state */
int ddp_count_cg_gl; /* The DD part. count for index_gl */
- int ncg_gl; /* The number of local charge groups */
- int *cg_gl; /* The global cg number of the local cgs */
- int cg_gl_nalloc; /* Allocation size of cg_gl; */
+ std::vector<int> cg_gl; /* The global cg number of the local cgs */
} t_state;
typedef struct t_extmass
void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int dfhistNumLambda);
-void done_state(t_state *state);
-
void comp_state(const t_state *st1, const t_state *st2, gmx_bool bRMSD, real ftol, real abstol);
+/*! \brief Allocate an rvec pointer and copy the contents of v to it */
+rvec *getRvecArrayFromPaddedRVecVector(const PaddedRVecVector *v,
+ unsigned int n);
+
#endif
{
fprintf(stderr, "Will write subset %s of original tpx containing %d "
"atoms\n", grpname, gnx);
- reduce_topology_x(gnx, index, &mtop, state.x, state.v);
+ reduce_topology_x(gnx, index, &mtop, as_rvec_array(state.x.data()), as_rvec_array(state.v.data()));
state.natoms = gnx;
}
else if (bZeroQ)
{
FILE *gp;
int indent, i, j, **gcount, atot;
- t_state state;
+ t_state state {};
t_inputrec *ir = nullptr;
t_tpxheader tpx;
gmx_mtop_t mtop;
pr_rvecs(stdout, indent, "svir_prev", tpx.bBox ? state.svir_prev : NULL, DIM);
pr_rvecs(stdout, indent, "fvir_prev", tpx.bBox ? state.fvir_prev : NULL, DIM);
/* leave nosehoover_xi in for now to match the tpr version */
- pr_doubles(stdout, indent, "nosehoover_xi", state.nosehoover_xi, state.ngtc);
+ pr_doubles(stdout, indent, "nosehoover_xi", state.nosehoover_xi.data(), state.ngtc);
/*pr_doubles(stdout,indent,"nosehoover_vxi",state.nosehoover_vxi,state.ngtc);*/
/*pr_doubles(stdout,indent,"therm_integral",state.therm_integral,state.ngtc);*/
- pr_rvecs(stdout, indent, "x", tpx.bX ? state.x : NULL, state.natoms);
- pr_rvecs(stdout, indent, "v", tpx.bV ? state.v : NULL, state.natoms);
+ pr_rvecs(stdout, indent, "x", tpx.bX ? as_rvec_array(state.x.data()) : NULL, state.natoms);
+ pr_rvecs(stdout, indent, "v", tpx.bV ? as_rvec_array(state.v.data()) : NULL, state.natoms);
}
groups = &mtop.groups;
}
sfree(gcount);
}
- done_state(&state);
}
static void list_top(const char *fn)
#include <stdlib.h>
#include <algorithm>
+#include <memory>
#include "thread_mpi/threads.h"
int nchkpt = 1;
gmx_localtop_t *top;
t_mdebin *mdebin = NULL;
- t_state *state = NULL;
gmx_enerdata_t *enerd;
- rvec *f = NULL;
+ PaddedRVecVector f {};
gmx_global_stat_t gstat;
gmx_update_t *upd = NULL;
t_graph *graph = NULL;
{
/* Initialize ion swapping code */
init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
- top_global, state_global->x, state_global->box, &state_global->swapstate, cr, oenv, Flags);
+ top_global, as_rvec_array(state_global->x.data()), state_global->box, &state_global->swapstate, cr, oenv, Flags);
}
/* Initial values */
- init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
+ init_md(fplog, cr, ir, oenv, &t, &t0, &state_global->lambda,
&(state_global->fep_state), lam0,
nrnb, top_global, &upd,
nfile, fnm, &outf, &mdebin,
snew(enerd, 1);
init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
enerd);
- if (DOMAINDECOMP(cr))
+ if (!DOMAINDECOMP(cr))
{
- f = NULL;
- }
- else
- {
- snew(f, top_global->natoms);
+ f.resize(top_global->natoms + 1);
}
/* Kinetic energy data */
}
}
+ std::unique_ptr<t_state> stateInstance;
+ t_state * state;
+
if (DOMAINDECOMP(cr))
{
top = dd_init_local_top(top_global);
- snew(state, 1);
+ stateInstance = std::unique_ptr<t_state>(new t_state {});
+ state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
}
else
mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
&graph, mdatoms, vsite, shellfc);
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
/* Set up interactive MD (IMD) */
- init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+ init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
nfile, fnm, oenv, imdport, Flags);
if (DOMAINDECOMP(cr))
vsite, constr,
nrnb, NULL, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
update_mdatoms(mdatoms, state->lambda[efptMASS]);
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, state->x, ir->delta_t, NULL,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, NULL,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
/* Following is necessary because the graph may get out of sync
* with the coordinates if we only have every N'th coordinate set
*/
- mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
- shift_self(graph, state->box, state->x);
+ mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
}
}
nrnb, wcycle,
do_verbose && !bPMETunePrinting);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
}
relax_shell_flexcon(fplog, cr, bVerbose, step,
ir, bNS, force_flags, top,
constr, enerd, fcd,
- state, f, force_vir, mdatoms,
+ state, &f, force_vir, mdatoms,
nrnb, wcycle, graph, groups,
shellfc, fr, bBornRadii, t, mu_tot,
vsite, mdoutf_get_fp_field(outf));
* Check comments in sim_util.c
*/
do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
- state->lambda, graph,
+ state->box, &state->x, &state->hist,
+ &f, force_vir, mdatoms, enerd, fcd,
+ &state->lambda, graph,
fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
(bNS ? GMX_FORCE_NS : 0) | force_flags);
}
* so that the input is actually the initial step.
*/
snew(vbuf, state->natoms);
- copy_rvecn(state->v, vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
}
else
{
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY1,
cr, constr);
{
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ state, fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
{
/* Need to unshift here if a do_force has been
called in the previous step */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
/* if it's the initial step, we performed this first step just to get the constraint virial */
if (ir->eI == eiVV && bInitStep)
{
- copy_rvecn(vbuf, state->v, 0, state->natoms);
+ copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
sfree(vbuf);
}
wallcycle_stop(wcycle, ewcUPDATE);
/* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
if (ir->efep != efepNO && !bRerunMD)
{
- sum_dhdl(enerd, state->lambda, ir->fepvals);
+ sum_dhdl(enerd, &state->lambda, ir->fepvals);
}
}
statistics, but if performing simulated tempering, we
do update the velocities and the tau_t. */
- lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v, mdatoms);
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
copy_df_history(state_global->dfhist, state->dfhist);
}
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
ir, state, state_global, top_global, fr,
- outf, mdebin, ekind, f,
+ outf, mdebin, ekind, &f,
&nchkpt,
bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
- bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
if (startingFromCheckpoint && bTrotter)
if (constr && bIfRandomize)
{
update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY2,
cr, constr);
}
cbuf_nalloc = state->natoms;
srenew(cbuf, cbuf_nalloc);
}
- copy_rvecn(state->x, cbuf, 0, state->natoms);
+ copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC, graph, f,
+ fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again again */
- copy_rvecn(cbuf, state->x, 0, state->natoms);
+ copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
/* are the small terms in the shake_vir here due
* to numerical errors, or are they important
* physically? I'm thinking they are just errors, but not completely sure.
* For now, will call without actually constraining, constr=NULL*/
update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
cr, nrnb, wcycle, upd, NULL,
FALSE, bCalcVir);
else if (graph)
{
/* Need to unshift here */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
if (vsite != NULL)
wallcycle_start(wcycle, ewcVSITECONSTR);
if (graph != NULL)
{
- shift_self(graph, state->box, state->x);
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph != NULL)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
{
/* Sum up the foreign energy and dhdl terms for md and sd.
Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, ir->fepvals);
+ sum_dhdl(enerd, &state->lambda, ir->fepvals);
}
- update_box(fplog, step, ir, mdatoms, state, f,
+ update_box(fplog, step, ir, mdatoms, state,
pcoupl_mu, nrnb, upd);
/* ################# END UPDATE STEP 2 ################# */
do_per_step(step, ir->swap->nstswap))
{
bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
- bRerunMD ? rerun_fr.x : state->x,
+ bRerunMD ? rerun_fr.x : as_rvec_array(state->x.data()),
bRerunMD ? rerun_fr.box : state->box,
top_global, MASTER(cr) && bVerbose, bRerunMD);
vsite, constr,
nrnb, wcycle, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
bFirstStep = FALSE;
if ( (membed != NULL) && (!bLastStep) )
{
- rescale_membed(step_rel, membed, state_global->x);
+ rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
}
if (bRerunMD)
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
mtop->mols.index = new_mols;
mtop->natoms -= n;
state->natoms -= n;
- state->nalloc = state->natoms;
- snew(x_tmp, state->nalloc);
- snew(v_tmp, state->nalloc);
+ snew(x_tmp, state->natoms);
+ snew(v_tmp, state->natoms);
for (i = 0; i < egcNR; i++)
{
}
}
}
- sfree(state->x);
- state->x = x_tmp;
- sfree(state->v);
- state->v = v_tmp;
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(x_tmp[i], state->x[i]);
+ }
+ sfree(x_tmp);
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(v_tmp[i], state->v[i]);
+ }
+ sfree(v_tmp);
for (i = 0; i < egcNR; i++)
{
/* Check that moleculetypes in insertion group are not part of the rest of the system */
check_types(ins_at, rest_at, mtop);
- init_mem_at(mem_p, mtop, state->x, state->box, pos_ins);
+ init_mem_at(mem_p, mtop, as_rvec_array(state->x.data()), state->box, pos_ins);
- prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p);
+ prot_area = est_prot_area(pos_ins, as_rvec_array(state->x.data()), ins_at, mem_p);
if ( (prot_area > prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX]) < box_vs_prot) )
{
warn++;
/* resize the protein by xy and by z if necessary*/
snew(r_ins, ins_at->nr);
- init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY);
+ init_resize(ins_at, r_ins, pos_ins, mem_p, as_rvec_array(state->x.data()), bALLOW_ASYMMETRY);
membed->fac[0] = membed->fac[1] = xy_fac;
membed->fac[2] = z_fac;
membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
membed->z_step = (z_max-z_fac)/(double)(it_z-1);
- resize(r_ins, state->x, pos_ins, membed->fac);
+ resize(r_ins, as_rvec_array(state->x.data()), pos_ins, membed->fac);
/* remove overlapping lipids and water from the membrane box*/
/*mark molecules to be removed*/
set_pbc(pbc, inputrec->ePBC, state->box);
snew(rm_p, 1);
- lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, mem_p, pos_ins,
+ lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, as_rvec_array(state->x.data()), mem_p, pos_ins,
probe_rad, low_up_rm, bALLOW_ASYMMETRY);
lip_rm -= low_up_rm;
exchange_rvecs(ms, b, state->svir_prev, DIM);
exchange_rvecs(ms, b, state->fvir_prev, DIM);
exchange_rvecs(ms, b, state->pres_prev, DIM);
- exchange_doubles(ms, b, state->nosehoover_xi, ngtc);
- exchange_doubles(ms, b, state->nosehoover_vxi, ngtc);
- exchange_doubles(ms, b, state->nhpres_xi, nnhpres);
- exchange_doubles(ms, b, state->nhpres_vxi, nnhpres);
- exchange_doubles(ms, b, state->therm_integral, state->ngtc);
- exchange_rvecs(ms, b, state->x, state->natoms);
- exchange_rvecs(ms, b, state->v, state->natoms);
+ exchange_doubles(ms, b, state->nosehoover_xi.data(), ngtc);
+ exchange_doubles(ms, b, state->nosehoover_vxi.data(), ngtc);
+ exchange_doubles(ms, b, state->nhpres_xi.data(), nnhpres);
+ exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres);
+ exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc);
+ exchange_rvecs(ms, b, as_rvec_array(state->x.data()), state->natoms);
+ exchange_rvecs(ms, b, as_rvec_array(state->v.data()), state->natoms);
}
static void copy_state_serial(const t_state *src, t_state *dest)
{
int i;
- if (state->v)
+ if (as_rvec_array(state->v.data()))
{
for (i = 0; i < state->natoms; i++)
{
real rlistWithReferenceNstlist, rlist_inc, rlist_ok, rlist_max;
real rlist_new, rlist_prev;
size_t nstlist_ind = 0;
- t_state state_tmp;
gmx_bool bBox, bDD, bCont;
const char *nstl_gpu = "\nFor optimal performance with a GPU nstlist (now %d) should be larger.\nThe optimum depends on your CPU and GPU resources.\nYou might want to try several nstlist values.\n";
const char *nve_err = "Can not increase nstlist because an NVE ensemble is used";
{
gmx_incons("Changing nstlist with domain decomposition and unbounded dimensions is not implemented yet");
}
+ t_state state_tmp {};
copy_mat(box, state_tmp.box);
bDD = change_dd_cutoff(cr, &state_tmp, ir, rlist_new);
}
{
gmx_bool bForceUseGPU, bTryUseGPU, bRerunMD;
t_inputrec *inputrec;
- t_state *state = NULL;
matrix box;
gmx_ddbox_t ddbox = {0};
int npme_major, npme_minor;
please_cite(fplog, "Berendsen95a");
}
- snew(state, 1);
+ std::unique_ptr<t_state> stateInstance = std::unique_ptr<t_state>(new t_state {});
+ t_state * state = stateInstance.get();
+
if (SIMMASTER(cr))
{
/* Read (nearly) all data required for the simulation */
/* This needs to be called before read_checkpoint to extend the state */
init_disres(fplog, mtop, inputrec, cr, fcd, state, repl_ex_nst > 0);
- init_orires(fplog, mtop, state->x, inputrec, cr, &(fcd->orires),
+ init_orires(fplog, mtop, as_rvec_array(state->x.data()), inputrec, cr, &(fcd->orires),
state);
if (inputrecDeform(inputrec))
dddlb_opt, dlb_scale,
ddcsx, ddcsy, ddcsz,
mtop, inputrec,
- box, state->x,
+ box, as_rvec_array(state->x.data()),
&ddbox, &npme_major, &npme_minor);
}
else
/* Make molecules whole at start of run */
if (fr->ePBC != epbcNONE)
{
- do_pbc_first_mtop(fplog, inputrec->ePBC, box, mtop, state->x);
+ do_pbc_first_mtop(fplog, inputrec->ePBC, box, mtop, as_rvec_array(state->x.data()));
}
if (vsite)
{
* for the initial distribution in the domain decomposition
* and for the initial shell prediction.
*/
- construct_vsites_mtop(vsite, mtop, state->x);
+ construct_vsites_mtop(vsite, mtop, as_rvec_array(state->x.data()));
}
}
/* This is a PME only node */
/* We don't need the state */
- done_state(state);
+ stateInstance.reset();
+ state = NULL;
ewaldcoeff_q = calc_ewaldcoeff_q(inputrec->rcoulomb, inputrec->ewald_rtol);
ewaldcoeff_lj = calc_ewaldcoeff_lj(inputrec->rvdw, inputrec->ewald_rtol_lj);
if (inputrec->bRot)
{
/* Initialize enforced rotation code */
- init_rot(fplog, inputrec, nfile, fnm, cr, state->x, state->box, mtop, oenv,
+ init_rot(fplog, inputrec, nfile, fnm, cr, as_rvec_array(state->x.data()), state->box, mtop, oenv,
bVerbose, Flags);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff