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 22 : keys.addInputKeyword("compulsory","ARG","matrix","three matrices containing the bond vectors of interest");
56 22 : keys.addInputKeyword("compulsory","WEIGHT","matrix","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 : {
65 6 : if( getNumberOfArguments()!=3 ) error("found wrong number of arguments in input");
66 12 : std::vector<Value*> wval; parseArgumentList("WEIGHT",wval);
67 6 : if( wval.size()!=1 ) error("keyword WEIGHT should be provided with the label of a single action");
68 :
69 24 : for(unsigned i=0; i<3; ++i) {
70 18 : if( getPntrToArgument(i)->getRank()!=2 ) error("input argument should be a matrix");
71 18 : if( wval[0]->getShape()[0]!=getPntrToArgument(i)->getShape()[0] || wval[0]->getShape()[1]!=getPntrToArgument(i)->getShape()[1] ) error("mismatched shapes of matrices in input");
72 : }
73 6 : log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() );
74 : // Rerequest the arguments
75 6 : std::vector<Value*> myargs( getArguments() ); myargs.push_back( wval[0] ); requestArguments( myargs );
76 30 : for(unsigned i=0; i<myargs.size(); ++i) myargs[i]->buildDataStore();
77 6 : std::vector<unsigned> shape(1); shape[0] = getPntrToArgument(0)->getShape()[0];
78 :
79 : // And now read the functions to compute
80 6 : for(int i=1;; i++) {
81 17 : std::string myfunc, mystr, lab, num; Tools::convert(i,num);
82 34 : if( !parseNumbered("FUNCTION",i,mystr ) ) break;
83 11 : std::vector<std::string> data=Tools::getWords(mystr);
84 22 : if( !Tools::parse(data,"LABEL",lab ) ) error("found no LABEL in FUNCTION" + num + " specification");
85 11 : addComponent( lab, shape ); componentIsNotPeriodic( lab );
86 22 : if( !Tools::parse(data,"FUNC",myfunc) ) error("found no FUNC in FUNCTION" + num + " specification");
87 11 : log.printf(" component labelled %s is computed using %s \n",lab.c_str(), myfunc.c_str() );
88 11 : functions.push_back( LeptonCall() ); std::vector<std::string> argnames(1); argnames[0]="ajik";
89 15 : if( myfunc.find("rij")!=std::string::npos ) argnames.push_back("rij");
90 11 : if( myfunc.find("rik")!=std::string::npos ) {
91 4 : if( argnames.size()<2 ) 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");
92 8 : argnames.push_back("rik");
93 : }
94 11 : if( myfunc.find("rjk")!=std::string::npos ) {
95 3 : if( argnames.size()<2 ) 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");
96 6 : argnames.push_back("rjk");
97 : }
98 11 : functions[i-1].set( myfunc, argnames, this, true );
99 22 : }
100 6 : checkRead();
101 6 : }
102 :
103 2 : std::string ThreeBodyGFunctions::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
104 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
105 6 : if( getConstPntrToComponent(i)->getName().find(cname)!=std::string::npos ) {
106 4 : std::string num; Tools::convert( i+1, num ); return "the function defined by the FUNCTION" + num + " keyword";
107 : }
108 : }
109 0 : plumed_error(); return "";
110 : }
111 :
112 36 : unsigned ThreeBodyGFunctions::getNumberOfDerivatives() {
113 36 : return 0;
114 : }
115 :
116 19 : void ThreeBodyGFunctions::calculate() {
117 19 : runAllTasks();
118 19 : }
119 :
120 472 : void ThreeBodyGFunctions::performTask( const unsigned& task_index, MultiValue& myvals ) const {
121 : const Value* wval = getPntrToArgument(3);
122 : const Value* xval = getPntrToArgument(0);
123 : const Value* yval = getPntrToArgument(1);
124 : const Value* zval = getPntrToArgument(2);
125 472 : Angle angle; Vector disti, distj; unsigned matsize = wval->getNumberOfValues();
126 472 : std::vector<double> values(4); std::vector<Vector> der_i(4), der_j(4);
127 472 : unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
128 13880 : for(unsigned i=0; i<nbonds; ++i) {
129 13408 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
130 13408 : double weighti = wval->get( ipos );
131 13408 : if( weighti<epsilon ) continue ;
132 6824 : disti[0] = xval->get( ipos );
133 6824 : disti[1] = yval->get( ipos );
134 6824 : disti[2] = zval->get( ipos );
135 6824 : values[1] = disti.modulo2(); der_i[1]=2*disti; der_i[2].zero();
136 189327 : for(unsigned j=0; j<i; ++j) {
137 182503 : unsigned jpos = ncols*task_index + wval->getRowIndex( task_index, j );
138 182503 : double weightj = wval->get( jpos );
139 182503 : if( weightj<epsilon ) continue ;
140 137009 : distj[0] = xval->get( jpos );
141 137009 : distj[1] = yval->get( jpos );
142 137009 : distj[2] = zval->get( jpos );
143 137009 : values[2] = distj.modulo2(); der_j[1].zero(); der_j[2]=2*distj;
144 137009 : der_i[3] = ( disti - distj ); values[3] = der_i[3].modulo2();
145 137009 : der_i[3] = 2*der_i[3]; der_j[3] = -der_i[3];
146 : // Compute angle between bonds
147 137009 : values[0] = angle.compute( disti, distj, der_i[0], der_j[0] );
148 : // Compute product of weights
149 137009 : double weightij = weighti*weightj;
150 : // Now compute all symmetry functions
151 414532 : for(unsigned n=0; n<functions.size(); ++n) {
152 277523 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
153 277523 : double nonweight = functions[n].evaluate( values ); myvals.addValue( ostrn, nonweight*weightij );
154 277523 : if( doNotCalculateDerivatives() ) continue;
155 :
156 567230 : for(unsigned m=0; m<functions[n].getNumberOfArguments(); ++m) {
157 296550 : double der = weightij*functions[n].evaluateDeriv( m, values );
158 296550 : myvals.addDerivative( ostrn, ipos, der*der_i[m][0] );
159 296550 : myvals.addDerivative( ostrn, matsize+ipos, der*der_i[m][1] );
160 296550 : myvals.addDerivative( ostrn, 2*matsize+ipos, der*der_i[m][2] );
161 296550 : myvals.addDerivative( ostrn, jpos, der*der_j[m][0] );
162 296550 : myvals.addDerivative( ostrn, matsize+jpos, der*der_j[m][1] );
163 296550 : myvals.addDerivative( ostrn, 2*matsize+jpos, der*der_j[m][2] );
164 : }
165 270680 : myvals.addDerivative( ostrn, 3*matsize+ipos, nonweight*weightj );
166 270680 : myvals.addDerivative( ostrn, 3*matsize+jpos, nonweight*weighti );
167 : }
168 : }
169 : }
170 472 : if( doNotCalculateDerivatives() ) return ;
171 :
172 : // And update the elements that have derivatives
173 : // Needs a separate loop here as there may be forces from j
174 8320 : for(unsigned i=0; i<nbonds; ++i) {
175 8192 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
176 8192 : double weighti = wval->get( ipos );
177 8192 : if( weighti<epsilon ) continue ;
178 :
179 16286 : for(unsigned n=0; n<functions.size(); ++n) {
180 11416 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
181 11416 : myvals.updateIndex( ostrn, ipos ); myvals.updateIndex( ostrn, matsize+ipos );
182 11416 : myvals.updateIndex( ostrn, 2*matsize+ipos ); myvals.updateIndex( ostrn, 3*matsize+ipos );
183 : }
184 : }
185 : }
186 :
187 : }
188 : }
|
