
1 
/* * mode: c; tabwidth: 4; indenttabsmode: nil; cbasicoffset: 4; cfilestyle: "stroustrup"; *


2 
*


3 
* This source code is part of


4 
*


5 
* G R O M A C S


6 
*


7 
* GROningen MAchine for Chemical Simulations


8 
*


9 
* VERSION 4.0.99


10 
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.


11 
* Copyright (c) 19912000, University of Groningen, The Netherlands.


12 
* Copyright (c) 20012008, The GROMACS development team,


13 
* check out http://www.gromacs.org for more information.


14 


15 
* This program is free software; you can redistribute it and/or


16 
* modify it under the terms of the GNU General Public License


17 
* as published by the Free Software Foundation; either version 2


18 
* of the License, or (at your option) any later version.


19 
*


20 
* If you want to redistribute modifications, please consider that


21 
* scientific software is very special. Version control is crucial 


22 
* bugs must be traceable. We will be happy to consider code for


23 
* inclusion in the official distribution, but derived work must not


24 
* be called official GROMACS. Details are found in the README & COPYING


25 
* files  if they are missing, get the official version at www.gromacs.org.


26 
*


27 
* To help us fund GROMACS development, we humbly ask that you cite


28 
* the papers on the package  you can find them in the top README file.


29 
*


30 
* For more info, check our website at http://www.gromacs.org


31 
*


32 
* And Hey:


33 
* Groningen Machine for Chemical Simulation


34 
*/


35 


36 
/*


37 
* gmx_correl.c (c) 2010 Alexey Shvetsov


38 
*/


39 


40 
#ifdef HAVE_CONFIG_H


41 
#include <config.h>


42 
#endif


43 
#include <math.h>


44 
#include <string.h>


45 
#include <time.h>


46 


47 
#if ((defined WIN32  defined _WIN32  defined WIN64  defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)


48 
#include <direct.h>


49 
#include <io.h>


50 
#endif


51 


52 
#include "statutil.h"


53 
#include "sysstuff.h"


54 
#include "typedefs.h"


55 
#include "smalloc.h"


56 
#include "macros.h"


57 
#include "vec.h"


58 
#include "pbc.h"


59 
#include "copyrite.h"


60 
#include "futil.h"


61 
#include "statutil.h"


62 
#include "index.h"


63 
#include "confio.h"


64 
#include "trnio.h"


65 
#include "mshift.h"


66 
#include "xvgr.h"


67 
#include "do_fit.h"


68 
#include "rmpbc.h"


69 
#include "txtdump.h"


70 
#include "matio.h"


71 
#include "physics.h"


72 
#include "gmx_ana.h"


73 


74 
int gmx_correl(int argc,char *argv[]) {


75 
const char *desc[] = {


76 
"[TT]g_correlation[tt] calculates and the correlation matrix.",


77 
"All structures are fitted to the structure in the structure file.",


78 
"When this is not a run input file periodicity will not be taken into",


79 
"account.",


80 
"[PAR]",


81 
"When the same atoms are used for the fit and the correlation analysis,",


82 
"the reference structure for the fit is written first with t=1.",


83 
"[PAR]",


84 
"Option [TT]ascii[tt] writes the whole correlation matrix to",


85 
"an ASCII file.",


86 
"[PAR]",


87 
"Option [TT]xpm[tt] writes the whole covariance matrix to an xpm file.",


88 
"[PAR]"


89 
};


90 
static bool bFit=TRUE,bRef=FALSE,bPBC=TRUE;


91 
t_pargs pa[] = {


92 
{ "fit", FALSE, etBOOL, {&bFit},


93 
"Fit to a reference structure"},


94 
{ "ref", FALSE, etBOOL, {&bRef},


95 
"Use the deviation from the conformation in the structure file instead of from the average" },


96 
{ "pbc", FALSE, etBOOL, {&bPBC},


97 
"Apply corrections for periodic boundary conditions" }


98 
};


99 
FILE *out;


100 
int status;


101 
t_topology top;


102 
int ePBC;


103 
t_atoms *atoms;


104 
rvec *xread,*xref;


105 
rvec *xa,*xava;


106 
rvec *xb,*xavb;


107 
matrix box,zerobox;


108 
real inv_nframes;


109 
real t,tstart,tend;


110 
real *w_rls=NULL;


111 
real min,max,*axisa,*axisb;


112 
real **mat;


113 
real *r2a,*r2b;


114 
int natoms,nat,nframes0,nframes,nlevels;


115 
int natomsa,natomsb;


116 
gmx_large_int_t i,j;


117 
const char *fitfile,*trxfile,*ndxfile;


118 
const char *asciifile,*xpmfile,*averfilea,*averfileb;


119 
char str[STRLEN],*fitname;


120 
char *ananamea,*ananameb;


121 
int d,nfit;


122 
atom_id *ifit;


123 
atom_id *indexa,*indexb;


124 
t_rgb rlo,rmi,rhi;


125 
output_env_t oenv;


126 


127 
t_filenm fnm[] = {


128 
{ efTRX, "f", NULL, ffREAD },


129 
{ efTPS, NULL, NULL, ffREAD },


130 
{ efNDX, NULL, NULL, ffOPTRD },


131 
{ efSTO, "ava", "averagea.pdb", ffWRITE },


132 
{ efSTO, "avb", "averageb.pdb", ffWRITE },


133 
{ efDAT, "ascii","correlation", ffOPTWR },


134 
{ efXPM, "xpm","correlation", ffOPTWR }


135 
};


136 
#define NFILE asize(fnm)


137 


138 
CopyRight(stderr,argv[0]);


139 
parse_common_args(&argc,argv,PCA_CAN_TIME  PCA_TIME_UNIT  PCA_BE_NICE,


140 
NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);


141 


142 
clear_mat(zerobox);


143 


144 
fitfile = ftp2fn(efTPS,NFILE,fnm);


145 
trxfile = ftp2fn(efTRX,NFILE,fnm);


146 
ndxfile = ftp2fn_null(efNDX,NFILE,fnm);


147 
averfilea = ftp2fn(efSTO,NFILE,fnm);


148 
averfileb = ftp2fn(efSTO,NFILE,fnm);


149 
asciifile = opt2fn_null("ascii",NFILE,fnm);


150 
xpmfile = opt2fn_null("xpm",NFILE,fnm);


151 


152 
read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE);


153 
atoms=&top.atoms;


154 


155 
if (bFit) {


156 
printf("\nChoose a group for the least squares fit\n");


157 
get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname);


158 
if (nfit < 3)


159 
gmx_fatal(FARGS,"Need >= 3 points to fit!\n");


160 
} else


161 
nfit=0;


162 
printf("\nChoose two groups for the correlation analysis\n");


163 
get_index(atoms,ndxfile,1,&natomsa,&indexa,&ananamea);


164 
get_index(atoms,ndxfile,1,&natomsb,&indexb,&ananameb);


165 


166 
if (bFit) {


167 
snew(w_rls,atoms>nr);


168 
for(i=0; (i<nfit); i++)


169 
w_rls[ifit[i]]=1.0;


170 
}


171 


172 
/* Prepare reference frame */


173 
if (bPBC)


174 
rm_pbc(&(top.idef),ePBC,atoms>nr,box,xref,xref);


175 
if (bFit)


176 
reset_x(nfit,ifit,atoms>nr,NULL,xref,w_rls);


177 


178 
snew(xa,natomsa);


179 
snew(xb,natomsb);


180 
snew(xava,natomsa);


181 
snew(xavb,natomsb);


182 


183 
fprintf(stderr,"Calculating the average structure ...\n");


184 
nframes0 = 0;


185 
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);


186 
if (nat != atoms>nr)


187 
fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);


188 
do {


189 
nframes0++;


190 
/* calculate x: a fitted struture of the selected atoms */


191 
if (bPBC)


192 
rm_pbc(&(top.idef),ePBC,nat,box,xread,xread);


193 
if (bFit) {


194 
reset_x(nfit,ifit,nat,NULL,xread,w_rls);


195 
do_fit(nat,w_rls,xref,xread);


196 
}


197 
for (i=0; i<natomsa; i++)


198 
rvec_inc(xava[i],xread[indexa[i]]);


199 
for (i=0; i<natomsb; i++)


200 
rvec_inc(xavb[i],xread[indexb[i]]);


201 
} while (read_next_x(oenv,status,&t,nat,xread,box));


202 
close_trj(status);


203 


204 
inv_nframes = 1.0/nframes0;


205 
for(i=0; i<natomsa; i++)


206 
for(d=0; d<DIM; d++) {


207 
xava[i][d] *= inv_nframes;


208 
xread[indexa[i]][d] = xava[i][d];


209 
}


210 
write_sto_conf_indexed(opt2fn("ava",NFILE,fnm),"Average structure",


211 
atoms,xread,NULL,epbcNONE,zerobox,natomsa,indexa);


212 
for(i=0; i<natomsb; i++)


213 
for(d=0; d<DIM; d++) {


214 
xavb[i][d] *= inv_nframes;


215 
xread[indexb[i]][d] = xavb[i][d];


216 
}


217 
write_sto_conf_indexed(opt2fn("avb",NFILE,fnm),"Average structure",


218 
atoms,xread,NULL,epbcNONE,zerobox,natomsb,indexb);


219 
sfree(xread);


220 


221 
fprintf(stderr,"Constructing correlation matrix (%dx%d) ...\n",(int)natomsa,(int)natomsb);


222 
nframes=0;


223 
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);


224 
tstart = t;


225 
snew(r2a,natomsa);


226 
snew(r2b,natomsb);


227 
/* Allocate 2d matrix */


228 
snew(mat,natomsa);


229 
for (i=0;i<natomsb;i++)


230 
snew(mat[i],natomsb);


231 
do {


232 
nframes++;


233 
tend = t;


234 
/* calculate x: a (fitted) structure of the selected atoms */


235 
if (bPBC)


236 
rm_pbc(&(top.idef),ePBC,nat,box,xread,xread);


237 
if (bFit) {


238 
reset_x(nfit,ifit,nat,NULL,xread,w_rls);


239 
do_fit(nat,w_rls,xref,xread);


240 
}


241 
if (bRef) {


242 
for (i=0; i<natomsa; i++)


243 
rvec_sub(xread[indexa[i]],xref[indexa[i]],xa[i]);


244 
for (i=0; i<natomsb;i++)


245 
rvec_sub(xread[indexb[i]],xref[indexb[i]],xb[i]);


246 
}


247 
else {


248 
for (i=0; i<natomsa; i++)


249 
rvec_sub(xread[indexa[i]],xava[i],xa[i]);


250 
for (i=0; i<natomsb; i++)


251 
rvec_sub(xread[indexb[i]],xavb[i],xb[i]);


252 
}


253 
/* find actual rmsf */


254 
for (i=0; i<natomsa; i++)


255 
r2a[i]+=norm2(xa[i]);


256 
for (i=0; i<natomsb; i++)


257 
r2b[i]+=norm2(xb[i]);


258 
/* TODO: Optimize calculation of correlation matrix */


259 
for (i=0; i<natomsa; i++) {


260 
for (j=0; j<natomsb; j++) {


261 
mat[i][j] += iprod(xa[i],xb[j]);


262 
}


263 
}


264 
} while (read_next_x(oenv,status,&t,nat,xread,box) &&


265 
(bRef  nframes < nframes0));


266 
close_trj(status);


267 
fprintf(stderr,"Read %d frames\n",nframes);


268 
/* Correct mat so it will be real correlation one */


269 
for (i=0;i<natomsa;i++) {


270 
for (j=0;j<natomsb;j++) {


271 
mat[i][j] = mat[i][j]/sqrt(r2a[i]*r2b[j]);


272 
}


273 
}


274 
if (asciifile) {


275 
out = ffopen(asciifile,"w");


276 
for (i=0; i<natomsa; i++) {


277 
for (j=0; j<natomsb; j++) {


278 
fprintf(out,"%g ",mat[i][j]);


279 
}


280 
fprintf(out,"\n");


281 
}


282 
ffclose(out);


283 
}


284 


285 
if (xpmfile) {


286 
min = 0;


287 
max = 0;


288 
for (i=0; i<natomsa; i++) {


289 
for (j=0; j<natomsb; j++) {


290 
if (mat[i][j] < min)


291 
min = mat[i][j];


292 
if (mat[i][j] > max)


293 
max = mat[i][j];


294 
}


295 
}


296 
snew(axisa,natomsa);


297 
snew(axisb,natomsb);


298 
for(i=0;i<natomsa;i++)


299 
axisa[i] = i+1;


300 
for(i=0;i<natomsb;i++)


301 
axisb[i] = i+1;


302 
rlo.r = 0; rlo.g = 0; rlo.b = 1;


303 
rmi.r = 1; rmi.g = 1; rmi.b = 1;


304 
rhi.r = 1; rhi.g = 0; rhi.b = 0;


305 
out = ffopen(xpmfile,"w");


306 
nlevels = 80;


307 
write_xpm3(out,0,"Correlation", "",


308 
"Residue Number","Residue Number",natomsa,natomsb,axisa,axisb,


309 
mat,min,0.0,max,rlo,rmi,rhi,&nlevels);


310 
ffclose(out);


311 
sfree(axisa);


312 
sfree(axisb);


313 
}


314 
for (i=0;i<natomsb;i++)


315 
sfree(mat[i]);


316 
sfree(mat);


317 
sfree(r2a);


318 
sfree(r2b);


319 


320 
thanx(stderr);


321 


322 
return 0;


323 
}

0 


