repl_ex.cpp

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team.
 * Copyright (c) 2011,2012,2013,2014,2015,2016,2017, 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.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */

#include "gmxpre.h"

#include "repl_ex.h"

#include "config.h"

#include <cmath>

#include <random>

#include "gromacs/domdec/domdec.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/main.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/random/threefry.h"
#include "gromacs/random/uniformintdistribution.h"
#include "gromacs/random/uniformrealdistribution.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"


/* PLUMED */
#include "../../../Plumed.h"
extern int    plumedswitch;
extern plumed plumedmain;
/* END PLUMED */

/* PLUMED HREX */
extern int plumed_hrex;
/* END PLUMED HREX */

#define PROBABILITYCUTOFF 100
/* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */

//! Rank in the multisimulaiton
#define MSRANK(ms, nodeid)  (nodeid)

enum {
    ereTEMP, ereLAMBDA, ereENDSINGLE, ereTL, ereNR
};
const char *erename[ereNR] = { "temperature", "lambda", "end_single_marker", "temperature and lambda"};
/* end_single_marker merely notes the end of single variable replica exchange. All types higher than
   it are multiple replica exchange methods */
/* Eventually, should add 'pressure', 'temperature and pressure', 'lambda_and_pressure', 'temperature_lambda_pressure'?;
   Let's wait until we feel better about the pressure control methods giving exact ensembles.  Right now, we assume constant pressure  */

typedef struct gmx_repl_ex
{
    int       repl;        /* replica ID */
    int       nrepl;       /* total number of replica */
    real      temp;        /* temperature */
    int       type;        /* replica exchange type from ere enum */
    real    **q;           /* quantity, e.g. temperature or lambda; first index is ere, second index is replica ID */
    gmx_bool  bNPT;        /* use constant pressure and temperature */
    real     *pres;        /* replica pressures */
    int      *ind;         /* replica indices */
    int      *allswaps;    /* used for keeping track of all the replica swaps */
    int       nst;         /* replica exchange interval (number of steps) */
    int       nex;         /* number of exchanges per interval */
    int       seed;        /* random seed */
    int       nattempt[2]; /* number of even and odd replica change attempts */
    real     *prob_sum;    /* sum of probabilities */
    int     **nmoves;      /* number of moves between replicas i and j */
    int      *nexchange;   /* i-th element of the array is the number of exchanges between replica i-1 and i */

    /* these are helper arrays for replica exchange; allocated here so they
       don't have to be allocated each time */
    int      *destinations;
    int     **cyclic;
    int     **order;
    int      *tmpswap;
    gmx_bool *incycle;
    gmx_bool *bEx;

    /* helper arrays to hold the quantities that are exchanged */
    real  *prob;
    real  *Epot;
    real  *beta;
    real  *Vol;
    real **de;

} t_gmx_repl_ex;

static gmx_bool repl_quantity(const gmx_multisim_t *ms,
                              struct gmx_repl_ex *re, int ere, real q)
{
    real    *qall;
    gmx_bool bDiff;
    int      s;

    snew(qall, ms->nsim);
    qall[re->repl] = q;
    gmx_sum_sim(ms->nsim, qall, ms);

    /* PLUMED */
    //bDiff = FALSE;
    //for (s = 1; s < ms->nsim; s++)
    //{
    //    if (qall[s] != qall[0])
    //    {
            bDiff = TRUE;
    //    }
    //}
    /* END PLUMED */

    if (bDiff)
    {
        /* Set the replica exchange type and quantities */
        re->type = ere;

        snew(re->q[ere], re->nrepl);
        for (s = 0; s < ms->nsim; s++)
        {
            re->q[ere][s] = qall[s];
        }
    }
    sfree(qall);
    return bDiff;
}

gmx_repl_ex_t
init_replica_exchange(FILE                            *fplog,
                      const gmx_multisim_t            *ms,
                      int                              numAtomsInSystem,
                      const t_inputrec                *ir,
                      const ReplicaExchangeParameters &replExParams)
{
    real                pres;
    int                 i, j, k;
    struct gmx_repl_ex *re;
    gmx_bool            bTemp;
    gmx_bool            bLambda = FALSE;

    fprintf(fplog, "\nInitializing Replica Exchange\n");

    if (ms == nullptr || ms->nsim == 1)
    {
        gmx_fatal(FARGS, "Nothing to exchange with only one replica, maybe you forgot to set the -multi option of mdrun?");
    }
    if (!EI_DYNAMICS(ir->eI))
    {
        gmx_fatal(FARGS, "Replica exchange is only supported by dynamical simulations");
        /* Note that PAR(cr) is defined by cr->nnodes > 1, which is
         * distinct from MULTISIM(cr). A multi-simulation only runs
         * with real MPI parallelism, but this does not imply PAR(cr)
         * is true!
         *
         * Since we are using a dynamical integrator, the only
         * decomposition is DD, so PAR(cr) and DOMAINDECOMP(cr) are
         * synonymous. The only way for cr->nnodes > 1 to be true is
         * if we are using DD. */
    }

    snew(re, 1);

    re->repl     = ms->sim;
    re->nrepl    = ms->nsim;