Merge release-4.5-patches into HEAD

[alexxy/gromacs.git] / src / tools / gmx_hbond.c

diff --git a/src/tools/gmx_hbond.c b/src/tools/gmx_hbond.c
index e7e31384f5798f484ed0ea94a15d32aaf069ef6d..375f6df2df610003d398e0c4cc818cbb01da364c 100644 (file)
--- a/src/tools/gmx_hbond.c
+++ b/src/tools/gmx_hbond.c
@@ -796,12 +796,15 @@ static void add_hbond(t_hbdata *hb,int d,int a,int h,int grpd,int grpa,
                   grpa,hb->a.grp[ia],a+1);
 
     if (bMerge)
-        if ((daSwap = isInterchangable(hb, d, a, grpd, grpa) || bContact) && d>a)
+    {
+        
+        if (isInterchangable(hb, d, a, grpd, grpa) && d>a)
             /* Then swap identity so that the id of d is lower then that of a.
              *
              * This should really be redundant by now, as is_hbond() now ought to return
              * hbNo in the cases where this conditional is TRUE. */
         {
+            daSwap = TRUE;
             k = d;
             d = a;
             a = k;
@@ -818,6 +821,7 @@ static void add_hbond(t_hbdata *hb,int d,int a,int h,int grpd,int grpa,
                 gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
                           grpa,hb->a.grp[ia],a+1);
         }
+    }
 
     if (hb->hbmap) {
         /* Loop over hydrogens to find which hydrogen is in this particular HB */
@@ -1437,12 +1441,13 @@ static int is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
     rvec_sub(x[d],x[a],r_da);
     /* Insert projection code here */
 
-    /* if (d>a && ((isInterchangable(hb, d, a, grpd, grpa) && bMerge) || bContact)) */
-/*         /\* Then this hbond will be found again, or it has already been found. *\/ */
-/*         return hbNo; */
-
+    if (bMerge && d>a && isInterchangable(hb, d, a, grpd, grpa))
+    {
+        /* Then this hbond/contact will be found again, or it has already been found. */
+        /*return hbNo;*/
+    }
     if (bBox){
-        if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
+        if (d>a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) { /* acceptor is also a donor and vice versa? */
             /* return hbNo; */
             daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
         }
@@ -1456,7 +1461,7 @@ static int is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
     rda2 = iprod(r_da,r_da);
   
     if (bContact) {
-        if (daSwap)
+        if (daSwap && grpa == grpd)
             return hbNo;
         if (rda2 <= rc2){
             if (hb->bGem){
@@ -2182,7 +2187,7 @@ static void parallel_print(int *data, int nThreads)
         fprintf(stderr, "%-7i",data[i]);
 }
 
-static void normalizeACF(real *ct, real *gt, int len)
+static void normalizeACF(real *ct, real *gt, int nhb, int len)
 {
     real ct_fac, gt_fac;
     int i;
@@ -2190,7 +2195,8 @@ static void normalizeACF(real *ct, real *gt, int len)
     /* Xu and Berne use the same normalization constant */
 
     ct_fac = 1.0/ct[0];
-    gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
+    gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
+    
     printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
     for (i=0; i<len; i++)
     {
@@ -2397,8 +2403,7 @@ static void do_hbac(const char *fn,t_hbdata *hb,
         {
             sfree(dondata);
         }
-
-        normalizeACF(ct, NULL, nn);
+        normalizeACF(ct, NULL, 0, nn);
         snew(ctdouble, nn);
         snew(timedouble, nn);
         for (j=0; j<nn; j++)
@@ -2618,7 +2623,7 @@ static void do_hbac(const char *fn,t_hbdata *hb,
             sfree(dondata);
         }
 
-        normalizeACF(ct, NULL, nn);
+        normalizeACF(ct, NULL, 0, nn);
 
         fprintf(stderr, "\n\nACF successfully calculated.\n");
 
@@ -2764,7 +2769,7 @@ static void do_hbac(const char *fn,t_hbdata *hb,
         fprintf(stderr,"\n");
         sfree(h);
         sfree(g);
-        normalizeACF(ct, gt, nn);
+        normalizeACF(ct, ght, nhb, nn);
 
         /* Determine tail value for statistics */
         tail  = 0;
@@ -3954,7 +3959,7 @@ int gmx_hbond(int argc,char *argv[])
             int id,ia,hh,x,y;
       
             mat.nx=nframes;
-            mat.ny=(bContact ? hb->nrdist : hb->nrhb);
+            mat.ny=hb->nrhb;
 
             snew(mat.matrix,mat.nx);
             for(x=0; (x<mat.nx); x++)