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 "core/ActionShortcut.h"
23 : #include "core/ActionWithArguments.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/ActionRegister.h"
27 : #include "tools/Tools.h"
28 :
29 : namespace PLMD {
30 : namespace bias {
31 :
32 : class Walls : public ActionShortcut {
33 : public:
34 : explicit Walls(const ActionOptions&);
35 : static void registerKeywords(Keywords& keys);
36 : };
37 :
38 : PLUMED_REGISTER_ACTION(Walls,"UPPER_WALLS")
39 : PLUMED_REGISTER_ACTION(Walls,"LOWER_WALLS")
40 :
41 112 : void Walls::registerKeywords(Keywords& keys) {
42 112 : ActionShortcut::registerKeywords(keys);
43 224 : keys.addInputKeyword("numbered","ARG","scalar/vector","the arguments on which the bias is acting");
44 112 : keys.add("compulsory","AT","the positions of the wall. The a_i in the expression for a wall.");
45 112 : keys.add("compulsory","KAPPA","the force constant for the wall. The k_i in the expression for a wall.");
46 112 : keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall. The o_i in the expression for a wall.");
47 112 : keys.add("compulsory","EXP","2.0","the powers for the walls. The e_i in the expression for a wall.");
48 112 : keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
49 112 : keys.add("hidden","STRIDE","the frequency with which the forces due to the bias should be calculated. This can be used to correctly set up multistep algorithms");
50 224 : keys.addOutputComponent("bias","default","scalar","the instantaneous value of the bias potential");
51 224 : keys.addOutputComponent("force2","default","scalar","the instantaneous value of the squared force due to this bias potential");
52 112 : keys.addActionNameSuffix("_SCALAR");
53 112 : keys.needsAction("COMBINE");
54 112 : keys.needsAction("CUSTOM");
55 112 : keys.needsAction("SUM");
56 112 : keys.needsAction("BIASVALUE");
57 112 : }
58 :
59 65 : Walls::Walls(const ActionOptions&ao):
60 : Action(ao),
61 65 : ActionShortcut(ao) {
62 : // Read the arguments
63 65 : std::vector<std::string> args;
64 130 : parseVector("ARG",args);
65 65 : if( args.size()==0 ) {
66 0 : error("found no input arguments");
67 : }
68 65 : std::string allargs=args[0];
69 69 : for(unsigned i=1; i<args.size(); ++i) {
70 8 : allargs += "," + args[i];
71 : }
72 65 : std::vector<Value*> vals;
73 65 : ActionWithArguments::interpretArgumentList( args, plumed.getActionSet(), this, vals );
74 65 : if( vals.size()==0 ) {
75 0 : error("found no input arguments");
76 : }
77 :
78 : // Find the rank
79 65 : unsigned rank=vals[0]->getRank();
80 134 : for(unsigned i=0; i<vals.size(); ++i) {
81 69 : if( vals[i]->getRank()>0 && vals[i]->hasDerivatives() ) {
82 0 : error("argument should not be function on grid");
83 : }
84 69 : if( vals[i]->getRank()!=rank ) {
85 0 : error("all arguments should have same rank");
86 : }
87 : }
88 65 : if( rank==0 ) {
89 55 : if( getName()=="UPPER_WALLS") {
90 144 : readInputLine( getShortcutLabel() + ": UPPER_WALLS_SCALAR ARG=" + allargs + " " + convertInputLineToString() );
91 7 : } else if( getName()=="LOWER_WALLS") {
92 21 : readInputLine( getShortcutLabel() + ": LOWER_WALLS_SCALAR ARG=" + allargs + " " + convertInputLineToString() );
93 : } else {
94 0 : plumed_merror( getName() + " is not valid");
95 : }
96 : return;
97 : }
98 :
99 : // Note : the sizes of these vectors are checked automatically by parseVector
100 10 : std::vector<std::string> kappa(args.size());
101 20 : parseVector("KAPPA",kappa);
102 10 : std::vector<std::string> offset(kappa.size());
103 20 : parseVector("OFFSET",offset);
104 10 : std::vector<std::string> eps(kappa.size());
105 20 : parseVector("EPS",eps);
106 10 : std::vector<std::string> exp(kappa.size());
107 20 : parseVector("EXP",exp);
108 10 : std::vector<std::string> at(kappa.size());
109 20 : parseVector("AT",at);
110 :
111 : std::string stride;
112 20 : parse("STRIDE",stride);
113 10 : if( stride.length()>0 ) {
114 0 : stride=" STRIDE=" + stride;
115 : }
116 : std::string biasinp, forceinp;
117 20 : for(unsigned i=0; i<args.size(); ++i) {
118 : std::string argn;
119 10 : std::size_t dot=args[i].find_first_of(".");
120 10 : if(dot!=std::string::npos) {
121 0 : argn = args[i].substr(0,dot) + "_" + args[i].substr(dot+1);
122 : } else {
123 : argn = args[i];
124 : }
125 : // The wall acts on the value of the argument in its domain, without periodic images
126 : // (see https://github.com/plumed/plumed2/issues/658). COMBINE cannot be used here
127 : // because it applies periodic boundary conditions to the difference: CUSTOM with PERIODIC=NO does not.
128 10 : std::vector<Value*> argvals;
129 10 : ActionWithArguments::interpretArgumentList( std::vector<std::string>(1,args[i]), plumed.getActionSet(), this, argvals );
130 20 : for(const auto & v : argvals) {
131 10 : if( v->isPeriodic() ) {
132 : double min, max, atv;
133 4 : v->getDomain( min, max );
134 4 : Tools::convert( at[i], atv );
135 4 : if( atv<min || atv>max ) {
136 0 : error("the position of the wall (AT) for periodic argument " + v->getName() + " should be within its periodic domain");
137 : }
138 4 : log.printf(" WARNING: argument %s is periodic. The wall acts on its value in the domain [%f,%f]\n"
139 : " without periodic images, so the bias is discontinuous at the boundary of the domain.\n"
140 : " Make sure that the walls keep the argument far from it (e.g. use both an upper and a lower wall).\n",
141 4 : v->getName().c_str(), min, max );
142 : }
143 : }
144 : // The scaled distance from the wall, s=(x-AT+OFFSET)/EPS (UPPER_WALLS) or s=(x-AT-OFFSET)/EPS (LOWER_WALLS).
145 : // The wall acts where s is positive (UPPER_WALLS) or negative (LOWER_WALLS): select() is used instead of
146 : // a product with step(), because a non-integer power of a negative number is nan and nan*0 is nan.
147 : // Energy and squared force of each element are computed by one CUSTOM each: the chain of one CUSTOM per
148 : // operation that was used before took about four times longer.
149 : std::string scaled, inactive, power;
150 10 : if( getName()=="UPPER_WALLS" ) {
151 28 : scaled = "(((x-(" + at[i] + "))+(" + offset[i] + "))/(" + eps[i] + "))";
152 14 : inactive = "step(-" + scaled + ")";
153 21 : power = scaled + "^(" + exp[i] + ")";
154 3 : } else if( getName()=="LOWER_WALLS" ) {
155 12 : scaled = "(((x-(" + at[i] + "))-(" + offset[i] + "))/(" + eps[i] + "))";
156 6 : inactive = "step(" + scaled + ")";
157 12 : power = "(-" + scaled + ")^(" + exp[i] + ")";
158 : }
159 50 : readInputLine( getShortcutLabel() + "_v_wall_" + argn + ": CUSTOM PERIODIC=NO FUNC=select(" + inactive + ",0,(" + kappa[i] + ")*" + power + ") ARG=" + args[i] );
160 50 : readInputLine( getShortcutLabel() + "_wall_" + argn + ": SUM ARG=" + getShortcutLabel() + "_v_wall_" + argn + " PERIODIC=NO");
161 : // the derivative of the wall with respect to x is KAPPA*EXP*s^EXP/(s*EPS)
162 80 : readInputLine( getShortcutLabel() + "_v_force2_" + argn + ": CUSTOM PERIODIC=NO FUNC=select(" + inactive + ",0,((" + kappa[i] + ")*(" + exp[i] + ")*" + power + "/(" + scaled + "*(" + eps[i] + ")))^2) ARG=" + args[i] );
163 40 : readInputLine( getShortcutLabel() + "_force2_" + argn + ": SUM ARG=" + getShortcutLabel() + "_v_force2_" + argn + " PERIODIC=NO");
164 10 : if(i==0) {
165 20 : biasinp = " ARG=" + getShortcutLabel() + "_wall_" + argn;
166 30 : forceinp = " ARG=" + getShortcutLabel() + "_force2_" + argn;
167 : } else {
168 0 : biasinp += "," + getShortcutLabel() + "_wall_" + argn;
169 0 : forceinp += "," + getShortcutLabel() + "_force2_" + argn;
170 : }
171 : }
172 10 : readInputLine( getShortcutLabel() + "_bias: COMBINE PERIODIC=NO " + biasinp );
173 20 : readInputLine( "BIASVALUE ARG=" + getShortcutLabel() + "_bias" + stride );
174 20 : readInputLine( getShortcutLabel() + "_force2: COMBINE PERIODIC=NO " + forceinp );
175 85 : }
176 :
177 : }
178 : }
|
