1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-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
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 lugano-2
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 \plumedfile
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 \endplumedfile
59
60 */
61 //+ENDPLUMEDOC
62
63 class Restraint : public Bias {
64 std::vector<double> at;
65 std::vector<double> kappa;
66 std::vector<double> slope;
67 Value* valueForce2;
68 public:
69 explicit Restraint(const ActionOptions&);
70 void calculate() override;
71 static void registerKeywords(Keywords& keys);
72 };
73
74 PLUMED_REGISTER_ACTION(Restraint,"RESTRAINT")
75
registerKeywords(Keywords & keys)76 void Restraint::registerKeywords(Keywords& keys) {
77 Bias::registerKeywords(keys);
78 keys.use("ARG");
79 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");
80 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");
81 keys.add("compulsory","AT","the position of the restraint");
82 keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
83 }
84
Restraint(const ActionOptions & ao)85 Restraint::Restraint(const ActionOptions&ao):
86 PLUMED_BIAS_INIT(ao),
87 at(getNumberOfArguments()),
88 kappa(getNumberOfArguments(),0.0),
89 slope(getNumberOfArguments(),0.0)
90 {
91 parseVector("SLOPE",slope);
92 parseVector("KAPPA",kappa);
93 parseVector("AT",at);
94 checkRead();
95
96 log.printf(" at");
97 for(unsigned i=0; i<at.size(); i++) log.printf(" %f",at[i]);
98 log.printf("\n");
99 log.printf(" with harmonic force constant");
100 for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
101 log.printf("\n");
102 log.printf(" and linear force constant");
103 for(unsigned i=0; i<slope.size(); i++) log.printf(" %f",slope[i]);
104 log.printf("\n");
105
106 addComponent("force2");
107 componentIsNotPeriodic("force2");
108 valueForce2=getPntrToComponent("force2");
109 }
110
111
calculate()112 void Restraint::calculate() {
113 double ene=0.0;
114 double totf2=0.0;
115 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
116 const double cv=difference(i,at[i],getArgument(i));
117 const double k=kappa[i];
118 const double m=slope[i];
119 const double f=-(k*cv+m);
120 ene+=0.5*k*cv*cv+m*cv;
121 setOutputForce(i,f);
122 totf2+=f*f;
123 }
124 setBias(ene);
125 valueForce2->set(totf2);
126 }
127
128 }
129
130
131 }
132