Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 :
25 :
26 : using namespace std;
27 :
28 :
29 : namespace PLMD {
30 : namespace bias {
31 :
32 : //+PLUMEDOC BIAS RESTRAINT
33 : /*
34 : Adds harmonic and/or linear restraints on one or more variables.
35 :
36 : Either or both
37 : of SLOPE and KAPPA must be present to specify the linear and harmonic force constants
38 : respectively. The resulting potential is given by:
39 : \f[
40 : \sum_i \frac{k_i}{2} (x_i-a_i)^2 + m_i*(x_i-a_i)
41 : \f].
42 :
43 : The number of components for any vector of force constants must be equal to the number
44 : of arguments to the action.
45 :
46 : Additional material and examples can be also found in the tutorial \ref belfast-4
47 :
48 : \par Examples
49 :
50 : The following input tells plumed to restrain the distance between atoms 3 and 5
51 : and the distance between atoms 2 and 4, at different equilibrium
52 : values, and to print the energy of the restraint
53 : \verbatim
54 : DISTANCE ATOMS=3,5 LABEL=d1
55 : DISTANCE ATOMS=2,4 LABEL=d2
56 : RESTRAINT ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 LABEL=restraint
57 : PRINT ARG=restraint.bias
58 : \endverbatim
59 : (See also \ref DISTANCE and \ref PRINT).
60 :
61 : */
62 : //+ENDPLUMEDOC
63 :
64 214 : class Restraint : public Bias {
65 : std::vector<double> at;
66 : std::vector<double> kappa;
67 : std::vector<double> slope;
68 : Value* valueForce2;
69 : public:
70 : explicit Restraint(const ActionOptions&);
71 : void calculate();
72 : static void registerKeywords(Keywords& keys);
73 : };
74 :
75 2630 : PLUMED_REGISTER_ACTION(Restraint,"RESTRAINT")
76 :
77 108 : void Restraint::registerKeywords(Keywords& keys) {
78 108 : Bias::registerKeywords(keys);
79 108 : keys.use("ARG");
80 108 : keys.add("compulsory","SLOPE","0.0","specifies that the restraint is linear and what the values of the force constants on each of the variables are");
81 108 : keys.add("compulsory","KAPPA","0.0","specifies that the restraint is harmonic and what the values of the force constants on each of the variables are");
82 108 : keys.add("compulsory","AT","the position of the restraint");
83 108 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
84 108 : }
85 :
86 107 : Restraint::Restraint(const ActionOptions&ao):
87 : PLUMED_BIAS_INIT(ao),
88 107 : at(getNumberOfArguments()),
89 107 : kappa(getNumberOfArguments(),0.0),
90 321 : slope(getNumberOfArguments(),0.0)
91 : {
92 107 : parseVector("SLOPE",slope);
93 107 : parseVector("KAPPA",kappa);
94 107 : parseVector("AT",at);
95 107 : checkRead();
96 :
97 107 : log.printf(" at");
98 107 : for(unsigned i=0; i<at.size(); i++) log.printf(" %f",at[i]);
99 107 : log.printf("\n");
100 107 : log.printf(" with harmonic force constant");
101 107 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
102 107 : log.printf("\n");
103 107 : log.printf(" and linear force constant");
104 107 : for(unsigned i=0; i<slope.size(); i++) log.printf(" %f",slope[i]);
105 107 : log.printf("\n");
106 :
107 107 : addComponent("force2");
108 107 : componentIsNotPeriodic("force2");
109 107 : valueForce2=getPntrToComponent("force2");
110 107 : }
111 :
112 :
113 1065 : void Restraint::calculate() {
114 1065 : double ene=0.0;
115 1065 : double totf2=0.0;
116 2296 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
117 1231 : const double cv=difference(i,at[i],getArgument(i));
118 1231 : const double k=kappa[i];
119 1231 : const double m=slope[i];
120 1231 : const double f=-(k*cv+m);
121 1231 : ene+=0.5*k*cv*cv+m*cv;
122 1231 : setOutputForce(i,f);
123 1231 : totf2+=f*f;
124 : }
125 1065 : setBias(ene);
126 1065 : valueForce2->set(totf2);
127 1065 : }
128 :
129 : }
130 :
131 :
132 2523 : }
|
