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