[gmx-users] problems with g_wham

[chris.neale at utoronto.ca]

Gavin, let's take g_wham out of the picture for now to simplify  
things. First, you umbrella is pretty strange. You want them to have a  
distance of 0.0 nm (pull_init1 = 0)???

So let's look at making a histogram from your raw data.

1. cat pullx.xvg |awk '{print sqrt($1,$5*$5+$6*$6+$7*$7)}' > my.dist
2. ./histo my.dist 0.01 > my.histo
** select the bin width 0.01 so that the curve is still smooth but  
also keep the number small for precision.
** compile the histo.c code below to get ./histo as follows (gcc  
histo.c -o histo -lm)
3. Make a plot of the histogram and post it on the internet (e.g.  
using photobucket website)
4. reply with a link to your histogram.

#### PS: i make no claim about the accuracy of the histo.c program so  
if you use it for other activities then its at your own risk. In fact,  
I think that there are cases where the very fist value is reported  
incorrectly, but that should not matter here.

#include <stdio.h>
#include <stdlib.h>
#include "math.h"

const int LINE_SIZE=1000;

void showUsage(const char *c){
   printf("Usage: %s <input_file> <binwidth> [<min> <max>]\n",c);
}

int main(int argn, char *args[]){
   char linein[LINE_SIZE];
   FILE *f;
   int *data;
   int nbin,tot,i;
   double min,max,bin,c;

   if ((argn!=3 &&argn!=5)|| sscanf(args[2],"%lf",&bin)!=1){
     showUsage(args[0]);
     exit(1);
   }

   if((f=fopen(args[1],"r"))==NULL){
     printf("Error: unable to open file\n");
     showUsage(args[0]);
     exit(1);
   }
   while(fgets(linein,LINE_SIZE,f)!=NULL){
     if((sscanf(linein,"%lf",&c))!=1)continue;
     if(c<min)min=c;
     if(c>max)max=c;
   }
   fclose(f);

   if (argn==5){
     if(sscanf(args[3],"%lf",&min)!=1 || sscanf(args[4],"%lf",&max)!=1){
       showUsage(args[0]);
       exit(1);
     }
   }

   min=floor(min/bin)*bin;
   max=ceil(max/bin)*bin;

   nbin=(int)ceil((max-min)/bin);
   data=(int *)malloc((nbin+1)*sizeof(int));
   if(data==NULL){
     printf("Memory allocation error\n");
     exit(1);
   }
   for(i=0; i<=nbin; ++i){
     data[i]=0;
   }

   tot=0;
   if((f=fopen(args[1],"r"))==NULL){
     printf("Error: unable to open file on second attempt ...  
programming error?\n");
     exit(1);
   }
   while(fgets(linein,LINE_SIZE,f)!=NULL){
     if((sscanf(linein,"%lf",&c))!=1)continue;
     ++data[(int)ceil((c-min)/bin)];
     ++tot;
   }
   fclose(f);

   //Bin 0 only holds the absolute minimum value
   i=1;
    
printf("%lf\t%lf\n",min+bin*((double)i-0.5),(double)(data[0]+data[1])/(double)tot);
   for(i=2; i<=nbin; ++i){
     printf("%lf\t%lf\n",min+bin*((double)i-0.5),(double)data[i]/(double)tot);
   }
}



-- original message --

Hi Guys

I am having great difficulty generating a histogram of the correct shape
from g_wham. I am running umbrella sampling on one configuration fo two
molecules at a fixed distance using the following
pull parameters. ( I am just running one configuration to begin with so
that can make sure that the shape of my distribution is correct).

pull        = umbrella
pull_geometry = distance
pull_dim = Y Y Y
pull_start = yes
pull_ngroups = 1
pull_group0 = cage_1
pull_group1 = cage_2
pull_init1 = 0
pull_rate1 = 0.0
pull_k1 = 2000
pull_nstxout = 100
pull_nstfout = 100

I have set the force constant quite high in the hope that the system
will sample the equilibrium distance much more often than the distances
accessible form the harmonic potential. In the plot of count vs distance
I get an upside down bell curve in which there is a larger distribution
at the extremeties of the curve. Does anyone know why this may be
happening? If so please let me know. Am I right in assuming that the x
axis in the histo.xvg file is the distance between the centres of masses
of the two molecules?

Cheers

Gavin