Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/Atoms.h"
25 : #include "ActionWithInputGrid.h"
26 :
27 : //+PLUMEDOC GRIDANALYSIS CONVERT_TO_FES
28 : /*
29 : Convert a histogram, \f$H(x)\f$, to a free energy surface using \f$F(x) = -k_B T \ln H(x)\f$.
30 :
31 : This action allows you to take a free energy surface that was calculated using the \ref HISTOGRAM
32 : action and to convert it to a free energy surface. This transformation performed by doing:
33 :
34 : \f[
35 : F(x) = -k_B T \ln H(x)
36 : \f]
37 :
38 : The free energy calculated on a grid is output by this action and can be printed using \ref DUMPGRID
39 :
40 : \par Examples
41 :
42 : This is a typical example showing how CONVERT_TO_FES might be used when postprocessing a trajectory.
43 : The input below calculates the free energy as a function of the distance between atom 1 and atom 2.
44 : This is done by accumulating a histogram as a function of this distance using kernel density estimation
45 : and the HISTOGRAM action. All the data within this trajectory is used in the construction of this
46 : HISTOGRAM. Finally, once all the data has been read in, the histogram is converted to a free energy
47 : using the formula above and the free energy is output to a file called fes.dat
48 :
49 : \plumedfile
50 : x: DISTANCE ATOMS=1,2
51 : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
52 : ff: CONVERT_TO_FES GRID=hA1 TEMP=300
53 : DUMPGRID GRID=ff FILE=fes.dat
54 : \endplumedfile
55 :
56 : */
57 : //+ENDPLUMEDOC
58 :
59 : namespace PLMD {
60 : namespace gridtools {
61 :
62 14 : class ConvertToFES : public ActionWithInputGrid {
63 : private:
64 : double simtemp;
65 : bool activated;
66 : public:
67 : static void registerKeywords( Keywords& keys );
68 : explicit ConvertToFES(const ActionOptions&ao);
69 : unsigned getNumberOfQuantities() const ;
70 12 : void prepare() { activated=true; }
71 15 : void prepareForAveraging() { ActionWithInputGrid::prepareForAveraging(); activated=false; }
72 : void compute( const unsigned& current, MultiValue& myvals ) const ;
73 0 : bool isPeriodic() { return false; }
74 72 : bool onStep() const { return activated; }
75 : void runFinalJobs();
76 : };
77 :
78 6459 : PLUMED_REGISTER_ACTION(ConvertToFES,"CONVERT_TO_FES")
79 :
80 8 : void ConvertToFES::registerKeywords( Keywords& keys ) {
81 8 : ActionWithInputGrid::registerKeywords( keys );
82 32 : keys.add("optional","TEMP","the temperature at which you are operating");
83 32 : keys.remove("STRIDE"); keys.remove("KERNEL"); keys.remove("BANDWIDTH");
84 32 : keys.remove("LOGWEIGHTS"); keys.remove("CLEAR"); keys.remove("NORMALIZATION");
85 8 : }
86 :
87 7 : ConvertToFES::ConvertToFES(const ActionOptions&ao):
88 : Action(ao),
89 : ActionWithInputGrid(ao),
90 7 : activated(false)
91 : {
92 7 : plumed_assert( ingrid->getNumberOfComponents()==1 );
93 :
94 : // Create a grid
95 49 : createGrid( "grid", "COMPONENTS=" + getLabel() + " " + ingrid->getInputString() );
96 14 : if( ingrid->noDerivatives() ) mygrid->setNoDerivatives();
97 : std::vector<double> fspacing;
98 14 : mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), ingrid->getNbin(), fspacing);
99 7 : setAveragingAction( mygrid, true );
100 :
101 14 : simtemp=0.; parse("TEMP",simtemp);
102 14 : if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
103 0 : else simtemp=plumed.getAtoms().getKbT();
104 7 : if( simtemp==0 ) error("TEMP not set - use keyword TEMP");
105 :
106 : // Now create task list
107 5579 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) addTaskToList(i);
108 : // And activate all tasks
109 7 : deactivateAllTasks();
110 16709 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) taskFlags[i]=1;
111 7 : lockContributors();
112 7 : }
113 :
114 60 : unsigned ConvertToFES::getNumberOfQuantities() const {
115 120 : if( mygrid->noDerivatives() ) return 2;
116 56 : return 2 + mygrid->getDimension();
117 : }
118 :
119 11293 : void ConvertToFES::compute( const unsigned& current, MultiValue& myvals ) const {
120 11293 : double val=getFunctionValue( current ); myvals.setValue(1, -simtemp*std::log(val) );
121 22586 : if( !mygrid->noDerivatives() && val>0 ) {
122 95070 : for(unsigned i=0; i<mygrid->getDimension(); ++i) myvals.setValue( 2+i, -(simtemp/val)*ingrid->getGridElement(current,i+1) );
123 : }
124 11293 : }
125 :
126 7 : void ConvertToFES::runFinalJobs() {
127 7 : activated=true; update();
128 7 : }
129 :
130 : }
131 4839 : }
|
