Merge branch 'release-4-5-patches' into release-4-6
authorSzilard Pall <pszilard@cbr.su.se>
Mon, 3 Sep 2012 10:40:56 +0000 (12:40 +0200)
committerSzilard Pall <pszilard@cbr.su.se>
Mon, 3 Sep 2012 10:40:56 +0000 (12:40 +0200)
Conflicts:
src/mdlib/tpi.c

Change-Id: I1b5d92291ece7cf3194b4e1d88cd686a5aa18696

14 files changed:
include/edsam.h
share/top/charmm27.ff/aminoacids.rtp
share/top/oplsaa.ff/aminoacids.n.tdb
src/kernel/gmxcheck.c
src/kernel/vsite_parm.c
src/mdlib/domdec.c
src/mdlib/edsam.c
src/mdlib/groupcoord.c
src/mdlib/pme.c
src/mdlib/sim_util.c
src/mdlib/tpi.c
src/tools/gmx_energy.c
src/tools/gmx_trjconv.c
src/tools/make_edi.c

index d33c52ad2446d2fb07e222055f75f74d9cb059c1..564d61128b09d99daa6df36b9b0b0eb33bb2618b 100644 (file)
@@ -59,7 +59,7 @@ void dd_make_local_ed_indices(gmx_domdec_t *dd, gmx_edsam_t ed);
  * Should be called at every domain decomposition. */
  
 void do_flood(FILE *log, t_commrec *cr, rvec x[],rvec force[], gmx_edsam_t ed,
-        matrix box, gmx_large_int_t step);
+        matrix box, gmx_large_int_t step, gmx_bool bNS);
 /* Flooding - called from do_force() */
 
 #ifdef __cplusplus
index 7eef837af42ea2867b6f1623e35cf152378fb15b..d24ef05e5173956f03aa56c6cba390f58c9a9345 100644 (file)
        OH2     H1
        OH2     H2
 
+[ HO4 ]
+; TIP4P
+ [ atoms ]
+       OW      OWT4    0.00    0
+       HW1     HWT4    0.52    0
+       HW2     HWT4    0.52    0
+       MW      MWT4    -1.04   0
+ [ bonds ]
+       OW      HW1
+       OW      HW2
+
 [ SOD ]
  [ atoms ]
        SOD     SOD     1.00    0
index 92a6153bb0fe084e37d1703f15383015c16fe210..a790503475dd83b47d3781830f59f07615ece033 100644 (file)
@@ -87,3 +87,6 @@ CA            opls_912B       12.011  0.12
 [ delete ]
 H
 
+[ PGLU-NH ]
+[ replace ]
+CA              opls_246        12.011  0.14
index 84e3880b0613775c6504cef18aa9c8ddbbd0e31d..0c4d26054d0587b10fc19608d4fd0d885cb47d9a 100644 (file)
@@ -157,10 +157,12 @@ static void chk_coords(int frame,int natoms,rvec *x,matrix box,real fac,real tol
        printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
               frame,i,x[i][j]);
     }
-    if ((fabs(x[j][XX]) < tol) && 
-       (fabs(x[j][YY]) < tol) && 
-       (fabs(x[j][ZZ]) < tol))
-      nNul++;
+    if ((fabs(x[i][XX]) < tol) && 
+        (fabs(x[i][YY]) < tol) && 
+        (fabs(x[i][ZZ]) < tol))
+    {
+        nNul++;
+    }
   }
   if (nNul > 0)
     printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
index d0688eff9012c264760e13454c2f910aefc7f771..bcceb218a52cd4dac3ed6c640f62cb80558aa77b 100644 (file)
@@ -38,6 +38,7 @@
 
 #include <stdio.h>
 #include <math.h>
+#include <assert.h>
 #include <string.h>
 #include "vsite_parm.h"
 #include "smalloc.h"
@@ -1020,13 +1021,11 @@ static void clean_vsite_bonds(t_params *plist, t_pindex pindex[],
     if ( bKeep ) {
       if(debug)fprintf(debug," keeping");
       /* now copy the bond to the new array */
-      memcpy(&(ps->param[kept_i]),
-            &(ps->param[i]),(size_t)sizeof(ps->param[0]));
+      ps->param[kept_i] = ps->param[i];
       kept_i++;
     } else if (IS_CHEMBOND(cftype)) {
       srenew(plist[F_CONNBONDS].param,plist[F_CONNBONDS].nr+1);
-      memcpy(&(plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr]),
-            &(ps->param[i]),(size_t)sizeof(plist[F_CONNBONDS].param[0]));
+      plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr] = ps->param[i];
       plist[F_CONNBONDS].nr++;
       nconverted++;
     } else
@@ -1134,10 +1133,9 @@ static void clean_vsite_angles(t_params *plist, t_pindex pindex[],
       }
     
     if ( bKeep ) {
-      /* now copy the angle to the new array */
-      memcpy(&(ps->param[kept_i]),
-            &(ps->param[i]),(size_t)sizeof(ps->param[0]));
-      kept_i++;
+        /* now copy the angle to the new array */
+        ps->param[kept_i] = ps->param[i];
+        kept_i++;
     }
   }
   
@@ -1152,7 +1150,7 @@ static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
 {
   int      ftype,i,parnr,k,l,m,n,nvsite,kept_i,vsnral;
   atom_id  atom,constr;
-  atom_id  vsiteatoms[3];
+  atom_id  vsiteatoms[4];
   gmx_bool     bKeep,bUsed,bPresent;
   t_params *ps;
   
@@ -1171,6 +1169,7 @@ static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
        if (nvsite==1) {
          /* store construction atoms of first vsite */
          vsnral=NRAL(pindex[atom].ftype)-1;
+         assert(vsnral<=4);
          for(m=0; (m<vsnral); m++)
            vsiteatoms[m]=
              plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
@@ -1218,9 +1217,8 @@ static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
     }
       
     if ( bKeep ) {
-      memcpy(&(ps->param[kept_i]),
-            &(ps->param[i]),(size_t)sizeof(ps->param[0]));
-      kept_i++;
+        ps->param[kept_i] = ps->param[i];
+        kept_i++;
     }
   }
 
index 5bfe1f2fa772e8b64dc960e1c184bf1e6b91bab1..472778af04349f2ac47f6e16b70e96bbbc5df8b5 100644 (file)
@@ -4092,7 +4092,7 @@ static void print_cg_move(FILE *fplog,
     fprintf(fplog,"\nStep %s:\n",gmx_step_str(step,buf));
     if (bHaveLimitdAndCMOld)
     {
-        fprintf(fplog,"The charge group starting at atom %d moved than the distance allowed by the domain decomposition (%f) in direction %c\n",
+        fprintf(fplog,"The charge group starting at atom %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
                 ddglatnr(dd,dd->cgindex[cg]),limitd,dim2char(dim));
     }
     else
@@ -4578,19 +4578,19 @@ static int dd_redistribute_cg(FILE *fplog,gmx_large_int_t step,
             if (dim >= npbcdim && dd->nc[dim] > 2)
             {
                 /* No pbc in this dim and more than one domain boundary.
-                 * We to a separate check if a charge did not move too far.
+                 * We do a separate check if a charge group didn't move too far.
                  */
                 if (((flag & DD_FLAG_FW(d)) &&
-                     comm->vbuf.v[buf_pos][d] > cell_x1[dim]) ||
+                     comm->vbuf.v[buf_pos][dim] > cell_x1[dim]) ||
                     ((flag & DD_FLAG_BW(d)) &&
-                     comm->vbuf.v[buf_pos][d] < cell_x0[dim]))
+                     comm->vbuf.v[buf_pos][dim] < cell_x0[dim]))
                 {
-                    cg_move_error(fplog,dd,step,cg,d,
+                    cg_move_error(fplog,dd,step,cg,dim,
                                   (flag & DD_FLAG_FW(d)) ? 1 : 0,
                                    FALSE,0,
                                    comm->vbuf.v[buf_pos],
                                    comm->vbuf.v[buf_pos],
-                                   comm->vbuf.v[buf_pos][d]);
+                                   comm->vbuf.v[buf_pos][dim]);
                 }
             }
 
index 3a0b21dc91703b83aeb11cf58af22e4d6c7f2f09..15fdd6873b416f23058f858bd719c56cc4e94f4d 100644 (file)
@@ -857,7 +857,8 @@ static void do_single_flood(
         t_edpar *edi,
         gmx_large_int_t step,
         matrix box,
-        t_commrec *cr)
+        t_commrec *cr,
+        gmx_bool bNS)       /* Are we in a neighbor searching step? */
 {
     int i;
     matrix  rotmat;         /* rotation matrix */
@@ -868,15 +869,16 @@ static void do_single_flood(
 
     buf=edi->buf->do_edsam;
 
+
     /* Broadcast the positions of the AVERAGE structure such that they are known on
      * every processor. Each node contributes its local positions x and stores them in
      * the collective ED array buf->xcoll */
-    communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, x,
+    communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
 
     /* Only assembly REFERENCE positions if their indices differ from the average ones */
     if (!edi->bRefEqAv)
-        communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, x,
+        communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
                 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
 
     /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
@@ -935,7 +937,8 @@ extern void do_flood(
         rvec            force[], /* forcefield forces, to these the flooding forces are added */
         gmx_edsam_t     ed,      /* ed data structure contains all ED and flooding datasets */
         matrix          box,     /* the box */
-        gmx_large_int_t step)    /* The relative time step since ir->init_step is already subtracted */
+        gmx_large_int_t step,    /* The relative time step since ir->init_step is already subtracted */
+        gmx_bool        bNS)     /* Are we in a neighbor searching step? */
 {
     t_edpar *edi;
 
@@ -948,7 +951,7 @@ extern void do_flood(
     {
         /* Call flooding for one matrix */
         if (edi->flood.vecs.neig)
-            do_single_flood(ed->edo,x,force,edi,step,box,cr);
+            do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
         edi = edi->next_edi;
     }
 }
@@ -1532,8 +1535,7 @@ static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int ed
     /* allocate space for reference positions and read them */
     snew(edi->sref.anrs,edi->sref.nr);
     snew(edi->sref.x   ,edi->sref.nr);
-    if (PAR(cr))
-        snew(edi->sref.x_old,edi->sref.nr);
+    snew(edi->sref.x_old,edi->sref.nr);
     edi->sref.sqrtm    =NULL;
     read_edx(in,edi->sref.nr,edi->sref.anrs,edi->sref.x);
 
@@ -1541,8 +1543,7 @@ static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int ed
     edi->sav.nr=read_checked_edint(in,"NAV");
     snew(edi->sav.anrs,edi->sav.nr);
     snew(edi->sav.x   ,edi->sav.nr);
-    if (PAR(cr))
-        snew(edi->sav.x_old,edi->sav.nr);
+    snew(edi->sav.x_old,edi->sav.nr);
     read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
 
     /* Check if the same atom indices are used for reference and average positions */
@@ -2118,13 +2119,16 @@ static int ed_constraints(gmx_bool edtype, t_edpar *edi)
  * umbrella sampling simulations. */
 static void copyEvecReference(t_eigvec* floodvecs)
 {
-       int i;
+    int i;
+
 
+    if (NULL==floodvecs->refproj0)
+        snew(floodvecs->refproj0, floodvecs->neig);
 
-       for (i=0; i<floodvecs->neig; i++)
-       {
-               floodvecs->refproj0[i] = floodvecs->refproj[i];
-       }
+    for (i=0; i<floodvecs->neig; i++)
+    {
+        floodvecs->refproj0[i] = floodvecs->refproj[i];
+    }
 }
 
 
@@ -2137,7 +2141,7 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
 {
     t_edpar *edi = NULL;    /* points to a single edi data set */
     int     numedis=0;      /* keep track of the number of ED data sets in edi file */
-    int     i,nr_edi;
+    int     i,nr_edi,avindex;
     rvec    *x_pbc  = NULL; /* positions of the whole MD system with pbc removed  */
     rvec    *xfit   = NULL; /* the positions which will be fitted to the reference structure  */
     rvec    *xstart = NULL; /* the positions which are subject to ED sampling */
@@ -2209,10 +2213,9 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             {
                 copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
 
-                /* Save the sref positions such that in the next time step the molecule can
-                 * be made whole again (in the parallel case) */
-                if (PAR(cr))
-                    copy_rvec(xfit[i], edi->sref.x_old[i]);
+                /* Save the sref positions such that in the next time step we can make the ED group whole
+                 * in case any of the atoms do not have the correct PBC representation */
+                copy_rvec(xfit[i], edi->sref.x_old[i]);
             }
 
             /* Extract the positions of the atoms subject to ED sampling */
@@ -2220,10 +2223,9 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             {
                 copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
 
-                /* Save the sav positions such that in the next time step the molecule can
-                 * be made whole again (in the parallel case) */
-                if (PAR(cr))
-                    copy_rvec(xstart[i], edi->sav.x_old[i]);
+                /* Save the sav positions such that in the next time step we can make the ED group whole
+                 * in case any of the atoms do not have the correct PBC representation */
+                copy_rvec(xstart[i], edi->sav.x_old[i]);
             }
 
             /* Make the fit to the REFERENCE structure, get translation and rotation */
@@ -2240,35 +2242,66 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             /* calculate initial projections */
             project(xstart, edi);
 
+            /* For the target and origin structure both a reference (fit) and an
+             * average structure can be provided in make_edi. If both structures
+             * are the same, make_edi only stores one of them in the .edi file.
+             * If they differ, first the fit and then the average structure is stored
+             * in star (or sor), thus the number of entries in star/sor is
+             * (n_fit + n_av) with n_fit the size of the fitting group and n_av
+             * the size of the average group. */
+
             /* process target structure, if required */
             if (edi->star.nr > 0)
             {
                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
+
                 /* get translation & rotation for fit of target structure to reference structure */
                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
                 /* do the fit */
-                translate_and_rotate(edi->star.x, edi->sav.nr, fit_transvec, fit_rotmat);
-                rad_project(edi, edi->star.x, &edi->vecs.radcon, cr);
+                translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
+                if (edi->star.nr == edi->sav.nr)
+                {
+                    avindex = 0;
+                }
+                else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
+                {
+                    /* The last sav.nr indices of the target structure correspond to
+                     * the average structure, which must be projected */
+                    avindex = edi->star.nr - edi->sav.nr;
+                }
+                rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon, cr);
             } else
                 rad_project(edi, xstart, &edi->vecs.radcon, cr);
 
             /* process structure that will serve as origin of expansion circle */
             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
+
             if (edi->sori.nr > 0)
             {
                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
+
                 /* fit this structure to reference structure */
                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
                 /* do the fit */
-                translate_and_rotate(edi->sori.x, edi->sav.nr, fit_transvec, fit_rotmat);
-                rad_project(edi, edi->sori.x, &edi->vecs.radacc, cr);
-                rad_project(edi, edi->sori.x, &edi->vecs.radfix, cr);
+                translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
+                if (edi->sori.nr == edi->sav.nr)
+                {
+                    avindex = 0;
+                }
+                else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
+                {
+                    /* For the projection, we need the last sav.nr indices of sori */
+                    avindex = edi->sori.nr - edi->sav.nr;
+                }
+
+                rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc, cr);
+                rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix, cr);
                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
                 {
                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
                     /* Set center of flooding potential to the ORIGIN structure */
-                    rad_project(edi, edi->sori.x, &edi->flood.vecs, cr);
+                    rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs, cr);
                     /* We already know that no (moving) reference position was provided,
                      * therefore we can overwrite refproj[0]*/
                     copyEvecReference(&edi->flood.vecs);
@@ -2470,7 +2503,7 @@ void do_edsam(t_inputrec  *ir,
              * the collective buf->xcoll array. Note that for edinr > 1
              * xs could already have been modified by an earlier ED */
 
-            communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, xs,
+            communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
 
 #ifdef DEBUG_ED
@@ -2478,11 +2511,12 @@ void do_edsam(t_inputrec  *ir,
 #endif
             /* Only assembly reference positions if their indices differ from the average ones */
             if (!edi->bRefEqAv)
-                communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, xs,
+                communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
                         edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
 
-            /* If bUpdateShifts was TRUE then the shifts have just been updated in get_positions.
-             * We do not need to uptdate the shifts until the next NS step */
+            /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
+             * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
+             * set bUpdateShifts=TRUE in the parallel case. */
             buf->bUpdateShifts = FALSE;
 
             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
index 6ec672badb0e2c175ed28e9336fb1b9ca746d263..2f46d5c60fbcd609ea19daba028dddf0e21b8b51 100644 (file)
@@ -195,7 +195,7 @@ extern void communicate_group_positions(
         rvec       *xcoll,        /* OUT: Collective array of positions */
         ivec       *shifts,       /* IN+OUT: Collective array of shifts for xcoll */
         ivec       *extra_shifts, /* BUF: Extra shifts since last time step */
-        const gmx_bool bNS,           /* IN:  NS step, the shifts have changed */
+        const gmx_bool bNS,       /* IN:  NS step, the shifts have changed */
         rvec       *x_loc,        /* IN:  Local positions on this node */ 
         const int  nr,            /* IN:  Total number of atoms in the group */
         const int  nr_loc,        /* IN:  Local number of atoms in the group */
@@ -221,35 +221,34 @@ extern void communicate_group_positions(
     {
         /* Add the arrays from all nodes together */
         gmx_sum(nr*3, xcoll[0], cr);
+    }
+    /* To make the group whole, start with a whole group and each
+     * step move the assembled positions at closest distance to the positions
+     * from the last step. First shift the positions with the saved shift
+     * vectors (these are 0 when this routine is called for the first time!) */
+    shift_positions_group(box, xcoll, shifts, nr);
+
+    /* Now check if some shifts changed since the last step.
+     * This only needs to be done when the shifts are expected to have changed,
+     * i.e. after neighboursearching */
+    if (bNS)
+    {
+        get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
 
-        /* To make the group whole, start with a whole group and each
-         * step move the assembled positions at closest distance to the positions 
-         * from the last step. First shift the positions with the saved shift 
-         * vectors (these are 0 when this routine is called for the first time!) */
-        shift_positions_group(box, xcoll, shifts, nr);
-        
-        /* Now check if some shifts changed since the last step.
-         * This only needs to be done when the shifts are expected to have changed,
-         * i.e. after neighboursearching */
-        if (bNS) 
+        /* Shift with the additional shifts such that we get a whole group now */
+        shift_positions_group(box, xcoll, extra_shifts, nr);
+
+        /* Add the shift vectors together for the next time step */
+        for (i=0; i<nr; i++)
         {
-            get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
-            
-            /* Shift with the additional shifts such that we get a whole group now */
-            shift_positions_group(box, xcoll, extra_shifts, nr);
-            
-            /* Add the shift vectors together for the next time step */
-            for (i=0; i<nr; i++)
-            {
-                shifts[i][XX] += extra_shifts[i][XX];
-                shifts[i][YY] += extra_shifts[i][YY];
-                shifts[i][ZZ] += extra_shifts[i][ZZ];
-            }
-            
-            /* Store current correctly-shifted positions for comparison in the next NS time step */
-            for (i=0; i<nr; i++)
-                copy_rvec(xcoll[i],xcoll_old[i]);   
+            shifts[i][XX] += extra_shifts[i][XX];
+            shifts[i][YY] += extra_shifts[i][YY];
+            shifts[i][ZZ] += extra_shifts[i][ZZ];
         }
+
+        /* Store current correctly-shifted positions for comparison in the next NS time step */
+        for (i=0; i<nr; i++)
+            copy_rvec(xcoll[i],xcoll_old[i]);
     }
     
     GMX_MPE_LOG(ev_get_group_x_finish);
index c39c978f605de808143b77c114da2bfe481a333c..8342621cad818ecbe0e5fd396580cdf47e2656fa 100644 (file)
@@ -2974,6 +2974,10 @@ int gmx_pme_init(gmx_pme_t *         pmedata,
 
     if (pme->nnodes == 1)
     {
+#ifdef GMX_MPI
+        pme->mpi_comm_d[0] = MPI_COMM_NULL;
+        pme->mpi_comm_d[1] = MPI_COMM_NULL;
+#endif
         pme->ndecompdim = 0;
         pme->nodeid_major = 0;
         pme->nodeid_minor = 0;
index bf6743544041b800f1be9c05cca044848f8bbfca..e06f3d88ce411a6dcdc26d2978d30fe31491b061 100644 (file)
@@ -796,7 +796,7 @@ void do_force(FILE *fplog,t_commrec *cr,
 
     if (ed)
     {
-        do_flood(fplog,cr,x,f,ed,box,step);
+        do_flood(fplog,cr,x,f,ed,box,step,bNS);
     }
 
     if (DOMAINDECOMP(cr))
index 9e0fe967c6ab7f1d7662c939925ce4b5f9ff1a8b..65e3fde04336ab0b28763c7ea27c64e532cc2571 100644 (file)
@@ -414,8 +414,9 @@ double do_tpi(FILE *fplog,t_commrec *cr,
         {
             copy_rvec(rerun_fr.x[i],state->x[i]);
         }
+        copy_mat(rerun_fr.box,state->box);
         
-        V = det(rerun_fr.box);
+        V = det(state->box);
         logV = log(V);
         
         bStateChanged = TRUE;
@@ -562,20 +563,19 @@ double do_tpi(FILE *fplog,t_commrec *cr,
                 cr->nnodes = 1;
                 do_force(fplog,cr,inputrec,
                          step,nrnb,wcycle,top,top_global,&top_global->groups,
-                         rerun_fr.box,state->x,&state->hist,
+                         state->box,state->x,&state->hist,
                          f,force_vir,mdatoms,enerd,fcd,
                          state->lambda,
                          NULL,fr,NULL,mu_tot,t,NULL,NULL,FALSE,
                          GMX_FORCE_NONBONDED |
-                         (bNS ? GMX_FORCE_NS | GMX_FORCE_DOLR : 0) |
+                         (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS | GMX_FORCE_DOLR : 0) |
                          (bStateChanged ? GMX_FORCE_STATECHANGED : 0)); 
                 cr->nnodes = nnodes;
                 bStateChanged = FALSE;
                 bNS = FALSE;
                 
                 /* Calculate long range corrections to pressure and energy */
-                calc_dispcorr(fplog,inputrec,fr,step,top_global->natoms,
-                              rerun_fr.box,
+                calc_dispcorr(fplog,inputrec,fr,step,top_global->natoms,state->box,
                               lambda,pres,vir,&prescorr,&enercorr,&dvdlcorr);
                 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
                 enerd->term[F_DISPCORR] = enercorr;
index 52ec7ffbd14d7aff53e061a015647f3a2aa4cbb7..83c00bf2a718afb8a1ac562dd21b997f33d42e61 100644 (file)
@@ -1606,7 +1606,8 @@ int gmx_energy(int argc,char *argv[])
     "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
 
     "Some fluctuation-dependent properties can be calculated provided",
-    "the correct energy terms are selected. The following properties",
+    "the correct energy terms are selected, and that the command line option",
+    "[TT]-fluct_props[tt] is given. The following properties",
     "will be computed:[BR]",
     "Property                        Energy terms needed[BR]",
     "---------------------------------------------------[BR]",
@@ -1670,7 +1671,7 @@ int gmx_energy(int argc,char *argv[])
 
   };
   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
-  static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
+  static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
   static int  skip=0,nmol=1,nbmin=5,nbmax=5;
   static real reftemp=300.0,ezero=0;
   t_pargs pa[] = {
@@ -1696,6 +1697,8 @@ int gmx_energy(int argc,char *argv[])
       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
     { "-nmol", FALSE, etINT,  {&nmol},
       "Number of molecules in your sample: the energies are divided by this number" },
+    { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
+      "Compute properties based on energy fluctuations, like heat capacity" },
     { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
       "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
     { "-fluc", FALSE, etBOOL, {&bFluct},
@@ -2420,8 +2423,9 @@ int gmx_energy(int argc,char *argv[])
                    time,reftemp,&edat,
                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
                    oenv);
-      calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
-                             nbmin,nbmax);
+      if (bFluctProps)
+          calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
+                                 nbmin,nbmax);
   }
   if (opt2bSet("-f2",NFILE,fnm)) {
       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
index 330a67766809387e604a768ad9b28968eb2e87a1..b7c323358176282300963d4d5e367bc100bdd22a 100644 (file)
@@ -1361,8 +1361,8 @@ int gmx_trjconv(int argc,char *argv[])
                         }
                         /* Copy the input trxframe struct to the output trxframe struct */
                         frout = fr;
-                       frout.bV    &= bVels;
-                       frout.bF    &= bForce;
+                       frout.bV = (frout.bV && bVels);
+                       frout.bF = (frout.bF && bForce);
                         frout.natoms = nout;
                         if (bNeedPrec && (bSetPrec || !fr.bPrec)) {
                             frout.bPrec = TRUE;
index d715409034828199dcc59a39065510bac422e65d..1a9b9ef622a1c860730b7c092047661aa9cd06e2 100644 (file)
@@ -443,7 +443,7 @@ void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
 
 void get_structure(t_atoms *atoms,const char *IndexFile,
                    const char *StructureFile,struct edix *edx,int nfit,
-                   atom_id ifit[],int natoms, atom_id index[])
+                   atom_id ifit[],int nav, atom_id index[])
 {
   atom_id *igro;  /*index corresponding to target or origin structure*/
   int ngro;
@@ -461,8 +461,10 @@ void get_structure(t_atoms *atoms,const char *IndexFile,
      gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
   init_edx(edx);
   filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
+
+  /* If average and reference/fitting structure differ, append the average structure as well */
   if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
-     filter2edx(edx,natoms,index,ngro,igro,xtar,StructureFile);
+     filter2edx(edx,nav,index,ngro,igro,xtar,StructureFile);
 }
 
 int main(int argc,char *argv[])
@@ -636,12 +638,12 @@ int main(int argc,char *argv[])
     int        nvec1,*eignr1=NULL;
     rvec       *xav1,**eigvec1=NULL;
     t_atoms    *atoms=NULL;
-    int natoms;
+    int        nav;  /* Number of atoms in the average structure */
     char       *grpname;
     const char *indexfile;
     int        i;
     atom_id    *index,*ifit;
-    int        nfit;
+    int        nfit; /* Number of atoms in the reference/fit structure */
     int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
     int nvecs;
     real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */
@@ -730,7 +732,7 @@ int main(int argc,char *argv[])
     EigvecFile=opt2fn("-f",NFILE,fnm);
 
     /*read eigenvectors from eigvec.trr*/
-    read_eigenvectors(EigvecFile,&natoms,&bFit1,
+    read_eigenvectors(EigvecFile,&nav,&bFit1,
                       &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);
 
     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
@@ -738,11 +740,11 @@ int main(int argc,char *argv[])
     atoms=&top.atoms;
 
 
-    printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
+    printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",nav);
     get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
-    if (i!=natoms) {
+    if (i!=nav) {
         gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
-                  i,natoms);
+                  i,nav);
     }
     printf("\n");
 
@@ -768,7 +770,7 @@ int main(int argc,char *argv[])
     }
     else
     {
-        nfit=natoms;
+        nfit=nav;
         ifit=index;
     }
 
@@ -804,15 +806,15 @@ int main(int argc,char *argv[])
         }
     }
 
-    edi_params.ned=natoms;
+    edi_params.ned=nav;
 
   /*number of system atoms  */
   edi_params.nini=atoms->nr;
 
 
   /*store reference and average structure in edi_params*/
-  make_t_edx(&edi_params.sref,nfit,xref1,ifit);
-  make_t_edx(&edi_params.sav,natoms,xav1,index);
+  make_t_edx(&edi_params.sref,nfit,xref1,ifit );
+  make_t_edx(&edi_params.sav ,nav ,xav1 ,index);
 
 
   /* Store target positions in edi_params */
@@ -823,7 +825,7 @@ int main(int argc,char *argv[])
           fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
                           "      You may want to use -ori to define the flooding potential center.\n\n");
       }
-      get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,natoms,index);
+      get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,nav,index);
   }
   else
   {
@@ -833,7 +835,7 @@ int main(int argc,char *argv[])
   /* Store origin positions */
   if (opt2bSet("-ori",NFILE,fnm))
   {
-      get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,natoms,index);
+      get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,nav,index);
   }
   else
   {