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 "Value.h"
23 : #include "ActionWithValue.h"
24 : #include "ActionAtomistic.h"
25 : #include "ActionWithArguments.h"
26 : #include "ActionWithVirtualAtom.h"
27 : #include "tools/Exception.h"
28 : #include "Atoms.h"
29 : #include "PlumedMain.h"
30 :
31 : namespace PLMD {
32 :
33 3561022 : Value::Value():
34 3561022 : action(NULL),
35 3561022 : value_set(false),
36 3561022 : value(0.0),
37 3561022 : inputForce(0.0),
38 3561022 : hasForce(false),
39 3561022 : hasDeriv(true),
40 3561022 : periodicity(unset),
41 3561022 : min(0.0),
42 3561022 : max(0.0),
43 3561022 : max_minus_min(0.0),
44 3561022 : inv_max_minus_min(0.0) {
45 3561022 : }
46 :
47 18 : Value::Value(const std::string& name):
48 18 : action(NULL),
49 18 : value_set(false),
50 18 : value(0.0),
51 18 : inputForce(0.0),
52 18 : hasForce(false),
53 18 : name(name),
54 18 : hasDeriv(true),
55 18 : periodicity(unset),
56 18 : min(0.0),
57 18 : max(0.0),
58 18 : max_minus_min(0.0),
59 18 : inv_max_minus_min(0.0) {
60 18 : }
61 :
62 49117 : Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv):
63 49117 : action(av),
64 49117 : value_set(false),
65 49117 : value(0.0),
66 49117 : inputForce(0.0),
67 49117 : hasForce(false),
68 49117 : name(name),
69 49117 : hasDeriv(withderiv),
70 49117 : periodicity(unset),
71 49117 : min(0.0),
72 49117 : max(0.0),
73 49117 : max_minus_min(0.0),
74 49117 : inv_max_minus_min(0.0) {
75 49117 : }
76 :
77 1184097 : void Value::setupPeriodicity() {
78 1184097 : if( min==0 && max==0 ) {
79 32591 : periodicity=notperiodic;
80 : } else {
81 1151506 : periodicity=periodic;
82 1151506 : max_minus_min=max-min;
83 1151506 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
84 1151506 : inv_max_minus_min=1.0/max_minus_min;
85 : }
86 1184097 : }
87 :
88 1674137 : bool Value::isPeriodic()const {
89 1674137 : plumed_massert(periodicity!=unset,"periodicity should be set");
90 1674137 : return periodicity==periodic;
91 : }
92 :
93 988341 : bool Value::applyForce(std::vector<double>& forces ) const {
94 988341 : if( !hasForce ) {
95 : return false;
96 : }
97 : plumed_dbg_massert( derivatives.size()==forces.size()," forces array has wrong size" );
98 120729 : const unsigned N=derivatives.size();
99 50118218 : for(unsigned i=0; i<N; ++i) {
100 49997489 : forces[i]=inputForce*derivatives[i];
101 : }
102 : return true;
103 : }
104 :
105 3529072 : void Value::setNotPeriodic() {
106 3529072 : min=0;
107 3529072 : max=0;
108 3529072 : periodicity=notperiodic;
109 3529072 : }
110 :
111 1151506 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
112 1151506 : str_min=pmin;
113 1151506 : if( !Tools::convertNoexcept(str_min,min) ) {
114 0 : action->error("could not convert period string " + str_min + " to real");
115 : }
116 1151506 : str_max=pmax;
117 1151506 : if( !Tools::convertNoexcept(str_max,max) ) {
118 0 : action->error("could not convert period string " + str_max + " to read");
119 : }
120 1151506 : setupPeriodicity();
121 1151506 : }
122 :
123 27988 : void Value::getDomain(std::string&minout,std::string&maxout) const {
124 27988 : plumed_massert(periodicity==periodic,"function should be periodic");
125 27988 : minout=str_min;
126 27988 : maxout=str_max;
127 27988 : }
128 :
129 34 : void Value::getDomain(double&minout,double&maxout) const {
130 34 : plumed_massert(periodicity==periodic,"function should be periodic");
131 34 : minout=min;
132 34 : maxout=max;
133 34 : }
134 :
135 183 : void Value::setGradients() {
136 : // Can't do gradients if we don't have derivatives
137 183 : if( !hasDeriv ) {
138 : return;
139 : }
140 : gradients.clear();
141 183 : ActionAtomistic*aa=dynamic_cast<ActionAtomistic*>(action);
142 183 : ActionWithArguments*aw=dynamic_cast<ActionWithArguments*>(action);
143 183 : if(aa) {
144 100 : const Atoms&atoms((aa->plumed).getAtoms());
145 8092 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
146 7992 : AtomNumber an=aa->getAbsoluteIndex(j);
147 7992 : if(atoms.isVirtualAtom(an)) {
148 : const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
149 33 : for(const auto & p : a->getGradients()) {
150 : // controllare l'ordine del matmul:
151 30 : gradients[p.first]+=matmul(Vector(derivatives[3*j],derivatives[3*j+1],derivatives[3*j+2]),p.second);
152 : }
153 : } else {
154 31956 : for(unsigned i=0; i<3; i++) {
155 23967 : gradients[an][i]+=derivatives[3*j+i];
156 : }
157 : }
158 : }
159 83 : } else if(aw) {
160 83 : std::vector<Value*> values=aw->getArguments();
161 249 : for(unsigned j=0; j<derivatives.size(); j++) {
162 8174 : for(const auto & p : values[j]->gradients) {
163 8008 : AtomNumber iatom=p.first;
164 8008 : gradients[iatom]+=p.second*derivatives[j];
165 : }
166 : }
167 : } else {
168 0 : plumed_error();
169 : }
170 : }
171 :
172 261 : double Value::projection(const Value& v1,const Value&v2) {
173 : double proj=0.0;
174 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
175 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
176 16167 : for(const auto & p1 : grad1) {
177 15906 : AtomNumber a=p1.first;
178 : const auto p2=grad2.find(a);
179 15906 : if(p2!=grad2.end()) {
180 8224 : proj+=dotProduct(p1.second,(*p2).second);
181 : }
182 : }
183 261 : return proj;
184 : }
185 :
186 4926 : ActionWithValue* Value::getPntrToAction() {
187 4926 : plumed_assert( action!=NULL );
188 4926 : return action;
189 : }
190 :
191 0 : void copy( const Value& val1, Value& val2 ) {
192 0 : unsigned nder=val1.getNumberOfDerivatives();
193 0 : if( nder!=val2.getNumberOfDerivatives() ) {
194 0 : val2.resizeDerivatives( nder );
195 : }
196 : val2.clearDerivatives();
197 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
198 : val2.addDerivative( i, val1.getDerivative(i) );
199 : }
200 : val2.set( val1.get() );
201 0 : }
202 :
203 0 : void copy( const Value& val1, Value* val2 ) {
204 0 : unsigned nder=val1.getNumberOfDerivatives();
205 0 : if( nder!=val2->getNumberOfDerivatives() ) {
206 0 : val2->resizeDerivatives( nder );
207 : }
208 : val2->clearDerivatives();
209 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
210 : val2->addDerivative( i, val1.getDerivative(i) );
211 : }
212 : val2->set( val1.get() );
213 0 : }
214 :
215 0 : void add( const Value& val1, Value* val2 ) {
216 0 : plumed_assert( val1.getNumberOfDerivatives()==val2->getNumberOfDerivatives() );
217 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
218 : val2->addDerivative( i, val1.getDerivative(i) );
219 : }
220 0 : val2->set( val1.get() + val2->get() );
221 0 : }
222 :
223 : }
|
