Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2018-2023 The plumed team 3 : (see the PEOPLE file at the root of the distribution for a list of names) 4 : 5 : See http://www.plumed.org for more information. 6 : 7 : This file is part of plumed, version 2. 8 : 9 : plumed is free software: you can redistribute it and/or modify 10 : it under the terms of the GNU Lesser General Public License as published by 11 : the Free Software Foundation, either version 3 of the License, or 12 : (at your option) any later version. 13 : 14 : plumed is distributed in the hope that it will be useful, 15 : but WITHOUT ANY WARRANTY; without even the implied warranty of 16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 : GNU Lesser General Public License for more details. 18 : 19 : You should have received a copy of the GNU Lesser General Public License 20 : along with plumed. If not, see <http://www.gnu.org/licenses/>. 21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ 22 : #include "core/ActionShortcut.h" 23 : #include "core/PlumedMain.h" 24 : #include "core/ActionSet.h" 25 : #include "core/ActionRegister.h" 26 : 27 : namespace PLMD { 28 : namespace wham { 29 : 30 : //+PLUMEDOC REWEIGHTING WHAM_HISTOGRAM 31 : /* 32 : This can be used to output the a histogram using the weighted histogram technique 33 : 34 : This shortcut action allows you to calculate a histogram using the weighted histogram 35 : analysis technique. For more detail on how this the weights for configurations are 36 : computed see \ref REWEIGHT_WHAM 37 : 38 : \par Examples 39 : 40 : The following input can be used to analyze the output from a series of umbrella sampling calculations. 41 : The trajectory from each of the simulations run with the different biases should be concatenated into a 42 : single trajectory before running the following analysis script on the concatenated trajectory using PLUMED 43 : driver. The umbrella sampling simulations that will be analyzed using the script below applied a harmonic 44 : restraint that restrained the torsional angle involving atoms 5, 7, 9 and 15 to particular values. The script 45 : below calculates the reweighting weights for each of the trajectories and then applies the binless WHAM algorithm 46 : to determine a weight for each configuration in the concatenated trajectory. A histogram is then constructed from 47 : the configurations visited and their weights. This histogram is then converted into a free energy surface and output 48 : to a file called fes.dat 49 : 50 : \plumedfile 51 : #SETTINGS NREPLICAS=4 52 : phi: TORSION ATOMS=5,7,9,15 53 : psi: TORSION ATOMS=7,9,15,17 54 : rp: RESTRAINT ARG=phi KAPPA=50.0 ... 55 : AT=@replicas:{ 56 : -3.00000000000000000000 57 : -1.45161290322580645168 58 : .09677419354838709664 59 : 1.64516129032258064496 60 : } 61 : ... 62 : 63 : hh: WHAM_HISTOGRAM ARG=phi BIAS=rp.bias TEMP=300 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 64 : fes: CONVERT_TO_FES GRID=hh TEMP=300 65 : DUMPGRID GRID=fes FILE=fes.dat 66 : \endplumedfile 67 : 68 : The script above must be run with multiple replicas using the following command: 69 : 70 : \verbatim 71 : mpirun -np 4 plumed driver --mf_xtc alltraj.xtc --multi 4 72 : \endverbatim 73 : 74 : */ 75 : //+ENDPLUMEDOC 76 : 77 : class WhamHistogram : public ActionShortcut { 78 : public: 79 : static void registerKeywords( Keywords& keys ); 80 : explicit WhamHistogram( const ActionOptions& ); 81 : }; 82 : 83 : PLUMED_REGISTER_ACTION(WhamHistogram,"WHAM_HISTOGRAM") 84 : 85 10 : void WhamHistogram::registerKeywords( Keywords& keys ) { 86 10 : ActionShortcut::registerKeywords( keys ); 87 20 : keys.add("compulsory","ARG","the arguments that you would like to make the histogram for"); 88 20 : keys.add("compulsory","BIAS","*.bias","the value of the biases to use when performing WHAM"); 89 20 : keys.add("compulsory","TEMP","the temperature at which the simulation was run"); 90 20 : keys.add("compulsory","STRIDE","1","the frequency with which the data should be stored to perform WHAM"); 91 20 : keys.add("compulsory","GRID_MIN","the minimum to use for the grid"); 92 20 : keys.add("compulsory","GRID_MAX","the maximum to use for the grid"); 93 20 : keys.add("compulsory","GRID_BIN","the number of bins to use for the grid"); 94 20 : keys.add("optional","BANDWIDTH","the bandwidth for kernel density estimation"); 95 10 : keys.setValueDescription("the histogram that was generated using the WHAM weights"); 96 20 : keys.needsAction("GATHER_REPLICAS"); keys.needsAction("CONCATENATE"); 97 30 : keys.needsAction("COLLECT"); keys.needsAction("WHAM"); keys.needsAction("KDE"); 98 10 : } 99 : 100 : 101 7 : WhamHistogram::WhamHistogram( const ActionOptions& ao ) : 102 : Action(ao), 103 7 : ActionShortcut(ao) 104 : { 105 : // Input for collection of weights for WHAM 106 14 : std::string bias; parse("BIAS",bias); 107 7 : std::string stride; parse("STRIDE",stride); 108 : // Input for GATHER_REPLICAS 109 14 : readInputLine( getShortcutLabel() + "_gather: GATHER_REPLICAS ARG=" + bias ); 110 : // Put all the replicas in a single vector 111 14 : readInputLine( getShortcutLabel() + "_gatherv: CONCATENATE ARG=" + getShortcutLabel() + "_gather.*"); 112 : // Input for COLLECT_FRAMES 113 14 : readInputLine( getShortcutLabel() + "_collect: COLLECT TYPE=vector ARG=" + getShortcutLabel() + "_gatherv STRIDE=" + stride); 114 : // Input for WHAM 115 21 : std::string temp, tempstr=""; parse("TEMP",temp); if( temp.length()>0 ) tempstr="TEMP=" + temp; 116 14 : readInputLine( getShortcutLabel() + "_wham: WHAM ARG=" + getShortcutLabel() + "_collect " + tempstr ); 117 : // Input for COLLECT_FRAMES 118 14 : std::vector<std::string> args; parseVector("ARG",args); std::string argstr; 119 14 : for(unsigned i=0; i<args.size(); ++i) { 120 14 : readInputLine( getShortcutLabel() + "_data_" + args[i] + ": COLLECT ARG=" + args[i] ); 121 7 : if( i==0 ) argstr = " ARG="; else argstr += ","; 122 14 : argstr += getShortcutLabel() + "_data_" + args[i]; 123 : } 124 : // Input for HISTOGRAM 125 14 : std::string histo_line, bw=""; parse("BANDWIDTH",bw); 126 7 : if( bw!="" ) histo_line += " BANDWIDTH=" + bw; 127 : else histo_line += " KERNEL=DISCRETE"; 128 21 : std::string min; parse("GRID_MIN",min); histo_line += " GRID_MIN=" + min; 129 21 : std::string max; parse("GRID_MAX",max); histo_line += " GRID_MAX=" + max; 130 14 : std::string bin; parse("GRID_BIN",bin); histo_line += " GRID_BIN=" + bin; 131 14 : readInputLine( getShortcutLabel() + ": KDE " + argstr + " HEIGHTS=" + getShortcutLabel() + "_wham" + histo_line ); 132 14 : } 133 : 134 : } 135 : }