Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 3553258 : Value::Value():
34 : action(NULL),
35 : value_set(false),
36 : value(0.0),
37 : inputForce(0.0),
38 : hasForce(false),
39 : hasDeriv(true),
40 : periodicity(unset),
41 : min(0.0),
42 : max(0.0),
43 : max_minus_min(0.0),
44 10659774 : inv_max_minus_min(0.0)
45 : {
46 3553258 : }
47 :
48 47797 : Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv):
49 : action(av),
50 : value_set(false),
51 : value(0.0),
52 : inputForce(0.0),
53 : hasForce(false),
54 : name(name),
55 : hasDeriv(withderiv),
56 : periodicity(unset),
57 : min(0.0),
58 : max(0.0),
59 : max_minus_min(0.0),
60 143391 : inv_max_minus_min(0.0)
61 : {
62 47797 : }
63 :
64 1175690 : void Value::setupPeriodicity() {
65 1175690 : if( min==0 && max==0 ) {
66 31454 : periodicity=notperiodic;
67 : } else {
68 1144236 : periodicity=periodic;
69 1144236 : max_minus_min=max-min;
70 1144236 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
71 1144236 : inv_max_minus_min=1.0/max_minus_min;
72 : }
73 1175690 : }
74 :
75 295394 : bool Value::isPeriodic()const {
76 295394 : plumed_massert(periodicity!=unset,"periodicity should be set");
77 295394 : return periodicity==periodic;
78 : }
79 :
80 612319 : bool Value::applyForce(std::vector<double>& forces ) const {
81 612319 : if( !hasForce ) return false;
82 : plumed_dbg_massert( derivatives.size()==forces.size()," forces array has wrong size" );
83 105607 : const unsigned N=derivatives.size();
84 92583969 : for(unsigned i=0; i<N; ++i) forces[i]=inputForce*derivatives[i];
85 : return true;
86 : }
87 :
88 2522320 : void Value::setNotPeriodic() {
89 2522320 : min=0; max=0; periodicity=notperiodic;
90 2522320 : }
91 :
92 1144235 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
93 1144235 : str_min=pmin;
94 1144233 : if( !Tools::convert(str_min,min) ) action->error("could not convert period string " + str_min + " to real");
95 1144242 : str_max=pmax;
96 1144239 : if( !Tools::convert(str_max,max) ) action->error("could not convert period string " + str_max + " to read");
97 1144240 : setupPeriodicity();
98 1144236 : }
99 :
100 25988 : void Value::getDomain(std::string&minout,std::string&maxout) const {
101 25988 : plumed_massert(periodicity==periodic,"function should be periodic");
102 25988 : minout=str_min;
103 25988 : maxout=str_max;
104 25988 : }
105 :
106 22 : void Value::getDomain(double&minout,double&maxout) const {
107 22 : plumed_massert(periodicity==periodic,"function should be periodic");
108 22 : minout=min;
109 22 : maxout=max;
110 22 : }
111 :
112 594 : void Value::setGradients() {
113 : // Can't do gradients if we don't have derivatives
114 1188 : if( !hasDeriv ) return;
115 : gradients.clear();
116 402 : ActionAtomistic*aa=dynamic_cast<ActionAtomistic*>(action);
117 402 : ActionWithArguments*aw=dynamic_cast<ActionWithArguments*>(action);
118 402 : if(aa) {
119 210 : Atoms&atoms((aa->plumed).getAtoms());
120 39972 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
121 39552 : AtomNumber an=aa->getAbsoluteIndex(j);
122 19776 : if(atoms.isVirtualAtom(an)) {
123 : const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
124 36 : for(const auto & p : a->getGradients()) {
125 : // controllare l'ordine del matmul:
126 120 : gradients[p.first]+=matmul(Vector(derivatives[3*j],derivatives[3*j+1],derivatives[3*j+2]),p.second);
127 : }
128 : } else {
129 118638 : for(unsigned i=0; i<3; i++) gradients[an][i]+=derivatives[3*j+i];
130 : }
131 : }
132 192 : } else if(aw) {
133 192 : std::vector<Value*> values=aw->getArguments();
134 960 : for(unsigned j=0; j<derivatives.size(); j++) {
135 21048 : for(const auto & p : values[j]->gradients) {
136 19896 : AtomNumber iatom=p.first;
137 19896 : gradients[iatom]+=p.second*derivatives[j];
138 : }
139 : }
140 0 : } else plumed_error();
141 : }
142 :
143 261 : double Value::projection(const Value& v1,const Value&v2) {
144 : double proj=0.0;
145 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
146 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
147 16428 : for(const auto & p1 : grad1) {
148 15906 : AtomNumber a=p1.first;
149 : const auto p2=grad2.find(a);
150 15906 : if(p2!=grad2.end()) {
151 8224 : proj+=dotProduct(p1.second,(*p2).second);
152 : }
153 : }
154 261 : return proj;
155 : }
156 :
157 4474 : ActionWithValue* Value::getPntrToAction() {
158 4474 : plumed_assert( action!=NULL );
159 4474 : return action;
160 : }
161 :
162 0 : void copy( const Value& val1, Value& val2 ) {
163 0 : unsigned nder=val1.getNumberOfDerivatives();
164 0 : if( nder!=val2.getNumberOfDerivatives() ) { val2.resizeDerivatives( nder ); }
165 : val2.clearDerivatives();
166 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2.addDerivative( i, val1.getDerivative(i) );
167 : val2.set( val1.get() );
168 0 : }
169 :
170 0 : void copy( const Value& val1, Value* val2 ) {
171 0 : unsigned nder=val1.getNumberOfDerivatives();
172 0 : if( nder!=val2->getNumberOfDerivatives() ) { val2->resizeDerivatives( nder ); }
173 : val2->clearDerivatives();
174 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2->addDerivative( i, val1.getDerivative(i) );
175 : val2->set( val1.get() );
176 0 : }
177 :
178 0 : void add( const Value& val1, Value* val2 ) {
179 0 : plumed_assert( val1.getNumberOfDerivatives()==val2->getNumberOfDerivatives() );
180 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2->addDerivative( i, val1.getDerivative(i) );
181 0 : val2->set( val1.get() + val2->get() );
182 0 : }
183 :
184 5874 : }
|
