include(ThreadMPI)
endif()
- # we need this linker flag in case if we have ld >= 2.22 (typically with gcc 4.5+),
- # but it's too cumbersome to check the ld version and the flag should not hurt
- if(CMAKE_COMPILER_IS_GNUCC)
- set(GROMACS_LINKER_FLAGS "-Wl,--add-needed ${GROMACS_LINKER_FLAGS}")
- endif()
endif()
# Annoyingly enough, FindCUDA leaves a few variables behind as non-advanced.
# We need to mark these advanced outside the conditional, otherwise, if the user
extern "C" {
#endif
+int
+gmx_anadock(int argc,char *argv[]);
+
int
gmx_analyze(int argc,char *argv[]);
int
gmx_dist(int argc,char *argv[]);
+int
+gmx_do_dssp(int argc,char *argv[]);
+
int
gmx_dos(int argc,char *argv[]);
int
gmx_kinetics(int argc,char *argv[]);
+int
+gmx_make_edi(int argc,char *argv[]);
+
+int
+gmx_make_ndx(int argc,char *argv[]);
+
int
gmx_mindist(int argc,char *argv[]);
+int
+gmx_mk_angndx(int argc,char *argv[]);
+
int
gmx_msd(int argc,char *argv[]);
int
gmx_sham(int argc,char *argv[]);
+int
+gmx_sigeps(int argc,char *argv[]);
+
int
gmx_sorient(int argc,char *argv[]);
set(GMX_OPENMM_LIBRARIES openmm_api_wrapper ${OpenMM_LIBRARIES})
endif(GMX_OPENMM)
+if(GMX_GPU OR GMX_OPENMM OR GMX_FORCE_CXX)
+ set_source_files_properties(main.c PROPERTIES LANGUAGE CXX)
+ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+ set_source_files_properties(main.c PROPERTIES COMPILE_FLAGS "-x c++")
+ endif()
+endif()
+
if(GMX_FAHCORE)
add_library(fahcore ${MDRUN_SOURCES})
else(GMX_FAHCORE)
grompp tpbconv pdb2gmx g_protonate g_luck gmxdump g_x2top gmxcheck)
foreach(PROGRAM ${GMX_KERNEL_PROGRAMS})
- add_executable(${PROGRAM} ${PROGRAM}.c)
+ add_executable(${PROGRAM} ${PROGRAM}.c main.c)
if (NOT ${PROGRAM} STREQUAL "g_luck")
gmx_add_man_page(${PROGRAM})
endif()
set_target_properties(${PROGRAM} PROPERTIES OUTPUT_NAME "${PROGRAM}${GMX_BINARY_SUFFIX}")
endforeach()
-add_executable(mdrun ${MDRUN_SOURCES})
+add_executable(mdrun ${MDRUN_SOURCES} main.c)
gmx_add_man_page(mdrun)
target_link_libraries(mdrun ${GMX_EXTRA_LIBRARIES} ${GMX_OPENMM_LIBRARIES})
set_target_properties(mdrun PROPERTIES OUTPUT_NAME "mdrun${GMX_BINARY_SUFFIX}" COMPILE_FLAGS "${OpenMP_C_FLAGS}")
#include <string.h>
#include "statutil.h"
-int main(int argc,char *argv[])
+int cmain(int argc,char *argv[])
{
char quote[256];
int i,bc=0,bb=0;
#include "vec.h"
#include "hackblock.h"
-int main (int argc,char *argv[])
+int cmain (int argc,char *argv[])
{
const char *desc[] = {
"[TT]g_protonate[tt] reads (a) conformation(s) and adds all missing",
gmx_fio_fclose(fp);
}
-int main(int argc, char *argv[])
+int cmain(int argc, char *argv[])
{
const char *desc[] = {
"[TT]g_x2top[tt] generates a primitive topology from a coordinate file.",
sfree(fr);
}
-int main(int argc,char *argv[])
+int cmain(int argc,char *argv[])
{
const char *desc[] = {
"[TT]gmxcheck[tt] reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
sfree(full);
}
-int main(int argc,char *argv[])
+int cmain(int argc,char *argv[])
{
const char *desc[] = {
"[TT]gmxdump[tt] reads a run input file ([TT].tpa[tt]/[TT].tpr[tt]/[TT].tpb[tt]),",
}
}
-int main (int argc, char *argv[])
+int cmain (int argc, char *argv[])
{
static const char *desc[] = {
"The gromacs preprocessor",
--- /dev/null
+/*
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * Copyright (c) 2001-2012, The GROMACS Development Team
+ *
+ * Gromacs is a library for molecular simulation and trajectory analysis,
+ * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
+ * a full list of developers and information, check out http://www.gromacs.org
+ *
+ * This program is free software; you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation; either version 2 of the License, or (at your option) any
+ * later version.
+ *
+ * To help fund GROMACS development, we humbly ask that you cite
+ * the papers people have written on it - you can find them on the website.
+ */
+#ifdef __cplusplus
+extern "C" {
+#endif
+int cmain(int argc,char *argv[]);
+#ifdef __cplusplus
+}
+#endif
+int main(int argc,char *argv[])
+{
+ return cmain(argc, argv);
+}
/* afm stuf */
#include "pull.h"
-int main(int argc,char *argv[])
+int cmain(int argc,char *argv[])
{
const char *desc[] = {
#ifdef GMX_OPENMM
rvec *x;
} t_chain;
-int main(int argc, char *argv[])
+int cmain(int argc, char *argv[])
{
const char *desc[] = {
"This program reads a [TT].pdb[tt] (or [TT].gro[tt]) file, reads",
}
}
-int main (int argc, char *argv[])
+int cmain (int argc, char *argv[])
{
const char *desc[] = {
"tpbconv can edit run input files in four ways.[PAR]",
addconf.c calcpot.c edittop.c gmx_bar.c
gmx_membed.c gmx_pme_error.c gmx_options.c gmx_dos.c
gmx_hydorder.c gmx_densorder.c powerspect.c dens_filter.c
- binsearch.c gmx_dyecoupl.c
+ binsearch.c gmx_dyecoupl.c gmx_make_edi.c gmx_sigeps.c
+ gmx_do_dssp.c gmx_anadock.c gmx_make_ndx.c gmx_mk_angndx.c
)
foreach(TOOL ${GMX_TOOLS_PROGRAMS} ${GMX_TOOLS_PROGRAMS_NOT_FOR_INSTALLATION})
+ if(GMX_GPU OR GMX_OPENMM OR GMX_FORCE_CXX)
+ set_source_files_properties(${TOOL}.c PROPERTIES LANGUAGE CXX)
+ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+ set_source_files_properties(${TOOL}.c PROPERTIES COMPILE_FLAGS "-x c++")
+ endif()
+ endif()
add_executable(${TOOL} ${TOOL}.c)
target_link_libraries(${TOOL} gmxana ${OpenMP_LINKER_FLAGS})
set_target_properties(${TOOL} PROPERTIES OUTPUT_NAME "${TOOL}${GMX_BINARY_SUFFIX}")
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
- *
+/*
*
* This source code is part of
*
#include <config.h>
#endif
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "string2.h"
-#include "strdb.h"
-#include "macros.h"
-#include "smalloc.h"
-#include "mshift.h"
-#include "statutil.h"
-#include "copyrite.h"
-#include "pdbio.h"
-#include "gmx_fatal.h"
-#include "xvgr.h"
-#include "matio.h"
-#include "index.h"
-#include "gstat.h"
-#include "tpxio.h"
-#include "viewit.h"
+#include <gmx_ana.h>
-static int strip_dssp(char *dsspfile,int nres,
- gmx_bool bPhobres[],real t,
- real *acc,FILE *fTArea,
- t_matrix *mat,int average_area[],
- const output_env_t oenv)
+/* This is just a wrapper binary.
+*/
+int
+main(int argc,char *argv[])
{
- static gmx_bool bFirst=TRUE;
- static char *ssbuf;
- FILE *tapeout;
- static int xsize,frame;
- char buf[STRLEN+1];
- char SSTP;
- int i,nr,iacc,nresidues;
- int naccf,naccb; /* Count hydrophobic and hydrophilic residues */
- real iaccf,iaccb;
- t_xpmelmt c;
-
- tapeout=ffopen(dsspfile,"r");
-
- /* Skip header */
- do {
- fgets2(buf,STRLEN,tapeout);
- } while (strstr(buf,"KAPPA") == NULL);
- if (bFirst) {
- /* Since we also have empty lines in the dssp output (temp) file,
- * and some line content is saved to the ssbuf variable,
- * we need more memory than just nres elements. To be shure,
- * we allocate 2*nres-1, since for each chain there is a
- * separating line in the temp file. (At most each residue
- * could have been defined as a separate chain.) */
- snew(ssbuf,2*nres-1);
- }
-
- iaccb=iaccf=0;
- nresidues = 0;
- naccf = 0;
- naccb = 0;
- for(nr=0; (fgets2(buf,STRLEN,tapeout) != NULL); nr++) {
- if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
- SSTP='='; /* Chain separator sign '=' */
- else
- SSTP=buf[16]==' ' ? '~' : buf[16];
- ssbuf[nr]=SSTP;
-
- buf[39]='\0';
-
- /* Only calculate solvent accessible area if needed */
- if ((NULL != acc) && (buf[13] != '!'))
- {
- sscanf(&(buf[34]),"%d",&iacc);
- acc[nr]=iacc;
- /* average_area and bPhobres are counted from 0...nres-1 */
- average_area[nresidues]+=iacc;
- if (bPhobres[nresidues])
- {
- naccb++;
- iaccb+=iacc;
- }
- else
- {
- naccf++;
- iaccf+=iacc;
- }
- /* Keep track of the residue number (do not count chain separator lines '!') */
- nresidues++;
- }
- }
- ssbuf[nr]='\0';
-
- if (bFirst) {
- if (0 != acc)
- fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
-
- sprintf(mat->title,"Secondary structure");
- mat->legend[0]=0;
- sprintf(mat->label_x,"%s",output_env_get_time_label(oenv));
- sprintf(mat->label_y,"Residue");
- mat->bDiscrete=TRUE;
- mat->ny=nr;
- snew(mat->axis_y,nr);
- for(i=0; i<nr; i++)
- mat->axis_y[i]=i+1;
- mat->axis_x=NULL;
- mat->matrix=NULL;
- xsize=0;
- frame=0;
- bFirst=FALSE;
- }
- if (frame>=xsize) {
- xsize+=10;
- srenew(mat->axis_x,xsize);
- srenew(mat->matrix,xsize);
- }
- mat->axis_x[frame]=t;
- snew(mat->matrix[frame],nr);
- c.c2=0;
- for(i=0; i<nr; i++) {
- c.c1=ssbuf[i];
- mat->matrix[frame][i] = max(0,searchcmap(mat->nmap,mat->map,c));
- }
- frame++;
- mat->nx=frame;
-
- if (fTArea)
- fprintf(fTArea,"%10g %10g %10g\n",t,0.01*iaccb,0.01*iaccf);
- ffclose(tapeout);
-
- /* Return the number of lines found in the dssp file (i.e. number
- * of redidues plus chain separator lines).
- * This is the number of y elements needed for the area xpm file */
- return nr;
-}
-
-static gmx_bool *bPhobics(t_atoms *atoms)
-{
- int i,nb;
- char **cb;
- gmx_bool *bb;
-
-
- nb = get_strings("phbres.dat",&cb);
- snew(bb,atoms->nres);
-
- for (i=0; (i<atoms->nres); i++)
- {
- if ( -1 != search_str(nb,cb,*atoms->resinfo[i].name) )
- {
- bb[i]=TRUE;
- }
- }
- return bb;
-}
-
-static void check_oo(t_atoms *atoms)
-{
- char *OOO;
-
- int i;
-
- OOO=strdup("O");
-
- for(i=0; (i<atoms->nr); i++) {
- if (strcmp(*(atoms->atomname[i]),"OXT")==0)
- *atoms->atomname[i]=OOO;
- else if (strcmp(*(atoms->atomname[i]),"O1")==0)
- *atoms->atomname[i]=OOO;
- else if (strcmp(*(atoms->atomname[i]),"OC1")==0)
- *atoms->atomname[i]=OOO;
- }
-}
-
-static void norm_acc(t_atoms *atoms, int nres,
- real av_area[], real norm_av_area[])
-{
- int i,n,n_surf;
-
- char surffn[]="surface.dat";
- char **surf_res, **surf_lines;
- double *surf;
-
- n_surf = get_lines(surffn, &surf_lines);
- snew(surf, n_surf);
- snew(surf_res, n_surf);
- for(i=0; (i<n_surf); i++) {
- snew(surf_res[i], 5);
- sscanf(surf_lines[i],"%s %lf",surf_res[i],&surf[i]);
- }
-
- for(i=0; (i<nres); i++) {
- n = search_str(n_surf,surf_res,*atoms->resinfo[i].name);
- if ( n != -1)
- norm_av_area[i] = av_area[i] / surf[n];
- else
- fprintf(stderr,"Residue %s not found in surface database (%s)\n",
- *atoms->resinfo[i].name,surffn);
- }
-}
-
-void prune_ss_legend(t_matrix *mat)
-{
- gmx_bool *present;
- int *newnum;
- int i,r,f,newnmap;
- t_mapping *newmap;
-
- snew(present,mat->nmap);
- snew(newnum,mat->nmap);
-
- for(f=0; f<mat->nx; f++)
- for(r=0; r<mat->ny; r++)
- present[mat->matrix[f][r]]=TRUE;
-
- newnmap=0;
- newmap=NULL;
- for (i=0; i<mat->nmap; i++) {
- newnum[i]=-1;
- if (present[i]) {
- newnum[i]=newnmap;
- newnmap++;
- srenew(newmap,newnmap);
- newmap[newnmap-1]=mat->map[i];
- }
- }
- if (newnmap!=mat->nmap) {
- mat->nmap=newnmap;
- mat->map=newmap;
- for(f=0; f<mat->nx; f++)
- for(r=0; r<mat->ny; r++)
- mat->matrix[f][r]=newnum[mat->matrix[f][r]];
- }
-}
-
-void write_sas_mat(const char *fn,real **accr,int nframe,int nres,t_matrix *mat)
-{
- real lo,hi;
- int i,j,nlev;
- t_rgb rlo={1,1,1}, rhi={0,0,0};
- FILE *fp;
-
- if(fn) {
- hi=lo=accr[0][0];
- for(i=0; i<nframe; i++)
- for(j=0; j<nres; j++) {
- lo=min(lo,accr[i][j]);
- hi=max(hi,accr[i][j]);
- }
- fp=ffopen(fn,"w");
- nlev=hi-lo+1;
- write_xpm(fp,0,"Solvent Accessible Surface","Surface (A^2)",
- "Time","Residue Index",nframe,nres,
- mat->axis_x,mat->axis_y,accr,lo,hi,rlo,rhi,&nlev);
- ffclose(fp);
- }
-}
-
-void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
- const output_env_t oenv)
-{
- FILE *fp;
- t_mapping *map;
- int f,r,*count,*total,ss_count,total_count;
- size_t s;
- const char** leg;
-
- map=mat->map;
- snew(count,mat->nmap);
- snew(total,mat->nmap);
- snew(leg,mat->nmap+1);
- leg[0]="Structure";
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- leg[s+1]=strdup(map[s].desc);
- }
-
- fp=xvgropen(outfile,"Secondary Structure",
- output_env_get_xvgr_tlabel(oenv),"Number of Residues",oenv);
- if (output_env_get_print_xvgr_codes(oenv))
- {
- fprintf(fp,"@ subtitle \"Structure = ");
- }
- for(s=0; s<strlen(ss_string); s++)
- {
- if (s>0)
- {
- fprintf(fp," + ");
- }
- for(f=0; f<mat->nmap; f++)
- {
- if (ss_string[s]==map[f].code.c1)
- {
- fprintf(fp,"%s",map[f].desc);
- }
- }
- }
- fprintf(fp,"\"\n");
- xvgr_legend(fp,mat->nmap+1,leg,oenv);
-
- total_count = 0;
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- total[s]=0;
- }
- for(f=0; f<mat->nx; f++)
- {
- ss_count=0;
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- count[s]=0;
- }
- for(r=0; r<mat->ny; r++)
- {
- count[mat->matrix[f][r]]++;
- total[mat->matrix[f][r]]++;
- }
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- if (strchr(ss_string,map[s].code.c1))
- {
- ss_count+=count[s];
- total_count += count[s];
- }
- }
- fprintf(fp,"%8g %5d",mat->axis_x[f],ss_count);
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- fprintf(fp," %5d",count[s]);
- }
- fprintf(fp,"\n");
- }
- /* now print column totals */
- fprintf(fp, "%-8s %5d", "# Totals", total_count);
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- fprintf(fp," %5d",total[s]);
- }
- fprintf(fp,"\n");
-
- /* now print percentages */
- fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));
- for(s=0; s<(size_t)mat->nmap; s++)
- {
- fprintf(fp," %5.2f",total[s] / (real) (mat->nx * mat->ny));
- }
- fprintf(fp,"\n");
-
- ffclose(fp);
- sfree(leg);
- sfree(count);
-}
-
-int main(int argc,char *argv[])
-{
- const char *desc[] = {
- "[TT]do_dssp[tt] ",
- "reads a trajectory file and computes the secondary structure for",
- "each time frame ",
- "calling the dssp program. If you do not have the dssp program,",
- "get it from http://swift.cmbi.ru.nl/gv/dssp. [TT]do_dssp[tt] assumes ",
- "that the dssp executable is located in ",
- "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
- "set an environment variable [TT]DSSP[tt] pointing to the dssp",
- "executable, e.g.: [PAR]",
- "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
- "Since version 2.0.0, dssp is invoked with a syntax that differs",
- "from earlier versions. If you have an older version of dssp,",
- "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
- "By default, do_dssp uses the syntax introduced with version 2.0.0.",
- "Even newer versions (which at the time of writing are not yet released)",
- "are assumed to have the same syntax as 2.0.0.[PAR]",
- "The structure assignment for each residue and time is written to an",
- "[TT].xpm[tt] matrix file. This file can be visualized with for instance",
- "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
- "Individual chains are separated by light grey lines in the [TT].xpm[tt] and",
- "postscript files.",
- "The number of residues with each secondary structure type and the",
- "total secondary structure ([TT]-sss[tt]) count as a function of",
- "time are also written to file ([TT]-sc[tt]).[PAR]",
- "Solvent accessible surface (SAS) per residue can be calculated, both in",
- "absolute values (A^2) and in fractions of the maximal accessible",
- "surface of a residue. The maximal accessible surface is defined as",
- "the accessible surface of a residue in a chain of glycines.",
- "[BB]Note[bb] that the program [TT]g_sas[tt] can also compute SAS",
- "and that is more efficient.[PAR]",
- "Finally, this program can dump the secondary structure in a special file",
- "[TT]ssdump.dat[tt] for usage in the program [TT]g_chi[tt]. Together",
- "these two programs can be used to analyze dihedral properties as a",
- "function of secondary structure type."
- };
- static gmx_bool bVerbose;
- static const char *ss_string="HEBT";
- static int dsspVersion=2;
- t_pargs pa[] = {
- { "-v", FALSE, etBOOL, {&bVerbose},
- "HIDDENGenerate miles of useless information" },
- { "-sss", FALSE, etSTR, {&ss_string},
- "Secondary structures for structure count"},
- { "-ver", FALSE, etINT, {&dsspVersion},
- "DSSP major version. Syntax changed with version 2"}
- };
-
- t_trxstatus *status;
- FILE *tapein;
- FILE *ss,*acc,*fTArea,*tmpf;
- const char *fnSCount,*fnArea,*fnTArea,*fnAArea;
- const char *leg[] = { "Phobic", "Phylic" };
- t_topology top;
- int ePBC;
- t_atoms *atoms;
- t_matrix mat;
- int nres,nr0,naccr,nres_plus_separators;
- gmx_bool *bPhbres,bDoAccSurf;
- real t;
- int i,j,natoms,nframe=0;
- matrix box;
- int gnx;
- char *grpnm,*ss_str;
- atom_id *index;
- rvec *xp,*x;
- int *average_area;
- real **accr,*accr_ptr=NULL,*av_area, *norm_av_area;
- char pdbfile[32],tmpfile[32],title[256];
- char dssp[256];
- const char *dptr;
- output_env_t oenv;
- gmx_rmpbc_t gpbc=NULL;
-
- t_filenm fnm[] = {
- { efTRX, "-f", NULL, ffREAD },
- { efTPS, NULL, NULL, ffREAD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efDAT, "-ssdump", "ssdump", ffOPTWR },
- { efMAP, "-map", "ss", ffLIBRD },
- { efXPM, "-o", "ss", ffWRITE },
- { efXVG, "-sc", "scount", ffWRITE },
- { efXPM, "-a", "area", ffOPTWR },
- { efXVG, "-ta", "totarea", ffOPTWR },
- { efXVG, "-aa", "averarea",ffOPTWR }
- };
-#define NFILE asize(fnm)
-
- CopyRight(stderr,argv[0]);
- parse_common_args(&argc,argv,
- PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
- NFILE,fnm, asize(pa),pa, asize(desc),desc,0,NULL,&oenv);
- fnSCount= opt2fn("-sc",NFILE,fnm);
- fnArea = opt2fn_null("-a", NFILE,fnm);
- fnTArea = opt2fn_null("-ta",NFILE,fnm);
- fnAArea = opt2fn_null("-aa",NFILE,fnm);
- bDoAccSurf=(fnArea || fnTArea || fnAArea);
-
- read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
- atoms=&(top.atoms);
- check_oo(atoms);
- bPhbres = bPhobics(atoms);
-
- get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
- nres=0;
- nr0=-1;
- for(i=0; (i<gnx); i++) {
- if (atoms->atom[index[i]].resind != nr0) {
- nr0=atoms->atom[index[i]].resind;
- nres++;
- }
- }
- fprintf(stderr,"There are %d residues in your selected group\n",nres);
-
- strcpy(pdbfile,"ddXXXXXX");
- gmx_tmpnam(pdbfile);
- if ((tmpf = fopen(pdbfile,"w")) == NULL) {
- sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
- gmx_tmpnam(pdbfile);
- if ((tmpf = fopen(pdbfile,"w")) == NULL)
- gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);
- }
- else
- fclose(tmpf);
-
- strcpy(tmpfile,"ddXXXXXX");
- gmx_tmpnam(tmpfile);
- if ((tmpf = fopen(tmpfile,"w")) == NULL) {
- sprintf(tmpfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
- gmx_tmpnam(tmpfile);
- if ((tmpf = fopen(tmpfile,"w")) == NULL)
- gmx_fatal(FARGS,"Can not open tmp file %s",tmpfile);
- }
- else
- fclose(tmpf);
-
- if ((dptr=getenv("DSSP")) == NULL)
- dptr="/usr/local/bin/dssp";
- if (!gmx_fexist(dptr))
- gmx_fatal(FARGS,"DSSP executable (%s) does not exist (use setenv DSSP)",
- dptr);
- if (dsspVersion >= 2)
- {
- if (dsspVersion > 2)
- {
- printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
- }
-
- sprintf(dssp,"%s -i %s -o %s > /dev/null %s",
- dptr,pdbfile,tmpfile,bVerbose?"":"2> /dev/null");
- }
- else
- {
- sprintf(dssp,"%s %s %s %s > /dev/null %s",
- dptr,bDoAccSurf?"":"-na",pdbfile,tmpfile,bVerbose?"":"2> /dev/null");
-
- }
- fprintf(stderr,"dssp cmd='%s'\n",dssp);
-
- if (fnTArea) {
- fTArea=xvgropen(fnTArea,"Solvent Accessible Surface Area",
- output_env_get_xvgr_tlabel(oenv),"Area (nm\\S2\\N)",oenv);
- xvgr_legend(fTArea,2,leg,oenv);
- } else
- fTArea=NULL;
-
- mat.map=NULL;
- mat.nmap=getcmap(libopen(opt2fn("-map",NFILE,fnm)),
- opt2fn("-map",NFILE,fnm),&(mat.map));
-
- natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
- if (natoms > atoms->nr)
- gmx_fatal(FARGS,"\nTrajectory does not match topology!");
- if (gnx > natoms)
- gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
-
- snew(average_area, atoms->nres);
- snew(av_area , atoms->nres);
- snew(norm_av_area, atoms->nres);
- accr=NULL;
- naccr=0;
-
- gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
- do {
- t = output_env_conv_time(oenv,t);
- if (bDoAccSurf && nframe>=naccr) {
- naccr+=10;
- srenew(accr,naccr);
- for(i=naccr-10; i<naccr; i++)
- snew(accr[i],2*atoms->nres-1);
- }
- gmx_rmpbc(gpbc,natoms,box,x);
- tapein=ffopen(pdbfile,"w");
- write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,' ',-1,gnx,index,NULL,TRUE);
- ffclose(tapein);
-
-#ifdef GMX_NO_SYSTEM
- printf("Warning-- No calls to system(3) supported on this platform.");
- printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
- exit(1);
-#else
- if(0 != system(dssp))
- {
- gmx_fatal(FARGS,"Failed to execute command: %s\n",
- "Try specifying your dssp version with the -ver option.",dssp);
- }
-#endif
-
- /* strip_dssp returns the number of lines found in the dssp file, i.e.
- * the number of residues plus the separator lines */
-
- if (bDoAccSurf)
- accr_ptr = accr[nframe];
-
- nres_plus_separators = strip_dssp(tmpfile,nres,bPhbres,t,
- accr_ptr,fTArea,&mat,average_area,oenv);
- remove(tmpfile);
- remove(pdbfile);
- nframe++;
- } while(read_next_x(oenv,status,&t,natoms,x,box));
- fprintf(stderr,"\n");
- close_trj(status);
- if (fTArea)
- ffclose(fTArea);
- gmx_rmpbc_done(gpbc);
-
- prune_ss_legend(&mat);
-
- ss=opt2FILE("-o",NFILE,fnm,"w");
- mat.flags = 0;
- write_xpm_m(ss,mat);
- ffclose(ss);
-
- if (opt2bSet("-ssdump",NFILE,fnm))
- {
- ss = opt2FILE("-ssdump",NFILE,fnm,"w");
- snew(ss_str,nres+1);
- fprintf(ss,"%d\n",nres);
- for(j=0; j<mat.nx; j++)
- {
- for(i=0; (i<mat.ny); i++)
- {
- ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
- }
- ss_str[i] = '\0';
- fprintf(ss,"%s\n",ss_str);
- }
- ffclose(ss);
- sfree(ss_str);
- }
- analyse_ss(fnSCount,&mat,ss_string,oenv);
-
- if (bDoAccSurf) {
- write_sas_mat(fnArea,accr,nframe,nres_plus_separators,&mat);
-
- for(i=0; i<atoms->nres; i++)
- av_area[i] = (average_area[i] / (real)nframe);
-
- norm_acc(atoms, nres, av_area, norm_av_area);
-
- if (fnAArea) {
- acc=xvgropen(fnAArea,"Average Accessible Area",
- "Residue","A\\S2",oenv);
- for(i=0; (i<nres); i++)
- fprintf(acc,"%5d %10g %10g\n",i+1,av_area[i], norm_av_area[i]);
- ffclose(acc);
- }
- }
-
- view_all(oenv, NFILE, fnm);
-
- thanx(stderr);
-
- return 0;
+ return gmx_do_dssp(argc,argv);
}
#include <config.h>
#endif
-#include "confio.h"
-#include "copyrite.h"
-#include "futil.h"
-#include "gmx_fatal.h"
-#include "smalloc.h"
-#include "string2.h"
-#include "vec.h"
-#include "gmx_statistics.h"
-#include "statutil.h"
-#include "typedefs.h"
-#include "xvgr.h"
-
-static const char *etitles[] = { "E-docked", "Free Energy" };
+#include <gmx_ana.h>
-typedef struct {
- real edocked,efree;
- int index,cluster_id;
- t_atoms atoms;
- rvec *x;
- int ePBC;
- matrix box;
-} t_pdbfile;
-static t_pdbfile *read_pdbf(const char *fn)
+/* This is just a wrapper binary.
+*/
+int
+main(int argc,char *argv[])
{
- t_pdbfile *pdbf;
- double e;
- char buf[256],*ptr;
- int natoms;
- FILE *fp;
-
- snew(pdbf,1);
- get_stx_coordnum (fn,&natoms);
- init_t_atoms(&(pdbf->atoms),natoms,FALSE);
- snew(pdbf->x,natoms);
- read_stx_conf(fn,buf,&pdbf->atoms,pdbf->x,NULL,&pdbf->ePBC,pdbf->box);
- fp = ffopen(fn,"r");
- do {
- ptr = fgets2(buf,255,fp);
- if (ptr) {
- if (strstr(buf,"Intermolecular") != NULL) {
- ptr = strchr(buf,'=');
- sscanf(ptr+1,"%lf",&e);
- pdbf->edocked = e;
- }
- else if (strstr(buf,"Estimated Free") != NULL) {
- ptr = strchr(buf,'=');
- sscanf(ptr+1,"%lf",&e);
- pdbf->efree = e;
- }
- }
- } while (ptr != NULL);
- ffclose(fp);
-
- return pdbf;
-}
-
-static t_pdbfile **read_em_all(const char *fn,int *npdbf)
-{
- t_pdbfile **pdbf=0;
- int i,maxpdbf;
- char buf[256],name[256];
- gmx_bool bExist;
-
- strcpy(buf,fn);
- buf[strlen(buf)-4] = '\0';
- maxpdbf = 100;
- snew(pdbf,maxpdbf);
- i=0;
- do {
- sprintf(name,"%s_%d.pdb",buf,i+1);
- if ((bExist = gmx_fexist(name)) == TRUE) {
- pdbf[i] = read_pdbf(name);
- pdbf[i]->index = i+1;
- i++;
- if (i >= maxpdbf) {
- maxpdbf += 100;
- srenew(pdbf,maxpdbf);
- }
- }
- } while (bExist);
-
- *npdbf = i;
-
- printf("Found %d pdb files\n",i);
-
- return pdbf;
-}
-
-static gmx_bool bFreeSort=FALSE;
-
-static int pdbf_comp(const void *a,const void *b)
-{
- t_pdbfile *pa,*pb;
- real x;
- int dc;
-
- pa = *(t_pdbfile **)a;
- pb = *(t_pdbfile **)b;
-
- dc = pa->cluster_id - pb->cluster_id;
-
- if (dc == 0) {
- if (bFreeSort)
- x = pa->efree - pb->efree;
- else
- x = pa->edocked - pb->edocked;
-
- if (x < 0)
- return -1;
- else if (x > 0)
- return 1;
- else
- return 0;
- }
- else
- return dc;
-}
-
-static void analyse_em_all(int npdb,t_pdbfile *pdbf[], const char *edocked,
- const char *efree, const output_env_t oenv)
-{
- FILE *fp;
- int i;
-
- for(bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++) {
- qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
- fp = xvgropen(bFreeSort ? efree : edocked,
- etitles[bFreeSort],"()","E (kJ/mol)",oenv);
- for(i=0; (i<npdb); i++)
- fprintf(fp,"%12lf\n",bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
- ffclose(fp);
- }
-}
-
-static void clust_stat(FILE *fp,int start,int end,t_pdbfile *pdbf[])
-{
- int i;
- gmx_stats_t ed,ef;
- real aver,sigma;
-
- ed = gmx_stats_init();
- ef = gmx_stats_init();
- for(i=start; (i<end); i++) {
- gmx_stats_add_point(ed,i-start,pdbf[i]->edocked,0,0);
- gmx_stats_add_point(ef,i-start,pdbf[i]->efree,0,0);
- }
- gmx_stats_get_ase(ed,&aver,&sigma,NULL);
- fprintf(fp," <%12s> = %8.3f (+/- %6.3f)\n",etitles[FALSE],aver,sigma);
- gmx_stats_get_ase(ef,&aver,&sigma,NULL);
- fprintf(fp," <%12s> = %8.3f (+/- %6.3f)\n",etitles[TRUE],aver,sigma);
- gmx_stats_done(ed);
- gmx_stats_done(ef);
- sfree(ed);
- sfree(ef);
-}
-
-static real rmsd_dist(t_pdbfile *pa,t_pdbfile *pb,gmx_bool bRMSD)
-{
- int i;
- real rmsd,dist;
- rvec acm,bcm,dx;
-
- if (bRMSD) {
- rmsd = 0;
- for(i=0; (i<pa->atoms.nr); i++) {
- rvec_sub(pa->x[i],pb->x[i],dx);
- rmsd += iprod(dx,dx);
- }
- rmsd = sqrt(rmsd/pa->atoms.nr);
- }
- else {
- dist = 0;
- clear_rvec(acm);
- clear_rvec(bcm);
- for(i=0; (i<pa->atoms.nr); i++) {
- rvec_inc(acm,pa->x[i]);
- rvec_inc(bcm,pb->x[i]);
- }
- rvec_sub(acm,bcm,dx);
- for(i=0; (i<DIM); i++)
- dx[i] /= pa->atoms.nr;
- rmsd = norm(dx);
- }
- return rmsd;
-}
-
-static void line(FILE *fp)
-{
- fprintf(fp," - - - - - - - - - -\n");
-}
-
-static void cluster_em_all(FILE *fp,int npdb,t_pdbfile *pdbf[],
- const char *pdbout,gmx_bool bFree,gmx_bool bRMSD,real cutoff)
-{
- int i,j,k;
- int *cndx,ncluster;
- real rmsd;
-
- bFreeSort = bFree;
- qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
-
- fprintf(fp,"Statistics over all structures:\n");
- clust_stat(fp,0,npdb,pdbf);
- line(fp);
-
- /* Index to first structure in a cluster */
- snew(cndx,npdb);
- ncluster=1;
-
- for(i=1; (i<npdb); i++) {
- for(j=0; (j<ncluster); j++) {
- rmsd = rmsd_dist(pdbf[cndx[j]],pdbf[i],bRMSD);
- if (rmsd <= cutoff) {
- /* Structure i is in cluster j */
- pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
- break;
- }
- }
- if (j == ncluster) {
- /* New cluster! Cool! */
- cndx[ncluster] = i;
- pdbf[i]->cluster_id = ncluster;
- ncluster++;
- }
- }
- cndx[ncluster] = npdb;
- qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
-
- j=0;
- cndx[j++] = 0;
- for(i=1; (i<npdb); i++) {
- if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id) {
- cndx[j++] = i;
- }
- }
- cndx[j] = npdb;
- if (j != ncluster)
- gmx_fatal(FARGS,"Consistency error: j = %d, ncluster = %d",j,ncluster);
-
- fprintf(fp,"I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
- ncluster,etitles[bFree],bRMSD ? "RMSD" : "distance",cutoff);
- line(fp);
- for(j=0; (j<ncluster); j++) {
- fprintf(fp,"Cluster: %3d %s: %10.5f kJ/mol %3d elements\n",
- j,etitles[bFree],
- bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
- cndx[j+1]-cndx[j]);
- clust_stat(fp,cndx[j],cndx[j+1],pdbf);
- for(k=cndx[j]; (k<cndx[j+1]); k++)
- fprintf(fp," %3d",pdbf[k]->index);
- fprintf(fp,"\n");
- line(fp);
- }
- sfree(cndx);
-}
-
-int main(int argc,char *argv[])
-{
- const char *desc[] = {
- "[TT]g_anadock[tt] analyses the results of an Autodock run and clusters the",
- "structures together, based on distance or RMSD. The docked energy",
- "and free energy estimates are analysed, and for each cluster the",
- "energy statistics are printed.[PAR]",
- "An alternative approach to this is to cluster the structures first",
- "using [TT]g_cluster[tt] and then sort the clusters on either lowest",
- "energy or average energy."
- };
- t_filenm fnm[] = {
- { efPDB, "-f", NULL, ffREAD },
- { efPDB, "-ox", "cluster", ffWRITE },
- { efXVG, "-od", "edocked", ffWRITE },
- { efXVG, "-of", "efree", ffWRITE },
- { efLOG, "-g", "anadock", ffWRITE }
- };
- output_env_t oenv;
-#define NFILE asize(fnm)
- static gmx_bool bFree=FALSE,bRMS=TRUE;
- static real cutoff = 0.2;
- t_pargs pa[] = {
- { "-free", FALSE, etBOOL, {&bFree},
- "Use Free energy estimate from autodock for sorting the classes" },
- { "-rms", FALSE, etBOOL, {&bRMS},
- "Cluster on RMS or distance" },
- { "-cutoff", FALSE, etREAL, {&cutoff},
- "Maximum RMSD/distance for belonging to the same cluster" }
- };
-#define NPA asize(pa)
-
- FILE *fp;
- t_pdbfile **pdbf=NULL;
- int npdbf;
-
- CopyRight(stderr,argv[0]);
- parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa, asize(desc),desc,0,
- NULL,&oenv);
-
- fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
- please_cite(stdout,"Hetenyi2002b");
- please_cite(fp,"Hetenyi2002b");
-
- pdbf = read_em_all(opt2fn("-f",NFILE,fnm),&npdbf);
-
- analyse_em_all(npdbf,pdbf,opt2fn("-od",NFILE,fnm),opt2fn("-of",NFILE,fnm),
- oenv);
-
- cluster_em_all(fp,npdbf,pdbf,opt2fn("-ox",NFILE,fnm),
- bFree,bRMS,cutoff);
-
- thanx(fp);
- ffclose(fp);
-
- thanx(stdout);
-
- return 0;
+ return gmx_anadock(argc,argv);
}
* For more info, check our website at http://www.gromacs.org
*
* And Hey:
- * Good gRace! Old Maple Actually Chews Slate
+ * Green Red Orange Magenta Azure Cyan Skyblue
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include <stdio.h>
-#include <math.h>
-#include "typedefs.h"
-#include "statutil.h"
-#include "copyrite.h"
-#include "gmx_fatal.h"
-#include "xvgr.h"
-#include "pdbio.h"
-#include "macros.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "pbc.h"
-#include "physics.h"
-#include "names.h"
-#include "txtdump.h"
-#include "trnio.h"
-#include "symtab.h"
-#include "confio.h"
+#include <gmx_ana.h>
-real pot(real x,real qq,real c6,real cn,int npow)
-{
- return cn*pow(x,-npow)-c6*pow(x,-6)+qq*ONE_4PI_EPS0/x;
-}
-real bhpot(real x,real qq,real A,real B,real C)
+/* This is just a wrapper binary.
+*/
+int
+main(int argc,char *argv[])
{
- return A*exp(-B*x) - C*pow(x,-6.0);
+ return gmx_sigeps(argc,argv);
}
-
-real dpot(real x,real qq,real c6,real cn,int npow)
-{
- return -(npow*cn*pow(x,-npow-1)-6*c6*pow(x,-7)+qq*ONE_4PI_EPS0/sqr(x));
-}
-
-int main(int argc,char *argv[])
-{
- const char *desc[] = {
- "[TT]g_sigeps[tt] is a simple utility that converts C6/C12 or C6/Cn combinations",
- "to [GRK]sigma[grk] and [GRK]epsilon[grk], or vice versa. It can also plot the potential",
- "in file. In addition, it makes an approximation of a Buckingham potential",
- "to a Lennard-Jones potential."
- };
- static real c6=1.0e-3,cn=1.0e-6,qi=0,qj=0,sig=0.3,eps=1,sigfac=0.7;
- static real Abh=1e5,Bbh=32,Cbh=1e-3;
- static int npow=12;
- t_pargs pa[] = {
- { "-c6", FALSE, etREAL, {&c6}, "C6" },
- { "-cn", FALSE, etREAL, {&cn}, "Constant for repulsion" },
- { "-pow", FALSE, etINT, {&npow},"Power of the repulsion term" },
- { "-sig", FALSE, etREAL, {&sig}, "[GRK]sigma[grk]" },
- { "-eps", FALSE, etREAL, {&eps}, "[GRK]epsilon[grk]" },
- { "-A", FALSE, etREAL, {&Abh}, "Buckingham A" },
- { "-B", FALSE, etREAL, {&Bbh}, "Buckingham B" },
- { "-C", FALSE, etREAL, {&Cbh}, "Buckingham C" },
- { "-qi", FALSE, etREAL, {&qi}, "qi" },
- { "-qj", FALSE, etREAL, {&qj}, "qj" },
- { "-sigfac", FALSE, etREAL, {&sigfac}, "Factor in front of [GRK]sigma[grk] for starting the plot" }
- };
- t_filenm fnm[] = {
- { efXVG, "-o", "potje", ffWRITE }
- };
- output_env_t oenv;
-#define NFILE asize(fnm)
- const char *legend[] = { "Lennard-Jones", "Buckingham" };
- FILE *fp;
- int i;
- gmx_bool bBham;
- real qq,x,oldx,minimum,mval,dp[2],pp[2];
- int cur=0;
-#define next (1-cur)
-
- /* CopyRight(stdout,argv[0]);*/
- parse_common_args(&argc,argv,PCA_CAN_VIEW,
- NFILE,fnm,asize(pa),pa,asize(desc),
- desc,0,NULL,&oenv);
-
- bBham = (opt2parg_bSet("-A",asize(pa),pa) ||
- opt2parg_bSet("-B",asize(pa),pa) ||
- opt2parg_bSet("-C",asize(pa),pa));
-
- if (bBham) {
- c6 = Cbh;
- sig = pow((6.0/npow)*pow(npow/Bbh,npow-6.0),1.0/(npow-6.0));
- eps = c6/(4*pow(sig,6.0));
- cn = 4*eps*pow(sig,npow);
- }
- else {
- if (opt2parg_bSet("-sig",asize(pa),pa) ||
- opt2parg_bSet("-eps",asize(pa),pa)) {
- c6 = 4*eps*pow(sig,6);
- cn = 4*eps*pow(sig,npow);
- }
- else if (opt2parg_bSet("-c6",asize(pa),pa) ||
- opt2parg_bSet("-cn",asize(pa),pa) ||
- opt2parg_bSet("-pow",asize(pa),pa)) {
- sig = pow(cn/c6,1.0/(npow-6.0));
- eps = 0.25*c6*pow(sig,-6.0);
- }
- else {
- sig = eps = 0;
- }
- printf("c6 = %12.5e, c%d = %12.5e\n",c6,npow,cn);
- printf("sigma = %12.5f, epsilon = %12.5f\n",sig,eps);
-
- minimum = pow(npow/6.0*pow(sig,npow-6.0),1.0/(npow-6));
- printf("Van der Waals minimum at %g, V = %g\n\n",
- minimum,pot(minimum,0,c6,cn,npow));
- printf("Fit of Lennard Jones (%d-6) to Buckingham:\n",npow);
- Bbh = npow/minimum;
- Cbh = c6;
- Abh = 4*eps*pow(sig/minimum,npow)*exp(npow);
- printf("A = %g, B = %g, C = %g\n",Abh,Bbh,Cbh);
- }
- qq = qi*qj;
-
- fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"Potential","r (nm)","E (kJ/mol)",
- oenv);
- xvgr_legend(fp,asize(legend),legend,
- oenv);
- if (sig == 0)
- sig=0.25;
- minimum = -1;
- mval = 0;
- oldx = 0;
- for(i=0; (i<100); i++) {
- x = sigfac*sig+sig*i*0.02;
- dp[next] = dpot(x,qq,c6,cn,npow);
- fprintf(fp,"%10g %10g %10g\n",x,pot(x,qq,c6,cn,npow),
- bhpot(x,qq,Abh,Bbh,Cbh));
- if (qq != 0) {
- if ((i > 0) && (dp[cur]*dp[next] < 0)) {
- minimum = oldx + dp[cur]*(x-oldx)/(dp[cur]-dp[next]);
- mval = pot(minimum,qq,c6,cn,npow);
- printf("Van der Waals + Coulomb minimum at r = %g (nm). Value = %g (kJ/mol)\n",
- minimum,mval);
- }
- }
- cur = next;
- oldx = x;
-
- }
- ffclose(fp);
-
- do_view(oenv,ftp2fn(efXVG,NFILE,fnm),NULL);
-
- thanx(stderr);
-
- return 0;
-}
-
-
--- /dev/null
+/*
+ *
+ * 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 "confio.h"
+#include "copyrite.h"
+#include "futil.h"
+#include "gmx_fatal.h"
+#include "smalloc.h"
+#include "string2.h"
+#include "vec.h"
+#include "gmx_statistics.h"
+#include "statutil.h"
+#include "typedefs.h"
+#include "xvgr.h"
+
+static const char *etitles[] = { "E-docked", "Free Energy" };
+
+typedef struct {
+ real edocked,efree;
+ int index,cluster_id;
+ t_atoms atoms;
+ rvec *x;
+ int ePBC;
+ matrix box;
+} t_pdbfile;
+
+static t_pdbfile *read_pdbf(const char *fn)
+{
+ t_pdbfile *pdbf;
+ double e;
+ char buf[256],*ptr;
+ int natoms;
+ FILE *fp;
+
+ snew(pdbf,1);
+ get_stx_coordnum (fn,&natoms);
+ init_t_atoms(&(pdbf->atoms),natoms,FALSE);
+ snew(pdbf->x,natoms);
+ read_stx_conf(fn,buf,&pdbf->atoms,pdbf->x,NULL,&pdbf->ePBC,pdbf->box);
+ fp = ffopen(fn,"r");
+ do {
+ ptr = fgets2(buf,255,fp);
+ if (ptr) {
+ if (strstr(buf,"Intermolecular") != NULL) {
+ ptr = strchr(buf,'=');
+ sscanf(ptr+1,"%lf",&e);
+ pdbf->edocked = e;
+ }
+ else if (strstr(buf,"Estimated Free") != NULL) {
+ ptr = strchr(buf,'=');
+ sscanf(ptr+1,"%lf",&e);
+ pdbf->efree = e;
+ }
+ }
+ } while (ptr != NULL);
+ ffclose(fp);
+
+ return pdbf;
+}
+
+static t_pdbfile **read_em_all(const char *fn,int *npdbf)
+{
+ t_pdbfile **pdbf=0;
+ int i,maxpdbf;
+ char buf[256],name[256];
+ gmx_bool bExist;
+
+ strcpy(buf,fn);
+ buf[strlen(buf)-4] = '\0';
+ maxpdbf = 100;
+ snew(pdbf,maxpdbf);
+ i=0;
+ do {
+ sprintf(name,"%s_%d.pdb",buf,i+1);
+ if ((bExist = gmx_fexist(name)) == TRUE) {
+ pdbf[i] = read_pdbf(name);
+ pdbf[i]->index = i+1;
+ i++;
+ if (i >= maxpdbf) {
+ maxpdbf += 100;
+ srenew(pdbf,maxpdbf);
+ }
+ }
+ } while (bExist);
+
+ *npdbf = i;
+
+ printf("Found %d pdb files\n",i);
+
+ return pdbf;
+}
+
+static gmx_bool bFreeSort=FALSE;
+
+static int pdbf_comp(const void *a,const void *b)
+{
+ t_pdbfile *pa,*pb;
+ real x;
+ int dc;
+
+ pa = *(t_pdbfile **)a;
+ pb = *(t_pdbfile **)b;
+
+ dc = pa->cluster_id - pb->cluster_id;
+
+ if (dc == 0) {
+ if (bFreeSort)
+ x = pa->efree - pb->efree;
+ else
+ x = pa->edocked - pb->edocked;
+
+ if (x < 0)
+ return -1;
+ else if (x > 0)
+ return 1;
+ else
+ return 0;
+ }
+ else
+ return dc;
+}
+
+static void analyse_em_all(int npdb,t_pdbfile *pdbf[], const char *edocked,
+ const char *efree, const output_env_t oenv)
+{
+ FILE *fp;
+ int i;
+
+ for(bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++) {
+ qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
+ fp = xvgropen(bFreeSort ? efree : edocked,
+ etitles[bFreeSort],"()","E (kJ/mol)",oenv);
+ for(i=0; (i<npdb); i++)
+ fprintf(fp,"%12lf\n",bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
+ ffclose(fp);
+ }
+}
+
+static void clust_stat(FILE *fp,int start,int end,t_pdbfile *pdbf[])
+{
+ int i;
+ gmx_stats_t ed,ef;
+ real aver,sigma;
+
+ ed = gmx_stats_init();
+ ef = gmx_stats_init();
+ for(i=start; (i<end); i++) {
+ gmx_stats_add_point(ed,i-start,pdbf[i]->edocked,0,0);
+ gmx_stats_add_point(ef,i-start,pdbf[i]->efree,0,0);
+ }
+ gmx_stats_get_ase(ed,&aver,&sigma,NULL);
+ fprintf(fp," <%12s> = %8.3f (+/- %6.3f)\n",etitles[FALSE],aver,sigma);
+ gmx_stats_get_ase(ef,&aver,&sigma,NULL);
+ fprintf(fp," <%12s> = %8.3f (+/- %6.3f)\n",etitles[TRUE],aver,sigma);
+ gmx_stats_done(ed);
+ gmx_stats_done(ef);
+ sfree(ed);
+ sfree(ef);
+}
+
+static real rmsd_dist(t_pdbfile *pa,t_pdbfile *pb,gmx_bool bRMSD)
+{
+ int i;
+ real rmsd,dist;
+ rvec acm,bcm,dx;
+
+ if (bRMSD) {
+ rmsd = 0;
+ for(i=0; (i<pa->atoms.nr); i++) {
+ rvec_sub(pa->x[i],pb->x[i],dx);
+ rmsd += iprod(dx,dx);
+ }
+ rmsd = sqrt(rmsd/pa->atoms.nr);
+ }
+ else {
+ dist = 0;
+ clear_rvec(acm);
+ clear_rvec(bcm);
+ for(i=0; (i<pa->atoms.nr); i++) {
+ rvec_inc(acm,pa->x[i]);
+ rvec_inc(bcm,pb->x[i]);
+ }
+ rvec_sub(acm,bcm,dx);
+ for(i=0; (i<DIM); i++)
+ dx[i] /= pa->atoms.nr;
+ rmsd = norm(dx);
+ }
+ return rmsd;
+}
+
+static void line(FILE *fp)
+{
+ fprintf(fp," - - - - - - - - - -\n");
+}
+
+static void cluster_em_all(FILE *fp,int npdb,t_pdbfile *pdbf[],
+ const char *pdbout,gmx_bool bFree,gmx_bool bRMSD,real cutoff)
+{
+ int i,j,k;
+ int *cndx,ncluster;
+ real rmsd;
+
+ bFreeSort = bFree;
+ qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
+
+ fprintf(fp,"Statistics over all structures:\n");
+ clust_stat(fp,0,npdb,pdbf);
+ line(fp);
+
+ /* Index to first structure in a cluster */
+ snew(cndx,npdb);
+ ncluster=1;
+
+ for(i=1; (i<npdb); i++) {
+ for(j=0; (j<ncluster); j++) {
+ rmsd = rmsd_dist(pdbf[cndx[j]],pdbf[i],bRMSD);
+ if (rmsd <= cutoff) {
+ /* Structure i is in cluster j */
+ pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
+ break;
+ }
+ }
+ if (j == ncluster) {
+ /* New cluster! Cool! */
+ cndx[ncluster] = i;
+ pdbf[i]->cluster_id = ncluster;
+ ncluster++;
+ }
+ }
+ cndx[ncluster] = npdb;
+ qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
+
+ j=0;
+ cndx[j++] = 0;
+ for(i=1; (i<npdb); i++) {
+ if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id) {
+ cndx[j++] = i;
+ }
+ }
+ cndx[j] = npdb;
+ if (j != ncluster)
+ gmx_fatal(FARGS,"Consistency error: j = %d, ncluster = %d",j,ncluster);
+
+ fprintf(fp,"I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
+ ncluster,etitles[bFree],bRMSD ? "RMSD" : "distance",cutoff);
+ line(fp);
+ for(j=0; (j<ncluster); j++) {
+ fprintf(fp,"Cluster: %3d %s: %10.5f kJ/mol %3d elements\n",
+ j,etitles[bFree],
+ bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
+ cndx[j+1]-cndx[j]);
+ clust_stat(fp,cndx[j],cndx[j+1],pdbf);
+ for(k=cndx[j]; (k<cndx[j+1]); k++)
+ fprintf(fp," %3d",pdbf[k]->index);
+ fprintf(fp,"\n");
+ line(fp);
+ }
+ sfree(cndx);
+}
+
+int gmx_anadock(int argc,char *argv[])
+{
+ const char *desc[] = {
+ "[TT]g_anadock[tt] analyses the results of an Autodock run and clusters the",
+ "structures together, based on distance or RMSD. The docked energy",
+ "and free energy estimates are analysed, and for each cluster the",
+ "energy statistics are printed.[PAR]",
+ "An alternative approach to this is to cluster the structures first",
+ "using [TT]g_cluster[tt] and then sort the clusters on either lowest",
+ "energy or average energy."
+ };
+ t_filenm fnm[] = {
+ { efPDB, "-f", NULL, ffREAD },
+ { efPDB, "-ox", "cluster", ffWRITE },
+ { efXVG, "-od", "edocked", ffWRITE },
+ { efXVG, "-of", "efree", ffWRITE },
+ { efLOG, "-g", "anadock", ffWRITE }
+ };
+ output_env_t oenv;
+#define NFILE asize(fnm)
+ static gmx_bool bFree=FALSE,bRMS=TRUE;
+ static real cutoff = 0.2;
+ t_pargs pa[] = {
+ { "-free", FALSE, etBOOL, {&bFree},
+ "Use Free energy estimate from autodock for sorting the classes" },
+ { "-rms", FALSE, etBOOL, {&bRMS},
+ "Cluster on RMS or distance" },
+ { "-cutoff", FALSE, etREAL, {&cutoff},
+ "Maximum RMSD/distance for belonging to the same cluster" }
+ };
+#define NPA asize(pa)
+
+ FILE *fp;
+ t_pdbfile **pdbf=NULL;
+ int npdbf;
+
+ CopyRight(stderr,argv[0]);
+ parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa, asize(desc),desc,0,
+ NULL,&oenv);
+
+ fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
+ please_cite(stdout,"Hetenyi2002b");
+ please_cite(fp,"Hetenyi2002b");
+
+ pdbf = read_em_all(opt2fn("-f",NFILE,fnm),&npdbf);
+
+ analyse_em_all(npdbf,pdbf,opt2fn("-od",NFILE,fnm),opt2fn("-of",NFILE,fnm),
+ oenv);
+
+ cluster_em_all(fp,npdbf,pdbf,opt2fn("-ox",NFILE,fnm),
+ bFree,bRMS,cutoff);
+
+ thanx(fp);
+ ffclose(fp);
+
+ thanx(stdout);
+
+ return 0;
+}
--- /dev/null
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * 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 "sysstuff.h"
+#include "typedefs.h"
+#include "string2.h"
+#include "strdb.h"
+#include "macros.h"
+#include "smalloc.h"
+#include "mshift.h"
+#include "statutil.h"
+#include "copyrite.h"
+#include "pdbio.h"
+#include "gmx_fatal.h"
+#include "xvgr.h"
+#include "matio.h"
+#include "index.h"
+#include "gstat.h"
+#include "tpxio.h"
+#include "viewit.h"
+
+
+static int strip_dssp(char *dsspfile,int nres,
+ gmx_bool bPhobres[],real t,
+ real *acc,FILE *fTArea,
+ t_matrix *mat,int average_area[],
+ const output_env_t oenv)
+{
+ static gmx_bool bFirst=TRUE;
+ static char *ssbuf;
+ FILE *tapeout;
+ static int xsize,frame;
+ char buf[STRLEN+1];
+ char SSTP;
+ int i,nr,iacc,nresidues;
+ int naccf,naccb; /* Count hydrophobic and hydrophilic residues */
+ real iaccf,iaccb;
+ t_xpmelmt c;
+
+ tapeout=ffopen(dsspfile,"r");
+
+ /* Skip header */
+ do {
+ fgets2(buf,STRLEN,tapeout);
+ } while (strstr(buf,"KAPPA") == NULL);
+ if (bFirst) {
+ /* Since we also have empty lines in the dssp output (temp) file,
+ * and some line content is saved to the ssbuf variable,
+ * we need more memory than just nres elements. To be shure,
+ * we allocate 2*nres-1, since for each chain there is a
+ * separating line in the temp file. (At most each residue
+ * could have been defined as a separate chain.) */
+ snew(ssbuf,2*nres-1);
+ }
+
+ iaccb=iaccf=0;
+ nresidues = 0;
+ naccf = 0;
+ naccb = 0;
+ for(nr=0; (fgets2(buf,STRLEN,tapeout) != NULL); nr++) {
+ if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
+ SSTP='='; /* Chain separator sign '=' */
+ else
+ SSTP=buf[16]==' ' ? '~' : buf[16];
+ ssbuf[nr]=SSTP;
+
+ buf[39]='\0';
+
+ /* Only calculate solvent accessible area if needed */
+ if ((NULL != acc) && (buf[13] != '!'))
+ {
+ sscanf(&(buf[34]),"%d",&iacc);
+ acc[nr]=iacc;
+ /* average_area and bPhobres are counted from 0...nres-1 */
+ average_area[nresidues]+=iacc;
+ if (bPhobres[nresidues])
+ {
+ naccb++;
+ iaccb+=iacc;
+ }
+ else
+ {
+ naccf++;
+ iaccf+=iacc;
+ }
+ /* Keep track of the residue number (do not count chain separator lines '!') */
+ nresidues++;
+ }
+ }
+ ssbuf[nr]='\0';
+
+ if (bFirst) {
+ if (0 != acc)
+ fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
+
+ sprintf(mat->title,"Secondary structure");
+ mat->legend[0]=0;
+ sprintf(mat->label_x,"%s",output_env_get_time_label(oenv));
+ sprintf(mat->label_y,"Residue");
+ mat->bDiscrete=TRUE;
+ mat->ny=nr;
+ snew(mat->axis_y,nr);
+ for(i=0; i<nr; i++)
+ mat->axis_y[i]=i+1;
+ mat->axis_x=NULL;
+ mat->matrix=NULL;
+ xsize=0;
+ frame=0;
+ bFirst=FALSE;
+ }
+ if (frame>=xsize) {
+ xsize+=10;
+ srenew(mat->axis_x,xsize);
+ srenew(mat->matrix,xsize);
+ }
+ mat->axis_x[frame]=t;
+ snew(mat->matrix[frame],nr);
+ c.c2=0;
+ for(i=0; i<nr; i++) {
+ c.c1=ssbuf[i];
+ mat->matrix[frame][i] = max(0,searchcmap(mat->nmap,mat->map,c));
+ }
+ frame++;
+ mat->nx=frame;
+
+ if (fTArea)
+ fprintf(fTArea,"%10g %10g %10g\n",t,0.01*iaccb,0.01*iaccf);
+ ffclose(tapeout);
+
+ /* Return the number of lines found in the dssp file (i.e. number
+ * of redidues plus chain separator lines).
+ * This is the number of y elements needed for the area xpm file */
+ return nr;
+}
+
+static gmx_bool *bPhobics(t_atoms *atoms)
+{
+ int i,nb;
+ char **cb;
+ gmx_bool *bb;
+
+
+ nb = get_strings("phbres.dat",&cb);
+ snew(bb,atoms->nres);
+
+ for (i=0; (i<atoms->nres); i++)
+ {
+ if ( -1 != search_str(nb,cb,*atoms->resinfo[i].name) )
+ {
+ bb[i]=TRUE;
+ }
+ }
+ return bb;
+}
+
+static void check_oo(t_atoms *atoms)
+{
+ char *OOO;
+
+ int i;
+
+ OOO=strdup("O");
+
+ for(i=0; (i<atoms->nr); i++) {
+ if (strcmp(*(atoms->atomname[i]),"OXT")==0)
+ *atoms->atomname[i]=OOO;
+ else if (strcmp(*(atoms->atomname[i]),"O1")==0)
+ *atoms->atomname[i]=OOO;
+ else if (strcmp(*(atoms->atomname[i]),"OC1")==0)
+ *atoms->atomname[i]=OOO;
+ }
+}
+
+static void norm_acc(t_atoms *atoms, int nres,
+ real av_area[], real norm_av_area[])
+{
+ int i,n,n_surf;
+
+ char surffn[]="surface.dat";
+ char **surf_res, **surf_lines;
+ double *surf;
+
+ n_surf = get_lines(surffn, &surf_lines);
+ snew(surf, n_surf);
+ snew(surf_res, n_surf);
+ for(i=0; (i<n_surf); i++) {
+ snew(surf_res[i], 5);
+ sscanf(surf_lines[i],"%s %lf",surf_res[i],&surf[i]);
+ }
+
+ for(i=0; (i<nres); i++) {
+ n = search_str(n_surf,surf_res,*atoms->resinfo[i].name);
+ if ( n != -1)
+ norm_av_area[i] = av_area[i] / surf[n];
+ else
+ fprintf(stderr,"Residue %s not found in surface database (%s)\n",
+ *atoms->resinfo[i].name,surffn);
+ }
+}
+
+void prune_ss_legend(t_matrix *mat)
+{
+ gmx_bool *present;
+ int *newnum;
+ int i,r,f,newnmap;
+ t_mapping *newmap;
+
+ snew(present,mat->nmap);
+ snew(newnum,mat->nmap);
+
+ for(f=0; f<mat->nx; f++)
+ for(r=0; r<mat->ny; r++)
+ present[mat->matrix[f][r]]=TRUE;
+
+ newnmap=0;
+ newmap=NULL;
+ for (i=0; i<mat->nmap; i++) {
+ newnum[i]=-1;
+ if (present[i]) {
+ newnum[i]=newnmap;
+ newnmap++;
+ srenew(newmap,newnmap);
+ newmap[newnmap-1]=mat->map[i];
+ }
+ }
+ if (newnmap!=mat->nmap) {
+ mat->nmap=newnmap;
+ mat->map=newmap;
+ for(f=0; f<mat->nx; f++)
+ for(r=0; r<mat->ny; r++)
+ mat->matrix[f][r]=newnum[mat->matrix[f][r]];
+ }
+}
+
+void write_sas_mat(const char *fn,real **accr,int nframe,int nres,t_matrix *mat)
+{
+ real lo,hi;
+ int i,j,nlev;
+ t_rgb rlo={1,1,1}, rhi={0,0,0};
+ FILE *fp;
+
+ if(fn) {
+ hi=lo=accr[0][0];
+ for(i=0; i<nframe; i++)
+ for(j=0; j<nres; j++) {
+ lo=min(lo,accr[i][j]);
+ hi=max(hi,accr[i][j]);
+ }
+ fp=ffopen(fn,"w");
+ nlev=hi-lo+1;
+ write_xpm(fp,0,"Solvent Accessible Surface","Surface (A^2)",
+ "Time","Residue Index",nframe,nres,
+ mat->axis_x,mat->axis_y,accr,lo,hi,rlo,rhi,&nlev);
+ ffclose(fp);
+ }
+}
+
+void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
+ const output_env_t oenv)
+{
+ FILE *fp;
+ t_mapping *map;
+ int f,r,*count,*total,ss_count,total_count;
+ size_t s;
+ const char** leg;
+
+ map=mat->map;
+ snew(count,mat->nmap);
+ snew(total,mat->nmap);
+ snew(leg,mat->nmap+1);
+ leg[0]="Structure";
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ leg[s+1]=strdup(map[s].desc);
+ }
+
+ fp=xvgropen(outfile,"Secondary Structure",
+ output_env_get_xvgr_tlabel(oenv),"Number of Residues",oenv);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp,"@ subtitle \"Structure = ");
+ }
+ for(s=0; s<strlen(ss_string); s++)
+ {
+ if (s>0)
+ {
+ fprintf(fp," + ");
+ }
+ for(f=0; f<mat->nmap; f++)
+ {
+ if (ss_string[s]==map[f].code.c1)
+ {
+ fprintf(fp,"%s",map[f].desc);
+ }
+ }
+ }
+ fprintf(fp,"\"\n");
+ xvgr_legend(fp,mat->nmap+1,leg,oenv);
+
+ total_count = 0;
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ total[s]=0;
+ }
+ for(f=0; f<mat->nx; f++)
+ {
+ ss_count=0;
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ count[s]=0;
+ }
+ for(r=0; r<mat->ny; r++)
+ {
+ count[mat->matrix[f][r]]++;
+ total[mat->matrix[f][r]]++;
+ }
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ if (strchr(ss_string,map[s].code.c1))
+ {
+ ss_count+=count[s];
+ total_count += count[s];
+ }
+ }
+ fprintf(fp,"%8g %5d",mat->axis_x[f],ss_count);
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ fprintf(fp," %5d",count[s]);
+ }
+ fprintf(fp,"\n");
+ }
+ /* now print column totals */
+ fprintf(fp, "%-8s %5d", "# Totals", total_count);
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ fprintf(fp," %5d",total[s]);
+ }
+ fprintf(fp,"\n");
+
+ /* now print percentages */
+ fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));
+ for(s=0; s<(size_t)mat->nmap; s++)
+ {
+ fprintf(fp," %5.2f",total[s] / (real) (mat->nx * mat->ny));
+ }
+ fprintf(fp,"\n");
+
+ ffclose(fp);
+ sfree(leg);
+ sfree(count);
+}
+
+int gmx_do_dssp(int argc,char *argv[])
+{
+ const char *desc[] = {
+ "[TT]do_dssp[tt] ",
+ "reads a trajectory file and computes the secondary structure for",
+ "each time frame ",
+ "calling the dssp program. If you do not have the dssp program,",
+ "get it from http://swift.cmbi.ru.nl/gv/dssp. [TT]do_dssp[tt] assumes ",
+ "that the dssp executable is located in ",
+ "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
+ "set an environment variable [TT]DSSP[tt] pointing to the dssp",
+ "executable, e.g.: [PAR]",
+ "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
+ "Since version 2.0.0, dssp is invoked with a syntax that differs",
+ "from earlier versions. If you have an older version of dssp,",
+ "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
+ "By default, do_dssp uses the syntax introduced with version 2.0.0.",
+ "Even newer versions (which at the time of writing are not yet released)",
+ "are assumed to have the same syntax as 2.0.0.[PAR]",
+ "The structure assignment for each residue and time is written to an",
+ "[TT].xpm[tt] matrix file. This file can be visualized with for instance",
+ "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
+ "Individual chains are separated by light grey lines in the [TT].xpm[tt] and",
+ "postscript files.",
+ "The number of residues with each secondary structure type and the",
+ "total secondary structure ([TT]-sss[tt]) count as a function of",
+ "time are also written to file ([TT]-sc[tt]).[PAR]",
+ "Solvent accessible surface (SAS) per residue can be calculated, both in",
+ "absolute values (A^2) and in fractions of the maximal accessible",
+ "surface of a residue. The maximal accessible surface is defined as",
+ "the accessible surface of a residue in a chain of glycines.",
+ "[BB]Note[bb] that the program [TT]g_sas[tt] can also compute SAS",
+ "and that is more efficient.[PAR]",
+ "Finally, this program can dump the secondary structure in a special file",
+ "[TT]ssdump.dat[tt] for usage in the program [TT]g_chi[tt]. Together",
+ "these two programs can be used to analyze dihedral properties as a",
+ "function of secondary structure type."
+ };
+ static gmx_bool bVerbose;
+ static const char *ss_string="HEBT";
+ static int dsspVersion=2;
+ t_pargs pa[] = {
+ { "-v", FALSE, etBOOL, {&bVerbose},
+ "HIDDENGenerate miles of useless information" },
+ { "-sss", FALSE, etSTR, {&ss_string},
+ "Secondary structures for structure count"},
+ { "-ver", FALSE, etINT, {&dsspVersion},
+ "DSSP major version. Syntax changed with version 2"}
+ };
+
+ t_trxstatus *status;
+ FILE *tapein;
+ FILE *ss,*acc,*fTArea,*tmpf;
+ const char *fnSCount,*fnArea,*fnTArea,*fnAArea;
+ const char *leg[] = { "Phobic", "Phylic" };
+ t_topology top;
+ int ePBC;
+ t_atoms *atoms;
+ t_matrix mat;
+ int nres,nr0,naccr,nres_plus_separators;
+ gmx_bool *bPhbres,bDoAccSurf;
+ real t;
+ int i,j,natoms,nframe=0;
+ matrix box;
+ int gnx;
+ char *grpnm,*ss_str;
+ atom_id *index;
+ rvec *xp,*x;
+ int *average_area;
+ real **accr,*accr_ptr=NULL,*av_area, *norm_av_area;
+ char pdbfile[32],tmpfile[32],title[256];
+ char dssp[256];
+ const char *dptr;
+ output_env_t oenv;
+ gmx_rmpbc_t gpbc=NULL;
+
+ t_filenm fnm[] = {
+ { efTRX, "-f", NULL, ffREAD },
+ { efTPS, NULL, NULL, ffREAD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efDAT, "-ssdump", "ssdump", ffOPTWR },
+ { efMAP, "-map", "ss", ffLIBRD },
+ { efXPM, "-o", "ss", ffWRITE },
+ { efXVG, "-sc", "scount", ffWRITE },
+ { efXPM, "-a", "area", ffOPTWR },
+ { efXVG, "-ta", "totarea", ffOPTWR },
+ { efXVG, "-aa", "averarea",ffOPTWR }
+ };
+#define NFILE asize(fnm)
+
+ CopyRight(stderr,argv[0]);
+ parse_common_args(&argc,argv,
+ PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
+ NFILE,fnm, asize(pa),pa, asize(desc),desc,0,NULL,&oenv);
+ fnSCount= opt2fn("-sc",NFILE,fnm);
+ fnArea = opt2fn_null("-a", NFILE,fnm);
+ fnTArea = opt2fn_null("-ta",NFILE,fnm);
+ fnAArea = opt2fn_null("-aa",NFILE,fnm);
+ bDoAccSurf=(fnArea || fnTArea || fnAArea);
+
+ read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
+ atoms=&(top.atoms);
+ check_oo(atoms);
+ bPhbres = bPhobics(atoms);
+
+ get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
+ nres=0;
+ nr0=-1;
+ for(i=0; (i<gnx); i++) {
+ if (atoms->atom[index[i]].resind != nr0) {
+ nr0=atoms->atom[index[i]].resind;
+ nres++;
+ }
+ }
+ fprintf(stderr,"There are %d residues in your selected group\n",nres);
+
+ strcpy(pdbfile,"ddXXXXXX");
+ gmx_tmpnam(pdbfile);
+ if ((tmpf = fopen(pdbfile,"w")) == NULL) {
+ sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
+ gmx_tmpnam(pdbfile);
+ if ((tmpf = fopen(pdbfile,"w")) == NULL)
+ gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);
+ }
+ else
+ fclose(tmpf);
+
+ strcpy(tmpfile,"ddXXXXXX");
+ gmx_tmpnam(tmpfile);
+ if ((tmpf = fopen(tmpfile,"w")) == NULL) {
+ sprintf(tmpfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
+ gmx_tmpnam(tmpfile);
+ if ((tmpf = fopen(tmpfile,"w")) == NULL)
+ gmx_fatal(FARGS,"Can not open tmp file %s",tmpfile);
+ }
+ else
+ fclose(tmpf);
+
+ if ((dptr=getenv("DSSP")) == NULL)
+ dptr="/usr/local/bin/dssp";
+ if (!gmx_fexist(dptr))
+ gmx_fatal(FARGS,"DSSP executable (%s) does not exist (use setenv DSSP)",
+ dptr);
+ if (dsspVersion >= 2)
+ {
+ if (dsspVersion > 2)
+ {
+ printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
+ }
+
+ sprintf(dssp,"%s -i %s -o %s > /dev/null %s",
+ dptr,pdbfile,tmpfile,bVerbose?"":"2> /dev/null");
+ }
+ else
+ {
+ sprintf(dssp,"%s %s %s %s > /dev/null %s",
+ dptr,bDoAccSurf?"":"-na",pdbfile,tmpfile,bVerbose?"":"2> /dev/null");
+
+ }
+ fprintf(stderr,"dssp cmd='%s'\n",dssp);
+
+ if (fnTArea) {
+ fTArea=xvgropen(fnTArea,"Solvent Accessible Surface Area",
+ output_env_get_xvgr_tlabel(oenv),"Area (nm\\S2\\N)",oenv);
+ xvgr_legend(fTArea,2,leg,oenv);
+ } else
+ fTArea=NULL;
+
+ mat.map=NULL;
+ mat.nmap=getcmap(libopen(opt2fn("-map",NFILE,fnm)),
+ opt2fn("-map",NFILE,fnm),&(mat.map));
+
+ natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
+ if (natoms > atoms->nr)
+ gmx_fatal(FARGS,"\nTrajectory does not match topology!");
+ if (gnx > natoms)
+ gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
+
+ snew(average_area, atoms->nres);
+ snew(av_area , atoms->nres);
+ snew(norm_av_area, atoms->nres);
+ accr=NULL;
+ naccr=0;
+
+ gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
+ do {
+ t = output_env_conv_time(oenv,t);
+ if (bDoAccSurf && nframe>=naccr) {
+ naccr+=10;
+ srenew(accr,naccr);
+ for(i=naccr-10; i<naccr; i++)
+ snew(accr[i],2*atoms->nres-1);
+ }
+ gmx_rmpbc(gpbc,natoms,box,x);
+ tapein=ffopen(pdbfile,"w");
+ write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,' ',-1,gnx,index,NULL,TRUE);
+ ffclose(tapein);
+
+#ifdef GMX_NO_SYSTEM
+ printf("Warning-- No calls to system(3) supported on this platform.");
+ printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
+ exit(1);
+#else
+ if(0 != system(dssp))
+ {
+ gmx_fatal(FARGS,"Failed to execute command: %s\n",
+ "Try specifying your dssp version with the -ver option.",dssp);
+ }
+#endif
+
+ /* strip_dssp returns the number of lines found in the dssp file, i.e.
+ * the number of residues plus the separator lines */
+
+ if (bDoAccSurf)
+ accr_ptr = accr[nframe];
+
+ nres_plus_separators = strip_dssp(tmpfile,nres,bPhbres,t,
+ accr_ptr,fTArea,&mat,average_area,oenv);
+ remove(tmpfile);
+ remove(pdbfile);
+ nframe++;
+ } while(read_next_x(oenv,status,&t,natoms,x,box));
+ fprintf(stderr,"\n");
+ close_trj(status);
+ if (fTArea)
+ ffclose(fTArea);
+ gmx_rmpbc_done(gpbc);
+
+ prune_ss_legend(&mat);
+
+ ss=opt2FILE("-o",NFILE,fnm,"w");
+ mat.flags = 0;
+ write_xpm_m(ss,mat);
+ ffclose(ss);
+
+ if (opt2bSet("-ssdump",NFILE,fnm))
+ {
+ ss = opt2FILE("-ssdump",NFILE,fnm,"w");
+ snew(ss_str,nres+1);
+ fprintf(ss,"%d\n",nres);
+ for(j=0; j<mat.nx; j++)
+ {
+ for(i=0; (i<mat.ny); i++)
+ {
+ ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
+ }
+ ss_str[i] = '\0';
+ fprintf(ss,"%s\n",ss_str);
+ }
+ ffclose(ss);
+ sfree(ss_str);
+ }
+ analyse_ss(fnSCount,&mat,ss_string,oenv);
+
+ if (bDoAccSurf) {
+ write_sas_mat(fnArea,accr,nframe,nres_plus_separators,&mat);
+
+ for(i=0; i<atoms->nres; i++)
+ av_area[i] = (average_area[i] / (real)nframe);
+
+ norm_acc(atoms, nres, av_area, norm_av_area);
+
+ if (fnAArea) {
+ acc=xvgropen(fnAArea,"Average Accessible Area",
+ "Residue","A\\S2",oenv);
+ for(i=0; (i<nres); i++)
+ fprintf(acc,"%5d %10g %10g\n",i+1,av_area[i], norm_av_area[i]);
+ ffclose(acc);
+ }
+ }
+
+ view_all(oenv, NFILE, fnm);
+
+ thanx(stderr);
+
+ return 0;
+}
--- /dev/null
+/*
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ *
+ * The make_edi program was generously contributed by Oliver Lange, based
+ * on the code from g_anaeig. You can reach him as olange@gwdg.de.
+ *
+ * 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:
+ * Gromacs Runs One Microsecond At Cannonball Speeds
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "gmx_header_config.h"
+
+#include <math.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "readinp.h"
+#include "statutil.h"
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "gmx_fatal.h"
+#include "vec.h"
+#include "pbc.h"
+#include "copyrite.h"
+#include "futil.h"
+#include "statutil.h"
+#include "pdbio.h"
+#include "confio.h"
+#include "tpxio.h"
+#include "matio.h"
+#include "mshift.h"
+#include "xvgr.h"
+#include "do_fit.h"
+#include "rmpbc.h"
+#include "txtdump.h"
+#include "eigio.h"
+#include "index.h"
+
+/* Suppress Cygwin compiler warnings from using newlib version of
+ * ctype.h */
+#ifdef GMX_CYGWIN
+#undef isdigit
+#endif
+
+typedef struct
+{
+ real deltaF0;
+ gmx_bool bHarmonic;
+ gmx_bool bConstForce; /* Do constant force flooding instead of
+ evaluating a flooding potential */
+ real tau;
+ real deltaF;
+ real kT;
+ real constEfl;
+ real alpha2;
+} t_edflood;
+
+
+/* This type is for the average, reference, target, and origin structure */
+typedef struct edix
+{
+ int nr; /* number of atoms this structure contains */
+ int *anrs; /* atom index numbers */
+ rvec *x; /* positions */
+ real *sqrtm; /* sqrt of the masses used for mass-
+ * weighting of analysis */
+} t_edix;
+
+
+typedef struct edipar
+{
+ int nini; /* total Nr of atoms */
+ gmx_bool fitmas; /* true if trans fit with cm */
+ gmx_bool pcamas; /* true if mass-weighted PCA */
+ int presteps; /* number of steps to run without any
+ * perturbations ... just monitoring */
+ int outfrq; /* freq (in steps) of writing to edo */
+ int maxedsteps; /* max nr of steps per cycle */
+ struct edix sref; /* reference positions, to these fitting
+ * will be done */
+ struct edix sav; /* average positions */
+ struct edix star; /* target positions */
+ struct edix sori; /* origin positions */
+ real slope; /* minimal slope in acceptance radexp */
+ int ned; /* Nr of atoms in essdyn buffer */
+ t_edflood flood; /* parameters especially for flooding */
+} t_edipar;
+
+
+
+void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
+ edx->nr=natoms;
+ edx->anrs=index;
+ edx->x=pos;
+}
+
+void write_t_edx(FILE *fp, struct edix edx, const char *comment) {
+ /*here we copy only the pointers into the t_edx struct
+ no data is copied and edx.box is ignored */
+ int i;
+ fprintf(fp,"#%s \n %d \n",comment,edx.nr);
+ for (i=0;i<edx.nr;i++) {
+ fprintf(fp,"%d %f %f %f\n",(edx.anrs)[i]+1,(edx.x)[i][XX],(edx.x)[i][YY],(edx.x)[i][ZZ]);
+ }
+}
+
+int sscan_list(int *list[], const char *str, const char *listname) {
+ /*this routine scans a string of the form 1,3-6,9 and returns the
+ selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
+ memory for this list will be allocated in this routine -- sscan_list expects *list to
+ be a NULL-Pointer
+
+ listname is a string used in the errormessage*/
+
+
+ int i,istep;
+ char c;
+ char *pos,*startpos,*step;
+ int n=strlen(str);
+
+ /*enums to define the different lexical stati */
+ enum { sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange };
+
+ int status=sBefore; /*status of the deterministic automat to scan str */
+ int number=0;
+ int end_number=0;
+
+ char *start=NULL; /*holds the string of the number behind a ','*/
+ char *end=NULL; /*holds the string of the number behind a '-' */
+
+ int nvecs=0; /* counts the number of vectors in the list*/
+
+ step=NULL;
+ snew(pos,n+4);
+ startpos=pos;
+ strcpy(pos,str);
+ pos[n]=',';
+ pos[n+1]='1';
+ pos[n+2]='\0';
+
+ *list = NULL;
+
+ while ((c=*pos)!=0) {
+ switch(status) {
+ /* expect a number */
+ case sBefore: if (isdigit(c)) {
+ start=pos;
+ status=sNumber;
+ break;
+ } else status=sError; break;
+
+ /* have read a number, expect ',' or '-' */
+ case sNumber: if (c==',') {
+ /*store number*/
+ srenew(*list,nvecs+1);
+ (*list)[nvecs++]=number=strtol(start,NULL,10);
+ status=sBefore;
+ if (number==0)
+ status=sZero;
+ break;
+ } else if (c=='-') { status=sMinus; break; }
+ else if (isdigit(c)) break;
+ else status=sError; break;
+
+ /* have read a '-' -> expect a number */
+ case sMinus:
+ if (isdigit(c)) {
+ end=pos;
+ status=sRange; break;
+ } else status=sError; break;
+
+ case sSteppedRange:
+ if (isdigit(c)) {
+ if (step) {
+ status=sError; break;
+ } else
+ step=pos;
+ status=sRange;
+ break;
+ } else status=sError; break;
+
+ /* have read the number after a minus, expect ',' or ':' */
+ case sRange:
+ if (c==',') {
+ /*store numbers*/
+ end_number=strtol(end,NULL,10);
+ number=strtol(start,NULL,10);
+ status=sBefore;
+ if (number==0) {
+ status=sZero; break;
+ }
+ if (end_number<=number) {
+ status=sSmaller; break;
+ }
+ srenew(*list,nvecs+end_number-number+1);
+ if (step) {
+ istep=strtol(step,NULL,10);
+ step=NULL;
+ } else istep=1;
+ for (i=number;i<=end_number;i+=istep)
+ (*list)[nvecs++]=i;
+ break;
+ } else if (c==':') {
+ status = sSteppedRange;
+ break;
+ } else if (isdigit(c)) break; else status=sError; break;
+
+ /* format error occured */
+ case sError:
+ gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));
+ break;
+ /* logical error occured */
+ case sZero:
+ gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid",listname,pos-startpos);
+ break;
+ case sSmaller:
+ gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: second index %d is not bigger than %d",listname,pos-startpos,end_number,number);
+ break;
+ }
+ ++pos; /* read next character */
+ } /*scanner has finished */
+
+ /* append zero to list of eigenvectors */
+ srenew(*list,nvecs+1);
+ (*list)[nvecs]=0;
+ sfree(startpos);
+ return nvecs;
+} /*sscan_list*/
+
+void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec, const char *grouptitle,real steps[]) {
+/* eig_list is a zero-terminated list of indices into the eigvecs array.
+ eigvecs are coordinates of eigenvectors
+ grouptitle to write in the comment line
+ steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC
+*/
+
+ int n=0,i; rvec x;
+ real sum;
+ while (eig_list[n++]); /*count selected eigenvecs*/
+
+ fprintf(fp,"# NUMBER OF EIGENVECTORS + %s\n %d\n",grouptitle,n-1);
+
+ /* write list of eigenvector indicess */
+ for(n=0;eig_list[n];n++) {
+ if (steps)
+ fprintf(fp,"%8d %g\n",eig_list[n],steps[n]);
+ else
+ fprintf(fp,"%8d %g\n",eig_list[n],1.0);
+ }
+ n=0;
+
+ /* dump coordinates of the selected eigenvectors */
+ while (eig_list[n]) {
+ sum=0;
+ for (i=0; i<natoms; i++) {
+ if (eig_list[n]>nvec)
+ gmx_fatal(FARGS,"Selected eigenvector %d is higher than maximum number %d of available eigenvectors",eig_list[n],nvec);
+ copy_rvec(eigvecs[eig_list[n]-1][i],x);
+ sum+=norm2(x);
+ fprintf(fp,"%8.5f %8.5f %8.5f\n",x[XX],x[YY],x[ZZ]);
+ }
+ n++;
+ }
+}
+
+
+/*enum referring to the different lists of eigenvectors*/
+enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON, evNr };
+#define oldMAGIC 666
+#define MAGIC 670
+
+
+void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
+ int nvec, int *eig_listen[], real* evStepList[])
+{
+/* write edi-file */
+
+ /*Header*/
+ fprintf(fp,"#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
+ MAGIC,edpars->nini,edpars->fitmas,edpars->pcamas);
+ fprintf(fp,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
+ edpars->outfrq,edpars->maxedsteps,edpars->slope);
+ fprintf(fp,"#PRESTEPS\n %d\n#DELTA_F0\n %f\n#INIT_DELTA_F\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n#KT\n %f\n#HARMONIC\n %d\n#CONST_FORCE_FLOODING\n %d\n",
+ edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,
+ edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic,edpars->flood.bConstForce);
+
+ /* Average and reference positions */
+ write_t_edx(fp,edpars->sref,"NREF, XREF");
+ write_t_edx(fp,edpars->sav,"NAV, XAV");
+
+ /*Eigenvectors */
+
+ write_eigvec(fp, edpars->ned, eig_listen[evMON],eigvecs,nvec,"COMPONENTS GROUP 1",NULL);
+ write_eigvec(fp, edpars->ned, eig_listen[evLINFIX],eigvecs,nvec,"COMPONENTS GROUP 2",evStepList[evLINFIX]);
+ write_eigvec(fp, edpars->ned, eig_listen[evLINACC],eigvecs,nvec,"COMPONENTS GROUP 3",evStepList[evLINACC]);
+ write_eigvec(fp, edpars->ned, eig_listen[evRADFIX],eigvecs,nvec,"COMPONENTS GROUP 4",evStepList[evRADFIX]);
+ write_eigvec(fp, edpars->ned, eig_listen[evRADACC],eigvecs,nvec,"COMPONENTS GROUP 5",NULL);
+ write_eigvec(fp, edpars->ned, eig_listen[evRADCON],eigvecs,nvec,"COMPONENTS GROUP 6",NULL);
+ write_eigvec(fp, edpars->ned, eig_listen[evFLOOD],eigvecs,nvec,"COMPONENTS GROUP 7",evStepList[evFLOOD]);
+
+
+ /*Target and Origin positions */
+ write_t_edx(fp,edpars->star,"NTARGET, XTARGET");
+ write_t_edx(fp,edpars->sori,"NORIGIN, XORIGIN");
+}
+
+int read_conffile(const char *confin,char *title,rvec *x[])
+{
+/* read coordinates out of STX file */
+ int natoms;
+ t_atoms confat;
+ matrix box;
+ printf("read coordnumber from file %s\n",confin);
+ get_stx_coordnum(confin,&natoms);
+ printf("number of coordinates in file %d\n",natoms);
+/* if (natoms != ncoords)
+ gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
+ " does not match topology (= %d)",
+ confin,natoms,ncoords);
+ else {*/
+ /* make space for coordinates and velocities */
+ init_t_atoms(&confat,natoms,FALSE);
+ snew(*x,natoms);
+ read_stx_conf(confin,title,&confat,*x,NULL,NULL,box);
+ return natoms;
+}
+
+
+void read_eigenvalues(int vecs[],const char *eigfile, real values[],
+ gmx_bool bHesse, real kT)
+{
+ int neig,nrow,i;
+ double **eigval;
+
+ neig = read_xvg(eigfile,&eigval,&nrow);
+
+ fprintf(stderr,"Read %d eigenvalues\n",neig);
+ for(i=bHesse ? 6 : 0 ; i<neig; i++) {
+ if (eigval[1][i] < -0.001 && bHesse)
+ fprintf(stderr,
+ "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n",eigval[1][i]);
+
+ if (eigval[1][i] < 0)
+ eigval[1][i] = 0;
+ }
+ if (bHesse)
+ for (i=0; vecs[i]; i++) {
+ if (vecs[i]<7)
+ gmx_fatal(FARGS,"ERROR: You have chosen one of the first 6 eigenvectors of the HESSE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
+ values[i]=eigval[1][vecs[i]-1]/kT;
+ }
+ else
+ for (i=0; vecs[i]; i++) {
+ if (vecs[i]>(neig-6))
+ gmx_fatal(FARGS,"ERROR: You have chosen one of the last 6 eigenvectors of the COVARIANCE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
+ values[i]=1/eigval[1][vecs[i]-1];
+ }
+ /* free memory */
+ for (i=0; i<nrow; i++)
+ sfree(eigval[i]);
+ sfree(eigval);
+}
+
+
+static real *scan_vecparams(const char *str,const char * par, int nvecs)
+{
+ char f0[256],f1[256]; /*format strings adapted every pass of the loop*/
+ double d,tcap=0;
+ int i;
+ real *vec_params;
+
+ snew(vec_params,nvecs);
+ if (str) {
+ f0[0] = '\0';
+ for(i=0; (i<nvecs); i++) {
+ strcpy(f1,f0); /*f0 is the format string for the "to-be-ignored" numbers*/
+ strcat(f1,"%lf"); /*and f1 to read the actual number in this pass of the loop*/
+ if (sscanf(str,f1,&d) != 1)
+ gmx_fatal(FARGS,"Not enough elements for %s parameter (I need %d)",par,nvecs);
+ vec_params[i] = d;
+ tcap += d;
+ strcat(f0,"%*s");
+ }
+ }
+ return vec_params;
+}
+
+
+void init_edx(struct edix *edx) {
+ edx->nr=0;
+ snew(edx->x,1);
+ snew(edx->anrs,1);
+}
+
+void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
+ atom_id igro[],rvec *x,const char* structure)
+{
+/* filter2edx copies coordinates from x to edx which are given in index
+*/
+
+ int pos,i;
+ int ix=edx->nr;
+ edx->nr+=nindex;
+ srenew(edx->x,edx->nr);
+ srenew(edx->anrs,edx->nr);
+ for (i=0;i<nindex;i++,ix++) {
+ for (pos=0; pos<ngro-1 && igro[pos]!=index[i] ; ++pos) {}; /*search element in igro*/
+ if (igro[pos]!=index[i])
+ gmx_fatal(FARGS,"Couldn't find atom with index %d in structure %s",index[i],structure);
+ edx->anrs[ix]=index[i];
+ copy_rvec(x[pos],edx->x[ix]);
+ }
+}
+
+void get_structure(t_atoms *atoms,const char *IndexFile,
+ const char *StructureFile,struct edix *edx,int nfit,
+ atom_id ifit[],int nav, atom_id index[])
+{
+ atom_id *igro; /*index corresponding to target or origin structure*/
+ int ngro;
+ int ntar;
+ rvec *xtar;
+ char title[STRLEN];
+ char* grpname;
+
+
+ ntar=read_conffile(StructureFile,title,&xtar);
+ printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
+ ntar,StructureFile);
+ get_index(atoms,IndexFile,1,&ngro,&igro,&grpname);
+ if (ngro!=ntar)
+ gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
+ init_edx(edx);
+ filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
+
+ /* If average and reference/fitting structure differ, append the average structure as well */
+ if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
+ filter2edx(edx,nav,index,ngro,igro,xtar,StructureFile);
+}
+
+int gmx_make_edi(int argc,char *argv[])
+{
+
+ static const char *desc[] = {
+ "[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
+ "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
+ "normal modes anaysis ([TT]g_nmeig[tt]).",
+ "ED sampling can be used to manipulate the position along collective coordinates",
+ "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
+ "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
+ "the system to explore new regions along these collective coordinates. A number",
+ "of different algorithms are implemented to drive the system along the eigenvectors",
+ "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
+ "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
+ "or to only monitor the projections of the positions onto",
+ "these coordinates ([TT]-mon[tt]).[PAR]",
+ "References:[BR]",
+ "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
+ "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
+ "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
+ "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
+ "Towards an exhaustive sampling of the configurational spaces of the ",
+ "two forms of the peptide hormone guanylin,",
+ "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
+ "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
+ "An extended sampling of the configurational space of HPr from E. coli",
+ "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
+ "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
+ "reference structure, target positions, etc.[PAR]",
+
+ "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
+ "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
+ "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
+ "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
+ "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
+ "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
+ "(steps in the desired direction will be accepted, others will be rejected).",
+ "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
+ "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
+ "to read in a structure file that defines an external origin.[PAR]",
+ "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
+ "towards a target structure specified with [TT]-tar[tt].[PAR]",
+ "NOTE: each eigenvector can be selected only once. [PAR]",
+ "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].edo[tt] file[PAR]",
+ "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
+ "cycle will be started if the spontaneous increase of the radius (in nm/step)",
+ "is less than the value specified.[PAR]",
+ "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
+ "before a new cycle is started.[PAR]",
+ "Note on the parallel implementation: since ED sampling is a 'global' thing",
+ "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
+ "is not very parallel-friendly from an implentation point of view. Because",
+ "parallel ED requires some extra communication, expect the performance to be",
+ "lower as in a free MD simulation, especially on a large number of nodes. [PAR]",
+ "All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a .edo file. In the output",
+ "file, per OUTFRQ step the following information is present: [PAR]",
+ "[TT]*[tt] the step number[BR]",
+ "[TT]*[tt] the number of the ED dataset. ([BB]Note[bb] that you can impose multiple ED constraints in",
+ "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
+ "first. The constraints are applied in the order they appear in the [TT].edi[tt] file.) [BR]",
+ "[TT]*[tt] RMSD (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
+ "* projections of the positions onto selected eigenvectors[BR]",
+ "[PAR][PAR]",
+ "FLOODING:[PAR]",
+ "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
+ "which will lead to extra forces expelling the structure out of the region described",
+ "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
+ "is kept in that region.",
+ "[PAR]",
+ "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
+ "It can be changed with [TT]-ori[tt] to an arbitrary position in configurational space.",
+ "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
+ "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
+ "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
+ "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
+ "To use constant Efl set [TT]-tau[tt] to zero.",
+ "[PAR]",
+ "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
+ "to give good results for most standard cases in flooding of proteins.",
+ "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
+ "increase, this is mimicked by [GRK]alpha[grk] > 1.",
+ "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
+ "[PAR]",
+ "RESTART and FLOODING:",
+ "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
+ "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
+ };
+
+ /* Save all the params in this struct and then save it in an edi file.
+ * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
+ */
+ static t_edipar edi_params;
+
+ enum { evStepNr = evRADFIX + 1};
+ static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
+ static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
+ static const char* evParams[evStepNr] ={NULL,NULL};
+ static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
+ static const char* ConstForceStr;
+ static real* evStepList[evStepNr];
+ static real radfix=0.0;
+ static real deltaF0=150;
+ static real deltaF=0;
+ static real tau=.1;
+ static real constEfl=0.0;
+ static real alpha=1;
+ static int eqSteps=0;
+ static int* listen[evNr];
+ static real T=300.0;
+ const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
+ static gmx_bool bRestrain = FALSE;
+ static gmx_bool bHesse=FALSE;
+ static gmx_bool bHarmonic=FALSE;
+ t_pargs pa[] = {
+ { "-mon", FALSE, etSTR, {&evSelections[evMON]},
+ "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
+ { "-linfix", FALSE, etSTR, {&evSelections[0]},
+ "Indices of eigenvectors for fixed increment linear sampling" },
+ { "-linacc", FALSE, etSTR, {&evSelections[1]},
+ "Indices of eigenvectors for acceptance linear sampling" },
+ { "-radfix", FALSE, etSTR, {&evSelections[3]},
+ "Indices of eigenvectors for fixed increment radius expansion" },
+ { "-radacc", FALSE, etSTR, {&evSelections[4]},
+ "Indices of eigenvectors for acceptance radius expansion" },
+ { "-radcon", FALSE, etSTR, {&evSelections[5]},
+ "Indices of eigenvectors for acceptance radius contraction" },
+ { "-flood", FALSE, etSTR, {&evSelections[2]},
+ "Indices of eigenvectors for flooding"},
+ { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
+ "Freqency (in steps) of writing output in [TT].edo[tt] file" },
+ { "-slope", FALSE, etREAL, { &edi_params.slope},
+ "Minimal slope in acceptance radius expansion"},
+ { "-linstep", FALSE, etSTR, {&evParams[0]},
+ "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
+ { "-accdir", FALSE, etSTR, {&evParams[1]},
+ "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
+ { "-radstep", FALSE, etREAL, {&radfix},
+ "Stepsize (nm/step) for fixed increment radius expansion"},
+ { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
+ "Maximum number of steps per cycle" },
+ { "-eqsteps", FALSE, etINT, {&eqSteps},
+ "Number of steps to run without any perturbations "},
+ { "-deltaF0", FALSE,etREAL, {&deltaF0},
+ "Target destabilization energy for flooding"},
+ { "-deltaF", FALSE, etREAL, {&deltaF},
+ "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
+ { "-tau", FALSE, etREAL, {&tau},
+ "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
+ { "-Eflnull", FALSE, etREAL, {&constEfl},
+ "The starting value of the flooding strength. The flooding strength is updated "
+ "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
+ { "-T", FALSE, etREAL, {&T},
+ "T is temperature, the value is needed if you want to do flooding "},
+ { "-alpha",FALSE,etREAL,{&alpha},
+ "Scale width of gaussian flooding potential with alpha^2 "},
+ { "-restrain",FALSE, etBOOL, {&bRestrain},
+ "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
+ { "-hessian",FALSE, etBOOL, {&bHesse},
+ "The eigenvectors and eigenvalues are from a Hessian matrix"},
+ { "-harmonic",FALSE, etBOOL, {&bHarmonic},
+ "The eigenvalues are interpreted as spring constant"},
+ { "-constF", FALSE, etSTR, {&ConstForceStr},
+ "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
+ "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
+ };
+#define NPA asize(pa)
+
+ rvec *xref1;
+ int nvec1,*eignr1=NULL;
+ rvec *xav1,**eigvec1=NULL;
+ t_atoms *atoms=NULL;
+ int nav; /* Number of atoms in the average structure */
+ char *grpname;
+ const char *indexfile;
+ int i;
+ atom_id *index,*ifit;
+ int nfit; /* Number of atoms in the reference/fit structure */
+ int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
+ int nvecs;
+ real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */
+
+ const char *EdiFile;
+ const char *TargetFile;
+ const char *OriginFile;
+ const char *EigvecFile;
+
+ output_env_t oenv;
+
+ /*to read topology file*/
+ t_topology top;
+ int ePBC;
+ char title[STRLEN];
+ matrix topbox;
+ rvec *xtop;
+ gmx_bool bTop, bFit1;
+
+ t_filenm fnm[] = {
+ { efTRN, "-f", "eigenvec", ffREAD },
+ { efXVG, "-eig", "eigenval", ffOPTRD },
+ { efTPS, NULL, NULL, ffREAD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efSTX, "-tar", "target", ffOPTRD},
+ { efSTX, "-ori", "origin", ffOPTRD},
+ { efEDI, "-o", "sam", ffWRITE }
+ };
+#define NFILE asize(fnm)
+ edi_params.outfrq=100; edi_params.slope=0.0; edi_params.maxedsteps=0;
+ CopyRight(stderr,argv[0]);
+ parse_common_args(&argc,argv, 0 ,
+ NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
+
+ indexfile=ftp2fn_null(efNDX,NFILE,fnm);
+ EdiFile=ftp2fn(efEDI,NFILE,fnm);
+ TargetFile = opt2fn_null("-tar",NFILE,fnm);
+ OriginFile = opt2fn_null("-ori",NFILE,fnm);
+
+
+ for (ev_class=0; ev_class<evNr; ++ev_class) {
+ if (opt2parg_bSet(evOptions[ev_class],NPA,pa)) {
+ /*get list of eigenvectors*/
+ nvecs=sscan_list(&(listen[ev_class]),opt2parg_str(evOptions[ev_class],NPA,pa),evOptions[ev_class]);
+ if (ev_class<evStepNr-2) {
+ /*if apropriate get list of stepsizes for these eigenvectors*/
+ if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
+ evStepList[ev_class]=
+ scan_vecparams(opt2parg_str(evStepOptions[ev_class],NPA,pa),evStepOptions[ev_class],nvecs);
+ else { /*if list is not given fill with zeros */
+ snew(evStepList[ev_class],nvecs);
+ for (i=0; i<nvecs; i++)
+ evStepList[ev_class][i]=0.0;
+ }
+ } else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
+ snew(evStepList[ev_class],nvecs);
+ for (i=0; i<nvecs; i++)
+ evStepList[ev_class][i]=radfix;
+ } else if (ev_class == evFLOOD) {
+ snew(evStepList[ev_class],nvecs);
+
+ /* Are we doing constant force flooding? In that case, we read in
+ * the fproj values from the command line */
+ if (opt2parg_bSet("-constF", NPA, pa))
+ {
+ evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
+ }
+ } else {}; /*to avoid ambiguity */
+ } else { /* if there are no eigenvectors for this option set list to zero */
+ listen[ev_class]=NULL;
+ snew(listen[ev_class],1);
+ listen[ev_class][0]=0;
+ }
+ }
+
+ /* print the interpreted list of eigenvectors - to give some feedback*/
+ for (ev_class=0; ev_class<evNr; ++ev_class) {
+ printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
+ i=0;
+ while(listen[ev_class][i])
+ printf("%d ",listen[ev_class][i++]);
+ printf("\n");
+ }
+
+ EigvecFile=NULL;
+ EigvecFile=opt2fn("-f",NFILE,fnm);
+
+ /*read eigenvectors from eigvec.trr*/
+ read_eigenvectors(EigvecFile,&nav,&bFit1,
+ &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);
+
+ bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
+ title,&top,&ePBC,&xtop,NULL,topbox,0);
+ atoms=&top.atoms;
+
+
+ printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",nav);
+ get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
+ if (i!=nav) {
+ gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
+ i,nav);
+ }
+ printf("\n");
+
+
+ if (xref1==NULL)
+ {
+ if (bFit1)
+ {
+ /* if g_covar used different coordinate groups to fit and to do the PCA */
+ printf("\nNote: the structure in %s should be the same\n"
+ " as the one used for the fit in g_covar\n",ftp2fn(efTPS,NFILE,fnm));
+ printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
+ }
+ else
+ {
+ printf("\nNote: Apparently no fitting was done in g_covar.\n"
+ " However, you need to select a reference group for fitting in mdrun\n");
+ }
+ get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
+ snew(xref1,nfit);
+ for (i=0;i<nfit;i++)
+ copy_rvec(xtop[ifit[i]],xref1[i]);
+ }
+ else
+ {
+ nfit=nav;
+ ifit=index;
+ }
+
+ if (opt2parg_bSet("-constF", NPA, pa))
+ {
+ /* Constant force flooding is special: Most of the normal flooding
+ * options are not needed. */
+ edi_params.flood.bConstForce = TRUE;
+ }
+ else
+ {
+ /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
+
+ if ( listen[evFLOOD][0]!=0)
+ read_eigenvalues(listen[evFLOOD],opt2fn("-eig",NFILE,fnm),evStepList[evFLOOD],bHesse,kB*T);
+
+ edi_params.flood.tau=tau;
+ edi_params.flood.deltaF0=deltaF0;
+ edi_params.flood.deltaF=deltaF;
+ edi_params.presteps=eqSteps;
+ edi_params.flood.kT=kB*T;
+ edi_params.flood.bHarmonic=bHarmonic;
+ if (bRestrain)
+ {
+ /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
+ edi_params.flood.constEfl=-constEfl;
+ edi_params.flood.alpha2=-sqr(alpha);
+ }
+ else
+ {
+ edi_params.flood.constEfl=constEfl;
+ edi_params.flood.alpha2=sqr(alpha);
+ }
+ }
+
+ edi_params.ned=nav;
+
+ /*number of system atoms */
+ edi_params.nini=atoms->nr;
+
+
+ /*store reference and average structure in edi_params*/
+ make_t_edx(&edi_params.sref,nfit,xref1,ifit );
+ make_t_edx(&edi_params.sav ,nav ,xav1 ,index);
+
+
+ /* Store target positions in edi_params */
+ if (opt2bSet("-tar",NFILE,fnm))
+ {
+ if (0 != listen[evFLOOD][0])
+ {
+ fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
+ " You may want to use -ori to define the flooding potential center.\n\n");
+ }
+ get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,nav,index);
+ }
+ else
+ {
+ make_t_edx(&edi_params.star,0,NULL,index);
+ }
+
+ /* Store origin positions */
+ if (opt2bSet("-ori",NFILE,fnm))
+ {
+ get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,nav,index);
+ }
+ else
+ {
+ make_t_edx(&edi_params.sori,0,NULL,index);
+ }
+
+ /* Write edi-file */
+ write_the_whole_thing(ffopen(EdiFile,"w"), &edi_params, eigvec1,nvec1, listen, evStepList);
+ thanx(stderr);
+
+ return 0;
+}
+
--- /dev/null
+/*
+ *
+ * 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 <ctype.h>
+#include "sysstuff.h"
+#include "strdb.h"
+#include "futil.h"
+#include "macros.h"
+#include "string2.h"
+#include "statutil.h"
+#include "confio.h"
+#include "copyrite.h"
+#include "typedefs.h"
+#include "index.h"
+#include "smalloc.h"
+#include "vec.h"
+#include "index.h"
+
+#define MAXNAMES 30
+#define NAME_LEN 30
+
+gmx_bool bCase=FALSE;
+
+static int or_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
+ atom_id *nr,atom_id *at)
+{
+ atom_id i1,i2,max=0;
+ gmx_bool bNotIncr;
+
+ *nr=0;
+
+ bNotIncr=FALSE;
+ for(i1=0; i1<nr1; i1++) {
+ if ((i1>0) && (at1[i1] <= max))
+ bNotIncr=TRUE;
+ max=at1[i1];
+ }
+ for(i1=0; i1<nr2; i1++) {
+ if ((i1>0) && (at2[i1] <= max))
+ bNotIncr=TRUE;
+ max=at2[i1];
+ }
+
+ if (bNotIncr)
+ printf("One of your groups is not ascending\n");
+ else {
+ i1=0;
+ i2=0;
+ *nr=0;
+ while ((i1 < nr1) || (i2 < nr2)) {
+ if ((i2 == nr2) || ((i1<nr1) && (at1[i1] < at2[i2]))) {
+ at[*nr]=at1[i1];
+ (*nr)++;
+ i1++;
+ }
+ else {
+ if ((i2<nr2) && ((i1==nr1) || (at1[i1] > at2[i2]))) {
+ at[*nr]=at2[i2];
+ (*nr)++;
+ }
+ i2++;
+ }
+ }
+
+ printf("Merged two groups with OR: %u %u -> %u\n",nr1,nr2,*nr);
+ }
+
+ return *nr;
+}
+
+static int and_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
+ atom_id *nr,atom_id *at)
+{
+ atom_id i1,i2;
+
+ *nr=0;
+ for (i1=0; i1<nr1; i1++) {
+ for (i2=0; i2<nr2; i2++) {
+ if (at1[i1]==at2[i2]) {
+ at[*nr]=at1[i1];
+ (*nr)++;
+ }
+ }
+ }
+
+ printf("Merged two groups with AND: %u %u -> %u\n",nr1,nr2,*nr);
+
+ return *nr;
+}
+
+static gmx_bool is_name_char(char c)
+{
+ /* This string should contain all characters that can not be
+ * the first letter of a name due to the make_ndx syntax.
+ */
+ const char *spec=" !&|";
+
+ return (c != '\0' && strchr(spec,c) == NULL);
+}
+
+static int parse_names(char **string,int *n_names,char **names)
+{
+ int i;
+
+ *n_names=0;
+ while (is_name_char((*string)[0]) || (*string)[0] == ' ') {
+ if (is_name_char((*string)[0])) {
+ if (*n_names >= MAXNAMES)
+ gmx_fatal(FARGS,"To many names: %d\n",*n_names+1);
+ i=0;
+ while (is_name_char((*string)[i])) {
+ names[*n_names][i]=(*string)[i];
+ i++;
+ if (i > NAME_LEN) {
+ printf("Name is too long, the maximum is %d characters\n",NAME_LEN);
+ return 0;
+ }
+ }
+ names[*n_names][i]='\0';
+ if (!bCase)
+ upstring(names[*n_names]);
+ *string += i;
+ (*n_names)++;
+ }
+ else
+ (*string)++;
+ }
+
+ return *n_names;
+}
+
+static gmx_bool parse_int_char(char **string,int *nr,char *c)
+{
+ char *orig;
+ gmx_bool bRet;
+
+ orig = *string;
+
+ while ((*string)[0]==' ')
+ (*string)++;
+
+ bRet=FALSE;
+
+ *c = ' ';
+
+ if (isdigit((*string)[0])) {
+ *nr=(*string)[0]-'0';
+ (*string)++;
+ while (isdigit((*string)[0])) {
+ *nr = (*nr)*10+(*string)[0]-'0';
+ (*string)++;
+ }
+ if (isalpha((*string)[0])) {
+ *c = (*string)[0];
+ (*string)++;
+ }
+ /* Check if there is at most one non-digit character */
+ if (!isalnum((*string)[0])) {
+ bRet = TRUE;
+ } else {
+ *string = orig;
+ }
+ }
+ else
+ *nr=NOTSET;
+
+ return bRet;
+}
+
+static gmx_bool parse_int(char **string,int *nr)
+{
+ char *orig,c;
+ gmx_bool bRet;
+
+ orig = *string;
+ bRet = parse_int_char(string,nr,&c);
+ if (bRet && c != ' ') {
+ *string = orig;
+ bRet = FALSE;
+ }
+
+ return bRet;
+}
+
+static gmx_bool isquote(char c)
+{
+ return (c == '\"');
+}
+
+static gmx_bool parse_string(char **string,int *nr, int ngrps, char **grpname)
+{
+ char *s, *sp;
+ char c;
+
+ while ((*string)[0]==' ')
+ (*string)++;
+
+ (*nr) = NOTSET;
+ if (isquote((*string)[0])) {
+ c=(*string)[0];
+ (*string)++;
+ s = strdup((*string));
+ sp = strchr(s, c);
+ if (sp!=NULL) {
+ (*string) += sp-s + 1;
+ sp[0]='\0';
+ (*nr) = find_group(s, ngrps, grpname);
+ }
+ }
+
+ return (*nr) != NOTSET;
+}
+
+static int select_atomnumbers(char **string,t_atoms *atoms,atom_id n1,
+ atom_id *nr,atom_id *index,char *gname)
+{
+ char buf[STRLEN];
+ int i,up;
+
+ *nr=0;
+ while ((*string)[0]==' ')
+ (*string)++;
+ if ((*string)[0]=='-') {
+ (*string)++;
+ parse_int(string,&up);
+ if ((n1<1) || (n1>atoms->nr) || (up<1) || (up>atoms->nr))
+ printf("Invalid atom range\n");
+ else {
+ for(i=n1-1; i<=up-1; i++) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ printf("Found %u atom%s in range %u-%d\n",*nr,(*nr==1)?"":"s",n1,up);
+ if (n1==up)
+ sprintf(buf,"a_%u",n1);
+ else
+ sprintf(buf,"a_%u-%d",n1,up);
+ strcpy(gname,buf);
+ }
+ }
+ else {
+ i=n1;
+ sprintf(gname,"a");
+ do {
+ if ((i-1>=0) && (i-1<atoms->nr)) {
+ index[*nr] = i-1;
+ (*nr)++;
+ sprintf(buf,"_%d",i);
+ strcat(gname,buf);
+ } else {
+ printf("Invalid atom number %d\n",i);
+ *nr = 0;
+ }
+ } while ((*nr!=0) && (parse_int(string,&i)));
+ }
+
+ return *nr;
+}
+
+static int select_residuenumbers(char **string,t_atoms *atoms,
+ atom_id n1,char c,
+ atom_id *nr,atom_id *index,char *gname)
+{
+ char buf[STRLEN];
+ int i,j,up;
+ t_resinfo *ri;
+
+ *nr=0;
+ while ((*string)[0]==' ')
+ (*string)++;
+ if ((*string)[0]=='-') {
+ /* Residue number range selection */
+ if (c != ' ') {
+ printf("Error: residue insertion codes can not be used with residue range selection\n");
+ return 0;
+ }
+ (*string)++;
+ parse_int(string,&up);
+
+ for(i=0; i<atoms->nr; i++) {
+ ri = &atoms->resinfo[atoms->atom[i].resind];
+ for(j=n1; (j<=up); j++) {
+ if (ri->nr == j && (c == ' ' || ri->ic == c)) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ }
+ printf("Found %u atom%s with res.nr. in range %u-%d\n",
+ *nr,(*nr==1)?"":"s",n1,up);
+ if (n1==up)
+ sprintf(buf,"r_%u",n1);
+ else
+ sprintf(buf,"r_%u-%d",n1,up);
+ strcpy(gname,buf);
+ }
+ else {
+ /* Individual residue number/insertion code selection */
+ j=n1;
+ sprintf(gname,"r");
+ do {
+ for(i=0; i<atoms->nr; i++) {
+ ri = &atoms->resinfo[atoms->atom[i].resind];
+ if (ri->nr == j && ri->ic == c) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ sprintf(buf,"_%d",j);
+ strcat(gname,buf);
+ } while (parse_int_char(string,&j,&c));
+ }
+
+ return *nr;
+}
+
+static int select_residueindices(char **string,t_atoms *atoms,
+ atom_id n1,char c,
+ atom_id *nr,atom_id *index,char *gname)
+{
+ /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
+ /*resind+1 for 1-indexing*/
+ char buf[STRLEN];
+ int i,j,up;
+ t_resinfo *ri;
+
+ *nr=0;
+ while ((*string)[0]==' ')
+ (*string)++;
+ if ((*string)[0]=='-') {
+ /* Residue number range selection */
+ if (c != ' ') {
+ printf("Error: residue insertion codes can not be used with residue range selection\n");
+ return 0;
+ }
+ (*string)++;
+ parse_int(string,&up);
+
+ for(i=0; i<atoms->nr; i++) {
+ ri = &atoms->resinfo[atoms->atom[i].resind];
+ for(j=n1; (j<=up); j++) {
+ if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c)) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ }
+ printf("Found %u atom%s with resind.+1 in range %u-%d\n",
+ *nr,(*nr==1)?"":"s",n1,up);
+ if (n1==up)
+ sprintf(buf,"r_%u",n1);
+ else
+ sprintf(buf,"r_%u-%d",n1,up);
+ strcpy(gname,buf);
+ }
+ else {
+ /* Individual residue number/insertion code selection */
+ j=n1;
+ sprintf(gname,"r");
+ do {
+ for(i=0; i<atoms->nr; i++) {
+ ri = &atoms->resinfo[atoms->atom[i].resind];
+ if (atoms->atom[i].resind+1 == j && ri->ic == c) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ sprintf(buf,"_%d",j);
+ strcat(gname,buf);
+ } while (parse_int_char(string,&j,&c));
+ }
+
+ return *nr;
+}
+
+
+static gmx_bool atoms_from_residuenumbers(t_atoms *atoms,int group,t_blocka *block,
+ atom_id *nr,atom_id *index,char *gname)
+{
+ int i,j,j0,j1,resnr,nres;
+
+ j0=block->index[group];
+ j1=block->index[group+1];
+ nres = atoms->nres;
+ for(j=j0; j<j1; j++)
+ if (block->a[j]>=nres) {
+ printf("Index %s contains number>nres (%d>%d)\n",
+ gname,block->a[j]+1,nres);
+ return FALSE;
+ }
+ for(i=0; i<atoms->nr; i++) {
+ resnr = atoms->resinfo[atoms->atom[i].resind].nr;
+ for (j=j0; j<j1; j++)
+ if (block->a[j]+1 == resnr) {
+ index[*nr]=i;
+ (*nr)++;
+ break;
+ }
+ }
+ printf("Found %u atom%s in %d residues from group %s\n",
+ *nr,(*nr==1)?"":"s",j1-j0,gname);
+ return *nr;
+}
+
+static gmx_bool comp_name(char *name,char *search)
+{
+ while (name[0] != '\0' && search[0] != '\0') {
+ switch (search[0]) {
+ case '?':
+ /* Always matches */
+ break;
+ case '*':
+ if (search[1] != '\0') {
+ printf("WARNING: Currently '*' is only supported at the end of an expression\n");
+ return FALSE;
+ } else {
+ return TRUE;
+ }
+ break;
+ default:
+ /* Compare a single character */
+ if (( bCase && strncmp(name,search,1)) ||
+ (!bCase && gmx_strncasecmp(name,search,1))) {
+ return FALSE;
+ }
+ }
+ name++;
+ search++;
+ }
+
+ return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
+}
+
+static int select_chainnames(t_atoms *atoms,int n_names,char **names,
+ atom_id *nr,atom_id *index)
+{
+ char name[2];
+ int j;
+ atom_id i;
+
+ name[1]=0;
+ *nr=0;
+ for(i=0; i<atoms->nr; i++) {
+ name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
+ j=0;
+ while (j<n_names && !comp_name(name,names[j]))
+ j++;
+ if (j<n_names) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ printf("Found %u atom%s with chain identifier%s",
+ *nr,(*nr==1)?"":"s",(n_names==1)?"":"s");
+ for(j=0; (j<n_names); j++)
+ printf(" %s",names[j]);
+ printf("\n");
+
+ return *nr;
+}
+
+static int select_atomnames(t_atoms *atoms,int n_names,char **names,
+ atom_id *nr,atom_id *index,gmx_bool bType)
+{
+ char *name;
+ int j;
+ atom_id i;
+
+ *nr=0;
+ for(i=0; i<atoms->nr; i++) {
+ if (bType)
+ name=*(atoms->atomtype[i]);
+ else
+ name=*(atoms->atomname[i]);
+ j=0;
+ while (j<n_names && !comp_name(name,names[j]))
+ j++;
+ if (j<n_names) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ printf("Found %u atoms with %s%s",
+ *nr,bType ? "type" : "name",(n_names==1)?"":"s");
+ for(j=0; (j<n_names); j++)
+ printf(" %s",names[j]);
+ printf("\n");
+
+ return *nr;
+}
+
+static int select_residuenames(t_atoms *atoms,int n_names,char **names,
+ atom_id *nr,atom_id *index)
+{
+ char *name;
+ int j;
+ atom_id i;
+
+ *nr=0;
+ for(i=0; i<atoms->nr; i++) {
+ name = *(atoms->resinfo[atoms->atom[i].resind].name);
+ j=0;
+ while (j<n_names && !comp_name(name,names[j]))
+ j++;
+ if (j<n_names) {
+ index[*nr]=i;
+ (*nr)++;
+ }
+ }
+ printf("Found %u atoms with residue name%s",*nr,(n_names==1)?"":"s");
+ for(j=0; (j<n_names); j++)
+ printf(" %s",names[j]);
+ printf("\n");
+
+ return *nr;
+}
+
+static void copy2block(int n,atom_id *index,t_blocka *block)
+{
+ int i,n0;
+
+ block->nr++;
+ n0=block->nra;
+ block->nra=n0+n;
+ srenew(block->index,block->nr+1);
+ block->index[block->nr]=n0+n;
+ srenew(block->a,n0+n);
+ for(i=0; (i<n); i++)
+ block->a[n0+i]=index[i];
+}
+
+static void make_gname(int n,char **names,char *gname)
+{
+ int i;
+
+ strcpy(gname,names[0]);
+ for (i=1; i<n; i++) {
+ strcat(gname,"_");
+ strcat(gname,names[i]);
+ }
+}
+
+static void copy_group(int g,t_blocka *block,atom_id *nr,atom_id *index)
+{
+ int i,i0;
+
+ i0=block->index[g];
+ *nr=block->index[g+1]-i0;
+ for (i=0; i<*nr; i++)
+ index[i]=block->a[i0+i];
+}
+
+static void remove_group(int nr,int nr2,t_blocka *block,char ***gn)
+{
+ int i,j,shift;
+ char *name;
+
+ if (nr2==NOTSET)
+ nr2=nr;
+
+ for(j=0; j<=nr2-nr; j++) {
+ if ((nr<0) || (nr>=block->nr))
+ printf("Group %d does not exist\n",nr+j);
+ else {
+ shift=block->index[nr+1]-block->index[nr];
+ for(i=block->index[nr+1]; i<block->nra; i++)
+ block->a[i-shift]=block->a[i];
+
+ for(i=nr; i<block->nr; i++)
+ block->index[i]=block->index[i+1]-shift;
+ name = strdup((*gn)[nr]);
+ sfree((*gn)[nr]);
+ for(i=nr; i<block->nr-1; i++) {
+ (*gn)[i]=(*gn)[i+1];
+ }
+ block->nr--;
+ block->nra=block->index[block->nr];
+ printf("Removed group %d '%s'\n",nr+j,name);
+ sfree(name);
+ }
+ }
+}
+
+static void split_group(t_atoms *atoms,int sel_nr,t_blocka *block,char ***gn,
+ gmx_bool bAtom)
+{
+ char buf[STRLEN],*name;
+ int i,resind;
+ atom_id a,n0,n1;
+
+ printf("Splitting group %d '%s' into %s\n",sel_nr,(*gn)[sel_nr],
+ bAtom ? "atoms" : "residues");
+
+ n0=block->index[sel_nr];
+ n1=block->index[sel_nr+1];
+ srenew(block->a,block->nra+n1-n0);
+ for (i=n0; i<n1; i++) {
+ a=block->a[i];
+ resind = atoms->atom[a].resind;
+ name = *(atoms->resinfo[resind].name);
+ if (bAtom || (i==n0) || (atoms->atom[block->a[i-1]].resind!=resind)) {
+ if (i>n0)
+ block->index[block->nr]=block->nra;
+ block->nr++;
+ srenew(block->index,block->nr+1);
+ srenew(*gn,block->nr);
+ if (bAtom)
+ sprintf(buf,"%s_%s_%u",(*gn)[sel_nr],*atoms->atomname[a],a+1);
+ else
+ sprintf(buf,"%s_%s_%d",(*gn)[sel_nr],name,atoms->resinfo[resind].nr);
+ (*gn)[block->nr-1]=strdup(buf);
+ }
+ block->a[block->nra]=a;
+ block->nra++;
+ }
+ block->index[block->nr]=block->nra;
+}
+
+static int split_chain(t_atoms *atoms,rvec *x,
+ int sel_nr,t_blocka *block,char ***gn)
+{
+ char buf[STRLEN];
+ int j,nchain;
+ atom_id i,a,natoms,*start=NULL,*end=NULL,ca_start,ca_end;
+ rvec vec;
+
+ natoms=atoms->nr;
+ nchain=0;
+ ca_start=0;
+
+ while (ca_start<natoms) {
+ while((ca_start<natoms) && strcmp(*atoms->atomname[ca_start],"CA"))
+ ca_start++;
+ if (ca_start<natoms) {
+ srenew(start,nchain+1);
+ srenew(end,nchain+1);
+ start[nchain]=ca_start;
+ while ((start[nchain]>0) &&
+ (atoms->atom[start[nchain]-1].resind ==
+ atoms->atom[ca_start].resind))
+ start[nchain]--;
+
+ i=ca_start;
+ do {
+ ca_end=i;
+ do {
+ i++;
+ } while ((i<natoms) && strcmp(*atoms->atomname[i],"CA"));
+ if (i<natoms)
+ rvec_sub(x[ca_end],x[i],vec);
+ } while ((i<natoms) && (norm(vec)<0.45));
+
+ end[nchain]=ca_end;
+ while ((end[nchain]+1<natoms) &&
+ (atoms->atom[end[nchain]+1].resind==atoms->atom[ca_end].resind))
+ end[nchain]++;
+ ca_start=end[nchain]+1;
+ nchain++;
+ }
+ }
+ if (nchain==1)
+ printf("Found 1 chain, will not split\n");
+ else
+ printf("Found %d chains\n",nchain);
+ for (j=0; j<nchain; j++)
+ printf("%d:%6u atoms (%u to %u)\n",
+ j+1,end[j]-start[j]+1,start[j]+1,end[j]+1);
+
+ if (nchain > 1) {
+ srenew(block->a,block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
+ for (j=0; j<nchain; j++) {
+ block->nr++;
+ srenew(block->index,block->nr+1);
+ srenew(*gn,block->nr);
+ sprintf(buf,"%s_chain%d",(*gn)[sel_nr],j+1);
+ (*gn)[block->nr-1]=strdup(buf);
+ for (i=block->index[sel_nr]; i<block->index[sel_nr+1]; i++) {
+ a=block->a[i];
+ if ((a>=start[j]) && (a<=end[j])) {
+ block->a[block->nra]=a;
+ block->nra++;
+ }
+ }
+ block->index[block->nr]=block->nra;
+ if (block->index[block->nr-1]==block->index[block->nr])
+ remove_group(block->nr-1,NOTSET,block,gn);
+ }
+ }
+ sfree(start);
+ sfree(end);
+
+ return nchain;
+}
+
+static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
+{
+ if ( atoms==NULL ) {
+ printf("Can not process '%s' without atom info, use option -f\n", string);
+ return FALSE;
+ } else
+ return TRUE;
+}
+
+static gmx_bool parse_entry(char **string,int natoms,t_atoms *atoms,
+ t_blocka *block,char ***gn,
+ atom_id *nr,atom_id *index,char *gname)
+{
+ static char **names, *ostring;
+ static gmx_bool bFirst=TRUE;
+ int j,n_names,sel_nr1;
+ atom_id i,nr1,*index1;
+ char c;
+ gmx_bool bRet,bCompl;
+
+ if (bFirst) {
+ bFirst=FALSE;
+ snew(names,MAXNAMES);
+ for (i=0; i<MAXNAMES; i++)
+ snew(names[i],NAME_LEN+1);
+ }
+
+ bRet=FALSE;
+ sel_nr1=NOTSET;
+
+ while(*string[0]==' ')
+ (*string)++;
+
+ if ((*string)[0]=='!') {
+ bCompl=TRUE;
+ (*string)++;
+ while(*string[0]==' ')
+ (*string)++;
+ } else
+ bCompl=FALSE;
+
+ ostring = *string;
+
+ if (parse_int(string,&sel_nr1) ||
+ parse_string(string,&sel_nr1,block->nr,*gn)) {
+ if ((sel_nr1>=0) && (sel_nr1<block->nr)) {
+ copy_group(sel_nr1,block,nr,index);
+ strcpy(gname,(*gn)[sel_nr1]);
+ printf("Copied index group %d '%s'\n",sel_nr1,(*gn)[sel_nr1]);
+ bRet=TRUE;
+ } else
+ printf("Group %d does not exist\n",sel_nr1);
+ }
+ else if ((*string)[0]=='a') {
+ (*string)++;
+ if (check_have_atoms(atoms, ostring)) {
+ if (parse_int(string,&sel_nr1)) {
+ bRet=select_atomnumbers(string,atoms,sel_nr1,nr,index,gname);
+ }
+ else if (parse_names(string,&n_names,names)) {
+ bRet=select_atomnames(atoms,n_names,names,nr,index,FALSE);
+ make_gname(n_names,names,gname);
+ }
+ }
+ }
+ else if ((*string)[0]=='t') {
+ (*string)++;
+ if (check_have_atoms(atoms, ostring) &&
+ parse_names(string,&n_names,names)) {
+ if (atoms->atomtype == NULL)
+ printf("Need a run input file to select atom types\n");
+ else {
+ bRet=select_atomnames(atoms,n_names,names,nr,index,TRUE);
+ make_gname(n_names,names,gname);
+ }
+ }
+ }
+ else if (strncmp(*string,"res",3)==0) {
+ (*string)+=3;
+ if ( check_have_atoms(atoms, ostring) &&
+ parse_int(string,&sel_nr1) &&
+ (sel_nr1>=0) && (sel_nr1<block->nr) ) {
+ bRet=atoms_from_residuenumbers(atoms,
+ sel_nr1,block,nr,index,(*gn)[sel_nr1]);
+ sprintf(gname,"atom_%s",(*gn)[sel_nr1]);
+ }
+ }
+ else if (strncmp(*string,"ri",2)==0) {
+ (*string)+=2;
+ if (check_have_atoms(atoms, ostring) &&
+ parse_int_char(string,&sel_nr1,&c)) {
+ bRet=select_residueindices(string,atoms,sel_nr1,c,nr,index,gname);
+ }
+ }
+ else if ((*string)[0]=='r') {
+ (*string)++;
+ if (check_have_atoms(atoms, ostring)) {
+ if (parse_int_char(string,&sel_nr1,&c)) {
+ bRet=select_residuenumbers(string,atoms,sel_nr1,c,nr,index,gname);
+ }
+ else if (parse_names(string,&n_names,names)) {
+ bRet=select_residuenames(atoms,n_names,names,nr,index);
+ make_gname(n_names,names,gname);
+ }
+ }
+ }
+ else if (strncmp(*string,"chain",5)==0) {
+ (*string)+=5;
+ if (check_have_atoms(atoms, ostring) &&
+ parse_names(string,&n_names,names)) {
+ bRet=select_chainnames(atoms,n_names,names,nr,index);
+ sprintf(gname,"ch%s",names[0]);
+ for (i=1; i<n_names; i++)
+ strcat(gname,names[i]);
+ }
+ }
+ if (bRet && bCompl) {
+ snew(index1,natoms-*nr);
+ nr1=0;
+ for(i=0; i<natoms; i++) {
+ j=0;
+ while ((j<*nr) && (index[j] != i))
+ j++;
+ if (j==*nr) {
+ if (nr1 >= natoms-*nr) {
+ printf("There are double atoms in your index group\n");
+ break;
+ }
+ index1[nr1]=i;
+ nr1++;
+ }
+ }
+ *nr=nr1;
+ for(i=0; i<nr1; i++)
+ index[i]=index1[i];
+ sfree(index1);
+
+ for (i=strlen(gname)+1; i>0; i--)
+ gname[i]=gname[i-1];
+ gname[0]='!';
+ printf("Complemented group: %u atoms\n",*nr);
+ }
+
+ return bRet;
+}
+
+static void list_residues(t_atoms *atoms)
+{
+ int i,j,start,end,prev_resind,resind;
+ gmx_bool bDiff;
+
+ /* Print all the residues, assuming continuous resnr count */
+ start = atoms->atom[0].resind;
+ prev_resind = start;
+ for(i=0; i<atoms->nr; i++) {
+ resind = atoms->atom[i].resind;
+ if ((resind != prev_resind) || (i==atoms->nr-1)) {
+ if ((bDiff=strcmp(*atoms->resinfo[resind].name,
+ *atoms->resinfo[start].name)) ||
+ (i==atoms->nr-1)) {
+ if (bDiff)
+ end = prev_resind;
+ else
+ end = resind;
+ if (end < start+3)
+ for(j=start; j<=end; j++)
+ printf("%4d %-5s",
+ j+1,*(atoms->resinfo[j].name));
+ else
+ printf(" %4d - %4d %-5s ",
+ start+1,end+1,*(atoms->resinfo[start].name));
+ start = resind;
+ }
+ }
+ prev_resind = resind;
+ }
+ printf("\n");
+}
+
+static void edit_index(int natoms, t_atoms *atoms,rvec *x,t_blocka *block, char ***gn, gmx_bool bVerbose)
+{
+ static char **atnames, *ostring;
+ static gmx_bool bFirst=TRUE;
+ char inp_string[STRLEN],*string;
+ char gname[STRLEN],gname1[STRLEN],gname2[STRLEN];
+ int i,i0,i1,sel_nr,sel_nr2,newgroup;
+ atom_id nr,nr1,nr2,*index,*index1,*index2;
+ gmx_bool bAnd,bOr,bPrintOnce;
+
+ if (bFirst) {
+ bFirst=FALSE;
+ snew(atnames,MAXNAMES);
+ for (i=0; i<MAXNAMES; i++)
+ snew(atnames[i],NAME_LEN+1);
+ }
+
+ string=NULL;
+
+ snew(index,natoms);
+ snew(index1,natoms);
+ snew(index2,natoms);
+
+ newgroup=NOTSET;
+ bPrintOnce = TRUE;
+ do {
+ gname1[0]='\0';
+ if (bVerbose || bPrintOnce || newgroup!=NOTSET) {
+ printf("\n");
+ if (bVerbose || bPrintOnce || newgroup==NOTSET) {
+ i0=0;
+ i1=block->nr;
+ } else {
+ i0=newgroup;
+ i1=newgroup+1;
+ }
+ for(i=i0; i<i1; i++)
+ printf("%3d %-20s: %5u atoms\n",i,(*gn)[i],
+ block->index[i+1]-block->index[i]);
+ newgroup=NOTSET;
+ }
+ if (bVerbose || bPrintOnce) {
+ printf("\n");
+ printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
+ printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
+ printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
+ printf(" 'r': residue 'res' nr 'chain' char\n");
+ printf(" \"name\": group 'case': case %s 'q': save and quit\n",
+ bCase ? "insensitive" : "sensitive ");
+ printf(" 'ri': residue index\n");
+ bPrintOnce = FALSE;
+ }
+ printf("\n");
+ printf("> ");
+ if(NULL==fgets(inp_string,STRLEN,stdin))
+ {
+ gmx_fatal(FARGS,"Error reading user input");
+ }
+ inp_string[strlen(inp_string)-1]=0;
+ printf("\n");
+ string=inp_string;
+ while (string[0]==' ')
+ string++;
+
+ ostring = string;
+ nr=0;
+ if (string[0] == 'h') {
+ printf(" nr : selects an index group by number or quoted string.\n");
+ printf(" The string is first matched against the whole group name,\n");
+ printf(" then against the beginning and finally against an\n");
+ printf(" arbitrary substring. A multiple match is an error.\n");
+
+ printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
+ printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
+ printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
+ printf(" wildcard '*' allowed at the end of a name.\n");
+ printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
+ printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
+ printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
+ printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
+ printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
+ printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
+ printf(" not available with a .gro file as input.\n");
+ printf(" ! : takes the complement of a group with respect to all\n");
+ printf(" the atoms in the input file.\n");
+ printf(" & | : AND and OR, can be placed between any of the options\n");
+ printf(" above, the input is processed from left to right.\n");
+ printf(" 'name' nr name : rename group nr to name.\n");
+ printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
+ printf(" 'keep' nr : deletes all groups except nr.\n");
+ printf(" 'case' : make all name compares case (in)sensitive.\n");
+ printf(" 'splitch' nr : split group into chains using CA distances.\n");
+ printf(" 'splitres' nr : split group into residues.\n");
+ printf(" 'splitat' nr : split group into atoms.\n");
+ printf(" 'res' nr : interpret numbers in group as residue numbers\n");
+ printf(" Enter : list the currently defined groups and commands\n");
+ printf(" 'l' : list the residues.\n");
+ printf(" 'h' : show this help.\n");
+ printf(" 'q' : save and quit.\n");
+ printf("\n");
+ printf(" Examples:\n");
+ printf(" > 2 | 4 & r 3-5\n");
+ printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
+ printf(" > a C* & !a C CA\n");
+ printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
+ printf(" > \"protein\" & ! \"backb\"\n");
+ printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
+ if (bVerbose) {
+ printf("\npress Enter ");
+ getchar();
+ }
+ }
+ else if (strncmp(string,"del",3)==0) {
+ string+=3;
+ if (parse_int(&string,&sel_nr)) {
+ while(string[0]==' ')
+ string++;
+ if (string[0]=='-') {
+ string++;
+ parse_int(&string,&sel_nr2);
+ } else
+ sel_nr2=NOTSET;
+ while(string[0]==' ')
+ string++;
+ if (string[0]=='\0')
+ remove_group(sel_nr,sel_nr2,block,gn);
+ else
+ printf("\nSyntax error: \"%s\"\n",string);
+ }
+ }
+ else if (strncmp(string,"keep",4)==0) {
+ string+=4;
+ if (parse_int(&string,&sel_nr)) {
+ remove_group(sel_nr+1,block->nr-1,block,gn);
+ remove_group(0,sel_nr-1,block,gn);
+ }
+ }
+ else if (strncmp(string,"name",4)==0) {
+ string+=4;
+ if (parse_int(&string,&sel_nr)) {
+ if ((sel_nr>=0) && (sel_nr<block->nr)) {
+ sscanf(string,"%s",gname);
+ sfree((*gn)[sel_nr]);
+ (*gn)[sel_nr]=strdup(gname);
+ }
+ }
+ }
+ else if (strncmp(string,"case",4)==0) {
+ bCase=!bCase;
+ printf("Switched to case %s\n",bCase ? "sensitive" : "insensitive");
+ }
+ else if (string[0] == 'v' ) {
+ bVerbose=!bVerbose;
+ printf("Turned verbose %s\n",bVerbose ? "on" : "off");
+ }
+ else if (string[0] == 'l') {
+ if ( check_have_atoms(atoms, ostring) )
+ list_residues(atoms);
+ }
+ else if (strncmp(string,"splitch",7)==0) {
+ string+=7;
+ if ( check_have_atoms(atoms, ostring) &&
+ parse_int(&string,&sel_nr) &&
+ (sel_nr>=0) && (sel_nr<block->nr))
+ split_chain(atoms,x,sel_nr,block,gn);
+ }
+ else if (strncmp(string,"splitres",8)==0 ) {
+ string+=8;
+ if ( check_have_atoms(atoms, ostring) &&
+ parse_int(&string,&sel_nr) &&
+ (sel_nr>=0) && (sel_nr<block->nr))
+ split_group(atoms,sel_nr,block,gn,FALSE);
+ }
+ else if (strncmp(string,"splitat",7)==0 ) {
+ string+=7;
+ if ( check_have_atoms(atoms, ostring) &&
+ parse_int(&string,&sel_nr) &&
+ (sel_nr>=0) && (sel_nr<block->nr))
+ split_group(atoms,sel_nr,block,gn,TRUE);
+ }
+ else if (string[0] == '\0') {
+ bPrintOnce = TRUE;
+ }
+ else if (string[0] != 'q') {
+ nr1=-1;
+ nr2=-1;
+ if (parse_entry(&string,natoms,atoms,block,gn,&nr,index,gname)) {
+ do {
+ while (string[0]==' ')
+ string++;
+
+ bAnd=FALSE;
+ bOr=FALSE;
+ if (string[0]=='&')
+ bAnd=TRUE;
+ else if (string[0]=='|')
+ bOr=TRUE;
+
+ if (bAnd || bOr) {
+ string++;
+ nr1=nr;
+ for(i=0; i<nr; i++)
+ index1[i]=index[i];
+ strcpy(gname1,gname);
+ if (parse_entry(&string,natoms,atoms,block,gn,&nr2,index2,gname2)) {
+ if (bOr) {
+ or_groups(nr1,index1,nr2,index2,&nr,index);
+ sprintf(gname,"%s_%s",gname1,gname2);
+ }
+ else {
+ and_groups(nr1,index1,nr2,index2,&nr,index);
+ sprintf(gname,"%s_&_%s",gname1,gname2);
+ }
+ }
+ }
+ } while (bAnd || bOr);
+ }
+ while(string[0]==' ')
+ string++;
+ if (string[0])
+ printf("\nSyntax error: \"%s\"\n",string);
+ else if (nr>0) {
+ copy2block(nr,index,block);
+ srenew(*gn,block->nr);
+ newgroup = block->nr-1;
+ (*gn)[newgroup]=strdup(gname);
+ }
+ else
+ printf("Group is empty\n");
+ }
+ } while (string[0]!='q');
+
+ sfree(index);
+ sfree(index1);
+ sfree(index2);
+}
+
+static int block2natoms(t_blocka *block)
+{
+ int i, natoms;
+
+ natoms = 0;
+ for(i=0; i<block->nra; i++)
+ natoms = max(natoms, block->a[i]);
+ natoms++;
+
+ return natoms;
+}
+
+void merge_blocks(t_blocka *dest, t_blocka *source)
+{
+ int i,nra0,i0;
+
+ /* count groups, srenew and fill */
+ i0 = dest->nr;
+ nra0 = dest->nra;
+ dest->nr+=source->nr;
+ srenew(dest->index, dest->nr+1);
+ for(i=0; i<source->nr; i++)
+ dest->index[i0+i] = nra0 + source->index[i];
+ /* count atoms, srenew and fill */
+ dest->nra+=source->nra;
+ srenew(dest->a, dest->nra);
+ for(i=0; i<source->nra; i++)
+ dest->a[nra0+i] = source->a[i];
+
+ /* terminate list */
+ dest->index[dest->nr]=dest->nra;
+
+}
+
+int gmx_make_ndx(int argc,char *argv[])
+{
+ const char *desc[] = {
+ "Index groups are necessary for almost every gromacs program.",
+ "All these programs can generate default index groups. You ONLY",
+ "have to use [TT]make_ndx[tt] when you need SPECIAL index groups.",
+ "There is a default index group for the whole system, 9 default",
+ "index groups for proteins, and a default index group",
+ "is generated for every other residue name.[PAR]",
+ "When no index file is supplied, also [TT]make_ndx[tt] will generate the",
+ "default groups.",
+ "With the index editor you can select on atom, residue and chain names",
+ "and numbers.",
+ "When a run input file is supplied you can also select on atom type.",
+ "You can use NOT, AND and OR, you can split groups",
+ "into chains, residues or atoms. You can delete and rename groups.[PAR]",
+ "The atom numbering in the editor and the index file starts at 1."
+ };
+
+ static int natoms=0;
+ static gmx_bool bVerbose=FALSE;
+ t_pargs pa[] = {
+ { "-natoms", FALSE, etINT, {&natoms},
+ "set number of atoms (default: read from coordinate or index file)" },
+ { "-verbose", FALSE, etBOOL, {&bVerbose},
+ "HIDDENVerbose output" }
+ };
+#define NPA asize(pa)
+
+ output_env_t oenv;
+ char title[STRLEN];
+ int nndxin;
+ const char *stxfile;
+ char **ndxinfiles;
+ const char *ndxoutfile;
+ gmx_bool bNatoms;
+ int i,j;
+ t_atoms *atoms;
+ rvec *x,*v;
+ int ePBC;
+ matrix box;
+ t_blocka *block,*block2;
+ char **gnames,**gnames2;
+ t_filenm fnm[] = {
+ { efSTX, "-f", NULL, ffOPTRD },
+ { efNDX, "-n", NULL, ffOPTRDMULT },
+ { efNDX, "-o", NULL, ffWRITE }
+ };
+#define NFILE asize(fnm)
+
+ CopyRight(stderr,argv[0]);
+
+ parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,asize(desc),desc,
+ 0,NULL,&oenv);
+
+ stxfile = ftp2fn_null(efSTX,NFILE,fnm);
+ if (opt2bSet("-n",NFILE,fnm)) {
+ nndxin = opt2fns(&ndxinfiles,"-n",NFILE,fnm);
+ } else {
+ nndxin = 0;
+ }
+ ndxoutfile = opt2fn("-o",NFILE,fnm);
+ bNatoms = opt2parg_bSet("-natoms",NPA,pa);
+
+ if (!stxfile && !nndxin)
+ gmx_fatal(FARGS,"No input files (structure or index)");
+
+ if (stxfile) {
+ snew(atoms,1);
+ get_stx_coordnum(stxfile,&(atoms->nr));
+ init_t_atoms(atoms,atoms->nr,TRUE);
+ snew(x,atoms->nr);
+ snew(v,atoms->nr);
+ fprintf(stderr,"\nReading structure file\n");
+ read_stx_conf(stxfile,title,atoms,x,v,&ePBC,box);
+ natoms = atoms->nr;
+ bNatoms=TRUE;
+ } else {
+ atoms = NULL;
+ x = NULL;
+ }
+
+ /* read input file(s) */
+ block = new_blocka();
+ gnames = NULL;
+ printf("Going to read %d old index file(s)\n",nndxin);
+ if (nndxin) {
+ for(i=0; i<nndxin; i++) {
+ block2 = init_index(ndxinfiles[i],&gnames2);
+ srenew(gnames, block->nr+block2->nr);
+ for(j=0; j<block2->nr; j++)
+ gnames[block->nr+j]=gnames2[j];
+ sfree(gnames2);
+ merge_blocks(block, block2);
+ sfree(block2->a);
+ sfree(block2->index);
+/* done_block(block2); */
+ sfree(block2);
+ }
+ }
+ else {
+ snew(gnames,1);
+ analyse(atoms,block,&gnames,FALSE,TRUE);
+ }
+
+ if (!bNatoms) {
+ natoms = block2natoms(block);
+ printf("Counted atom numbers up to %d in index file\n", natoms);
+ }
+
+ edit_index(natoms,atoms,x,block,&gnames,bVerbose);
+
+ write_index(ndxoutfile,block,gnames);
+
+ thanx(stderr);
+
+ return 0;
+}
--- /dev/null
+/*
+ *
+ * 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 "typedefs.h"
+#include "smalloc.h"
+#include "copyrite.h"
+#include "statutil.h"
+#include "macros.h"
+#include "string2.h"
+#include "futil.h"
+#include "gmx_fatal.h"
+
+static int calc_ntype(int nft,int *ft,t_idef *idef)
+{
+ int i,f,nf=0;
+
+ for(i=0; (i<idef->ntypes); i++) {
+ for(f=0; f<nft; f++) {
+ if (idef->functype[i] == ft[f])
+ nf++;
+ }
+ }
+
+ return nf;
+}
+
+static void fill_ft_ind(int nft,int *ft,t_idef *idef,
+ int ft_ind[],char *grpnames[])
+{
+ char buf[125];
+ int i,f,ftype,ind=0;
+
+ /* Loop over all the function types in the topology */
+ for(i=0; (i<idef->ntypes); i++) {
+ ft_ind[i] = -1;
+ /* Check all the selected function types */
+ for(f=0; f<nft; f++) {
+ ftype = ft[f];
+ if (idef->functype[i] == ftype) {
+ ft_ind[i] = ind;
+ switch (ftype) {
+ case F_ANGLES:
+ sprintf(buf,"Theta=%.1f_%.2f",idef->iparams[i].harmonic.rA,
+ idef->iparams[i].harmonic.krA);
+ break;
+ case F_G96ANGLES:
+ sprintf(buf,"Cos_th=%.1f_%.2f",idef->iparams[i].harmonic.rA,
+ idef->iparams[i].harmonic.krA);
+ break;
+ case F_UREY_BRADLEY:
+ sprintf(buf,"UB_th=%.1f_%.2f2f",idef->iparams[i].u_b.thetaA,
+ idef->iparams[i].u_b.kthetaA);
+ break;
+ case F_QUARTIC_ANGLES:
+ sprintf(buf,"Q_th=%.1f_%.2f_%.2f",idef->iparams[i].qangle.theta,
+ idef->iparams[i].qangle.c[0],idef->iparams[i].qangle.c[1]);
+ break;
+ case F_TABANGLES:
+ sprintf(buf,"Table=%d_%.2f",idef->iparams[i].tab.table,
+ idef->iparams[i].tab.kA);
+ break;
+ case F_PDIHS:
+ sprintf(buf,"Phi=%.1f_%d_%.2f",idef->iparams[i].pdihs.phiA,
+ idef->iparams[i].pdihs.mult,idef->iparams[i].pdihs.cpA);
+ break;
+ case F_IDIHS:
+ sprintf(buf,"Xi=%.1f_%.2f",idef->iparams[i].harmonic.rA,
+ idef->iparams[i].harmonic.krA);
+ break;
+ case F_RBDIHS:
+ sprintf(buf,"RB-A1=%.2f",idef->iparams[i].rbdihs.rbcA[1]);
+ break;
+ default:
+ gmx_fatal(FARGS,"Unsupported function type '%s' selected",
+ interaction_function[ftype].longname);
+ }
+ grpnames[ind]=strdup(buf);
+ ind++;
+ }
+ }
+ }
+}
+
+static void fill_ang(int nft,int *ft,int fac,
+ int nr[],int *index[],int ft_ind[],t_topology *top,
+ gmx_bool bNoH,real hq)
+{
+ int f,ftype,i,j,indg,nr_fac;
+ gmx_bool bUse;
+ t_idef *idef;
+ t_atom *atom;
+ t_iatom *ia;
+
+
+ idef = &top->idef;
+ atom = top->atoms.atom;
+
+ for(f=0; f<nft; f++) {
+ ftype = ft[f];
+ ia = idef->il[ftype].iatoms;
+ for(i=0; (i<idef->il[ftype].nr); ) {
+ indg = ft_ind[ia[0]];
+ if (indg == -1)
+ gmx_incons("Routine fill_ang");
+ bUse = TRUE;
+ if (bNoH) {
+ for(j=0; j<fac; j++) {
+ if (atom[ia[1+j]].m < 1.5)
+ bUse = FALSE;
+ }
+ }
+ if (hq) {
+ for(j=0; j<fac; j++) {
+ if (atom[ia[1+j]].m < 1.5 && fabs(atom[ia[1+j]].q) < hq)
+ bUse = FALSE;
+ }
+ }
+ if (bUse) {
+ if (nr[indg] % 1000 == 0) {
+ srenew(index[indg],fac*(nr[indg]+1000));
+ }
+ nr_fac = fac*nr[indg];
+ for(j=0; (j<fac); j++)
+ index[indg][nr_fac+j] = ia[j+1];
+ nr[indg]++;
+ }
+ ia += interaction_function[ftype].nratoms+1;
+ i += interaction_function[ftype].nratoms+1;
+ }
+ }
+}
+
+static int *select_ftype(const char *opt,int *nft,int *mult)
+{
+ int *ft=NULL,ftype;
+
+ if (opt[0] == 'a') {
+ *mult = 3;
+ for(ftype=0; ftype<F_NRE; ftype++) {
+ if ((interaction_function[ftype].flags & IF_ATYPE) ||
+ ftype == F_TABANGLES) {
+ (*nft)++;
+ srenew(ft,*nft);
+ ft[*nft-1] = ftype;
+ }
+ }
+ } else {
+ *mult = 4;
+ *nft = 1;
+ snew(ft,*nft);
+ switch(opt[0]) {
+ case 'd':
+ ft[0] = F_PDIHS;
+ break;
+ case 'i':
+ ft[0] = F_IDIHS;
+ break;
+ case 'r':
+ ft[0] = F_RBDIHS;
+ break;
+ default:
+ break;
+ }
+ }
+
+ return ft;
+}
+
+int gmx_mk_angndx(int argc,char *argv[])
+{
+ static const char *desc[] = {
+ "[TT]mk_angndx[tt] makes an index file for calculation of",
+ "angle distributions etc. It uses a run input file ([TT].tpx[tt]) for the",
+ "definitions of the angles, dihedrals etc."
+ };
+ static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
+ static gmx_bool bH=TRUE;
+ static real hq=-1;
+ t_pargs pa[] = {
+ { "-type", FALSE, etENUM, {opt},
+ "Type of angle" },
+ { "-hyd", FALSE, etBOOL, {&bH},
+ "Include angles with atoms with mass < 1.5" },
+ { "-hq", FALSE, etREAL, {&hq},
+ "Ignore angles with atoms with mass < 1.5 and magnitude of their charge less than this value" }
+ };
+
+ output_env_t oenv;
+ FILE *out;
+ t_topology *top;
+ int i,j,ntype;
+ int nft=0,*ft,mult=0;
+ int **index;
+ int *ft_ind;
+ int *nr;
+ char **grpnames;
+ t_filenm fnm[] = {
+ { efTPX, NULL, NULL, ffREAD },
+ { efNDX, NULL, "angle", ffWRITE }
+ };
+#define NFILE asize(fnm)
+
+ CopyRight(stderr,argv[0]);
+ parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
+ asize(desc),desc,0,NULL,&oenv);
+
+
+ ft = select_ftype(opt[0],&nft,&mult);
+
+ top = read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
+
+ ntype = calc_ntype(nft,ft,&(top->idef));
+ snew(grpnames,ntype);
+ snew(ft_ind,top->idef.ntypes);
+ fill_ft_ind(nft,ft,&top->idef,ft_ind,grpnames);
+
+ snew(nr,ntype);
+ snew(index,ntype);
+ fill_ang(nft,ft,mult,nr,index,ft_ind,top,!bH,hq);
+
+ out=ftp2FILE(efNDX,NFILE,fnm,"w");
+ for(i=0; (i<ntype); i++) {
+ if (nr[i] > 0) {
+ fprintf(out,"[ %s ]\n",grpnames[i]);
+ for(j=0; (j<nr[i]*mult); j++) {
+ fprintf(out," %5d",index[i][j]+1);
+ if ((j % 12) == 11)
+ fprintf(out,"\n");
+ }
+ fprintf(out,"\n");
+ }
+ }
+ ffclose(out);
+
+ thanx(stderr);
+
+ return 0;
+}
--- /dev/null
+/*
+ *
+ * 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:
+ * Good gRace! Old Maple Actually Chews Slate
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include "typedefs.h"
+#include "statutil.h"
+#include "copyrite.h"
+#include "gmx_fatal.h"
+#include "xvgr.h"
+#include "pdbio.h"
+#include "macros.h"
+#include "smalloc.h"
+#include "vec.h"
+#include "pbc.h"
+#include "physics.h"
+#include "names.h"
+#include "txtdump.h"
+#include "trnio.h"
+#include "symtab.h"
+#include "confio.h"
+
+real pot(real x,real qq,real c6,real cn,int npow)
+{
+ return cn*pow(x,-npow)-c6*pow(x,-6)+qq*ONE_4PI_EPS0/x;
+}
+
+real bhpot(real x,real qq,real A,real B,real C)
+{
+ return A*exp(-B*x) - C*pow(x,-6.0);
+}
+
+real dpot(real x,real qq,real c6,real cn,int npow)
+{
+ return -(npow*cn*pow(x,-npow-1)-6*c6*pow(x,-7)+qq*ONE_4PI_EPS0/sqr(x));
+}
+
+int gmx_sigeps(int argc,char *argv[])
+{
+ const char *desc[] = {
+ "[TT]g_sigeps[tt] is a simple utility that converts C6/C12 or C6/Cn combinations",
+ "to [GRK]sigma[grk] and [GRK]epsilon[grk], or vice versa. It can also plot the potential",
+ "in file. In addition, it makes an approximation of a Buckingham potential",
+ "to a Lennard-Jones potential."
+ };
+ static real c6=1.0e-3,cn=1.0e-6,qi=0,qj=0,sig=0.3,eps=1,sigfac=0.7;
+ static real Abh=1e5,Bbh=32,Cbh=1e-3;
+ static int npow=12;
+ t_pargs pa[] = {
+ { "-c6", FALSE, etREAL, {&c6}, "C6" },
+ { "-cn", FALSE, etREAL, {&cn}, "Constant for repulsion" },
+ { "-pow", FALSE, etINT, {&npow},"Power of the repulsion term" },
+ { "-sig", FALSE, etREAL, {&sig}, "[GRK]sigma[grk]" },
+ { "-eps", FALSE, etREAL, {&eps}, "[GRK]epsilon[grk]" },
+ { "-A", FALSE, etREAL, {&Abh}, "Buckingham A" },
+ { "-B", FALSE, etREAL, {&Bbh}, "Buckingham B" },
+ { "-C", FALSE, etREAL, {&Cbh}, "Buckingham C" },
+ { "-qi", FALSE, etREAL, {&qi}, "qi" },
+ { "-qj", FALSE, etREAL, {&qj}, "qj" },
+ { "-sigfac", FALSE, etREAL, {&sigfac}, "Factor in front of [GRK]sigma[grk] for starting the plot" }
+ };
+ t_filenm fnm[] = {
+ { efXVG, "-o", "potje", ffWRITE }
+ };
+ output_env_t oenv;
+#define NFILE asize(fnm)
+ const char *legend[] = { "Lennard-Jones", "Buckingham" };
+ FILE *fp;
+ int i;
+ gmx_bool bBham;
+ real qq,x,oldx,minimum,mval,dp[2],pp[2];
+ int cur=0;
+#define next (1-cur)
+
+ /* CopyRight(stdout,argv[0]);*/
+ parse_common_args(&argc,argv,PCA_CAN_VIEW,
+ NFILE,fnm,asize(pa),pa,asize(desc),
+ desc,0,NULL,&oenv);
+
+ bBham = (opt2parg_bSet("-A",asize(pa),pa) ||
+ opt2parg_bSet("-B",asize(pa),pa) ||
+ opt2parg_bSet("-C",asize(pa),pa));
+
+ if (bBham) {
+ c6 = Cbh;
+ sig = pow((6.0/npow)*pow(npow/Bbh,npow-6.0),1.0/(npow-6.0));
+ eps = c6/(4*pow(sig,6.0));
+ cn = 4*eps*pow(sig,npow);
+ }
+ else {
+ if (opt2parg_bSet("-sig",asize(pa),pa) ||
+ opt2parg_bSet("-eps",asize(pa),pa)) {
+ c6 = 4*eps*pow(sig,6);
+ cn = 4*eps*pow(sig,npow);
+ }
+ else if (opt2parg_bSet("-c6",asize(pa),pa) ||
+ opt2parg_bSet("-cn",asize(pa),pa) ||
+ opt2parg_bSet("-pow",asize(pa),pa)) {
+ sig = pow(cn/c6,1.0/(npow-6.0));
+ eps = 0.25*c6*pow(sig,-6.0);
+ }
+ else {
+ sig = eps = 0;
+ }
+ printf("c6 = %12.5e, c%d = %12.5e\n",c6,npow,cn);
+ printf("sigma = %12.5f, epsilon = %12.5f\n",sig,eps);
+
+ minimum = pow(npow/6.0*pow(sig,npow-6.0),1.0/(npow-6));
+ printf("Van der Waals minimum at %g, V = %g\n\n",
+ minimum,pot(minimum,0,c6,cn,npow));
+ printf("Fit of Lennard Jones (%d-6) to Buckingham:\n",npow);
+ Bbh = npow/minimum;
+ Cbh = c6;
+ Abh = 4*eps*pow(sig/minimum,npow)*exp(npow);
+ printf("A = %g, B = %g, C = %g\n",Abh,Bbh,Cbh);
+ }
+ qq = qi*qj;
+
+ fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"Potential","r (nm)","E (kJ/mol)",
+ oenv);
+ xvgr_legend(fp,asize(legend),legend,
+ oenv);
+ if (sig == 0)
+ sig=0.25;
+ minimum = -1;
+ mval = 0;
+ oldx = 0;
+ for(i=0; (i<100); i++) {
+ x = sigfac*sig+sig*i*0.02;
+ dp[next] = dpot(x,qq,c6,cn,npow);
+ fprintf(fp,"%10g %10g %10g\n",x,pot(x,qq,c6,cn,npow),
+ bhpot(x,qq,Abh,Bbh,Cbh));
+ if (qq != 0) {
+ if ((i > 0) && (dp[cur]*dp[next] < 0)) {
+ minimum = oldx + dp[cur]*(x-oldx)/(dp[cur]-dp[next]);
+ mval = pot(minimum,qq,c6,cn,npow);
+ printf("Van der Waals + Coulomb minimum at r = %g (nm). Value = %g (kJ/mol)\n",
+ minimum,mval);
+ }
+ }
+ cur = next;
+ oldx = x;
+
+ }
+ ffclose(fp);
+
+ do_view(oenv,ftp2fn(efXVG,NFILE,fnm),NULL);
+
+ thanx(stderr);
+
+ return 0;
+}
+
+
/*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
- *
- * The make_edi program was generously contributed by Oliver Lange, based
- * on the code from g_anaeig. You can reach him as olange@gwdg.de.
- *
+ * 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:
- * Gromacs Runs One Microsecond At Cannonball Speeds
+ * Green Red Orange Magenta Azure Cyan Skyblue
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include "gmx_header_config.h"
-
-#include <math.h>
-#include <stdlib.h>
-#include <ctype.h>
-#include <string.h>
-#include "readinp.h"
-#include "statutil.h"
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "gmx_fatal.h"
-#include "vec.h"
-#include "pbc.h"
-#include "copyrite.h"
-#include "futil.h"
-#include "statutil.h"
-#include "pdbio.h"
-#include "confio.h"
-#include "tpxio.h"
-#include "matio.h"
-#include "mshift.h"
-#include "xvgr.h"
-#include "do_fit.h"
-#include "rmpbc.h"
-#include "txtdump.h"
-#include "eigio.h"
-#include "index.h"
-
-/* Suppress Cygwin compiler warnings from using newlib version of
- * ctype.h */
-#ifdef GMX_CYGWIN
-#undef isdigit
-#endif
-
-typedef struct
-{
- real deltaF0;
- gmx_bool bHarmonic;
- gmx_bool bConstForce; /* Do constant force flooding instead of
- evaluating a flooding potential */
- real tau;
- real deltaF;
- real kT;
- real constEfl;
- real alpha2;
-} t_edflood;
-
-
-/* This type is for the average, reference, target, and origin structure */
-typedef struct edix
-{
- int nr; /* number of atoms this structure contains */
- int *anrs; /* atom index numbers */
- rvec *x; /* positions */
- real *sqrtm; /* sqrt of the masses used for mass-
- * weighting of analysis */
-} t_edix;
-
-
-typedef struct edipar
-{
- int nini; /* total Nr of atoms */
- gmx_bool fitmas; /* true if trans fit with cm */
- gmx_bool pcamas; /* true if mass-weighted PCA */
- int presteps; /* number of steps to run without any
- * perturbations ... just monitoring */
- int outfrq; /* freq (in steps) of writing to edo */
- int maxedsteps; /* max nr of steps per cycle */
- struct edix sref; /* reference positions, to these fitting
- * will be done */
- struct edix sav; /* average positions */
- struct edix star; /* target positions */
- struct edix sori; /* origin positions */
- real slope; /* minimal slope in acceptance radexp */
- int ned; /* Nr of atoms in essdyn buffer */
- t_edflood flood; /* parameters especially for flooding */
-} t_edipar;
-
-
-
-void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
- edx->nr=natoms;
- edx->anrs=index;
- edx->x=pos;
-}
-
-void write_t_edx(FILE *fp, struct edix edx, const char *comment) {
- /*here we copy only the pointers into the t_edx struct
- no data is copied and edx.box is ignored */
- int i;
- fprintf(fp,"#%s \n %d \n",comment,edx.nr);
- for (i=0;i<edx.nr;i++) {
- fprintf(fp,"%d %f %f %f\n",(edx.anrs)[i]+1,(edx.x)[i][XX],(edx.x)[i][YY],(edx.x)[i][ZZ]);
- }
-}
-
-int sscan_list(int *list[], const char *str, const char *listname) {
- /*this routine scans a string of the form 1,3-6,9 and returns the
- selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
- memory for this list will be allocated in this routine -- sscan_list expects *list to
- be a NULL-Pointer
-
- listname is a string used in the errormessage*/
-
-
- int i,istep;
- char c;
- char *pos,*startpos,*step;
- int n=strlen(str);
-
- /*enums to define the different lexical stati */
- enum { sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange };
-
- int status=sBefore; /*status of the deterministic automat to scan str */
- int number=0;
- int end_number=0;
-
- char *start=NULL; /*holds the string of the number behind a ','*/
- char *end=NULL; /*holds the string of the number behind a '-' */
-
- int nvecs=0; /* counts the number of vectors in the list*/
-
- step=NULL;
- snew(pos,n+4);
- startpos=pos;
- strcpy(pos,str);
- pos[n]=',';
- pos[n+1]='1';
- pos[n+2]='\0';
-
- *list = NULL;
-
- while ((c=*pos)!=0) {
- switch(status) {
- /* expect a number */
- case sBefore: if (isdigit(c)) {
- start=pos;
- status=sNumber;
- break;
- } else status=sError; break;
-
- /* have read a number, expect ',' or '-' */
- case sNumber: if (c==',') {
- /*store number*/
- srenew(*list,nvecs+1);
- (*list)[nvecs++]=number=strtol(start,NULL,10);
- status=sBefore;
- if (number==0)
- status=sZero;
- break;
- } else if (c=='-') { status=sMinus; break; }
- else if (isdigit(c)) break;
- else status=sError; break;
-
- /* have read a '-' -> expect a number */
- case sMinus:
- if (isdigit(c)) {
- end=pos;
- status=sRange; break;
- } else status=sError; break;
-
- case sSteppedRange:
- if (isdigit(c)) {
- if (step) {
- status=sError; break;
- } else
- step=pos;
- status=sRange;
- break;
- } else status=sError; break;
-
- /* have read the number after a minus, expect ',' or ':' */
- case sRange:
- if (c==',') {
- /*store numbers*/
- end_number=strtol(end,NULL,10);
- number=strtol(start,NULL,10);
- status=sBefore;
- if (number==0) {
- status=sZero; break;
- }
- if (end_number<=number) {
- status=sSmaller; break;
- }
- srenew(*list,nvecs+end_number-number+1);
- if (step) {
- istep=strtol(step,NULL,10);
- step=NULL;
- } else istep=1;
- for (i=number;i<=end_number;i+=istep)
- (*list)[nvecs++]=i;
- break;
- } else if (c==':') {
- status = sSteppedRange;
- break;
- } else if (isdigit(c)) break; else status=sError; break;
-
- /* format error occured */
- case sError:
- gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));
- break;
- /* logical error occured */
- case sZero:
- gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid",listname,pos-startpos);
- break;
- case sSmaller:
- gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: second index %d is not bigger than %d",listname,pos-startpos,end_number,number);
- break;
- }
- ++pos; /* read next character */
- } /*scanner has finished */
-
- /* append zero to list of eigenvectors */
- srenew(*list,nvecs+1);
- (*list)[nvecs]=0;
- sfree(startpos);
- return nvecs;
-} /*sscan_list*/
-
-void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec, const char *grouptitle,real steps[]) {
-/* eig_list is a zero-terminated list of indices into the eigvecs array.
- eigvecs are coordinates of eigenvectors
- grouptitle to write in the comment line
- steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC
-*/
-
- int n=0,i; rvec x;
- real sum;
- while (eig_list[n++]); /*count selected eigenvecs*/
-
- fprintf(fp,"# NUMBER OF EIGENVECTORS + %s\n %d\n",grouptitle,n-1);
-
- /* write list of eigenvector indicess */
- for(n=0;eig_list[n];n++) {
- if (steps)
- fprintf(fp,"%8d %g\n",eig_list[n],steps[n]);
- else
- fprintf(fp,"%8d %g\n",eig_list[n],1.0);
- }
- n=0;
-
- /* dump coordinates of the selected eigenvectors */
- while (eig_list[n]) {
- sum=0;
- for (i=0; i<natoms; i++) {
- if (eig_list[n]>nvec)
- gmx_fatal(FARGS,"Selected eigenvector %d is higher than maximum number %d of available eigenvectors",eig_list[n],nvec);
- copy_rvec(eigvecs[eig_list[n]-1][i],x);
- sum+=norm2(x);
- fprintf(fp,"%8.5f %8.5f %8.5f\n",x[XX],x[YY],x[ZZ]);
- }
- n++;
- }
-}
-
-
-/*enum referring to the different lists of eigenvectors*/
-enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON, evNr };
-#define oldMAGIC 666
-#define MAGIC 670
-
-
-void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
- int nvec, int *eig_listen[], real* evStepList[])
-{
-/* write edi-file */
-
- /*Header*/
- fprintf(fp,"#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
- MAGIC,edpars->nini,edpars->fitmas,edpars->pcamas);
- fprintf(fp,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
- edpars->outfrq,edpars->maxedsteps,edpars->slope);
- fprintf(fp,"#PRESTEPS\n %d\n#DELTA_F0\n %f\n#INIT_DELTA_F\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n#KT\n %f\n#HARMONIC\n %d\n#CONST_FORCE_FLOODING\n %d\n",
- edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,
- edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic,edpars->flood.bConstForce);
-
- /* Average and reference positions */
- write_t_edx(fp,edpars->sref,"NREF, XREF");
- write_t_edx(fp,edpars->sav,"NAV, XAV");
-
- /*Eigenvectors */
-
- write_eigvec(fp, edpars->ned, eig_listen[evMON],eigvecs,nvec,"COMPONENTS GROUP 1",NULL);
- write_eigvec(fp, edpars->ned, eig_listen[evLINFIX],eigvecs,nvec,"COMPONENTS GROUP 2",evStepList[evLINFIX]);
- write_eigvec(fp, edpars->ned, eig_listen[evLINACC],eigvecs,nvec,"COMPONENTS GROUP 3",evStepList[evLINACC]);
- write_eigvec(fp, edpars->ned, eig_listen[evRADFIX],eigvecs,nvec,"COMPONENTS GROUP 4",evStepList[evRADFIX]);
- write_eigvec(fp, edpars->ned, eig_listen[evRADACC],eigvecs,nvec,"COMPONENTS GROUP 5",NULL);
- write_eigvec(fp, edpars->ned, eig_listen[evRADCON],eigvecs,nvec,"COMPONENTS GROUP 6",NULL);
- write_eigvec(fp, edpars->ned, eig_listen[evFLOOD],eigvecs,nvec,"COMPONENTS GROUP 7",evStepList[evFLOOD]);
-
-
- /*Target and Origin positions */
- write_t_edx(fp,edpars->star,"NTARGET, XTARGET");
- write_t_edx(fp,edpars->sori,"NORIGIN, XORIGIN");
-}
-
-int read_conffile(const char *confin,char *title,rvec *x[])
-{
-/* read coordinates out of STX file */
- int natoms;
- t_atoms confat;
- matrix box;
- printf("read coordnumber from file %s\n",confin);
- get_stx_coordnum(confin,&natoms);
- printf("number of coordinates in file %d\n",natoms);
-/* if (natoms != ncoords)
- gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
- " does not match topology (= %d)",
- confin,natoms,ncoords);
- else {*/
- /* make space for coordinates and velocities */
- init_t_atoms(&confat,natoms,FALSE);
- snew(*x,natoms);
- read_stx_conf(confin,title,&confat,*x,NULL,NULL,box);
- return natoms;
-}
+#include <gmx_ana.h>
-void read_eigenvalues(int vecs[],const char *eigfile, real values[],
- gmx_bool bHesse, real kT)
-{
- int neig,nrow,i;
- double **eigval;
-
- neig = read_xvg(eigfile,&eigval,&nrow);
-
- fprintf(stderr,"Read %d eigenvalues\n",neig);
- for(i=bHesse ? 6 : 0 ; i<neig; i++) {
- if (eigval[1][i] < -0.001 && bHesse)
- fprintf(stderr,
- "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n",eigval[1][i]);
-
- if (eigval[1][i] < 0)
- eigval[1][i] = 0;
- }
- if (bHesse)
- for (i=0; vecs[i]; i++) {
- if (vecs[i]<7)
- gmx_fatal(FARGS,"ERROR: You have chosen one of the first 6 eigenvectors of the HESSE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
- values[i]=eigval[1][vecs[i]-1]/kT;
- }
- else
- for (i=0; vecs[i]; i++) {
- if (vecs[i]>(neig-6))
- gmx_fatal(FARGS,"ERROR: You have chosen one of the last 6 eigenvectors of the COVARIANCE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
- values[i]=1/eigval[1][vecs[i]-1];
- }
- /* free memory */
- for (i=0; i<nrow; i++)
- sfree(eigval[i]);
- sfree(eigval);
-}
-
-
-static real *scan_vecparams(const char *str,const char * par, int nvecs)
-{
- char f0[256],f1[256]; /*format strings adapted every pass of the loop*/
- double d,tcap=0;
- int i;
- real *vec_params;
- snew(vec_params,nvecs);
- if (str) {
- f0[0] = '\0';
- for(i=0; (i<nvecs); i++) {
- strcpy(f1,f0); /*f0 is the format string for the "to-be-ignored" numbers*/
- strcat(f1,"%lf"); /*and f1 to read the actual number in this pass of the loop*/
- if (sscanf(str,f1,&d) != 1)
- gmx_fatal(FARGS,"Not enough elements for %s parameter (I need %d)",par,nvecs);
- vec_params[i] = d;
- tcap += d;
- strcat(f0,"%*s");
- }
- }
- return vec_params;
-}
-
-
-void init_edx(struct edix *edx) {
- edx->nr=0;
- snew(edx->x,1);
- snew(edx->anrs,1);
-}
-
-void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
- atom_id igro[],rvec *x,const char* structure)
-{
-/* filter2edx copies coordinates from x to edx which are given in index
+/* This is just a wrapper binary.
*/
-
- int pos,i;
- int ix=edx->nr;
- edx->nr+=nindex;
- srenew(edx->x,edx->nr);
- srenew(edx->anrs,edx->nr);
- for (i=0;i<nindex;i++,ix++) {
- for (pos=0; pos<ngro-1 && igro[pos]!=index[i] ; ++pos) {}; /*search element in igro*/
- if (igro[pos]!=index[i])
- gmx_fatal(FARGS,"Couldn't find atom with index %d in structure %s",index[i],structure);
- edx->anrs[ix]=index[i];
- copy_rvec(x[pos],edx->x[ix]);
- }
-}
-
-void get_structure(t_atoms *atoms,const char *IndexFile,
- const char *StructureFile,struct edix *edx,int nfit,
- atom_id ifit[],int nav, atom_id index[])
+int
+main(int argc,char *argv[])
{
- atom_id *igro; /*index corresponding to target or origin structure*/
- int ngro;
- int ntar;
- rvec *xtar;
- char title[STRLEN];
- char* grpname;
-
-
- ntar=read_conffile(StructureFile,title,&xtar);
- printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
- ntar,StructureFile);
- get_index(atoms,IndexFile,1,&ngro,&igro,&grpname);
- if (ngro!=ntar)
- gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
- init_edx(edx);
- filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
-
- /* If average and reference/fitting structure differ, append the average structure as well */
- if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
- filter2edx(edx,nav,index,ngro,igro,xtar,StructureFile);
+ return gmx_make_edi(argc,argv);
}
-
-int main(int argc,char *argv[])
-{
-
- static const char *desc[] = {
- "[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
- "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
- "normal modes anaysis ([TT]g_nmeig[tt]).",
- "ED sampling can be used to manipulate the position along collective coordinates",
- "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
- "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
- "the system to explore new regions along these collective coordinates. A number",
- "of different algorithms are implemented to drive the system along the eigenvectors",
- "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
- "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
- "or to only monitor the projections of the positions onto",
- "these coordinates ([TT]-mon[tt]).[PAR]",
- "References:[BR]",
- "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
- "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
- "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
- "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
- "Towards an exhaustive sampling of the configurational spaces of the ",
- "two forms of the peptide hormone guanylin,",
- "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
- "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
- "An extended sampling of the configurational space of HPr from E. coli",
- "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
- "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
- "reference structure, target positions, etc.[PAR]",
-
- "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
- "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
- "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
- "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
- "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
- "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
- "(steps in the desired direction will be accepted, others will be rejected).",
- "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
- "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
- "to read in a structure file that defines an external origin.[PAR]",
- "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
- "towards a target structure specified with [TT]-tar[tt].[PAR]",
- "NOTE: each eigenvector can be selected only once. [PAR]",
- "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].edo[tt] file[PAR]",
- "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
- "cycle will be started if the spontaneous increase of the radius (in nm/step)",
- "is less than the value specified.[PAR]",
- "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
- "before a new cycle is started.[PAR]",
- "Note on the parallel implementation: since ED sampling is a 'global' thing",
- "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
- "is not very parallel-friendly from an implentation point of view. Because",
- "parallel ED requires some extra communication, expect the performance to be",
- "lower as in a free MD simulation, especially on a large number of nodes. [PAR]",
- "All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a .edo file. In the output",
- "file, per OUTFRQ step the following information is present: [PAR]",
- "[TT]*[tt] the step number[BR]",
- "[TT]*[tt] the number of the ED dataset. ([BB]Note[bb] that you can impose multiple ED constraints in",
- "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
- "first. The constraints are applied in the order they appear in the [TT].edi[tt] file.) [BR]",
- "[TT]*[tt] RMSD (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
- "* projections of the positions onto selected eigenvectors[BR]",
- "[PAR][PAR]",
- "FLOODING:[PAR]",
- "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
- "which will lead to extra forces expelling the structure out of the region described",
- "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
- "is kept in that region.",
- "[PAR]",
- "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
- "It can be changed with [TT]-ori[tt] to an arbitrary position in configurational space.",
- "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
- "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
- "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
- "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
- "To use constant Efl set [TT]-tau[tt] to zero.",
- "[PAR]",
- "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
- "to give good results for most standard cases in flooding of proteins.",
- "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
- "increase, this is mimicked by [GRK]alpha[grk] > 1.",
- "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
- "[PAR]",
- "RESTART and FLOODING:",
- "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
- "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
- };
-
- /* Save all the params in this struct and then save it in an edi file.
- * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
- */
- static t_edipar edi_params;
-
- enum { evStepNr = evRADFIX + 1};
- static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
- static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
- static const char* evParams[evStepNr] ={NULL,NULL};
- static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
- static const char* ConstForceStr;
- static real* evStepList[evStepNr];
- static real radfix=0.0;
- static real deltaF0=150;
- static real deltaF=0;
- static real tau=.1;
- static real constEfl=0.0;
- static real alpha=1;
- static int eqSteps=0;
- static int* listen[evNr];
- static real T=300.0;
- const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
- static gmx_bool bRestrain = FALSE;
- static gmx_bool bHesse=FALSE;
- static gmx_bool bHarmonic=FALSE;
- t_pargs pa[] = {
- { "-mon", FALSE, etSTR, {&evSelections[evMON]},
- "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
- { "-linfix", FALSE, etSTR, {&evSelections[0]},
- "Indices of eigenvectors for fixed increment linear sampling" },
- { "-linacc", FALSE, etSTR, {&evSelections[1]},
- "Indices of eigenvectors for acceptance linear sampling" },
- { "-radfix", FALSE, etSTR, {&evSelections[3]},
- "Indices of eigenvectors for fixed increment radius expansion" },
- { "-radacc", FALSE, etSTR, {&evSelections[4]},
- "Indices of eigenvectors for acceptance radius expansion" },
- { "-radcon", FALSE, etSTR, {&evSelections[5]},
- "Indices of eigenvectors for acceptance radius contraction" },
- { "-flood", FALSE, etSTR, {&evSelections[2]},
- "Indices of eigenvectors for flooding"},
- { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
- "Freqency (in steps) of writing output in [TT].edo[tt] file" },
- { "-slope", FALSE, etREAL, { &edi_params.slope},
- "Minimal slope in acceptance radius expansion"},
- { "-linstep", FALSE, etSTR, {&evParams[0]},
- "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
- { "-accdir", FALSE, etSTR, {&evParams[1]},
- "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
- { "-radstep", FALSE, etREAL, {&radfix},
- "Stepsize (nm/step) for fixed increment radius expansion"},
- { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
- "Maximum number of steps per cycle" },
- { "-eqsteps", FALSE, etINT, {&eqSteps},
- "Number of steps to run without any perturbations "},
- { "-deltaF0", FALSE,etREAL, {&deltaF0},
- "Target destabilization energy for flooding"},
- { "-deltaF", FALSE, etREAL, {&deltaF},
- "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
- { "-tau", FALSE, etREAL, {&tau},
- "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
- { "-Eflnull", FALSE, etREAL, {&constEfl},
- "The starting value of the flooding strength. The flooding strength is updated "
- "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
- { "-T", FALSE, etREAL, {&T},
- "T is temperature, the value is needed if you want to do flooding "},
- { "-alpha",FALSE,etREAL,{&alpha},
- "Scale width of gaussian flooding potential with alpha^2 "},
- { "-restrain",FALSE, etBOOL, {&bRestrain},
- "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
- { "-hessian",FALSE, etBOOL, {&bHesse},
- "The eigenvectors and eigenvalues are from a Hessian matrix"},
- { "-harmonic",FALSE, etBOOL, {&bHarmonic},
- "The eigenvalues are interpreted as spring constant"},
- { "-constF", FALSE, etSTR, {&ConstForceStr},
- "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
- "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
- };
-#define NPA asize(pa)
-
- rvec *xref1;
- int nvec1,*eignr1=NULL;
- rvec *xav1,**eigvec1=NULL;
- t_atoms *atoms=NULL;
- int nav; /* Number of atoms in the average structure */
- char *grpname;
- const char *indexfile;
- int i;
- atom_id *index,*ifit;
- int nfit; /* Number of atoms in the reference/fit structure */
- int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
- int nvecs;
- real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */
-
- const char *EdiFile;
- const char *TargetFile;
- const char *OriginFile;
- const char *EigvecFile;
-
- output_env_t oenv;
-
- /*to read topology file*/
- t_topology top;
- int ePBC;
- char title[STRLEN];
- matrix topbox;
- rvec *xtop;
- gmx_bool bTop, bFit1;
-
- t_filenm fnm[] = {
- { efTRN, "-f", "eigenvec", ffREAD },
- { efXVG, "-eig", "eigenval", ffOPTRD },
- { efTPS, NULL, NULL, ffREAD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efSTX, "-tar", "target", ffOPTRD},
- { efSTX, "-ori", "origin", ffOPTRD},
- { efEDI, "-o", "sam", ffWRITE }
- };
-#define NFILE asize(fnm)
- edi_params.outfrq=100; edi_params.slope=0.0; edi_params.maxedsteps=0;
- CopyRight(stderr,argv[0]);
- parse_common_args(&argc,argv, 0 ,
- NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
-
- indexfile=ftp2fn_null(efNDX,NFILE,fnm);
- EdiFile=ftp2fn(efEDI,NFILE,fnm);
- TargetFile = opt2fn_null("-tar",NFILE,fnm);
- OriginFile = opt2fn_null("-ori",NFILE,fnm);
-
-
- for (ev_class=0; ev_class<evNr; ++ev_class) {
- if (opt2parg_bSet(evOptions[ev_class],NPA,pa)) {
- /*get list of eigenvectors*/
- nvecs=sscan_list(&(listen[ev_class]),opt2parg_str(evOptions[ev_class],NPA,pa),evOptions[ev_class]);
- if (ev_class<evStepNr-2) {
- /*if apropriate get list of stepsizes for these eigenvectors*/
- if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
- evStepList[ev_class]=
- scan_vecparams(opt2parg_str(evStepOptions[ev_class],NPA,pa),evStepOptions[ev_class],nvecs);
- else { /*if list is not given fill with zeros */
- snew(evStepList[ev_class],nvecs);
- for (i=0; i<nvecs; i++)
- evStepList[ev_class][i]=0.0;
- }
- } else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
- snew(evStepList[ev_class],nvecs);
- for (i=0; i<nvecs; i++)
- evStepList[ev_class][i]=radfix;
- } else if (ev_class == evFLOOD) {
- snew(evStepList[ev_class],nvecs);
-
- /* Are we doing constant force flooding? In that case, we read in
- * the fproj values from the command line */
- if (opt2parg_bSet("-constF", NPA, pa))
- {
- evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
- }
- } else {}; /*to avoid ambiguity */
- } else { /* if there are no eigenvectors for this option set list to zero */
- listen[ev_class]=NULL;
- snew(listen[ev_class],1);
- listen[ev_class][0]=0;
- }
- }
-
- /* print the interpreted list of eigenvectors - to give some feedback*/
- for (ev_class=0; ev_class<evNr; ++ev_class) {
- printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
- i=0;
- while(listen[ev_class][i])
- printf("%d ",listen[ev_class][i++]);
- printf("\n");
- }
-
- EigvecFile=NULL;
- EigvecFile=opt2fn("-f",NFILE,fnm);
-
- /*read eigenvectors from eigvec.trr*/
- read_eigenvectors(EigvecFile,&nav,&bFit1,
- &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);
-
- bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
- title,&top,&ePBC,&xtop,NULL,topbox,0);
- atoms=&top.atoms;
-
-
- printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",nav);
- get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
- if (i!=nav) {
- gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
- i,nav);
- }
- printf("\n");
-
-
- if (xref1==NULL)
- {
- if (bFit1)
- {
- /* if g_covar used different coordinate groups to fit and to do the PCA */
- printf("\nNote: the structure in %s should be the same\n"
- " as the one used for the fit in g_covar\n",ftp2fn(efTPS,NFILE,fnm));
- printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
- }
- else
- {
- printf("\nNote: Apparently no fitting was done in g_covar.\n"
- " However, you need to select a reference group for fitting in mdrun\n");
- }
- get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
- snew(xref1,nfit);
- for (i=0;i<nfit;i++)
- copy_rvec(xtop[ifit[i]],xref1[i]);
- }
- else
- {
- nfit=nav;
- ifit=index;
- }
-
- if (opt2parg_bSet("-constF", NPA, pa))
- {
- /* Constant force flooding is special: Most of the normal flooding
- * options are not needed. */
- edi_params.flood.bConstForce = TRUE;
- }
- else
- {
- /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
-
- if ( listen[evFLOOD][0]!=0)
- read_eigenvalues(listen[evFLOOD],opt2fn("-eig",NFILE,fnm),evStepList[evFLOOD],bHesse,kB*T);
-
- edi_params.flood.tau=tau;
- edi_params.flood.deltaF0=deltaF0;
- edi_params.flood.deltaF=deltaF;
- edi_params.presteps=eqSteps;
- edi_params.flood.kT=kB*T;
- edi_params.flood.bHarmonic=bHarmonic;
- if (bRestrain)
- {
- /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
- edi_params.flood.constEfl=-constEfl;
- edi_params.flood.alpha2=-sqr(alpha);
- }
- else
- {
- edi_params.flood.constEfl=constEfl;
- edi_params.flood.alpha2=sqr(alpha);
- }
- }
-
- edi_params.ned=nav;
-
- /*number of system atoms */
- edi_params.nini=atoms->nr;
-
-
- /*store reference and average structure in edi_params*/
- make_t_edx(&edi_params.sref,nfit,xref1,ifit );
- make_t_edx(&edi_params.sav ,nav ,xav1 ,index);
-
-
- /* Store target positions in edi_params */
- if (opt2bSet("-tar",NFILE,fnm))
- {
- if (0 != listen[evFLOOD][0])
- {
- fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
- " You may want to use -ori to define the flooding potential center.\n\n");
- }
- get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,nav,index);
- }
- else
- {
- make_t_edx(&edi_params.star,0,NULL,index);
- }
-
- /* Store origin positions */
- if (opt2bSet("-ori",NFILE,fnm))
- {
- get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,nav,index);
- }
- else
- {
- make_t_edx(&edi_params.sori,0,NULL,index);
- }
-
- /* Write edi-file */
- write_the_whole_thing(ffopen(EdiFile,"w"), &edi_params, eigvec1,nvec1, listen, evStepList);
- thanx(stderr);
-
- return 0;
-}
-
#include <config.h>
#endif
-#include <ctype.h>
-#include "sysstuff.h"
-#include "strdb.h"
-#include "futil.h"
-#include "macros.h"
-#include "string2.h"
-#include "statutil.h"
-#include "confio.h"
-#include "copyrite.h"
-#include "typedefs.h"
-#include "index.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "index.h"
+#include <gmx_ana.h>
-#define MAXNAMES 30
-#define NAME_LEN 30
-gmx_bool bCase=FALSE;
-
-static int or_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
- atom_id *nr,atom_id *at)
-{
- atom_id i1,i2,max=0;
- gmx_bool bNotIncr;
-
- *nr=0;
-
- bNotIncr=FALSE;
- for(i1=0; i1<nr1; i1++) {
- if ((i1>0) && (at1[i1] <= max))
- bNotIncr=TRUE;
- max=at1[i1];
- }
- for(i1=0; i1<nr2; i1++) {
- if ((i1>0) && (at2[i1] <= max))
- bNotIncr=TRUE;
- max=at2[i1];
- }
-
- if (bNotIncr)
- printf("One of your groups is not ascending\n");
- else {
- i1=0;
- i2=0;
- *nr=0;
- while ((i1 < nr1) || (i2 < nr2)) {
- if ((i2 == nr2) || ((i1<nr1) && (at1[i1] < at2[i2]))) {
- at[*nr]=at1[i1];
- (*nr)++;
- i1++;
- }
- else {
- if ((i2<nr2) && ((i1==nr1) || (at1[i1] > at2[i2]))) {
- at[*nr]=at2[i2];
- (*nr)++;
- }
- i2++;
- }
- }
-
- printf("Merged two groups with OR: %u %u -> %u\n",nr1,nr2,*nr);
- }
-
- return *nr;
-}
-
-static int and_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
- atom_id *nr,atom_id *at)
-{
- atom_id i1,i2;
-
- *nr=0;
- for (i1=0; i1<nr1; i1++) {
- for (i2=0; i2<nr2; i2++) {
- if (at1[i1]==at2[i2]) {
- at[*nr]=at1[i1];
- (*nr)++;
- }
- }
- }
-
- printf("Merged two groups with AND: %u %u -> %u\n",nr1,nr2,*nr);
-
- return *nr;
-}
-
-static gmx_bool is_name_char(char c)
-{
- /* This string should contain all characters that can not be
- * the first letter of a name due to the make_ndx syntax.
- */
- const char *spec=" !&|";
-
- return (c != '\0' && strchr(spec,c) == NULL);
-}
-
-static int parse_names(char **string,int *n_names,char **names)
-{
- int i;
-
- *n_names=0;
- while (is_name_char((*string)[0]) || (*string)[0] == ' ') {
- if (is_name_char((*string)[0])) {
- if (*n_names >= MAXNAMES)
- gmx_fatal(FARGS,"To many names: %d\n",*n_names+1);
- i=0;
- while (is_name_char((*string)[i])) {
- names[*n_names][i]=(*string)[i];
- i++;
- if (i > NAME_LEN) {
- printf("Name is too long, the maximum is %d characters\n",NAME_LEN);
- return 0;
- }
- }
- names[*n_names][i]='\0';
- if (!bCase)
- upstring(names[*n_names]);
- *string += i;
- (*n_names)++;
- }
- else
- (*string)++;
- }
-
- return *n_names;
-}
-
-static gmx_bool parse_int_char(char **string,int *nr,char *c)
-{
- char *orig;
- gmx_bool bRet;
-
- orig = *string;
-
- while ((*string)[0]==' ')
- (*string)++;
-
- bRet=FALSE;
-
- *c = ' ';
-
- if (isdigit((*string)[0])) {
- *nr=(*string)[0]-'0';
- (*string)++;
- while (isdigit((*string)[0])) {
- *nr = (*nr)*10+(*string)[0]-'0';
- (*string)++;
- }
- if (isalpha((*string)[0])) {
- *c = (*string)[0];
- (*string)++;
- }
- /* Check if there is at most one non-digit character */
- if (!isalnum((*string)[0])) {
- bRet = TRUE;
- } else {
- *string = orig;
- }
- }
- else
- *nr=NOTSET;
-
- return bRet;
-}
-
-static gmx_bool parse_int(char **string,int *nr)
-{
- char *orig,c;
- gmx_bool bRet;
-
- orig = *string;
- bRet = parse_int_char(string,nr,&c);
- if (bRet && c != ' ') {
- *string = orig;
- bRet = FALSE;
- }
-
- return bRet;
-}
-
-static gmx_bool isquote(char c)
-{
- return (c == '\"');
-}
-
-static gmx_bool parse_string(char **string,int *nr, int ngrps, char **grpname)
-{
- char *s, *sp;
- char c;
-
- while ((*string)[0]==' ')
- (*string)++;
-
- (*nr) = NOTSET;
- if (isquote((*string)[0])) {
- c=(*string)[0];
- (*string)++;
- s = strdup((*string));
- sp = strchr(s, c);
- if (sp!=NULL) {
- (*string) += sp-s + 1;
- sp[0]='\0';
- (*nr) = find_group(s, ngrps, grpname);
- }
- }
-
- return (*nr) != NOTSET;
-}
-
-static int select_atomnumbers(char **string,t_atoms *atoms,atom_id n1,
- atom_id *nr,atom_id *index,char *gname)
-{
- char buf[STRLEN];
- int i,up;
-
- *nr=0;
- while ((*string)[0]==' ')
- (*string)++;
- if ((*string)[0]=='-') {
- (*string)++;
- parse_int(string,&up);
- if ((n1<1) || (n1>atoms->nr) || (up<1) || (up>atoms->nr))
- printf("Invalid atom range\n");
- else {
- for(i=n1-1; i<=up-1; i++) {
- index[*nr]=i;
- (*nr)++;
- }
- printf("Found %u atom%s in range %u-%d\n",*nr,(*nr==1)?"":"s",n1,up);
- if (n1==up)
- sprintf(buf,"a_%u",n1);
- else
- sprintf(buf,"a_%u-%d",n1,up);
- strcpy(gname,buf);
- }
- }
- else {
- i=n1;
- sprintf(gname,"a");
- do {
- if ((i-1>=0) && (i-1<atoms->nr)) {
- index[*nr] = i-1;
- (*nr)++;
- sprintf(buf,"_%d",i);
- strcat(gname,buf);
- } else {
- printf("Invalid atom number %d\n",i);
- *nr = 0;
- }
- } while ((*nr!=0) && (parse_int(string,&i)));
- }
-
- return *nr;
-}
-
-static int select_residuenumbers(char **string,t_atoms *atoms,
- atom_id n1,char c,
- atom_id *nr,atom_id *index,char *gname)
-{
- char buf[STRLEN];
- int i,j,up;
- t_resinfo *ri;
-
- *nr=0;
- while ((*string)[0]==' ')
- (*string)++;
- if ((*string)[0]=='-') {
- /* Residue number range selection */
- if (c != ' ') {
- printf("Error: residue insertion codes can not be used with residue range selection\n");
- return 0;
- }
- (*string)++;
- parse_int(string,&up);
-
- for(i=0; i<atoms->nr; i++) {
- ri = &atoms->resinfo[atoms->atom[i].resind];
- for(j=n1; (j<=up); j++) {
- if (ri->nr == j && (c == ' ' || ri->ic == c)) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- }
- printf("Found %u atom%s with res.nr. in range %u-%d\n",
- *nr,(*nr==1)?"":"s",n1,up);
- if (n1==up)
- sprintf(buf,"r_%u",n1);
- else
- sprintf(buf,"r_%u-%d",n1,up);
- strcpy(gname,buf);
- }
- else {
- /* Individual residue number/insertion code selection */
- j=n1;
- sprintf(gname,"r");
- do {
- for(i=0; i<atoms->nr; i++) {
- ri = &atoms->resinfo[atoms->atom[i].resind];
- if (ri->nr == j && ri->ic == c) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- sprintf(buf,"_%d",j);
- strcat(gname,buf);
- } while (parse_int_char(string,&j,&c));
- }
-
- return *nr;
-}
-
-static int select_residueindices(char **string,t_atoms *atoms,
- atom_id n1,char c,
- atom_id *nr,atom_id *index,char *gname)
-{
- /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
- /*resind+1 for 1-indexing*/
- char buf[STRLEN];
- int i,j,up;
- t_resinfo *ri;
-
- *nr=0;
- while ((*string)[0]==' ')
- (*string)++;
- if ((*string)[0]=='-') {
- /* Residue number range selection */
- if (c != ' ') {
- printf("Error: residue insertion codes can not be used with residue range selection\n");
- return 0;
- }
- (*string)++;
- parse_int(string,&up);
-
- for(i=0; i<atoms->nr; i++) {
- ri = &atoms->resinfo[atoms->atom[i].resind];
- for(j=n1; (j<=up); j++) {
- if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c)) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- }
- printf("Found %u atom%s with resind.+1 in range %u-%d\n",
- *nr,(*nr==1)?"":"s",n1,up);
- if (n1==up)
- sprintf(buf,"r_%u",n1);
- else
- sprintf(buf,"r_%u-%d",n1,up);
- strcpy(gname,buf);
- }
- else {
- /* Individual residue number/insertion code selection */
- j=n1;
- sprintf(gname,"r");
- do {
- for(i=0; i<atoms->nr; i++) {
- ri = &atoms->resinfo[atoms->atom[i].resind];
- if (atoms->atom[i].resind+1 == j && ri->ic == c) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- sprintf(buf,"_%d",j);
- strcat(gname,buf);
- } while (parse_int_char(string,&j,&c));
- }
-
- return *nr;
-}
-
-
-static gmx_bool atoms_from_residuenumbers(t_atoms *atoms,int group,t_blocka *block,
- atom_id *nr,atom_id *index,char *gname)
-{
- int i,j,j0,j1,resnr,nres;
-
- j0=block->index[group];
- j1=block->index[group+1];
- nres = atoms->nres;
- for(j=j0; j<j1; j++)
- if (block->a[j]>=nres) {
- printf("Index %s contains number>nres (%d>%d)\n",
- gname,block->a[j]+1,nres);
- return FALSE;
- }
- for(i=0; i<atoms->nr; i++) {
- resnr = atoms->resinfo[atoms->atom[i].resind].nr;
- for (j=j0; j<j1; j++)
- if (block->a[j]+1 == resnr) {
- index[*nr]=i;
- (*nr)++;
- break;
- }
- }
- printf("Found %u atom%s in %d residues from group %s\n",
- *nr,(*nr==1)?"":"s",j1-j0,gname);
- return *nr;
-}
-
-static gmx_bool comp_name(char *name,char *search)
-{
- while (name[0] != '\0' && search[0] != '\0') {
- switch (search[0]) {
- case '?':
- /* Always matches */
- break;
- case '*':
- if (search[1] != '\0') {
- printf("WARNING: Currently '*' is only supported at the end of an expression\n");
- return FALSE;
- } else {
- return TRUE;
- }
- break;
- default:
- /* Compare a single character */
- if (( bCase && strncmp(name,search,1)) ||
- (!bCase && gmx_strncasecmp(name,search,1))) {
- return FALSE;
- }
- }
- name++;
- search++;
- }
-
- return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
-}
-
-static int select_chainnames(t_atoms *atoms,int n_names,char **names,
- atom_id *nr,atom_id *index)
-{
- char name[2];
- int j;
- atom_id i;
-
- name[1]=0;
- *nr=0;
- for(i=0; i<atoms->nr; i++) {
- name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
- j=0;
- while (j<n_names && !comp_name(name,names[j]))
- j++;
- if (j<n_names) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- printf("Found %u atom%s with chain identifier%s",
- *nr,(*nr==1)?"":"s",(n_names==1)?"":"s");
- for(j=0; (j<n_names); j++)
- printf(" %s",names[j]);
- printf("\n");
-
- return *nr;
-}
-
-static int select_atomnames(t_atoms *atoms,int n_names,char **names,
- atom_id *nr,atom_id *index,gmx_bool bType)
-{
- char *name;
- int j;
- atom_id i;
-
- *nr=0;
- for(i=0; i<atoms->nr; i++) {
- if (bType)
- name=*(atoms->atomtype[i]);
- else
- name=*(atoms->atomname[i]);
- j=0;
- while (j<n_names && !comp_name(name,names[j]))
- j++;
- if (j<n_names) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- printf("Found %u atoms with %s%s",
- *nr,bType ? "type" : "name",(n_names==1)?"":"s");
- for(j=0; (j<n_names); j++)
- printf(" %s",names[j]);
- printf("\n");
-
- return *nr;
-}
-
-static int select_residuenames(t_atoms *atoms,int n_names,char **names,
- atom_id *nr,atom_id *index)
-{
- char *name;
- int j;
- atom_id i;
-
- *nr=0;
- for(i=0; i<atoms->nr; i++) {
- name = *(atoms->resinfo[atoms->atom[i].resind].name);
- j=0;
- while (j<n_names && !comp_name(name,names[j]))
- j++;
- if (j<n_names) {
- index[*nr]=i;
- (*nr)++;
- }
- }
- printf("Found %u atoms with residue name%s",*nr,(n_names==1)?"":"s");
- for(j=0; (j<n_names); j++)
- printf(" %s",names[j]);
- printf("\n");
-
- return *nr;
-}
-
-static void copy2block(int n,atom_id *index,t_blocka *block)
-{
- int i,n0;
-
- block->nr++;
- n0=block->nra;
- block->nra=n0+n;
- srenew(block->index,block->nr+1);
- block->index[block->nr]=n0+n;
- srenew(block->a,n0+n);
- for(i=0; (i<n); i++)
- block->a[n0+i]=index[i];
-}
-
-static void make_gname(int n,char **names,char *gname)
-{
- int i;
-
- strcpy(gname,names[0]);
- for (i=1; i<n; i++) {
- strcat(gname,"_");
- strcat(gname,names[i]);
- }
-}
-
-static void copy_group(int g,t_blocka *block,atom_id *nr,atom_id *index)
-{
- int i,i0;
-
- i0=block->index[g];
- *nr=block->index[g+1]-i0;
- for (i=0; i<*nr; i++)
- index[i]=block->a[i0+i];
-}
-
-static void remove_group(int nr,int nr2,t_blocka *block,char ***gn)
-{
- int i,j,shift;
- char *name;
-
- if (nr2==NOTSET)
- nr2=nr;
-
- for(j=0; j<=nr2-nr; j++) {
- if ((nr<0) || (nr>=block->nr))
- printf("Group %d does not exist\n",nr+j);
- else {
- shift=block->index[nr+1]-block->index[nr];
- for(i=block->index[nr+1]; i<block->nra; i++)
- block->a[i-shift]=block->a[i];
-
- for(i=nr; i<block->nr; i++)
- block->index[i]=block->index[i+1]-shift;
- name = strdup((*gn)[nr]);
- sfree((*gn)[nr]);
- for(i=nr; i<block->nr-1; i++) {
- (*gn)[i]=(*gn)[i+1];
- }
- block->nr--;
- block->nra=block->index[block->nr];
- printf("Removed group %d '%s'\n",nr+j,name);
- sfree(name);
- }
- }
-}
-
-static void split_group(t_atoms *atoms,int sel_nr,t_blocka *block,char ***gn,
- gmx_bool bAtom)
-{
- char buf[STRLEN],*name;
- int i,resind;
- atom_id a,n0,n1;
-
- printf("Splitting group %d '%s' into %s\n",sel_nr,(*gn)[sel_nr],
- bAtom ? "atoms" : "residues");
-
- n0=block->index[sel_nr];
- n1=block->index[sel_nr+1];
- srenew(block->a,block->nra+n1-n0);
- for (i=n0; i<n1; i++) {
- a=block->a[i];
- resind = atoms->atom[a].resind;
- name = *(atoms->resinfo[resind].name);
- if (bAtom || (i==n0) || (atoms->atom[block->a[i-1]].resind!=resind)) {
- if (i>n0)
- block->index[block->nr]=block->nra;
- block->nr++;
- srenew(block->index,block->nr+1);
- srenew(*gn,block->nr);
- if (bAtom)
- sprintf(buf,"%s_%s_%u",(*gn)[sel_nr],*atoms->atomname[a],a+1);
- else
- sprintf(buf,"%s_%s_%d",(*gn)[sel_nr],name,atoms->resinfo[resind].nr);
- (*gn)[block->nr-1]=strdup(buf);
- }
- block->a[block->nra]=a;
- block->nra++;
- }
- block->index[block->nr]=block->nra;
-}
-
-static int split_chain(t_atoms *atoms,rvec *x,
- int sel_nr,t_blocka *block,char ***gn)
-{
- char buf[STRLEN];
- int j,nchain;
- atom_id i,a,natoms,*start=NULL,*end=NULL,ca_start,ca_end;
- rvec vec;
-
- natoms=atoms->nr;
- nchain=0;
- ca_start=0;
-
- while (ca_start<natoms) {
- while((ca_start<natoms) && strcmp(*atoms->atomname[ca_start],"CA"))
- ca_start++;
- if (ca_start<natoms) {
- srenew(start,nchain+1);
- srenew(end,nchain+1);
- start[nchain]=ca_start;
- while ((start[nchain]>0) &&
- (atoms->atom[start[nchain]-1].resind ==
- atoms->atom[ca_start].resind))
- start[nchain]--;
-
- i=ca_start;
- do {
- ca_end=i;
- do {
- i++;
- } while ((i<natoms) && strcmp(*atoms->atomname[i],"CA"));
- if (i<natoms)
- rvec_sub(x[ca_end],x[i],vec);
- } while ((i<natoms) && (norm(vec)<0.45));
-
- end[nchain]=ca_end;
- while ((end[nchain]+1<natoms) &&
- (atoms->atom[end[nchain]+1].resind==atoms->atom[ca_end].resind))
- end[nchain]++;
- ca_start=end[nchain]+1;
- nchain++;
- }
- }
- if (nchain==1)
- printf("Found 1 chain, will not split\n");
- else
- printf("Found %d chains\n",nchain);
- for (j=0; j<nchain; j++)
- printf("%d:%6u atoms (%u to %u)\n",
- j+1,end[j]-start[j]+1,start[j]+1,end[j]+1);
-
- if (nchain > 1) {
- srenew(block->a,block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
- for (j=0; j<nchain; j++) {
- block->nr++;
- srenew(block->index,block->nr+1);
- srenew(*gn,block->nr);
- sprintf(buf,"%s_chain%d",(*gn)[sel_nr],j+1);
- (*gn)[block->nr-1]=strdup(buf);
- for (i=block->index[sel_nr]; i<block->index[sel_nr+1]; i++) {
- a=block->a[i];
- if ((a>=start[j]) && (a<=end[j])) {
- block->a[block->nra]=a;
- block->nra++;
- }
- }
- block->index[block->nr]=block->nra;
- if (block->index[block->nr-1]==block->index[block->nr])
- remove_group(block->nr-1,NOTSET,block,gn);
- }
- }
- sfree(start);
- sfree(end);
-
- return nchain;
-}
-
-static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
-{
- if ( atoms==NULL ) {
- printf("Can not process '%s' without atom info, use option -f\n", string);
- return FALSE;
- } else
- return TRUE;
-}
-
-static gmx_bool parse_entry(char **string,int natoms,t_atoms *atoms,
- t_blocka *block,char ***gn,
- atom_id *nr,atom_id *index,char *gname)
-{
- static char **names, *ostring;
- static gmx_bool bFirst=TRUE;
- int j,n_names,sel_nr1;
- atom_id i,nr1,*index1;
- char c;
- gmx_bool bRet,bCompl;
-
- if (bFirst) {
- bFirst=FALSE;
- snew(names,MAXNAMES);
- for (i=0; i<MAXNAMES; i++)
- snew(names[i],NAME_LEN+1);
- }
-
- bRet=FALSE;
- sel_nr1=NOTSET;
-
- while(*string[0]==' ')
- (*string)++;
-
- if ((*string)[0]=='!') {
- bCompl=TRUE;
- (*string)++;
- while(*string[0]==' ')
- (*string)++;
- } else
- bCompl=FALSE;
-
- ostring = *string;
-
- if (parse_int(string,&sel_nr1) ||
- parse_string(string,&sel_nr1,block->nr,*gn)) {
- if ((sel_nr1>=0) && (sel_nr1<block->nr)) {
- copy_group(sel_nr1,block,nr,index);
- strcpy(gname,(*gn)[sel_nr1]);
- printf("Copied index group %d '%s'\n",sel_nr1,(*gn)[sel_nr1]);
- bRet=TRUE;
- } else
- printf("Group %d does not exist\n",sel_nr1);
- }
- else if ((*string)[0]=='a') {
- (*string)++;
- if (check_have_atoms(atoms, ostring)) {
- if (parse_int(string,&sel_nr1)) {
- bRet=select_atomnumbers(string,atoms,sel_nr1,nr,index,gname);
- }
- else if (parse_names(string,&n_names,names)) {
- bRet=select_atomnames(atoms,n_names,names,nr,index,FALSE);
- make_gname(n_names,names,gname);
- }
- }
- }
- else if ((*string)[0]=='t') {
- (*string)++;
- if (check_have_atoms(atoms, ostring) &&
- parse_names(string,&n_names,names)) {
- if (atoms->atomtype == NULL)
- printf("Need a run input file to select atom types\n");
- else {
- bRet=select_atomnames(atoms,n_names,names,nr,index,TRUE);
- make_gname(n_names,names,gname);
- }
- }
- }
- else if (strncmp(*string,"res",3)==0) {
- (*string)+=3;
- if ( check_have_atoms(atoms, ostring) &&
- parse_int(string,&sel_nr1) &&
- (sel_nr1>=0) && (sel_nr1<block->nr) ) {
- bRet=atoms_from_residuenumbers(atoms,
- sel_nr1,block,nr,index,(*gn)[sel_nr1]);
- sprintf(gname,"atom_%s",(*gn)[sel_nr1]);
- }
- }
- else if (strncmp(*string,"ri",2)==0) {
- (*string)+=2;
- if (check_have_atoms(atoms, ostring) &&
- parse_int_char(string,&sel_nr1,&c)) {
- bRet=select_residueindices(string,atoms,sel_nr1,c,nr,index,gname);
- }
- }
- else if ((*string)[0]=='r') {
- (*string)++;
- if (check_have_atoms(atoms, ostring)) {
- if (parse_int_char(string,&sel_nr1,&c)) {
- bRet=select_residuenumbers(string,atoms,sel_nr1,c,nr,index,gname);
- }
- else if (parse_names(string,&n_names,names)) {
- bRet=select_residuenames(atoms,n_names,names,nr,index);
- make_gname(n_names,names,gname);
- }
- }
- }
- else if (strncmp(*string,"chain",5)==0) {
- (*string)+=5;
- if (check_have_atoms(atoms, ostring) &&
- parse_names(string,&n_names,names)) {
- bRet=select_chainnames(atoms,n_names,names,nr,index);
- sprintf(gname,"ch%s",names[0]);
- for (i=1; i<n_names; i++)
- strcat(gname,names[i]);
- }
- }
- if (bRet && bCompl) {
- snew(index1,natoms-*nr);
- nr1=0;
- for(i=0; i<natoms; i++) {
- j=0;
- while ((j<*nr) && (index[j] != i))
- j++;
- if (j==*nr) {
- if (nr1 >= natoms-*nr) {
- printf("There are double atoms in your index group\n");
- break;
- }
- index1[nr1]=i;
- nr1++;
- }
- }
- *nr=nr1;
- for(i=0; i<nr1; i++)
- index[i]=index1[i];
- sfree(index1);
-
- for (i=strlen(gname)+1; i>0; i--)
- gname[i]=gname[i-1];
- gname[0]='!';
- printf("Complemented group: %u atoms\n",*nr);
- }
-
- return bRet;
-}
-
-static void list_residues(t_atoms *atoms)
-{
- int i,j,start,end,prev_resind,resind;
- gmx_bool bDiff;
-
- /* Print all the residues, assuming continuous resnr count */
- start = atoms->atom[0].resind;
- prev_resind = start;
- for(i=0; i<atoms->nr; i++) {
- resind = atoms->atom[i].resind;
- if ((resind != prev_resind) || (i==atoms->nr-1)) {
- if ((bDiff=strcmp(*atoms->resinfo[resind].name,
- *atoms->resinfo[start].name)) ||
- (i==atoms->nr-1)) {
- if (bDiff)
- end = prev_resind;
- else
- end = resind;
- if (end < start+3)
- for(j=start; j<=end; j++)
- printf("%4d %-5s",
- j+1,*(atoms->resinfo[j].name));
- else
- printf(" %4d - %4d %-5s ",
- start+1,end+1,*(atoms->resinfo[start].name));
- start = resind;
- }
- }
- prev_resind = resind;
- }
- printf("\n");
-}
-
-static void edit_index(int natoms, t_atoms *atoms,rvec *x,t_blocka *block, char ***gn, gmx_bool bVerbose)
-{
- static char **atnames, *ostring;
- static gmx_bool bFirst=TRUE;
- char inp_string[STRLEN],*string;
- char gname[STRLEN],gname1[STRLEN],gname2[STRLEN];
- int i,i0,i1,sel_nr,sel_nr2,newgroup;
- atom_id nr,nr1,nr2,*index,*index1,*index2;
- gmx_bool bAnd,bOr,bPrintOnce;
-
- if (bFirst) {
- bFirst=FALSE;
- snew(atnames,MAXNAMES);
- for (i=0; i<MAXNAMES; i++)
- snew(atnames[i],NAME_LEN+1);
- }
-
- string=NULL;
-
- snew(index,natoms);
- snew(index1,natoms);
- snew(index2,natoms);
-
- newgroup=NOTSET;
- bPrintOnce = TRUE;
- do {
- gname1[0]='\0';
- if (bVerbose || bPrintOnce || newgroup!=NOTSET) {
- printf("\n");
- if (bVerbose || bPrintOnce || newgroup==NOTSET) {
- i0=0;
- i1=block->nr;
- } else {
- i0=newgroup;
- i1=newgroup+1;
- }
- for(i=i0; i<i1; i++)
- printf("%3d %-20s: %5u atoms\n",i,(*gn)[i],
- block->index[i+1]-block->index[i]);
- newgroup=NOTSET;
- }
- if (bVerbose || bPrintOnce) {
- printf("\n");
- printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
- printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
- printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
- printf(" 'r': residue 'res' nr 'chain' char\n");
- printf(" \"name\": group 'case': case %s 'q': save and quit\n",
- bCase ? "insensitive" : "sensitive ");
- printf(" 'ri': residue index\n");
- bPrintOnce = FALSE;
- }
- printf("\n");
- printf("> ");
- if(NULL==fgets(inp_string,STRLEN,stdin))
- {
- gmx_fatal(FARGS,"Error reading user input");
- }
- inp_string[strlen(inp_string)-1]=0;
- printf("\n");
- string=inp_string;
- while (string[0]==' ')
- string++;
-
- ostring = string;
- nr=0;
- if (string[0] == 'h') {
- printf(" nr : selects an index group by number or quoted string.\n");
- printf(" The string is first matched against the whole group name,\n");
- printf(" then against the beginning and finally against an\n");
- printf(" arbitrary substring. A multiple match is an error.\n");
-
- printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
- printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
- printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
- printf(" wildcard '*' allowed at the end of a name.\n");
- printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
- printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
- printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
- printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
- printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
- printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
- printf(" not available with a .gro file as input.\n");
- printf(" ! : takes the complement of a group with respect to all\n");
- printf(" the atoms in the input file.\n");
- printf(" & | : AND and OR, can be placed between any of the options\n");
- printf(" above, the input is processed from left to right.\n");
- printf(" 'name' nr name : rename group nr to name.\n");
- printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
- printf(" 'keep' nr : deletes all groups except nr.\n");
- printf(" 'case' : make all name compares case (in)sensitive.\n");
- printf(" 'splitch' nr : split group into chains using CA distances.\n");
- printf(" 'splitres' nr : split group into residues.\n");
- printf(" 'splitat' nr : split group into atoms.\n");
- printf(" 'res' nr : interpret numbers in group as residue numbers\n");
- printf(" Enter : list the currently defined groups and commands\n");
- printf(" 'l' : list the residues.\n");
- printf(" 'h' : show this help.\n");
- printf(" 'q' : save and quit.\n");
- printf("\n");
- printf(" Examples:\n");
- printf(" > 2 | 4 & r 3-5\n");
- printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
- printf(" > a C* & !a C CA\n");
- printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
- printf(" > \"protein\" & ! \"backb\"\n");
- printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
- if (bVerbose) {
- printf("\npress Enter ");
- getchar();
- }
- }
- else if (strncmp(string,"del",3)==0) {
- string+=3;
- if (parse_int(&string,&sel_nr)) {
- while(string[0]==' ')
- string++;
- if (string[0]=='-') {
- string++;
- parse_int(&string,&sel_nr2);
- } else
- sel_nr2=NOTSET;
- while(string[0]==' ')
- string++;
- if (string[0]=='\0')
- remove_group(sel_nr,sel_nr2,block,gn);
- else
- printf("\nSyntax error: \"%s\"\n",string);
- }
- }
- else if (strncmp(string,"keep",4)==0) {
- string+=4;
- if (parse_int(&string,&sel_nr)) {
- remove_group(sel_nr+1,block->nr-1,block,gn);
- remove_group(0,sel_nr-1,block,gn);
- }
- }
- else if (strncmp(string,"name",4)==0) {
- string+=4;
- if (parse_int(&string,&sel_nr)) {
- if ((sel_nr>=0) && (sel_nr<block->nr)) {
- sscanf(string,"%s",gname);
- sfree((*gn)[sel_nr]);
- (*gn)[sel_nr]=strdup(gname);
- }
- }
- }
- else if (strncmp(string,"case",4)==0) {
- bCase=!bCase;
- printf("Switched to case %s\n",bCase ? "sensitive" : "insensitive");
- }
- else if (string[0] == 'v' ) {
- bVerbose=!bVerbose;
- printf("Turned verbose %s\n",bVerbose ? "on" : "off");
- }
- else if (string[0] == 'l') {
- if ( check_have_atoms(atoms, ostring) )
- list_residues(atoms);
- }
- else if (strncmp(string,"splitch",7)==0) {
- string+=7;
- if ( check_have_atoms(atoms, ostring) &&
- parse_int(&string,&sel_nr) &&
- (sel_nr>=0) && (sel_nr<block->nr))
- split_chain(atoms,x,sel_nr,block,gn);
- }
- else if (strncmp(string,"splitres",8)==0 ) {
- string+=8;
- if ( check_have_atoms(atoms, ostring) &&
- parse_int(&string,&sel_nr) &&
- (sel_nr>=0) && (sel_nr<block->nr))
- split_group(atoms,sel_nr,block,gn,FALSE);
- }
- else if (strncmp(string,"splitat",7)==0 ) {
- string+=7;
- if ( check_have_atoms(atoms, ostring) &&
- parse_int(&string,&sel_nr) &&
- (sel_nr>=0) && (sel_nr<block->nr))
- split_group(atoms,sel_nr,block,gn,TRUE);
- }
- else if (string[0] == '\0') {
- bPrintOnce = TRUE;
- }
- else if (string[0] != 'q') {
- nr1=-1;
- nr2=-1;
- if (parse_entry(&string,natoms,atoms,block,gn,&nr,index,gname)) {
- do {
- while (string[0]==' ')
- string++;
-
- bAnd=FALSE;
- bOr=FALSE;
- if (string[0]=='&')
- bAnd=TRUE;
- else if (string[0]=='|')
- bOr=TRUE;
-
- if (bAnd || bOr) {
- string++;
- nr1=nr;
- for(i=0; i<nr; i++)
- index1[i]=index[i];
- strcpy(gname1,gname);
- if (parse_entry(&string,natoms,atoms,block,gn,&nr2,index2,gname2)) {
- if (bOr) {
- or_groups(nr1,index1,nr2,index2,&nr,index);
- sprintf(gname,"%s_%s",gname1,gname2);
- }
- else {
- and_groups(nr1,index1,nr2,index2,&nr,index);
- sprintf(gname,"%s_&_%s",gname1,gname2);
- }
- }
- }
- } while (bAnd || bOr);
- }
- while(string[0]==' ')
- string++;
- if (string[0])
- printf("\nSyntax error: \"%s\"\n",string);
- else if (nr>0) {
- copy2block(nr,index,block);
- srenew(*gn,block->nr);
- newgroup = block->nr-1;
- (*gn)[newgroup]=strdup(gname);
- }
- else
- printf("Group is empty\n");
- }
- } while (string[0]!='q');
-
- sfree(index);
- sfree(index1);
- sfree(index2);
-}
-
-static int block2natoms(t_blocka *block)
+/* This is just a wrapper binary.
+*/
+int
+main(int argc,char *argv[])
{
- int i, natoms;
-
- natoms = 0;
- for(i=0; i<block->nra; i++)
- natoms = max(natoms, block->a[i]);
- natoms++;
-
- return natoms;
-}
-
-void merge_blocks(t_blocka *dest, t_blocka *source)
-{
- int i,nra0,i0;
-
- /* count groups, srenew and fill */
- i0 = dest->nr;
- nra0 = dest->nra;
- dest->nr+=source->nr;
- srenew(dest->index, dest->nr+1);
- for(i=0; i<source->nr; i++)
- dest->index[i0+i] = nra0 + source->index[i];
- /* count atoms, srenew and fill */
- dest->nra+=source->nra;
- srenew(dest->a, dest->nra);
- for(i=0; i<source->nra; i++)
- dest->a[nra0+i] = source->a[i];
-
- /* terminate list */
- dest->index[dest->nr]=dest->nra;
-
-}
-
-int main(int argc,char *argv[])
-{
- const char *desc[] = {
- "Index groups are necessary for almost every gromacs program.",
- "All these programs can generate default index groups. You ONLY",
- "have to use [TT]make_ndx[tt] when you need SPECIAL index groups.",
- "There is a default index group for the whole system, 9 default",
- "index groups for proteins, and a default index group",
- "is generated for every other residue name.[PAR]",
- "When no index file is supplied, also [TT]make_ndx[tt] will generate the",
- "default groups.",
- "With the index editor you can select on atom, residue and chain names",
- "and numbers.",
- "When a run input file is supplied you can also select on atom type.",
- "You can use NOT, AND and OR, you can split groups",
- "into chains, residues or atoms. You can delete and rename groups.[PAR]",
- "The atom numbering in the editor and the index file starts at 1."
- };
-
- static int natoms=0;
- static gmx_bool bVerbose=FALSE;
- t_pargs pa[] = {
- { "-natoms", FALSE, etINT, {&natoms},
- "set number of atoms (default: read from coordinate or index file)" },
- { "-verbose", FALSE, etBOOL, {&bVerbose},
- "HIDDENVerbose output" }
- };
-#define NPA asize(pa)
-
- output_env_t oenv;
- char title[STRLEN];
- int nndxin;
- const char *stxfile;
- char **ndxinfiles;
- const char *ndxoutfile;
- gmx_bool bNatoms;
- int i,j;
- t_atoms *atoms;
- rvec *x,*v;
- int ePBC;
- matrix box;
- t_blocka *block,*block2;
- char **gnames,**gnames2;
- t_filenm fnm[] = {
- { efSTX, "-f", NULL, ffOPTRD },
- { efNDX, "-n", NULL, ffOPTRDMULT },
- { efNDX, "-o", NULL, ffWRITE }
- };
-#define NFILE asize(fnm)
-
- CopyRight(stderr,argv[0]);
-
- parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,asize(desc),desc,
- 0,NULL,&oenv);
-
- stxfile = ftp2fn_null(efSTX,NFILE,fnm);
- if (opt2bSet("-n",NFILE,fnm)) {
- nndxin = opt2fns(&ndxinfiles,"-n",NFILE,fnm);
- } else {
- nndxin = 0;
- }
- ndxoutfile = opt2fn("-o",NFILE,fnm);
- bNatoms = opt2parg_bSet("-natoms",NPA,pa);
-
- if (!stxfile && !nndxin)
- gmx_fatal(FARGS,"No input files (structure or index)");
-
- if (stxfile) {
- snew(atoms,1);
- get_stx_coordnum(stxfile,&(atoms->nr));
- init_t_atoms(atoms,atoms->nr,TRUE);
- snew(x,atoms->nr);
- snew(v,atoms->nr);
- fprintf(stderr,"\nReading structure file\n");
- read_stx_conf(stxfile,title,atoms,x,v,&ePBC,box);
- natoms = atoms->nr;
- bNatoms=TRUE;
- } else {
- atoms = NULL;
- x = NULL;
- }
-
- /* read input file(s) */
- block = new_blocka();
- gnames = NULL;
- printf("Going to read %d old index file(s)\n",nndxin);
- if (nndxin) {
- for(i=0; i<nndxin; i++) {
- block2 = init_index(ndxinfiles[i],&gnames2);
- srenew(gnames, block->nr+block2->nr);
- for(j=0; j<block2->nr; j++)
- gnames[block->nr+j]=gnames2[j];
- sfree(gnames2);
- merge_blocks(block, block2);
- sfree(block2->a);
- sfree(block2->index);
-/* done_block(block2); */
- sfree(block2);
- }
- }
- else {
- snew(gnames,1);
- analyse(atoms,block,&gnames,FALSE,TRUE);
- }
-
- if (!bNatoms) {
- natoms = block2natoms(block);
- printf("Counted atom numbers up to %d in index file\n", natoms);
- }
-
- edit_index(natoms,atoms,x,block,&gnames,bVerbose);
-
- write_index(ndxoutfile,block,gnames);
-
- thanx(stderr);
-
- return 0;
+ return gmx_make_ndx(argc,argv);
}
#include <config.h>
#endif
-#include <math.h>
-#include "typedefs.h"
-#include "smalloc.h"
-#include "copyrite.h"
-#include "statutil.h"
-#include "macros.h"
-#include "string2.h"
-#include "futil.h"
-#include "gmx_fatal.h"
+#include <gmx_ana.h>
-static int calc_ntype(int nft,int *ft,t_idef *idef)
-{
- int i,f,nf=0;
-
- for(i=0; (i<idef->ntypes); i++) {
- for(f=0; f<nft; f++) {
- if (idef->functype[i] == ft[f])
- nf++;
- }
- }
-
- return nf;
-}
-
-static void fill_ft_ind(int nft,int *ft,t_idef *idef,
- int ft_ind[],char *grpnames[])
-{
- char buf[125];
- int i,f,ftype,ind=0;
-
- /* Loop over all the function types in the topology */
- for(i=0; (i<idef->ntypes); i++) {
- ft_ind[i] = -1;
- /* Check all the selected function types */
- for(f=0; f<nft; f++) {
- ftype = ft[f];
- if (idef->functype[i] == ftype) {
- ft_ind[i] = ind;
- switch (ftype) {
- case F_ANGLES:
- sprintf(buf,"Theta=%.1f_%.2f",idef->iparams[i].harmonic.rA,
- idef->iparams[i].harmonic.krA);
- break;
- case F_G96ANGLES:
- sprintf(buf,"Cos_th=%.1f_%.2f",idef->iparams[i].harmonic.rA,
- idef->iparams[i].harmonic.krA);
- break;
- case F_UREY_BRADLEY:
- sprintf(buf,"UB_th=%.1f_%.2f2f",idef->iparams[i].u_b.thetaA,
- idef->iparams[i].u_b.kthetaA);
- break;
- case F_QUARTIC_ANGLES:
- sprintf(buf,"Q_th=%.1f_%.2f_%.2f",idef->iparams[i].qangle.theta,
- idef->iparams[i].qangle.c[0],idef->iparams[i].qangle.c[1]);
- break;
- case F_TABANGLES:
- sprintf(buf,"Table=%d_%.2f",idef->iparams[i].tab.table,
- idef->iparams[i].tab.kA);
- break;
- case F_PDIHS:
- sprintf(buf,"Phi=%.1f_%d_%.2f",idef->iparams[i].pdihs.phiA,
- idef->iparams[i].pdihs.mult,idef->iparams[i].pdihs.cpA);
- break;
- case F_IDIHS:
- sprintf(buf,"Xi=%.1f_%.2f",idef->iparams[i].harmonic.rA,
- idef->iparams[i].harmonic.krA);
- break;
- case F_RBDIHS:
- sprintf(buf,"RB-A1=%.2f",idef->iparams[i].rbdihs.rbcA[1]);
- break;
- default:
- gmx_fatal(FARGS,"Unsupported function type '%s' selected",
- interaction_function[ftype].longname);
- }
- grpnames[ind]=strdup(buf);
- ind++;
- }
- }
- }
-}
-
-static void fill_ang(int nft,int *ft,int fac,
- int nr[],int *index[],int ft_ind[],t_topology *top,
- gmx_bool bNoH,real hq)
-{
- int f,ftype,i,j,indg,nr_fac;
- gmx_bool bUse;
- t_idef *idef;
- t_atom *atom;
- t_iatom *ia;
-
-
- idef = &top->idef;
- atom = top->atoms.atom;
- for(f=0; f<nft; f++) {
- ftype = ft[f];
- ia = idef->il[ftype].iatoms;
- for(i=0; (i<idef->il[ftype].nr); ) {
- indg = ft_ind[ia[0]];
- if (indg == -1)
- gmx_incons("Routine fill_ang");
- bUse = TRUE;
- if (bNoH) {
- for(j=0; j<fac; j++) {
- if (atom[ia[1+j]].m < 1.5)
- bUse = FALSE;
- }
- }
- if (hq) {
- for(j=0; j<fac; j++) {
- if (atom[ia[1+j]].m < 1.5 && fabs(atom[ia[1+j]].q) < hq)
- bUse = FALSE;
- }
- }
- if (bUse) {
- if (nr[indg] % 1000 == 0) {
- srenew(index[indg],fac*(nr[indg]+1000));
- }
- nr_fac = fac*nr[indg];
- for(j=0; (j<fac); j++)
- index[indg][nr_fac+j] = ia[j+1];
- nr[indg]++;
- }
- ia += interaction_function[ftype].nratoms+1;
- i += interaction_function[ftype].nratoms+1;
- }
- }
-}
-
-static int *select_ftype(const char *opt,int *nft,int *mult)
-{
- int *ft=NULL,ftype;
-
- if (opt[0] == 'a') {
- *mult = 3;
- for(ftype=0; ftype<F_NRE; ftype++) {
- if ((interaction_function[ftype].flags & IF_ATYPE) ||
- ftype == F_TABANGLES) {
- (*nft)++;
- srenew(ft,*nft);
- ft[*nft-1] = ftype;
- }
- }
- } else {
- *mult = 4;
- *nft = 1;
- snew(ft,*nft);
- switch(opt[0]) {
- case 'd':
- ft[0] = F_PDIHS;
- break;
- case 'i':
- ft[0] = F_IDIHS;
- break;
- case 'r':
- ft[0] = F_RBDIHS;
- break;
- default:
- break;
- }
- }
-
- return ft;
-}
-
-int main(int argc,char *argv[])
+/* This is just a wrapper binary.
+*/
+int
+main(int argc,char *argv[])
{
- static const char *desc[] = {
- "[TT]mk_angndx[tt] makes an index file for calculation of",
- "angle distributions etc. It uses a run input file ([TT].tpx[tt]) for the",
- "definitions of the angles, dihedrals etc."
- };
- static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
- static gmx_bool bH=TRUE;
- static real hq=-1;
- t_pargs pa[] = {
- { "-type", FALSE, etENUM, {opt},
- "Type of angle" },
- { "-hyd", FALSE, etBOOL, {&bH},
- "Include angles with atoms with mass < 1.5" },
- { "-hq", FALSE, etREAL, {&hq},
- "Ignore angles with atoms with mass < 1.5 and magnitude of their charge less than this value" }
- };
-
- output_env_t oenv;
- FILE *out;
- t_topology *top;
- int i,j,ntype;
- int nft=0,*ft,mult=0;
- int **index;
- int *ft_ind;
- int *nr;
- char **grpnames;
- t_filenm fnm[] = {
- { efTPX, NULL, NULL, ffREAD },
- { efNDX, NULL, "angle", ffWRITE }
- };
-#define NFILE asize(fnm)
-
- CopyRight(stderr,argv[0]);
- parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
- asize(desc),desc,0,NULL,&oenv);
-
-
- ft = select_ftype(opt[0],&nft,&mult);
-
- top = read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
-
- ntype = calc_ntype(nft,ft,&(top->idef));
- snew(grpnames,ntype);
- snew(ft_ind,top->idef.ntypes);
- fill_ft_ind(nft,ft,&top->idef,ft_ind,grpnames);
-
- snew(nr,ntype);
- snew(index,ntype);
- fill_ang(nft,ft,mult,nr,index,ft_ind,top,!bH,hq);
-
- out=ftp2FILE(efNDX,NFILE,fnm,"w");
- for(i=0; (i<ntype); i++) {
- if (nr[i] > 0) {
- fprintf(out,"[ %s ]\n",grpnames[i]);
- for(j=0; (j<nr[i]*mult); j++) {
- fprintf(out," %5d",index[i][j]+1);
- if ((j % 12) == 11)
- fprintf(out,"\n");
- }
- fprintf(out,"\n");
- }
- }
- ffclose(out);
-
- thanx(stderr);
-
- return 0;
+ return gmx_mk_angndx(argc,argv);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff