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