HI Carla,
I wrote a similar code, see attached. But it is written for my
condition. You should modify it accordingly.
regards,
Baofu Qiao
On 07/12/2011 02:04 PM, Carla Jamous wrote:
Dear gmx-users,
I used trjorder in order to study the water molecules that are closer
than 5 A from my protein.
trjorder -s structure.tpr -f traj.xtc -n prot_water.ndx -o ordered.pdb
-nshell nshell_.xvg -r 0.5 -b 0 -e 5000
But now I need to analyse the residence time of a water molecule, I
mean the number of times a single water molecule stays in a radius of
5 A of the protein and divide this number by the total number of
conformations, in order to have a pourcentage value.
Please is there any gromacs tool able to do this calculation or else
does anyone have an idea of how to do that?
Thank you
Carla
/*
* $Id: gmx_density.c,v 1.9.2.1 2009/01/18 23:06:27 lindahl Exp $
*
* This source code is part of
*
* G R O M A C S
*
* GROningen MAchine for Chemical Simulations
*
* VERSION 3.2.0
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team,
* check out http://www.gromacs.org for more information.
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* If you want to redistribute modifications, please consider that
* scientific software is very special. Version control is crucial -
* bugs must be traceable. We will be happy to consider code for
* inclusion in the official distribution, but derived work must not
* be called official GROMACS. Details are found in the README & COPYING
* files - if they are missing, get the official version at www.gromacs.org.
*
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
*
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
* Green Red Orange Magenta Azure Cyan Skyblue
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
#include <ctype.h>
#include "sysstuff.h"
#include "string.h"
#include "string2.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
#include "gstat.h"
#include "vec.h"
#include "xvgr.h"
#include "pbc.h"
#include "copyrite.h"
#include "futil.h"
#include "statutil.h"
#include "index.h"
#include "tpxio.h"
#include "physics.h"
int main(int argc,char *argv[])
{
const char *desc[] = {
"This is to calculate the residence time of particles in certain slab region."
};
output_env_t oenv;
static int axis = 2; /* normal to memb. default z */
static const char *axtitle="Z";
static real ZMin=-1.0, ZMax=-1.0;
t_pargs pa[] = {
{ "-d", FALSE, etSTR, {&axtitle},
"Dimension for Residence Time calculation: X, Y or Z." },
{ "-z0", FALSE, etREAL, {&ZMin},
"Minimum of Z coordinate (nm). <=0 means 0."},
{ "-z1", FALSE, etREAL, {&ZMax},
"Maximum of Z coordinate (nm). <=0 measn BOX[Z][Z]."}
};
#define NPA asize(pa)
gmx_bool bTop;
t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
{ efTPX, "-s", NULL, ffREAD },
{ efNDX, "-n", NULL, ffOPTRD },
{ efXVG,"-o","residence",ffWRITE },
};
#define NFILE asize(fnm)
FILE *fp_resid;
t_topology *top;
atom_id *index;
rvec *x;
matrix box;
t_trxstatus *status;
gmx_rmpbc_t gpbc=NULL;
char *grpname,title[256],*ylabel=NULL;
int ePBC,natoms,gnx,nr_frames=0;
int ax1=0, ax2=0;
real z;
int i,*lifetime,*Cr,length=10000;
real t,t0,t1,dt=0;
time_t tstart,tend;
int all_time,hour,min,sec;
//CopyRight(stderr,argv[0]);
parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
/* Calculate axis */
axis = toupper(axtitle[0]) - 'X';
switch(axis) {
case 0: ax1=1; ax2=2; break;
case 1: ax1=0; ax2=2; break;
case 2: ax1=0; ax2=1; break;
default: gmx_fatal(FARGS,"\nInvalid axes. Terminating...\n\n");
}
//bTop = read_tps_conf(ftp2fn(efTPX,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
fprintf(stdout,"\nPlease select (only) one group:\n");
get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
if ((natoms = read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box)) == 0)
gmx_fatal(FARGS,"\nCould not read coordinates from statusfile!\n\n");
if (ZMin<=0) ZMin = 0;
if (ZMax<=0 || ZMax>box[axis][axis]) ZMax = box[axis][axis];
fprintf(stdout,"\nRegime of [%g,%g]nm employed!\n",ZMin,ZMax);
snew(lifetime,gnx);
for (i=0; i<gnx; i++) lifetime[i]=0;
snew(Cr,length);
for (i=0; i<length; i++) Cr[i]=0;
time(&tstart);
t0 = t;
gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
/*********** Start processing trajectory ***********/
do {
gmx_rmpbc(gpbc,natoms,box,x);
for (i=0; i<gnx; i++) {
z = x[index[i]][axis];
if (z < 0) //This is not used in the system with "nwall=2"
z += box[axis][axis];
else if (z > box[axis][axis])
z -= box[axis][axis];
if (z>=ZMin && z<=ZMax) {
lifetime[i]++ ;
if (lifetime[i]>=length) {
length+=100;
srenew(Cr,length);
}
Cr[lifetime[i]-1]++ ;
} else {
lifetime[i] = 0;
}
}
if (nr_frames == 1) {
t1 = t;
dt = t1-t0;
}
nr_frames++;
} while (read_next_x(oenv,status,&t,natoms,x,box));
gmx_rmpbc_done(gpbc);
close_trj(status);
sfree(x);
fprintf(stdout,"\nRead %d frames from trajectory. Printing......\n",nr_frames);
/*Print out the results */
sprintf(title,"Residence probability: %s",grpname);
fp_resid=xvgropen(opt2fn("-o",NFILE,fnm),title,"t (ps)","C\\sR\\N(t) (%)",oenv);
fprintf(fp_resid, "#Coordinate range in %s Dimension: %g-%g (nm)\n",axtitle,ZMin,ZMax);
for (i=0; i<length; i++) {
fprintf(fp_resid,"\n%6g %8g",i*dt,Cr[i]*1.0/Cr[0]);
if (Cr[i] < Cr[0]/100.0) break;
}
fprintf(fp_resid,"\n");
fclose(fp_resid);
time(&tend);
all_time = (int)(tend-tstart);
hour = (int) (all_time/3600);
min = (int)(( all_time-hour*3600) /60);
sec = all_time - hour*3600 - min*60;
fprintf(stderr,"\nIt runs %d:%d:%d (hh:mm:ss) totally\n",hour,min,sec);
thanx(stderr);
return 0;
}
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
