/* This type is for the average, reference, target, and origin structure */
typedef struct edix
{
- int nr; /* number of atoms this structure contains */
+ 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-
-void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
- edx->nr=natoms;
- edx->anrs=index;
- edx->x=pos;
+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]);
- }
+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;
+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[]) {
+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*/
+ 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);
+ }
+ 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]);
+ /* 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++;
}
- n++;
- }
}
/*enum referring to the different lists of eigenvectors*/
-enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON, evNr };
+enum {
+ evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON, evMON, evNr
+};
#define oldMAGIC 666
#define MAGIC 670
/* 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);
+ 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");
+ 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]);
+ 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");
+ 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[])
+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);
+ 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 {*/
+ 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);
+ 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[],
+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;
+ 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];
+ }
}
- 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]);
}
- /* free memory */
- for (i=0; i<nrow; i++)
- sfree(eigval[i]);
- sfree(eigval);
+ sfree(eigval);
}
-static real *scan_vecparams(const char *str,const char * par, int nvecs)
+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");
+ 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;
+ return vec_params;
}
-void init_edx(struct edix *edx) {
- edx->nr=0;
- snew(edx->x,1);
- snew(edx->anrs,1);
+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)
+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]);
- }
+ */
+
+ 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[])
+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);
+ 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[])
+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 analysis ([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].xvg[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 implementation 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 and/or",
- "when the ED group contains a lot of atoms. [PAR]",
- "Please also note that if your ED group contains more than a single protein,",
- "then the [TT].tpr[tt] file must contain the correct PBC representation of the ED group.",
- "Take a look on the initial RMSD from the reference structure, which is printed",
- "out at the start of the simulation; if this is much higher than expected, one",
- "of the ED molecules might be shifted by a box vector. [PAR]",
- "All ED-related output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a [TT].xvg[tt] file",
- "as a function of time in intervals of OUTFRQ steps.[PAR]",
- "[BB]Note[bb] that you can impose multiple ED constraints and flooding potentials 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. ",
- "Depending on what was specified in the [TT].edi[tt] input file, the output file contains for each ED dataset[PAR]",
- "[TT]*[tt] the RMSD of the fitted molecule to the reference structure (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
- "[TT]*[tt] 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 configuration 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."
- };
+ 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 analysis ([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].xvg[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 implementation 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 and/or",
+ "when the ED group contains a lot of atoms. [PAR]",
+ "Please also note that if your ED group contains more than a single protein,",
+ "then the [TT].tpr[tt] file must contain the correct PBC representation of the ED group.",
+ "Take a look on the initial RMSD from the reference structure, which is printed",
+ "out at the start of the simulation; if this is much higher than expected, one",
+ "of the ED molecules might be shifted by a box vector. [PAR]",
+ "All ED-related output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a [TT].xvg[tt] file",
+ "as a function of time in intervals of OUTFRQ steps.[PAR]",
+ "[BB]Note[bb] that you can impose multiple ED constraints and flooding potentials 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. ",
+ "Depending on what was specified in the [TT].edi[tt] input file, the output file contains for each ED dataset[PAR]",
+ "[TT]*[tt] the RMSD of the fitted molecule to the reference structure (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
+ "[TT]*[tt] 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 configuration 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
- */
+ * 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"};
+ 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].xvg[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."}
+ 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].xvg[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;
+ 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;
+ 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 }
+ 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;
- parse_common_args(&argc,argv, 0 ,
- NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
+ edi_params.outfrq = 100; edi_params.slope = 0.0; edi_params.maxedsteps = 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);
+ 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)) {
+ 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) {
+ 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;
+ 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 (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);
+ 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);
+ 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;
+ }
+ 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++]);
+ 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);
+ 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);
+ 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;
+ 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("\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 (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));
+ " 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");
+ " 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]);
}
- 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;
+ nfit = nav;
+ ifit = index;
}
if (opt2parg_bSet("-constF", NPA, pa))
{
/* 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);
+ 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;
+ 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);
+ 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.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;
-}
+ 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;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff