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 1135539 : 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 3406617 : inv_max_minus_min(0.0)
45 : {
46 1135539 : }
47 :
48 46671 : 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 186684 : inv_max_minus_min(0.0)
61 : {
62 46671 : }
63 :
64 1127152 : void Value::setupPeriodicity() {
65 1127152 : if( min==0 && max==0 ) {
66 30595 : periodicity=notperiodic;
67 : } else {
68 1096557 : periodicity=periodic;
69 1096557 : max_minus_min=max-min;
70 1096557 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
71 1096557 : inv_max_minus_min=1.0/max_minus_min;
72 : }
73 1127152 : }
74 :
75 145612 : bool Value::isPeriodic()const {
76 145612 : plumed_massert(periodicity!=unset,"periodicity should be set");
77 145612 : return periodicity==periodic;
78 : }
79 :
80 419938 : bool Value::applyForce(std::vector<double>& forces ) const {
81 419938 : if( !hasForce ) return false;
82 : plumed_dbg_massert( derivatives.size()==forces.size()," forces array has wrong size" );
83 81272 : const unsigned N=derivatives.size();
84 133373711 : for(unsigned i=0; i<N; ++i) forces[i]=inputForce*derivatives[i];
85 : return true;
86 : }
87 :
88 61309 : void Value::setNotPeriodic() {
89 61309 : min=0; max=0; periodicity=notperiodic;
90 61309 : }
91 :
92 1096557 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
93 1096557 : str_min=pmin;
94 1096557 : if( !Tools::convert(str_min,min) ) action->error("could not convert period string " + str_min + " to real");
95 1096557 : str_max=pmax;
96 1096557 : if( !Tools::convert(str_max,max) ) action->error("could not convert period string " + str_max + " to read");
97 1096557 : setupPeriodicity();
98 1096557 : }
99 :
100 18067 : void Value::getDomain(std::string&minout,std::string&maxout) const {
101 18067 : plumed_massert(periodicity==periodic,"function should be periodic");
102 18067 : minout=str_min;
103 18067 : maxout=str_max;
104 18067 : }
105 :
106 12 : void Value::getDomain(double&minout,double&maxout) const {
107 12 : plumed_massert(periodicity==periodic,"function should be periodic");
108 12 : minout=min;
109 12 : maxout=max;
110 12 : }
111 :
112 266 : void Value::setGradients() {
113 266 : gradients.clear();
114 266 : ActionAtomistic*aa=dynamic_cast<ActionAtomistic*>(action);
115 266 : ActionWithArguments*aw=dynamic_cast<ActionWithArguments*>(action);
116 266 : if(aa) {
117 100 : Atoms&atoms((aa->plumed).getAtoms());
118 16084 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
119 15984 : AtomNumber an=aa->getAbsoluteIndex(j);
120 7992 : if(atoms.isVirtualAtom(an)) {
121 : const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
122 33 : for(const auto & p : a->getGradients()) {
123 : // controllare l'ordine del matmul:
124 150 : gradients[p.first]+=matmul(Vector(derivatives[3*j],derivatives[3*j+1],derivatives[3*j+2]),p.second);
125 : }
126 : } else {
127 55923 : for(unsigned i=0; i<3; i++) gradients[an][i]+=derivatives[3*j+i];
128 : }
129 : }
130 166 : } else if(aw) {
131 166 : std::vector<Value*> values=aw->getArguments();
132 830 : for(unsigned j=0; j<derivatives.size(); j++) {
133 8340 : for(const auto & p : values[j]->gradients) {
134 8008 : AtomNumber iatom=p.first;
135 8008 : gradients[iatom]+=p.second*derivatives[j];
136 : }
137 : }
138 0 : } else plumed_error();
139 266 : }
140 :
141 261 : double Value::projection(const Value& v1,const Value&v2) {
142 : double proj=0.0;
143 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
144 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
145 16167 : for(const auto & p1 : grad1) {
146 15906 : AtomNumber a=p1.first;
147 : const auto p2=grad2.find(a);
148 15906 : if(p2!=grad2.end()) {
149 8224 : proj+=dotProduct(p1.second,(*p2).second);
150 : }
151 : }
152 261 : return proj;
153 : }
154 :
155 4197 : ActionWithValue* Value::getPntrToAction() {
156 4197 : plumed_assert( action!=NULL );
157 4197 : return action;
158 : }
159 :
160 0 : void copy( const Value& val1, Value& val2 ) {
161 0 : unsigned nder=val1.getNumberOfDerivatives();
162 0 : if( nder!=val2.getNumberOfDerivatives() ) { val2.resizeDerivatives( nder ); }
163 : val2.clearDerivatives();
164 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2.addDerivative( i, val1.getDerivative(i) );
165 : val2.set( val1.get() );
166 0 : }
167 :
168 0 : void copy( const Value& val1, Value* val2 ) {
169 0 : unsigned nder=val1.getNumberOfDerivatives();
170 0 : if( nder!=val2->getNumberOfDerivatives() ) { val2->resizeDerivatives( nder ); }
171 : val2->clearDerivatives();
172 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2->addDerivative( i, val1.getDerivative(i) );
173 : val2->set( val1.get() );
174 0 : }
175 :
176 0 : void add( const Value& val1, Value* val2 ) {
177 0 : plumed_assert( val1.getNumberOfDerivatives()==val2->getNumberOfDerivatives() );
178 0 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) val2->addDerivative( i, val1.getDerivative(i) );
179 0 : val2->set( val1.get() + val2->get() );
180 0 : }
181 :
182 4839 : }
|
