changed number of iatoms for F_SETTLE from 1 to 3

Change-Id: Ia0cb1f277b5d3e8988c6fd1fb26896fa68c2fb3c

diff --git a/src/gmxlib/ifunc.c b/src/gmxlib/ifunc.c
index c3b62831cdaba02490fb3105d3bf3c2223c96001..316d01f6a308ce4b75eacfcdaf4b9069e4b88f68 100644 (file)
--- a/src/gmxlib/ifunc.c
+++ b/src/gmxlib/ifunc.c
@@ -147,7 +147,7 @@ const t_interaction_function interaction_function[F_NRE]=
   def_nofc    ("DIHVIOL",  "Dih. Rest. viol."                                     ),    
   def_shkcb   ("CONSTR",   "Constraint",      2, 1, 1                             ),
   def_shk     ("CONSTRNC", "Constr. No Conn.",2, 1, 1                             ),
-  def_shkcb   ("SETTLE",   "Settle",          1, 2, 0                             ),
+  def_shkcb   ("SETTLE",   "Settle",          3, 2, 0                             ),
   def_vsite   ("VSITE2",   "Virtual site 2",  3, 1                                ),
   def_vsite   ("VSITE3",   "Virtual site 3",  4, 2                                ),
   def_vsite   ("VSITE3FD", "Virtual site 3fd",4, 2                                ),

diff --git a/src/gmxlib/mshift.c b/src/gmxlib/mshift.c
index 0e9e102539a6f3da260b9fb48b37ce2d7e4052e9..7fb2314fdb64c2939a82096fc22b29540a0d91ed 100644 (file)
--- a/src/gmxlib/mshift.c
+++ b/src/gmxlib/mshift.c
@@ -98,7 +98,7 @@ static void mk_igraph(t_graph *g,int ftype,t_ilist *il,
     np = interaction_function[ftype].nratoms;
     
     if (ia[1] >= at_start && ia[1] < at_end) {
-      if (ia[np] >= at_end || (ftype == F_SETTLE && ia[1]+2 >= at_end))
+      if (ia[np] >= at_end)
        gmx_fatal(FARGS,
                  "Molecule in topology has atom numbers below and "
                  "above natoms (%d).\n"
@@ -108,8 +108,8 @@ static void mk_igraph(t_graph *g,int ftype,t_ilist *il,
                  at_end,at_end);
       if (ftype == F_SETTLE) {
        /* Bond all the atoms in the settle */
-       add_gbond(g,ia[1],ia[1]+1);
-       add_gbond(g,ia[1],ia[1]+2);
+       add_gbond(g,ia[1],ia[2]);
+       add_gbond(g,ia[1],ia[3]);
       } else if (part == NULL) {
        /* Simply add this bond */
        for(j=1; j<np; j++) {
@@ -179,8 +179,8 @@ static void calc_1se(t_graph *g,int ftype,t_ilist *il,
       iaa          = ia[1];
       if (iaa >= at_start && iaa < at_end) {
        nbond[iaa]   += 2;
-       nbond[iaa+1] += 1;
-       nbond[iaa+2] += 1;
+       nbond[ia[2]] += 1;
+       nbond[ia[3]] += 1;
        g->at_start   = min(g->at_start,iaa);
        g->at_end     = max(g->at_end,iaa+2+1);
       }

diff --git a/src/gmxlib/splitter.c b/src/gmxlib/splitter.c
index 284aa1c3b9f59ff3a0274e38215cef67b85b4a43..ce0be7646de3b56e19d1119b9d1b8faf6c6c6cd9 100644 (file)
--- a/src/gmxlib/splitter.c
+++ b/src/gmxlib/splitter.c
@@ -188,10 +188,11 @@ static void split_force2(t_inputrec *ir, int nnodes,int hid[],int ftype,t_ilist
       /* Only the first particle is stored for settles ... */
       ai=ilist->iatoms[i+1];
       nodeid=hid[ai];
-      if ((nodeid != hid[ai+1]) ||
-         (nodeid != hid[ai+2]))
+      if (nodeid != hid[ilist->iatoms[i+2]] ||
+          nodeid != hid[ilist->iatoms[i+3]])
        gmx_fatal(FARGS,"Settle block crossing node boundaries\n"
-                 "constraint between atoms (%d-%d)",ai+1,ai+2+1);
+                 "constraint between atoms %d, %d, %d)",
+              ai,ilist->iatoms[i+2],ilist->iatoms[i+3]);
     }
     else if(interaction_function[ftype].flags & IF_VSITE) 
        {

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index c640bb5dac989eefa2a9ec9c00a2396ce6aaae40..88fb4e335c57841378008dd84632a676b896bb2e 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -73,7 +73,7 @@
 static const char *tpx_tag = TPX_TAG_RELEASE;
 
 /* This number should be increased whenever the file format changes! */
-static const int tpx_version = 77;
+static const int tpx_version = 78;
 
 /* This number should only be increased when you edit the TOPOLOGY section
  * of the tpx format. This way we can maintain forward compatibility too
@@ -84,7 +84,7 @@ static const int tpx_version = 77;
  * to the end of the tpx file, so we can just skip it if we only
  * want the topology.
  */
-static const int tpx_generation = 23;
+static const int tpx_generation = 24;
 
 /* 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.
@@ -1360,6 +1360,22 @@ static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams,
   }
 }
 
+static void add_settle_atoms(t_ilist *ilist)
+{
+    int i;
+
+    /* Settle used to only store the first atom: add the other two */
+    srenew(ilist->iatoms,2*ilist->nr);
+    for(i=ilist->nr/2-1; i>=0; i--)
+    {
+        ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
+        ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
+        ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
+        ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
+    }
+    ilist->nr = 2*ilist->nr;
+}
+
 static void do_ilists(t_fileio *fio, t_ilist *ilist,gmx_bool bRead, 
                       int file_version)
 {
@@ -1378,6 +1394,10 @@ static void do_ilists(t_fileio *fio, t_ilist *ilist,gmx_bool bRead,
       ilist[j].iatoms = NULL;
     } else {
       do_ilist(fio, &ilist[j],bRead,file_version,j);
+      if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
+      {
+          add_settle_atoms(&ilist[j]);
+      }
     }
     /*
     if (bRead && gmx_debug_at)

diff --git a/src/kernel/grompp.c b/src/kernel/grompp.c
index d3b43bbcb268638472fe755a94608cfe043f0aec..bde3490087fa92a622e5c9ef8e5cd27cf8caebad 100644 (file)
--- a/src/kernel/grompp.c
+++ b/src/kernel/grompp.c
@@ -282,10 +282,10 @@ static void check_bonds_timestep(gmx_mtop_t *mtop,double dt,warninp_t wi)
                                 bFound = TRUE;
                             }
                         }
-                    for(j=0; j<ils->nr; j+=2)
+                    for(j=0; j<ils->nr; j+=4)
                     {
-                        if ((a1 >= ils->iatoms[j+1] && a1 < ils->iatoms[j+1]+3) &&
-                            (a2 >= ils->iatoms[j+1] && a2 < ils->iatoms[j+1]+3))
+                        if ((a1 == ils->iatoms[j+1] || a1 == ils->iatoms[j+2] || a1 == ils->iatoms[j+3]) &&
+                            (a2 == ils->iatoms[j+1] || a2 == ils->iatoms[j+2] || a2 == ils->iatoms[j+3]))
                         {
                             bFound = TRUE;
                         }

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index 3b6e4e7f4b68e1829aa959b45dc8125a886ffb65..9c1d3da3d8efb5968430eee0f18886de86848408 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -1601,14 +1601,15 @@ static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
       ia = molt->ilist[F_SETTLE].iatoms;
       for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
        /* Subtract 1 dof from every atom in the SETTLE */
-       for(ai=as+ia[1]; ai<as+ia[1]+3; ai++) {
+       for(j=0; j<3; j++) {
+      ai = as + ia[1+j];
          imin = min(2,nrdf2[ai]);
          nrdf2[ai] -= imin;
          nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
          nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
        }
-       ia += 2;
-       i  += 2;
+       ia += 4;
+       i  += 4;
       }
       as += molt->atoms.nr;
     }

diff --git a/src/kernel/toppush.c b/src/kernel/toppush.c
index eabb245ce43dd30c30a587403745c5be2b94d8f8..1ce51d0500ab8de032377ead1ede47cbf1bd5d13 100644 (file)
--- a/src/kernel/toppush.c
+++ b/src/kernel/toppush.c
@@ -1487,7 +1487,7 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
     "%*s%*s%*s%*s%*s%*s%*s"
   };
   const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
-  int      nr,i,j,nral,nread,ftype;
+  int      nr,i,j,nral,nral_fmt,nread,ftype;
   char     format[STRLEN];
   /* One force parameter more, so we can check if we read too many */
   double   cc[MAXFORCEPARAM+1];
@@ -1499,36 +1499,65 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
 
   nparam_defA=nparam_defB=0;
        
-  ftype = ifunc_index(d,1);
-  nral  = NRAL(ftype);
-  for(j=0; j<MAXATOMLIST; j++)
-    aa[j]=NOTSET;
-  bDef = (NRFP(ftype) > 0);
-  
-  nread = sscanf(line,aaformat[nral-1],
-                &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
-  if (nread < nral) {
-    too_few(wi);
-    return;
-  } else if (nread == nral) 
     ftype = ifunc_index(d,1);
-  else {
-    /* this is a hack to allow for virtual sites with swapped parity */
-    bSwapParity = (aa[nral]<0);
-    if (bSwapParity)
-      aa[nral] = -aa[nral];
-    ftype = ifunc_index(d,aa[nral]);
-    if (bSwapParity)
-      switch(ftype) {
-      case F_VSITE3FAD:
-      case F_VSITE3OUT:
-       break;
-      default:
-       gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
-                   interaction_function[F_VSITE3FAD].longname,
-                   interaction_function[F_VSITE3OUT].longname);
-      }
-  }
+    nral  = NRAL(ftype);
+    for(j=0; j<MAXATOMLIST; j++)
+    {
+        aa[j] = NOTSET;
+    }
+    bDef = (NRFP(ftype) > 0);
+
+    if (ftype == F_SETTLE)
+    {
+        /* SETTLE acts on 3 atoms, but the topology format only specifies
+         * the first atom (for historical reasons).
+         */
+        nral_fmt = 1;
+    }
+    else
+    {
+        nral_fmt = nral;
+    }
+
+    nread = sscanf(line,aaformat[nral_fmt-1],
+                   &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
+
+    if (ftype == F_SETTLE)
+    {
+        aa[3] = aa[1];
+        aa[1] = aa[0] + 1;
+        aa[2] = aa[0] + 2;
+    }
+
+    if (nread < nral_fmt)
+    {
+        too_few(wi);
+        return;
+    }
+    else if (nread > nral_fmt)
+    {
+        /* this is a hack to allow for virtual sites with swapped parity */
+        bSwapParity = (aa[nral]<0);
+        if (bSwapParity)
+        {
+            aa[nral] = -aa[nral];
+        }
+        ftype = ifunc_index(d,aa[nral]);
+        if (bSwapParity)
+        {
+            switch(ftype)
+            {
+            case F_VSITE3FAD:
+            case F_VSITE3OUT:
+                break;
+            default:
+                gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
+                          interaction_function[F_VSITE3FAD].longname,
+                          interaction_function[F_VSITE3OUT].longname);
+            }
+        }
+    }
+
   
   /* Check for double atoms and atoms out of bounds */
   for(i=0; (i<nral); i++) {
@@ -1546,14 +1575,6 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
        warning(wi,errbuf);
       }
   }
-  if (ftype == F_SETTLE)
-    if (aa[0]+2 > at->nr)
-      gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
-                 "             Atom index (%d) in %s out of bounds (1-%d)\n"
-                 "             Settle works on atoms %d, %d and %d",
-                 get_warning_file(wi),get_warning_line(wi),
-                 aa[0],dir2str(d),at->nr,
-                 aa[0],aa[0]+1,aa[0]+2);
   
   /* default force parameters  */
   for(j=0; (j<MAXATOMLIST); j++)
@@ -1597,10 +1618,10 @@ void push_bond(directive d,t_params bondtype[],t_params bond[],
     gmx_incons("Unknown function type in push_bond");
   }
        
-  if (nread > nral) {  
+  if (nread > nral_fmt) {  
          /* Manually specified parameters - in this case we discard multiple torsion info! */
-         
-    strcpy(format,asformat[nral-1]);
+
+    strcpy(format,asformat[nral_fmt-1]);
     strcat(format,ccformat);
     
     nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],

diff --git a/src/mdlib/constr.c b/src/mdlib/constr.c
index 94abc70af7d517fe6330ee33dcf85caa4c5d86a5..6997e13c31e4184645521004b1d4e3cfe2aeaad2 100644 (file)
--- a/src/mdlib/constr.c
+++ b/src/mdlib/constr.c
@@ -406,7 +406,7 @@ gmx_bool constrain(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
     settle  = &idef->il[F_SETTLE];
     if (settle->nr > 0)
     {
-        nsettle = settle->nr/2;
+        nsettle = settle->nr/4;
         
         switch (econq)
         {
@@ -431,7 +431,7 @@ gmx_bool constrain(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
                 sprintf(buf,
                         "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
                         "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
-                        step,ddglatnr(cr->dd,settle->iatoms[error*2+1]));
+                        step,ddglatnr(cr->dd,settle->iatoms[error*4+1]));
                 if (fplog)
                 {
                     fprintf(fplog,"%s",buf);
@@ -780,7 +780,7 @@ void set_constraints(struct gmx_constr *constr,
     {
         settle = &idef->il[F_SETTLE];
         iO = settle->iatoms[1];
-        iH = settle->iatoms[1]+1;
+        iH = settle->iatoms[2];
         constr->settled =
             settle_init(md->massT[iO],md->massT[iH],
                         md->invmass[iO],md->invmass[iH],
@@ -1046,7 +1046,7 @@ gmx_constr_t init_constraints(FILE *fplog,
         iloop = gmx_mtop_ilistloop_init(mtop);
         while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol)) 
         {
-            for (i=0; i<ilist[F_SETTLE].nr; i+=2) 
+            for (i=0; i<ilist[F_SETTLE].nr; i+=4) 
             {
                 if (settle_type == -1) 
                 {

diff --git a/src/mdlib/csettle.c b/src/mdlib/csettle.c
index dbab4f2da1f68623492aa19e5273e473cd5e398c..eea04e2e3c12217e9e4183e401385e4d35f0bfd5 100644 (file)
--- a/src/mdlib/csettle.c
+++ b/src/mdlib/csettle.c
@@ -213,9 +213,9 @@ void settle_proj(FILE *fp,
     
     for (i=0; i<nsettle; i++)
     {
-        ow1 = iatoms[i*2+1];
-        hw2 = ow1 + 1;
-        hw3 = ow1 + 2;
+        ow1 = iatoms[i*4+1];
+        hw2 = iatoms[i*4+2];
+        hw3 = iatoms[i*4+3];
 
 
         for(m=0; m<DIM; m++)
@@ -441,9 +441,9 @@ void csettle(gmx_settledata_t settled,
   for (i = 0; i < nsettle; ++i) {
     doshake = 0;
     /*    --- Step1  A1' ---      */
-    ow1 = iatoms[i*2+1] * 3;
-    hw2 = ow1 + 3;
-    hw3 = ow1 + 6;
+    ow1 = iatoms[i*4+1] * 3;
+    hw2 = iatoms[i*4+2] * 3;
+    hw3 = iatoms[i*4+3] * 3;
     xb0 = b4[hw2    ] - b4[ow1];
     yb0 = b4[hw2 + 1] - b4[ow1 + 1];
     zb0 = b4[hw2 + 2] - b4[ow1 + 2];

diff --git a/src/ngmx/manager.c b/src/ngmx/manager.c
index b18ece2402fcbdbf270bdad12b3c0d675d0ecc40..23f883189f79e38f87f7a52b341b293171884140 100644 (file)
--- a/src/ngmx/manager.c
+++ b/src/ngmx/manager.c
@@ -87,8 +87,8 @@ static void add_bonds(t_manager *man,t_functype func[],
     delta = interaction_function[ftype].nratoms;
     
     if (ftype == F_SETTLE) {
-      aj=ai+1;
-      ak=ai+2;
+      aj = ia[2];
+      ak = ia[3];
       bB[ai]=bB[aj]=bB[ak]=TRUE;
       add_object(man,eOHBond,ai,aj);
       add_object(man,eOHBond,ai,ak);

diff --git a/src/tools/gmx_hbond.c b/src/tools/gmx_hbond.c
index e0c77591a7a49a3ee630dc89a5603766d8f67af3..3193e67588548f1ed55b7e1c0ce8ac66ce492ad5 100644 (file)
--- a/src/tools/gmx_hbond.c
+++ b/src/tools/gmx_hbond.c
@@ -1022,7 +1022,7 @@ static void search_donors(t_topology *top, int isize, atom_id *index,
     int        i,j,nra,n;
     t_functype func_type;
     t_ilist    *interaction;
-    atom_id    nr1,nr2;
+    atom_id    nr1,nr2,nr3;
     gmx_bool       stop;
 
     if (!ddd->dptr) {
@@ -1057,14 +1057,16 @@ static void search_donors(t_topology *top, int isize, atom_id *index,
        
                 /* check out this functype */
                 if (func_type == F_SETTLE) {
-                    nr1=interaction->iatoms[i+1];
+                    nr1 = interaction->iatoms[i+1];
+                    nr2 = interaction->iatoms[i+2];
+                    nr3 = interaction->iatoms[i+3];
          
                     if (ISINGRP(datable[nr1])) {
-                        if (ISINGRP(datable[nr1+1])) {
+                        if (ISINGRP(datable[nr2])) {
                             datable[nr1] |= DON;
                             add_dh(ddd,nr1,nr1+1,grp,datable);
                         }
-                        if (ISINGRP(datable[nr1+2])) {
+                        if (ISINGRP(datable[nr3])) {
                             datable[nr1] |= DON;
                             add_dh(ddd,nr1,nr1+2,grp,datable);
                         }