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