Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2017 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/ActionWithVector.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/LeptonCall.h"
25 : #include "tools/Angle.h"
26 :
27 : namespace PLMD {
28 : namespace symfunc {
29 :
30 : //+PLUMEDOC COLVAR GSYMFUNC_THREEBODY
31 : /*
32 : Calculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom
33 :
34 : \par Examples
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : class ThreeBodyGFunctions : public ActionWithVector {
40 : private:
41 : std::vector<LeptonCall> functions;
42 : public:
43 : static void registerKeywords( Keywords& keys );
44 : explicit ThreeBodyGFunctions(const ActionOptions&);
45 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
46 : void calculate() override ;
47 : unsigned getNumberOfDerivatives() override;
48 : void performTask( const unsigned& task_index, MultiValue& myvals ) const override ;
49 : };
50 :
51 : PLUMED_REGISTER_ACTION(ThreeBodyGFunctions,"GSYMFUNC_THREEBODY")
52 :
53 11 : void ThreeBodyGFunctions::registerKeywords( Keywords& keys ) {
54 11 : ActionWithVector::registerKeywords( keys );
55 11 : keys.use("ARG");
56 22 : keys.add("compulsory","WEIGHT","the matrix that contains the weights that should be used for each connection");
57 22 : keys.add("numbered","FUNCTION","the parameters of the function you would like to compute");
58 11 : ActionWithValue::useCustomisableComponents( keys );
59 11 : }
60 :
61 6 : ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao):
62 : Action(ao),
63 6 : ActionWithVector(ao) {
64 6 : if( getNumberOfArguments()!=3 ) {
65 0 : error("found wrong number of arguments in input");
66 : }
67 : std::vector<Value*> wval;
68 12 : parseArgumentList("WEIGHT",wval);
69 6 : if( wval.size()!=1 ) {
70 0 : error("keyword WEIGHT should be provided with the label of a single action");
71 : }
72 :
73 24 : for(unsigned i=0; i<3; ++i) {
74 18 : if( getPntrToArgument(i)->getRank()!=2 ) {
75 0 : error("input argument should be a matrix");
76 : }
77 18 : if( wval[0]->getShape()[0]!=getPntrToArgument(i)->getShape()[0] || wval[0]->getShape()[1]!=getPntrToArgument(i)->getShape()[1] ) {
78 0 : error("mismatched shapes of matrices in input");
79 : }
80 : }
81 6 : log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() );
82 : // Rerequest the arguments
83 6 : std::vector<Value*> myargs( getArguments() );
84 6 : myargs.push_back( wval[0] );
85 6 : requestArguments( myargs );
86 30 : for(unsigned i=0; i<myargs.size(); ++i) {
87 24 : myargs[i]->buildDataStore();
88 : }
89 6 : std::vector<unsigned> shape(1);
90 6 : shape[0] = getPntrToArgument(0)->getShape()[0];
91 :
92 : // And now read the functions to compute
93 6 : for(int i=1;; i++) {
94 : std::string myfunc, mystr, lab, num;
95 17 : Tools::convert(i,num);
96 34 : if( !parseNumbered("FUNCTION",i,mystr ) ) {
97 : break;
98 : }
99 11 : std::vector<std::string> data=Tools::getWords(mystr);
100 22 : if( !Tools::parse(data,"LABEL",lab ) ) {
101 0 : error("found no LABEL in FUNCTION" + num + " specification");
102 : }
103 11 : addComponent( lab, shape );
104 11 : componentIsNotPeriodic( lab );
105 22 : if( !Tools::parse(data,"FUNC",myfunc) ) {
106 0 : error("found no FUNC in FUNCTION" + num + " specification");
107 : }
108 11 : log.printf(" component labelled %s is computed using %s \n",lab.c_str(), myfunc.c_str() );
109 11 : functions.push_back( LeptonCall() );
110 11 : std::vector<std::string> argnames(1);
111 : argnames[0]="ajik";
112 11 : if( myfunc.find("rij")!=std::string::npos ) {
113 8 : argnames.push_back("rij");
114 : }
115 11 : if( myfunc.find("rik")!=std::string::npos ) {
116 4 : if( argnames.size()<2 ) {
117 0 : error("if you have a function of rik it must also be a function of rij -- email gareth.tribello@gmail.com if this is a problem");
118 : }
119 8 : argnames.push_back("rik");
120 : }
121 11 : if( myfunc.find("rjk")!=std::string::npos ) {
122 3 : if( argnames.size()<2 ) {
123 0 : error("if you have a function of rjk it must also be a function of rij and rik -- email gareth.tribello@gmail.com if this is a problem");
124 : }
125 6 : argnames.push_back("rjk");
126 : }
127 11 : functions[i-1].set( myfunc, argnames, this, true );
128 22 : }
129 6 : checkRead();
130 6 : }
131 :
132 2 : std::string ThreeBodyGFunctions::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
133 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
134 6 : if( getConstPntrToComponent(i)->getName().find(cname)!=std::string::npos ) {
135 : std::string num;
136 2 : Tools::convert( i+1, num );
137 4 : return "the function defined by the FUNCTION" + num + " keyword";
138 : }
139 : }
140 0 : plumed_error();
141 : return "";
142 : }
143 :
144 36 : unsigned ThreeBodyGFunctions::getNumberOfDerivatives() {
145 36 : return 0;
146 : }
147 :
148 19 : void ThreeBodyGFunctions::calculate() {
149 19 : runAllTasks();
150 19 : }
151 :
152 472 : void ThreeBodyGFunctions::performTask( const unsigned& task_index, MultiValue& myvals ) const {
153 : const Value* wval = getPntrToArgument(3);
154 : const Value* xval = getPntrToArgument(0);
155 : const Value* yval = getPntrToArgument(1);
156 : const Value* zval = getPntrToArgument(2);
157 : Angle angle;
158 472 : Vector disti, distj;
159 472 : unsigned matsize = wval->getNumberOfValues();
160 472 : std::vector<double> values(4);
161 472 : std::vector<Vector> der_i(4), der_j(4);
162 472 : unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
163 13880 : for(unsigned i=0; i<nbonds; ++i) {
164 13408 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
165 13408 : double weighti = wval->get( ipos );
166 13408 : if( weighti<epsilon ) {
167 6584 : continue ;
168 : }
169 6824 : disti[0] = xval->get( ipos );
170 6824 : disti[1] = yval->get( ipos );
171 6824 : disti[2] = zval->get( ipos );
172 6824 : values[1] = disti.modulo2();
173 6824 : der_i[1]=2*disti;
174 6824 : der_i[2].zero();
175 189327 : for(unsigned j=0; j<i; ++j) {
176 182503 : unsigned jpos = ncols*task_index + wval->getRowIndex( task_index, j );
177 182503 : double weightj = wval->get( jpos );
178 182503 : if( weightj<epsilon ) {
179 45494 : continue ;
180 : }
181 137009 : distj[0] = xval->get( jpos );
182 137009 : distj[1] = yval->get( jpos );
183 137009 : distj[2] = zval->get( jpos );
184 137009 : values[2] = distj.modulo2();
185 137009 : der_j[1].zero();
186 137009 : der_j[2]=2*distj;
187 137009 : der_i[3] = ( disti - distj );
188 137009 : values[3] = der_i[3].modulo2();
189 137009 : der_i[3] = 2*der_i[3];
190 137009 : der_j[3] = -der_i[3];
191 : // Compute angle between bonds
192 137009 : values[0] = angle.compute( disti, distj, der_i[0], der_j[0] );
193 : // Compute product of weights
194 137009 : double weightij = weighti*weightj;
195 : // Now compute all symmetry functions
196 414532 : for(unsigned n=0; n<functions.size(); ++n) {
197 277523 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
198 277523 : double nonweight = functions[n].evaluate( values );
199 277523 : myvals.addValue( ostrn, nonweight*weightij );
200 277523 : if( doNotCalculateDerivatives() ) {
201 6843 : continue;
202 : }
203 :
204 567230 : for(unsigned m=0; m<functions[n].getNumberOfArguments(); ++m) {
205 296550 : double der = weightij*functions[n].evaluateDeriv( m, values );
206 296550 : myvals.addDerivative( ostrn, ipos, der*der_i[m][0] );
207 296550 : myvals.addDerivative( ostrn, matsize+ipos, der*der_i[m][1] );
208 296550 : myvals.addDerivative( ostrn, 2*matsize+ipos, der*der_i[m][2] );
209 296550 : myvals.addDerivative( ostrn, jpos, der*der_j[m][0] );
210 296550 : myvals.addDerivative( ostrn, matsize+jpos, der*der_j[m][1] );
211 296550 : myvals.addDerivative( ostrn, 2*matsize+jpos, der*der_j[m][2] );
212 : }
213 270680 : myvals.addDerivative( ostrn, 3*matsize+ipos, nonweight*weightj );
214 270680 : myvals.addDerivative( ostrn, 3*matsize+jpos, nonweight*weighti );
215 : }
216 : }
217 : }
218 472 : if( doNotCalculateDerivatives() ) {
219 : return ;
220 : }
221 :
222 : // And update the elements that have derivatives
223 : // Needs a separate loop here as there may be forces from j
224 8320 : for(unsigned i=0; i<nbonds; ++i) {
225 8192 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
226 8192 : double weighti = wval->get( ipos );
227 8192 : if( weighti<epsilon ) {
228 3322 : continue ;
229 : }
230 :
231 16286 : for(unsigned n=0; n<functions.size(); ++n) {
232 11416 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
233 11416 : myvals.updateIndex( ostrn, ipos );
234 11416 : myvals.updateIndex( ostrn, matsize+ipos );
235 11416 : myvals.updateIndex( ostrn, 2*matsize+ipos );
236 11416 : myvals.updateIndex( ostrn, 3*matsize+ipos );
237 : }
238 : }
239 : }
240 :
241 : }
242 : }
|
