implemented free-energy coupling options in grompp and added dG/dlambda term for...

diff --git a/include/grompp.h b/include/grompp.h
index f84ba2ecb070d0694dd33642ca5662aca2e69c73..7ce3d7ab28fa3ce646dc22ff3098f630e29ae3b1 100644 (file)
--- a/include/grompp.h
+++ b/include/grompp.h
@@ -55,9 +55,7 @@ typedef struct {
  * non-bonded parameter combinations, which will be copied to t_params.
  */
 
-#define DEF_PARAM_TYPE NOTSET
 typedef struct {
-  t_functype type;
   atom_id    a[MAXATOMLIST];   /* The atom list (eg. bonds: particle   */
                                /* i = a[0] (AI), j = a[1] (AJ))        */
   real              c[MAXFORCEPARAM];  /* Force parameters (eg. b0 = c[0])     */
@@ -105,6 +103,7 @@ typedef enum {
   d_bonds,
   d_exclusions,
   d_pairs,
+  d_pairs_nb,
   d_angles,
   d_dihedrals,
   d_constraints,
@@ -144,6 +143,7 @@ static char *ds[d_maxdir+1] = {
   "bonds",
   "exclusions",
   "pairs",
+  "pairs_nb",
   "angles",
   "dihedrals",
   "constraints",

diff --git a/include/types/idef.h b/include/types/idef.h
index 8ad9a27be63bb322ff7b229a0d9f1260d70fb014..ffa06c9f37eac3e4846ca00d43ad867ea8882dc4 100644 (file)
--- a/include/types/idef.h
+++ b/include/types/idef.h
@@ -76,8 +76,8 @@ enum {
   F_TABDIHS,
   F_LJ14,
   F_COUL14,
-  F_LJC14_A,
-  F_LJC_PAIRS_A,
+  F_LJC14_Q,
+  F_LJC_PAIRS_NB,
   F_LJ,
   F_BHAM,
   F_LJ_LR,
@@ -120,7 +120,8 @@ enum {
   F_TEMP,
   F_PRES,
   F_DVDL,
-  F_DVDLKIN,
+  F_DKDL,
+  F_DGDL_CON,
   F_NRE                /* This number is for the total number of energies      */
 };
   
@@ -146,6 +147,8 @@ typedef union
   struct {real a,alpha1,alpha2,rfac;                       } thole;
   struct {real c6,c12;                                    } lj;
   struct {real c6A,c12A,c6B,c12B;                         } lj14;
+  struct {real fqq,qi,qj,c6,c12;                          } ljc14;
+  struct {real qi,qj,c6,c12;                              } ljcnb;
   /* Proper dihedrals can not have different multiplicity when
    * doing free energy calculations, because the potential would not
    * be periodic anymore.
@@ -204,6 +207,7 @@ typedef struct
   int atnr;
   t_functype *functype;
   t_iparams  *iparams;
+  real fudgeQQ;
 
   t_ilist il[F_NRE];
 } t_idef;

diff --git a/include/types/inputrec.h b/include/types/inputrec.h
index 6c641b64f7859d3274c789196a1373155ab302b1..b975f26b04750a25dba79ba5935f7a14feece249 100644 (file)
--- a/include/types/inputrec.h
+++ b/include/types/inputrec.h
@@ -184,7 +184,6 @@ typedef struct {
   real tabext;          /* Extension of the table beyond the cut-off,   *
                         * as well as the table length for 1-4 interac. */
   real shake_tol;      /* tolerance for shake                          */
-  real fudgeQQ;                /* Id. for 1-4 coulomb interactions             */
   int  efep;                   /* free energy interpolation no/yes             */
   real init_lambda;    /* initial value for perturbation variable      */
   real delta_lambda;   /* change of lambda per time step (1/dt)        */

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index 790f2c4d06c5f856d312e4737efbe30c1782b369..37efe262d0f327f033499517a093b3d87c42da5b 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -50,7 +50,7 @@ IF YOU'RE NOT SURE ABOUT WHAT YOU'RE DOING, DON'T DO IT!
 wall_density, wall_ewald_zfac)
 <li><A HREF="#pull"><b>COM pulling</b></A> (pull, ...)
 <li><A HREF="#nmr"><b>NMR refinement</b></A> (disre, disre_weighting, disre_mixed, disre_fc, disre_tau, nstdisreout, orire, orire_fc, orire_tau, orire_fitgrp, nstorireout)
-<li><A HREF="#free"><b>Free Energy calculations</b></A> (free_energy, init_lambda, delta_lambda, sc_alpha, sc_power, sc_sigma)
+<li><A HREF="#free"><b>Free Energy calculations</b></A> (free_energy, init_lambda, delta_lambda, sc_alpha, sc_power, sc_sigma, couple-moltype, couple-lambda0, couple-lambda1, couple-intramol)
 <li><A HREF="#neq"><b>Non-equilibrium MD</b></A> (acc_grps, accelerate, freezegrps, freezedim, cos_acceleration, deform)
 <li><A HREF="#ef"><b>Electric fields</b></A> (E_x, E_xt, E_y, E_yt, E_z, E_zt )
 <li><A HREF="#qmmm"><b>Mixed quantum/classical dynamics</b></A> (QMMM, QMMM-grps, QMMMscheme, QMmethod, QMbasis, QMcharge, Qmmult, CASorbitals, CASelectrons, SH)
@@ -1276,6 +1276,33 @@ only the values 1 and 2 are supported</dd>
 <dt><h4>sc_sigma: (0.3) [nm]</h4></dt>
 <dd>the soft-core sigma for particles which have a C6 or C12 parameter equal
 to zero</dd>
+<dt><h4>couple-moltype:</h4></dt>
+<dd>Here one can supply a molecule type (as defined in the topology)
+for calculating solvation or coupling free energies.
+<b>free_energy</b> has to be turned on.
+The Van der Waals interactions and/or charges in this molecule type can be
+turned on or off between lambda=0 and lambda=1, depending on the settings
+of <b>couple-lambda0</b> and <b>couple-lambda1</b>. If you want to decouple
+one of several copies of a molecule, you need to copy and rename
+the molecule definition in the topology.</dd>
+<dt><h4>couple-lambda0:</h4></dt>
+<dd><dl compact>
+<dt><b>vdw-q</b></dt>
+<dd>all interactions are on at lambda=0
+<dt><b>vdw</b></dt>
+<dd>the charges are zero (no Coulomb interactions) at lambda=0
+<dt><b>none</b></dt>
+<dd>the Van der Waals interactions are turned off and the charges are zero at lambda=0; soft-core interactions will be required to avoid singularities
+</dl>
+<dt><h4>couple-lambda1:</h4></dt>
+<dd> analogous to <b>couple-lambda1</b>, but for lambda=1
+<dt><h4>couple-intramol:</h4></dt>
+<dd><dl compact>
+<dt><b>no</b></dt>
+<dd>All intra-molecular non-bonded interactions for moleculetype <b>couple-moltype</b> are replaced by exclusions and explicit pair interactions. In this manner the decoupled state of the molecule corresponds to the proper vacuum state without periodicity effects.
+<dt><b>yes</b></dt>
+<dd>The intra-molecular Van der Waals and Coulomb interactions are also turned on/off. This can be useful for partitioning free-energies of relatively large molecules, where the intra-molecular non-bonded interactions might lead to kinetically trapped vacuum conformations.
+</dl>
 </dl>
 
 
@@ -1466,6 +1493,10 @@ reals to your subroutine. Check the inputrec definition in
 <A HREF="#bond">constraints</A><br>
 <A HREF="#neq">cos_acceleration</A><br>
 <A HREF="#el">coulombtype</A><br>
+<A HREF="#free">couple-intramol</A><br>
+<A HREF="#free">couple-lambda0</A><br>
+<A HREF="#free">couple-lambda1</A><br>
+<A HREF="#free">couple-moltype</A><br>
 <A HREF="#pp">define</A><br>
 <A HREF="#neq">deform</A><br>
 <A HREF="#free">delta_lambda</A><br>

diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c
index ee7755c8461aa37ae2aeaa2434841a3142a846fb..7634e242b21701220ad639bf91e2e0e2de263db6 100644 (file)
--- a/src/gmxlib/bondfree.c
+++ b/src/gmxlib/bondfree.c
@@ -1950,7 +1950,7 @@ void calc_bonds(FILE *fplog,const gmx_multisim_t *ms,
        ind = interaction_function[ftype].nrnb_ind;
        nat = interaction_function[ftype].nratoms+1;
        dvdl = 0;
-       if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_A) {
+       if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB) {
          v =
            interaction_function[ftype].ifunc(nbonds,idef->il[ftype].iatoms,
                                              idef->iparams,

diff --git a/src/gmxlib/ifunc.c b/src/gmxlib/ifunc.c
index 1d56123e989868730146a48bed493bd3af2f333d..66dbd70b3f170f642c35d26536f1f4637790856c 100644 (file)
--- a/src/gmxlib/ifunc.c
+++ b/src/gmxlib/ifunc.c
@@ -111,8 +111,8 @@ const t_interaction_function interaction_function[F_NRE]=
   def_bondedt ("TABDIHS", "Tab. Dih.",        4, 2, 2,  eNR_TABDIHS, tab_dihs     ),
   def_bondedz ("LJ14",     "LJ-14",           2, 2, 2,  eNR_NB14,   unimplemented ),
   def_nofc    ("COUL14",   "Coulomb-14"                                           ),
-  def_bondedz ("LJC14_A",  "LJC-14 A",        2, 2, 0,  eNR_NB14,   unimplemented ),
-  def_bondedz ("LJC_P_A",  "LJC Pairs A",     2, 0, 0,  eNR_NB14,   unimplemented ),
+  def_bondedz ("LJC14_Q",  "LJC-14 q",        2, 5, 0,  eNR_NB14,   unimplemented ),
+  def_bondedz ("LJC_NB",   "LJC Pairs NB",    2, 4, 0,  eNR_NB14,   unimplemented ),
   def_nb      ("LJ_SR",    "LJ (SR)",         2, 2                                ),
   def_nb      ("BHAM",     "Buck.ham (SR)",   2, 3                                ),
   def_nofc    ("LJ_LR",    "LJ (LR)"                                              ),
@@ -155,7 +155,8 @@ const t_interaction_function interaction_function[F_NRE]=
   def_nofc    ("TEMP",     "Temperature"      ),
   def_nofc    ("PRES",     "Pressure (bar)"   ),
   def_nofc    ("DV/DL",    "dVpot/dlambda"    ),
-  def_nofc    ("DK/DL",    "dEkin/dlambda"    )
+  def_nofc    ("DK/DL",    "dEkin/dlambda"    ),
+  def_nofc    ("DG/DL_CON","dG/dl constr."    )
 };
 
 bool have_interaction(t_idef *idef,int ftype)

diff --git a/src/gmxlib/names.c b/src/gmxlib/names.c
index 0522b3b7994e2afd382ab3b745bea6976913eb75..38e250282f761a75a6e527929da8025e5cac2b15 100644 (file)
--- a/src/gmxlib/names.c
+++ b/src/gmxlib/names.c
@@ -54,7 +54,7 @@ const char *ens_names[ensNR+1]=
 
 const char *ei_names[eiNR+1]=
 {
-  "md", "steep", "cg", "bd", "sd2", "nm", "l-bfgs", "tpi", "tpic", "sd", NULL 
+  "md", "steep", "cg", "bd", "sd", "nm", "l-bfgs", "tpi", "tpic", "sd1", NULL 
 };
 
 const char *bool_names[BOOL_NR+1]=

diff --git a/src/gmxlib/nonbonded/nonbonded.c b/src/gmxlib/nonbonded/nonbonded.c
index b69cb84d644bede843b32325b0aaeb067f8bae04..b96486f5559510e36eb43d34fbdb1228e3da2a40 100644 (file)
--- a/src/gmxlib/nonbonded/nonbonded.c
+++ b/src/gmxlib/nonbonded/nonbonded.c
@@ -527,7 +527,7 @@ do_nonbonded14(int ftype,int nbonds,
     rvec      dx,x14[2],f14[2];
     int       i,ai,aj,itype;
     int       typeA[2]={0,0},typeB[2]={0,1};
-    real      chargeA[2],chargeB[2];
+    real      chargeA[2]={0,0},chargeB[2];
     int       gid,shift_vir,shift_f;
     int       j_index[] = { 0, 1 };
     int       i0=0,i1=1,i2=2;
@@ -536,7 +536,7 @@ do_nonbonded14(int ftype,int nbonds,
     int       nthreads = 1;
     int       count;
     real      krf,crf,tabscale;
-    int       ntype;
+    int       ntype=0;
     real      *nbfp=NULL;
     real      *egnb=NULL,*egcoul=NULL;
     t_nblist  tmplist;
@@ -558,16 +558,15 @@ do_nonbonded14(int ftype,int nbonds,
 
     switch (ftype) {
     case F_LJ14:
-    case F_LJC14_A:
-      ntype = 1;
+    case F_LJC14_Q:
       eps = fr->epsfac*fr->fudgeQQ;
+      ntype  = 1;
       egnb   = gener->ee[egLJ14];
       egcoul = gener->ee[egCOUL14];
       break;
-    case F_LJC_PAIRS_A:
-      ntype = fr->ntype;
-      nbfp = fr->nbfp;
+    case F_LJC_PAIRS_NB:
       eps = fr->epsfac;
+      ntype  = 1;
       egnb   = gener->ee[egLJSR];
       egcoul = gener->ee[egCOULSR];
       break;
@@ -581,7 +580,7 @@ do_nonbonded14(int ftype,int nbonds,
 
     krf = fr->k_rf;
     crf = fr->c_rf;
-    
+
     /* Determine the values for icoul/ivdw. */
     if (fr->bEwald) {
       icoul = 1;
@@ -632,13 +631,20 @@ do_nonbonded14(int ftype,int nbonds,
             ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
              iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
              iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
-       case F_LJC14_A:
+         chargeA[0] = md->chargeA[ai];
+         chargeA[1] = md->chargeA[aj];
          nbfp = (real *)&(iparams[itype].lj14.c6A);
          break;
-       case F_LJC_PAIRS_A:
-         /* We could also consider doing these with the real atom numbers */
-         typeA[0] = md->typeA[ai];
-         typeA[1] = md->typeA[aj];
+       case F_LJC14_Q:
+         eps = fr->epsfac*iparams[itype].ljc14.fqq;
+         chargeA[0] = iparams[itype].ljc14.qi;
+         chargeA[1] = iparams[itype].ljc14.qj;
+         nbfp = (real *)&(iparams[itype].ljc14.c6);
+         break;
+       case F_LJC_PAIRS_NB:
+         chargeA[0] = iparams[itype].ljcnb.qi;
+         chargeA[1] = iparams[itype].ljcnb.qj;
+         nbfp = (real *)&(iparams[itype].ljcnb.c6);
          break;
        }
 
@@ -675,8 +681,6 @@ do_nonbonded14(int ftype,int nbonds,
         }
         else 
         {
-            chargeA[0] = md->chargeA[ai];
-            chargeA[1] = md->chargeA[aj];
             copy_rvec(x[ai],x14[0]);
             copy_rvec(x[aj],x14[1]);
             clear_rvec(f14[0]);
@@ -720,7 +724,7 @@ do_nonbonded14(int ftype,int nbonds,
                                           eps,
                                           krf,
                                           crf,
-                                                                                 fr->ewaldcoeff,
+                                         fr->ewaldcoeff,
                                           egcoul,
                                           typeA,
                                           typeB,
@@ -732,7 +736,7 @@ do_nonbonded14(int ftype,int nbonds,
                                           lambda,
                                           dvdlambda,
                                           fr->sc_alpha,
-                                                                                 fr->sc_power,
+                                         fr->sc_power,
                                           fr->sc_sigma6,
                                           &outeriter,
                                           &inneriter);

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index ab959563d07a66fe370bc3f8ce23843266a45432..b95597b49c6a108c005772fea84914c82314d9fb 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -63,7 +63,7 @@
 #endif
 
 /* This number should be increased whenever the file format changes! */
-static const int tpx_version = 53;
+static const int tpx_version = 54;
 
 /* This number should only be increased when you edit the TOPOLOGY section
  * of the tpx format. This way we can maintain forward compatibility too
@@ -74,7 +74,7 @@ static const int tpx_version = 53;
  * to the end of the tpx file, so we can just skip it if we only
  * want the topology.
  */
-static const int tpx_generation = 15;
+static const int tpx_generation = 16;
 
 /* This number should be the most recent backwards incompatible version 
  * I.e., if this number is 9, we cannot read tpx version 9 with this code.
@@ -132,8 +132,8 @@ static const t_ftupd ftupd[] = {
   { 26, F_FOURDIHS          },
   { 26, F_PIDIHS            },
   { 43, F_TABDIHS           },
-  { 41, F_LJC14_A           },
-  { 41, F_LJC_PAIRS_A       },
+  { 41, F_LJC14_Q           },
+  { 41, F_LJC_PAIRS_NB      },
   { 32, F_BHAM_LR           },
   { 32, F_RF_EXCL           },
   { 32, F_COUL_RECIP        },
@@ -149,7 +149,8 @@ static const t_ftupd ftupd[] = {
   { 50, F_VSITEN            },
   { 46, F_COM_PULL          },
   { 20, F_EQM               },
-  { 46, F_ECONSERVED        }
+  { 46, F_ECONSERVED        },
+  { 54, F_DGDL_CON          }
 };
 #define NFTUPD asize(ftupd)
 
@@ -231,7 +232,8 @@ static void do_pull(t_pull *pull,bool bRead, int file_version)
     do_pullgrp(&pull->grp[g],bRead,file_version);
 }
 
-static void do_inputrec(t_inputrec *ir,bool bRead, int file_version)
+static void do_inputrec(t_inputrec *ir,bool bRead, int file_version,
+                       real *fudgeQQ)
 {
   int  i,j,k,*tmp,idum=0; 
   bool bDum=TRUE;
@@ -457,7 +459,8 @@ static void do_inputrec(t_inputrec *ir,bool bRead, int file_version)
       do_real(rdum); 
 
     do_real(ir->shake_tol);
-    do_real(ir->fudgeQQ);
+    if (file_version < 54)
+      do_real(*fudgeQQ);
     do_int(ir->efep);
     if (file_version <= 14 && ir->efep > efepNO)
       ir->efep = efepYES;
@@ -834,11 +837,18 @@ void do_iparams(t_functype ftype,t_iparams *iparams,bool bRead, int file_version
     do_real(iparams->lj14.c6B);
     do_real(iparams->lj14.c12B);
     break;
-  case F_LJC14_A:
-    do_real(iparams->lj14.c6A);
-    do_real(iparams->lj14.c12A);
+  case F_LJC14_Q:
+    do_real(iparams->ljc14.fqq);
+    do_real(iparams->ljc14.qi);
+    do_real(iparams->ljc14.qj);
+    do_real(iparams->ljc14.c6);
+    do_real(iparams->ljc14.c12);
     break;
-  case F_LJC_PAIRS_A:
+  case F_LJC_PAIRS_NB:
+    do_real(iparams->ljcnb.qi);
+    do_real(iparams->ljcnb.qj);
+    do_real(iparams->ljcnb.c6);
+    do_real(iparams->ljcnb.c12);
     break;
   case F_PDIHS:
   case F_PIDIHS:
@@ -1004,6 +1014,8 @@ static void do_idef(t_idef *idef,bool bRead, int file_version)
     if (bRead && debug)
       pr_iparams(debug,idef->functype[i],&idef->iparams[i]);
   }
+  if (file_version >= 54)
+    do_real(idef->fudgeQQ);
   
   for(j=0; (j<F_NRE); j++) {
     bClear = FALSE;
@@ -1488,12 +1500,12 @@ static int do_tpx(int fp,bool bRead,int *step,real *t,
     do_section(eitemIR,bRead);
     if (tpx.bIr) {
       if (ir) {
-       do_inputrec(ir,bRead,file_version);
+       do_inputrec(ir,bRead,file_version,top ? &top->idef.fudgeQQ : NULL);
        if (bRead && debug) 
          pr_inputrec(debug,0,"inputrec",ir,FALSE);
       }
       else {
-       do_inputrec  (&dum_ir,bRead,file_version);
+       do_inputrec(&dum_ir,bRead,file_version,top ? &top->idef.fudgeQQ :NULL);
        if (bRead && debug) 
          pr_inputrec(debug,0,"inputrec",&dum_ir,FALSE);
        done_inputrec(&dum_ir);
@@ -1561,7 +1573,7 @@ static int do_tpx(int fp,bool bRead,int *step,real *t,
        do_int(bPeriodicMols);
       }
       if (file_generation <= tpx_generation && ir) {
-       do_inputrec(ir,bRead,file_version);
+       do_inputrec(ir,bRead,file_version,top ? &top->idef.fudgeQQ : NULL);
        if (bRead && debug) 
          pr_inputrec(debug,0,"inputrec",ir,FALSE);
        if (file_version < 51)

diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c
index 53553f372a3b779dd0ee6816b7e6f76b03d062a8..7c3902bdcb250a9b5635f20d676ac6382295ca73 100644 (file)
--- a/src/gmxlib/txtdump.c
+++ b/src/gmxlib/txtdump.c
@@ -513,7 +513,6 @@ void pr_inputrec(FILE *fp,int indent,const char *title,t_inputrec *ir,
     PR("gb_saltconc",ir->gb_saltconc);
     PS("implicit_solvent",EIMPLICITSOL(ir->implicit_solvent));
     PS("DispCorr",EDISPCORR(ir->eDispCorr));
-    PR("fudgeQQ",ir->fudgeQQ);
     PS("free_energy",EFEPTYPE(ir->efep));
     PR("init_lambda",ir->init_lambda);
     PR("sc_alpha",ir->sc_alpha);
@@ -679,12 +678,16 @@ void pr_iparams(FILE *fp,t_functype ftype,t_iparams *iparams)
            iparams->lj14.c6A,iparams->lj14.c12A,
            iparams->lj14.c6B,iparams->lj14.c12B);
     break;
-  case F_LJC14_A:
-    fprintf(fp,"c6=%15.8e, c12=%15.8e\n",
-           iparams->lj14.c6A,iparams->lj14.c12A);
+  case F_LJC14_Q:
+    fprintf(fp,"fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n",
+           iparams->ljc14.fqq,
+           iparams->ljc14.qi,iparams->ljc14.qj,
+           iparams->ljc14.c6,iparams->ljc14.c12);
     break;
-  case F_LJC_PAIRS_A:
-    fprintf(fp,"\n");
+  case F_LJC_PAIRS_NB:
+    fprintf(fp,"qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n",
+           iparams->ljcnb.qi,iparams->ljcnb.qj,
+           iparams->ljcnb.c6,iparams->ljcnb.c12);
     break;
   case F_PDIHS:
   case F_ANGRES:
@@ -839,6 +842,8 @@ void pr_idef(FILE *fp,int indent,const char *title,t_idef *idef, bool bShowNumbe
                     interaction_function[idef->functype[i]].name);
       pr_iparams(fp,idef->functype[i],&idef->iparams[i]);
     }
+    (void) pr_real(fp,indent,"fudgeQQ",idef->fudgeQQ);
+
     for(j=0; (j<F_NRE); j++)
       pr_ilist(fp,indent,interaction_function[j].longname,
               idef,&idef->il[j],bShowNumbers);

diff --git a/src/kernel/convparm.c b/src/kernel/convparm.c
index 81599eded07802699b384e0b3d3ed615347f46eb..0d2a79c89ee30c0976c00924f41d7714a8f904ea 100644 (file)
--- a/src/kernel/convparm.c
+++ b/src/kernel/convparm.c
@@ -71,6 +71,23 @@ static int round_check(real r,int limit,int ftype,char *name)
   return i;
 }
 
+static void set_ljparams(int comb,real reppow,real v,real w,real *c6,real *c12)
+{
+  if (comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) {
+    if (v >= 0) {
+      *c6  = 4*w*pow(v,6.0);
+      *c12 = 4*w*pow(v,reppow);
+    } else {
+      /* Interpret negative sigma as c6=0 and c12 with -sigma */
+      *c6  = 0;
+      *c12 = 4*w*pow(-v,reppow);
+    }
+  } else {
+    *c6  = v;
+    *c12 = w;
+  }
+}
+
 static void assign_param(t_functype ftype,t_iparams *new,
                         real old[MAXFORCEPARAM],int comb,real reppow)
 {
@@ -176,48 +193,22 @@ static void assign_param(t_functype ftype,t_iparams *new,
     new->bham.c = old[2];
     break;
   case F_LJ14:
-  case F_LJC14_A:
-    if (comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) {
-      if (old[0] >= 0) {
-       new->lj14.c6A  = 4*old[1]*pow(old[0],6.0);
-       new->lj14.c12A = 4*old[1]*pow(old[0],reppow);
-      } else {
-       new->lj14.c6A  = 0;
-       new->lj14.c12A = 4*old[1]*pow(-old[0],reppow);
-      }
-      if (ftype == F_LJ14) {
-       if (old[2] >= 0) {
-         new->lj14.c6B  = 4*old[3]*pow(old[2],6.0);
-         new->lj14.c12B = 4*old[3]*pow(old[2],reppow);
-       } else {
-         new->lj14.c6B  = 0;
-         new->lj14.c12B = 4*old[3]*pow(-old[2],reppow);
-       }
-      }
-    } else {
-      new->lj14.c6A  = old[0]; 
-      new->lj14.c12A = old[1];
-      if (ftype == F_LJ14) {
-       new->lj14.c6B  = old[2]; 
-       new->lj14.c12B = old[3];
-      }
-    }
+    set_ljparams(comb,reppow,old[0],old[1],&new->lj14.c6A,&new->lj14.c12A);
+    set_ljparams(comb,reppow,old[2],old[3],&new->lj14.c6B,&new->lj14.c12B);
+    break;
+  case F_LJC14_Q:
+    new->ljc14.fqq = old[0];
+    new->ljc14.qi  = old[1];
+    new->ljc14.qj  = old[2];
+    set_ljparams(comb,reppow,old[3],old[4],&new->ljc14.c6,&new->ljc14.c12);
     break;
-  case F_LJC_PAIRS_A:
+  case F_LJC_PAIRS_NB:
+    new->ljcnb.qi = old[0];
+    new->ljcnb.qj = old[1];
+    set_ljparams(comb,reppow,old[2],old[3],&new->ljcnb.c6,&new->ljcnb.c12);
     break;
   case F_LJ:
-    if (comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) {
-      if (old[0] >= 0) {
-       new->lj.c6  = 4*old[1]*pow(old[0],6);
-       new->lj.c12 = 4*old[1]*pow(old[0],reppow);
-      } else {
-       new->lj.c6  = 0;
-       new->lj.c12 = 4*old[1]*pow(-old[0],reppow);
-      }
-    } else {
-      new->lj.c6  = old[0]; 
-      new->lj.c12 = old[1];
-    }
+    set_ljparams(comb,reppow,old[0],old[1],&new->lj.c6,&new->lj.c12);
     break;
   case F_PDIHS:
   case F_ANGRES:
@@ -410,7 +401,8 @@ static void new_interaction_list(t_ilist *ilist)
 }
 
 void convert_params(int atnr,t_params nbtypes[],
-                   t_params plist[],int comb,real reppow,t_idef *idef)
+                   t_params plist[],int comb,real reppow,real fudgeQQ,
+                   t_idef *idef)
 {
   int    i,j,maxtypes;
   unsigned long  flags;
@@ -448,4 +440,5 @@ void convert_params(int atnr,t_params nbtypes[],
       printf("# %10s:   %d\n",
             interaction_function[j].name,idef->il[j].nr/(1+NRAL(j)));
   }
+  idef->fudgeQQ = fudgeQQ;
 }

diff --git a/src/kernel/convparm.h b/src/kernel/convparm.h
index 4b8f6c82200b011394932e521d61a73ee8a72471..90c3e3e6c2f098dc356ca424fcfaf465596b7bcb 100644 (file)
--- a/src/kernel/convparm.h
+++ b/src/kernel/convparm.h
@@ -41,6 +41,7 @@
 
 extern void convert_params(int atnr,t_params plist[],
                           t_params nbtypes[],int comb,real reppow,
+                          real fudgeQQ,
                           t_idef *idef);
 
 #endif /* _convparm_h */

diff --git a/src/kernel/grompp.c b/src/kernel/grompp.c
index 9f8d063ccd049267ca46d985abe69ebb53ecc51f..7312bfc9eef9293881f3970858eaf82dc920fe56 100644 (file)
--- a/src/kernel/grompp.c
+++ b/src/kernel/grompp.c
@@ -194,7 +194,8 @@ new_status(char *topfile,char *topppfile,char *confin,
           t_gromppopts *opts,t_inputrec *ir,bool bZero,
           bool bGenVel,bool bVerbose,t_state *state,
           t_atomtype atype,t_topology *sys,
-          t_molinfo *msys,t_params plist[],int *comb,real *reppow,
+          t_molinfo *msys,t_params plist[],
+          int *comb,real *reppow,real *fudgeQQ,
           bool bEnsemble,bool bMorse,
           int *nerror)
 {
@@ -209,7 +210,8 @@ new_status(char *topfile,char *topppfile,char *confin,
   
   /* TOPOLOGY processing */
   msys->name=do_top(bVerbose,topfile,topppfile,opts,bZero,&(sys->symtab),
-                   plist,comb,reppow,atype,&nrmols,&molinfo,ir,&Nsim,&Sims);
+                   plist,comb,reppow,fudgeQQ,
+                   atype,&nrmols,&molinfo,ir,&Nsim,&Sims);
   
   if (bMorse)
     convert_harmonics(nrmols,molinfo,atype);
@@ -570,7 +572,7 @@ int main (int argc, char *argv[])
   t_params     *plist;
   t_state      state;
   matrix       box;
-  real         max_spacing,reppow;
+  real         max_spacing,reppow,fudgeQQ;
   char         fn[STRLEN],fnB[STRLEN],*mdparin;
   int          nerror,ntype;
   bool         bNeedVel,bGenVel;
@@ -655,7 +657,7 @@ int main (int argc, char *argv[])
     gmx_fatal(FARGS,"%s does not exist",fn);
   new_status(fn,opt2fn_null("-pp",NFILE,fnm),opt2fn("-c",NFILE,fnm),
             opts,ir,bZero,bGenVel,bVerbose,&state,
-            atype,sys,&msys,plist,&comb,&reppow,
+            atype,sys,&msys,plist,&comb,&reppow,&fudgeQQ,
             (opts->eDisre==edrEnsemble),opts->bMorse,
             &nerror);
   
@@ -767,7 +769,7 @@ int main (int argc, char *argv[])
 
   if (bVerbose) 
     fprintf(stderr,"converting bonded parameters...\n");
-  convert_params(ntype, plist, msys.plist, comb, reppow, &sys->idef);
+  convert_params(ntype, plist, msys.plist, comb, reppow, fudgeQQ, &sys->idef);
   
   if (debug)
     pr_symtab(debug,0,"After convert_params",&sys->symtab);

diff --git a/src/kernel/md.c b/src/kernel/md.c
index eb6258b71db7ae487bb7f90f4275e27eb3694aac..705b5b11d91cdbdad927ad275bd2ca0d94b8498c 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -960,7 +960,7 @@ time_t do_md(FILE *fplog,t_commrec *cr,int nfile,t_filenm fnm[],
             bNEMD,bDoBerendsenCoupl,bFirstStep,pres);
       if (fr->bSepDVDL && fplog && do_log)
        fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
-      ener[F_DVDL] += dvdl;
+      ener[F_DGDL_CON] += dvdl;
       wallcycle_stop(wcycle,ewcUPDATE);
     } else if (graph) {
       /* Need to unshift here */
@@ -1107,7 +1107,7 @@ time_t do_md(FILE *fplog,t_commrec *cr,int nfile,t_filenm fnm[],
 
     /* Sum the kinetic energies of the groups & calc temp */
     ener[F_TEMP] = sum_ekin(bRerunMD,&(ir->opts),grps,ekin,
-                           &(ener[F_DVDLKIN]));
+                           &(ener[F_DKDL]));
     ener[F_EKIN] = trace(ekin);
 
     /* Calculate Temperature coupling parameters lambda and adjust

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index 001711e4969a6c2d04a47eff77bde8bc6bfbc07c..c8573cb274f236ef94b73db304d434c1626b6a66 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -71,7 +71,7 @@
 static char tcgrps[STRLEN],tau_t[STRLEN],ref_t[STRLEN],
   acc[STRLEN],accgrps[STRLEN],freeze[STRLEN],frdim[STRLEN],
   energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
-  orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
+  couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
   wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
 static char **pull_grp;
 static char anneal[STRLEN],anneal_npoints[STRLEN],
@@ -756,6 +756,10 @@ void get_ir(char *mdparin,char *mdparout,
   RTYPE ("sc-alpha",ir->sc_alpha,0.0);
   ITYPE ("sc-power",ir->sc_power,0);
   RTYPE ("sc-sigma",ir->sc_sigma,0.3);
+  STYPE ("couple-moltype",  couple_moltype,  NULL);
+  EETYPE("couple-lambda0", opts->couple_lam0, couple_lam, nerror, TRUE);
+  EETYPE("couple-lambda1", opts->couple_lam1, couple_lam, nerror, TRUE);
+  EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names, nerror, TRUE);
 
   /* Non-equilibrium MD stuff */  
   CCTYPE("Non-equilibrium MD stuff");
@@ -864,6 +868,19 @@ void get_ir(char *mdparin,char *mdparout,
   if (ir->comm_mode == ecmNO)
     ir->nstcomm = 0;
 
+  opts->couple_moltype = NULL;
+  if (strlen(couple_moltype) > 0) {
+    if (ir->efep != efepNO) {
+      opts->couple_moltype = strdup(couple_moltype);
+      if (opts->couple_lam0 == opts->couple_lam1)
+       warning("The lambda=0 and lambda=1 states for coupling are identical");
+      if (ir->eI == eiMD)
+       warning("For proper sampling with coupling, stochastic dynamics should be used");
+    } else {
+      warning("Can not couple a molecule with free_energy = no");
+    }
+  }
+
   do_wall_params(ir,wall_atomtype,wall_density,opts);
   
   if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {

diff --git a/src/kernel/readir.h b/src/kernel/readir.h
index f9b1d6d8b3d9b0706e057c8c761566d4621a4aac..9772ab80d7dad5e67134fea24bb4bb7805d3b7b6 100644 (file)
--- a/src/kernel/readir.h
+++ b/src/kernel/readir.h
@@ -45,15 +45,19 @@
 enum { eshNONE, eshHBONDS, eshALLBONDS, eshHANGLES, eshALLANGLES, eshNR };
 
 static const char *constraints[eshNR+1]    = { 
-    "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL 
-  };
+  "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL 
+};
+
+enum { ecouplamVDWQ, ecouplamVDW, ecouplamNONE, ecouplamNR };
+
+static const char *couple_lam[ecouplamNR+1]    = { 
+  "vdw-q", "vdw", "none", NULL 
+};
 
 typedef struct {
   int warnings;
   int nshake;
   real fourierspacing;
-  int nprocs;
-  int splitalg;
   char *include;
   char *define;
   bool bGenVel;
@@ -66,6 +70,10 @@ typedef struct {
   bool bMorse;
   char *wall_atomtype[2];
   bool pull_start;
+  char *couple_moltype;
+  int  couple_lam0;
+  int  couple_lam1;
+  bool bCoupleIntra;
 } t_gromppopts;
 
 

diff --git a/src/kernel/topdirs.c b/src/kernel/topdirs.c
index 89df9a37f0bf6df68757be2fc7cc1e187482cd4b..fe0bc5a63cae1cea5ca41c254a80e279f458016d 100644 (file)
--- a/src/kernel/topdirs.c
+++ b/src/kernel/topdirs.c
@@ -94,13 +94,11 @@ int ifunc_index(directive d,int type)
     if (type == 1 || (d == d_pairtypes && type == 2))
       return F_LJ14;
     else if (type == 2)
-      return F_LJC14_A;
-    else if (type == 3) {
-      if (d == d_pairtypes)
-       gmx_fatal(FARGS,"Can not have parameters for pair type %d",type);
-      return F_LJC_PAIRS_A;
-    } else
-      gmx_fatal(FARGS,"Invalid pair type %d",type);
+      return F_LJC14_Q;
+    else
+      gmx_fatal(FARGS,"Invalid pairs type %d",type);
+  case d_pairs_nb:
+    return F_LJC_PAIRS_NB;
   case d_dihedrals:
   case d_dihedraltypes:
     switch (type) {
@@ -264,6 +262,7 @@ void DS_Init(DirStack **DS)
     set_nec(&(necessary[d_bonds]),d_atoms,d_none);
     set_nec(&(necessary[d_exclusions]),d_bonds,d_constraints,d_settles,d_none);
     set_nec(&(necessary[d_pairs]),d_atoms,d_none);
+    set_nec(&(necessary[d_pairs_nb]),d_atoms,d_none);
     set_nec(&(necessary[d_angles]),d_atoms,d_none);
     set_nec(&(necessary[d_polarization]),d_atoms,d_none);
     set_nec(&(necessary[d_water_polarization]),d_atoms,d_none);

diff --git a/src/kernel/topexcl.h b/src/kernel/topexcl.h
index d62b5ff33e87ba544b94ba8aab98bc959c719db2..09febce62cfbcc7547176db2c850b45c634de68b 100644 (file)
--- a/src/kernel/topexcl.h
+++ b/src/kernel/topexcl.h
@@ -78,4 +78,5 @@ extern void generate_excl (int nrexcl, int nratoms,
 /* Generate an exclusion block from bonds and constraints in
  * plist.
  */
+
 #endif /* _topexcl_h */

diff --git a/src/kernel/topio.c b/src/kernel/topio.c
index 6741ba46d12af9845e9a3288f0d599b7f4c2ce1c..09214059ac8fcb1bd979bab4e238b28c2408e9f4 100644 (file)
--- a/src/kernel/topio.c
+++ b/src/kernel/topio.c
@@ -286,7 +286,7 @@ static char **read_topol(char *infile,char *outfile,
                         t_params    plist[],
                         int         *combination_rule,
                         real        *reppow,
-                        int         nshake,
+                        t_gromppopts *opts,
                         real        *fudgeQQ,
                         int         *nsim,
                         t_simsystem **sims,
@@ -314,6 +314,7 @@ static char **read_topol(char *infile,char *outfile,
   double     qt=0,qBt=0; /* total charge */
   t_bond_atomtype batype;
   int        lastcg=-1;
+  int        dcatt=-1,nmol_couple;
   /* File handling variables */
   int        status,done;
   gmx_cpp_t  handle;
@@ -345,6 +346,7 @@ static char **read_topol(char *infile,char *outfile,
   bReadDefaults = FALSE;
   bGenPairs     = FALSE;
   bReadMolType  = FALSE;
+  nmol_couple = 0;
   
   do {
     status = cpp_read_line(&handle,STRLEN,line);
@@ -485,7 +487,13 @@ static char **read_topol(char *infile,char *outfile,
            */
          case d_moleculetype: {
            if (!bReadMolType) {
-             int ntype = get_atomtype_ntypes(atype);
+             int ntype;
+             if (opts->couple_moltype &&
+                 (opts->couple_lam0 == ecouplamNONE ||
+                  opts->couple_lam1 == ecouplamNONE))
+               dcatt = add_atomtype_decoupled(symtab,atype,
+                                              &nbparam,bGenPairs?&pair:NULL);
+             ntype = get_atomtype_ntypes(atype);
              ncombs = (ntype*(ntype+1))/2;
              generate_nbparams(comb,nb_funct,&(plist[nb_funct]),atype);
              ncopy = copy_nbparams(nbparam,nb_funct,&(plist[nb_funct]),
@@ -516,7 +524,7 @@ static char **read_topol(char *infile,char *outfile,
            
          case d_pairs: 
            push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,FALSE,
-                     bGenPairs,bZero,&bWarn_copy_A_B);
+                     bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B);
            break;
            
          case d_vsites2:
@@ -537,7 +545,7 @@ static char **read_topol(char *infile,char *outfile,
          case d_water_polarization:
          case d_thole_polarization:
            push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,TRUE,
-                     bGenPairs,bZero,&bWarn_copy_A_B);
+                     bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B);
            break;
          case d_vsitesn:
            push_vsitesn(d,plist,mi0->plist,&(mi0->atoms),atype,pline);
@@ -552,7 +560,8 @@ static char **read_topol(char *infile,char *outfile,
            title=put_symtab(symtab,pline);
            break;
          case d_molecules: {
-           int whichmol;
+           int  whichmol;
+           bool bCouple;
            
            push_mol(nmol,*molinfo,pline,&whichmol,&nrcopies);
            mi0=&((*molinfo)[whichmol]);
@@ -560,6 +569,12 @@ static char **read_topol(char *infile,char *outfile,
            Sims[Nsim].whichmol=whichmol;
            Sims[Nsim].nrcopies=nrcopies;
            Nsim++;
+           
+           bCouple = (opts->couple_moltype &&
+                      strcmp(*(mi0->name),opts->couple_moltype) == 0);
+           if (bCouple)
+             nmol_couple += nrcopies;
+
            if (mi0->atoms.nr == 0)
              gmx_fatal(FARGS,"Moleculetype %s contains no atoms",*mi0->name);
            fprintf(stderr,"Excluding %d bonded neighbours for %s\n",
@@ -572,7 +587,13 @@ static char **read_topol(char *infile,char *outfile,
                            &(mi0->excls));
              merge_excl(&(mi0->excls),&(block2[whichmol]));
              done_block2(&(block2[whichmol]));
-             make_shake(mi0->plist,&mi0->atoms,atype,nshake); 
+             make_shake(mi0->plist,&mi0->atoms,atype,opts->nshake);
+             if (bCouple) {
+               convert_moltype_couple(mi0,dcatt,*fudgeQQ,
+                                      opts->couple_lam0,opts->couple_lam1,
+                                      opts->bCoupleIntra,
+                                      nb_funct,&(plist[nb_funct]));
+             }
              stupid_fill_block(&mi0->mols,mi0->atoms.nr,TRUE);
              mi0->bProcessed=TRUE;
            }
@@ -593,6 +614,16 @@ static char **read_topol(char *infile,char *outfile,
     gmx_fatal(FARGS,cpp_error(&handle,status));
   if (out)
     fclose(out);
+
+  if (opts->couple_moltype) {
+    if (nmol_couple == 0) {
+      gmx_fatal(FARGS,"Did not find any molecules of type '%s' for coupling",
+               opts->couple_moltype);
+    } else {
+      fprintf(stderr,"Found %d copies of molecule type '%s' for coupling\n",
+             nmol_couple,opts->couple_moltype);
+    }
+  }
   
   /* this is not very clean, but fixes core dump on empty system name */
   if(!title)
@@ -625,6 +656,7 @@ char **do_top(bool         bVerbose,
              t_params     plist[],
              int          *combination_rule,
              real         *repulsion_power,
+             real         *fudgeQQ,
              t_atomtype   atype,
              int          *nrmols,
              t_molinfo    **molinfo,
@@ -645,8 +677,8 @@ char **do_top(bool         bVerbose,
   title=read_topol(topfile,tmpfile,opts->define,opts->include,
                   symtab,atype,nrmols,molinfo,
                   plist,combination_rule,repulsion_power,
-                  opts->nshake,&ir->fudgeQQ,nsim,sims,ir->efep!=efepNO,
-                  bZero,bVerbose);
+                  opts,fudgeQQ,nsim,sims,
+                  ir->efep!=efepNO,bZero,bVerbose);
   if ((*combination_rule != eCOMB_GEOMETRIC) && 
       (ir->vdwtype == evdwUSER)) {
     warning("Using sigma/epsilon based combination rules with"

diff --git a/src/kernel/topio.h b/src/kernel/topio.h
index 925320886be05be502576f602b978627dfbbeaba..2f81ac36e0268771025fcca04fab5d83e4346cc1 100644 (file)
--- a/src/kernel/topio.h
+++ b/src/kernel/topio.h
@@ -59,6 +59,7 @@ extern char **do_top(bool         bVerbose,
                     t_params     plist[],
                     int          *combination_rule,
                     real         *repulsion_power,
+                    real         *fudgeQQ,
                     t_atomtype   atype,
                     int          *nrmols,
                     t_molinfo    **molinfo,

diff --git a/src/kernel/toppush.c b/src/kernel/toppush.c
index 2e0d1a2a1a6a37f4e651d62b65041e52b09a94c9..2111d8264b0ab3cd33ee44b129ef344806880274 100644 (file)
--- a/src/kernel/toppush.c
+++ b/src/kernel/toppush.c
@@ -48,6 +48,7 @@
 #include "toputil.h"
 #include "toppush.h"
 #include "topdirs.h"
+#include "readir.h"
 #include "symtab.h"
 #include "gmx_fatal.h"
 #include "gpp_atomtype.h"
@@ -143,6 +144,19 @@ void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype atype)
   }
 }
 
+static void realloc_nb_params(t_atomtype at,
+                             t_nbparam ***nbparam,t_nbparam ***pair)
+{
+  /* Add space in the non-bonded parameters matrix */
+  int atnr = get_atomtype_ntypes(at);
+  srenew(*nbparam,atnr);
+  snew((*nbparam)[atnr-1],atnr);
+  if (pair) {
+    srenew(*pair,atnr);
+    snew((*pair)[atnr-1],atnr);
+  }
+}
+
 void push_at (t_symtab *symtab, t_atomtype at, t_bond_atomtype bat,
              char *line,int nb_funct,
              t_nbparam ***nbparam,t_nbparam ***pair)
@@ -428,13 +442,7 @@ void push_at (t_symtab *symtab, t_atomtype at, t_bond_atomtype bat,
     gmx_fatal(FARGS,"Adding atomtype %s failed",type);
   else {  
     /* Add space in the non-bonded parameters matrix */
-    int atnr = get_atomtype_ntypes(at);
-    srenew(*nbparam,atnr);
-    snew((*nbparam)[atnr-1],atnr);
-    if (pair) {
-      srenew(*pair,atnr);
-      snew((*pair)[atnr-1],atnr);
-    }
+    realloc_nb_params(at,nbparam,pair);
   }
   sfree(atom);
   sfree(param);
@@ -729,7 +737,7 @@ void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
   
   /* Get the force parameters */
   nrfp = NRFP(ftype);
-  if (ftype == F_LJ14 || ftype == F_LJC14_A) {
+  if (ftype == F_LJ14) {
     n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
     if (n < 2) {
       too_few();
@@ -741,6 +749,13 @@ void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
     for(i=n; i<nrfp; i++)
       c[i] = c[i-2];
   }
+  else if (ftype == F_LJC14_Q) {
+    n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
+    if (n < 4) {
+      too_few();
+      return;
+    }
+  }
   else if (nrfp == 2) {
     if (sscanf(pline,form2,&c[0],&c[1]) != 2) {
       too_few();
@@ -945,7 +960,7 @@ void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line)
 }
 
 static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
-                             t_param *p,bool bB,bool bGenPairs)
+                             t_param *p,int c_start,bool bB,bool bGenPairs)
 {
   int      i,j,ti,tj,ntype;
   bool     bFound;
@@ -964,7 +979,7 @@ static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
      * time when we have 1000*1000 entries for e.g. OPLS...
      */
     ntype=sqrt(nr);
-    if(bB) {
+    if (bB) {
       ti=at->atom[p->a[0]].typeB;
       tj=at->atom[p->a[1]].typeB;
     } else {
@@ -994,15 +1009,15 @@ static bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
       if (nrfp+nrfpB > MAXFORCEPARAM) {
        gmx_incons("Too many force parameters");
       }
-      for (j=0; (j < nrfpB); j++)
+      for (j=c_start; (j < nrfpB); j++)
        p->c[nrfp+j] = pi->c[j];
     }
     else
-      for (j=0; (j < nrfp); j++)
+      for (j=c_start; (j < nrfp); j++)
        p->c[j] = pi->c[j];
   }
   else {
-    for (j=0; (j < nrfp); j++)
+    for (j=c_start; (j < nrfp); j++)
       p->c[j] = 0.0;
   }
   return bFound;
@@ -1108,7 +1123,8 @@ void push_bondnow(t_params *bond, t_param *b)
 
 void push_bond(directive d,t_params bondtype[],t_params bond[],
               t_atoms *at,t_atomtype atype,char *line,
-              bool bBonded,bool bGenPairs,bool bZero,bool *bWarn_copy_A_B)
+              bool bBonded,bool bGenPairs,real fudgeQQ,
+              bool bZero,bool *bWarn_copy_A_B)
 {
   const char *aaformat[MAXATOMLIST]= {
     "%d%d",
@@ -1131,7 +1147,7 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
   double   cc[MAXFORCEPARAM+1];
   int      aa[MAXATOMLIST+1];
   t_param  param,paramB,*param_def;
-  bool     bFoundA,bFoundB,bDef,bPert,bSwapParity=FALSE;
+  bool     bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
   char  errbuf[256];
 
   ftype = ifunc_index(d,1);
@@ -1199,7 +1215,7 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
   /* Get force params for normal and free energy perturbation
    * studies, as determined by types!
    */
-  if(bBonded) {
+  if (bBonded) {
     bFoundA = default_params(ftype,bondtype,at,atype,&param,FALSE,&param_def);
     if (bFoundA) {
       /* Copy the A-state and B-state default parameters */
@@ -1212,10 +1228,23 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
       for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
        param.c[j] = param_def->c[j];
     }
+  } else if (ftype == F_LJ14) {
+    bFoundA = default_nb_params(ftype, bondtype,at,&param,0,FALSE,bGenPairs);
+    bFoundB = default_nb_params(ftype, bondtype,at,&param,0,TRUE, bGenPairs);
+  } else if (ftype == F_LJC14_Q) {
+    param.c[0] = fudgeQQ;
+    /* Fill in the A-state charges as default parameters */
+    param.c[1] = at->atom[param.a[0]].q;
+    param.c[2] = at->atom[param.a[1]].q;
+    /* The default LJ parameters are the standard 1-4 parameters */
+    bFoundA = default_nb_params(F_LJ14,bondtype,at,&param,3,FALSE,bGenPairs);
+    bFoundB = TRUE;
+  } else if (ftype == F_LJC_PAIRS_NB) {
+    /* Defaults are not supported here */
+    bFoundA = FALSE;
+    bFoundB = TRUE;
   } else {
-    bFoundA = default_nb_params(ftype == F_LJC14_A ? F_LJ14 : ftype,
-                               bondtype,at,&param,FALSE,bGenPairs);
-    bFoundB = default_nb_params(ftype,bondtype,at,&param,TRUE,bGenPairs);
+    gmx_incons("Unknown function type in push_bond");
   }
   
   if (nread > nral) {  
@@ -1251,9 +1280,11 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
     /* If nread was 0 or EOF, no parameters were read => use defaults.
      * If nread was nrfpA we copied above so nread=nrfp.
      * If nread was nrfp we are cool.
+     * For F_LJC14_Q we allow supplying fudgeQQ only.
      * Anything else is an error!
      */        
-    if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)))
+    if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
+       !(ftype == F_LJC14_Q && nread == 1))
       {
        gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
                  nread,NRFPA(ftype),NRFP(ftype),
@@ -1449,7 +1480,7 @@ void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
 }
 
 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
-                 int *nrcopies)
+             int *nrcopies)
 {
   int  i,copies;
   char type[STRLEN];
@@ -1606,3 +1637,139 @@ void merge_excl(t_blocka *excl, t_block2 *b2)
   b2_to_b(b2,excl);
 }
 
+int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
+                          t_nbparam ***nbparam,t_nbparam ***pair)
+{
+  t_atom  atom;
+  t_param param;
+  int     i,nr;
+
+  /* Define an atom type with all parameters set to zero (no interactions) */
+  atom.q = 0.0;
+  atom.m = 0.0;
+  /* Type for decoupled atoms could be anything,
+   * this should be changed automatically later when required.
+   */
+  atom.ptype = eptAtom;
+  for (i=0; (i<MAXFORCEPARAM); i++)
+    param.c[i] = 0.0;
+
+  nr = add_atomtype(at,symtab,&atom,"decoupled",&param,-1,0.0,0.0,0.0,0);
+
+  /* Add space in the non-bonded parameters matrix */
+  realloc_nb_params(at,nbparam,pair);
+
+  return nr;
+}
+
+static void convert_pairs_to_pairsQ(t_params *plist,
+                                   real fudgeQQ,t_atoms *atoms)
+{
+  t_param *param;
+  int  i;
+  real v,w;
+
+  /* Copy the pair list to the pairQ list */
+  plist[F_LJC14_Q] = plist[F_LJ14];
+  /* Empty the pair list */
+  plist[F_LJ14].nr    = 0;
+  plist[F_LJ14].param = NULL;
+  param = plist[F_LJC14_Q].param;
+  for(i=0; i<plist[F_LJC14_Q].nr; i++) {
+    v = param[i].c[0];
+    w = param[i].c[1];
+    param[i].c[0] = fudgeQQ;
+    param[i].c[1] = atoms->atom[param[i].a[0]].q;
+    param[i].c[2] = atoms->atom[param[i].a[1]].q;
+    param[i].c[3] = v;
+    param[i].c[4] = w;
+  }
+}
+
+static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
+{
+  int n,ntype,i,j,k;
+  t_atom *atom;
+  t_blocka *excl;
+  bool bExcl;
+  t_param param;
+
+  n = mol->atoms.nr;
+  atom = mol->atoms.atom;
+  
+  ntype = sqrt(nbp->nr);
+
+  for (i=0; i<MAXATOMLIST; i++) 
+    param.a[i] = NOTSET;
+  for (i=0; i<MAXFORCEPARAM; i++) 
+    param.c[i] = NOTSET;
+
+  /* Add a pair interaction for all non-excluded atom pairs */
+  excl = &mol->excls;
+  for(i=0; i<n; i++) {
+    for(j=i+1; j<n; j++) {
+      bExcl = FALSE;
+      for(k=excl->index[i]; k<excl->index[i+1]; k++) {
+       if (excl->a[k] == j)
+         bExcl = TRUE;
+      }
+      if (!bExcl) {
+       if (nb_funct != F_LJ)
+         gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
+       param.a[0] = i;
+       param.a[1] = j;
+       param.c[0] = atom[i].q;
+       param.c[1] = atom[j].q;
+       param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
+       param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
+       push_bondnow(&mol->plist[F_LJC_PAIRS_NB],&param);
+      }
+    }
+  }
+}
+
+static void set_excl_all(t_blocka *excl)
+{
+  int nat,i,j,k;
+
+  /* Get rid of the current exclusions and exclude all atom pairs */
+  nat = excl->nr;
+  excl->nra = nat*nat;
+  srenew(excl->a,excl->nra);
+  k = 0;
+  for(i=0; i<nat; i++) {
+    excl->index[i] = k;
+    for(j=0; j<nat; j++)
+      excl->a[k++] = j;
+  }
+  excl->index[nat] = k;
+}
+
+static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
+                          int couple_lam0,int couple_lam1)
+{
+  int i;
+
+  for(i=0; i<atoms->nr; i++) {
+    if (couple_lam0 != ecouplamVDWQ)
+      atoms->atom[i].q     = 0.0;
+    if (couple_lam0 == ecouplamNONE)
+      atoms->atom[i].type  = atomtype_decouple;
+    if (couple_lam1 != ecouplamVDWQ)
+      atoms->atom[i].qB    = 0.0;
+    if (couple_lam1 == ecouplamNONE)
+      atoms->atom[i].typeB = atomtype_decouple;
+  }
+}
+
+void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
+                           int couple_lam0,int couple_lam1,
+                           bool bCoupleIntra,int nb_funct,t_params *nbp)
+{
+  convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
+  if (!bCoupleIntra) {
+    generate_LJCpairsNB(mol,nb_funct,nbp);
+    set_excl_all(&mol->excls);
+  }
+  decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);
+}

diff --git a/src/kernel/toppush.h b/src/kernel/toppush.h
index e9643a14be67ece05a7a705d763ce3473695f6cd..e57521b7031d706cebd4dc04ffd9ce0831cba215 100644 (file)
--- a/src/kernel/toppush.h
+++ b/src/kernel/toppush.h
@@ -77,7 +77,7 @@ extern void push_bondnow (t_params *bond, t_param *b);
 
 extern void push_bond(directive d,t_params bondtype[],t_params bond[],
                      t_atoms *at,t_atomtype atype,char *line,
-                     bool bBonded,bool bGenPairs,
+                     bool bBonded,bool bGenPairs,real fudgeQQ,
                      bool bZero,bool *bWarn_copy_A_B);
 
 extern void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
@@ -100,4 +100,19 @@ extern void b_to_b2(t_blocka *b, t_block2 *b2);
 
 extern void b2_to_b(t_block2 *b2, t_blocka *b);
 
+extern int add_atomtype_decoupled(t_symtab *symtab,t_atomtype at,
+                                 t_nbparam ***nbparam,t_nbparam ***pair);
+/* Add an atom type with all parameters set to zero (no interactions).
+ * Returns the atom type number.
+ */
+
+extern void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,
+                                  real fudgeQQ,
+                                  int couple_lam0,int couple_lam1,
+                                  bool bCoupleIntra,
+                                  int nb_funct,t_params *nbp);
+/* Setup mol such that the B-state has no interaction with the rest
+ * of the system, but full interaction with itself.
+ */
+
 #endif /* _toppush_h */

diff --git a/src/kernel/topshake.c b/src/kernel/topshake.c
index d9d7e1a6445e1e0e7f83ed54048ef9b2f7ff5b43..b91a683c0d0b0f7ac3fcb57ae9fa969fc28b3ca5 100644 (file)
--- a/src/kernel/topshake.c
+++ b/src/kernel/topshake.c
@@ -86,7 +86,6 @@ void make_shake (t_params plist[],t_atoms *atoms,t_atomtype at,int nshake)
   int          nb,b,i,j,ftype,ftype_a;
   bool         bFound;
   
-  p.type = DEF_PARAM_TYPE;
   if (nshake != eshNONE) {
     switch (nshake) {
     case eshHBONDS:

diff --git a/src/kernel/tpbcmp.c b/src/kernel/tpbcmp.c
index 80d492cb11f5fdb38a305623552be40589aaab09..7d738fb93b3e7469f06165395469182135d53cac 100644 (file)
--- a/src/kernel/tpbcmp.c
+++ b/src/kernel/tpbcmp.c
@@ -220,6 +220,7 @@ static void cmp_idef(FILE *fp,t_idef *id1,t_idef *id2,real ftol)
       cmp_iparm(fp,buf2,id1->functype[i],
                id1->iparams[i],id2->iparams[i],ftol);
     }
+    cmp_real(fp,"fudgeQQ",-1,id1->fudgeQQ,id2->fudgeQQ,ftol);
     for(i=0; (i<F_NRE); i++)
       cmp_ilist(fp,i,&(id1->il[i]),&(id2->il[i]));
   } else {
@@ -450,7 +451,6 @@ static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2,real ftol)
   cmp_int(fp,"inputrec->implicit_solvent",-1,ir1->implicit_solvent,ir2->implicit_solvent);
   cmp_int(fp,"inputrec->eDispCorr",-1,ir1->eDispCorr,ir2->eDispCorr);
   cmp_real(fp,"inputrec->shake_tol",-1,ir1->shake_tol,ir2->shake_tol,ftol);
-  cmp_real(fp,"inputrec->fudgeQQ",-1,ir1->fudgeQQ,ir2->fudgeQQ,ftol);
   cmp_int(fp,"inputrec->efep",-1,ir1->efep,ir2->efep);
   cmp_real(fp,"inputrec->init_lambda",-1,ir1->init_lambda,ir2->init_lambda,ftol);
   cmp_real(fp,"inputrec->delta_lambda",-1,ir1->delta_lambda,ir2->delta_lambda,ftol);

diff --git a/src/mdlib/force.c b/src/mdlib/force.c
index dd85e5b4e9e72a3296eb7dd00d89f798d70f72a0..074971699c77dd36f915d59ce770276e08446b34 100644 (file)
--- a/src/mdlib/force.c
+++ b/src/mdlib/force.c
@@ -1030,7 +1030,7 @@ void init_forcerec(FILE *fp,
   /* Electrostatics */
   fr->epsilon_r  = ir->epsilon_r;
   fr->epsilon_rf = ir->epsilon_rf;
-  fr->fudgeQQ    = ir->fudgeQQ;
+  fr->fudgeQQ    = idef->fudgeQQ;
   fr->rcoulomb_switch = ir->rcoulomb_switch;
   fr->rcoulomb        = ir->rcoulomb;
 
@@ -1204,8 +1204,8 @@ void init_forcerec(FILE *fp,
 
   bTab = fr->bcoultab || fr->bvdwtab;
   bSep14tab = ((top->idef.il[F_LJ14].nr > 0 ||
-               top->idef.il[F_LJC14_A].nr > 0 ||
-               top->idef.il[F_LJC_PAIRS_A].nr > 0) &&
+               top->idef.il[F_LJC14_Q].nr > 0 ||
+               top->idef.il[F_LJC_PAIRS_NB].nr > 0) &&
               (!bTab || fr->eeltype!=eelCUT || fr->vdwtype!=evdwCUT));
   
   negp_pp = ir->opts.ngener - ir->nwall;

diff --git a/src/mdlib/mdebin.c b/src/mdlib/mdebin.c
index 2d6830927cc8a2bf204067a0ad283b0b677ce121..0773430f933205fd36eafecd1881e53facc175d0 100644 (file)
--- a/src/mdlib/mdebin.c
+++ b/src/mdlib/mdebin.c
@@ -81,7 +81,7 @@ static char *pcouplmu_nm[] = {
 #define NBOXS asize(boxs_nm)
 #define NTRICLBOXS asize(tricl_boxs_nm)
 
-static bool bConstr,bTricl,bDynBox;
+static bool bConstr,bConstrVir,bTricl,bDynBox;
 static int  f_nre=0,epc,etc,nCrmsd;
 
 t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
@@ -129,7 +129,25 @@ t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
   bool     bBHAM,b14;
   
   bBHAM = (idef->functype[0] == F_BHAM);
-  b14   = (idef->il[F_LJ14].nr > 0 || idef->il[F_LJC14_A].nr > 0);
+  b14   = (idef->il[F_LJ14].nr > 0 || idef->il[F_LJC14_Q].nr > 0);
+
+  nc[0] = idef->il[F_CONSTR].nr/3;
+  nc[1] = idef->il[F_SETTLE].nr*3/2;
+  if (PARTDECOMP(cr))
+    gmx_sumi(2,nc,cr);
+  bConstr    = (nc[0]+nc[1] > 0);
+  bConstrVir = FALSE;
+  if (bConstr) {
+    if (nc[0] > 0 && ir->eConstrAlg == estLINCS) {
+      if (ir->eI == eiSD2)
+       nCrmsd = 2;
+      else
+       nCrmsd = 1;
+    }
+    bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
+  } else {
+    nCrmsd = 0;
+  }
 
   for(i=0; i<F_NRE; i++) {
     bEner[i] = FALSE;
@@ -153,8 +171,12 @@ t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
       bEner[i] = b14;
     else if (i == F_COUL14)
       bEner[i] = b14;
-    else if ((i == F_DVDL) || (i == F_DVDLKIN))
+    else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
+      bEner[i] = FALSE;
+    else if ((i == F_DVDL) || (i == F_DKDL))
       bEner[i] = (ir->efep != efepNO);
+    else if (i == F_DGDL_CON)
+      bEner[i] = (ir->efep != efepNO && bConstr);
     else if ((interaction_function[i].flags & IF_VSITE) ||
             (i == F_CONSTR) || (i == F_SETTLE))
       bEner[i] = FALSE;
@@ -167,9 +189,7 @@ t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
       bEner[i] = (idef->il[F_DISRES].nr > 0);
     else if (i == F_ORIRESDEV)
       bEner[i] = (idef->il[F_ORIRES].nr > 0);
-    else if (i == F_CONNBONDS || 
-            i == F_LJC14_A ||
-            i == F_LJC_PAIRS_A)
+    else if (i == F_CONNBONDS)
       bEner[i] = FALSE;
     else if (i == F_COM_PULL)
       bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F);
@@ -187,22 +207,6 @@ t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
       f_nre++;
     }
 
-  nc[0] = idef->il[F_CONSTR].nr/3;
-  nc[1] = idef->il[F_SETTLE].nr*3/2;
-  if (PARTDECOMP(cr))
-    gmx_sumi(2,nc,cr);
-  bConstr = (nc[0]+nc[1] > 0);
-  if (bConstr) {
-    if (nc[0] > 0 && ir->eConstrAlg == estLINCS) {
-      if (ir->eI == eiSD2)
-       nCrmsd = 2;
-      else
-       nCrmsd = 1;
-    }
-    bConstr = (getenv("GMX_CONSTRAINTVIR") != NULL);
-  } else {
-    nCrmsd = 0;
-  }
   epc = ir->epc;
   bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
   bDynBox = DYNAMIC_BOX(*ir);
@@ -221,7 +225,7 @@ t_mdebin *init_mdebin(int fp_ene,const t_groups *grps,const t_atoms *atoms,
   if (bDynBox)
     md->ib    = get_ebin_space(md->ebin, bTricl ? NTRICLBOXS :
                               NBOXS, bTricl ? tricl_boxs_nm : boxs_nm);
-  if (bConstr) {
+  if (bConstrVir) {
     md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm);
     md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm);
   }
@@ -418,7 +422,7 @@ void upd_mdebin(t_mdebin *md,FILE *fp_dgdl,
       add_ebin(md->ebin,md->ib,NBOXS,bs,bSum,step);
     }
   }  
-  if (bConstr) {
+  if (bConstrVir) {
     add_ebin(md->ebin,md->isvir,9,svir[0],bSum,step);
     add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum,step);
   }
@@ -493,7 +497,8 @@ void upd_mdebin(t_mdebin *md,FILE *fp_dgdl,
     add_ebin(md->ebin,md->iu,3*md->nU,uuu[0],bSum,step);
   }
   if (fp_dgdl)
-    fprintf(fp_dgdl,"%.4f %g\n",time,ener[F_DVDL]+ener[F_DVDLKIN]);
+    fprintf(fp_dgdl,"%.4f %g\n",
+           time,ener[F_DVDL]+ener[F_DKDL]+ener[F_DGDL_CON]);
 }
 
 static void npr(FILE *log,int n,char c)
@@ -595,13 +600,13 @@ void print_ebin(int fp_ene,bool bEne,bool bDR,bool bOR,bool bDihR,
                nsteps,TRUE);      
        fprintf(log,"\n");
       }
-      if (bConstr) {
+      if (bConstrVir) {
        fprintf(log,"   Constraint Virial %s\n",kjm);
        pr_ebin(log,md->ebin,md->isvir,9,3,mode,nsteps,FALSE);  
        fprintf(log,"\n");
        fprintf(log,"   Force Virial %s\n",kjm);
        pr_ebin(log,md->ebin,md->ifvir,9,3,mode,nsteps,FALSE);  
-      fprintf(log,"\n");
+       fprintf(log,"\n");
       }
       fprintf(log,"   Total Virial %s\n",kjm);
       pr_ebin(log,md->ebin,md->ivir,9,3,mode,nsteps,FALSE);   

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 722c26c817a3a7c076f2779e0be559e587975021..d4c2fd422612254edd573c74c565f7bf3569e106 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -173,8 +173,9 @@ static void reset_energies(t_grpopts *opts,t_groups *grp,
   /* Normal potential energy components */
   for(i=0; (i<=F_EPOT); i++)
     epot[i] = 0.0;
-  epot[F_DVDL]    = 0.0;
-  epot[F_DVDLKIN] = 0.0;
+  epot[F_DVDL]     = 0.0;
+  epot[F_DKDL]     = 0.0;
+  epot[F_DGDL_CON] = 0.0;
 }
 
 /*