/* -*- 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.
* 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
*/
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 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 eigenvector */
+ 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 kT;
real Vfl;
real dt;
real constEfl;
- real alpha2;
+ real alpha2;
int flood_id;
- rvec *forces_cartesian;
+ rvec *forces_cartesian;
t_eigvec vecs; /* use flooding for these */
} t_edflood;
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
+ 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; /* used to keep track of the shift vectors
- such that the ED molecule can always be
+ such that the ED molecule can always be
made whole in the parallel case */
real *m; /* masses */
real mtot; /* total mass (only used in sref) */
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
+ 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 */
matrix old_rotmat;
real oldrad;
rvec old_transvec,older_transvec,transvec_compact;
- rvec *xcoll; /* Positions from all nodes, this is the
+ rvec *xcoll; /* Positions from all nodes, this is the
collective set we work on.
- These are the positions of atoms with
+ 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 *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 dataset need to
+ ED shifts for this ED dataset need to
be updated */
};
rad += pow((vec->refproj[i]-vec->xproj[i]),2);
}
vec->radius=sqrt(rad);
-
+
/* Add average positions */
- for (i = 0; i < edi->sav.nr; i++)
+ 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
+/* 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 */
+static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
t_eigvec *vec, /* The eigenvectors */
t_edpar *edi)
{
if (!vec->neig) return;
/* Subtract average positions */
- for (i=0; i<edi->sav.nr; i++)
+ 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++)
+ 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 */
+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
real rad=0.0;
- for (i=0; i<vec->neig; i++)
+ 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,
+static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
int step)
{
int i;
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,
+ 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],
+ shifts[i][XX] , shifts[i][YY] , shifts[i][ZZ],
eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
fclose(fp);
{
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++)
+ 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)
{
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",
+ 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_rvec(FILE *out, int dim, rvec *x)
+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]);
}
{
int i,j;
-
+
fprintf(out,"MATRIX:\n");
for (i=0;i<dim;i++)
{
#endif
-struct t_do_edfit {
+struct t_do_edfit {
double **omega;
double **om;
};
}
loc = edi->buf->do_edfit;
- if (bFirst)
+ if (bFirst)
{
snew(loc->omega,2*DIM);
snew(loc->om,2*DIM);
}
-static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
+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 */
/**********************************************************************************
******************** 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).
+
+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
+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
+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.
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.
+ 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.
+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.
+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
+ 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);
+ 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
+ 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
*/
char buf[22];
gmx_bool bOutputRef=FALSE;
-
+
fprintf(fp,"%d.th FL: %s %12.5e %12.5e %12.5e\n",
edi->flood.flood_id, gmx_step_str(step,buf),
edi->flood.Efl, edi->flood.Vfl, edi->flood.deltaF);
}
}
fprintf(fp,"FL_FORCES: ");
-
+
for (i=0; i<edi->flood.vecs.neig; i++)
fprintf(fp," %12.5e",edi->flood.vecs.fproj[i]);
-
+
fprintf(fp,"\n");
}
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)
{
{
/* 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;
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
+ 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 momentarily we want to compute the forces separately
+ 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++)
+ 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++)
{
/* 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,
+ /* 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)
gmx_large_int_t step,
matrix box,
t_commrec *cr)
-{
+{
int i;
matrix rotmat; /* rotation matrix */
matrix tmat; /* inverse rotation */
rvec transvec; /* translation vector */
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, buf->bUpdateShifts, x,
- edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
-
+ * the collective ED array buf->xcoll */
+ communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, 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, buf->bUpdateShifts, x,
+ communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, 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.
/* 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);
-
- /* Compute Vfl(x) from flood.xproj */
- edi->flood.Vfl = flood_energy(edi, step);
-
- update_adaption(edi);
+ 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);
- /* Compute the flooding forces */
- flood_forces(edi);
+ update_adaption(edi);
+
+ /* Compute the flooding forces */
+ flood_forces(edi);
+ }
/* Translate them into cartesian positions */
flood_blowup(edi, edi->flood.forces_cartesian);
{
t_edpar *edi;
-
+
if (ed->eEDtype != eEDflood)
return;
-
+
edi = ed->edpar;
while (edi)
{
}
-/* Called by init_edi, configure some flooding related variables and structures,
+/* 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, t_commrec *cr)
{
+ int i;
+
+
edi->flood.Efl = edi->flood.constEfl;
edi->flood.Vfl = 0;
edi->flood.dt = dt;
{
/* If in any of the datasets we find a flooding vector, flooding is turned on */
ed->eEDtype = eEDflood;
-
+
+ fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
+
+ 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]);
+ }
+ }
fprintf(ed->edo,"FL_HEADER: Flooding of matrix %d is switched on! The flooding output will have the following format:\n",
edi->flood.flood_id);
- fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
- if (edi->flood.flood_id<1)
- fprintf(ed->edo,"FL_HEADER: Step Efl Vfl deltaF\n");
+ fprintf(ed->edo,"FL_HEADER: Step Efl Vfl deltaF\n");
}
}
t_edpar *actual;
int count;
-
+
actual=edi;
count = 1;
while (actual)
actual=actual->next_edi;
count++;
}
- if (nnames!=count-1)
+ if (nnames!=count-1)
gmx_fatal(FARGS,"Number of energies is not consistent with t_edi structure");
}
/************* END of FLOODING IMPLEMENTATION ****************************/
gmx_edsam_t ed_open(int nfile,const t_filenm fnm[],unsigned long Flags,t_commrec *cr)
-{
+{
gmx_edsam_t ed;
-
-
+
+
/* 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))
+ if (MASTER(cr))
{
/* Open .edi input file: */
ed->edinam=ftp2fn(efEDI,nfile,fnm);
/* The master opens the .edo output file */
- fprintf(stderr,"ED sampling will be performed!\n");
+ fprintf(stderr,"ED sampling will be performed!\n");
ed->edonam = ftp2fn(efEDO,nfile,fnm);
ed->edo = gmx_fio_fopen(ed->edonam,(Flags & MD_APPENDFILES)? "a+" : "w+");
ed->bStartFromCpt = Flags & MD_STARTFROMCPT;
{
/* 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.
+ /* 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, ... */
+ 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 */
}
{
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);
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 dataset */
if (edi->next_edi)
{
/* init-routine called for every *.edi-cycle, initialises t_edpar structure */
static void init_edi(gmx_mtop_t *mtop,t_inputrec *ir,
- t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
+ t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
{
int i;
real totalmass = 0.0;
rvec com;
t_atom *atom;
- /* NOTE Init_edi is executed on the master process only
+ /* 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 */
- edi->bNeedDoEdsam = edi->vecs.mon.neig
- || edi->vecs.linfix.neig
- || edi->vecs.linacc.neig
- || edi->vecs.radfix.neig
- || edi->vecs.radacc.neig
+ 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;
/* evaluate masses (reference structure) */
{
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
+ /* 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 );
{
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);
+ }
}
/* put reference structure in origin */
static void check(const char *line, const char *label)
{
- if (!strstr(line,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);
}
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)
int idum;
char *eof;
-
+
eof=fgets2 (line,STRLEN,file);
if (eof==NULL)
{
*bEOF = TRUE;
return -1;
}
- eof=fgets2 (line,STRLEN,file);
+ eof=fgets2 (line,STRLEN,file);
if (eof==NULL)
{
*bEOF = TRUE;
char line[STRLEN+1];
double rdum;
-
+
fgets2 (line,STRLEN,file);
check(line,label);
fgets2 (line,STRLEN,file);
char line[STRLEN+1];
double d[3];
-
+
for(i=0; i<number; i++)
{
fgets2 (line,STRLEN,file);
int i;
double x,y,z;
-
+
for(i=0; (i < nr); i++)
{
fgets2 (line,STRLEN,in);
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)
{
{
int i;
-
+
/* If the number of atoms differs between the two structures,
* they cannot be identical */
if (sref.nr != sav.nr)
static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
{
int readmagic;
- const int magic=669;
+ const int magic=670;
gmx_bool bEOF;
-
+
/* the edi file is not free format, so expect problems if the input is corrupt. */
/* check the magic number */
{
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)
+ ;
else
gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
}
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)",
- ed->edinam,edi->nini,nr_mdatoms);
+ ed->edinam,edi->nini,nr_mdatoms);
/* Done checking. For the rest we blindly trust the input */
edi->fitmas = read_checked_edint(in,"FITMAS");
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->flood.flood_id = edi_nr;
edi->sref.nr = read_checked_edint(in,"NREF");
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 void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
+static void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
{
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(ed->edinam,"r");
+ in = gmx_fio_fopen(ed->edinam,"r");
fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
/* Now read a sequence of ED input parameter sets from the edi file */
edi_nr++;
/* Make shure that the number of atoms in each dataset is the same as in the tpr file */
if (edi->nini != nr_mdatoms)
- gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
+ gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
ed->edinam, edi_nr, edi->nini, nr_mdatoms);
/* 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: */
last_edi = curr_edi;
/* Let's prepare to read in the next edi data set: */
curr_edi = edi_read;
- }
- if (edi_nr == 0)
+ }
+ if (edi_nr == 0)
gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
/* Terminate the edi dataset list with a NULL pointer: */
last_edi->next_edi = NULL;
fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
-
+
/* Close the .edi file again */
gmx_fio_fclose(in);
}
rvec *xcopy; /* Working copy of the positions in fit_to_reference */
};
-/* Fit the current positions to the reference positions
+/* 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
+ * 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 */
+ 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;
-
+
GMX_MPE_LOG(ev_fit_to_reference_start);
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]);
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 */
+ rvec transvec, /* The translation vector */
matrix rotmat) /* The rotation matrix */
{
/* Translation */
/* 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 */
+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;
{
t_edpar *edi;
-
+
if (ed->eEDtype != eEDnone)
{
/* Loop over ED datasets (usually there is just one dataset, though) */
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
+ /* 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 dataset (if any) */
edi=edi->next_edi;
}
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];
rvec vec_dum;
- if (edi->vecs.radfix.neig == 0)
+ if (edi->vecs.radfix.neig == 0)
return;
snew(proj, edi->vecs.radfix.neig);
rvec vec_dum;
- if (edi->vecs.radacc.neig == 0)
+ if (edi->vecs.radacc.neig == 0)
return;
snew(proj,edi->vecs.radacc.neig);
ratio = edi->vecs.radacc.radius/rad - 1.0;
rad = edi->vecs.radacc.radius;
}
- else
+ else
edi->vecs.radacc.radius = rad;
/* loop over radacc vectors */
svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
rvec_inc(xcoll[j], vec_dum);
}
- }
+ }
sfree(proj);
}
{
bFirst = FALSE;
loc = edi->buf->do_radcon;
- }
+ }
else
{
bFirst = TRUE;
}
loc = edi->buf->do_radcon;
- if (edi->vecs.radcon.neig == 0)
+ if (edi->vecs.radcon.neig == 0)
return;
if (bFirst)
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_dec(xcoll[i], edi->sav.x[i]);
/* apply the constraints */
- if (step >= 0)
+ if (step >= 0)
do_linfix(xcoll, edi, step, cr);
do_linacc(xcoll, edi, cr);
- if (step >= 0)
+ if (step >= 0)
do_radfix(xcoll, edi, step, cr);
do_radacc(xcoll, edi, cr);
do_radcon(xcoll, edi, cr);
/* add back the average positions */
- for (i=0; i<edi->sav.nr; i++)
+ for (i=0; i<edi->sav.nr; i++)
rvec_inc(xcoll[i], edi->sav.x[i]);
GMX_MPE_LOG(ev_ed_apply_cons_finish);
/* Write out the projections onto the eigenvectors */
static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
-{
+{
int i;
char buf[22];
-
+
if (edi->bNeedDoEdsam)
{
if (step == -1)
{
fprintf(ed->edo,"Step %s, ED #%d ", gmx_step_str(step, buf), nr_edi);
fprintf(ed->edo," RMSD %f nm\n",rmsd);
- if (ed->eEDtype == eEDflood)
- fprintf(ed->edo, " Efl=%f deltaF=%f Vfl=%f\n",edi->flood.Efl,edi->flood.deltaF,edi->flood.Vfl);
}
if (edi->vecs.mon.neig)
/* Returns if any constraints are switched on */
static int ed_constraints(gmx_bool edtype, t_edpar *edi)
-{
- if (edtype == eEDedsam || edtype == eEDflood)
+{
+ if (edtype == eEDedsam || edtype == eEDflood)
{
- return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
- edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
+ return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
+ edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
edi->vecs.radcon.neig);
- }
+ }
return 0;
}
rvec *xstart = NULL; /* the positions which are subject to ED sampling */
rvec fit_transvec; /* translation ... */
matrix fit_rotmat; /* ... and rotation from fit to reference structure */
-
+
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");
-
+
/* Needed for initializing radacc radius in do_edsam */
ed->bFirst = 1;
/* Read the whole edi file at once: */
read_edi_file(ed,ed->edpar,mtop->natoms,cr);
- /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
+ /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset 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
+ * 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,ir,cr,ed,edi);
-
+
/* Init flooding parameters if needed */
init_flood(edi,ed,ir->delta_t,cr);
}
}
- /* The master does the work here. The other nodes get the positions
+ /* 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))
{
/* 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 <= numedis; nr_edi++)
{
- /* We use srenew to allocate memory since the size of the buffers
+ /* We use srenew to allocate memory since the size of the buffers
* is likely to change with every ED dataset */
srenew(xfit , edi->sref.nr );
srenew(xstart, edi->sav.nr );
for (i=0; i < edi->sref.nr; i++)
{
copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
-
- /* Save the sref positions such that in the next time step the molecule can
+
+ /* Save the sref positions such that in the next time step the molecule can
* be made whole again (in the parallel case) */
if (PAR(cr))
copy_rvec(xfit[i], edi->sref.x_old[i]);
}
-
+
/* Extract the positions of the atoms subject to ED sampling */
for (i=0; i < edi->sav.nr; i++)
{
copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
- /* Save the sav positions such that in the next time step the molecule can
+ /* Save the sav positions such that in the next time step the molecule can
* be made whole again (in the parallel case) */
if (PAR(cr))
copy_rvec(xstart[i], edi->sav.x_old[i]);
}
-
+
/* 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);
+ translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm (dataset #%d)\n",
rmsd_from_structure(xfit, &edi->sref), nr_edi);
rad_project(edi, xstart, &edi->vecs.radcon, cr);
/* process structure that will serve as origin of expansion circle */
- if (eEDflood == ed->eEDtype)
+ 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)
{
translate_and_rotate(edi->sori.x, edi->sav.nr, fit_transvec, fit_rotmat);
rad_project(edi, edi->sori.x, &edi->vecs.radacc, cr);
rad_project(edi, edi->sori.x, &edi->vecs.radfix, cr);
- if (ed->eEDtype == eEDflood)
+ 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, xstart, &edi->vecs.radacc, cr);
rad_project(edi, xstart, &edi->vecs.radfix, cr);
- if (ed->eEDtype == eEDflood)
+ if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
{
if (edi->flood.bHarmonic)
{
}
}
/* For convenience, output the center of the flooding potential for the eigenvectors */
- if (eEDflood == ed->eEDtype)
+ 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 (adding %11.4e/timestep)\n",
- i, edi->flood.vecs.refproj[i], edi->flood.vecs.refprojslope[i]);
+ fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", i, edi->flood.vecs.refproj[i]);
+ if (edi->flood.bHarmonic)
+ fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
+ fprintf(stdout, "\n");
}
}
write_edo(nr_edi, edi, ed, -1, 0);
/* Prepare for the next edi data set: */
- edi=edi->next_edi;
+ edi=edi->next_edi;
}
/* Cleaning up on the master node: */
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 */
}
else
{
- /* In the single-CPU case, point the local atom numbers pointers to the global
+ /* 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 <= numedis; nr_edi++)
/* 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);
+ 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->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);
/* An on we go to the next edi dataset */
edi=edi->next_edi;
}
-
- /* Flush the edo file so that the user can check some things
+
+ /* Flush the edo file so that the user can check some things
* when the simulation has started */
if (ed->edo)
fflush(ed->edo);
* 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 datasets (usually one) */
- edi = ed->edpar;
+ edi = ed->edpar;
edinr = 0;
while (edi != NULL)
{
/* initialise radacc radius for slope criterion */
buf->oldrad=calc_radius(&edi->vecs.radacc);
- /* Copy the positions into buf->xc* arrays and after ED
+ /* 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, buf->bUpdateShifts, xs,
- edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
+
+ communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, xs,
+ edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
#ifdef DEBUG_ED
dump_xcoll(edi, buf, cr, step);
#endif
/* 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, buf->bUpdateShifts, xs,
+ communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, 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 get_positions.
* We do not need to uptdate the shifts until the next NS step */
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_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);
{
if (edi->bRefEqAv)
{
- /* Indices of reference and average structures are identical,
+ /* 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);
}
if (MASTER(cr) && !bSuppress)
write_edo(edinr, edi, ed, step, rmsdev);
}
-
+
/* Copy back the positions unless monitoring only */
if (ed_constraints(ed->eEDtype, edi))
{
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
+ /* 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]],
+ 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: */
}
}
} /* END of if (edi->bNeedDoEdsam) */
-
+
/* Prepare for the next ED dataset */
edi = edi->next_edi;
-
+
} /* END of loop over ED datasets */
ed->bFirst = FALSE;
* G R O M A C S
*
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
*
* The make_edi program was generously contributed by Oliver Lange, based
* 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
#include "rmpbc.h"
#include "txtdump.h"
#include "eigio.h"
-#include "index.h"
+#include "index.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 kT;
real constEfl;
- real alpha2;
+ real alpha2;
} t_edflood;
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
+ 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 */
/*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]=',';
pos[n+1]='1';
pos[n+2]='\0';
*list = NULL;
-
+
while ((c=*pos)!=0) {
switch(status) {
/* expect a number */
else status=sError; break;
/* have read a '-' -> expect a number */
- case sMinus:
+ case sMinus:
if (isdigit(c)) {
end=pos;
status=sRange; break;
} else status=sError; break;
-
+
case sSteppedRange:
if (isdigit(c)) {
if (step) {
- status=sError; break;
+ 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==',') {
/* 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));
+ gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));
/* logical error occured */
case sZero:
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
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],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;
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]);
- }
+ 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, evEND };
+enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON, evNr };
#define oldMAGIC 666
-#define MAGIC 669
+#define MAGIC 670
-void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
- int nvec, int *eig_listen[], real* evStepList[])
+void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
+ int nvec, int *eig_listen[], real* evStepList[])
{
/* write edi-file */
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",
- edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic);
-
+ 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]);
/*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[])
+int read_conffile(const char *confin,char *title,rvec *x[])
{
/* read coordinates out of STX file */
int natoms;
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)
+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++) {
+ 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;
}
tcap += d;
strcat(f0,"%*s");
}
- }
+ }
return vec_params;
-}
+}
void init_edx(struct edix *edx) {
snew(edx->anrs,1);
}
-void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
- atom_id igro[],rvec *x,const char* structure)
+void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
+ atom_id igro[],rvec *x,const char* structure)
{
/* filter2edx copies coordinates from x to edx which are given in index
*/
-
+
int pos,i;
int ix=edx->nr;
edx->nr+=nindex;
void get_structure(t_atoms *atoms,const char *IndexFile,
const char *StructureFile,struct edix *edx,int nfit,
- atom_id ifit[],int natoms, atom_id index[])
+ atom_id ifit[],int natoms, atom_id index[])
{
atom_id *igro; /*index corresponding to target or origin structure*/
int ngro;
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",
static const char *desc[] = {
"[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
- "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
+ "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
"normal modes anaysis ([TT]g_nmeig[tt]).",
"ED sampling can be used to manipulate the position along collective coordinates",
"(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
"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.",
"Note on the parallel implementation: since ED sampling is a 'global' thing",
"(collective coordinates etc.), at least on the 'protein' side, ED sampling",
"is not very parallel-friendly from an implentation point of view. Because",
- "parallel ED requires much extra communication, expect the performance to be",
+ "parallel ED requires some extra communication, expect the performance to be",
"lower as in a free MD simulation, especially on a large number of nodes. [PAR]",
"All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a .edo file. In the output",
"file, per OUTFRQ step the following information is present: [PAR]",
"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 { evSTEPEND = evRADFIX + 1};
- static const char* evSelections[evEND]= {NULL,NULL,NULL,NULL,NULL,NULL};
- static const char* evOptions[evEND] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
- static const char* evParams[evSTEPEND] ={NULL,NULL};
- static const char* evStepOptions[evSTEPEND] = {"-linstep","-accdir","-not_used","-radstep"};
- static real* evStepList[evSTEPEND];
+ static t_edipar edi_params;
+
+ enum { evStepNr = evRADFIX + 1};
+ static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
+ static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
+ static const char* evParams[evStepNr] ={NULL,NULL};
+ static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
+ static const char* ConstForceStr;
+ static real* evStepList[evStepNr];
static real radfix=0.0;
static real deltaF0=150;
static real deltaF=0;
static real constEfl=0.0;
static real alpha=1;
static int eqSteps=0;
- static int* listen[evEND];
+ 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;
"Indices of eigenvectors for fixed increment linear sampling" },
{ "-linacc", FALSE, etSTR, {&evSelections[1]},
"Indices of eigenvectors for acceptance linear sampling" },
- { "-flood", FALSE, etSTR, {&evSelections[2]},
- "Indices of eigenvectors for flooding"},
{ "-radfix", FALSE, etSTR, {&evSelections[3]},
"Indices of eigenvectors for fixed increment radius expansion" },
{ "-radacc", FALSE, etSTR, {&evSelections[4]},
"Indices of eigenvectors for acceptance radius expansion" },
{ "-radcon", FALSE, etSTR, {&evSelections[5]},
"Indices of eigenvectors for acceptance radius contraction" },
+ { "-flood", FALSE, etSTR, {&evSelections[2]},
+ "Indices of eigenvectors for flooding"},
{ "-outfrq", FALSE, etINT, {&edi_params.outfrq},
"Freqency (in steps) of writing output in [TT].edo[tt] file" },
{ "-slope", FALSE, etREAL, { &edi_params.slope},
"Minimal slope in acceptance radius expansion"},
+ { "-linstep", FALSE, etSTR, {&evParams[0]},
+ "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
+ { "-accdir", FALSE, etSTR, {&evParams[1]},
+ "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
+ { "-radstep", FALSE, etREAL, {&radfix},
+ "Stepsize (nm/step) for fixed increment radius expansion"},
{ "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
- "Max nr of steps per cycle" },
+ "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 - used for flooding"},
+ "Target destabilization energy for flooding"},
{ "-deltaF", FALSE, etREAL, {&deltaF},
- "Start deltaF with this parameter - default 0, i.e. nonzero values only needed for restart"},
- { "-tau", FALSE, etREAL, {&tau},
+ "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"},
- { "-eqsteps", FALSE, etINT, {&eqSteps},
- "Number of steps to run without any perturbations "},
{ "-Eflnull", FALSE, etREAL, {&constEfl},
- "This is the starting value of the flooding strength. The flooding strength is updated according to the adaptive flooding scheme. To use a constant flooding strength use [TT]-tau[tt] 0. "},
+ "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 "},
- { "-linstep", FALSE, etSTR, {&evParams[0]},
- "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
- { "-accdir", FALSE, etSTR, {&evParams[1]},
- "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
- { "-radstep", FALSE, etREAL, {&radfix},
- "Stepsize (nm/step) for fixed increment radius expansion"},
{ "-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},
+ { "-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;
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;
matrix topbox;
rvec *xtop;
gmx_bool bTop, bFit1;
-
+
t_filenm fnm[] = {
{ efTRN, "-f", "eigenvec", ffREAD },
- { efXVG, "-eig", "eigenval", ffOPTRD },
+ { efXVG, "-eig", "eigenval", ffOPTRD },
{ efTPS, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
{ efSTX, "-tar", "target", ffOPTRD},
CopyRight(stderr,argv[0]);
parse_common_args(&argc,argv, 0 ,
NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
-
+
indexfile=ftp2fn_null(efNDX,NFILE,fnm);
EdiFile=ftp2fn(efEDI,NFILE,fnm);
TargetFile = opt2fn_null("-tar",NFILE,fnm);
OriginFile = opt2fn_null("-ori",NFILE,fnm);
-
-
- for (ev_class=0; ev_class<evEND; ++ev_class) {
+
+
+ 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<evSTEPEND-2) {
+ if (ev_class<evStepNr-2) {
/*if apropriate get list of stepsizes for these eigenvectors*/
- if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
+ 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++)
+ for (i=0; i<nvecs; i++)
evStepList[ev_class][i]=0.0;
}
} else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
snew(evStepList[ev_class],nvecs);
for (i=0; i<nvecs; i++)
evStepList[ev_class][i]=radfix;
-
} else if (ev_class == evFLOOD) {
snew(evStepList[ev_class],nvecs);
+
+ /* Are we doing constant force flooding? In that case, we read in
+ * the fproj values from the command line */
+ if (opt2parg_bSet("-constF", NPA, pa))
+ {
+ evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
+ }
} else {}; /*to avoid ambiguity */
} else { /* if there are no eigenvectors for this option set list to zero */
listen[ev_class]=NULL;
snew(listen[ev_class],1);
listen[ev_class][0]=0;
- };
- };
-
+ }
+ }
+
/* print the interpreted list of eigenvectors - to give some feedback*/
- for (ev_class=0; ev_class<evEND; ++ev_class) {
+ for (ev_class=0; ev_class<evNr; ++ev_class) {
printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
i=0;
while(listen[ev_class][i])
printf("%d ",listen[ev_class][i++]);
printf("\n");
}
-
- EigvecFile=NULL;
- EigvecFile=opt2fn("-f",NFILE,fnm);
-
-
+
+ EigvecFile=NULL;
+ EigvecFile=opt2fn("-f",NFILE,fnm);
+
/*read eigenvectors from eigvec.trr*/
read_eigenvectors(EigvecFile,&natoms,&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",natoms);
get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
if (i!=natoms) {
i,natoms);
}
printf("\n");
-
-
+
+
if (xref1==NULL)
{
if (bFit1)
ifit=index;
}
- /* if -flood 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);
-
- /*number of eigenvectors*/
- edi_params.ned=natoms;
- 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);
- }
-
+ 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=natoms;
+
/*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,natoms,xav1,index);
-
+
/* Store target positions in edi_params */
if (opt2bSet("-tar",NFILE,fnm))
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff