Merge branch 'release-5-0'
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / readir.c
index 4b53e5b1e1d47e1e2add3b095c1aacbc4aa394b5..2d86536b178e5ddbd4077a3c550f8c50bca99fb9 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "readir.h"
 
 #include <ctype.h>
-#include <stdlib.h>
 #include <limits.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
-#include "typedefs.h"
-#include "physics.h"
-#include "names.h"
-#include "gmx_fatal.h"
-#include "macros.h"
-#include "index.h"
-#include "symtab.h"
+#include <stdlib.h>
+
+#include "gromacs/gmxpreprocess/toputil.h"
+#include "gromacs/legacyheaders/chargegroup.h"
+#include "gromacs/legacyheaders/inputrec.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/readinp.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/warninp.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/calc_verletbuf.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/symtab.h"
 #include "gromacs/utility/cstringutil.h"
-#include "readinp.h"
-#include "warninp.h"
-#include "readir.h"
-#include "toputil.h"
-#include "index.h"
-#include "network.h"
-#include "vec.h"
-#include "pbc.h"
-#include "mtop_util.h"
-#include "chargegroup.h"
-#include "inputrec.h"
-#include "calc_verletbuf.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
 #define MAXPTR 254
 #define NOGID  255
@@ -160,7 +159,7 @@ static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lamb
         }
         else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
         {
-            simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
+            simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(gmx_expm1(temperature_lambdas[i])/gmx_expm1(1.0));
         }
         else
         {
@@ -270,6 +269,8 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
     {
         warning_error(wi, "rlist should be >= 0");
     }
+    sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
+    CHECK(ir->nstlist < 0);
 
     process_interaction_modifier(ir, &ir->coulomb_modifier);
     process_interaction_modifier(ir, &ir->vdw_modifier);
@@ -306,9 +307,9 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         {
             warning_error(wi, "rlistlong can not be shorter than rlist");
         }
-        if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
+        if (IR_TWINRANGE(*ir) && ir->nstlist == 0)
         {
-            warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
+            warning_error(wi, "Can not have nstlist == 0 with twin-range interactions");
         }
     }
 
@@ -880,10 +881,6 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
               (ir->ePBC     != epbcNONE) ||
               (ir->rcoulomb != 0.0)      || (ir->rvdw != 0.0));
 
-        if (ir->nstlist < 0)
-        {
-            warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
-        }
         if (ir->nstlist > 0)
         {
             warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
@@ -1285,8 +1282,7 @@ nd %s",
     if (ir->cutoff_scheme == ecutsGROUP)
     {
         if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
-             (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
-            ir->nstlist != 1)
+             (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
         {
             warning_note(wi, "With exact cut-offs, rlist should be "
                          "larger than rcoulomb and rvdw, so that there "
@@ -1313,14 +1309,6 @@ nd %s",
         warning_note(wi, "You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
     }
 
-    if (ir->nstlist == -1)
-    {
-        sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
-        CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
-    }
-    sprintf(err_buf, "nstlist can not be smaller than -1");
-    CHECK(ir->nstlist < -1);
-
     if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
         && ir->rvdw != 0)
     {
@@ -1448,7 +1436,7 @@ int str_nelem(const char *str, int maxptr, char *ptr[])
     int   np = 0;
     char *copy0, *copy;
 
-    copy0 = strdup(str);
+    copy0 = gmx_strdup(str);
     copy  = copy0;
     ltrim(copy);
     while (*copy != '\0')
@@ -1703,7 +1691,7 @@ static void do_wall_params(t_inputrec *ir,
         }
         for (i = 0; i < ir->nwall; i++)
         {
-            opts->wall_atomtype[i] = strdup(names[i]);
+            opts->wall_atomtype[i] = gmx_strdup(names[i]);
         }
 
         if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
@@ -2094,12 +2082,11 @@ void get_ir(const char *mdparin, const char *mdparout,
 
     /* COM pulling */
     CCTYPE("COM PULLING");
-    CTYPE("Pull type: no, umbrella, constraint or constant-force");
-    EETYPE("pull",          ir->ePull, epull_names);
-    if (ir->ePull != epullNO)
+    EETYPE("pull",          ir->bPull, yesno_names);
+    if (ir->bPull)
     {
         snew(ir->pull, 1);
-        is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
+        is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
     }
 
     /* Enforced rotation */
@@ -2369,7 +2356,7 @@ void get_ir(const char *mdparin, const char *mdparout,
     {
         if (ir->efep != efepNO)
         {
-            opts->couple_moltype = strdup(is->couple_moltype);
+            opts->couple_moltype = gmx_strdup(is->couple_moltype);
             if (opts->couple_lam0 == opts->couple_lam1)
             {
                 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
@@ -2547,6 +2534,9 @@ static int search_QMstring(const char *s, int ng, const char *gn[])
 } /* search_QMstring */
 
 /* We would like gn to be const as well, but C doesn't allow this */
+/* TODO this is utility functionality (search for the index of a
+   string in a collection), so should be refactored and located more
+   centrally. */
 int search_string(const char *s, int ng, char *gn[])
 {
     int i;
@@ -2852,7 +2842,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
         }
     }
 
-    if (ir->ePull == epullCONSTRAINT)
+    if (ir->bPull)
     {
         /* Correct nrdf for the COM constraints.
          * We correct using the TC and VCM group of the first atom
@@ -2864,6 +2854,11 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
 
         for (i = 0; i < pull->ncoord; i++)
         {
+            if (pull->coord[i].eType != epullCONSTRAINT)
+            {
+                continue;
+            }
+
             imin = 1;
 
             for (j = 0; j < 2; j++)
@@ -2977,7 +2972,7 @@ static void decode_cos(char *s, t_cosines *cosine)
     double  a, phi;
     int     i;
 
-    t = strdup(s);
+    t = gmx_strdup(s);
     trim(t);
 
     cosine->n   = 0;
@@ -3158,7 +3153,7 @@ static void make_swap_groups(
 
 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
 {
-    int      ig = -1, i;
+    int      ig, i;
 
 
     ig            = search_string(IMDgname, grps->nr, gnames);
@@ -3489,7 +3484,7 @@ void do_index(const char* mdparin, const char *ndx,
         }
     }
 
-    if (ir->ePull != epullNO)
+    if (ir->bPull)
     {
         make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
 
@@ -3865,7 +3860,15 @@ static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
                         case efbposresSPHERE:
                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
                             break;
+                        case efbposresCYLINDERX:
+                            AbsRef[YY] = AbsRef[ZZ] = 1;
+                            break;
+                        case efbposresCYLINDERY:
+                            AbsRef[XX] = AbsRef[ZZ] = 1;
+                            break;
                         case efbposresCYLINDER:
+                        /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
+                        case efbposresCYLINDERZ:
                             AbsRef[XX] = AbsRef[YY] = 1;
                             break;
                         case efbposresX: /* d=XX */
@@ -4239,45 +4242,42 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
     }
 
-    if (ir->ePull != epullNO)
+    if (ir->bPull)
     {
-        gmx_bool bPullAbsoluteRef;
+        gmx_bool bWarned;
 
-        bPullAbsoluteRef = FALSE;
-        for (i = 0; i < ir->pull->ncoord; i++)
+        bWarned = FALSE;
+        for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
         {
-            bPullAbsoluteRef = bPullAbsoluteRef ||
-                ir->pull->coord[i].group[0] == 0 ||
-                ir->pull->coord[i].group[1] == 0;
-        }
-        if (bPullAbsoluteRef)
-        {
-            absolute_reference(ir, sys, FALSE, AbsRef);
-            for (m = 0; m < DIM; m++)
+            if (ir->pull->coord[i].group[0] == 0 ||
+                ir->pull->coord[i].group[1] == 0)
             {
-                if (ir->pull->dim[m] && !AbsRef[m])
+                absolute_reference(ir, sys, FALSE, AbsRef);
+                for (m = 0; m < DIM; m++)
                 {
-                    warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
-                    break;
+                    if (ir->pull->coord[i].dim[m] && !AbsRef[m])
+                    {
+                        warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
+                        bWarned = TRUE;
+                        break;
+                    }
                 }
             }
         }
 
-        if (ir->pull->eGeom == epullgDIRPBC)
+        for (i = 0; i < 3; i++)
         {
-            for (i = 0; i < 3; i++)
+            for (m = 0; m <= i; m++)
             {
-                for (m = 0; m <= i; m++)
+                if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
+                    ir->deform[i][m] != 0)
                 {
-                    if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
-                        ir->deform[i][m] != 0)
+                    for (c = 0; c < ir->pull->ncoord; c++)
                     {
-                        for (c = 0; c < ir->pull->ncoord; c++)
+                        if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
+                            ir->pull->coord[c].vec[m] != 0)
                         {
-                            if (ir->pull->coord[c].vec[m] != 0)
-                            {
-                                gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
-                            }
+                            gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
                         }
                     }
                 }
@@ -4288,7 +4288,10 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
     check_disre(sys);
 }
 
-void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
+void double_check(t_inputrec *ir, matrix box,
+                  gmx_bool bHasNormalConstraints,
+                  gmx_bool bHasAnyConstraints,
+                  warninp_t wi)
 {
     real        min_size;
     gmx_bool    bTWIN;
@@ -4301,7 +4304,7 @@ void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
         warning_error(wi, ptr);
     }
 
-    if (bConstr && ir->eConstrAlg == econtSHAKE)
+    if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
     {
         if (ir->shake_tol <= 0.0)
         {
@@ -4324,7 +4327,7 @@ void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
         }
     }
 
-    if ( (ir->eConstrAlg == econtLINCS) && bConstr)
+    if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
     {
         /* If we have Lincs constraints: */
         if (ir->eI == eiMD && ir->etc == etcNO &&
@@ -4345,9 +4348,9 @@ void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
         }
     }
 
-    if (bConstr && ir->epc == epcMTTK)
+    if (bHasAnyConstraints && ir->epc == epcMTTK)
     {
-        warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
+        warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
     }
 
     if (ir->LincsWarnAngle > 90.0)