3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
11 * The make_edi program was generously contributed by Oliver Lange, based
12 * on the code from g_anaeig. You can reach him as olange@gwdg.de.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
32 * Gromacs Runs One Microsecond At Cannonball Speeds
37 #include "gmx_header_config.h"
49 #include "gmx_fatal.h"
68 /* Suppress Cygwin compiler warnings from using newlib version of
78 gmx_bool bConstForce; /* Do constant force flooding instead of
79 evaluating a flooding potential */
88 /* This type is for the average, reference, target, and origin structure */
91 int nr; /* number of atoms this structure contains */
92 int *anrs; /* atom index numbers */
93 rvec *x; /* positions */
94 real *sqrtm; /* sqrt of the masses used for mass-
95 * weighting of analysis */
101 int nini; /* total Nr of atoms */
102 gmx_bool fitmas; /* true if trans fit with cm */
103 gmx_bool pcamas; /* true if mass-weighted PCA */
104 int presteps; /* number of steps to run without any
105 * perturbations ... just monitoring */
106 int outfrq; /* freq (in steps) of writing to edo */
107 int maxedsteps; /* max nr of steps per cycle */
108 struct edix sref; /* reference positions, to these fitting
110 struct edix sav; /* average positions */
111 struct edix star; /* target positions */
112 struct edix sori; /* origin positions */
113 real slope; /* minimal slope in acceptance radexp */
114 int ned; /* Nr of atoms in essdyn buffer */
115 t_edflood flood; /* parameters especially for flooding */
120 void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
126 void write_t_edx(FILE *fp, struct edix edx, const char *comment) {
127 /*here we copy only the pointers into the t_edx struct
128 no data is copied and edx.box is ignored */
130 fprintf(fp,"#%s \n %d \n",comment,edx.nr);
131 for (i=0;i<edx.nr;i++) {
132 fprintf(fp,"%d %f %f %f\n",(edx.anrs)[i]+1,(edx.x)[i][XX],(edx.x)[i][YY],(edx.x)[i][ZZ]);
136 int sscan_list(int *list[], const char *str, const char *listname) {
137 /*this routine scans a string of the form 1,3-6,9 and returns the
138 selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
139 memory for this list will be allocated in this routine -- sscan_list expects *list to
142 listname is a string used in the errormessage*/
147 char *pos,*startpos,*step;
150 /*enums to define the different lexical stati */
151 enum { sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange };
153 int status=sBefore; /*status of the deterministic automat to scan str */
157 char *start=NULL; /*holds the string of the number behind a ','*/
158 char *end=NULL; /*holds the string of the number behind a '-' */
160 int nvecs=0; /* counts the number of vectors in the list*/
172 while ((c=*pos)!=0) {
174 /* expect a number */
175 case sBefore: if (isdigit(c)) {
179 } else status=sError; break;
181 /* have read a number, expect ',' or '-' */
182 case sNumber: if (c==',') {
184 srenew(*list,nvecs+1);
185 (*list)[nvecs++]=number=strtol(start,NULL,10);
190 } else if (c=='-') { status=sMinus; break; }
191 else if (isdigit(c)) break;
192 else status=sError; break;
194 /* have read a '-' -> expect a number */
198 status=sRange; break;
199 } else status=sError; break;
204 status=sError; break;
209 } else status=sError; break;
211 /* have read the number after a minus, expect ',' or ':' */
215 end_number=strtol(end,NULL,10);
216 number=strtol(start,NULL,10);
221 if (end_number<=number) {
222 status=sSmaller; break;
224 srenew(*list,nvecs+end_number-number+1);
226 istep=strtol(step,NULL,10);
229 for (i=number;i<=end_number;i+=istep)
233 status = sSteppedRange;
235 } else if (isdigit(c)) break; else status=sError; break;
237 /* format error occured */
239 gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));
241 /* logical error occured */
243 gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid",listname,pos-startpos);
246 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);
249 ++pos; /* read next character */
250 } /*scanner has finished */
252 /* append zero to list of eigenvectors */
253 srenew(*list,nvecs+1);
259 void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec, const char *grouptitle,real steps[]) {
260 /* eig_list is a zero-terminated list of indices into the eigvecs array.
261 eigvecs are coordinates of eigenvectors
262 grouptitle to write in the comment line
263 steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC
268 while (eig_list[n++]); /*count selected eigenvecs*/
270 fprintf(fp,"# NUMBER OF EIGENVECTORS + %s\n %d\n",grouptitle,n-1);
272 /* write list of eigenvector indicess */
273 for(n=0;eig_list[n];n++) {
275 fprintf(fp,"%8d %g\n",eig_list[n],steps[n]);
277 fprintf(fp,"%8d %g\n",eig_list[n],1.0);
281 /* dump coordinates of the selected eigenvectors */
282 while (eig_list[n]) {
284 for (i=0; i<natoms; i++) {
285 if (eig_list[n]>nvec)
286 gmx_fatal(FARGS,"Selected eigenvector %d is higher than maximum number %d of available eigenvectors",eig_list[n],nvec);
287 copy_rvec(eigvecs[eig_list[n]-1][i],x);
289 fprintf(fp,"%8.5f %8.5f %8.5f\n",x[XX],x[YY],x[ZZ]);
296 /*enum referring to the different lists of eigenvectors*/
297 enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON, evNr };
302 void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
303 int nvec, int *eig_listen[], real* evStepList[])
308 fprintf(fp,"#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
309 MAGIC,edpars->nini,edpars->fitmas,edpars->pcamas);
310 fprintf(fp,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
311 edpars->outfrq,edpars->maxedsteps,edpars->slope);
312 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",
313 edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,
314 edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic,edpars->flood.bConstForce);
316 /* Average and reference positions */
317 write_t_edx(fp,edpars->sref,"NREF, XREF");
318 write_t_edx(fp,edpars->sav,"NAV, XAV");
322 write_eigvec(fp, edpars->ned, eig_listen[evMON],eigvecs,nvec,"COMPONENTS GROUP 1",NULL);
323 write_eigvec(fp, edpars->ned, eig_listen[evLINFIX],eigvecs,nvec,"COMPONENTS GROUP 2",evStepList[evLINFIX]);
324 write_eigvec(fp, edpars->ned, eig_listen[evLINACC],eigvecs,nvec,"COMPONENTS GROUP 3",evStepList[evLINACC]);
325 write_eigvec(fp, edpars->ned, eig_listen[evRADFIX],eigvecs,nvec,"COMPONENTS GROUP 4",evStepList[evRADFIX]);
326 write_eigvec(fp, edpars->ned, eig_listen[evRADACC],eigvecs,nvec,"COMPONENTS GROUP 5",NULL);
327 write_eigvec(fp, edpars->ned, eig_listen[evRADCON],eigvecs,nvec,"COMPONENTS GROUP 6",NULL);
328 write_eigvec(fp, edpars->ned, eig_listen[evFLOOD],eigvecs,nvec,"COMPONENTS GROUP 7",evStepList[evFLOOD]);
331 /*Target and Origin positions */
332 write_t_edx(fp,edpars->star,"NTARGET, XTARGET");
333 write_t_edx(fp,edpars->sori,"NORIGIN, XORIGIN");
336 int read_conffile(const char *confin,char *title,rvec *x[])
338 /* read coordinates out of STX file */
342 printf("read coordnumber from file %s\n",confin);
343 get_stx_coordnum(confin,&natoms);
344 printf("number of coordinates in file %d\n",natoms);
345 /* if (natoms != ncoords)
346 gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
347 " does not match topology (= %d)",
348 confin,natoms,ncoords);
350 /* make space for coordinates and velocities */
351 init_t_atoms(&confat,natoms,FALSE);
353 read_stx_conf(confin,title,&confat,*x,NULL,NULL,box);
358 void read_eigenvalues(int vecs[],const char *eigfile, real values[],
359 gmx_bool bHesse, real kT)
364 neig = read_xvg(eigfile,&eigval,&nrow);
366 fprintf(stderr,"Read %d eigenvalues\n",neig);
367 for(i=bHesse ? 6 : 0 ; i<neig; i++) {
368 if (eigval[1][i] < -0.001 && bHesse)
370 "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n",eigval[1][i]);
372 if (eigval[1][i] < 0)
376 for (i=0; vecs[i]; i++) {
378 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");
379 values[i]=eigval[1][vecs[i]-1]/kT;
382 for (i=0; vecs[i]; i++) {
383 if (vecs[i]>(neig-6))
384 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");
385 values[i]=1/eigval[1][vecs[i]-1];
388 for (i=0; i<nrow; i++)
394 static real *scan_vecparams(const char *str,const char * par, int nvecs)
396 char f0[256],f1[256]; /*format strings adapted every pass of the loop*/
401 snew(vec_params,nvecs);
404 for(i=0; (i<nvecs); i++) {
405 strcpy(f1,f0); /*f0 is the format string for the "to-be-ignored" numbers*/
406 strcat(f1,"%lf"); /*and f1 to read the actual number in this pass of the loop*/
407 if (sscanf(str,f1,&d) != 1)
408 gmx_fatal(FARGS,"Not enough elements for %s parameter (I need %d)",par,nvecs);
418 void init_edx(struct edix *edx) {
424 void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
425 atom_id igro[],rvec *x,const char* structure)
427 /* filter2edx copies coordinates from x to edx which are given in index
433 srenew(edx->x,edx->nr);
434 srenew(edx->anrs,edx->nr);
435 for (i=0;i<nindex;i++,ix++) {
436 for (pos=0; pos<ngro-1 && igro[pos]!=index[i] ; ++pos) {}; /*search element in igro*/
437 if (igro[pos]!=index[i])
438 gmx_fatal(FARGS,"Couldn't find atom with index %d in structure %s",index[i],structure);
439 edx->anrs[ix]=index[i];
440 copy_rvec(x[pos],edx->x[ix]);
444 void get_structure(t_atoms *atoms,const char *IndexFile,
445 const char *StructureFile,struct edix *edx,int nfit,
446 atom_id ifit[],int nav, atom_id index[])
448 atom_id *igro; /*index corresponding to target or origin structure*/
456 ntar=read_conffile(StructureFile,title,&xtar);
457 printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
459 get_index(atoms,IndexFile,1,&ngro,&igro,&grpname);
461 gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
463 filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
465 /* If average and reference/fitting structure differ, append the average structure as well */
466 if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
467 filter2edx(edx,nav,index,ngro,igro,xtar,StructureFile);
470 int gmx_make_edi(int argc,char *argv[])
473 static const char *desc[] = {
474 "[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
475 "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
476 "normal modes analysis ([TT]g_nmeig[tt]).",
477 "ED sampling can be used to manipulate the position along collective coordinates",
478 "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
479 "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
480 "the system to explore new regions along these collective coordinates. A number",
481 "of different algorithms are implemented to drive the system along the eigenvectors",
482 "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
483 "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
484 "or to only monitor the projections of the positions onto",
485 "these coordinates ([TT]-mon[tt]).[PAR]",
487 "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
488 "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
489 "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
490 "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
491 "Towards an exhaustive sampling of the configurational spaces of the ",
492 "two forms of the peptide hormone guanylin,",
493 "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
494 "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
495 "An extended sampling of the configurational space of HPr from E. coli",
496 "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
497 "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
498 "reference structure, target positions, etc.[PAR]",
500 "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
501 "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
502 "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
503 "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
504 "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
505 "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
506 "(steps in the desired direction will be accepted, others will be rejected).",
507 "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
508 "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
509 "to read in a structure file that defines an external origin.[PAR]",
510 "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
511 "towards a target structure specified with [TT]-tar[tt].[PAR]",
512 "NOTE: each eigenvector can be selected only once. [PAR]",
513 "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].edo[tt] file[PAR]",
514 "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
515 "cycle will be started if the spontaneous increase of the radius (in nm/step)",
516 "is less than the value specified.[PAR]",
517 "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
518 "before a new cycle is started.[PAR]",
519 "Note on the parallel implementation: since ED sampling is a 'global' thing",
520 "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
521 "is not very parallel-friendly from an implementation point of view. Because",
522 "parallel ED requires some extra communication, expect the performance to be",
523 "lower as in a free MD simulation, especially on a large number of nodes and/or",
524 "when the ED group contains a lot of atoms. [PAR]",
525 "Please also note that if your ED group contains more than a single protein,",
526 "then the [TT].tpr[tt] file must contain the correct PBC representation of the ED group.",
527 "Take a look on the initial RMSD from the reference structure, which is printed",
528 "out at the start of the simulation; if this is much higher than expected, one",
529 "of the ED molecules might be shifted by a box vector. [PAR]",
530 "All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a [TT].edo[tt] file. In the output",
531 "file, per OUTFRQ step the following information is present: [PAR]",
532 "[TT]*[tt] the step number[BR]",
533 "[TT]*[tt] the number of the ED dataset. ([BB]Note[bb] that you can impose multiple ED constraints in",
534 "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
535 "first. The constraints are applied in the order they appear in the [TT].edi[tt] file.) [BR]",
536 "[TT]*[tt] RMSD (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
537 "* projections of the positions onto selected eigenvectors[BR]",
540 "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
541 "which will lead to extra forces expelling the structure out of the region described",
542 "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
543 "is kept in that region.",
545 "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
546 "It can be changed with [TT]-ori[tt] to an arbitrary position in configuration space.",
547 "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
548 "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
549 "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
550 "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
551 "To use constant Efl set [TT]-tau[tt] to zero.",
553 "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
554 "to give good results for most standard cases in flooding of proteins.",
555 "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
556 "increase, this is mimicked by [GRK]alpha[grk] > 1.",
557 "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
559 "RESTART and FLOODING:",
560 "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
561 "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
564 /* Save all the params in this struct and then save it in an edi file.
565 * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
567 static t_edipar edi_params;
569 enum { evStepNr = evRADFIX + 1};
570 static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
571 static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
572 static const char* evParams[evStepNr] ={NULL,NULL};
573 static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
574 static const char* ConstForceStr;
575 static real* evStepList[evStepNr];
576 static real radfix=0.0;
577 static real deltaF0=150;
578 static real deltaF=0;
580 static real constEfl=0.0;
582 static int eqSteps=0;
583 static int* listen[evNr];
585 const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
586 static gmx_bool bRestrain = FALSE;
587 static gmx_bool bHesse=FALSE;
588 static gmx_bool bHarmonic=FALSE;
590 { "-mon", FALSE, etSTR, {&evSelections[evMON]},
591 "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
592 { "-linfix", FALSE, etSTR, {&evSelections[0]},
593 "Indices of eigenvectors for fixed increment linear sampling" },
594 { "-linacc", FALSE, etSTR, {&evSelections[1]},
595 "Indices of eigenvectors for acceptance linear sampling" },
596 { "-radfix", FALSE, etSTR, {&evSelections[3]},
597 "Indices of eigenvectors for fixed increment radius expansion" },
598 { "-radacc", FALSE, etSTR, {&evSelections[4]},
599 "Indices of eigenvectors for acceptance radius expansion" },
600 { "-radcon", FALSE, etSTR, {&evSelections[5]},
601 "Indices of eigenvectors for acceptance radius contraction" },
602 { "-flood", FALSE, etSTR, {&evSelections[2]},
603 "Indices of eigenvectors for flooding"},
604 { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
605 "Freqency (in steps) of writing output in [TT].edo[tt] file" },
606 { "-slope", FALSE, etREAL, { &edi_params.slope},
607 "Minimal slope in acceptance radius expansion"},
608 { "-linstep", FALSE, etSTR, {&evParams[0]},
609 "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
610 { "-accdir", FALSE, etSTR, {&evParams[1]},
611 "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
612 { "-radstep", FALSE, etREAL, {&radfix},
613 "Stepsize (nm/step) for fixed increment radius expansion"},
614 { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
615 "Maximum number of steps per cycle" },
616 { "-eqsteps", FALSE, etINT, {&eqSteps},
617 "Number of steps to run without any perturbations "},
618 { "-deltaF0", FALSE,etREAL, {&deltaF0},
619 "Target destabilization energy for flooding"},
620 { "-deltaF", FALSE, etREAL, {&deltaF},
621 "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
622 { "-tau", FALSE, etREAL, {&tau},
623 "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
624 { "-Eflnull", FALSE, etREAL, {&constEfl},
625 "The starting value of the flooding strength. The flooding strength is updated "
626 "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
627 { "-T", FALSE, etREAL, {&T},
628 "T is temperature, the value is needed if you want to do flooding "},
629 { "-alpha",FALSE,etREAL,{&alpha},
630 "Scale width of gaussian flooding potential with alpha^2 "},
631 { "-restrain",FALSE, etBOOL, {&bRestrain},
632 "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
633 { "-hessian",FALSE, etBOOL, {&bHesse},
634 "The eigenvectors and eigenvalues are from a Hessian matrix"},
635 { "-harmonic",FALSE, etBOOL, {&bHarmonic},
636 "The eigenvalues are interpreted as spring constant"},
637 { "-constF", FALSE, etSTR, {&ConstForceStr},
638 "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
639 "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
641 #define NPA asize(pa)
644 int nvec1,*eignr1=NULL;
645 rvec *xav1,**eigvec1=NULL;
647 int nav; /* Number of atoms in the average structure */
649 const char *indexfile;
651 atom_id *index,*ifit;
652 int nfit; /* Number of atoms in the reference/fit structure */
653 int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
655 real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */
658 const char *TargetFile;
659 const char *OriginFile;
660 const char *EigvecFile;
664 /*to read topology file*/
670 gmx_bool bTop, bFit1;
673 { efTRN, "-f", "eigenvec", ffREAD },
674 { efXVG, "-eig", "eigenval", ffOPTRD },
675 { efTPS, NULL, NULL, ffREAD },
676 { efNDX, NULL, NULL, ffOPTRD },
677 { efSTX, "-tar", "target", ffOPTRD},
678 { efSTX, "-ori", "origin", ffOPTRD},
679 { efEDI, "-o", "sam", ffWRITE }
681 #define NFILE asize(fnm)
682 edi_params.outfrq=100; edi_params.slope=0.0; edi_params.maxedsteps=0;
683 CopyRight(stderr,argv[0]);
684 parse_common_args(&argc,argv, 0 ,
685 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
687 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
688 EdiFile=ftp2fn(efEDI,NFILE,fnm);
689 TargetFile = opt2fn_null("-tar",NFILE,fnm);
690 OriginFile = opt2fn_null("-ori",NFILE,fnm);
693 for (ev_class=0; ev_class<evNr; ++ev_class) {
694 if (opt2parg_bSet(evOptions[ev_class],NPA,pa)) {
695 /*get list of eigenvectors*/
696 nvecs=sscan_list(&(listen[ev_class]),opt2parg_str(evOptions[ev_class],NPA,pa),evOptions[ev_class]);
697 if (ev_class<evStepNr-2) {
698 /*if apropriate get list of stepsizes for these eigenvectors*/
699 if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
700 evStepList[ev_class]=
701 scan_vecparams(opt2parg_str(evStepOptions[ev_class],NPA,pa),evStepOptions[ev_class],nvecs);
702 else { /*if list is not given fill with zeros */
703 snew(evStepList[ev_class],nvecs);
704 for (i=0; i<nvecs; i++)
705 evStepList[ev_class][i]=0.0;
707 } else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
708 snew(evStepList[ev_class],nvecs);
709 for (i=0; i<nvecs; i++)
710 evStepList[ev_class][i]=radfix;
711 } else if (ev_class == evFLOOD) {
712 snew(evStepList[ev_class],nvecs);
714 /* Are we doing constant force flooding? In that case, we read in
715 * the fproj values from the command line */
716 if (opt2parg_bSet("-constF", NPA, pa))
718 evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
720 } else {}; /*to avoid ambiguity */
721 } else { /* if there are no eigenvectors for this option set list to zero */
722 listen[ev_class]=NULL;
723 snew(listen[ev_class],1);
724 listen[ev_class][0]=0;
728 /* print the interpreted list of eigenvectors - to give some feedback*/
729 for (ev_class=0; ev_class<evNr; ++ev_class) {
730 printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
732 while(listen[ev_class][i])
733 printf("%d ",listen[ev_class][i++]);
738 EigvecFile=opt2fn("-f",NFILE,fnm);
740 /*read eigenvectors from eigvec.trr*/
741 read_eigenvectors(EigvecFile,&nav,&bFit1,
742 &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);
744 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
745 title,&top,&ePBC,&xtop,NULL,topbox,0);
749 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",nav);
750 get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
752 gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
762 /* if g_covar used different coordinate groups to fit and to do the PCA */
763 printf("\nNote: the structure in %s should be the same\n"
764 " as the one used for the fit in g_covar\n",ftp2fn(efTPS,NFILE,fnm));
765 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
769 printf("\nNote: Apparently no fitting was done in g_covar.\n"
770 " However, you need to select a reference group for fitting in mdrun\n");
772 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
775 copy_rvec(xtop[ifit[i]],xref1[i]);
783 if (opt2parg_bSet("-constF", NPA, pa))
785 /* Constant force flooding is special: Most of the normal flooding
786 * options are not needed. */
787 edi_params.flood.bConstForce = TRUE;
791 /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
793 if ( listen[evFLOOD][0]!=0)
794 read_eigenvalues(listen[evFLOOD],opt2fn("-eig",NFILE,fnm),evStepList[evFLOOD],bHesse,kB*T);
796 edi_params.flood.tau=tau;
797 edi_params.flood.deltaF0=deltaF0;
798 edi_params.flood.deltaF=deltaF;
799 edi_params.presteps=eqSteps;
800 edi_params.flood.kT=kB*T;
801 edi_params.flood.bHarmonic=bHarmonic;
804 /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
805 edi_params.flood.constEfl=-constEfl;
806 edi_params.flood.alpha2=-sqr(alpha);
810 edi_params.flood.constEfl=constEfl;
811 edi_params.flood.alpha2=sqr(alpha);
817 /*number of system atoms */
818 edi_params.nini=atoms->nr;
821 /*store reference and average structure in edi_params*/
822 make_t_edx(&edi_params.sref,nfit,xref1,ifit );
823 make_t_edx(&edi_params.sav ,nav ,xav1 ,index);
826 /* Store target positions in edi_params */
827 if (opt2bSet("-tar",NFILE,fnm))
829 if (0 != listen[evFLOOD][0])
831 fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
832 " You may want to use -ori to define the flooding potential center.\n\n");
834 get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,nav,index);
838 make_t_edx(&edi_params.star,0,NULL,index);
841 /* Store origin positions */
842 if (opt2bSet("-ori",NFILE,fnm))
844 get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,nav,index);
848 make_t_edx(&edi_params.sori,0,NULL,index);
852 write_the_whole_thing(ffopen(EdiFile,"w"), &edi_params, eigvec1,nvec1, listen, evStepList);