Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 "core/ActionRegister.h"
24 :
25 : namespace PLMD {
26 : namespace bias {
27 :
28 : //+PLUMEDOC BIAS RESTRAINT
29 : /*
30 : Adds harmonic and/or linear restraints on one or more variables.
31 :
32 : Either or both of SLOPE and KAPPA must be present to specify the linear and harmonic force constants
33 : respectively. If the input arguments are scalars the resulting potential is given by:
34 :
35 : $$
36 : \sum_i \frac{k_i}{2} (x_i-a_i)^2 + m_i*(x_i-a_i)
37 : $$.
38 :
39 : The following example input illustrates how you can apply a potential of this type. This input
40 : tells plumed to restrain the distance between atoms 3 and 5 and the distance between atoms 2 and 4, at different equilibrium
41 : values, and to print the energy of the restraint
42 :
43 : ```plumed
44 : d1: DISTANCE ATOMS=3,5
45 : d2: DISTANCE ATOMS=2,4
46 : restraint: RESTRAINT ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0
47 : PRINT ARG=restraint.bias
48 : ```
49 :
50 : If by contrast the inputs are vectors as in the example below:
51 :
52 : ```plumed
53 : d1: DISTANCE ATOMS1=3,5 ATOMS2=5,6 ATOMS3=6,7
54 : d2: DISTANCE ATOMS1=2,4 ATOMS2=4,6 ATOMS3=6,8
55 : restraint: RESTRAINT ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0
56 : PRINT ARG=restraint.bias
57 : ```
58 :
59 : Then the resulting potential is given by:
60 :
61 : $$
62 : \sum_i \sum_j \frac{k_i}{2} (x_{ij}-a_i)^2 + m_i*(x_{ij}-a_i)
63 : $$
64 :
65 : The sum over $i$ here is runs over the two arguments, while the sum over $j$ runs over the three components of the input vectors.
66 : Notice, that regardless of whether the input is a scalar, vector or matrix the number of $k_i$, $a_i$ and $m_i$ values must be
67 : equal to the number of arguments to the action.
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 : class Restraint : public Bias {
73 : std::vector<double> at;
74 : std::vector<double> kappa;
75 : std::vector<double> slope;
76 : Value* valueForce2;
77 : public:
78 : explicit Restraint(const ActionOptions&);
79 : void calculate() override;
80 : static void registerKeywords(Keywords& keys);
81 : };
82 :
83 : PLUMED_REGISTER_ACTION(Restraint,"RESTRAINT_SCALAR")
84 :
85 782 : void Restraint::registerKeywords(Keywords& keys) {
86 782 : Bias::registerKeywords(keys);
87 1564 : keys.setDisplayName("RESTRAINT");
88 782 : 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");
89 782 : 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");
90 782 : keys.add("compulsory","AT","the position of the restraint");
91 1564 : keys.addOutputComponent("force2","default","scalar","the instantaneous value of the squared force due to this bias potential");
92 782 : }
93 :
94 389 : Restraint::Restraint(const ActionOptions&ao):
95 : PLUMED_BIAS_INIT(ao),
96 776 : at(getNumberOfArguments()),
97 388 : kappa(getNumberOfArguments(),0.0),
98 777 : slope(getNumberOfArguments(),0.0) {
99 388 : parseVector("SLOPE",slope);
100 388 : parseVector("KAPPA",kappa);
101 388 : parseVector("AT",at);
102 388 : checkRead();
103 :
104 388 : log.printf(" at");
105 892 : for(unsigned i=0; i<at.size(); i++) {
106 504 : log.printf(" %f",at[i]);
107 : }
108 388 : log.printf("\n");
109 388 : log.printf(" with harmonic force constant");
110 892 : for(unsigned i=0; i<kappa.size(); i++) {
111 504 : log.printf(" %f",kappa[i]);
112 : }
113 388 : log.printf("\n");
114 388 : log.printf(" and linear force constant");
115 892 : for(unsigned i=0; i<slope.size(); i++) {
116 504 : log.printf(" %f",slope[i]);
117 : }
118 388 : log.printf("\n");
119 :
120 776 : addComponent("force2");
121 388 : componentIsNotPeriodic("force2");
122 388 : valueForce2=getPntrToComponent("force2");
123 389 : }
124 :
125 :
126 5393 : void Restraint::calculate() {
127 : double ene=0.0;
128 : double totf2=0.0;
129 11634 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
130 6241 : const double cv=difference(i,at[i],getArgument(i));
131 6241 : const double k=kappa[i];
132 6241 : const double m=slope[i];
133 6241 : const double force=-(k*cv+m);
134 6241 : ene+=0.5*k*cv*cv+m*cv;
135 6241 : setOutputForce(i,force);
136 6241 : totf2+=force*force;
137 : }
138 5393 : setBias(ene);
139 5393 : valueForce2->set(totf2);
140 5393 : }
141 :
142 : }
143 :
144 :
145 : }
|
