Added PBC checking patch to g_cluster from Alexey Shvetsov (Bugzilla 526)

diff --git a/src/tools/gmx_cluster.c b/src/tools/gmx_cluster.c
index 3f2fd734f819d059dfdb16c66de1d398c545a88d..c393379d3513d39af44139564242249b74e0673a 100644 (file)
--- a/src/tools/gmx_cluster.c
+++ b/src/tools/gmx_cluster.c
@@ -50,6 +50,7 @@
 #include "index.h"
 #include "random.h"
 #include "pbc.h"
+#include "rmpbc.h"
 #include "xvgr.h"
 #include "futil.h"
 #include "matio.h"
@@ -239,7 +240,7 @@ void gather(t_mat *m,real cutoff,t_clusters *clust)
   t_clustid *c;
   t_dist    *d;
   int       i,j,k,nn,cid,n1,diff;
-  gmx_bool      bChange;
+  gmx_bool  bChange;
   
   /* First we sort the entries in the RMSD matrix */
   n1 = m->nn;
@@ -333,7 +334,7 @@ static void jarvis_patrick(int n1,real **mat,int M,int P,
   t_clustid *c;
   int       **nnb;
   int       i,j,k,cid,diff,max;
-  gmx_bool      bChange;
+  gmx_bool  bChange;
   real      **mcpy=NULL;
 
   if (rmsdcut < 0)
@@ -390,7 +391,7 @@ static void jarvis_patrick(int n1,real **mat,int M,int P,
 
   c = new_clustid(n1);
   fprintf(stderr,"Linking structures ");
-  /* Use mcpy for temporary storage of gmx_booleans */
+  /* Use mcpy for temporary storage of booleans */
   mcpy = mk_matrix(n1,n1,FALSE);
   for(i=0; i<n1; i++)
     for(j=i+1; j<n1; j++)
@@ -549,7 +550,7 @@ static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
 }
 
 rvec **read_whole_trj(const char *fn,int isize,atom_id index[],int skip,
-                      int *nframe, real **time,const output_env_t oenv)
+                      int *nframe, real **time,const output_env_t oenv,gmx_bool bPBC, gmx_rmpbc_t gpbc)
 {
   rvec   **xx,*x;
   matrix box;
@@ -566,6 +567,9 @@ rvec **read_whole_trj(const char *fn,int isize,atom_id index[],int skip,
   i  = 0;
   i0 = 0;
   do {
+    if (bPBC) {
+      gmx_rmpbc(gpbc,natom,box,x);
+    }
     if (i0 >= max_nf) {
       max_nf += 10;
       srenew(xx,max_nf);
@@ -1027,7 +1031,7 @@ int gmx_cluster(int argc,char *argv[])
   char     *grpname;
   real     rmsd,**d1,**d2,*time=NULL,time_invfac,*mass=NULL;
   char     buf[STRLEN],buf1[80],title[STRLEN];
-  gmx_bool     bAnalyze,bUseRmsdCut,bJP_RMSD=FALSE,bReadMat,bReadTraj;
+  gmx_bool bAnalyze,bUseRmsdCut,bJP_RMSD=FALSE,bReadMat,bReadTraj,bPBC=TRUE;
 
   int method,ncluster=0;  
   static const char *methodname[] = { 
@@ -1043,11 +1047,13 @@ int gmx_cluster(int argc,char *argv[])
   static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
   static int  nlevels=40,skip=1;
   static real scalemax=-1.0,rmsdcut=0.1,rmsmin=0.0;
-  static gmx_bool bRMSdist=FALSE,bBinary=FALSE,bAverage=FALSE,bFit=TRUE;
+  gmx_bool bRMSdist=FALSE,bBinary=FALSE,bAverage=FALSE,bFit=TRUE;
   static int  niter=10000,seed=1993,write_ncl=0,write_nst=1,minstruct=1;
   static real kT=1e-3;
   static int  M=10,P=3;
   output_env_t oenv;
+  gmx_rmpbc_t gpbc=NULL;
+  
   t_pargs pa[] = {
     { "-dista", FALSE, etBOOL, {&bRMSdist},
       "Use RMSD of distances instead of RMS deviation" },
@@ -1087,7 +1093,9 @@ int gmx_cluster(int argc,char *argv[])
       "Number of iterations for MC" },
     { "-kT",    FALSE, etREAL, {&kT},
       "Boltzmann weighting factor for Monte Carlo optimization "
-      "(zero turns off uphill steps)" }
+      "(zero turns off uphill steps)" },
+    { "-pbc", FALSE, etBOOL,
+      { &bPBC }, "PBC check" }
   };
   t_filenm fnm[] = {
     { efTRX, "-f",     NULL,        ffOPTRD },
@@ -1182,6 +1190,9 @@ int gmx_cluster(int argc,char *argv[])
     /* don't read mass-database as masses (and top) are not used */
     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),buf,&top,&ePBC,&xtps,NULL,box,
                  bAnalyze);
+    if(bPBC) {
+       gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
+    }
     
     fprintf(stderr,"\nSelect group for least squares fit%s:\n",
            bReadMat?"":" and RMSD calculation");
@@ -1233,8 +1244,8 @@ int gmx_cluster(int argc,char *argv[])
   if (bReadTraj) {
     /* Loop over first coordinate file */
     fn = opt2fn("-f",NFILE,fnm);
-    
-    xx = read_whole_trj(fn,isize,index,skip,&nf,&time,oenv);
+
+    xx = read_whole_trj(fn,isize,index,skip,&nf,&time,oenv,&bPBC,&gpbc);
     output_env_conv_times(oenv, nf, time);
     if (!bRMSdist || bAnalyze) {
       /* Center all frames on zero */
@@ -1246,6 +1257,10 @@ int gmx_cluster(int argc,char *argv[])
        reset_x(ifsize,fitidx,isize,NULL,xx[i],mass);
     }
   }
+  if (bPBC) {
+    gmx_rmpbc_done(gpbc);
+  }
+
   if (bReadMat) {
     fprintf(stderr,"Reading rms distance matrix ");
     read_xpm_matrix(opt2fn("-dm",NFILE,fnm),&readmat);