--- /dev/null
- static real radfix = 0.0;
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS 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.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, 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 http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/* 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. He
+ * probably also holds copyright to the following code.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "readinp.h"
+#include "gromacs/commandline/pargs.h"
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "gmx_fatal.h"
+#include "vec.h"
+#include "pbc.h"
+#include "gromacs/fileio/futil.h"
+#include "gromacs/fileio/pdbio.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/matio.h"
+#include "mshift.h"
+#include "xvgr.h"
+#include "do_fit.h"
+#include "rmpbc.h"
+#include "txtdump.h"
+#include "eigio.h"
+#include "index.h"
+#include "string2.h"
+
+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[] = {
+ "[THISMODULE] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
+ "based on eigenvectors of a covariance matrix ([gmx-covar]) or from a",
+ "normal modes analysis ([gmx-nmeig]).",
+ "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
+ */
+ 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];
- { "-radstep", FALSE, etREAL, {&radfix},
++ static real radstep = 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\")"},
- else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class], NPA, pa))
++ { "-radstep", FALSE, etREAL, {&radstep},
+ "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;
+ if (!parse_common_args(&argc, argv, 0,
+ NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+
+ 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;
+ }
+ }
+ }
- evStepList[ev_class][i] = radfix;
++ else if (ev_class == evRADFIX)
+ {
+ snew(evStepList[ev_class], nvecs);
+ for (i = 0; i < nvecs; i++)
+ {
++ evStepList[ev_class][i] = radstep;
+ }
+ }
+ 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);
+
+ return 0;
+}
--- /dev/null
- static int get_nthreads_hw_avail(const t_commrec gmx_unused *cr)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS 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.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, 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 http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+#ifdef HAVE_UNISTD_H
+/* For sysconf */
+#include <unistd.h>
+#endif
+
+#include "types/enums.h"
+#include "types/hw_info.h"
+#include "types/commrec.h"
+#include "gmx_fatal.h"
+#include "gmx_fatal_collective.h"
+#include "smalloc.h"
+#include "gpu_utils.h"
+#include "copyrite.h"
+#include "gmx_detect_hardware.h"
+#include "main.h"
+#include "md_logging.h"
++#include "gromacs/utility/gmxomp.h"
+
+#include "thread_mpi/threads.h"
+
+#if ((defined(WIN32) || defined( _WIN32 ) || defined(WIN64) || defined( _WIN64 )) && !(defined (__CYGWIN__) || defined (__CYGWIN32__)))
+#include "windows.h"
+#endif
+
+#ifdef GMX_GPU
+const gmx_bool bGPUBinary = TRUE;
+#else
+const gmx_bool bGPUBinary = FALSE;
+#endif
+
+static const char * invalid_gpuid_hint =
+ "A delimiter-free sequence of valid numeric IDs of available GPUs is expected.";
+
+/* The globally shared hwinfo structure. */
+static gmx_hw_info_t *hwinfo_g;
+/* A reference counter for the hwinfo structure */
+static int n_hwinfo = 0;
+/* A lock to protect the hwinfo structure */
+static tMPI_Thread_mutex_t hw_info_lock = TMPI_THREAD_MUTEX_INITIALIZER;
+
+
+/* FW decl. */
+static void limit_num_gpus_used(gmx_gpu_opt_t *gpu_opt, int count);
+static int gmx_count_gpu_dev_unique(const gmx_gpu_info_t *gpu_info,
+ const gmx_gpu_opt_t *gpu_opt);
+
+static void sprint_gpus(char *sbuf, const gmx_gpu_info_t *gpu_info)
+{
+ int i, ndev;
+ char stmp[STRLEN];
+
+ ndev = gpu_info->ncuda_dev;
+
+ sbuf[0] = '\0';
+ for (i = 0; i < ndev; i++)
+ {
+ get_gpu_device_info_string(stmp, gpu_info, i);
+ strcat(sbuf, " ");
+ strcat(sbuf, stmp);
+ if (i < ndev - 1)
+ {
+ strcat(sbuf, "\n");
+ }
+ }
+}
+
+static void print_gpu_detection_stats(FILE *fplog,
+ const gmx_gpu_info_t *gpu_info,
+ const t_commrec *cr)
+{
+ char onhost[266], stmp[STRLEN];
+ int ngpu;
+
+ if (!gpu_info->bDetectGPUs)
+ {
+ /* We skipped the detection, so don't print detection stats */
+ return;
+ }
+
+ ngpu = gpu_info->ncuda_dev;
+
+#if defined GMX_MPI && !defined GMX_THREAD_MPI
+ /* We only print the detection on one, of possibly multiple, nodes */
+ strncpy(onhost, " on host ", 10);
+ gmx_gethostname(onhost+9, 256);
+#else
+ /* We detect all relevant GPUs */
+ strncpy(onhost, "", 1);
+#endif
+
+ if (ngpu > 0)
+ {
+ sprint_gpus(stmp, gpu_info);
+ md_print_warn(cr, fplog, "%d GPU%s detected%s:\n%s\n",
+ ngpu, (ngpu > 1) ? "s" : "", onhost, stmp);
+ }
+ else
+ {
+ md_print_warn(cr, fplog, "No GPUs detected%s\n", onhost);
+ }
+}
+
+static void print_gpu_use_stats(FILE *fplog,
+ const gmx_gpu_info_t *gpu_info,
+ const gmx_gpu_opt_t *gpu_opt,
+ const t_commrec *cr)
+{
+ char sbuf[STRLEN], stmp[STRLEN];
+ int i, ngpu_comp, ngpu_use;
+
+ ngpu_comp = gpu_info->ncuda_dev_compatible;
+ ngpu_use = gpu_opt->ncuda_dev_use;
+
+ /* Issue a note if GPUs are available but not used */
+ if (ngpu_comp > 0 && ngpu_use < 1)
+ {
+ sprintf(sbuf,
+ "%d compatible GPU%s detected in the system, but none will be used.\n"
+ "Consider trying GPU acceleration with the Verlet scheme!",
+ ngpu_comp, (ngpu_comp > 1) ? "s" : "");
+ }
+ else
+ {
+ int ngpu_use_uniq;
+
+ ngpu_use_uniq = gmx_count_gpu_dev_unique(gpu_info, gpu_opt);
+
+ sprintf(sbuf, "%d GPU%s %sselected for this run.\n"
+ "Mapping of GPU%s to the %d PP rank%s in this node: ",
+ ngpu_use_uniq, (ngpu_use_uniq > 1) ? "s" : "",
+ gpu_opt->bUserSet ? "user-" : "auto-",
+ (ngpu_use > 1) ? "s" : "",
+ cr->nrank_pp_intranode,
+ (cr->nrank_pp_intranode > 1) ? "s" : "");
+
+ for (i = 0; i < ngpu_use; i++)
+ {
+ sprintf(stmp, "#%d", get_gpu_device_id(gpu_info, gpu_opt, i));
+ if (i < ngpu_use - 1)
+ {
+ strcat(stmp, ", ");
+ }
+ strcat(sbuf, stmp);
+ }
+ }
+ md_print_info(cr, fplog, "%s\n\n", sbuf);
+}
+
+/* Parse a "plain" GPU ID string which contains a sequence of digits corresponding
+ * to GPU IDs; the order will indicate the process/tMPI thread - GPU assignment. */
+static void parse_gpu_id_plain_string(const char *idstr, int *nid, int **idlist)
+{
+ int i;
+
+ *nid = strlen(idstr);
+
+ snew(*idlist, *nid);
+
+ for (i = 0; i < *nid; i++)
+ {
+ if (idstr[i] < '0' || idstr[i] > '9')
+ {
+ gmx_fatal(FARGS, "Invalid character in GPU ID string: '%c'\n%s\n",
+ idstr[i], invalid_gpuid_hint);
+ }
+ (*idlist)[i] = idstr[i] - '0';
+ }
+}
+
+static void parse_gpu_id_csv_string(const char gmx_unused *idstr, int gmx_unused *nid, int gmx_unused *idlist)
+{
+ /* XXX implement cvs format to support more than 10 different GPUs in a box. */
+ gmx_incons("Not implemented yet");
+}
+
+void gmx_check_hw_runconf_consistency(FILE *fplog,
+ const gmx_hw_info_t *hwinfo,
+ const t_commrec *cr,
+ const gmx_hw_opt_t *hw_opt,
+ gmx_bool bUseGPU)
+{
+ int npppn, ntmpi_pp;
+ char sbuf[STRLEN], th_or_proc[STRLEN], th_or_proc_plural[STRLEN], pernode[STRLEN];
+ gmx_bool btMPI, bMPI, bMaxMpiThreadsSet, bNthreadsAuto, bEmulateGPU;
+
+ assert(hwinfo);
+ assert(cr);
+
+ /* Below we only do consistency checks for PP and GPUs,
+ * this is irrelevant for PME only nodes, so in that case we return
+ * here.
+ */
+ if (!(cr->duty & DUTY_PP))
+ {
+ return;
+ }
+
+ btMPI = bMPI = FALSE;
+ bNthreadsAuto = FALSE;
+#if defined(GMX_THREAD_MPI)
+ btMPI = TRUE;
+ bNthreadsAuto = (hw_opt->nthreads_tmpi < 1);
+#elif defined(GMX_LIB_MPI)
+ bMPI = TRUE;
+#endif
+
+ /* GPU emulation detection is done later, but we need here as well
+ * -- uncool, but there's no elegant workaround */
+ bEmulateGPU = (getenv("GMX_EMULATE_GPU") != NULL);
+ bMaxMpiThreadsSet = (getenv("GMX_MAX_MPI_THREADS") != NULL);
+
+ /* check the acceleration mdrun is compiled with against hardware
+ capabilities */
+ /* TODO: Here we assume homogeneous hardware which is not necessarily
+ the case! Might not hurt to add an extra check over MPI. */
+ gmx_cpuid_acceleration_check(hwinfo->cpuid_info, fplog, SIMMASTER(cr));
+
+ /* NOTE: this print is only for and on one physical node */
+ print_gpu_detection_stats(fplog, &hwinfo->gpu_info, cr);
+
+ if (hwinfo->gpu_info.ncuda_dev_compatible > 0)
+ {
+ /* NOTE: this print is only for and on one physical node */
+ print_gpu_use_stats(fplog, &hwinfo->gpu_info, &hw_opt->gpu_opt, cr);
+ }
+
+ /* Need to ensure that we have enough GPUs:
+ * - need one GPU per PP node
+ * - no GPU oversubscription with tMPI
+ * */
+ /* number of PP processes per node */
+ npppn = cr->nrank_pp_intranode;
+
+ pernode[0] = '\0';
+ th_or_proc_plural[0] = '\0';
+ if (btMPI)
+ {
+ sprintf(th_or_proc, "thread-MPI thread");
+ if (npppn > 1)
+ {
+ sprintf(th_or_proc_plural, "s");
+ }
+ }
+ else if (bMPI)
+ {
+ sprintf(th_or_proc, "MPI process");
+ if (npppn > 1)
+ {
+ sprintf(th_or_proc_plural, "es");
+ }
+ sprintf(pernode, " per node");
+ }
+ else
+ {
+ /* neither MPI nor tMPI */
+ sprintf(th_or_proc, "process");
+ }
+
+ if (bUseGPU && hwinfo->gpu_info.ncuda_dev_compatible > 0 &&
+ !bEmulateGPU)
+ {
+ int ngpu_comp, ngpu_use;
+ char gpu_comp_plural[2], gpu_use_plural[2];
+
+ ngpu_comp = hwinfo->gpu_info.ncuda_dev_compatible;
+ ngpu_use = hw_opt->gpu_opt.ncuda_dev_use;
+
+ sprintf(gpu_comp_plural, "%s", (ngpu_comp> 1) ? "s" : "");
+ sprintf(gpu_use_plural, "%s", (ngpu_use > 1) ? "s" : "");
+
+ /* number of tMPI threads auto-adjusted */
+ if (btMPI && bNthreadsAuto)
+ {
+ if (hw_opt->gpu_opt.bUserSet && npppn < ngpu_use)
+ {
+ /* The user manually provided more GPUs than threads we
+ could automatically start. */
+ gmx_fatal(FARGS,
+ "%d GPU%s provided, but only %d PP thread-MPI thread%s coud be started.\n"
+ "%s requires one PP tread-MPI thread per GPU; use fewer GPUs%s.",
+ ngpu_use, gpu_use_plural,
+ npppn, th_or_proc_plural,
+ ShortProgram(), bMaxMpiThreadsSet ? "\nor allow more threads to be used" : "");
+ }
+
+ if (!hw_opt->gpu_opt.bUserSet && npppn < ngpu_comp)
+ {
+ /* There are more GPUs than tMPI threads; we have
+ limited the number GPUs used. */
+ md_print_warn(cr, fplog,
+ "NOTE: %d GPU%s were detected, but only %d PP thread-MPI thread%s can be started.\n"
+ " %s can use one GPU per PP tread-MPI thread, so only %d GPU%s will be used.%s\n",
+ ngpu_comp, gpu_comp_plural,
+ npppn, th_or_proc_plural,
+ ShortProgram(), npppn,
+ npppn > 1 ? "s" : "",
+ bMaxMpiThreadsSet ? "\n Also, you can allow more threads to be used by increasing GMX_MAX_MPI_THREADS" : "");
+ }
+ }
+
+ if (hw_opt->gpu_opt.bUserSet)
+ {
+ if (ngpu_use != npppn)
+ {
+ gmx_fatal(FARGS,
+ "Incorrect launch configuration: mismatching number of PP %s%s and GPUs%s.\n"
+ "%s was started with %d PP %s%s%s, but you provided %d GPU%s.",
+ th_or_proc, btMPI ? "s" : "es", pernode,
+ ShortProgram(), npppn, th_or_proc,
+ th_or_proc_plural, pernode,
+ ngpu_use, gpu_use_plural);
+ }
+ }
+ else
+ {
+ if (ngpu_comp > npppn)
+ {
+ md_print_warn(cr, fplog,
+ "NOTE: potentially sub-optimal launch configuration, %s started with less\n"
+ " PP %s%s%s than GPU%s available.\n"
+ " Each PP %s can use only one GPU, %d GPU%s%s will be used.\n",
+ ShortProgram(), th_or_proc,
+ th_or_proc_plural, pernode, gpu_comp_plural,
+ th_or_proc, npppn, gpu_use_plural, pernode);
+ }
+
+ if (ngpu_use != npppn)
+ {
+ /* Avoid duplicate error messages.
+ * Unfortunately we can only do this at the physical node
+ * level, since the hardware setup and MPI process count
+ * might differ between physical nodes.
+ */
+ if (cr->rank_pp_intranode == 0)
+ {
+ gmx_fatal(FARGS,
+ "Incorrect launch configuration: mismatching number of PP %s%s and GPUs%s.\n"
+ "%s was started with %d PP %s%s%s, but only %d GPU%s were detected.",
+ th_or_proc, btMPI ? "s" : "es", pernode,
+ ShortProgram(), npppn, th_or_proc,
+ th_or_proc_plural, pernode,
+ ngpu_use, gpu_use_plural);
+ }
+ }
+ }
+
+ {
+ int same_count;
+
+ same_count = gmx_count_gpu_dev_shared(&hw_opt->gpu_opt);
+
+ if (same_count > 0)
+ {
+ md_print_info(cr, fplog,
+ "NOTE: You assigned %s to multiple %s%s.\n",
+ same_count > 1 ? "GPUs" : "a GPU", th_or_proc, btMPI ? "s" : "es");
+ }
+ }
+ }
+
+#ifdef GMX_MPI
+ if (PAR(cr))
+ {
+ /* Avoid other ranks to continue after
+ inconsistency */
+ MPI_Barrier(cr->mpi_comm_mygroup);
+ }
+#endif
+
+}
+
+/* Return 0 if none of the GPU (per node) are shared among PP ranks.
+ *
+ * Sharing GPUs among multiple PP ranks is possible when the user passes
+ * GPU IDs. Here we check for sharing and return a non-zero value when
+ * this is detected. Note that the return value represents the number of
+ * PP rank pairs that share a device.
+ */
+int gmx_count_gpu_dev_shared(const gmx_gpu_opt_t *gpu_opt)
+{
+ int same_count = 0;
+ int ngpu = gpu_opt->ncuda_dev_use;
+
+ if (gpu_opt->bUserSet)
+ {
+ int i, j;
+
+ for (i = 0; i < ngpu - 1; i++)
+ {
+ for (j = i + 1; j < ngpu; j++)
+ {
+ same_count += (gpu_opt->cuda_dev_use[i] ==
+ gpu_opt->cuda_dev_use[j]);
+ }
+ }
+ }
+
+ return same_count;
+}
+
+/* Count and return the number of unique GPUs (per node) selected.
+ *
+ * As sharing GPUs among multiple PP ranks is possible when the user passes
+ * GPU IDs, the number of GPUs user (per node) can be different from the
+ * number of GPU IDs selected.
+ */
+static int gmx_count_gpu_dev_unique(const gmx_gpu_info_t *gpu_info,
+ const gmx_gpu_opt_t *gpu_opt)
+{
+ int i, uniq_count, ngpu;
+ int *uniq_ids;
+
+ assert(gpu_info);
+ assert(gpu_opt);
+
+ ngpu = gpu_info->ncuda_dev;
+ uniq_count = 0;
+
+ snew(uniq_ids, ngpu);
+
+ /* Each element in uniq_ids will be set to 0 or 1. The n-th element set
+ * to 1 indicates that the respective GPU was selected to be used. */
+ for (i = 0; i < gpu_opt->ncuda_dev_use; i++)
+ {
+ uniq_ids[get_gpu_device_id(gpu_info, gpu_opt, i)] = 1;
+ }
+ /* Count the devices used. */
+ for (i = 0; i < ngpu; i++)
+ {
+ uniq_count += uniq_ids[i];
+ }
+
+ sfree(uniq_ids);
+
+ return uniq_count;
+}
+
+
+/* Return the number of hardware threads supported by the current CPU.
+ * We assume that this is equal with the number of CPUs reported to be
+ * online by the OS at the time of the call.
+ */
- #ifdef GMX_OMPENMP
++static int get_nthreads_hw_avail(FILE gmx_unused *fplog, const t_commrec gmx_unused *cr)
+{
+ int ret = 0;
+
+#if ((defined(WIN32) || defined( _WIN32 ) || defined(WIN64) || defined( _WIN64 )) && !(defined (__CYGWIN__) || defined (__CYGWIN32__)))
+ /* Windows */
+ SYSTEM_INFO sysinfo;
+ GetSystemInfo( &sysinfo );
+ ret = sysinfo.dwNumberOfProcessors;
+#elif defined HAVE_SYSCONF
+ /* We are probably on Unix.
+ * Now check if we have the argument to use before executing the call
+ */
+#if defined(_SC_NPROCESSORS_ONLN)
+ ret = sysconf(_SC_NPROCESSORS_ONLN);
+#elif defined(_SC_NPROC_ONLN)
+ ret = sysconf(_SC_NPROC_ONLN);
+#elif defined(_SC_NPROCESSORS_CONF)
+ ret = sysconf(_SC_NPROCESSORS_CONF);
+#elif defined(_SC_NPROC_CONF)
+ ret = sysconf(_SC_NPROC_CONF);
+#else
+#warning "No valid sysconf argument value found. Executables will not be able to determine the number of hardware threads: mdrun will use 1 thread by default!"
+#endif /* End of check for sysconf argument values */
+
+#else
+ /* Neither windows nor Unix. No fscking idea how many CPUs we have! */
+ ret = -1;
+#endif
+
+ if (debug)
+ {
+ fprintf(debug, "Detected %d processors, will use this as the number "
+ "of supported hardware threads.\n", ret);
+ }
+
- hwinfo_g->nthreads_hw_avail = get_nthreads_hw_avail(cr);
++#ifdef GMX_OPENMP
+ if (ret != gmx_omp_get_num_procs())
+ {
+ md_print_warn(cr, fplog,
+ "Number of CPUs detected (%d) does not match the number reported by OpenMP (%d).\n"
+ "Consider setting the launch configuration manually!",
+ ret, gmx_omp_get_num_procs());
+ }
+#endif
+
+ return ret;
+}
+
+static void gmx_detect_gpus(FILE *fplog, const t_commrec *cr)
+{
+#ifdef GMX_LIB_MPI
+ int rank_world;
+ MPI_Comm physicalnode_comm;
+#endif
+ int rank_local;
+
+ /* Under certain circumstances MPI ranks on the same physical node
+ * can not simultaneously access the same GPU(s). Therefore we run
+ * the detection only on one MPI rank per node and broadcast the info.
+ * Note that with thread-MPI only a single thread runs this code.
+ *
+ * TODO: We should also do CPU hardware detection only once on each
+ * physical node and broadcast it, instead of do it on every MPI rank.
+ */
+#ifdef GMX_LIB_MPI
+ /* A split of MPI_COMM_WORLD over physical nodes is only required here,
+ * so we create and destroy it locally.
+ */
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank_world);
+ MPI_Comm_split(MPI_COMM_WORLD, gmx_physicalnode_id_hash(),
+ rank_world, &physicalnode_comm);
+ MPI_Comm_rank(physicalnode_comm, &rank_local);
+#else
+ /* Here there should be only one process, check this */
+ assert(cr->nnodes == 1 && cr->sim_nodeid == 0);
+
+ rank_local = 0;
+#endif
+
+ if (rank_local == 0)
+ {
+ char detection_error[STRLEN], sbuf[STRLEN];
+
+ if (detect_cuda_gpus(&hwinfo_g->gpu_info, detection_error) != 0)
+ {
+ if (detection_error != NULL && detection_error[0] != '\0')
+ {
+ sprintf(sbuf, ":\n %s\n", detection_error);
+ }
+ else
+ {
+ sprintf(sbuf, ".");
+ }
+ md_print_warn(cr, fplog,
+ "NOTE: Error occurred during GPU detection%s"
+ " Can not use GPU acceleration, will fall back to CPU kernels.\n",
+ sbuf);
+ }
+ }
+
+#ifdef GMX_LIB_MPI
+ /* Broadcast the GPU info to the other ranks within this node */
+ MPI_Bcast(&hwinfo_g->gpu_info.ncuda_dev, 1, MPI_INT, 0, physicalnode_comm);
+
+ if (hwinfo_g->gpu_info.ncuda_dev > 0)
+ {
+ int cuda_dev_size;
+
+ cuda_dev_size = hwinfo_g->gpu_info.ncuda_dev*sizeof_cuda_dev_info();
+
+ if (rank_local > 0)
+ {
+ hwinfo_g->gpu_info.cuda_dev =
+ (cuda_dev_info_ptr_t)malloc(cuda_dev_size);
+ }
+ MPI_Bcast(hwinfo_g->gpu_info.cuda_dev, cuda_dev_size, MPI_BYTE,
+ 0, physicalnode_comm);
+ MPI_Bcast(&hwinfo_g->gpu_info.ncuda_dev_compatible, 1, MPI_INT,
+ 0, physicalnode_comm);
+ }
+
+ MPI_Comm_free(&physicalnode_comm);
+#endif
+}
+
+gmx_hw_info_t *gmx_detect_hardware(FILE *fplog, const t_commrec *cr,
+ gmx_bool bDetectGPUs)
+{
+ gmx_hw_info_t *hw;
+ int ret;
+
+ /* make sure no one else is doing the same thing */
+ ret = tMPI_Thread_mutex_lock(&hw_info_lock);
+ if (ret != 0)
+ {
+ gmx_fatal(FARGS, "Error locking hwinfo mutex: %s", strerror(errno));
+ }
+
+ /* only initialize the hwinfo structure if it is not already initalized */
+ if (n_hwinfo == 0)
+ {
+ snew(hwinfo_g, 1);
+
+ /* detect CPUID info; no fuss, we don't detect system-wide
+ * -- sloppy, but that's it for now */
+ if (gmx_cpuid_init(&hwinfo_g->cpuid_info) != 0)
+ {
+ gmx_fatal_collective(FARGS, cr, NULL, "CPUID detection failed!");
+ }
+
+ /* detect number of hardware threads */
++ hwinfo_g->nthreads_hw_avail = get_nthreads_hw_avail(fplog, cr);
+
+ /* detect GPUs */
+ hwinfo_g->gpu_info.ncuda_dev = 0;
+ hwinfo_g->gpu_info.cuda_dev = NULL;
+ hwinfo_g->gpu_info.ncuda_dev_compatible = 0;
+
+ /* Run the detection if the binary was compiled with GPU support
+ * and we requested detection.
+ */
+ hwinfo_g->gpu_info.bDetectGPUs =
+ (bGPUBinary && bDetectGPUs &&
+ getenv("GMX_DISABLE_GPU_DETECTION") == NULL);
+ if (hwinfo_g->gpu_info.bDetectGPUs)
+ {
+ gmx_detect_gpus(fplog, cr);
+ }
+ }
+ /* increase the reference counter */
+ n_hwinfo++;
+
+ ret = tMPI_Thread_mutex_unlock(&hw_info_lock);
+ if (ret != 0)
+ {
+ gmx_fatal(FARGS, "Error unlocking hwinfo mutex: %s", strerror(errno));
+ }
+
+ return hwinfo_g;
+}
+
+void gmx_parse_gpu_ids(gmx_gpu_opt_t *gpu_opt)
+{
+ char *env;
+
+ if (gpu_opt->gpu_id != NULL && !bGPUBinary)
+ {
+ gmx_fatal(FARGS, "GPU ID string set, but %s was compiled without GPU support!", ShortProgram());
+ }
+
+ env = getenv("GMX_GPU_ID");
+ if (env != NULL && gpu_opt->gpu_id != NULL)
+ {
+ gmx_fatal(FARGS, "GMX_GPU_ID and -gpu_id can not be used at the same time");
+ }
+ if (env == NULL)
+ {
+ env = gpu_opt->gpu_id;
+ }
+
+ /* parse GPU IDs if the user passed any */
+ if (env != NULL)
+ {
+ parse_gpu_id_plain_string(env,
+ &gpu_opt->ncuda_dev_use,
+ &gpu_opt->cuda_dev_use);
+
+ if (gpu_opt->ncuda_dev_use == 0)
+ {
+ gmx_fatal(FARGS, "Empty GPU ID string encountered.\n%s\n",
+ invalid_gpuid_hint);
+ }
+
+ gpu_opt->bUserSet = TRUE;
+ }
+}
+
+void gmx_select_gpu_ids(FILE *fplog, const t_commrec *cr,
+ const gmx_gpu_info_t *gpu_info,
+ gmx_bool bForceUseGPU,
+ gmx_gpu_opt_t *gpu_opt)
+{
+ int i;
+ const char *env;
+ char sbuf[STRLEN], stmp[STRLEN];
+
+ /* Bail if binary is not compiled with GPU acceleration, but this is either
+ * explicitly (-nb gpu) or implicitly (gpu ID passed) requested. */
+ if (bForceUseGPU && !bGPUBinary)
+ {
+ gmx_fatal(FARGS, "GPU acceleration requested, but %s was compiled without GPU support!", ShortProgram());
+ }
+
+ if (gpu_opt->bUserSet)
+ {
+ /* Check the GPU IDs passed by the user.
+ * (GPU IDs have been parsed by gmx_parse_gpu_ids before)
+ */
+ int *checkres;
+ int res;
+
+ snew(checkres, gpu_opt->ncuda_dev_use);
+
+ res = check_selected_cuda_gpus(checkres, gpu_info, gpu_opt);
+
+ if (!res)
+ {
+ print_gpu_detection_stats(fplog, gpu_info, cr);
+
+ sprintf(sbuf, "Some of the requested GPUs do not exist, behave strangely, or are not compatible:\n");
+ for (i = 0; i < gpu_opt->ncuda_dev_use; i++)
+ {
+ if (checkres[i] != egpuCompatible)
+ {
+ sprintf(stmp, " GPU #%d: %s\n",
+ gpu_opt->cuda_dev_use[i],
+ gpu_detect_res_str[checkres[i]]);
+ strcat(sbuf, stmp);
+ }
+ }
+ gmx_fatal(FARGS, "%s", sbuf);
+ }
+
+ sfree(checkres);
+ }
+ else
+ {
+ pick_compatible_gpus(&hwinfo_g->gpu_info, gpu_opt);
+
+ if (gpu_opt->ncuda_dev_use > cr->nrank_pp_intranode)
+ {
+ /* We picked more GPUs than we can use: limit the number.
+ * We print detailed messages about this later in
+ * gmx_check_hw_runconf_consistency.
+ */
+ limit_num_gpus_used(gpu_opt, cr->nrank_pp_intranode);
+ }
+
+ gpu_opt->bUserSet = FALSE;
+ }
+
+ /* If the user asked for a GPU, check whether we have a GPU */
+ if (bForceUseGPU && gpu_info->ncuda_dev_compatible == 0)
+ {
+ gmx_fatal(FARGS, "GPU acceleration requested, but no compatible GPUs were detected.");
+ }
+}
+
+static void limit_num_gpus_used(gmx_gpu_opt_t *gpu_opt, int count)
+{
+ int ndev_use;
+
+ assert(gpu_opt);
+
+ ndev_use = gpu_opt->ncuda_dev_use;
+
+ if (count > ndev_use)
+ {
+ /* won't increase the # of GPUs */
+ return;
+ }
+
+ if (count < 1)
+ {
+ char sbuf[STRLEN];
+ sprintf(sbuf, "Limiting the number of GPUs to <1 doesn't make sense (detected %d, %d requested)!",
+ ndev_use, count);
+ gmx_incons(sbuf);
+ }
+
+ /* TODO: improve this implementation: either sort GPUs or remove the weakest here */
+ gpu_opt->ncuda_dev_use = count;
+}
+
+void gmx_hardware_info_free(gmx_hw_info_t *hwinfo)
+{
+ int ret;
+
+ ret = tMPI_Thread_mutex_lock(&hw_info_lock);
+ if (ret != 0)
+ {
+ gmx_fatal(FARGS, "Error locking hwinfo mutex: %s", strerror(errno));
+ }
+
+ /* decrease the reference counter */
+ n_hwinfo--;
+
+
+ if (hwinfo != hwinfo_g)
+ {
+ gmx_incons("hwinfo < hwinfo_g");
+ }
+
+ if (n_hwinfo < 0)
+ {
+ gmx_incons("n_hwinfo < 0");
+ }
+
+ if (n_hwinfo == 0)
+ {
+ gmx_cpuid_done(hwinfo_g->cpuid_info);
+ free_gpu_info(&hwinfo_g->gpu_info);
+ sfree(hwinfo_g);
+ }
+
+ ret = tMPI_Thread_mutex_unlock(&hw_info_lock);
+ if (ret != 0)
+ {
+ gmx_fatal(FARGS, "Error unlocking hwinfo mutex: %s", strerror(errno));
+ }
+}
--- /dev/null
- real *refproj; /* starting or target projecions */
+/* -*- 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:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <time.h>
+#include "typedefs.h"
+#include "string2.h"
+#include "smalloc.h"
+#include "names.h"
+#include "gromacs/fileio/confio.h"
+#include "mvdata.h"
+#include "txtdump.h"
+#include "vec.h"
+#include <time.h>
+#include "nrnb.h"
+#include "mshift.h"
+#include "mdrun.h"
+#include "update.h"
+#include "physics.h"
+#include "nrjac.h"
+#include "mtop_util.h"
+#include "edsam.h"
+#include "gromacs/fileio/gmxfio.h"
+#include "xvgr.h"
+#include "groupcoord.h"
+
+
+/* We use the same defines as in mvdata.c here */
+#define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
+#define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
+#define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
+
+/* These macros determine the column width in the output file */
+#define EDcol_sfmt "%17s"
+#define EDcol_efmt "%17.5e"
+#define EDcol_ffmt "%17f"
+
+/* enum to identify the type of ED: none, normal ED, flooding */
+enum {
+ eEDnone, eEDedsam, eEDflood, eEDnr
+};
+
+/* enum to identify operations on reference, average, origin, target structures */
+enum {
+ eedREF, eedAV, eedORI, eedTAR, eedNR
+};
+
+
+typedef struct
+{
+ int neig; /* nr of eigenvectors */
+ int *ieig; /* index nrs of eigenvectors */
+ real *stpsz; /* stepsizes (per eigenvector) */
+ rvec **vec; /* eigenvector components */
+ real *xproj; /* instantaneous x projections */
+ real *fproj; /* instantaneous f projections */
+ real radius; /* instantaneous radius */
- gmx_bool bNeedDoEdsam; /* if any of the options mon, linfix, ...
- * is used (i.e. apart from flooding) */
++ real *refproj; /* starting or target projections */
+ /* When using flooding as harmonic restraint: The current reference projection
+ * is at each step calculated from the initial refproj0 and the slope. */
+ real *refproj0, *refprojslope;
+} t_eigvec;
+
+
+typedef struct
+{
+ t_eigvec mon; /* only monitored, no constraints */
+ t_eigvec linfix; /* fixed linear constraints */
+ t_eigvec linacc; /* acceptance linear constraints */
+ t_eigvec radfix; /* fixed radial constraints (exp) */
+ t_eigvec radacc; /* acceptance radial constraints (exp) */
+ t_eigvec radcon; /* acceptance rad. contraction constr. */
+} t_edvecs;
+
+
+typedef struct
+{
+ real deltaF0;
+ gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
+ the eigenvector */
+ gmx_bool bConstForce; /* Do not calculate a flooding potential,
+ instead flood with a constant force */
+ real tau;
+ real deltaF;
+ real Efl;
+ real kT;
+ real Vfl;
+ real dt;
+ real constEfl;
+ real alpha2;
+ rvec *forces_cartesian;
+ t_eigvec vecs; /* use flooding for these */
+} t_edflood;
+
+
+/* This type is for the average, reference, target, and origin structure */
+typedef struct gmx_edx
+{
+ int nr; /* number of atoms this structure contains */
+ int nr_loc; /* number of atoms on local node */
+ int *anrs; /* atom index numbers */
+ int *anrs_loc; /* local atom index numbers */
+ int nalloc_loc; /* allocation size of anrs_loc */
+ int *c_ind; /* at which position of the whole anrs
+ * array is a local atom?, i.e.
+ * c_ind[0...nr_loc-1] gives the atom index
+ * with respect to the collective
+ * anrs[0...nr-1] array */
+ rvec *x; /* positions for this structure */
+ rvec *x_old; /* Last positions which have the correct PBC
+ representation of the ED group. In
+ combination with keeping track of the
+ shift vectors, the ED group can always
+ be made whole */
+ real *m; /* masses */
+ real mtot; /* total mass (only used in sref) */
+ real *sqrtm; /* sqrt of the masses used for mass-
+ * weighting of analysis (only used in sav) */
+} t_gmx_edx;
+
+
+typedef struct edpar
+{
+ 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 */
+
+ /* all gmx_edx datasets are copied to all nodes in the parallel case */
+ struct gmx_edx sref; /* reference positions, to these fitting
+ * will be done */
+ gmx_bool bRefEqAv; /* If true, reference & average indices
+ * are the same. Used for optimization */
+ struct gmx_edx sav; /* average positions */
+ struct gmx_edx star; /* target positions */
+ struct gmx_edx sori; /* origin positions */
+
+ t_edvecs vecs; /* eigenvectors */
+ real slope; /* minimal slope in acceptance radexp */
+
- * subroutine again, because these routines are rarely used simultanely */
+ t_edflood flood; /* parameters especially for flooding */
+ struct t_ed_buffer *buf; /* handle to local buffers */
+ struct edpar *next_edi; /* Pointer to another ED group */
+} t_edpar;
+
+
+typedef struct gmx_edsam
+{
+ int eEDtype; /* Type of ED: see enums above */
+ FILE *edo; /* output file pointer */
+ t_edpar *edpar;
+ gmx_bool bFirst;
+} t_gmx_edsam;
+
+
+struct t_do_edsam
+{
+ matrix old_rotmat;
+ real oldrad;
+ rvec old_transvec, older_transvec, transvec_compact;
+ rvec *xcoll; /* Positions from all nodes, this is the
+ collective set we work on.
+ These are the positions of atoms with
+ average structure indices */
+ rvec *xc_ref; /* same but with reference structure indices */
+ ivec *shifts_xcoll; /* Shifts for xcoll */
+ ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
+ ivec *shifts_xc_ref; /* Shifts for xc_ref */
+ ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
+ gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
+ ED shifts for this ED group need to
+ be updated */
+};
+
+
+/* definition of ED buffer structure */
+struct t_ed_buffer
+{
+ struct t_fit_to_ref * fit_to_ref;
+ struct t_do_edfit * do_edfit;
+ struct t_do_edsam * do_edsam;
+ struct t_do_radcon * do_radcon;
+};
+
+
+/* Function declarations */
+static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
+static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
+static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
+static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
+static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate);
+static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate);
+static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv);
+/* End function declarations */
+
+
++/* Do we have to perform essential dynamics constraints or possibly only flooding
++ * for any of the ED groups? */
++static gmx_bool bNeedDoEdsam(t_edpar *edi)
++{
++ return edi->vecs.mon.neig
++ || edi->vecs.linfix.neig
++ || edi->vecs.linacc.neig
++ || edi->vecs.radfix.neig
++ || edi->vecs.radacc.neig
++ || edi->vecs.radcon.neig;
++}
++
++
+/* Multiple ED groups will be labeled with letters instead of numbers
+ * to avoid confusion with eigenvector indices */
+static char get_EDgroupChar(int nr_edi, int nED)
+{
+ if (nED == 1)
+ {
+ return ' ';
+ }
+
+ /* nr_edi = 1 -> A
+ * nr_edi = 2 -> B ...
+ */
+ return 'A' + nr_edi - 1;
+}
+
+
+/* Does not subtract average positions, projection on single eigenvector is returned
+ * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
+ * Average position is subtracted in ed_apply_constraints prior to calling projectx
+ */
+static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
+{
+ int i;
+ real proj = 0.0;
+
+
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
+ }
+
+ return proj;
+}
+
+
+/* Specialized: projection is stored in vec->refproj
+ * -> used for radacc, radfix, radcon and center of flooding potential
+ * subtracts average positions, projects vector x */
+static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
+{
+ int i;
+ real rad = 0.0;
+
+ /* Subtract average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_dec(x[i], edi->sav.x[i]);
+ }
+
+ for (i = 0; i < vec->neig; i++)
+ {
+ vec->refproj[i] = projectx(edi, x, vec->vec[i]);
+ rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
+ }
+ vec->radius = sqrt(rad);
+
+ /* Add average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_inc(x[i], edi->sav.x[i]);
+ }
+}
+
+
+/* Project vector x, subtract average positions prior to projection and add
+ * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
+ * is applied. */
+static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
+ t_eigvec *vec, /* The eigenvectors */
+ t_edpar *edi)
+{
+ int i;
+
+
+ if (!vec->neig)
+ {
+ return;
+ }
+
+ /* Subtract average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_dec(x[i], edi->sav.x[i]);
+ }
+
+ for (i = 0; i < vec->neig; i++)
+ {
+ vec->xproj[i] = projectx(edi, x, vec->vec[i]);
+ }
+
+ /* Add average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_inc(x[i], edi->sav.x[i]);
+ }
+}
+
+
+/* Project vector x onto all edi->vecs (mon, linfix,...) */
+static void project(rvec *x, /* positions to project */
+ t_edpar *edi) /* edi data set */
+{
+ /* It is not more work to subtract the average position in every
- /* Remember for each ED group whether we have to do essential dynamics
- * constraints or possibly only flooding */
- edi->bNeedDoEdsam = edi->vecs.mon.neig
- || edi->vecs.linfix.neig
- || edi->vecs.linacc.neig
- || edi->vecs.radfix.neig
- || edi->vecs.radacc.neig
- || edi->vecs.radcon.neig;
-
++ * subroutine again, because these routines are rarely used simultaneously */
+ project_to_eigvectors(x, &edi->vecs.mon, edi);
+ project_to_eigvectors(x, &edi->vecs.linfix, edi);
+ project_to_eigvectors(x, &edi->vecs.linacc, edi);
+ project_to_eigvectors(x, &edi->vecs.radfix, edi);
+ project_to_eigvectors(x, &edi->vecs.radacc, edi);
+ project_to_eigvectors(x, &edi->vecs.radcon, edi);
+}
+
+
+static real calc_radius(t_eigvec *vec)
+{
+ int i;
+ real rad = 0.0;
+
+
+ for (i = 0; i < vec->neig; i++)
+ {
+ rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
+ }
+
+ return rad = sqrt(rad);
+}
+
+
+/* Debug helper */
+#ifdef DEBUGHELPERS
+static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
+ int step)
+{
+ int i;
+ FILE *fp;
+ char fn[STRLEN];
+ rvec *xcoll;
+ ivec *shifts, *eshifts;
+
+
+ if (!MASTER(cr))
+ {
+ return;
+ }
+
+ xcoll = buf->xcoll;
+ shifts = buf->shifts_xcoll;
+ eshifts = buf->extra_shifts_xcoll;
+
+ sprintf(fn, "xcolldump_step%d.txt", step);
+ fp = fopen(fn, "w");
+
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
+ edi->sav.anrs[i]+1,
+ xcoll[i][XX], xcoll[i][YY], xcoll[i][ZZ],
+ shifts[i][XX], shifts[i][YY], shifts[i][ZZ],
+ eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
+ }
+
+ fclose(fp);
+}
+
+
+/* Debug helper */
+static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
+{
+ int i;
+
+
+ fprintf(out, "#%s positions:\n%d\n", name, s->nr);
+ if (s->nr == 0)
+ {
+ return;
+ }
+
+ fprintf(out, "#index, x, y, z");
+ if (s->sqrtm)
+ {
+ fprintf(out, ", sqrt(m)");
+ }
+ for (i = 0; i < s->nr; i++)
+ {
+ fprintf(out, "\n%6d %11.6f %11.6f %11.6f", s->anrs[i], s->x[i][XX], s->x[i][YY], s->x[i][ZZ]);
+ if (s->sqrtm)
+ {
+ fprintf(out, "%9.3f", s->sqrtm[i]);
+ }
+ }
+ fprintf(out, "\n");
+}
+
+
+/* Debug helper */
+static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
+ const char name[], int length)
+{
+ int i, j;
+
+
+ fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
+ /* Dump the data for every eigenvector: */
+ for (i = 0; i < ev->neig; i++)
+ {
+ fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
+ ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
+ for (j = 0; j < length; j++)
+ {
+ fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
+ }
+ }
+}
+
+
+/* Debug helper */
+static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
+{
+ FILE *out;
+ char fn[STRLEN];
+
+
+ sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi);
+ out = ffopen(fn, "w");
+
+ fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
+ edpars->nini, edpars->fitmas, edpars->pcamas);
+ fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
+ edpars->outfrq, edpars->maxedsteps, edpars->slope);
+ fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
+ edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
+ edpars->flood.constEfl, edpars->flood.alpha2);
+
+ /* Dump reference, average, target, origin positions */
+ dump_edi_positions(out, &edpars->sref, "REFERENCE");
+ dump_edi_positions(out, &edpars->sav, "AVERAGE" );
+ dump_edi_positions(out, &edpars->star, "TARGET" );
+ dump_edi_positions(out, &edpars->sori, "ORIGIN" );
+
+ /* Dump eigenvectors */
+ dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
+ dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
+ dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
+ dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
+ dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
+ dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
+
+ /* Dump flooding eigenvectors */
+ dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
+
+ /* Dump ed local buffer */
+ fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
+ fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
+ fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
+
+ ffclose(out);
+}
+
+
+/* Debug helper */
+static void dump_rotmat(FILE* out, matrix rotmat)
+{
+ fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[XX][XX], rotmat[XX][YY], rotmat[XX][ZZ]);
+ fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[YY][XX], rotmat[YY][YY], rotmat[YY][ZZ]);
+ fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[ZZ][XX], rotmat[ZZ][YY], rotmat[ZZ][ZZ]);
+}
+
+
+/* Debug helper */
+static void dump_rvec(FILE *out, int dim, rvec *x)
+{
+ int i;
+
+
+ for (i = 0; i < dim; i++)
+ {
+ fprintf(out, "%4d %f %f %f\n", i, x[i][XX], x[i][YY], x[i][ZZ]);
+ }
+}
+
+
+/* Debug helper */
+static void dump_mat(FILE* out, int dim, double** mat)
+{
+ int i, j;
+
+
+ fprintf(out, "MATRIX:\n");
+ for (i = 0; i < dim; i++)
+ {
+ for (j = 0; j < dim; j++)
+ {
+ fprintf(out, "%f ", mat[i][j]);
+ }
+ fprintf(out, "\n");
+ }
+}
+#endif
+
+
+struct t_do_edfit {
+ double **omega;
+ double **om;
+};
+
+static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
+{
+ /* this is a copy of do_fit with some modifications */
+ int c, r, n, j, i, irot;
+ double d[6], xnr, xpc;
+ matrix vh, vk, u;
+ int index;
+ real max_d;
+
+ struct t_do_edfit *loc;
+ gmx_bool bFirst;
+
+ if (edi->buf->do_edfit != NULL)
+ {
+ bFirst = FALSE;
+ }
+ else
+ {
+ bFirst = TRUE;
+ snew(edi->buf->do_edfit, 1);
+ }
+ loc = edi->buf->do_edfit;
+
+ if (bFirst)
+ {
+ snew(loc->omega, 2*DIM);
+ snew(loc->om, 2*DIM);
+ for (i = 0; i < 2*DIM; i++)
+ {
+ snew(loc->omega[i], 2*DIM);
+ snew(loc->om[i], 2*DIM);
+ }
+ }
+
+ for (i = 0; (i < 6); i++)
+ {
+ d[i] = 0;
+ for (j = 0; (j < 6); j++)
+ {
+ loc->omega[i][j] = 0;
+ loc->om[i][j] = 0;
+ }
+ }
+
+ /* calculate the matrix U */
+ clear_mat(u);
+ for (n = 0; (n < natoms); n++)
+ {
+ for (c = 0; (c < DIM); c++)
+ {
+ xpc = xp[n][c];
+ for (r = 0; (r < DIM); r++)
+ {
+ xnr = x[n][r];
+ u[c][r] += xnr*xpc;
+ }
+ }
+ }
+
+ /* construct loc->omega */
+ /* loc->omega is symmetric -> loc->omega==loc->omega' */
+ for (r = 0; (r < 6); r++)
+ {
+ for (c = 0; (c <= r); c++)
+ {
+ if ((r >= 3) && (c < 3))
+ {
+ loc->omega[r][c] = u[r-3][c];
+ loc->omega[c][r] = u[r-3][c];
+ }
+ else
+ {
+ loc->omega[r][c] = 0;
+ loc->omega[c][r] = 0;
+ }
+ }
+ }
+
+ /* determine h and k */
+#ifdef DEBUG
+ {
+ int i;
+ dump_mat(stderr, 2*DIM, loc->omega);
+ for (i = 0; i < 6; i++)
+ {
+ fprintf(stderr, "d[%d] = %f\n", i, d[i]);
+ }
+ }
+#endif
+ jacobi(loc->omega, 6, d, loc->om, &irot);
+
+ if (irot == 0)
+ {
+ fprintf(stderr, "IROT=0\n");
+ }
+
+ index = 0; /* For the compiler only */
+
+ for (j = 0; (j < 3); j++)
+ {
+ max_d = -1000;
+ for (i = 0; (i < 6); i++)
+ {
+ if (d[i] > max_d)
+ {
+ max_d = d[i];
+ index = i;
+ }
+ }
+ d[index] = -10000;
+ for (i = 0; (i < 3); i++)
+ {
+ vh[j][i] = M_SQRT2*loc->om[i][index];
+ vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
+ }
+ }
+
+ /* determine R */
+ for (c = 0; (c < 3); c++)
+ {
+ for (r = 0; (r < 3); r++)
+ {
+ R[c][r] = vk[0][r]*vh[0][c]+
+ vk[1][r]*vh[1][c]+
+ vk[2][r]*vh[2][c];
+ }
+ }
+ if (det(R) < 0)
+ {
+ for (c = 0; (c < 3); c++)
+ {
+ for (r = 0; (r < 3); r++)
+ {
+ R[c][r] = vk[0][r]*vh[0][c]+
+ vk[1][r]*vh[1][c]-
+ vk[2][r]*vh[2][c];
+ }
+ }
+ }
+}
+
+
+static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
+{
+ rvec vec;
+ matrix tmat;
+
+
+ /* Remove rotation.
+ * The inverse rotation is described by the transposed rotation matrix */
+ transpose(rotmat, tmat);
+ rotate_x(xcoll, nat, tmat);
+
+ /* Remove translation */
+ vec[XX] = -transvec[XX];
+ vec[YY] = -transvec[YY];
+ vec[ZZ] = -transvec[ZZ];
+ translate_x(xcoll, nat, vec);
+}
+
+
+/**********************************************************************************
+ ******************** FLOODING ****************************************************
+ **********************************************************************************
+
+ The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
+ The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
+ read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
+
+ do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
+ the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
+
+ since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
+ edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
+
+ do_flood makes a copy of the positions,
+ fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
+ space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
+ fit. Then do_flood adds these forces to the forcefield-forces
+ (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
+
+ To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
+ structure is projected to the system of eigenvectors and then this position in the subspace is used as
+ center of the flooding potential. If the option is not used, the center will be zero in the subspace,
+ i.e. the average structure as given in the make_edi file.
+
+ To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
+ signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
+ Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
+ so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
+ further adaption is applied, Efl will stay constant at zero.
+
+ To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
+ used as spring constants for the harmonic potential.
+ Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
+ as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
+
+ To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
+ the routine read_edi_file reads all of theses flooding files.
+ The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
+ calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
+ edi there is no interdependence whatsoever. The forces are added together.
+
+ To write energies into the .edr file, call the function
+ get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
+ and call
+ get_flood_energies(real Vfl[],int nnames);
+
+ TODO:
+ - one could program the whole thing such that Efl, Vfl and deltaF is written to the .edr file. -- i dont know how to do that, yet.
+
+ Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
+ two edsam files from two peptide chains
+ */
+
+static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
+{
+ int i;
+
+
+ /* Output how well we fit to the reference structure */
+ fprintf(fp, EDcol_ffmt, rmsd);
+
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
+
+ /* Check whether the reference projection changes with time (this can happen
+ * in case flooding is used as harmonic restraint). If so, output the
+ * current reference projection */
+ if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
+ {
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
+ }
+
+ /* Output Efl if we are doing adaptive flooding */
+ if (0 != edi->flood.tau)
+ {
+ fprintf(fp, EDcol_efmt, edi->flood.Efl);
+ }
+ fprintf(fp, EDcol_efmt, edi->flood.Vfl);
+
+ /* Output deltaF if we are doing adaptive flooding */
+ if (0 != edi->flood.tau)
+ {
+ fprintf(fp, EDcol_efmt, edi->flood.deltaF);
+ }
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
+ }
+}
+
+
+/* From flood.xproj compute the Vfl(x) at this point */
+static real flood_energy(t_edpar *edi, gmx_int64_t step)
+{
+ /* compute flooding energy Vfl
+ Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
+ \lambda_i is the reciprocal eigenvalue 1/\sigma_i
+ it is already computed by make_edi and stored in stpsz[i]
+ bHarmonic:
+ Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
+ */
+ real sum;
+ real Vfl;
+ int i;
+
+
+ /* Each time this routine is called (i.e. each time step), we add a small
+ * value to the reference projection. This way a harmonic restraint towards
+ * a moving reference is realized. If no value for the additive constant
+ * is provided in the edi file, the reference will not change. */
+ if (edi->flood.bHarmonic)
+ {
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
+ }
+ }
+
+ sum = 0.0;
+ /* Compute sum which will be the exponent of the exponential */
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
+ sum += edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i])*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
+ }
+
+ /* Compute the Gauss function*/
+ if (edi->flood.bHarmonic)
+ {
+ Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
+ }
+ else
+ {
+ Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
+ }
+
+ return Vfl;
+}
+
+
+/* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
+static void flood_forces(t_edpar *edi)
+{
+ /* compute the forces in the subspace of the flooding eigenvectors
+ * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
+
+ int i;
+ real energy = edi->flood.Vfl;
+
+
+ if (edi->flood.bHarmonic)
+ {
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
+ }
+ }
+ else
+ {
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ /* if Efl is zero the forces are zero if not use the formula */
+ edi->flood.vecs.fproj[i] = edi->flood.Efl != 0 ? edi->flood.kT/edi->flood.Efl/edi->flood.alpha2*energy*edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]) : 0;
+ }
+ }
+}
+
+
+/* Raise forces from subspace into cartesian space */
+static void flood_blowup(t_edpar *edi, rvec *forces_cart)
+{
+ /* this function lifts the forces from the subspace to the cartesian space
+ all the values not contained in the subspace are assumed to be zero and then
+ a coordinate transformation from eigenvector to cartesian vectors is performed
+ The nonexistent values don't have to be set to zero explicitly, they would occur
+ as zero valued summands, hence we just stop to compute this part of the sum.
+
+ for every atom we add all the contributions to this atom from all the different eigenvectors.
+
+ NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
+ field forces_cart prior the computation, but we compute the forces separately
+ to have them accessible for diagnostics
+ */
+ int j, eig;
+ rvec dum;
+ real *forces_sub;
+
+
+ forces_sub = edi->flood.vecs.fproj;
+
+
+ /* Calculate the cartesian forces for the local atoms */
+
+ /* Clear forces first */
+ for (j = 0; j < edi->sav.nr_loc; j++)
+ {
+ clear_rvec(forces_cart[j]);
+ }
+
+ /* Now compute atomwise */
+ for (j = 0; j < edi->sav.nr_loc; j++)
+ {
+ /* Compute forces_cart[edi->sav.anrs[j]] */
+ for (eig = 0; eig < edi->flood.vecs.neig; eig++)
+ {
+ /* Force vector is force * eigenvector (compute only atom j) */
+ svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
+ /* Add this vector to the cartesian forces */
+ rvec_inc(forces_cart[j], dum);
+ }
+ }
+}
+
+
+/* Update the values of Efl, deltaF depending on tau and Vfl */
+static void update_adaption(t_edpar *edi)
+{
+ /* this function updates the parameter Efl and deltaF according to the rules given in
+ * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
+ * J. chem Phys. */
+
+ if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
+ {
+ edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
+ /* check if restrain (inverted flooding) -> don't let EFL become positive */
+ if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
+ {
+ edi->flood.Efl = 0;
+ }
+
+ edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
+ }
+}
+
+
+static void do_single_flood(
+ FILE *edo,
+ rvec x[],
+ rvec force[],
+ t_edpar *edi,
+ gmx_int64_t step,
+ matrix box,
+ t_commrec *cr,
+ gmx_bool bNS) /* Are we in a neighbor searching step? */
+{
+ int i;
+ matrix rotmat; /* rotation matrix */
+ matrix tmat; /* inverse rotation */
+ rvec transvec; /* translation vector */
+ real rmsdev;
+ struct t_do_edsam *buf;
+
+
+ buf = edi->buf->do_edsam;
+
+
+ /* Broadcast the positions of the AVERAGE structure such that they are known on
+ * every processor. Each node contributes its local positions x and stores them in
+ * the collective ED array buf->xcoll */
+ communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
+ edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
+
+ /* Only assembly REFERENCE positions if their indices differ from the average ones */
+ if (!edi->bRefEqAv)
+ {
+ communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
+ edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
+ }
+
+ /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
+ * We do not need to update the shifts until the next NS step */
+ buf->bUpdateShifts = FALSE;
+
+ /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
+ * as well as the indices in edi->sav.anrs */
+
+ /* Fit the reference indices to the reference structure */
+ if (edi->bRefEqAv)
+ {
+ fit_to_reference(buf->xcoll, transvec, rotmat, edi);
+ }
+ else
+ {
+ fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
+ }
+
+ /* Now apply the translation and rotation to the ED structure */
+ translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
+
+ /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
+ project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
+
+ if (FALSE == edi->flood.bConstForce)
+ {
+ /* Compute Vfl(x) from flood.xproj */
+ edi->flood.Vfl = flood_energy(edi, step);
+
+ update_adaption(edi);
+
+ /* Compute the flooding forces */
+ flood_forces(edi);
+ }
+
+ /* Translate them into cartesian positions */
+ flood_blowup(edi, edi->flood.forces_cartesian);
+
+ /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
+ /* Each node rotates back its local forces */
+ transpose(rotmat, tmat);
+ rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
+
+ /* Finally add forces to the main force variable */
+ for (i = 0; i < edi->sav.nr_loc; i++)
+ {
+ rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
+ }
+
+ /* Output is written by the master process */
+ if (do_per_step(step, edi->outfrq) && MASTER(cr))
+ {
+ /* Output how well we fit to the reference */
+ if (edi->bRefEqAv)
+ {
+ /* Indices of reference and average structures are identical,
+ * thus we can calculate the rmsd to SREF using xcoll */
+ rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
+ }
+ else
+ {
+ /* We have to translate & rotate the reference atoms first */
+ translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
+ rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
+ }
+
+ write_edo_flood(edi, edo, rmsdev);
+ }
+}
+
+
+/* Main flooding routine, called from do_force */
+extern void do_flood(
+ t_commrec *cr, /* Communication record */
+ t_inputrec *ir, /* Input record */
+ rvec x[], /* Positions on the local processor */
+ rvec force[], /* forcefield forces, to these the flooding forces are added */
+ gmx_edsam_t ed, /* ed data structure contains all ED and flooding groups */
+ matrix box, /* the box */
+ gmx_int64_t step, /* The relative time step since ir->init_step is already subtracted */
+ gmx_bool bNS) /* Are we in a neighbor searching step? */
+{
+ t_edpar *edi;
+
+
+ edi = ed->edpar;
+
+ /* Write time to edo, when required. Output the time anyhow since we need
+ * it in the output file for ED constraints. */
+ if (MASTER(cr) && do_per_step(step, edi->outfrq))
+ {
+ fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
+ }
+
+ if (ed->eEDtype != eEDflood)
+ {
+ return;
+ }
+
+ while (edi)
+ {
+ /* Call flooding for one matrix */
+ if (edi->flood.vecs.neig)
+ {
+ do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
+ }
+ edi = edi->next_edi;
+ }
+}
+
+
+/* Called by init_edi, configure some flooding related variables and structures,
+ * print headers to output files */
+static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
+{
+ int i;
+
+
+ edi->flood.Efl = edi->flood.constEfl;
+ edi->flood.Vfl = 0;
+ edi->flood.dt = dt;
+
+ if (edi->flood.vecs.neig)
+ {
+ /* If in any of the ED groups we find a flooding vector, flooding is turned on */
+ ed->eEDtype = eEDflood;
+
+ fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
+
+ if (edi->flood.bConstForce)
+ {
+ /* We have used stpsz as a vehicle to carry the fproj values for constant
+ * force flooding. Now we copy that to flood.vecs.fproj. Note that
+ * in const force flooding, fproj is never changed. */
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
+
+ fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
+ edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
+ }
+ }
+ }
+}
+
+
+#ifdef DEBUGHELPERS
+/*********** Energy book keeping ******/
+static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
+{
+ t_edpar *actual;
+ int count;
+ char buf[STRLEN];
+ actual = edi;
+ count = 1;
+ while (actual)
+ {
+ srenew(names, count);
+ sprintf(buf, "Vfl_%d", count);
+ names[count-1] = strdup(buf);
+ actual = actual->next_edi;
+ count++;
+ }
+ *nnames = count-1;
+}
+
+
+static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
+{
+ /*fl has to be big enough to capture nnames-many entries*/
+ t_edpar *actual;
+ int count;
+
+
+ actual = edi;
+ count = 1;
+ while (actual)
+ {
+ Vfl[count-1] = actual->flood.Vfl;
+ actual = actual->next_edi;
+ count++;
+ }
+ if (nnames != count-1)
+ {
+ gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
+ }
+}
+/************* END of FLOODING IMPLEMENTATION ****************************/
+#endif
+
+
+gmx_edsam_t ed_open(int natoms, edsamstate_t *EDstate, int nfile, const t_filenm fnm[], unsigned long Flags, const output_env_t oenv, t_commrec *cr)
+{
+ gmx_edsam_t ed;
+ int nED;
+
+
+ /* Allocate space for the ED data structure */
+ snew(ed, 1);
+
+ /* We want to perform ED (this switch might later be upgraded to eEDflood) */
+ ed->eEDtype = eEDedsam;
+
+ if (MASTER(cr))
+ {
+ fprintf(stderr, "ED sampling will be performed!\n");
+ snew(ed->edpar, 1);
+
+ /* Read the edi input file: */
+ nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms);
+
+ /* Make sure the checkpoint was produced in a run using this .edi file */
+ if (EDstate->bFromCpt)
+ {
+ crosscheck_edi_file_vs_checkpoint(ed, EDstate);
+ }
+ else
+ {
+ EDstate->nED = nED;
+ }
+ init_edsamstate(ed, EDstate);
+
+ /* The master opens the ED output file */
+ if (Flags & MD_APPENDFILES)
+ {
+ ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+");
+ }
+ else
+ {
+ ed->edo = xvgropen(opt2fn("-eo", nfile, fnm),
+ "Essential dynamics / flooding output",
+ "Time (ps)",
+ "RMSDs (nm), projections on EVs (nm), ...", oenv);
+
+ /* Make a descriptive legend */
+ write_edo_legend(ed, EDstate->nED, oenv);
+ }
+ }
+ return ed;
+}
+
+
+/* Broadcasts the structure data */
+static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
+{
+ snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
+ snew_bc(cr, s->x, s->nr ); /* Positions */
+ nblock_bc(cr, s->nr, s->anrs );
+ nblock_bc(cr, s->nr, s->x );
+
+ /* For the average & reference structures we need an array for the collective indices,
+ * and we need to broadcast the masses as well */
+ if (stype == eedAV || stype == eedREF)
+ {
+ /* We need these additional variables in the parallel case: */
+ snew(s->c_ind, s->nr ); /* Collective indices */
+ /* Local atom indices get assigned in dd_make_local_group_indices.
+ * There, also memory is allocated */
+ s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
+ snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
+ nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
+ }
+
+ /* broadcast masses for the reference structure (for mass-weighted fitting) */
+ if (stype == eedREF)
+ {
+ snew_bc(cr, s->m, s->nr);
+ nblock_bc(cr, s->nr, s->m);
+ }
+
+ /* For the average structure we might need the masses for mass-weighting */
+ if (stype == eedAV)
+ {
+ snew_bc(cr, s->sqrtm, s->nr);
+ nblock_bc(cr, s->nr, s->sqrtm);
+ snew_bc(cr, s->m, s->nr);
+ nblock_bc(cr, s->nr, s->m);
+ }
+}
+
+
+/* Broadcasts the eigenvector data */
+static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
+{
+ int i;
+
+ snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
+ snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
+ snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
+ snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
+ snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
+
+ nblock_bc(cr, ev->neig, ev->ieig );
+ nblock_bc(cr, ev->neig, ev->stpsz );
+ nblock_bc(cr, ev->neig, ev->xproj );
+ nblock_bc(cr, ev->neig, ev->fproj );
+ nblock_bc(cr, ev->neig, ev->refproj);
+
+ snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
+ for (i = 0; i < ev->neig; i++)
+ {
+ snew_bc(cr, ev->vec[i], length);
+ nblock_bc(cr, length, ev->vec[i]);
+ }
+
+ /* For harmonic restraints the reference projections can change with time */
+ if (bHarmonic)
+ {
+ snew_bc(cr, ev->refproj0, ev->neig);
+ snew_bc(cr, ev->refprojslope, ev->neig);
+ nblock_bc(cr, ev->neig, ev->refproj0 );
+ nblock_bc(cr, ev->neig, ev->refprojslope);
+ }
+}
+
+
+/* Broadcasts the ED / flooding data to other nodes
+ * and allocates memory where needed */
+static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
+{
+ int nr;
+ t_edpar *edi;
+
+
+ /* Master lets the other nodes know if its ED only or also flooding */
+ gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
+
+ snew_bc(cr, ed->edpar, 1);
+ /* Now transfer the ED data set(s) */
+ edi = ed->edpar;
+ for (nr = 0; nr < numedis; nr++)
+ {
+ /* Broadcast a single ED data set */
+ block_bc(cr, *edi);
+
+ /* Broadcast positions */
+ bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
+ bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
+ bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
+ bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
+
+ /* Broadcast eigenvectors */
+ bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
+ /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
+ bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
+
+ /* Set the pointer to the next ED group */
+ if (edi->next_edi)
+ {
+ snew_bc(cr, edi->next_edi, 1);
+ edi = edi->next_edi;
+ }
+ }
+}
+
+
+/* init-routine called for every *.edi-cycle, initialises t_edpar structure */
+static void init_edi(gmx_mtop_t *mtop, t_edpar *edi)
+{
+ int i;
+ real totalmass = 0.0;
+ rvec com;
+ gmx_mtop_atomlookup_t alook = NULL;
+ t_atom *atom;
+
+ /* NOTE Init_edi is executed on the master process only
+ * The initialized data sets are then transmitted to the
+ * other nodes in broadcast_ed_data */
+
+ alook = gmx_mtop_atomlookup_init(mtop);
+
+ /* evaluate masses (reference structure) */
+ snew(edi->sref.m, edi->sref.nr);
+ for (i = 0; i < edi->sref.nr; i++)
+ {
+ if (edi->fitmas)
+ {
+ gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom);
+ edi->sref.m[i] = atom->m;
+ }
+ else
+ {
+ edi->sref.m[i] = 1.0;
+ }
+
+ /* Check that every m > 0. Bad things will happen otherwise. */
+ if (edi->sref.m[i] <= 0.0)
+ {
+ gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
+ "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
+ "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
+ "atoms from the reference structure by creating a proper index group.\n",
+ i, edi->sref.anrs[i]+1, edi->sref.m[i]);
+ }
+
+ totalmass += edi->sref.m[i];
+ }
+ edi->sref.mtot = totalmass;
+
+ /* Masses m and sqrt(m) for the average structure. Note that m
+ * is needed if forces have to be evaluated in do_edsam */
+ snew(edi->sav.sqrtm, edi->sav.nr );
+ snew(edi->sav.m, edi->sav.nr );
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom);
+ edi->sav.m[i] = atom->m;
+ if (edi->pcamas)
+ {
+ edi->sav.sqrtm[i] = sqrt(atom->m);
+ }
+ else
+ {
+ edi->sav.sqrtm[i] = 1.0;
+ }
+
+ /* Check that every m > 0. Bad things will happen otherwise. */
+ if (edi->sav.sqrtm[i] <= 0.0)
+ {
+ gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
+ "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
+ "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
+ "atoms from the average structure by creating a proper index group.\n",
+ i, edi->sav.anrs[i]+1, atom->m);
+ }
+ }
+
+ gmx_mtop_atomlookup_destroy(alook);
+
+ /* put reference structure in origin */
+ get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
+ com[XX] = -com[XX];
+ com[YY] = -com[YY];
+ com[ZZ] = -com[ZZ];
+ translate_x(edi->sref.x, edi->sref.nr, com);
+
+ /* Init ED buffer */
+ snew(edi->buf, 1);
+}
+
+
+static void check(const char *line, const char *label)
+{
+ if (!strstr(line, label))
+ {
+ gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
+ }
+}
+
+
+static int read_checked_edint(FILE *file, const char *label)
+{
+ char line[STRLEN+1];
+ int idum;
+
+
+ fgets2 (line, STRLEN, file);
+ check(line, label);
+ fgets2 (line, STRLEN, file);
+ sscanf (line, "%d", &idum);
+ return idum;
+}
+
+
+static int read_edint(FILE *file, gmx_bool *bEOF)
+{
+ char line[STRLEN+1];
+ int idum;
+ char *eof;
+
+
+ eof = fgets2 (line, STRLEN, file);
+ if (eof == NULL)
+ {
+ *bEOF = TRUE;
+ return -1;
+ }
+ eof = fgets2 (line, STRLEN, file);
+ if (eof == NULL)
+ {
+ *bEOF = TRUE;
+ return -1;
+ }
+ sscanf (line, "%d", &idum);
+ *bEOF = FALSE;
+ return idum;
+}
+
+
+static real read_checked_edreal(FILE *file, const char *label)
+{
+ char line[STRLEN+1];
+ double rdum;
+
+
+ fgets2 (line, STRLEN, file);
+ check(line, label);
+ fgets2 (line, STRLEN, file);
+ sscanf (line, "%lf", &rdum);
+ return (real) rdum; /* always read as double and convert to single */
+}
+
+
+static void read_edx(FILE *file, int number, int *anrs, rvec *x)
+{
+ int i, j;
+ char line[STRLEN+1];
+ double d[3];
+
+
+ for (i = 0; i < number; i++)
+ {
+ fgets2 (line, STRLEN, file);
+ sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]);
+ anrs[i]--; /* we are reading FORTRAN indices */
+ for (j = 0; j < 3; j++)
+ {
+ x[i][j] = d[j]; /* always read as double and convert to single */
+ }
+ }
+}
+
+
+static void scan_edvec(FILE *in, int nr, rvec *vec)
+{
+ char line[STRLEN+1];
+ int i;
+ double x, y, z;
+
+
+ for (i = 0; (i < nr); i++)
+ {
+ fgets2 (line, STRLEN, in);
+ sscanf (line, "%le%le%le", &x, &y, &z);
+ vec[i][XX] = x;
+ vec[i][YY] = y;
+ vec[i][ZZ] = z;
+ }
+}
+
+
+static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
+{
+ int i, idum, nscan;
+ double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
+ char line[STRLEN+1];
+
+
+ tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
+ if (tvec->neig > 0)
+ {
+ snew(tvec->ieig, tvec->neig);
+ snew(tvec->stpsz, tvec->neig);
+ snew(tvec->vec, tvec->neig);
+ snew(tvec->xproj, tvec->neig);
+ snew(tvec->fproj, tvec->neig);
+ snew(tvec->refproj, tvec->neig);
+ if (bReadRefproj)
+ {
+ snew(tvec->refproj0, tvec->neig);
+ snew(tvec->refprojslope, tvec->neig);
+ }
+
+ for (i = 0; (i < tvec->neig); i++)
+ {
+ fgets2 (line, STRLEN, in);
+ if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
+ {
+ nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum);
+ /* Zero out values which were not scanned */
+ switch (nscan)
+ {
+ case 4:
+ /* Every 4 values read, including reference position */
+ *bHaveReference = TRUE;
+ break;
+ case 3:
+ /* A reference position is provided */
+ *bHaveReference = TRUE;
+ /* No value for slope, set to 0 */
+ refprojslope_dum = 0.0;
+ break;
+ case 2:
+ /* No values for reference projection and slope, set to 0 */
+ refproj_dum = 0.0;
+ refprojslope_dum = 0.0;
+ break;
+ default:
+ gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
+ break;
+ }
+ tvec->refproj[i] = refproj_dum;
+ tvec->refproj0[i] = refproj_dum;
+ tvec->refprojslope[i] = refprojslope_dum;
+ }
+ else /* Normal flooding */
+ {
+ nscan = sscanf(line, "%d%lf", &idum, &rdum);
+ if (nscan != 2)
+ {
+ gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
+ }
+ }
+ tvec->ieig[i] = idum;
+ tvec->stpsz[i] = rdum;
+ } /* end of loop over eigenvectors */
+
+ for (i = 0; (i < tvec->neig); i++)
+ {
+ snew(tvec->vec[i], nr);
+ scan_edvec(in, nr, tvec->vec[i]);
+ }
+ }
+}
+
+
+/* calls read_edvec for the vector groups, only for flooding there is an extra call */
+static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
+{
+ gmx_bool bHaveReference = FALSE;
+
+
+ read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
+ read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
+ read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
+ read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
+ read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
+ read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
+}
+
+
+/* Check if the same atom indices are used for reference and average positions */
+static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
+{
+ int i;
+
+
+ /* If the number of atoms differs between the two structures,
+ * they cannot be identical */
+ if (sref.nr != sav.nr)
+ {
+ return FALSE;
+ }
+
+ /* Now that we know that both stuctures have the same number of atoms,
+ * check if also the indices are identical */
+ for (i = 0; i < sav.nr; i++)
+ {
+ if (sref.anrs[i] != sav.anrs[i])
+ {
+ return FALSE;
+ }
+ }
+ fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
+
+ return TRUE;
+}
+
+
+static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
+{
+ int readmagic;
+ const int magic = 670;
+ gmx_bool bEOF;
+
+ /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
+ gmx_bool bHaveReference = FALSE;
+
+
+ /* the edi file is not free format, so expect problems if the input is corrupt. */
+
+ /* check the magic number */
+ readmagic = read_edint(in, &bEOF);
+ /* Check whether we have reached the end of the input file */
+ if (bEOF)
+ {
+ return 0;
+ }
+
+ if (readmagic != magic)
+ {
+ if (readmagic == 666 || readmagic == 667 || readmagic == 668)
+ {
+ gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
+ }
+ else if (readmagic != 669)
+ {
+ gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
+ }
+ }
+
+ /* check the number of atoms */
+ edi->nini = read_edint(in, &bEOF);
+ if (edi->nini != nr_mdatoms)
+ {
+ gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
+ }
+
+ /* Done checking. For the rest we blindly trust the input */
+ edi->fitmas = read_checked_edint(in, "FITMAS");
+ edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
+ edi->outfrq = read_checked_edint(in, "OUTFRQ");
+ edi->maxedsteps = read_checked_edint(in, "MAXLEN");
+ edi->slope = read_checked_edreal(in, "SLOPECRIT");
+
+ edi->presteps = read_checked_edint(in, "PRESTEPS");
+ edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
+ edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
+ edi->flood.tau = read_checked_edreal(in, "TAU");
+ edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
+ edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
+ edi->flood.kT = read_checked_edreal(in, "KT");
+ edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
+ if (readmagic > 669)
+ {
+ edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
+ }
+ else
+ {
+ edi->flood.bConstForce = FALSE;
+ }
+ edi->sref.nr = read_checked_edint(in, "NREF");
+
+ /* allocate space for reference positions and read them */
+ snew(edi->sref.anrs, edi->sref.nr);
+ snew(edi->sref.x, edi->sref.nr);
+ snew(edi->sref.x_old, edi->sref.nr);
+ edi->sref.sqrtm = NULL;
+ read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
+
+ /* average positions. they define which atoms will be used for ED sampling */
+ edi->sav.nr = read_checked_edint(in, "NAV");
+ snew(edi->sav.anrs, edi->sav.nr);
+ snew(edi->sav.x, edi->sav.nr);
+ snew(edi->sav.x_old, edi->sav.nr);
+ read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
+
+ /* Check if the same atom indices are used for reference and average positions */
+ edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
+
+ /* eigenvectors */
+ read_edvecs(in, edi->sav.nr, &edi->vecs);
+ read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
+
+ /* target positions */
+ edi->star.nr = read_edint(in, &bEOF);
+ if (edi->star.nr > 0)
+ {
+ snew(edi->star.anrs, edi->star.nr);
+ snew(edi->star.x, edi->star.nr);
+ edi->star.sqrtm = NULL;
+ read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
+ }
+
+ /* positions defining origin of expansion circle */
+ edi->sori.nr = read_edint(in, &bEOF);
+ if (edi->sori.nr > 0)
+ {
+ if (bHaveReference)
+ {
+ /* Both an -ori structure and a at least one manual reference point have been
+ * specified. That's ambiguous and probably not intentional. */
+ gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
+ " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
+ }
+ snew(edi->sori.anrs, edi->sori.nr);
+ snew(edi->sori.x, edi->sori.nr);
+ edi->sori.sqrtm = NULL;
+ read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
+ }
+
+ /* all done */
+ return 1;
+}
+
+
+
+/* Read in the edi input file. Note that it may contain several ED data sets which were
+ * achieved by concatenating multiple edi files. The standard case would be a single ED
+ * data set, though. */
+static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
+{
+ FILE *in;
+ t_edpar *curr_edi, *last_edi;
+ t_edpar *edi_read;
+ int edi_nr = 0;
+
+
+ /* This routine is executed on the master only */
+
+ /* Open the .edi parameter input file */
+ in = gmx_fio_fopen(fn, "r");
+ fprintf(stderr, "ED: Reading edi file %s\n", fn);
+
+ /* Now read a sequence of ED input parameter sets from the edi file */
+ curr_edi = edi;
+ last_edi = edi;
+ while (read_edi(in, curr_edi, nr_mdatoms, fn) )
+ {
+ edi_nr++;
+
+ /* Since we arrived within this while loop we know that there is still another data set to be read in */
+ /* We need to allocate space for the data: */
+ snew(edi_read, 1);
+ /* Point the 'next_edi' entry to the next edi: */
+ curr_edi->next_edi = edi_read;
+ /* Keep the curr_edi pointer for the case that the next group is empty: */
+ last_edi = curr_edi;
+ /* Let's prepare to read in the next edi data set: */
+ curr_edi = edi_read;
+ }
+ if (edi_nr == 0)
+ {
+ gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
+ }
+
+ /* Terminate the edi group list with a NULL pointer: */
+ last_edi->next_edi = NULL;
+
+ fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
+
+ /* Close the .edi file again */
+ gmx_fio_fclose(in);
+
+ return edi_nr;
+}
+
+
+struct t_fit_to_ref {
+ rvec *xcopy; /* Working copy of the positions in fit_to_reference */
+};
+
+/* Fit the current positions to the reference positions
+ * Do not actually do the fit, just return rotation and translation.
+ * Note that the COM of the reference structure was already put into
+ * the origin by init_edi. */
+static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
+ rvec transvec, /* The translation vector */
+ matrix rotmat, /* The rotation matrix */
+ t_edpar *edi) /* Just needed for do_edfit */
+{
+ rvec com; /* center of mass */
+ int i;
+ struct t_fit_to_ref *loc;
+
+
+ /* Allocate memory the first time this routine is called for each edi group */
+ if (NULL == edi->buf->fit_to_ref)
+ {
+ snew(edi->buf->fit_to_ref, 1);
+ snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
+ }
+ loc = edi->buf->fit_to_ref;
+
+ /* We do not touch the original positions but work on a copy. */
+ for (i = 0; i < edi->sref.nr; i++)
+ {
+ copy_rvec(xcoll[i], loc->xcopy[i]);
+ }
+
+ /* Calculate the center of mass */
+ get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
+
+ transvec[XX] = -com[XX];
+ transvec[YY] = -com[YY];
+ transvec[ZZ] = -com[ZZ];
+
+ /* Subtract the center of mass from the copy */
+ translate_x(loc->xcopy, edi->sref.nr, transvec);
+
+ /* Determine the rotation matrix */
+ do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
+}
+
+
+static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
+ int nat, /* How many positions are there? */
+ rvec transvec, /* The translation vector */
+ matrix rotmat) /* The rotation matrix */
+{
+ /* Translation */
+ translate_x(x, nat, transvec);
+
+ /* Rotation */
+ rotate_x(x, nat, rotmat);
+}
+
+
+/* Gets the rms deviation of the positions to the structure s */
+/* fit_to_structure has to be called before calling this routine! */
+static real rmsd_from_structure(rvec *x, /* The positions under consideration */
+ struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
+{
+ real rmsd = 0.0;
+ int i;
+
+
+ for (i = 0; i < s->nr; i++)
+ {
+ rmsd += distance2(s->x[i], x[i]);
+ }
+
+ rmsd /= (real) s->nr;
+ rmsd = sqrt(rmsd);
+
+ return rmsd;
+}
+
+
+void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
+{
+ t_edpar *edi;
+
+
+ if (ed->eEDtype != eEDnone)
+ {
+ /* Loop over ED groups */
+ edi = ed->edpar;
+ while (edi)
+ {
+ /* Local atoms of the reference structure (for fitting), need only be assembled
+ * if their indices differ from the average ones */
+ if (!edi->bRefEqAv)
+ {
+ dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
+ &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
+ }
+
+ /* Local atoms of the average structure (on these ED will be performed) */
+ dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
+ &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
+
+ /* Indicate that the ED shift vectors for this structure need to be updated
+ * at the next call to communicate_group_positions, since obviously we are in a NS step */
+ edi->buf->do_edsam->bUpdateShifts = TRUE;
+
+ /* Set the pointer to the next ED group (if any) */
+ edi = edi->next_edi;
+ }
+ }
+}
+
+
+static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
+{
+ int tx, ty, tz;
+
+
+ tx = is[XX];
+ ty = is[YY];
+ tz = is[ZZ];
+
+ if (TRICLINIC(box))
+ {
+ xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
+ xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
+ xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
+ }
+ else
+ {
+ xu[XX] = x[XX]-tx*box[XX][XX];
+ xu[YY] = x[YY]-ty*box[YY][YY];
+ xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
+ }
+}
+
+
+static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
+{
+ int i, j;
+ real proj, add;
+ rvec vec_dum;
+
+
+ /* loop over linfix vectors */
+ for (i = 0; i < edi->vecs.linfix.neig; i++)
+ {
+ /* calculate the projection */
+ proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
+
+ /* calculate the correction */
+ add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
+
+ /* apply the correction */
+ add /= edi->sav.sqrtm[i];
+ for (j = 0; j < edi->sav.nr; j++)
+ {
+ svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
+ rvec_inc(xcoll[j], vec_dum);
+ }
+ }
+}
+
+
+static void do_linacc(rvec *xcoll, t_edpar *edi)
+{
+ int i, j;
+ real proj, add;
+ rvec vec_dum;
+
+
+ /* loop over linacc vectors */
+ for (i = 0; i < edi->vecs.linacc.neig; i++)
+ {
+ /* calculate the projection */
+ proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
+
+ /* calculate the correction */
+ add = 0.0;
+ if (edi->vecs.linacc.stpsz[i] > 0.0)
+ {
+ if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
+ {
+ add = edi->vecs.linacc.refproj[i] - proj;
+ }
+ }
+ if (edi->vecs.linacc.stpsz[i] < 0.0)
+ {
+ if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
+ {
+ add = edi->vecs.linacc.refproj[i] - proj;
+ }
+ }
+
+ /* apply the correction */
+ add /= edi->sav.sqrtm[i];
+ for (j = 0; j < edi->sav.nr; j++)
+ {
+ svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
+ rvec_inc(xcoll[j], vec_dum);
+ }
+
+ /* new positions will act as reference */
+ edi->vecs.linacc.refproj[i] = proj + add;
+ }
+}
+
+
+static void do_radfix(rvec *xcoll, t_edpar *edi)
+{
+ int i, j;
+ real *proj, rad = 0.0, ratio;
+ rvec vec_dum;
+
+
+ if (edi->vecs.radfix.neig == 0)
+ {
+ return;
+ }
+
+ snew(proj, edi->vecs.radfix.neig);
+
+ /* loop over radfix vectors */
+ for (i = 0; i < edi->vecs.radfix.neig; i++)
+ {
+ /* calculate the projections, radius */
+ proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
+ rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
+ }
+
+ rad = sqrt(rad);
+ ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
+ edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
+
+ /* loop over radfix vectors */
+ for (i = 0; i < edi->vecs.radfix.neig; i++)
+ {
+ proj[i] -= edi->vecs.radfix.refproj[i];
+
+ /* apply the correction */
+ proj[i] /= edi->sav.sqrtm[i];
+ proj[i] *= ratio;
+ for (j = 0; j < edi->sav.nr; j++)
+ {
+ svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
+ rvec_inc(xcoll[j], vec_dum);
+ }
+ }
+
+ sfree(proj);
+}
+
+
+static void do_radacc(rvec *xcoll, t_edpar *edi)
+{
+ int i, j;
+ real *proj, rad = 0.0, ratio = 0.0;
+ rvec vec_dum;
+
+
+ if (edi->vecs.radacc.neig == 0)
+ {
+ return;
+ }
+
+ snew(proj, edi->vecs.radacc.neig);
+
+ /* loop over radacc vectors */
+ for (i = 0; i < edi->vecs.radacc.neig; i++)
+ {
+ /* calculate the projections, radius */
+ proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
+ rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
+ }
+ rad = sqrt(rad);
+
+ /* only correct when radius decreased */
+ if (rad < edi->vecs.radacc.radius)
+ {
+ ratio = edi->vecs.radacc.radius/rad - 1.0;
+ rad = edi->vecs.radacc.radius;
+ }
+ else
+ {
+ edi->vecs.radacc.radius = rad;
+ }
+
+ /* loop over radacc vectors */
+ for (i = 0; i < edi->vecs.radacc.neig; i++)
+ {
+ proj[i] -= edi->vecs.radacc.refproj[i];
+
+ /* apply the correction */
+ proj[i] /= edi->sav.sqrtm[i];
+ proj[i] *= ratio;
+ for (j = 0; j < edi->sav.nr; j++)
+ {
+ svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
+ rvec_inc(xcoll[j], vec_dum);
+ }
+ }
+ sfree(proj);
+}
+
+
+struct t_do_radcon {
+ real *proj;
+};
+
+static void do_radcon(rvec *xcoll, t_edpar *edi)
+{
+ int i, j;
+ real rad = 0.0, ratio = 0.0;
+ struct t_do_radcon *loc;
+ gmx_bool bFirst;
+ rvec vec_dum;
+
+
+ if (edi->buf->do_radcon != NULL)
+ {
+ bFirst = FALSE;
+ loc = edi->buf->do_radcon;
+ }
+ else
+ {
+ bFirst = TRUE;
+ snew(edi->buf->do_radcon, 1);
+ }
+ loc = edi->buf->do_radcon;
+
+ if (edi->vecs.radcon.neig == 0)
+ {
+ return;
+ }
+
+ if (bFirst)
+ {
+ snew(loc->proj, edi->vecs.radcon.neig);
+ }
+
+ /* loop over radcon vectors */
+ for (i = 0; i < edi->vecs.radcon.neig; i++)
+ {
+ /* calculate the projections, radius */
+ loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
+ rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
+ }
+ rad = sqrt(rad);
+ /* only correct when radius increased */
+ if (rad > edi->vecs.radcon.radius)
+ {
+ ratio = edi->vecs.radcon.radius/rad - 1.0;
+
+ /* loop over radcon vectors */
+ for (i = 0; i < edi->vecs.radcon.neig; i++)
+ {
+ /* apply the correction */
+ loc->proj[i] -= edi->vecs.radcon.refproj[i];
+ loc->proj[i] /= edi->sav.sqrtm[i];
+ loc->proj[i] *= ratio;
+
+ for (j = 0; j < edi->sav.nr; j++)
+ {
+ svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
+ rvec_inc(xcoll[j], vec_dum);
+ }
+ }
+ }
+ else
+ {
+ edi->vecs.radcon.radius = rad;
+ }
+
+ if (rad != edi->vecs.radcon.radius)
+ {
+ rad = 0.0;
+ for (i = 0; i < edi->vecs.radcon.neig; i++)
+ {
+ /* calculate the projections, radius */
+ loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
+ rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
+ }
+ rad = sqrt(rad);
+ }
+}
+
+
+static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
+{
+ int i;
+
+
+ /* subtract the average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_dec(xcoll[i], edi->sav.x[i]);
+ }
+
+ /* apply the constraints */
+ if (step >= 0)
+ {
+ do_linfix(xcoll, edi, step);
+ }
+ do_linacc(xcoll, edi);
+ if (step >= 0)
+ {
+ do_radfix(xcoll, edi);
+ }
+ do_radacc(xcoll, edi);
+ do_radcon(xcoll, edi);
+
+ /* add back the average positions */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ rvec_inc(xcoll[i], edi->sav.x[i]);
+ }
+}
+
+
+/* Write out the projections onto the eigenvectors. The order of output
+ * corresponds to ed_output_legend() */
+static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
+{
+ int i;
+
+
+ /* Output how well we fit to the reference structure */
+ fprintf(fp, EDcol_ffmt, rmsd);
+
+ for (i = 0; i < edi->vecs.mon.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
+ }
+
+ for (i = 0; i < edi->vecs.linfix.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
+ }
+
+ for (i = 0; i < edi->vecs.linacc.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
+ }
+
+ for (i = 0; i < edi->vecs.radfix.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
+ }
+ if (edi->vecs.radfix.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
+ }
+
+ for (i = 0; i < edi->vecs.radacc.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
+ }
+ if (edi->vecs.radacc.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
+ }
+
+ for (i = 0; i < edi->vecs.radcon.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
+ }
+ if (edi->vecs.radcon.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
+ }
+}
+
+/* Returns if any constraints are switched on */
+static int ed_constraints(gmx_bool edtype, t_edpar *edi)
+{
+ if (edtype == eEDedsam || edtype == eEDflood)
+ {
+ return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
+ edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
+ edi->vecs.radcon.neig);
+ }
+ return 0;
+}
+
+
+/* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
+ * umbrella sampling simulations. */
+static void copyEvecReference(t_eigvec* floodvecs)
+{
+ int i;
+
+
+ if (NULL == floodvecs->refproj0)
+ {
+ snew(floodvecs->refproj0, floodvecs->neig);
+ }
+
+ for (i = 0; i < floodvecs->neig; i++)
+ {
+ floodvecs->refproj0[i] = floodvecs->refproj[i];
+ }
+}
+
+
+/* Call on MASTER only. Check whether the essential dynamics / flooding
+ * groups of the checkpoint file are consistent with the provided .edi file. */
+static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
+{
+ t_edpar *edi = NULL; /* points to a single edi data set */
+ int edinum;
+
+
+ if (NULL == EDstate->nref || NULL == EDstate->nav)
+ {
+ gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
+ "start of a new simulation. If a simulation runs with/without ED constraints,\n"
+ "it must also continue with/without ED constraints when checkpointing.\n"
+ "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
+ "from without a checkpoint.\n");
+ }
+
+ edi = ed->edpar;
+ edinum = 0;
+ while (edi != NULL)
+ {
+ /* Check number of atoms in the reference and average structures */
+ if (EDstate->nref[edinum] != edi->sref.nr)
+ {
+ gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
+ "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
+ get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
+ }
+ if (EDstate->nav[edinum] != edi->sav.nr)
+ {
+ gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
+ "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
+ get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
+ }
+ edi = edi->next_edi;
+ edinum++;
+ }
+
+ if (edinum != EDstate->nED)
+ {
+ gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
+ "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
+ "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
+ }
+}
+
+
+/* The edsamstate struct stores the information we need to make the ED group
+ * whole again after restarts from a checkpoint file. Here we do the following:
+ * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
+ * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
+ * c) in any case, for subsequent checkpoint writing, we set the pointers in
+ * edsamstate to the x_old arrays, which contain the correct PBC representation of
+ * all ED structures at the last time step. */
+static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
+{
+ int i, nr_edi;
+ t_edpar *edi;
+
+
+ snew(EDstate->old_sref_p, EDstate->nED);
+ snew(EDstate->old_sav_p, EDstate->nED);
+
+ /* If we did not read in a .cpt file, these arrays are not yet allocated */
+ if (!EDstate->bFromCpt)
+ {
+ snew(EDstate->nref, EDstate->nED);
+ snew(EDstate->nav, EDstate->nED);
+ }
+
+ /* Loop over all ED/flooding data sets (usually only one, though) */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
+ {
+ /* We always need the last reference and average positions such that
+ * in the next time step we can make the ED group whole again
+ * if the atoms do not have the correct PBC representation */
+ if (EDstate->bFromCpt)
+ {
+ /* Copy the last whole positions of reference and average group from .cpt */
+ for (i = 0; i < edi->sref.nr; i++)
+ {
+ copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
+ }
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
+ }
+ }
+ else
+ {
+ EDstate->nref[nr_edi-1] = edi->sref.nr;
+ EDstate->nav [nr_edi-1] = edi->sav.nr;
+ }
+
+ /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
+ EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
+ EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
+
+ edi = edi->next_edi;
+ }
+}
+
+
+/* Adds 'buf' to 'str' */
+static void add_to_string(char **str, char *buf)
+{
+ int len;
+
+
+ len = strlen(*str) + strlen(buf) + 1;
+ srenew(*str, len);
+ strcat(*str, buf);
+}
+
+
+static void add_to_string_aligned(char **str, char *buf)
+{
+ char buf_aligned[STRLEN];
+
+ sprintf(buf_aligned, EDcol_sfmt, buf);
+ add_to_string(str, buf_aligned);
+}
+
+
+static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
+{
+ char tmp[STRLEN], tmp2[STRLEN];
+
+
+ sprintf(tmp, "%c %s", EDgroupchar, value);
+ add_to_string_aligned(LegendStr, tmp);
+ sprintf(tmp2, "%s (%s)", tmp, unit);
+ (*setname)[*nsets] = strdup(tmp2);
+ (*nsets)++;
+}
+
+
+static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
+{
+ int i;
+ char tmp[STRLEN];
+
+
+ for (i = 0; i < evec->neig; i++)
+ {
+ sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
+ nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
+ }
+}
+
+
+/* Makes a legend for the xvg output file. Call on MASTER only! */
+static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
+{
+ t_edpar *edi = NULL;
+ int i;
+ int nr_edi, nsets, n_flood, n_edsam;
+ const char **setname;
+ char buf[STRLEN];
+ char *LegendStr = NULL;
+
+
+ edi = ed->edpar;
+
+ fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
+
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
- if (edi->bNeedDoEdsam) /* Only print ED legend if at least one ED option is on */
+ fprintf(ed->edo, "#\n");
+ fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
+ fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
+ fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
+ fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
+
+ if (edi->flood.vecs.neig)
+ {
+ /* If in any of the groups we find a flooding vector, flooding is turned on */
+ ed->eEDtype = eEDflood;
+
+ /* Print what flavor of flooding we will do */
+ if (0 == edi->flood.tau) /* constant flooding strength */
+ {
+ fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
+ if (edi->flood.bHarmonic)
+ {
+ fprintf(ed->edo, ", harmonic");
+ }
+ }
+ else /* adaptive flooding */
+ {
+ fprintf(ed->edo, ", adaptive");
+ }
+ }
+ fprintf(ed->edo, "\n");
+
+ edi = edi->next_edi;
+ }
+
+ /* Print a nice legend */
+ snew(LegendStr, 1);
+ LegendStr[0] = '\0';
+ sprintf(buf, "# %6s", "time");
+ add_to_string(&LegendStr, buf);
+
+ /* Calculate the maximum number of columns we could end up with */
+ edi = ed->edpar;
+ nsets = 0;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ nsets += 5 +edi->vecs.mon.neig
+ +edi->vecs.linfix.neig
+ +edi->vecs.linacc.neig
+ +edi->vecs.radfix.neig
+ +edi->vecs.radacc.neig
+ +edi->vecs.radcon.neig
+ + 6*edi->flood.vecs.neig;
+ edi = edi->next_edi;
+ }
+ snew(setname, nsets);
+
+ /* In the mdrun time step in a first function call (do_flood()) the flooding
+ * forces are calculated and in a second function call (do_edsam()) the
+ * ED constraints. To get a corresponding legend, we need to loop twice
+ * over the edi groups and output first the flooding, then the ED part */
+
+ /* The flooding-related legend entries, if flooding is done */
+ nsets = 0;
+ if (eEDflood == ed->eEDtype)
+ {
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ /* Always write out the projection on the flooding EVs. Of course, this can also
+ * be achieved with the monitoring option in do_edsam() (if switched on by the
+ * user), but in that case the positions need to be communicated in do_edsam(),
+ * which is not necessary when doing flooding only. */
+ nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
+
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
+
+ /* Output the current reference projection if it changes with time;
+ * this can happen when flooding is used as harmonic restraint */
+ if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
+ {
+ sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
+ }
+
+ /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
+ if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
+ {
+ sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+ }
+
+ sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+
+ if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
+ {
+ sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+ }
+
+ sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
+ }
+
+ edi = edi->next_edi;
+ } /* End of flooding-related legend entries */
+ }
+ n_flood = nsets;
+
+ /* Now the ED-related entries, if essential dynamics is done */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
- /* Remove pbc, make molecule whole.
- * When ir->bContinuation=TRUE this has already been done, but ok.
- */
- snew(x_pbc, mtop->natoms);
- m_rveccopy(mtop->natoms, x, x_pbc);
- do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
-
++ if ( bNeedDoEdsam(edi) ) /* Only print ED legend if at least one ED option is on */
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
+
+ /* Essential dynamics, projections on eigenvectors */
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
+ if (edi->vecs.radfix.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
+ if (edi->vecs.radacc.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
+ if (edi->vecs.radcon.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+ }
+ edi = edi->next_edi;
+ } /* end of 'pure' essential dynamics legend entries */
+ n_edsam = nsets - n_flood;
+
+ xvgr_legend(ed->edo, nsets, setname, oenv);
+ sfree(setname);
+
+ fprintf(ed->edo, "#\n"
+ "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
+ n_flood, 1 == n_flood ? "" : "s",
+ n_edsam, 1 == n_edsam ? "" : "s");
+ fprintf(ed->edo, "%s", LegendStr);
+ sfree(LegendStr);
+
+ fflush(ed->edo);
+}
+
+
+void init_edsam(gmx_mtop_t *mtop, /* global topology */
+ t_inputrec *ir, /* input record */
+ t_commrec *cr, /* communication record */
+ gmx_edsam_t ed, /* contains all ED data */
+ rvec x[], /* positions of the whole MD system */
+ matrix box, /* the box */
+ edsamstate_t *EDstate)
+{
+ t_edpar *edi = NULL; /* points to a single edi data set */
+ int i, nr_edi, avindex;
+ rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
+ rvec *xfit = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs */
+ rvec fit_transvec; /* translation ... */
+ matrix fit_rotmat; /* ... and rotation from fit to reference structure */
++ rvec *ref_x_old = NULL; /* helper pointer */
+
+
+ if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
+ {
+ gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
+ }
+
+ if (MASTER(cr))
+ {
+ fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
+
+ if (NULL == ed)
+ {
+ gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
+ "flooding simulation. Please also provide the correct .edi file with -ei.\n");
+ }
+ }
+
+ /* Needed for initializing radacc radius in do_edsam */
+ ed->bFirst = TRUE;
+
+ /* The input file is read by the master and the edi structures are
+ * initialized here. Input is stored in ed->edpar. Then the edi
+ * structures are transferred to the other nodes */
+ if (MASTER(cr))
+ {
+ /* Initialization for every ED/flooding group. Flooding uses one edi group per
+ * flooding vector, Essential dynamics can be applied to more than one structure
+ * as well, but will be done in the order given in the edi file, so
+ * expect different results for different order of edi file concatenation! */
+ edi = ed->edpar;
+ while (edi != NULL)
+ {
+ init_edi(mtop, edi);
+ init_flood(edi, ed, ir->delta_t);
+ edi = edi->next_edi;
+ }
+ }
+
+ /* The master does the work here. The other nodes get the positions
+ * not before dd_partition_system which is called after init_edsam */
+ if (MASTER(cr))
+ {
- copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
++ if (!EDstate->bFromCpt)
++ {
++ /* Remove PBC, make molecule(s) subject to ED whole. */
++ snew(x_pbc, mtop->natoms);
++ m_rveccopy(mtop->natoms, x, x_pbc);
++ do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
++ }
+ /* Reset pointer to first ED data set which contains the actual ED data */
+ edi = ed->edpar;
+ /* Loop over all ED/flooding data sets (usually only one, though) */
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
+ {
+ /* For multiple ED groups we use the output frequency that was specified
+ * in the first set */
+ if (nr_edi > 1)
+ {
+ edi->outfrq = ed->edpar->outfrq;
+ }
+
+ /* Extract the initial reference and average positions. When starting
+ * from .cpt, these have already been read into sref.x_old
+ * in init_edsamstate() */
+ if (!EDstate->bFromCpt)
+ {
+ /* If this is the first run (i.e. no checkpoint present) we assume
+ * that the starting positions give us the correct PBC representation */
+ for (i = 0; i < edi->sref.nr; i++)
+ {
+ copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
+ }
+
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
+ }
+ }
+
+ /* Now we have the PBC-correct start positions of the reference and
+ average structure. We copy that over to dummy arrays on which we
+ can apply fitting to print out the RMSD. We srenew the memory since
+ the size of the buffers is likely different for every ED group */
+ srenew(xfit, edi->sref.nr );
+ srenew(xstart, edi->sav.nr );
- sfree(x_pbc);
++ if (edi->bRefEqAv)
++ {
++ /* Reference indices are the same as average indices */
++ ref_x_old = edi->sav.x_old;
++ }
++ else
++ {
++ ref_x_old = edi->sref.x_old;
++ }
++ copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
+ copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
+
+ /* Make the fit to the REFERENCE structure, get translation and rotation */
+ fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
+
+ /* Output how well we fit to the reference at the start */
+ translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
+ fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
+ rmsd_from_structure(xfit, &edi->sref));
+ if (EDstate->nED > 1)
+ {
+ fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
+ }
+ fprintf(stderr, "\n");
+
+ /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
+ translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
+
+ /* calculate initial projections */
+ project(xstart, edi);
+
+ /* For the target and origin structure both a reference (fit) and an
+ * average structure can be provided in make_edi. If both structures
+ * are the same, make_edi only stores one of them in the .edi file.
+ * If they differ, first the fit and then the average structure is stored
+ * in star (or sor), thus the number of entries in star/sor is
+ * (n_fit + n_av) with n_fit the size of the fitting group and n_av
+ * the size of the average group. */
+
+ /* process target structure, if required */
+ if (edi->star.nr > 0)
+ {
+ fprintf(stderr, "ED: Fitting target structure to reference structure\n");
+
+ /* get translation & rotation for fit of target structure to reference structure */
+ fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
+ /* do the fit */
+ translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
+ if (edi->star.nr == edi->sav.nr)
+ {
+ avindex = 0;
+ }
+ else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
+ {
+ /* The last sav.nr indices of the target structure correspond to
+ * the average structure, which must be projected */
+ avindex = edi->star.nr - edi->sav.nr;
+ }
+ rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
+ }
+ else
+ {
+ rad_project(edi, xstart, &edi->vecs.radcon);
+ }
+
+ /* process structure that will serve as origin of expansion circle */
+ if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
+ {
+ fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
+ }
+
+ if (edi->sori.nr > 0)
+ {
+ fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
+
+ /* fit this structure to reference structure */
+ fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
+ /* do the fit */
+ translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
+ if (edi->sori.nr == edi->sav.nr)
+ {
+ avindex = 0;
+ }
+ else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
+ {
+ /* For the projection, we need the last sav.nr indices of sori */
+ avindex = edi->sori.nr - edi->sav.nr;
+ }
+
+ rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
+ rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
+ if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
+ {
+ fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
+ /* Set center of flooding potential to the ORIGIN structure */
+ rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
+ /* We already know that no (moving) reference position was provided,
+ * therefore we can overwrite refproj[0]*/
+ copyEvecReference(&edi->flood.vecs);
+ }
+ }
+ else /* No origin structure given */
+ {
+ rad_project(edi, xstart, &edi->vecs.radacc);
+ rad_project(edi, xstart, &edi->vecs.radfix);
+ if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
+ {
+ if (edi->flood.bHarmonic)
+ {
+ fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
+ }
+ }
+ else
+ {
+ fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
+ /* Set center of flooding potential to the center of the covariance matrix,
+ * i.e. the average structure, i.e. zero in the projected system */
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ edi->flood.vecs.refproj[i] = 0.0;
+ }
+ }
+ }
+ }
+ /* For convenience, output the center of the flooding potential for the eigenvectors */
+ if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
+ {
+ for (i = 0; i < edi->flood.vecs.neig; i++)
+ {
+ fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
+ if (edi->flood.bHarmonic)
+ {
+ fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
+ }
+ fprintf(stdout, "\n");
+ }
+ }
+
+ /* set starting projections for linsam */
+ rad_project(edi, xstart, &edi->vecs.linacc);
+ rad_project(edi, xstart, &edi->vecs.linfix);
+
+ /* Prepare for the next edi data set: */
+ edi = edi->next_edi;
+ }
+ /* Cleaning up on the master node: */
- /* Allocate space for ED buffer */
- snew(edi->buf, 1);
++ if (!EDstate->bFromCpt)
++ {
++ sfree(x_pbc);
++ }
+ sfree(xfit);
+ sfree(xstart);
+
+ } /* end of MASTER only section */
+
+ if (PAR(cr))
+ {
+ /* First let everybody know how many ED data sets to expect */
+ gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
+ /* Broadcast the essential dynamics / flooding data to all nodes */
+ broadcast_ed_data(cr, ed, EDstate->nED);
+ }
+ else
+ {
+ /* In the single-CPU case, point the local atom numbers pointers to the global
+ * one, so that we can use the same notation in serial and parallel case: */
+
+ /* Loop over all ED data sets (usually only one, though) */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
+ {
+ edi->sref.anrs_loc = edi->sref.anrs;
+ edi->sav.anrs_loc = edi->sav.anrs;
+ edi->star.anrs_loc = edi->star.anrs;
+ edi->sori.anrs_loc = edi->sori.anrs;
+ /* For the same reason as above, make a dummy c_ind array: */
+ snew(edi->sav.c_ind, edi->sav.nr);
+ /* Initialize the array */
+ for (i = 0; i < edi->sav.nr; i++)
+ {
+ edi->sav.c_ind[i] = i;
+ }
+ /* In the general case we will need a different-sized array for the reference indices: */
+ if (!edi->bRefEqAv)
+ {
+ snew(edi->sref.c_ind, edi->sref.nr);
+ for (i = 0; i < edi->sref.nr; i++)
+ {
+ edi->sref.c_ind[i] = i;
+ }
+ }
+ /* Point to the very same array in case of other structures: */
+ edi->star.c_ind = edi->sav.c_ind;
+ edi->sori.c_ind = edi->sav.c_ind;
+ /* In the serial case, the local number of atoms is the global one: */
+ edi->sref.nr_loc = edi->sref.nr;
+ edi->sav.nr_loc = edi->sav.nr;
+ edi->star.nr_loc = edi->star.nr;
+ edi->sori.nr_loc = edi->sori.nr;
+
+ /* An on we go to the next ED group */
+ edi = edi->next_edi;
+ }
+ }
+
+ /* Allocate space for ED buffer variables */
+ /* Again, loop over ED data sets */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
+ {
- if (edi->bNeedDoEdsam)
++ /* Allocate space for ED buffer variables */
++ snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
+ snew(edi->buf->do_edsam, 1);
+
+ /* Space for collective ED buffer variables */
+
+ /* Collective positions of atoms with the average indices */
+ snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
+ snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
+ snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
+ /* Collective positions of atoms with the reference indices */
+ if (!edi->bRefEqAv)
+ {
+ snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
+ snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
+ snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
+ }
+
+ /* Get memory for flooding forces */
+ snew(edi->flood.forces_cartesian, edi->sav.nr);
+
+#ifdef DUMPEDI
+ /* Dump it all into one file per process */
+ dump_edi(edi, cr, nr_edi);
+#endif
+
+ /* Next ED group */
+ edi = edi->next_edi;
+ }
+
+ /* Flush the edo file so that the user can check some things
+ * when the simulation has started */
+ if (ed->edo)
+ {
+ fflush(ed->edo);
+ }
+}
+
+
+void do_edsam(t_inputrec *ir,
+ gmx_int64_t step,
+ t_commrec *cr,
+ rvec xs[], /* The local current positions on this processor */
+ rvec v[], /* The velocities */
+ matrix box,
+ gmx_edsam_t ed)
+{
+ int i, edinr, iupdate = 500;
+ matrix rotmat; /* rotation matrix */
+ rvec transvec; /* translation vector */
+ rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
+ real dt_1; /* 1/dt */
+ struct t_do_edsam *buf;
+ t_edpar *edi;
+ real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
+ gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
+
+
+ /* Check if ED sampling has to be performed */
+ if (ed->eEDtype == eEDnone)
+ {
+ return;
+ }
+
+ /* Suppress output on first call of do_edsam if
+ * two-step sd2 integrator is used */
+ if ( (ir->eI == eiSD2) && (v != NULL) )
+ {
+ bSuppress = TRUE;
+ }
+
+ dt_1 = 1.0/ir->delta_t;
+
+ /* Loop over all ED groups (usually one) */
+ edi = ed->edpar;
+ edinr = 0;
+ while (edi != NULL)
+ {
+ edinr++;
- /* initialise radacc radius for slope criterion */
++ if ( bNeedDoEdsam(edi) )
+ {
+
+ buf = edi->buf->do_edsam;
+
+ if (ed->bFirst)
+ {
- } /* END of if (edi->bNeedDoEdsam) */
++ /* initialize radacc radius for slope criterion */
+ buf->oldrad = calc_radius(&edi->vecs.radacc);
+ }
+
+ /* Copy the positions into buf->xc* arrays and after ED
+ * feed back corrections to the official positions */
+
+ /* Broadcast the ED positions such that every node has all of them
+ * Every node contributes its local positions xs and stores it in
+ * the collective buf->xcoll array. Note that for edinr > 1
+ * xs could already have been modified by an earlier ED */
+
+ communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
+ edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
+
+ /* Only assembly reference positions if their indices differ from the average ones */
+ if (!edi->bRefEqAv)
+ {
+ communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
+ edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
+ }
+
+ /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
+ * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
+ * set bUpdateShifts=TRUE in the parallel case. */
+ buf->bUpdateShifts = FALSE;
+
+ /* Now all nodes have all of the ED positions in edi->sav->xcoll,
+ * as well as the indices in edi->sav.anrs */
+
+ /* Fit the reference indices to the reference structure */
+ if (edi->bRefEqAv)
+ {
+ fit_to_reference(buf->xcoll, transvec, rotmat, edi);
+ }
+ else
+ {
+ fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
+ }
+
+ /* Now apply the translation and rotation to the ED structure */
+ translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
+
+ /* Find out how well we fit to the reference (just for output steps) */
+ if (do_per_step(step, edi->outfrq) && MASTER(cr))
+ {
+ if (edi->bRefEqAv)
+ {
+ /* Indices of reference and average structures are identical,
+ * thus we can calculate the rmsd to SREF using xcoll */
+ rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
+ }
+ else
+ {
+ /* We have to translate & rotate the reference atoms first */
+ translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
+ rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
+ }
+ }
+
+ /* update radsam references, when required */
+ if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
+ {
+ project(buf->xcoll, edi);
+ rad_project(edi, buf->xcoll, &edi->vecs.radacc);
+ rad_project(edi, buf->xcoll, &edi->vecs.radfix);
+ buf->oldrad = -1.e5;
+ }
+
+ /* update radacc references, when required */
+ if (do_per_step(step, iupdate) && step >= edi->presteps)
+ {
+ edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
+ if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
+ {
+ project(buf->xcoll, edi);
+ rad_project(edi, buf->xcoll, &edi->vecs.radacc);
+ buf->oldrad = 0.0;
+ }
+ else
+ {
+ buf->oldrad = edi->vecs.radacc.radius;
+ }
+ }
+
+ /* apply the constraints */
+ if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
+ {
+ /* ED constraints should be applied already in the first MD step
+ * (which is step 0), therefore we pass step+1 to the routine */
+ ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
+ }
+
+ /* write to edo, when required */
+ if (do_per_step(step, edi->outfrq))
+ {
+ project(buf->xcoll, edi);
+ if (MASTER(cr) && !bSuppress)
+ {
+ write_edo(edi, ed->edo, rmsdev);
+ }
+ }
+
+ /* Copy back the positions unless monitoring only */
+ if (ed_constraints(ed->eEDtype, edi))
+ {
+ /* remove fitting */
+ rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
+
+ /* Copy the ED corrected positions into the coordinate array */
+ /* Each node copies its local part. In the serial case, nat_loc is the
+ * total number of ED atoms */
+ for (i = 0; i < edi->sav.nr_loc; i++)
+ {
+ /* Unshift local ED coordinate and store in x_unsh */
+ ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
+ buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
+
+ /* dx is the ED correction to the positions: */
+ rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
+
+ if (v != NULL)
+ {
+ /* dv is the ED correction to the velocity: */
+ svmul(dt_1, dx, dv);
+ /* apply the velocity correction: */
+ rvec_inc(v[edi->sav.anrs_loc[i]], dv);
+ }
+ /* Finally apply the position correction due to ED: */
+ copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
+ }
+ }
++ } /* END of if ( bNeedDoEdsam(edi) ) */
+
+ /* Prepare for the next ED group */
+ edi = edi->next_edi;
+
+ } /* END of loop over ED groups */
+
+ ed->bFirst = FALSE;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff