Merge branch release-2016
[alexxy/gromacs.git] / src / gromacs / mdlib / md_support.cpp
index 2663e635d1bf4c6bc0f7490ee3b7db1302e9069f..39c7706ee5beef3f98e24c52dd9db8f7124b4201 100644 (file)
@@ -44,7 +44,6 @@
 #include <algorithm>
 
 #include "gromacs/domdec/domdec.h"
-#include "gromacs/gmxlib/md_logging.h"
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/math/vec.h"
@@ -69,6 +68,7 @@
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/logger.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/snprintf.h"
 
@@ -146,93 +146,6 @@ int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
     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;
-    switch (ir->etc)
-    {
-        case etcNO:
-            break;
-        case etcBERENDSEN:
-            break;
-        case etcNOSEHOOVER:
-            quantity = NPT_energy(ir, state, MassQ);
-            break;
-        case etcVRESCALE:
-            quantity = vrescale_energy(&(ir->opts), state->therm_integral);
-            break;
-        default:
-            break;
-    }
-    return quantity;
-}
-
 /* TODO Specialize this routine into init-time and loop-time versions?
    e.g. bReadEkin is only true when restoring from checkpoint */
 void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
@@ -291,7 +204,7 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input
     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)
@@ -311,7 +224,7 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input
             {
                 wallcycle_start(wcycle, ewcMoveE);
                 global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
-                            ir, ekind, constr, bStopCM ? vcm : NULL,
+                            ir, ekind, constr, bStopCM ? vcm : nullptr,
                             signalBuffer.size(), signalBuffer.data(),
                             totalNumberOfBondedInteractions,
                             *bSumEkinhOld, flags);
@@ -326,7 +239,7 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input
     {
         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);
     }
 
@@ -335,7 +248,7 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input
     {
         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);
     }
 
@@ -400,16 +313,18 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input
     }
 }
 
-void check_nst_param(FILE *fplog, t_commrec *cr,
-                     const char *desc_nst, int nst,
-                     const char *desc_p, int *p)
+/* check whether an 'nst'-style parameter p is a multiple of nst, and
+   set it to be one if not, with a warning. */
+static void check_nst_param(const gmx::MDLogger &mdlog,
+                            const char *desc_nst, int nst,
+                            const char *desc_p, int *p)
 {
     if (*p > 0 && *p % nst != 0)
     {
         /* Round up to the next multiple of nst */
         *p = ((*p)/nst + 1)*nst;
-        md_print_warn(cr, fplog,
-                      "NOTE: %s changes %s to %d\n", desc_nst, desc_p, *p);
+        GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+                "NOTE: %s changes %s to %d", desc_nst, desc_p, *p);
     }
 }
 
@@ -537,8 +452,7 @@ static int lcd4(int i1, int i2, int i3, int i4)
     return nst;
 }
 
-int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
-                        int nstglobalcomm, t_inputrec *ir)
+int check_nstglobalcomm(const gmx::MDLogger &mdlog, int nstglobalcomm, t_inputrec *ir)
 {
     if (!EI_DYNAMICS(ir->eI))
     {
@@ -586,42 +500,43 @@ int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
         {
             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
-            md_print_warn(cr, fplog, "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n", nstglobalcomm);
+            GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+                    "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d",
+                    nstglobalcomm);
         }
         if (ir->nstcalcenergy > 0)
         {
-            check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
+            check_nst_param(mdlog, "-gcom", nstglobalcomm,
                             "nstcalcenergy", &ir->nstcalcenergy);
         }
         if (ir->etc != etcNO && ir->nsttcouple > 0)
         {
-            check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
+            check_nst_param(mdlog, "-gcom", nstglobalcomm,
                             "nsttcouple", &ir->nsttcouple);
         }
         if (ir->epc != epcNO && ir->nstpcouple > 0)
         {
-            check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
+            check_nst_param(mdlog, "-gcom", nstglobalcomm,
                             "nstpcouple", &ir->nstpcouple);
         }
 
-        check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
+        check_nst_param(mdlog, "-gcom", nstglobalcomm,
                         "nstenergy", &ir->nstenergy);
 
-        check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
+        check_nst_param(mdlog, "-gcom", nstglobalcomm,
                         "nstlog", &ir->nstlog);
     }
 
     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
     {
-        md_print_warn(cr, fplog, "WARNING: Changing nstcomm from %d to %d\n",
-                      ir->nstcomm, nstglobalcomm);
+        GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+                "WARNING: Changing nstcomm from %d to %d",
+                ir->nstcomm, nstglobalcomm);
         ir->nstcomm = nstglobalcomm;
     }
 
-    if (fplog)
-    {
-        fprintf(fplog, "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm);
-    }
+    GMX_LOG(mdlog.info).appendTextFormatted(
+            "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm);
     return nstglobalcomm;
 }
 
@@ -656,33 +571,17 @@ void set_state_entries(t_state *state, const t_inputrec *ir)
         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);
-    }
+    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;
@@ -696,23 +595,17 @@ void set_state_entries(t_state *state, const t_inputrec *ir)
         if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
         {
             state->flags |= (1<<estBOXV);
+            state->flags |= (1<<estPRES_PREV);
         }
-        if (ir->epc != epcNO)
+        if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
         {
-            if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
-            {
-                state->nnhpres = 1;
-                state->flags  |= (1<<estNHPRES_XI);
-                state->flags  |= (1<<estNHPRES_VXI);
-                state->flags  |= (1<<estSVIR_PREV);
-                state->flags  |= (1<<estFVIR_PREV);
-                state->flags  |= (1<<estVETA);
-                state->flags  |= (1<<estVOL0);
-            }
-            else
-            {
-                state->flags |= (1<<estPRES_PREV);
-            }
+            state->nnhpres = 1;
+            state->flags  |= (1<<estNHPRES_XI);
+            state->flags  |= (1<<estNHPRES_VXI);
+            state->flags  |= (1<<estSVIR_PREV);
+            state->flags  |= (1<<estFVIR_PREV);
+            state->flags  |= (1<<estVETA);
+            state->flags  |= (1<<estVOL0);
         }
     }
 
@@ -729,8 +622,18 @@ void set_state_entries(t_state *state, const t_inputrec *ir)
 
     init_gtc_state(state, state->ngtc, state->nnhpres, ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
     init_ekinstate(&state->ekinstate, ir);
-    snew(state->enerhist, 1);
-    init_energyhistory(state->enerhist);
-    init_df_history(&state->dfhist, ir->fepvals->n_lambda);
-    state->swapstate.eSwapCoords = ir->eSwapCoords;
+
+    if (ir->bExpanded)
+    {
+        snew(state->dfhist, 1);
+        init_df_history(state->dfhist, ir->fepvals->n_lambda);
+    }
+    if (ir->eSwapCoords != eswapNO)
+    {
+        if (state->swapstate == nullptr)
+        {
+            snew(state->swapstate, 1);
+        }
+        state->swapstate->eSwapCoords = ir->eSwapCoords;
+    }
 }