patch-v2019.4

patch -p 1 << EOF
diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 6b82ca99b2..6855a15494 100644
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -1026,7 +1026,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
     {
         if (ir->epc > epcNO)
         {
-            if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
+            if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK) && (ir->epc != epcCRESCALE))
             {
                 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
             }
diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index 5927ace633..31ac1264a9 100644
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -736,6 +736,242 @@ void berendsen_pscale(const t_inputrec *ir, const matrix mu,
     inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
 }
 
+void crescale_pcoupl(FILE *fplog, int64_t step,
+                      const t_inputrec *ir, real dt,
+                      const tensor pres, const matrix box,
+                      const matrix force_vir, const matrix constraint_vir,
+                      matrix mu, double *baros_integral)
+{
+    int     d, n;
+    real    scalar_pressure, xy_pressure, p_corr_z;
+    char    buf[STRLEN];
+
+    /*
+     *  Calculate the scaling matrix mu
+     */
+    scalar_pressure = 0;
+    xy_pressure     = 0;
+    for (d = 0; d < DIM; d++)
+    {
+        scalar_pressure += pres[d][d]/DIM;
+        if (d != ZZ)
+        {
+            xy_pressure += pres[d][d]/(DIM-1);
+        }
+    }
+    /* Pressure is now in bar, everywhere. */
+#define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
+
+    /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
+     * necessary for triclinic scaling
+     */
+    clear_mat(mu);
+
+    /* C-RESCALE */
+    // seed hardcoded to zero now
+    gmx::ThreeFry2x64<64>          rng(0, gmx::RandomDomain::Barostat);
+    gmx::NormalDistribution<real>  normalDist;
+    rng.restart(step, 0);
+    real vol=1.0;
+    for (d = 0; d < DIM; d++) vol*=box[d][d];
+    real gauss;
+    real gauss2;
+    real kt;
+    real depsilon, dsigma;
+    kt=ir->opts.ref_t[0]*BOLTZ;
+    if(kt<0.0) kt=0.0;
+    /* END C-RESCALE */
+
+    switch (ir->epct)
+    {
+        case epctISOTROPIC:
+            /* C-RESCALE */
+            gauss=normalDist(rng);
+            vol=1.0; for (d = 0; d < DIM; d++) vol*=box[d][d];
+            /* END C-RESCALE */
+            for (d = 0; d < DIM; d++)
+            {
+                /* Berendsen: */
+                /* mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM; */
+                /* C-RESCALE */
+                mu[d][d] = exp(- factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM +
+                            sqrt(2.0*kt*factor(d, d)*PRESFAC/vol)*gauss/DIM);
+            }
+            break;
+        case epctSEMIISOTROPIC:
+            /* C-RESCALE */
+            gauss=normalDist(rng);
+            gauss2=normalDist(rng);
+            vol=1.0; for (d = 0; d < DIM; d++) vol*=box[d][d];
+            for (d = 0; d < ZZ; d++)
+            {
+                /* Berendsen: */
+                /* mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM; */
+                /* C-RESCALE */
+                mu[d][d] = exp(- factor(d, d)*(ir->ref_p[d][d] - xy_pressure) /DIM +
+                            sqrt((DIM-1)*2.0*kt*factor(d, d)*PRESFAC/vol/DIM)/(DIM-1)*gauss);
+            }
+            /* Notice: we here allow ir->ref_p[ZZ][ZZ] != ir->ref_p[d][d] although this is not well defined */
+            /* Berendsen: 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM; */
+            mu[ZZ][ZZ] = exp(- factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) /DIM +
+                            sqrt(2.0*kt*factor(d, d)*PRESFAC/vol/DIM)*gauss2);
+            break;
+#if 0
+        case epctANISOTROPIC:
+            for (d = 0; d < DIM; d++)
+            {
+                for (n = 0; n < DIM; n++)
+                {
+                    mu[d][n] = (d == n ? 1.0 : 0.0)
+                        -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
+                }
+            }
+            break;
+#endif
+        case epctSURFACETENSION:
+            /* C-RESCALE */
+            gauss=normalDist(rng);
+            gauss2=normalDist(rng);
+            vol=1.0; for (d = 0; d < DIM; d++) vol*=box[d][d];
+            for (d = 0; d < ZZ; d++)
+            {
+                /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
+                mu[d][d] = exp(- factor(d, d)*(ir->ref_p[ZZ][ZZ]-ir->ref_p[d][d]/box[ZZ][ZZ] - xy_pressure) /DIM +
+                            sqrt(4.0/3.0*kt*factor(d, d)*PRESFAC/vol)/(DIM-1)*gauss);
+            }
+            mu[ZZ][ZZ] = exp(- factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) /DIM +
+                            sqrt(2.0/3.0*kt*factor(d, d)*PRESFAC/vol)*gauss2);
+            break;
+        default:
+            gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
+                      EPCOUPLTYPETYPE(ir->epct));
+    }
+    /* To fullfill the orientation restrictions on triclinic boxes
+     * we will set mu_yx, mu_zx and mu_zy to 0 and correct
+     * the other elements of mu to first order.
+     */
+    mu[YY][XX] += mu[XX][YY];
+    mu[ZZ][XX] += mu[XX][ZZ];
+    mu[ZZ][YY] += mu[YY][ZZ];
+    mu[XX][YY]  = 0;
+    mu[XX][ZZ]  = 0;
+    mu[YY][ZZ]  = 0;
+
+    /* Keep track of the work the barostat applies on the system.
+     * Without constraints force_vir tells us how Epot changes when scaling.
+     * With constraints constraint_vir gives us the constraint contribution
+     * to both Epot and Ekin. Although we are not scaling velocities, scaling
+     * the coordinates leads to scaling of distances involved in constraints.
+     * This in turn changes the angular momentum (even if the constrained
+     * distances are corrected at the next step). The kinetic component
+     * of the constraint virial captures the angular momentum change.
+     */
+    for (int d = 0; d < DIM; d++)
+    {
+        for (int n = 0; n <= d; n++)
+        {
+            *baros_integral -= 2*(mu[d][n] - (n == d ? 1 : 0))*(force_vir[d][n] + constraint_vir[d][n]);
+        }
+    }
+
+    if (debug)
+    {
+        pr_rvecs(debug, 0, "PC: pres ", pres, 3);
+        pr_rvecs(debug, 0, "PC: mu   ", mu, 3);
+    }
+
+    if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
+        mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
+        mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
+    {
+        char buf2[22];
+        sprintf(buf, "\nStep %s  Warning: pressure scaling more than 1%%, "
+                "mu: %g %g %g\n",
+                gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
+        if (fplog)
+        {
+            fprintf(fplog, "%s", buf);
+        }
+        fprintf(stderr, "%s", buf);
+    }
+}
+
+void crescale_pscale(const t_inputrec *ir, const matrix mu,
+                      matrix box, matrix box_rel,
+                      int start, int nr_atoms,
+                      rvec x[], rvec v[], const unsigned short cFREEZE[],
+                      t_nrnb *nrnb)
+{
+    ivec   *nFreeze = ir->opts.nFreeze;
+    int     n, d;
+    int     nthreads gmx_unused;
+    matrix  inv_mu;
+
+#ifndef __clang_analyzer__
+    nthreads = gmx_omp_nthreads_get(emntUpdate);
+#endif
+
+    /* find the inverse matrix */
+    inv_mu[XX][XX]=1.0/mu[XX][XX];
+    inv_mu[YY][YY]=1.0/mu[YY][YY];
+    inv_mu[ZZ][ZZ]=1.0/mu[ZZ][ZZ];
+
+    inv_mu[YY][XX]=-mu[YY][XX]/(mu[XX][XX]*mu[YY][YY]);
+    inv_mu[YY][ZZ]=-mu[YY][ZZ]/(mu[YY][YY]*mu[ZZ][ZZ]);
+    inv_mu[ZZ][XX]=(mu[YY][XX]*mu[YY][ZZ] - mu[YY][YY]*mu[ZZ][XX]) / (mu[XX][XX]*mu[YY][YY]*mu[ZZ][ZZ]);
+
+    inv_mu[XX][YY]=0.0;
+    inv_mu[ZZ][YY]=0.0;
+    inv_mu[ZZ][XX]=0.0;
+
+    /* Scale the positions and the velocities */
+#pragma omp parallel for num_threads(nthreads) schedule(static)
+    for (n = start; n < start+nr_atoms; n++)
+    {
+        // Trivial OpenMP region that does not throw
+        int g;
+
+        if (cFREEZE == nullptr)
+        {
+            g = 0;
+        }
+        else
+        {
+            g = cFREEZE[n];
+        }
+
+        if (!nFreeze[g][XX])
+        {
+            x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
+            v[n][XX] = inv_mu[XX][XX]*v[n][XX]+inv_mu[YY][XX]*v[n][YY]+inv_mu[ZZ][XX]*v[n][ZZ];
+        }
+        if (!nFreeze[g][YY])
+        {
+            x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
+            v[n][YY] = inv_mu[YY][YY]*v[n][YY]+inv_mu[ZZ][YY]*v[n][ZZ];
+        }
+        if (!nFreeze[g][ZZ])
+        {
+            x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
+            v[n][ZZ] = inv_mu[ZZ][ZZ]*v[n][ZZ];
+        }
+    }
+    /* compute final boxlengths */
+    for (d = 0; d < DIM; d++)
+    {
+        box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
+        box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
+        box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
+    }
+
+    preserve_box_shape(ir, box_rel, box);
+
+    /* (un)shifting should NOT be done after this,
+     * since the box vectors might have changed
+     */
+    inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
+}
+
 void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt,
                       std::vector<double> &therm_integral)
 {
@@ -1405,6 +1641,9 @@ real NPT_energy(const t_inputrec *ir, const t_state *state, const t_extmass *Mas
             case epcBERENDSEN:
                 energyNPT += state->baros_integral;
                 break;
+            case epcCRESCALE:
+                energyNPT += state->baros_integral;
+                break;
             default:
                 GMX_RELEASE_ASSERT(false, "Conserved energy quantity for pressure coupling is not handled. A case should be added with either the conserved quantity added or nothing added and an exclusion added to integratorHasConservedEnergyQuantity().");
         }
diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp
index af0e485c75..ad3295efe6 100644
--- a/src/gromacs/mdlib/update.cpp
+++ b/src/gromacs/mdlib/update.cpp
@@ -1721,6 +1721,21 @@ void update_pcouple_after_coordinates(FILE             *fplog,
                                  md->cFREEZE, nrnb);
             }
             break;
+        case (epcCRESCALE):
+            if (isPressureCouplingStep(step, inputrec))
+            {
+                real   dtpc = inputrec->nstpcouple*dt;
+                matrix mu;
+                crescale_pcoupl(fplog, step, inputrec, dtpc,
+                                 pressure, state->box,
+                                 forceVirial, constraintVirial,
+                                 mu, &state->baros_integral);
+                crescale_pscale(inputrec, mu, state->box, state->box_rel,
+                                 start, homenr, state->x.rvec_array(),
+                                 state->v.rvec_array(),
+                                 md->cFREEZE, nrnb);
+            }
+            break;
         case (epcPARRINELLORAHMAN):
             if (isPressureCouplingStep(step, inputrec))
             {
diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h
index 81b1ace278..fa905622ce 100644
--- a/src/gromacs/mdlib/update.h
+++ b/src/gromacs/mdlib/update.h
@@ -278,6 +278,18 @@ void berendsen_pscale(const t_inputrec *ir, const matrix mu,
                       rvec x[], const unsigned short cFREEZE[],
                       t_nrnb *nrnb);
 
+void crescale_pcoupl(FILE *fplog, int64_t step,
+                     const t_inputrec *ir, real dt,
+                     const tensor pres, const matrix box,
+                     const matrix force_vir, const matrix constraint_vir,
+                     matrix mu, double *baros_integral);
+
+void crescale_pscale(const t_inputrec *ir, const matrix mu,
+                     matrix box, matrix box_rel,
+                     int start, int nr_atoms,
+                     rvec x[], rvec v[], const unsigned short cFREEZE[],
+                     t_nrnb *nrnb);
+
 // TODO: This doesn't seem to be used or implemented anywhere
 void correct_ekin(FILE *log, int start, int end, rvec v[],
                   rvec vcm, real mass[], real tmass, tensor ekin);
diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp
index 1afe02a840..596d8b5bf3 100644
--- a/src/gromacs/mdtypes/inputrec.cpp
+++ b/src/gromacs/mdtypes/inputrec.cpp
@@ -183,7 +183,8 @@ int pcouple_min_integration_steps(int epc)
         case epcNO:
             n = 0;
             break;
-        case etcBERENDSEN:
+        case epcBERENDSEN:
+        case epcCRESCALE:
         case epcISOTROPIC:
             n = nstmin_berendsen_pcouple;
             break;
diff --git a/src/gromacs/mdtypes/md_enums.cpp b/src/gromacs/mdtypes/md_enums.cpp
index 3b03645e3f..2d8308ea19 100644
--- a/src/gromacs/mdtypes/md_enums.cpp
+++ b/src/gromacs/mdtypes/md_enums.cpp
@@ -108,7 +108,7 @@ const char *etcoupl_names[etcNR+1] = {
 }; /* yes is alias for berendsen */
 
 const char *epcoupl_names[epcNR+1] = {
-    "No", "Berendsen", "Parrinello-Rahman", "Isotropic", "MTTK", nullptr
+    "No", "Berendsen", "Parrinello-Rahman", "Isotropic", "MTTK", "C-rescale", nullptr
 }; /* isotropic is alias for berendsen */
 
 const char *epcoupltype_names[epctNR+1] = {
diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h
index f5cf68e68e..df0bb57d12 100644
--- a/src/gromacs/mdtypes/md_enums.h
+++ b/src/gromacs/mdtypes/md_enums.h
@@ -94,7 +94,7 @@ extern const char *etcoupl_names[etcNR+1];
  * isotropic is an alias for berendsen
  */
 enum {
-    epcNO, epcBERENDSEN, epcPARRINELLORAHMAN, epcISOTROPIC, epcMTTK, epcNR
+    epcNO, epcBERENDSEN, epcPARRINELLORAHMAN, epcISOTROPIC, epcMTTK, epcCRESCALE, epcNR
 };
 //! String corresponding to pressure coupling algorithm
 extern const char *epcoupl_names[epcNR+1];
EOF