From 021c8b41cf3b0d47a8adb995d7ce53bcb51188f3 Mon Sep 17 00:00:00 2001 From: Alexey Shvetsov Date: Tue, 15 Jun 2010 03:20:47 +0400 Subject: [PATCH] Add new tool g_correlation. This tool computes Pearsons coefficients for atomic displacements and constructs {cross,auto}correlation matrix from them. --- src/tools/CMakeLists.txt | 4 +- src/tools/Makefile.am | 2 +- src/tools/g_correlation.c | 55 ++++++++ src/tools/gmx_correlation.c | 322 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 380 insertions(+), 3 deletions(-) create mode 100644 src/tools/g_correlation.c create mode 100644 src/tools/gmx_correlation.c diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 5fc666c..d11cbd5 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -24,7 +24,7 @@ add_library(gmxana gmx_sgangle.c gmx_sorient.c gmx_spol.c gmx_tcaf.c gmx_traj.c gmx_velacc.c gmx_helixorient.c gmx_clustsize.c gmx_mdmat.c gmx_wham.c - correl.c gmx_sham.c gmx_nmtraj.c + correl.c gmx_sham.c gmx_nmtraj.c gmx_correlation.c gmx_trjconv.c gmx_trjcat.c gmx_trjorder.c gmx_xpm2ps.c gmx_editconf.c gmx_genbox.c gmx_genion.c gmx_genconf.c gmx_genpr.c gmx_eneconv.c gmx_vanhove.c gmx_wheel.c @@ -43,7 +43,7 @@ set(GMX_TOOLS_PROGRAMS g_angle g_bond g_bundle g_chi g_cluster g_confrms g_covar g_current g_density g_densmap g_dih g_dielectric g_helixorient g_principal g_dipoles g_disre g_dist - g_dyndom g_enemat g_energy g_lie g_filter g_gyrate + g_dyndom g_enemat g_energy g_lie g_filter g_gyrate g_correlation g_h2order g_hbond g_helix g_mindist g_msd g_morph g_nmeig g_nmens g_order g_polystat g_potential g_rama g_rdf g_rms g_rmsf g_rotacf g_saltbr g_sas g_select g_sgangle g_sham g_sorient diff --git a/src/tools/Makefile.am b/src/tools/Makefile.am index 609667b..fc6e262 100644 --- a/src/tools/Makefile.am +++ b/src/tools/Makefile.am @@ -44,7 +44,7 @@ libgmxana@LIBSUFFIX@_la_SOURCES = \ gmx_editconf.c gmx_genbox.c gmx_genion.c gmx_genconf.c \ gmx_genpr.c gmx_eneconv.c gmx_vanhove.c gmx_wheel.c \ addconf.c addconf.h gmx_tune_pme.c \ - calcpot.c calcpot.h edittop.c + calcpot.c calcpot.h edittop.c gmx_correlation.c bin_PROGRAMS = \ do_dssp editconf eneconv \ diff --git a/src/tools/g_correlation.c b/src/tools/g_correlation.c new file mode 100644 index 0000000..9c3b124 --- /dev/null +++ b/src/tools/g_correlation.c @@ -0,0 +1,55 @@ +/* + * + * 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 +#endif + +#include + + + +/* This is just a wrapper binary. +* The code that used to be in g_covar.c is now in gmx_covar.c, +* where the old main function is called gmx_covar(). +*/ +int +main(int argc, char *argv[]) +{ + gmx_correlation(argc,argv); + return 0; +} + + + diff --git a/src/tools/gmx_correlation.c b/src/tools/gmx_correlation.c new file mode 100644 index 0000000..ea405ac --- /dev/null +++ b/src/tools/gmx_correlation.c @@ -0,0 +1,322 @@ +/* + * + * 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 +#endif +#include +#include +#include + +#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__) +#include +#include +#endif + +#include "statutil.h" +#include "sysstuff.h" +#include "typedefs.h" +#include "smalloc.h" +#include "macros.h" +#include "vec.h" +#include "pbc.h" +#include "copyrite.h" +#include "futil.h" +#include "statutil.h" +#include "index.h" +#include "confio.h" +#include "trnio.h" +#include "mshift.h" +#include "xvgr.h" +#include "do_fit.h" +#include "rmpbc.h" +#include "txtdump.h" +#include "matio.h" +#include "physics.h" +#include "gmx_ana.h" + +int gmx_correlation(int argc,char *argv[]) +{ + const char *desc[] = { + "[TT]g_correlation[tt] calculates and the correlation matrix.", + "All structures are fitted to the structure in the structure file.", + "When this is not a run input file periodicity will not be taken into", + "account.", + "[PAR]", + "When the same atoms are used for the fit and the correlation analysis,", + "the reference structure for the fit is written first with t=-1.", + "[PAR]", + "Option [TT]-ascii[tt] writes the whole correlation matrix to", + "an ASCII file.", + "[PAR]", + "Option [TT]-xpm[tt] writes the whole covariance matrix to an xpm file.", + "[PAR]" + }; + static bool bFit=TRUE,bRef=FALSE,bPBC=TRUE; + t_pargs pa[] = { + { "-fit", FALSE, etBOOL, {&bFit}, + "Fit to a reference structure"}, + { "-ref", FALSE, etBOOL, {&bRef}, + "Use the deviation from the conformation in the structure file instead of from the average" }, + { "-pbc", FALSE, etBOOL, {&bPBC}, + "Apply corrections for periodic boundary conditions" } + }; + FILE *out; + int status; + t_topology top; + int ePBC; + t_atoms *atoms; + rvec *xread,*xref; + rvec *xa,*xava; + rvec *xb,*xavb; + matrix box,zerobox; + real inv_nframes; + real t,tstart,tend; + real *w_rls=NULL; + real min,max,*axisa,*axisb; + real **mat; + real *r2a,*r2b; + int natoms,nat,nframes0,nframes,nlevels; + int natomsa,natomsb; + gmx_large_int_t i,j; + const char *fitfile,*trxfile,*ndxfile; + const char *asciifile,*xpmfile,*averfilea,*averfileb; + char str[STRLEN],*fitname; + char *ananamea,*ananameb; + int d,nfit; + atom_id *ifit; + atom_id *indexa,*indexb; + t_rgb rlo,rmi,rhi; + output_env_t oenv; + + t_filenm fnm[] = { + { efTRX, "-f", NULL, ffREAD }, + { efTPS, NULL, NULL, ffREAD }, + { efNDX, NULL, NULL, ffOPTRD }, + { efSTO, "-ava", "averagea.pdb", ffWRITE }, + { efSTO, "-avb", "averageb.pdb", ffWRITE }, + { efDAT, "-ascii","correlation", ffOPTWR }, + { efXPM, "-xpm","correlation", ffOPTWR } + }; +#define NFILE asize(fnm) + + CopyRight(stderr,argv[0]); + parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, + NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv); + + clear_mat(zerobox); + + fitfile = ftp2fn(efTPS,NFILE,fnm); + trxfile = ftp2fn(efTRX,NFILE,fnm); + ndxfile = ftp2fn_null(efNDX,NFILE,fnm); + averfilea = ftp2fn(efSTO,NFILE,fnm); + averfileb = ftp2fn(efSTO,NFILE,fnm); + asciifile = opt2fn_null("-ascii",NFILE,fnm); + xpmfile = opt2fn_null("-xpm",NFILE,fnm); + + read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE); + atoms=&top.atoms; + + if (bFit) { + printf("\nChoose a group for the least squares fit\n"); + get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname); + if (nfit < 3) + gmx_fatal(FARGS,"Need >= 3 points to fit!\n"); + } else + nfit=0; + printf("\nChoose two groups for the correlation analysis\n"); + get_index(atoms,ndxfile,1,&natomsa,&indexa,&ananamea); + get_index(atoms,ndxfile,1,&natomsb,&indexb,&ananameb); + + if (bFit) { + snew(w_rls,atoms->nr); + for(i=0; (inr,box,xref,xref); + if (bFit) + reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls); + + snew(xa,natomsa); + snew(xb,natomsb); + snew(xava,natomsa); + snew(xavb,natomsb); +/* if (sqrt(LARGE_INT_MAX)nr) + fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat); + do { + nframes0++; + /* calculate x: a fitted struture of the selected atoms */ + if (bPBC) + rm_pbc(&(top.idef),ePBC,nat,box,xread,xread); + if (bFit) { + reset_x(nfit,ifit,nat,NULL,xread,w_rls); + do_fit(nat,w_rls,xref,xread); + } + for (i=0; i max) + max = mat[i][j]; + } + } + snew(axisa,natomsa); + snew(axisb,natomsb); + for(i=0;i