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