[gmx-users] problems with g_wham

[Justin A. Lemkul]

Gavin Melaugh wrote:
> Hi Chris
> 
> Thanks very much for your help. The reason I used pull_init1 = 0 is
> because this is the value that is used in the tutorial. Maybe I am
> misinterpreting its meaning. In the tutorial it says that you do not
> have to specify a value separately for each configuration. I am now
> thinking that in my case this does not apply ? Should I set pull_init to
> the distance between the molecules at the first configuration?
> 

To clarify, the tutorial sets "pull_init1 = 0" because it is used in conjunction 
with "pull_start = yes," in order to take the starting COM distance as the 
reference.  Check the grompp output to verify that you're getting what you think 
you should be.  Otherwise, you'd have to set "pull_start = no" and modify your 
.mdp file to reflect the proper value of pull_init1 for every window.

-Justin

> Many thanks
> 
> Gavin
> chris.neale at utoronto.ca wrote:
>> 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
>>
>>
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================