1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2012-2021 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 "Bias.h"
23 #include "ActionRegister.h"
24 #include "tools/Grid.h"
25 #include "tools/Exception.h"
26 #include "tools/File.h"
27 #include <memory>
28
29
30 using namespace std;
31
32
33 namespace PLMD {
34 namespace bias {
35
36 //+PLUMEDOC BIAS EXTERNAL
37 /*
38 Calculate a restraint that is defined on a grid that is read during start up
39
40 \par Examples
41
42 The following is an input for a calculation with an external potential that is
43 defined in the file bias.dat and that acts on the distance between atoms 3 and 5.
44 \plumedfile
45 DISTANCE ATOMS=3,5 LABEL=d1
46 EXTERNAL ARG=d1 FILE=bias.grid LABEL=external
47 \endplumedfile
48
49 The bias.grid will then look something like this:
50 \auxfile{bias.grid}
51 #! FIELDS d1 external.bias der_d1
52 #! SET min_d1 1.14
53 #! SET max_d1 1.32
54 #! SET nbins_d1 6
55 #! SET periodic_d1 false
56 1.1400 0.0031 0.1101
57 1.1700 0.0086 0.2842
58 1.2000 0.0222 0.6648
59 1.2300 0.0521 1.4068
60 1.2600 0.1120 2.6873
61 1.2900 0.2199 4.6183
62 1.3200 0.3948 7.1055
63 \endauxfile
64
65 This should then be followed by the value of the potential and its derivative
66 at 100 equally spaced points along the distance between 0 and 1.
67
68 You can also include grids that are a function of more than one collective
69 variable. For instance the following would be the input for an external
70 potential acting on two torsional angles:
71 \plumedfile
72 TORSION ATOMS=4,5,6,7 LABEL=t1
73 TORSION ATOMS=6,7,8,9 LABEL=t2
74 EXTERNAL ARG=t1,t2 FILE=bias2.grid LABEL=ext
75 \endplumedfile
76
77 The file bias2.grid for this calculation would need to look something like this:
78 \auxfile{bias2.grid}
79 #! FIELDS t1 t2 ext.bias der_t1 der_t2
80 #! SET min_t1 -pi
81 #! SET max_t1 pi
82 #! SET nbins_t1 3
83 #! SET periodic_t1 true
84 #! SET min_t2 -pi
85 #! SET max_t2 pi
86 #! SET nbins_t2 3
87 #! SET periodic_t2 true
88 -3.141593 -3.141593 0.000000 -0.000000 -0.000000
89 -1.047198 -3.141593 0.000000 0.000000 -0.000000
90 1.047198 -3.141593 0.000000 -0.000000 -0.000000
91
92 -3.141593 -1.047198 0.000000 -0.000000 0.000000
93 -1.047198 -1.047198 0.007922 0.033185 0.033185
94 1.047198 -1.047198 0.007922 -0.033185 0.033185
95
96 -3.141593 1.047198 0.000000 -0.000000 -0.000000
97 -1.047198 1.047198 0.007922 0.033185 -0.033185
98 1.047198 1.047198 0.007922 -0.033185 -0.033185
99 \endauxfile
100
101 This would be then followed by 100 blocks of data. In the first block of data the
102 value of t1 (the value in the first column) is kept fixed and the value of
103 the function is given at 100 equally spaced values for t2 between \f$-pi\f$ and \f$+pi\f$. In the
104 second block of data t1 is fixed at \f$-pi + \frac{2pi}{100}\f$ and the value of the function is
105 given at 100 equally spaced values for t2 between \f$-pi\f$ and \f$+pi\f$. In the third block of
106 data the same is done but t1 is fixed at \f$-pi + \frac{4pi}{100}\f$ and so on until you get to
107 the one hundredth block of data where t1 is fixed at \f$+pi\f$.
108
109 Please note the order that the order of arguments in the plumed.dat file must be the same as
110 the order of arguments in the header of the grid file.
111 */
112 //+ENDPLUMEDOC
113
114 class External : public Bias {
115
116 private:
117 std::unique_ptr<GridBase> BiasGrid_;
118 double scale_;
119
120 public:
121 explicit External(const ActionOptions&);
122 void calculate() override;
123 static void registerKeywords(Keywords& keys);
124 };
125
126 PLUMED_REGISTER_ACTION(External,"EXTERNAL")
127
registerKeywords(Keywords & keys)128 void External::registerKeywords(Keywords& keys) {
129 Bias::registerKeywords(keys);
130 keys.use("ARG");
131 keys.add("compulsory","FILE","the name of the file containing the external potential.");
132 keys.addFlag("NOSPLINE",false,"specifies that no spline interpolation is to be used when calculating the energy and forces due to the external potential");
133 keys.addFlag("SPARSE",false,"specifies that the external potential uses a sparse grid");
134 keys.add("compulsory","SCALE","1.0","a factor that multiplies the external potential, useful to invert free energies");
135 }
136
External(const ActionOptions & ao)137 External::External(const ActionOptions& ao):
138 PLUMED_BIAS_INIT(ao)
139 {
140 string filename;
141 parse("FILE",filename);
142 if( filename.length()==0 ) error("No external potential file was specified");
143 bool sparsegrid=false;
144 parseFlag("SPARSE",sparsegrid);
145 bool nospline=false;
146 parseFlag("NOSPLINE",nospline);
147 bool spline=!nospline;
148 parse("SCALE",scale_);
149
150 checkRead();
151
152 log.printf(" External potential from file %s\n",filename.c_str());
153 log.printf(" Multiplied by %lf\n",scale_);
154 if(spline) {log.printf(" External potential uses spline interpolation\n");}
155 if(sparsegrid) {log.printf(" External potential uses sparse grid\n");}
156
157 // read grid
158 IFile gridfile; gridfile.open(filename);
159 std::string funcl=getLabel() + ".bias";
160 BiasGrid_=GridBase::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
161 if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
162 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
163 if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
164 }
165 }
166
calculate()167 void External::calculate()
168 {
169 unsigned ncv=getNumberOfArguments();
170 vector<double> cv(ncv), der(ncv);
171
172 for(unsigned i=0; i<ncv; ++i) {cv[i]=getArgument(i);}
173
174 double ene=scale_*BiasGrid_->getValueAndDerivatives(cv,der);
175
176 setBias(ene);
177
178 for(unsigned i=0; i<ncv; ++i) {
179 const double f=-scale_*der[i];
180 setOutputForce(i,f);
181 }
182 }
183
184 }
185 }
186