Merge "[tools] g_sans - splited out fixes and added citation" into release-4-6

diff --git a/src/gmxlib/copyrite.c b/src/gmxlib/copyrite.c
index 00d326a4aa131bd1d0cbce9a6bcaa154361910b6..44bd448c335b11b82a5b92fbd66cab6628cce7c9 100644 (file)
--- a/src/gmxlib/copyrite.c
+++ b/src/gmxlib/copyrite.c
@@ -557,7 +557,12 @@ void please_cite(FILE *fp,const char *key)
       "V. Ballenegger, J.J. Cerda, and C. Holm",
       "How to Convert SPME to P3M: Influence Functions and Error Estimates",
       "J. Chem. Theory Comput.",
-      8, 2012, "936-947" }
+      8, 2012, "936-947" },
+    { "Garmay2012",
+      "Garmay Yu, Shvetsov A, Karelov D, Lebedev D, Radulescu A, Petukhov M, Isaev-Ivanov V",
+      "Correlated motion of protein subdomains and large-scale conformational flexibility of RecA protein filament",
+      "Journal of Physics: Conference Series",
+      340, 2012, "012094" }
   };
 #define NSTR (int)asize(citedb)
   

diff --git a/src/tools/gmx_sans.c b/src/tools/gmx_sans.c
index eac6ce2b0a2c9d0dbb2241a477d7953e590bf26c..311fa9c3a4faf43dc943bd57582cc0ba670c834d 100644 (file)
--- a/src/tools/gmx_sans.c
+++ b/src/tools/gmx_sans.c
@@ -78,7 +78,7 @@ int gmx_sans(int argc,char *argv[])
     static gmx_bool bPBC=TRUE;
     static real binwidth=0.2,grid=0.05; /* bins shouldnt be smaller then bond (~0.1nm) length */
     static real start_q=0.0, end_q=2.0, q_step=0.01;
-    static gmx_large_int_t  nmc=1048576;
+    static real mcover=-1;
     static unsigned int  seed=0;
 
     static const char *emode[]= { NULL, "direct", "mc", NULL };
@@ -94,8 +94,8 @@ int gmx_sans(int argc,char *argv[])
           "[HIDDEN]Binwidth (nm)" },
         { "-mode", FALSE, etENUM, {emode},
           "Mode for sans spectra calculation" },
-        { "-nmc", FALSE, etINT, {&nmc},
-          "Number of iterations for Monte-Carlo run"},
+        { "-mcover", FALSE, etREAL, {&mcover},
+          "Monte-Carlo coverage"},
         { "-method", FALSE, etENUM, {emethod},
           "[HIDDEN]Method for sans spectra calculation" },
         { "-pbc", FALSE, etBOOL, {&bPBC},
@@ -119,7 +119,7 @@ int gmx_sans(int argc,char *argv[])
   gmx_rmpbc_t  gpbc=NULL;
   gmx_bool  bTPX;
   gmx_bool  bFFT=FALSE, bDEBYE=FALSE;
-  gmx_bool  bMC=FALSE, bDIRECT=FALSE;
+  gmx_bool  bMC=FALSE;
   int        ePBC=-1;
   matrix     box;
   char       title[STRLEN];
@@ -148,18 +148,19 @@ int gmx_sans(int argc,char *argv[])
   parse_common_args(&argc,argv,PCA_BE_NICE,
                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 
+  /* check that binwidth not smaller than smallers distance */
+  check_binwidth(binwidth);
+
   /* Now try to parse opts for modes */
   switch(emethod[0][0]) {
   case 'd':
       bDEBYE=TRUE;
       switch(emode[0][0]) {
       case 'd':
-          bDIRECT=TRUE;
-          fprintf(stderr,"Using direct Debye method to calculate spectrum\n");
+          bMC=FALSE;
           break;
       case 'm':
           bMC=TRUE;
-          fprintf(stderr,"Using Monte Carlo Debye method to calculate spectrum\n");
           break;
       default:
           break;
@@ -167,7 +168,6 @@ int gmx_sans(int argc,char *argv[])
       break;
   case 'f':
       bFFT=TRUE;
-      fprintf(stderr,"Using FFT method\n");
       break;
   default:
       break;
@@ -175,8 +175,6 @@ int gmx_sans(int argc,char *argv[])
 
   if (!bDEBYE && !bFFT)
       gmx_fatal(FARGS,"Unknown method. Set pr or fft!\n");
-  if (!bDIRECT && !bMC)
-      gmx_fatal(FARGS,"Unknown mode for g(r) method set to direct or mc!");
   /* Try to read files */
   fnDAT = ftp2fn(efDAT,NFILE,fnm);
   fnTPX = ftp2fn(efTPX,NFILE,fnm);
@@ -202,28 +200,22 @@ int gmx_sans(int argc,char *argv[])
   }
 
   natoms=top->atoms.nr;
+
   if (bDEBYE) {
-      if (bDIRECT) {
-          /* calc pr */
-          pr = calc_radial_distribution_histogram(gsans,x,index,isize,binwidth,bMC,nmc,seed);
-      } else if (bMC) {
-          if (nmc>(gmx_large_int_t)floor(0.5*isize*(isize-1))) {
-              fprintf(stderr,"Number of mc iteration larger then number of pairs in index group. Switching to direct method!\n");
-              bMC=FALSE;
-              bDIRECT=TRUE;
-              pr = calc_radial_distribution_histogram(gsans,x,index,isize,binwidth,bMC,nmc,seed);
-          } else {
-              pr = calc_radial_distribution_histogram(gsans,x,index,isize,binwidth,bMC,nmc,seed);
-          }
+      if (bMC) {
+          fprintf(stderr,"Using Monte Carlo Debye method to calculate spectrum\n");
       } else {
-          gmx_fatal(FARGS,"Unknown method!\n");
+          fprintf(stderr,"Using direct Debye method to calculate spectrum\n");
       }
   } else if (bFFT) {
-      gmx_fatal(FARGS,"Not implented!\n");
+      gmx_fatal(FARGS,"Not implented!");
   } else {
-      gmx_fatal(FARGS,"Whats this!?\n");
+      gmx_fatal(FARGS,"Whats this!");
   }
 
+  /*  realy calc p(r) */
+  pr = calc_radial_distribution_histogram(gsans,x,box,index,isize,binwidth,bMC,mcover,seed);
+
   /* prepare pr.xvg */
   fp = xvgropen(opt2fn_null("-pr",NFILE,fnm),"G(r)","Distance (nm)","Probability",oenv);
   for(i=0;i<pr->grn;i++)
@@ -240,7 +232,8 @@ int gmx_sans(int argc,char *argv[])
 
   sfree(pr);
 
+  please_cite(stdout,"Garmay2012");
   thanx(stderr);
-  
+
   return 0;
 }

diff --git a/src/tools/nsfactor.c b/src/tools/nsfactor.c
index 0393249f1eb59c3268fa4a1a860ed2084adaaffd..79ebdbcebcdbfe37a066ec07cb502c9b2e5076af 100644 (file)
--- a/src/tools/nsfactor.c
+++ b/src/tools/nsfactor.c
@@ -46,6 +46,12 @@
 #include "vec.h"
 #include "nsfactor.h"
 
+void check_binwidth(real binwidth) {
+    real smallest_bin=0.1;
+    if (binwidth<smallest_bin)
+        gmx_fatal(FARGS,"Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
+}
+
 void normalize_probability(int n,double *a){
     int i;
     double norm=0.0;
@@ -131,12 +137,21 @@ gmx_sans_t *gmx_sans_init (t_topology *top, gmx_nentron_atomic_structurefactors_
     return (gmx_sans_t *) gsans;
 }
 
-gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (gmx_sans_t *gsans, rvec *x, atom_id *index, int isize, double binwidth, gmx_bool bMC, gmx_large_int_t nmc, unsigned int seed) {
+gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
+                            gmx_sans_t *gsans,
+                            rvec *x,
+                            matrix box,
+                            atom_id *index,
+                            int isize,
+                            double binwidth,
+                            gmx_bool bMC,
+                            real mcover,
+                            unsigned int seed) {
     gmx_radial_distribution_histogram_t    *pr=NULL;
-    rvec        xmin, xmax;
+    rvec        dist;
     double      rmax;
-    int         i,j,d;
-    int         mc;
+    int         i,j;
+    gmx_large_int_t   mc=0,max;
     gmx_rng_t   rng=NULL;
 
     /* allocate memory for pr */
@@ -144,18 +159,14 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (gmx_san
     /* set some fields */
     pr->binwidth=binwidth;
 
-    /* Lets try to find min and max distance */
-    for(d=0;d<3;d++) {
-        xmax[d]=x[index[0]][d];
-        xmin[d]=x[index[0]][d];
-    }
-
-    for(i=1;i<isize;i++)
-        for(d=0;d<3;d++)
-            if (xmax[d]<x[index[i]][d]) xmax[d]=x[index[i]][d]; else
-                if (xmin[d]>x[index[i]][d]) xmin[d]=x[index[i]][d];
+    /*
+    * create max dist rvec
+    * dist = box[xx] + box[yy] + box[zz]
+    */
+    rvec_add(box[XX],box[YY],dist);
+    rvec_add(box[ZZ],dist,dist);
 
-    rmax=sqrt(distance2(xmax,xmin));
+    rmax=norm(dist);
 
     pr->grn=(int)floor(rmax/pr->binwidth)+1;
     rmax=pr->grn*pr->binwidth;
@@ -163,17 +174,21 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (gmx_san
     snew(pr->gr,pr->grn);
 
     if(bMC) {
-        /* Use several independent mc runs to collect better statistics */
-        for(d=0;d<(int)floor(nmc/524288);d++) {
-            rng=gmx_rng_init(seed);
-            for(mc=0;mc<524288;mc++) {
-                i=(int)floor(gmx_rng_uniform_real(rng)*isize);
-                j=(int)floor(gmx_rng_uniform_real(rng)*isize);
-                if(i!=j)
-                    pr->gr[(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
-            }
-            gmx_rng_destroy(rng);
+        /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
+        if (mcover==-1) {
+            max=(gmx_large_int_t)floor(0.5*0.01*isize*(isize-1));
+        } else {
+            max=(gmx_large_int_t)floor(0.5*mcover*isize*(isize-1));
+        }
+        rng=gmx_rng_init(seed);
+        while(mc<max) {
+            i=(int)floor(gmx_rng_uniform_real(rng)*isize);
+            j=(int)floor(gmx_rng_uniform_real(rng)*isize);
+            if(i!=j)
+                pr->gr[(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
+            mc++;
         }
+        gmx_rng_destroy(rng);
     } else {
         for(i=0;i<isize;i++)
             for(j=0;j<i;j++)

diff --git a/src/tools/nsfactor.h b/src/tools/nsfactor.h
index 6563ea3706f1bb9fe8b4254c9486e0fd0c8db550..7f2d9cc59e78c051f3d7a5559bcbee5e8dbc1ef1 100644 (file)
--- a/src/tools/nsfactor.h
+++ b/src/tools/nsfactor.h
@@ -72,6 +72,8 @@ typedef struct gmx_static_structurefator_t {
     double  qstep; /* q increment */
 } gmx_static_structurefator_t;
 
+void check_binwidth(real binwidth);
+
 void normalize_probability(int n, double *a);
 
 gmx_nentron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn);
@@ -79,11 +81,13 @@ gmx_nentron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const ch
 gmx_sans_t *gmx_sans_init(t_topology *top, gmx_nentron_atomic_structurefactors_t *gnsf);
 
 gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram  ( gmx_sans_t *gsans,
-                            rvec *x, atom_id *index,
+                            rvec *x,
+                            matrix box,
+                            atom_id *index,
                             int isize,
                             double binwidth,
                             gmx_bool bMC,
-                            gmx_large_int_t nmc,
+                            real mcover,
                             unsigned int seed);
 
 gmx_static_structurefator_t *convert_histogram_to_intensity_curve (gmx_radial_distribution_histogram_t *pr, double start_q, double end_q, double q_step);