Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 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 UPPER_WALLS
33 : /*
34 : Defines a wall for the value of one or more collective variables,
35 : which limits the region of the phase space accessible during the simulation.
36 :
37 : The restraining potential starts acting on the system when the value of the CV is greater
38 : (in the case of UPPER_WALLS) or lower (in the case of LOWER_WALLS) than a certain limit \f$a_i\f$ (AT)
39 : minus an offset \f$o_i\f$ (OFFSET).
40 : The expression for the bias due to the wall is given by:
41 :
42 : \f$
43 : \sum_i {k_i}((x_i-a_i+o_i)/s_i)^e_i
44 : \f$
45 :
46 : \f$k_i\f$ (KAPPA) is an energy constant in internal unit of the code, \f$s_i\f$ (EPS) a rescaling factor and
47 : \f$e_i\f$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFFSET = 0.
48 :
49 :
50 : \par Examples
51 :
52 : The following input tells plumed to add both a lower and an upper walls on the distance
53 : between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
54 : are defined at different values. The strength of the walls is the same for the four cases.
55 : It also tells plumed to print the energy of the walls.
56 : \plumedfile
57 : DISTANCE ATOMS=3,5 LABEL=d1
58 : DISTANCE ATOMS=2,4 LABEL=d2
59 : UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
60 : LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
61 : PRINT ARG=uwall.bias,lwall.bias
62 : \endplumedfile
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 51 : class UWalls : public Bias {
68 : std::vector<double> at;
69 : std::vector<double> kappa;
70 : std::vector<double> exp;
71 : std::vector<double> eps;
72 : std::vector<double> offset;
73 : public:
74 : explicit UWalls(const ActionOptions&);
75 : void calculate();
76 : static void registerKeywords(Keywords& keys);
77 : };
78 :
79 7390 : PLUMED_REGISTER_ACTION(UWalls,"UPPER_WALLS")
80 :
81 18 : void UWalls::registerKeywords(Keywords& keys) {
82 18 : Bias::registerKeywords(keys);
83 36 : keys.use("ARG");
84 72 : keys.add("compulsory","AT","the positions of the wall. The a_i in the expression for a wall.");
85 72 : keys.add("compulsory","KAPPA","the force constant for the wall. The k_i in the expression for a wall.");
86 90 : keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall. The o_i in the expression for a wall.");
87 90 : keys.add("compulsory","EXP","2.0","the powers for the walls. The e_i in the expression for a wall.");
88 90 : keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
89 72 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
90 18 : }
91 :
92 17 : UWalls::UWalls(const ActionOptions&ao):
93 : PLUMED_BIAS_INIT(ao),
94 : at(getNumberOfArguments(),0),
95 : kappa(getNumberOfArguments(),0.0),
96 : exp(getNumberOfArguments(),2.0),
97 : eps(getNumberOfArguments(),1.0),
98 102 : offset(getNumberOfArguments(),0.0)
99 : {
100 : // Note : the sizes of these vectors are checked automatically by parseVector
101 34 : parseVector("OFFSET",offset);
102 34 : parseVector("EPS",eps);
103 34 : parseVector("EXP",exp);
104 34 : parseVector("KAPPA",kappa);
105 34 : parseVector("AT",at);
106 17 : checkRead();
107 :
108 17 : log.printf(" at");
109 85 : for(unsigned i=0; i<at.size(); i++) log.printf(" %f",at[i]);
110 17 : log.printf("\n");
111 17 : log.printf(" with an offset");
112 85 : for(unsigned i=0; i<offset.size(); i++) log.printf(" %f",offset[i]);
113 17 : log.printf("\n");
114 17 : log.printf(" with force constant");
115 85 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
116 17 : log.printf("\n");
117 17 : log.printf(" and exponent");
118 85 : for(unsigned i=0; i<exp.size(); i++) log.printf(" %f",exp[i]);
119 17 : log.printf("\n");
120 17 : log.printf(" rescaled");
121 85 : for(unsigned i=0; i<eps.size(); i++) log.printf(" %f",eps[i]);
122 17 : log.printf("\n");
123 :
124 51 : addComponent("force2"); componentIsNotPeriodic("force2");
125 17 : }
126 :
127 28015 : void UWalls::calculate() {
128 : double ene=0.0;
129 : double totf2=0.0;
130 84045 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
131 : double f = 0.0;
132 28015 : const double cv=difference(i,at[i],getArgument(i));
133 28015 : const double off=offset[i];
134 28015 : const double epsilon=eps[i];
135 28015 : const double uscale = (cv+off)/epsilon;
136 28015 : if( uscale > 0.) {
137 6 : const double k=kappa[i];
138 6 : const double exponent=exp[i];
139 6 : double power = pow( uscale, exponent );
140 6 : f = -( k / epsilon ) * exponent * power / uscale;
141 6 : ene += k * power;
142 6 : totf2 += f * f;
143 : }
144 : setOutputForce(i,f);
145 : }
146 : setBias(ene);
147 56030 : getPntrToComponent("force2")->set(totf2);
148 28015 : }
149 :
150 : }
151 5517 : }
|
