From 89ccf9a78a1091ed34b439adcd53bbf88e3b4d01 Mon Sep 17 00:00:00 2001 From: Gareth Tribello <gareth.tribello@gmail.com> Date: Fri, 31 Aug 2012 17:57:49 +0200 Subject: [PATCH] Many improvements to the user manual --- regtest/rt19-mpi/plumed.dat | 12 +-- regtest/rt19/plumed.dat | 12 +-- regtest/rt22/plumed.dat | 8 +- src/ColvarContactMap.cpp | 6 +- src/ColvarCoordination.cpp | 14 +-- src/ColvarDRMSD.cpp | 26 +++++- src/ColvarRMSD.cpp | 40 +++++++- src/DistributionCVdens.cpp | 39 ++++++-- src/DistributionLessThan.cpp | 6 +- src/DistributionMoreThan.cpp | 7 +- src/DistributionWithin.cpp | 9 +- src/FunctionTarget.cpp | 6 +- src/HistogramBead.cpp | 109 ++++++++++++++++++---- src/HistogramBead.h | 12 +-- src/MultiColvarCoordination.cpp | 24 +++-- src/MultiColvarDensity.cpp | 5 +- src/MultiColvarDistance.cpp | 34 ++++--- src/MultiColvarSecondaryStructureRMSD.cpp | 4 +- src/PDB.cpp | 9 -- src/PDB.h | 2 - src/SetupMolInfo.cpp | 5 +- src/SwitchingFunction.cpp | 65 ++++++++++--- src/SwitchingFunction.h | 1 - user-doc/Colvar.txt | 13 +-- user-doc/extract | 15 ++- 25 files changed, 338 insertions(+), 145 deletions(-) diff --git a/regtest/rt19-mpi/plumed.dat b/regtest/rt19-mpi/plumed.dat index fa64e0c13..e58f9306a 100644 --- a/regtest/rt19-mpi/plumed.dat +++ b/regtest/rt19-mpi/plumed.dat @@ -1,9 +1,9 @@ -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={NBINS=3 LOWER=0.5 UPPER=2.0} LABEL=d1 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={NBINS=3 LOWER=0.5 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d1n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN NBINS=3 LOWER=0.5 UPPER=2.0} LABEL=d1 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN NBINS=3 LOWER=0.5 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d1n DUMPDERIVATIVES ARG=d1.*,d1n.* FILE=derivatives1 FMT=%8.5f -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={LOWER=0.5 UPPER=1.5} LABEL=d2 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={LOWER=0.5 UPPER=1.5} NUMERICAL_DERIVATIVES LABEL=d2n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN LOWER=0.5 UPPER=1.5} LABEL=d2 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN LOWER=0.5 UPPER=1.5} NUMERICAL_DERIVATIVES LABEL=d2n DUMPDERIVATIVES ARG=d2.*,d2n.* FILE=derivatives2 FMT=%8.5f -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={LOWER=0.5 UPPER=1.0} WITHIN2={LOWER=1.0 UPPER=2.0} LABEL=d3 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={LOWER=0.5 UPPER=1.0} WITHIN2={LOWER=1.0 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d3n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={GAUSSIAN LOWER=0.5 UPPER=1.0} WITHIN2={GAUSSIAN LOWER=1.0 UPPER=2.0} LABEL=d3 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={GAUSSIAN LOWER=0.5 UPPER=1.0} WITHIN2={GAUSSIAN LOWER=1.0 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d3n DUMPDERIVATIVES ARG=d3.*,d3n.* FILE=derivatives3 FMT=%8.5f diff --git a/regtest/rt19/plumed.dat b/regtest/rt19/plumed.dat index fa64e0c13..e58f9306a 100644 --- a/regtest/rt19/plumed.dat +++ b/regtest/rt19/plumed.dat @@ -1,9 +1,9 @@ -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={NBINS=3 LOWER=0.5 UPPER=2.0} LABEL=d1 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={NBINS=3 LOWER=0.5 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d1n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN NBINS=3 LOWER=0.5 UPPER=2.0} LABEL=d1 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN NBINS=3 LOWER=0.5 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d1n DUMPDERIVATIVES ARG=d1.*,d1n.* FILE=derivatives1 FMT=%8.5f -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={LOWER=0.5 UPPER=1.5} LABEL=d2 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={LOWER=0.5 UPPER=1.5} NUMERICAL_DERIVATIVES LABEL=d2n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN LOWER=0.5 UPPER=1.5} LABEL=d2 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN={GAUSSIAN LOWER=0.5 UPPER=1.5} NUMERICAL_DERIVATIVES LABEL=d2n DUMPDERIVATIVES ARG=d2.*,d2n.* FILE=derivatives2 FMT=%8.5f -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={LOWER=0.5 UPPER=1.0} WITHIN2={LOWER=1.0 UPPER=2.0} LABEL=d3 -DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={LOWER=0.5 UPPER=1.0} WITHIN2={LOWER=1.0 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d3n +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={GAUSSIAN LOWER=0.5 UPPER=1.0} WITHIN2={GAUSSIAN LOWER=1.0 UPPER=2.0} LABEL=d3 +DISTANCES ATOMS1=1,2 ATOMS2=5,4 WITHIN1={GAUSSIAN LOWER=0.5 UPPER=1.0} WITHIN2={GAUSSIAN LOWER=1.0 UPPER=2.0} NUMERICAL_DERIVATIVES LABEL=d3n DUMPDERIVATIVES ARG=d3.*,d3n.* FILE=derivatives3 FMT=%8.5f diff --git a/regtest/rt22/plumed.dat b/regtest/rt22/plumed.dat index 5a534421f..6fe7298a3 100644 --- a/regtest/rt22/plumed.dat +++ b/regtest/rt22/plumed.dat @@ -1,9 +1,9 @@ -DENSITY SPECIES=1-100 SUBCELL={XLOWER=0.0 XUPPER=0.5 XSMEAR=0.6} LABEL=d1 -DENSITY SPECIES=1-100 SUBCELL={XLOWER=0.0 XUPPER=0.5 XSMEAR=0.6} NUMERICAL_DERIVATIVES LABEL=d1n +DENSITY SPECIES=1-100 SUBCELL={GAUSSIAN XLOWER=0.0 XUPPER=0.5 XSMEAR=0.6} LABEL=d1 +DENSITY SPECIES=1-100 SUBCELL={GAUSSIAN XLOWER=0.0 XUPPER=0.5 XSMEAR=0.6} NUMERICAL_DERIVATIVES LABEL=d1n PRINT ARG=d1.* FILE=colvar1 FMT=%8.4f DUMPDERIVATIVES ARG=d1.*,d1n.* FILE=derivatives1 FMT=%8.4f -COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0} SUBCELL={YLOWER=0.0 YUPPER=0.5 YSMEAR=0.1} LABEL=c1 -COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0} SUBCELL={YLOWER=0.0 YUPPER=0.5 YSMEAR=0.1} NUMERICAL_DERIVATIVES LABEL=c1n +COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0} SUBCELL={GAUSSIAN YLOWER=0.0 YUPPER=0.5 YSMEAR=0.1} LABEL=c1 +COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0} SUBCELL={GAUSSIAN YLOWER=0.0 YUPPER=0.5 YSMEAR=0.1} NUMERICAL_DERIVATIVES LABEL=c1n PRINT ARG=c1.* FILE=colvar2 FMT=%8.4f DUMPDERIVATIVES ARG=c1.*,c1n.* FILE=derivatives2 FMT=%8.4f diff --git a/src/ColvarContactMap.cpp b/src/ColvarContactMap.cpp index 3c53baa09..9a9a0e4ce 100644 --- a/src/ColvarContactMap.cpp +++ b/src/ColvarContactMap.cpp @@ -74,7 +74,11 @@ void ColvarContactMap::registerKeywords( Keywords& keys ){ "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one contact will be " "calculated for each ATOM keyword you specify."); keys.reset_style("ATOMS","atoms"); - keys.add("numbered","SWITCH","The switching functions to use for each of the contacts in your map. You can either specify a global switching function using SWITCH or one switching function for each contact"); keys.reset_style("SWITCH","compulsory"); + keys.add("numbered","SWITCH","The switching functions to use for each of the contacts in your map. " + "You can either specify a global switching function using SWITCH or one " + "switching function for each contact. Details of the various switching " + "functions you can use are provided on \\ref switchingfunction."); + keys.reset_style("SWITCH","compulsory"); keys.addFlag("SUM",false,"calculate the sum of all the contacts in the input"); } diff --git a/src/ColvarCoordination.cpp b/src/ColvarCoordination.cpp index f11e0405b..364c6adf6 100644 --- a/src/ColvarCoordination.cpp +++ b/src/ColvarCoordination.cpp @@ -49,7 +49,7 @@ relevant subset of the pairwise distance are calculated at every step. The following example instructs plumed to the average coordination number of the atoms in group 1-10 with the atoms in group 20-100. These coordination numbers count the number of atoms within the group that are within 0.3 nm of the central atom. A neighbour list is used to make this calculation faster, this neighbour list is updated every 100 steps. \verbatim -COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NL_CUTOFF=0.5 UPDATE=100 +COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100 \endverbatim */ @@ -81,11 +81,13 @@ void ColvarCoordination::registerKeywords( Keywords& keys ){ keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation"); keys.add("atoms","GROUPA","The list of central atoms for which we are calculating our coordination numbers"); keys.add("atoms","GROUPB","The list of neighbourhood atoms for which we are using to calculate coordination numbers"); - keys.add("optional","NN","The n parameter of the switching function "); - keys.add("optional","MM","The m parameter of the switching function "); - keys.add("optional","D_0","The d_0 parameter of the switching function"); - keys.add("optional","R_0","The r_0 parameter of the switching function"); - keys.add("optional","SWITCH","A generic switching function"); + keys.add("compulsory","NN","6","The n parameter of the switching function "); + keys.add("compulsory","MM","12","The m parameter of the switching function "); + keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function"); + keys.add("compulsory","R_0","The r_0 parameter of the switching function"); + keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. " + "The following provides information on the \\ref switchingfunction that are available. " + "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords."); keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list"); keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list"); } diff --git a/src/ColvarDRMSD.cpp b/src/ColvarDRMSD.cpp index 51bfd837a..f9272bf1d 100644 --- a/src/ColvarDRMSD.cpp +++ b/src/ColvarDRMSD.cpp @@ -34,9 +34,27 @@ namespace PLMD{ /* Calculate the distance RMSD with respect to a reference structure. - -Only pairs of atoms whose distance in the reference structure is within -LOWER_CUTOFF and UPPER_CUTOFF are considered. +To calculate the root-mean-square deviation between the atoms in two configurations +you must first superimpose the two structures in some ways. Obviously, it is the internal vibrational +motions of the structure - i.e. not the translations and rotations - that are interesting. However, +aligning two structures by removing the translational and rotational motions is not easy. Furthermore, +in some cases there can be alignment issues caused by so-called frame-fitting problems. It is thus +often cheaper and easier to calculate the distances between all the pairs of atoms. The distance +between the two structures, \f$\mathbf{X}^a\f$ and \f$\mathbf{X}^b\f$ can then be measured as: + +\f[ +d(\mathbf{X}^A, \mathbf{X}^B) = \frac{1}{N(N-1)} \sum_{i \ne j} [ d(\mathbf{x}_i^a,\mathbf{x}_j^a) - d(\mathbf{x}_i^b,\mathbf{x}_j^b) ]^2 +\f] + +where \f$N\f$ is the number of atoms and \f$d(\mathbf{x}_i,\mathbf{x}_j)\f$ represents the distance between +atoms \f$i\f$ and \f$j\f$. Clearly, this representation of the configuration is invariant to translation and rotation. +However, it can become expensive to calculate when the number of atoms is large. This can be resolved +within the DRMSD colvar by setting LOWER_CUTOFF and UPPER_CUTOFF. These keywords ensure that only +pairs of atoms that are within a certain range are incorporated into the above sum. + +In PDB files the atomic coordinates and box lengths should be in Angstroms unless +you are working with natural units. If you are working with natural units then the coordinates +should be in your natural length unit. For more details on the PDB file format visit http://www.wwpdb.org/docs.html \par Examples @@ -71,7 +89,7 @@ PLUMED_REGISTER_ACTION(ColvarDRMSD,"DRMSD") void ColvarDRMSD::registerKeywords(Keywords& keys){ Colvar::registerKeywords(keys); - keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV. " + PDB::documentation() ); + keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV."); keys.add("compulsory","LOWER_CUTOFF","only pairs of atoms further than LOWER_CUTOFF are considered in the calculation."); keys.add("compulsory","UPPER_CUTOFF","only pairs of atoms closer than UPPER_CUTOFF are considered in the calculation."); } diff --git a/src/ColvarRMSD.cpp b/src/ColvarRMSD.cpp index ccd926cbb..d912f81ae 100644 --- a/src/ColvarRMSD.cpp +++ b/src/ColvarRMSD.cpp @@ -56,10 +56,40 @@ namespace PLMD{ /* Calculate the RMSD with respect to a reference structure. -To perform an ??optimal?? (what does this mean algorithmical speed wise?) alignment -using the Kearsley algorithm then use TYPE=OPTIMAL. Otherwise -use TYPE=SIMPLE, which will not perform optimal alignment and will only -remove the translation of the center of mass. +To calculate the root-mean-square deviation between the atoms in two configurations +you must first superimpose the two structures in some ways. Obviously, it is the internal vibrational +motions of the structure - i.e. not the translations and rotations - that are interesting. It is +possible to align two structures (i.e. remove the translational and rotational motions of the +system) using the Kearsley \cite kearsley algorithm. This algorithm calculates a +roto-translation matrix from a set of \e alignment atoms \f${\cal A}\f$. The amount by which +a \e displacement set of atoms, \f${\cal B}\f$, has been moved can then be calculated using: + +\f[ +d(X_j,X_i) = \sum_{a=1}^{N_{\cal B}} (X^{(j)}_a - M_{ij} X^{(i)}_a)^2 +\f] + +where \f$M_{ij}\f$ is the roto-translation matrix translated by the Kearsley algorithm. + +When a pdb input file is required to specify the structure to align to in an RMSD calculation +the occupancy and beta columns are used to specify which atoms form part of the alignment set +and which atoms form part of the displacement set. Values of 1.0 and 0.0 indicate that the +atom is to be used in the \e alignment set only, while values of 0.0 and 1.0 indicates that +the atom is to be used in the \e displacement set only. Values of 1.0 and 1.0 indicate that +the atom is to be used in both the \e alignment and \e displacement sets. Users can also +use fractional values for beta and the occupancy values. We recommend you only do this when +you really know what you are doing however as the results can be rather strange. +When this form of RMSD is used to calculate the secondary structure variables (\ref ALPHARMSD, +\ref ANTIBETARMSD and \ref PARABETARMSD) all the atoms in the segment are assumed to be part of +both the \e alignment and \e displacement sets. + +The \e OPTIMAL type of rmsd calculation removes both translational and rotational motions. +There may, however, be ocasions where you do not need to remove the rotational motions. That is +to say there may be times when you only need to align the centers of mass of the two structures. +This can be done by performing the \e SIMPLE form of RMSD calculation. + +In PDB files the atomic coordinates and box lengths should be in Angstroms unless +you are working with natural units. If you are working with natural units then the coordinates +should be in your natural length unit. For more details on the PDB file format visit http://www.wwpdb.org/docs.html \attention The documentation here needs some work as it is not very clear to me @@ -84,7 +114,7 @@ PLUMED_REGISTER_ACTION(ColvarRMSD,"RMSD") void ColvarRMSD::registerKeywords(Keywords& keys){ Colvar::registerKeywords(keys); - keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV. " + PDB::documentation() ); + keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV."); keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE."); keys.addFlag("SQUARED",false," This should be setted if you want MSD instead of RMSD "); } diff --git a/src/DistributionCVdens.cpp b/src/DistributionCVdens.cpp index 7925b9d4e..fcb912f24 100644 --- a/src/DistributionCVdens.cpp +++ b/src/DistributionCVdens.cpp @@ -26,6 +26,37 @@ namespace PLMD { +//+PLUMEDOC INTERNAL subcell +/* +Imagine we have a collection of collective variables that can all be assigned to a particular point in three +dimensional space. For example, we could have the values of the coordination numbers for all the atoms in the +system. Because each CV value can be assigned to a particular point in space we can calculate the average +value of the cv in a particular part of the box using: + +\f[ +\overline{s}_{\tau} = \frac{ \sum_i s_i w(x_i)w(y_i)w(z_i) }{ \sum_i w(x_i)w(y_i)w(z_i) } +\f] + +where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (x_i,y_i,z_i)\f$. +Because we want to calculate the average cv in the subregion \f$\tau\f$ we introduce the three \f$w()\f$ to tell us +whether or not \f$i\f$ is inside \f$\tau\f$. These functions are given by: + +\f[ +w(x_i) = \int_{a_x}^{b_x} \textrm{d}x K\left( \frac{x - x_i}{w_x} \right) +\f] + +where \f$K\f$ is one of the kernel functions described on \ref histogrambead. + +All the input to calculate these quantities is provided through a single keyword that will have the following form + +KEYWORD={TYPE XLOWER=\f$a_x\f$ XUPPER=\f$b_x\f$ XSMEAR=\f$\frac{w_x}{b_x-a_x}\f$ YLOWER=\f$a_y\f$ YUPPER=\f$b_y\f$ YSMEAR=\f$\frac{w_y}{b_y-a_y}\f$ ZLOWER=\f$a_z\f$ ZUPPER=\f$b_z\f$ ZSMEAR=\f$\frac{w_z}{b_z-a_z}\f$} + +The \f$a_x\f$, \f$b_x\f$, \f$a_y\f$ ... parameters are specified as a fraction of the box length. In directions for which upper and lower bounds are not +specified the \f$a\f$ and \f$b\f$ values are assumed equal to 0 and 1. If SMEAR values are not specified \f$\frac{w_x}{b_x-a_x}=0.5\f$. +*/ +//+ENDPLUMEDOC + + class cvdens : public NormedSumVessel { private: bool isDensity; @@ -47,12 +78,8 @@ public: PLUMED_REGISTER_VESSEL(cvdens,"SUBCELL") void cvdens::reserveKeyword( Keywords& keys ){ - keys.reserve("numbered","SUBCELL","calculate the average value of the CV within a portion of the box and store it in a value called subcell. " - "To make this quantity continuous it is calculated using \\f$ a = \\frac{\\sum_{i=1}^N s_i w(x_i)w(y_i)w(z_i)}{\\sum_{i=1}^N w(x_i)w(y_i)w(z_i)}\\f$ " - "where the sum is over the collective variables \\f$ s_i \\f$ which can be thought to be at \\f$ (x_i,y_i,z_i)\\f$. The three \\f$w(x)\\f$ functions " - "specify the extent of the subregion in which we are calculating the average CV in each direction. To make these assignment functions continuous they " - "are calculated using " + HistogramBead::documentation(true) + ". If no parameters for the \\f$ w(d) \\f$ function in a particular direction is specified " - "then the subcell is assumed to incorporate the entirity of the box in that direction."); + keys.reserve("numbered","SUBCELL","calculate the average value of the CV within a portion of the box and store it in a value called subcell. " + "For more details on how this quantity is calculated see \\ref subcell."); } cvdens::cvdens( const VesselOptions& da ) : diff --git a/src/DistributionLessThan.cpp b/src/DistributionLessThan.cpp index 5a27d9045..3da18f7b3 100644 --- a/src/DistributionLessThan.cpp +++ b/src/DistributionLessThan.cpp @@ -39,8 +39,10 @@ public: PLUMED_REGISTER_VESSEL(less_than,"LESS_THAN") void less_than::reserveKeyword( Keywords& keys ){ - keys.reserve("optional","LESS_THAN", "take the number of variables less than the specified target and " - "store it in a value called lt<target>. " + SwitchingFunction::documentation() ); + keys.reserve("optional","LESS_THAN","calculate the number of variables less than a certain target value. " + "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ " + "is a \\ref switchingfunction. The final value can be referenced using " + "\\e label.lt\\f$r_0\\f$."); } less_than::less_than( const VesselOptions& da ) : diff --git a/src/DistributionMoreThan.cpp b/src/DistributionMoreThan.cpp index 258b3fde6..a5c5ee955 100644 --- a/src/DistributionMoreThan.cpp +++ b/src/DistributionMoreThan.cpp @@ -39,9 +39,10 @@ public: PLUMED_REGISTER_VESSEL(more_than,"MORE_THAN") void more_than::reserveKeyword( Keywords& keys ){ - keys.reserve("numbered","MORE_THAN","take the number of variables more than the specified target and store it in a value called gt<target>. " - "This is calculated using \\f$1.0 - s(r)\\f$, where \\f$s(r)\\f$ is a switching function. " + - SwitchingFunction::documentation() ); + keys.reserve("optional","MORE_THAN","calculate the number of variables more than a certain target value. " + "This quantity is calculated using \\f$\\sum_i 1.0 - \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ " + "is a \\ref switchingfunction. The final value can be referenced using " + "\\e label.gt\\f$r_0\\f$."); } more_than::more_than( const VesselOptions& da ) : diff --git a/src/DistributionWithin.cpp b/src/DistributionWithin.cpp index 8b8f016e7..0855442b3 100644 --- a/src/DistributionWithin.cpp +++ b/src/DistributionWithin.cpp @@ -19,11 +19,10 @@ public: PLUMED_REGISTER_VESSEL(within,"WITHIN") void within::reserveKeyword( Keywords& keys ){ - keys.reserve("numbered","WITHIN", "calculate the number of variables that are within a certain range and store it in a value called between<lowerbound>&<upperbound> " - "or create a discretized histogram of the distribution for a particular range. To make these quantities continuous they are " - "calculated using " + HistogramBead::documentation(false) + " If you add the NBINS keyword the range between your upper and " - "lower bounds is divided into a discrete number of bins. Adding the NORM flag will calculate the fraction of colvars in " - "range of interest rather than the total number"); + keys.reserve("numbered","WITHIN", "calculate the number of values that are within a certain range or create a discretized " + "histogram of the distribution of cvs within a particular rage by adding the NBINS " + "to the function specifier. These quantities are described using kernel density estimation as described on " + "\\ref histogrambead. The final values can be referenced using \\e label.between\\f$a\\f$&\\f$b\\f$."); } within::within( const VesselOptions& da ) : diff --git a/src/FunctionTarget.cpp b/src/FunctionTarget.cpp index 30a0a5e95..57cb62579 100644 --- a/src/FunctionTarget.cpp +++ b/src/FunctionTarget.cpp @@ -56,7 +56,11 @@ PLUMED_REGISTER_ACTION(FunctionTarget,"TARGET") void FunctionTarget::registerKeywords(Keywords& keys){ Function::registerKeywords(keys); keys.use("ARG"); - keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure. " + PDB::documentation() ); + keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure. In the PDB file the atomic " + "coordinates and box lengths should be in Angstroms unless you are working with natural units. " + "If you are working with natural units then the coordinates should be in your natural length unit. " + "The charges and masses of the atoms (if required) should be inserted in the beta and occupancy " + "columns respectively. For more details on the PDB file format visit http://www.wwpdb.org/docs.html"); keys.add("optional","REFERENCE_VEC","the vector of values for the CVs at the reference point (if you use this you don't need REFERENCE)"); } diff --git a/src/HistogramBead.cpp b/src/HistogramBead.cpp index 55de76715..588825238 100644 --- a/src/HistogramBead.cpp +++ b/src/HistogramBead.cpp @@ -27,21 +27,59 @@ using namespace PLMD; -std::string HistogramBead::documentation( bool dir ) { - std::ostringstream ostr; - if(dir){ - ostr<<"\\f$ w(x)=\\int_a^b \\frac{1}{\\sqrt{2\\pi}\\sigma_d} \\exp\\left( -\\frac{(x'-x)^2}{2\\sigma_x^2} \\right) \\textrm{d}x' \\f$"; - ostr<<"where \\f$ \\sigma_x=(b_x-a_x)k_x \\f$. The parameters of the functions are specifed in fractional coordinates using "; - ostr<<"{XLOWER=\\f$a_x\\f$ XUPPER=\\f$b_x\\f$ XSMEAR=\\f$k_x\\f$ YLOWER=\\f$a_y\\f$ YUPPER=\\f$b_y\\f$ YSMEAR=\\f$k_y\\f$ ZLOWER=\\f$a_z\\f$ ZUPPER=\\f$b_z\\f$ ZSMEAR=\\f$k_z\\f$}."; - ostr<<"If any of the SMEAR keywords are not present then the default \\f$k_x=0.5\\f$ is used in that direction. "; - } else { - ostr<<"\\f$ w(r)=\\int_a^b \\frac{1}{\\sqrt{2\\pi}\\sigma} \\exp\\left( -\\frac{(r - r')^2}{2\\sigma^2} \\right) \\textrm{d}r' \\f$"; - ostr<<"where \\f$ \\sigma=k\\frac{b-a}{n} \\f$ and \\f$n\\f$ is the number of bins you are dividing the range into. "; - ostr<<"The parameters of the function are specified using {LOWER=\\f$a\\f$ UPPER=\\f$b\\f$ SMEAR=\\f$k\\f$}. "; - ostr<<"If the SMEAR keyword is not present then by default \\f$k=0.5\\f$."; - } - return ostr.str(); -} +//+PLUMEDOC INTERNAL histogrambead +/* +If we have multiple instances of a variable we can estimate the probability distribution (pdf) +for that variable using a process called kernel density estimation: + +\f[ +P(s) = \sum_i K\left( \frac{s - s_i}{w} \right) +\f] + +In this equation \f$K\f$ is a symmetric funciton that must integrate to one that is often +called a kernel function and \f$w\f$ is a smearing parameter. From a pdf calculated using +kernel density estimation we can calculate the number/fraction of values between an upper and lower +bound using: + +\f[ +w(s) = \int_a^b \sum_i K\left( \frac{s - s_i}{w} \right) +\f] + +All the input to calculate a quantity like \f$w(s)\f$ is generally provided through a single +keyword that will have the following form: + +KEYWORD={TYPE UPPER=\f$a\f$ LOWER=\f$b\f$ SMEAR=\f$\frac{w}{b-a}\f$} + +This will calculate the number of values between \f$a\f$ and \f$b\f$. To calculate +the fraction of values you add the word NORM to the input specification. If the +function keyword SMEAR is not present \f$w\f$ is set equal to \f$0.5(b-a)\f$. Finally, +type should specify one of the kernel types that is present in plumed. These are listed +in the table below: + +<table align=center frame=void width=95%% cellpadding=5%%> +<tr> +<td> TYPE </td> <td> FUNCTION </td> +</tr> <tr> +<td> GAUSSIAN </td> <td> \f$\frac{1}{\sqrt{2\pi}w} \exp\left( -\frac{(s-s_i)^2}{2w^2} \right)\f$ </td> +<td> TRIANGULAR </td> <td> \f$ \frac{1}{2w} \left( 1. - \left| \frac{s-s_i}{w} \right| \right) \quad \frac{s-s_i}{w}<1 \f$ </td> +</tr> +</table> + +Some keywords can also be used to calculate a descretized version of the histogram. That +is to say the number of values between \f$a\f$ and \f$b\f$, the number of values between +\f$b\f$ and \f$c\f$ and so on. A keyword that specifies this sort of calculation would look +something like + +KEYWORD={TYPE UPPER=\f$a\f$ LOWER=\f$b\f$ NBINS=\f$n\f$ SMEAR=\f$\frac{w}{n(b-a)}\f$} + +This specification would calculate the following vector of quantities: + +\f[ +w_j(s) \int_{a + \frac{j-1}{n}(b-a)}^{a + \frac{j}{n}(b-a)} \sum_i K\left( \frac{s - s_i}{w} \right) +\f] + +*/ +//+ENDPLUMEDOC std::string HistogramBead::description() const { std::ostringstream ostr; @@ -54,6 +92,8 @@ void HistogramBead::generateBins( const std::string& params, const std::string& std::vector<std::string> data=Tools::getWords(params); plumed_assert(data.size()>=1); + std::string name=data[0]; + unsigned nbins; std::vector<double> range(2); std::string smear; bool found_nb=Tools::parse(data,dd+"NBINS",nbins); plumed_massert(found_nb,"Number of bins in histogram not found"); @@ -72,10 +112,9 @@ void HistogramBead::generateBins( const std::string& params, const std::string& for(unsigned i=0;i<nbins;++i){ Tools::convert( range[0]+i*delr, lb ); Tools::convert( range[0]+(i+1)*delr, ub ); - bins.push_back( dd + "LOWER=" + lb + " " + dd + "UPPER=" + ub + " " + dd + "SMEAR=" + smear + " " + normstr ); + bins.push_back( name + " " + dd + "LOWER=" + lb + " " + dd + "UPPER=" + ub + " " + dd + "SMEAR=" + smear + " " + normstr ); } plumed_assert(bins.size()==nbins); - if( dd.empty() ) plumed_massert(data.empty(),"Error reading histogram"); } void HistogramBead::set( const std::string& params, const std::string& dd, std::string& errormsg ){ @@ -83,6 +122,11 @@ void HistogramBead::set( const std::string& params, const std::string& dd, std:: std::vector<std::string> data=Tools::getWords(params); plumed_assert(data.size()>=1); + std::string name=data[0]; + if(name=="GAUSSIAN") type=gaussian; + else if(name=="TRIANGULAR") type=triangular; + else plumed_merror("cannot understand kernel type " + name ); + double smear; bool found_r=Tools::parse(data,dd+"LOWER",lowb); if( !found_r ) errormsg="Lower bound has not been specified use LOWER"; @@ -93,7 +137,6 @@ void HistogramBead::set( const std::string& params, const std::string& dd, std:: smear=0.5; bool found_b=Tools::parse(data,dd+"SMEAR",smear); width=smear*(highb-lowb); init=true; bool usenorm; bool found_n=Tools::parseFlag(data,dd+"NORM",usenorm); - if( dd.empty() ){ if( !data.empty()) errormsg="Error reading within"; } } void HistogramBead::set( double l, double h, double w){ @@ -108,3 +151,33 @@ void HistogramBead::printKeywords( Log& log ) const { hkeys.addFlag("NORM",false,"normalize the histogram according to the number of values we are histograming"); hkeys.print( log ); } + +double HistogramBead::calculate( double x, double& df ) const { + const double pi=3.141592653589793238462643383279502884197169399375105820974944592307; + assert(init && periodicity!=unset ); + double lowB, upperB, f; + if( type==gaussian ){ + lowB = difference( x, lowb ) / ( sqrt(2.0) * width ); + upperB = difference( x, highb ) / ( sqrt(2.0) * width ) ; + df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width ); + f = 0.5*( erf( upperB ) - erf( lowB ) ); + } else if( type==triangular ){ + lowB = ( difference( x, lowb ) / width ); + upperB = ( difference( x, highb ) / width ); + df=0; + if( fabs(lowB)<1. ) df = 1 - fabs(lowB) / width; + if( fabs(upperB)<1. ) df -= fabs(upperB) / width; + double ia, ib; + if (upperB<=-1. || lowB >=1.){ + f=0.; + } else { + if( lowB>-1.0 ){ ia=lowB; }else{ ia=-1.0; } + if( upperB<1.0 ){ ib=upperB; } else{ ib=1.0; } + f = (ib*(2.-fabs(ib))-ia*(2.-fabs(ia)))*0.5; + } + } else { + plumed_merror("function type does not exist"); + } + return f; +} + diff --git a/src/HistogramBead.h b/src/HistogramBead.h index 35daa4015..c50079343 100644 --- a/src/HistogramBead.h +++ b/src/HistogramBead.h @@ -44,11 +44,11 @@ private: double lowb; double highb; double width; + enum {gaussian,triangular} type; enum {unset,periodic,notperiodic} periodicity; double min, max, max_minus_min, inv_max_minus_min; double difference( const double& d1, const double& d2 ) const ; public: - static std::string documentation( bool dir ); static void generateBins( const std::string& params, const std::string& dd, std::vector<std::string>& bins ); HistogramBead(); std::string description() const ; @@ -104,16 +104,6 @@ double HistogramBead::difference( const double& d1, const double& d2 ) const { } else plumed_assert(0); return 0; } - -inline -double HistogramBead::calculate( double x, double& df ) const { - const double pi=3.141592653589793238462643383279502884197169399375105820974944592307; - assert(init && periodicity!=unset ); double lowB, upperB; - lowB = difference( x, lowb ) / ( sqrt(2.0) * width ); - upperB = difference( x, highb ) / ( sqrt(2.0) * width ) ; - df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width ); - return 0.5*( erf( upperB ) - erf( lowB ) ); -} } diff --git a/src/MultiColvarCoordination.cpp b/src/MultiColvarCoordination.cpp index 9ff837806..54f18cb47 100644 --- a/src/MultiColvarCoordination.cpp +++ b/src/MultiColvarCoordination.cpp @@ -34,22 +34,28 @@ namespace PLMD{ //+PLUMEDOC MCOLVAR COORDINATIONNUMBER /* -Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of +Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of coordination numbers such as the minimum, the number less than a certain quantity and so on. +To make the calculation of coordination numbers differentiable the following function is used: + +\f[ +s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m } +\f] + \par Examples The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves. The minimum coordination number is then calculated. \verbatim -COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN=0.1 +COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1} \endverbatim The following input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The number of coordination numbers more than 6 is then computed. \verbatim -COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN=6.0 +COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0} \endverbatim */ @@ -77,11 +83,13 @@ void MultiColvarCoordination::registerKeywords( Keywords& keys ){ MultiColvar::registerKeywords( keys ); ActionWithDistribution::autoParallelize( keys ); keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB"); - keys.add("optional","SWITCH","A switching function that defines which atoms are in the first coordination sphere. " + SwitchingFunction::documentation() ); - keys.add("optional","NN","The n parameter of the switching function "); - keys.add("optional","MM","The m parameter of the switching function "); - keys.add("optional","D_0","The d_0 parameter of the switching function"); - keys.add("optional","R_0","The r_0 parameter of the switching function"); + keys.add("compulsory","NN","6","The n parameter of the switching function "); + keys.add("compulsory","MM","12","The m parameter of the switching function "); + keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function"); + keys.add("compulsory","R_0","The r_0 parameter of the switching function"); + keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. " + "The following provides information on the \\ref switchingfunction that are available. " + "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords."); // Use actionWithDistributionKeywords keys.use("AVERAGE"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MIN"); keys.use("WITHIN"); keys.use("MOMENTS"); diff --git a/src/MultiColvarDensity.cpp b/src/MultiColvarDensity.cpp index 3d958cef4..0b61fe024 100644 --- a/src/MultiColvarDensity.cpp +++ b/src/MultiColvarDensity.cpp @@ -33,16 +33,15 @@ namespace PLMD{ //+PLUMEDOC MCOLVAR DENSITY /* Calculate functions of the density of atoms as a function of the box. This allows one to calculate -density gradients, number of atoms in half the box and so on. +the number of atoms in half the box. \par Examples The following example calculates the number of atoms in one half of the simulation box. \verbatim -DENSITY SPECIES=1-100 SUBCELL=(XLOWER=0.0 XUPPER=0.5) LABEL=d1 +DENSITY SPECIES=1-100 SUBCELL={XLOWER=0.0 XUPPER=0.5} LABEL=d1 PRINT ARG=d1.* FILE=colvar1 FMT=%8.4f - \endverbatim */ diff --git a/src/MultiColvarDistance.cpp b/src/MultiColvarDistance.cpp index 4cbdd9032..7cc1b8cc3 100644 --- a/src/MultiColvarDistance.cpp +++ b/src/MultiColvarDistance.cpp @@ -37,17 +37,11 @@ distances such as the minimum, the number less than a certain quantity and so on \par Examples -The following input tells plumed to print the distance between atoms 3 and 5, +The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2 and to +print the minimum for these two distances. \verbatim -DISTANCES ATOMS=3,5 LABEL=d1 -PRINT ARG=d1.* -\endverbatim -(See also \ref PRINT). - -The following input tells plumed to print the distances between atoms 3 and 5 and between atoms 1 and 2 -\verbatim -DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 -PRINT ARG=d1.* +DISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1 +PRINT ARG=d1.min \endverbatim (See also \ref PRINT). @@ -55,11 +49,27 @@ The following input tells plumed to calculate the distances between atoms 3 and and then to calculate the number of these distances that are less than 0.1 nm. The number of distances less than 0.1nm is then printed to a file. \verbatim -DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN=0.1 +DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1} PRINT ARG=d1.lt0.1 \endverbatim -(See also \ref PRINT). +(See also \ref PRINT \ref switchingfunction). +The following input tells plumed to calculate all the distances between atoms 1, 2 and 3 (i.e. the distances between atoms +1 and 2, atoms 1 and 3 and atoms 2 and 3). The average of these distances is then calculated. +\verbatim +DISTANCES GROUP=1-3 AVERAGE LABEL=d1 +PRINT ARG=d1.average +\endverbatim +(See also \ref PRINT) + +The following input tells plumed to calculate all the distances between the atoms in GROUPA and the atoms in GROUPB. +In other words the distances between atoms 1 and 2 and the distance between atoms 1 and 3. The number of distances +more than 0.1 is then printed to a file. +\verbatim +DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1} +PRINT ARG=d1.gt0.1 +\endverbatim +(See also \ref PRINT \ref switchingfunction) */ //+ENDPLUMEDOC diff --git a/src/MultiColvarSecondaryStructureRMSD.cpp b/src/MultiColvarSecondaryStructureRMSD.cpp index d285fed87..3667445e1 100644 --- a/src/MultiColvarSecondaryStructureRMSD.cpp +++ b/src/MultiColvarSecondaryStructureRMSD.cpp @@ -37,7 +37,9 @@ void MultiColvarSecondaryStructureRMSD::registerKeywords( Keywords& keys ){ "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain " "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. " "As such if you define portions of the chain with fewer than N residues the code will crash."); - keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD."); + keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. " + "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD. For more details on the " + "DRMSD method see \\ref DRMSD."); keys.use("LESS_THAN"); keys.use("MIN"); keys.use("AVERAGE"); } diff --git a/src/PDB.cpp b/src/PDB.cpp index e73f9e54a..9d9980fb2 100644 --- a/src/PDB.cpp +++ b/src/PDB.cpp @@ -28,15 +28,6 @@ using namespace std; namespace PLMD{ -std::string PDB::documentation(){ - std::ostringstream ostr; - ostr<<"In PDB files the atomic coordinates and box lengths should be in Angstroms unless you are working with natural units. "; - ostr<<"If you are working with natural units then the coordinates should be in your natural length unit. In the PDB files used "; - ostr<<"by plumed the beta column is used to specify the charges on the atoms and the occupancy column is used to specify the atomic masses. "; - ostr<<"For more details on the PDB file format visit http://www.wwpdb.org/docs.html"; - return ostr.str(); -} - const std::vector<Vector> & PDB::getPositions()const{ return positions; } diff --git a/src/PDB.h b/src/PDB.h index 58876820e..4a6dfd794 100644 --- a/src/PDB.h +++ b/src/PDB.h @@ -40,8 +40,6 @@ class PDB{ std::vector<double> beta; std::vector<AtomNumber> numbers; public: -/// The documentation for PDB file parser - static std::string documentation(); /// Read the pdb from a file, scaling positions by a factor scale bool read(const std::string&file,bool naturalUnits,double scale); /// Access to the position array diff --git a/src/SetupMolInfo.cpp b/src/SetupMolInfo.cpp index aec4677ed..5f096ac65 100644 --- a/src/SetupMolInfo.cpp +++ b/src/SetupMolInfo.cpp @@ -63,8 +63,9 @@ PLUMED_REGISTER_ACTION(SetupMolInfo,"MOLINFO") void SetupMolInfo::registerKeywords( Keywords& keys ){ ActionSetup::registerKeywords(keys); keys.add("compulsory","STRUCTURE","a file in pdb format containing a reference structure. " - "This is used to defines the atoms in the various residues, chains, etc . " + PDB::documentation() ); - keys.add("atoms","CHAIN","(for masochists ( a.k.a. Davide Branduardi ) ) The atoms involved in each of the chains of interest in the structure."); + "This is used to defines the atoms in the various residues, chains, etc . " + "For more details on the PDB file format visit http://www.wwpdb.org/docs.html"); + keys.add("atoms","CHAIN","(for masochists ( mostly Davide Branduardi ) ) The atoms involved in each of the chains of interest in the structure."); } SetupMolInfo::SetupMolInfo( const ActionOptions&ao ): diff --git a/src/SwitchingFunction.cpp b/src/SwitchingFunction.cpp index 83afc1114..d4b12b2b7 100644 --- a/src/SwitchingFunction.cpp +++ b/src/SwitchingFunction.cpp @@ -28,19 +28,58 @@ using namespace std; using namespace PLMD; -std::string SwitchingFunction::documentation(){ - std::ostringstream ostr; - ostr<<"Within plumed you can use the following switching functions "; - ostr<<"\\f$s(r)=\\frac{ 1 - \\left(\\frac{ r - d_0 }{ r_0 }\\right)^{n} }{ 1 - \\left(\\frac{ r - d_0 }{ r_0 }\\right)^{m} } \\f$, "; - ostr<<"\\f$s(r)=\\exp\\left(-\\frac{ r - d_0 }{ r_0 }\\right)\\f$ or using \\f$s(r)=\\exp\\left(-\\frac{ (r - d_0)^2 }{ 2r_0^2 }\\right)\\f$. "; - ostr<<"The first of these options is specified using the syntax {RATIONAL R_0=\\f$r_0\\f$ D_0=\\f$d_0\\f$ NN=\\f$n\\f$ MM=\\f$m\\f$} and if "; - ostr<<"the D_0, NN and MM keywords are missing they are assumed equal to 0, 6 and 12 respectively. The second form is specified using "; - ostr<<"{EXP R_0=\\f$r_0\\f$ D_0=\\f$d_0\\f$} and if the D_0 is missing it is assumed equal to 0. The third form is specified using "; - ostr<<"{GAUSSIAN R_0=\\f$r_0\\f$ D_0=\\f$d_0\\f$} and if the D_0 is missing it is assumed equal to 0. You can add the D_MAX flag to "; - ostr<<"all switching function definitions to specify that if the input \\f$r\\f$ is greater than this value the value of the switching "; - ostr<<"function is very small and can thus be assumed to equal zero."; - return ostr.str(); -} +//+PLUMEDOC INTERNAL switchingfunction +/* + +Switching functions \f$s(r)\f$ take a minimum of one input parameter \f$d_0\f$. +For \f$r \le d_0 \quad s(r)=1.0\f$ while for \f$r > d_0\f$ the function decays smoothly to 0. +The various switching functions available in plumed differ in terms of how this decay is performed. + +Where there is an accepted convention in the literature (e.g. \ref COORDINATION) on the form of the +switching function we use the convention as the default. However, the flexibility to use different +switching functions is always present generally through a single keyword. This keyword generally +takes an input with the following form: + +\verbatim +KEYWORD={TYPE <list of parameters>} +\endverbatim + +The following table contains a list of the various switching functions that are available in plumed 2 +together with an example input. + +<table align=center frame=void width=95%% cellpadding=5%%> +<tr> +<td> TYPE </td> <td> FUNCTION </td> <td> EXAMPLE INPUT </td> <td> DEFAULT PARAMETERS </td> +</tr> <tr> <td>RATIONAL </td> <td> +\f$ +s(r)=\frac{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{n} }{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{m} } +\f$ +</td> <td> +{RATIONAL R_0=\f$r_0\f$ D_0=\f$d_0\f$ NN=\f$n\f$ MM=\f$m\f$} +</td> <td> \f$d_0=0.0\f$, \f$n=6\f$, \f$m=12\f$ </td> +</tr> <tr> +<td> EXP </td> <td> +\f$ +s(r)=\exp\left(-\frac{ r - d_0 }{ r_0 }\right) +\f$ +</td> <td> +{EXP R_0=\f$r_0\f$ D_0=\f$d_0\f$} +</td> <td> \f$d_0=0.0\f$ </td> +</tr> <tr> +<td> GAUSSIAN </td> <td> +\f$ +s(r)=\exp\left(-\frac{ (r - d_0)^2 }{ 2r_0^2 }\right) +\f$ +</td> <td> +{GAUSSIAN R_0=\f$r_0\f$ D_0=\f$d_0\f$} +</td> <td> \f$d_0=0.0\f$ </td> +</tr> +</table> + +For all the switching functions in the above table one can also specify a further (optional) parameter using the parameter +keyword D_MAX to assert that for \f$r>d_{\textrm{max}}\f$ the switching function can be assumed equal to zero. +*/ +//+ENDPLUMEDOC void SwitchingFunction::set(const std::string & definition,std::string& errormsg){ vector<string> data=Tools::getWords(definition); diff --git a/src/SwitchingFunction.h b/src/SwitchingFunction.h index 07ebcb842..6ea12ca00 100644 --- a/src/SwitchingFunction.h +++ b/src/SwitchingFunction.h @@ -41,7 +41,6 @@ class SwitchingFunction{ int nn,mm; double invr0,d0,dmax; public: - static std::string documentation(); SwitchingFunction(); void set(int nn,int mm,double r_0,double d_0); void set(const std::string& definition, std::string& errormsg); diff --git a/user-doc/Colvar.txt b/user-doc/Colvar.txt index 16650547a..0c8755863 100644 --- a/user-doc/Colvar.txt +++ b/user-doc/Colvar.txt @@ -69,18 +69,7 @@ periodic boundaries, a fact which plumed could only deal with were the topology work would involve a lot laborious programming and goes against our original aim of having a general patch that can be implemented in a wide variety of MD codes. Consequentially, we have implemented a more pragmatic solution to this probem - the user specifies in input any molecules (or parts of molecules) that must be kept in tact throughout the simulation run. In plumed 1.0 this was done -using the ALIGN_ATOMS keyword. In plumed 2.0 the same effect can be acchieved using any one from the following list of commands - -<table align=center frame=void width=95%% cellpadding=5%%> -<tr> -<td width=5%> \subpage WHOLEMOLECULES </td> <td> With this command you can specify the atoms in each of your molecules </td> -</tr> -</table> -@TOPOLOGY@ - -The most general of these is the \ref WHOLEMOLECULES keywords. However, there are advantages in using the others as oftentimes -by defining a topology using a pdb file you can take advantage of shortcuts when specifying other colvars. For example, if you -use the \ref PROTEIN_TOPOLOGY keyword you can then later refer to the atoms in a specific residue by the residue number. +using the ALIGN_ATOMS keyword. In plumed 2.0 the same effect can be acchieved using the \subpage WHOLEMOLECULES command. \section cvavail CV Documentation diff --git a/user-doc/extract b/user-doc/extract index 1c1235a20..3ae7ff600 100755 --- a/user-doc/extract +++ b/user-doc/extract @@ -37,16 +37,16 @@ awk 'BEGIN{gfile="automatic/GLOSSARY1.list"}{ inside=0; print "*/" >output } - if(inside==2 && NF==0){ + if(inside>=2 && NF==0){ print "</td> </tr>" > lfile - printf "</td> </tr>\n" > gfile + if(inside==2){ printf "</td> </tr>\n" > gfile; } inside=1 } if(inside==1 && $1!="/*" && $1!="*/") print $0 >output - if(inside==2 && $1!="/*" && $1!="*/" ){ + if(inside>=2 && $1!="/*" && $1!="*/" ){ print $0 > output print $0 > lfile - printf "%s", $0 > gfile + if(inside==2){ printf "%s", $0 > gfile; } } if($1=="//+PLUMEDOC"){ if( $2=="TOPOLOGY" || $2=="COLVAR" || $2=="MCOLVAR" || $2=="FUNCTION" || $2=="ANALYSIS" || $2=="BIAS" || $2=="GENERIC" || $2=="VATOM" || $2=="TOOLS" ){ @@ -61,6 +61,13 @@ awk 'BEGIN{gfile="automatic/GLOSSARY1.list"}{ print "\\page "$3 > output print "\\section "$3 >output inside=2; + } else if ( $2=="INTERNAL" ){ + if(output!="")close(output); + output="automatic/"$3".txt"; + print "/**" > output + print "\\page "$3 > output + print "\\section "$3 > output + inside=3; } } }' -- GitLab