Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 "ActionAtomistic.h"
23 : #include "PlumedMain.h"
24 : #include "ActionSet.h"
25 : #include "GenericMolInfo.h"
26 : #include "PbcAction.h"
27 : #include <vector>
28 : #include <string>
29 : #include "ActionWithValue.h"
30 : #include "ActionShortcut.h"
31 : #include "Group.h"
32 : #include "ActionWithVirtualAtom.h"
33 : #include "tools/Exception.h"
34 : #include "tools/Pbc.h"
35 : #include "tools/PDB.h"
36 : #include "tools/Tree.h"
37 :
38 : namespace PLMD {
39 :
40 16699 : ActionAtomistic::~ActionAtomistic() {
41 : // empty destructor to delete unique_ptr
42 33398 : }
43 :
44 16699 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
45 : Action(ao),
46 16699 : unique_local_needs_update(true),
47 16699 : boxValue(NULL),
48 16699 : lockRequestAtoms(false),
49 16699 : donotretrieve(false),
50 16699 : donotforce(false),
51 16699 : massesWereSet(false),
52 16699 : chargesWereSet(false) {
53 16699 : ActionWithValue* bv = plumed.getActionSet().selectWithLabel<ActionWithValue*>("Box");
54 16699 : actionMoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
55 16699 : if( bv ) {
56 16547 : boxValue=bv->copyOutput(0);
57 : }
58 : // We now get all the information about atoms that are lying about
59 16699 : getAtomValuesFromPlumedObject( plumed, xpos, ypos, zpos, masv, chargev );
60 16699 : if( xpos.size()!=ypos.size() || xpos.size()!=zpos.size() || xpos.size()!=masv.size() || xpos.size()!=chargev.size() ) {
61 0 : error("mismatch between value arrays");
62 : }
63 16699 : }
64 :
65 16708 : void ActionAtomistic::getAtomValuesFromPlumedObject( const PlumedMain& plumed, std::vector<Value*>& xpos, std::vector<Value*>& ypos, std::vector<Value*>& zpos, std::vector<Value*>& masv, std::vector<Value*>& chargev ) {
66 16708 : std::vector<ActionShortcut*> shortcuts = plumed.getActionSet().select<ActionShortcut*>();
67 : bool foundpdb=false;
68 6230159 : for(const auto & ss : shortcuts ) {
69 6213451 : if( ss->getName()=="READMASSCHARGE" ) {
70 : foundpdb=true;
71 2 : ActionWithValue* mv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_mass");
72 2 : plumed_assert( mv );
73 2 : masv.push_back( mv->copyOutput(0) );
74 2 : ActionWithValue* qv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_charges");
75 2 : plumed_assert( qv );
76 2 : chargev.push_back( qv->copyOutput(0) );
77 : }
78 : }
79 16708 : std::vector<ActionWithValue*> vatoms = plumed.getActionSet().select<ActionWithValue*>();
80 6920956 : for(const auto & vv : vatoms ) {
81 6904248 : plumed_assert(vv); // needed for following calls, see #1046
82 6904248 : ActionToPutData* ap = vv->castToActionToPutData();
83 6904248 : if( ap ) {
84 231810 : if( ap->getRole()=="x" ) {
85 16554 : xpos.push_back( ap->copyOutput(0) );
86 : }
87 231810 : if( ap->getRole()=="y" ) {
88 16554 : ypos.push_back( ap->copyOutput(0) );
89 : }
90 231810 : if( ap->getRole()=="z" ) {
91 16554 : zpos.push_back( ap->copyOutput(0) );
92 : }
93 347687 : if( !foundpdb && ap->getRole()=="m" ) {
94 16552 : masv.push_back( ap->copyOutput(0) );
95 : }
96 347687 : if( !foundpdb && ap->getRole()=="q" ) {
97 16552 : chargev.push_back( ap->copyOutput(0) );
98 : }
99 : }
100 6904248 : ActionWithVirtualAtom* av = vv->castToActionWithVirtualAtom();
101 6904248 : if( av || vv->getName()=="ARGS2VATOM" ) {
102 6673573 : xpos.push_back( vv->copyOutput( vv->getLabel() + ".x") );
103 6673573 : ypos.push_back( vv->copyOutput( vv->getLabel() + ".y") );
104 6673573 : zpos.push_back( vv->copyOutput( vv->getLabel() + ".z") );
105 6673573 : masv.push_back( vv->copyOutput( vv->getLabel() + ".mass") );
106 6673573 : chargev.push_back( vv->copyOutput( vv->getLabel() + ".charge") );
107 : }
108 : }
109 16708 : }
110 :
111 29654 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
112 : (void) keys; // avoid warning
113 29654 : }
114 :
115 :
116 14173 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
117 14173 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
118 : // this makes the EMST invalid so that it will be regenerated at next request
119 : tree.reset();
120 14173 : int nat=a.size();
121 14173 : indexes=a;
122 14173 : actionPositions.resize(nat);
123 14173 : forces.resize(nat);
124 14173 : masses.resize(nat);
125 14173 : charges.resize(nat);
126 14173 : atom_value_ind.resize( a.size() );
127 14173 : unsigned n=getTotAtoms();
128 14173 : if(clearDep) {
129 13528 : clearDependencies();
130 : }
131 14173 : unique.clear();
132 14173 : std::vector<bool> requirements( xpos.size(), false );
133 14173 : if( boxValue ) {
134 14173 : addDependency( boxValue->getPntrToAction() );
135 : }
136 1703091 : for(unsigned i=0; i<indexes.size(); i++) {
137 1688919 : if(indexes[i].index()>=n) {
138 : std::string num;
139 1 : Tools::convert( indexes[i].serial(),num );
140 3 : error("atom " + num + " out of range");
141 : }
142 1688918 : atom_value_ind[i] = getValueIndices( indexes[i] );
143 1688918 : requirements[atom_value_ind[i].first] = true;
144 1688918 : if( atom_value_ind[i].first==0 ) {
145 1673705 : unique.push_back(indexes[i]);
146 15213 : } else if( atom_value_ind[i].second>0 ) {
147 0 : error("action atomistic is not set up to deal with multiple vectors in position input");
148 : }
149 : }
150 :
151 14172 : atom_value_ind_grouped.clear();
152 :
153 14172 : if(atom_value_ind.size()>0) {
154 14130 : auto prev_nn=atom_value_ind[0].first;
155 : atom_value_ind_grouped.push_back(
156 14130 : std::pair<std::size_t,std::vector<std::size_t>>(
157 : prev_nn, {}));
158 14130 : atom_value_ind_grouped.back().second.push_back(atom_value_ind[0].second);
159 1688911 : for(unsigned i=1; i<atom_value_ind.size(); i++) {
160 1674781 : auto nn = atom_value_ind[i].first;
161 1674781 : auto kk = atom_value_ind[i].second;
162 1674781 : if(nn!=prev_nn)
163 16887 : atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
164 1674781 : atom_value_ind_grouped.back().second.push_back(kk);
165 : prev_nn=nn;
166 : }
167 : }
168 :
169 : // Add the dependencies to the actions that we require
170 14172 : Tools::removeDuplicates(unique);
171 14172 : value_depends.resize(0);
172 6633192 : for(unsigned i=0; i<requirements.size(); ++i ) {
173 6619020 : if( !requirements[i] ) {
174 6590939 : continue;
175 : }
176 28082 : value_depends.push_back( i );
177 28081 : addDependency( xpos[i]->getPntrToAction() );
178 28081 : addDependency( ypos[i]->getPntrToAction() );
179 28081 : addDependency( zpos[i]->getPntrToAction() );
180 28081 : addDependency( masv[i]->getPntrToAction() );
181 28081 : addDependency( chargev[i]->getPntrToAction() );
182 : }
183 14172 : unique_local_needs_update=true;
184 14172 : }
185 :
186 0 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
187 0 : actionPbc.apply(dlist, max_index);
188 0 : }
189 :
190 161 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
191 161 : calculateAtomicNumericalDerivatives( a, 0 );
192 161 : }
193 :
194 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
195 0 : actionPbc.setBox( newbox );
196 0 : }
197 :
198 246 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
199 246 : if(!a) {
200 246 : a=castToActionWithValue();
201 246 : plumed_massert(a,"only Actions with a value can be differentiated");
202 : }
203 :
204 : const size_t nval=a->getNumberOfComponents();
205 246 : const size_t natoms=getNumberOfAtoms();
206 246 : std::vector<Vector> value(nval*natoms);
207 246 : std::vector<Tensor> valuebox(nval);
208 246 : std::vector<Vector> savedPositions(natoms);
209 : const double delta=std::sqrt(epsilon);
210 :
211 21325 : for(unsigned i=0; i<natoms; i++)
212 84316 : for(int k=0; k<3; k++) {
213 63237 : savedPositions[i][k]=actionPositions[i][k];
214 63237 : actionPositions[i][k]=actionPositions[i][k]+delta;
215 63237 : a->calculate();
216 63237 : actionPositions[i][k]=savedPositions[i][k];
217 197478 : for(unsigned j=0; j<nval; j++) {
218 134241 : value[j*natoms+i][k]=a->getOutputQuantity(j);
219 : }
220 : }
221 246 : Tensor box(actionPbc.getBox());
222 984 : for(int i=0; i<3; i++)
223 2952 : for(int k=0; k<3; k++) {
224 2214 : double arg0=box(i,k);
225 191925 : for(unsigned j=0; j<natoms; j++) {
226 189711 : actionPositions[j]=actionPbc.realToScaled(actionPositions[j]);
227 : }
228 2214 : box(i,k)=box(i,k)+delta;
229 2214 : actionPbc.setBox(box);
230 191925 : for(unsigned j=0; j<natoms; j++) {
231 189711 : actionPositions[j]=actionPbc.scaledToReal(actionPositions[j]);
232 : }
233 2214 : a->calculate();
234 2214 : box(i,k)=arg0;
235 2214 : actionPbc.setBox(box);
236 191925 : for(unsigned j=0; j<natoms; j++) {
237 189711 : actionPositions[j]=savedPositions[j];
238 : }
239 15111 : for(unsigned j=0; j<nval; j++) {
240 12897 : valuebox[j](i,k)=a->getOutputQuantity(j);
241 : }
242 : }
243 :
244 246 : a->calculate();
245 246 : a->clearDerivatives();
246 1679 : for(unsigned j=0; j<nval; j++) {
247 1433 : Value* v=a->copyOutput(j);
248 1433 : double ref=v->get();
249 1433 : if(v->hasDerivatives()) {
250 23482 : for(unsigned i=0; i<natoms; i++)
251 91948 : for(int k=0; k<3; k++) {
252 68961 : double d=(value[j*natoms+i][k]-ref)/delta;
253 68961 : v->addDerivative(startnum+3*i+k,d);
254 : }
255 : Tensor virial;
256 1980 : for(int i=0; i<3; i++)
257 5940 : for(int k=0; k<3; k++) {
258 4455 : virial(i,k)= (valuebox[j](i,k)-ref)/delta;
259 : }
260 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
261 495 : virial=-matmul(box.transpose(),virial);
262 1980 : for(int i=0; i<3; i++)
263 5940 : for(int k=0; k<3; k++) {
264 4455 : v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
265 : }
266 : }
267 : }
268 246 : }
269 :
270 105884 : bool ActionAtomistic::actionHasForces() {
271 105884 : ActionWithValue* av = castToActionWithValue();
272 105884 : if( av ) {
273 105884 : return !av->doNotCalculateDerivatives();
274 : }
275 0 : if( indexes.size()>0 ) {
276 0 : plumed_merror("you have to overwrite the function actionHasForce to tell plumed if you method applies forces");
277 : }
278 : return true;
279 : }
280 :
281 13706 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
282 13706 : parseAtomList(key,-1,t);
283 13706 : }
284 :
285 212188 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
286 214534 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden") || keywords.style(key,"deprecated"), "keyword " + key + " should be registered as atoms");
287 : std::vector<std::string> strings;
288 212188 : if( num<0 ) {
289 17267 : parseVector(key,strings);
290 17267 : if(strings.empty()) {
291 : return;
292 : }
293 : } else {
294 194921 : if ( !parseNumberedVector(key,num,strings) ) {
295 : return;
296 : }
297 : }
298 207963 : t.resize(0);
299 207963 : interpretAtomList( strings, xpos, this, t );
300 212188 : }
301 :
302 33 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
303 33 : interpretAtomList( strings, xpos, this, t );
304 33 : }
305 :
306 208582 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t) {
307 208582 : Tools::interpretRanges(strings);
308 1717071 : for(unsigned i=0; i<strings.size(); ++i) {
309 : AtomNumber atom;
310 1508489 : bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
311 1508489 : if(ok) {
312 1482890 : t.push_back(atom);
313 : }
314 : // here we check if this is a special symbol for MOLINFO
315 1508489 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
316 837 : std::string symbol=strings[i].substr(1);
317 837 : if(symbol=="allatoms") {
318 : unsigned n=0;
319 24 : for(unsigned ii=0; ii<xpos.size(); ++ii) {
320 14 : n += xpos[ii]->getNumberOfValues();
321 : }
322 10 : t.reserve(n);
323 365 : for(unsigned ii=0; ii<n; ii++) {
324 355 : t.push_back(AtomNumber::index(ii));
325 : }
326 : ok=true;
327 827 : } else if(symbol=="mdatoms") {
328 15 : const auto n=xpos[0]->getNumberOfValues();
329 15 : t.reserve(t.size()+n);
330 1181 : for(unsigned ii=0; ii<n; ii++) {
331 1166 : t.push_back(AtomNumber::index(ii));
332 : }
333 : ok=true;
334 812 : } else if(Tools::startWith(symbol,"ndx:")) {
335 50 : auto words=Tools::getWords(symbol.substr(4));
336 : std::string ndxfile,ndxgroup;
337 25 : if(words.size()==1) {
338 : ndxfile=words[0];
339 23 : } else if(words.size()==2) {
340 : ndxfile=words[0];
341 : ndxgroup=words[1];
342 : } else {
343 0 : plumed_error()<<"Cannot intepret selection "<<symbol;
344 : }
345 :
346 25 : if(ndxgroup.size()>0) {
347 46 : action->log<<" importing group '"+ndxgroup+"'";
348 : } else {
349 2 : action->log<<" importing first group";
350 : }
351 25 : action->log<<" from index file "<<ndxfile<<"\n";
352 :
353 25 : IFile ifile;
354 25 : ifile.open(ndxfile);
355 : std::string line;
356 : std::string groupname;
357 : bool firstgroup=true;
358 : bool groupfound=false;
359 5799 : while(ifile.getline(line)) {
360 5774 : std::vector<std::string> groupWords=Tools::getWords(line);
361 11680 : if(groupWords.size()>=3 && groupWords[0]=="[" && groupWords[2]=="]") {
362 218 : if(groupname.length()>0) {
363 : firstgroup=false;
364 : }
365 : groupname=groupWords[1];
366 218 : if(groupname==ndxgroup || ndxgroup.length()==0) {
367 : groupfound=true;
368 : }
369 5556 : } else if(groupname==ndxgroup || (firstgroup && ndxgroup.length()==0)) {
370 9401 : for(unsigned ii=0; ii<groupWords.size(); ii++) {
371 : AtomNumber at;
372 8792 : Tools::convert(groupWords[ii],at);
373 8792 : t.push_back(at);
374 : }
375 : }
376 5774 : }
377 25 : if(!groupfound) {
378 0 : plumed_error()<<"group has not been found in index file";
379 : }
380 : ok=true;
381 50 : } else {
382 787 : auto* moldat=action->plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
383 787 : if( moldat ) {
384 : std::vector<AtomNumber> atom_list;
385 787 : moldat->interpretSymbol( symbol, atom_list );
386 787 : if( atom_list.size()>0 ) {
387 : ok=true;
388 787 : t.insert(t.end(),atom_list.begin(),atom_list.end());
389 : } else {
390 0 : action->error(strings[i] + " is not a label plumed knows");
391 : }
392 : } else {
393 0 : action->error("atoms specified using @ symbol but no MOLINFO was available");
394 : }
395 : }
396 : }
397 : // here we check if the atom name is the name of a group
398 1507652 : if(!ok) {
399 24762 : Group* mygrp=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]);
400 24762 : if(mygrp) {
401 466 : std::vector<std::string> grp_str( mygrp->getGroupAtoms() );
402 466 : interpretAtomList( grp_str, xpos, action, t );
403 : ok=true;
404 466 : } else {
405 24296 : Group* mygrp2=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]+"_grp");
406 24296 : if(mygrp2) {
407 112 : std::vector<std::string> grp_str( mygrp2->getGroupAtoms() );
408 112 : interpretAtomList( grp_str, xpos, action, t );
409 : ok=true;
410 112 : }
411 : }
412 : }
413 : // here we check if the atom name is the name of an added virtual atom
414 1507911 : if(!ok) {
415 : unsigned ind = 0;
416 13304972 : for(unsigned j=0; j<xpos.size(); ++j) {
417 13304972 : if( xpos[j]->getPntrToAction()->getLabel()==strings[i] ) {
418 24184 : t.push_back( AtomNumber::index(ind) );
419 : ok=true;
420 : break;
421 : }
422 13280788 : ind = ind + xpos[j]->getNumberOfValues();
423 : }
424 : }
425 1484305 : if(!ok) {
426 0 : action->error("it was not possible to interpret atom name " + strings[i]);
427 : }
428 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
429 : }
430 208582 : }
431 :
432 2174594 : std::pair<std::size_t, std::size_t> ActionAtomistic::getValueIndices( const AtomNumber& i ) const {
433 2174594 : std::size_t valno=0, k = i.index();
434 15620054 : for(unsigned j=0; j<xpos.size(); ++j) {
435 15620054 : if( k<xpos[j]->getNumberOfValues() ) {
436 : valno=j;
437 : break;
438 : }
439 13445460 : k = k - xpos[j]->getNumberOfValues();
440 : }
441 2174594 : return std::pair<std::size_t, std::size_t>( valno, k );
442 : }
443 :
444 448128 : void ActionAtomistic::retrieveAtoms( const bool& force ) {
445 448128 : if( boxValue ) {
446 431486 : auto* ptr=boxValue->getPntrToAction();
447 431486 : plumed_assert(ptr); // needed for following calls, see #1046
448 431486 : PbcAction* pbca = ptr->castToPbcAction();
449 431486 : plumed_assert( pbca );
450 431486 : actionPbc=pbca->pbc;
451 : }
452 448128 : if( donotretrieve || indexes.size()==0 ) {
453 : return;
454 : }
455 236100 : auto * mtr=masv[0]->getPntrToAction();
456 236100 : plumed_assert(mtr); // needed for following calls, see #1046
457 236100 : ActionToPutData* mv = mtr->castToActionToPutData();
458 236100 : if(mv) {
459 236098 : massesWereSet=mv->hasBeenSet();
460 2 : } else if( (masv[0]->getPntrToAction())->getName()=="CONSTANT" ) {
461 2 : massesWereSet=true; // Read masses from PDB file
462 : }
463 236100 : auto * ptr=chargev[0]->getPntrToAction();
464 236100 : plumed_assert(ptr); // needed for following calls, see #1046
465 236100 : ActionToPutData* cv = ptr->castToActionToPutData();
466 236100 : if(cv) {
467 236098 : chargesWereSet=cv->hasBeenSet();
468 2 : } else if( (chargev[0]->getPntrToAction())->getName()=="CONSTANT" ) {
469 2 : chargesWereSet=true; // Read masses from PDB file
470 : }
471 : unsigned j = 0;
472 :
473 : // for(const auto & a : atom_value_ind) {
474 : // std::size_t nn = a.first, kk = a.second;
475 : // actionPositions[j][0] = xpos[nn]->data[kk];
476 : // actionPositions[j][1] = ypos[nn]->data[kk];
477 : // actionPositions[j][2] = zpos[nn]->data[kk];
478 : // charges[j] = chargev[nn]->data[kk];
479 : // masses[j] = masv[nn]->data[kk];
480 : // j++;
481 : // }
482 :
483 874903 : for(const auto & a : atom_value_ind_grouped) {
484 638803 : const auto nn=a.first;
485 638803 : auto & xp=xpos[nn]->data;
486 638803 : auto & yp=ypos[nn]->data;
487 638803 : auto & zp=zpos[nn]->data;
488 638803 : auto & ch=chargev[nn]->data;
489 638803 : auto & ma=masv[nn]->data;
490 9497887 : for(const auto & kk : a.second) {
491 8859084 : actionPositions[j][0] = xp[kk];
492 8859084 : actionPositions[j][1] = yp[kk];
493 8859084 : actionPositions[j][2] = zp[kk];
494 8859084 : charges[j] = ch[kk];
495 8859084 : masses[j] = ma[kk];
496 8859084 : j++;
497 : }
498 : }
499 :
500 : }
501 :
502 304082 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned& ind) {
503 304082 : if( donotforce || (indexes.size()==0 && getName()!="FIXEDATOM") ) {
504 154561 : return;
505 : }
506 344296 : for(unsigned i=0; i<value_depends.size(); ++i) {
507 194775 : xpos[value_depends[i]]->hasForce = true;
508 194775 : ypos[value_depends[i]]->hasForce = true;
509 194775 : zpos[value_depends[i]]->hasForce = true;
510 : }
511 :
512 : // for(const auto & a : atom_value_ind) {
513 : // plumed_dbg_massert( ind<forcesToApply.size(), "problem setting forces in " + getLabel() );
514 : // std::size_t nn = a.first, kk = a.second;
515 : // xpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
516 : // ypos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
517 : // zpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
518 : // }
519 :
520 654684 : for(const auto & a : atom_value_ind_grouped) {
521 505163 : const auto nn=a.first;
522 : plumed_dbg_assert(ind<forcesToApply.size()) << "problem setting forces in " << getLabel();
523 505163 : auto & xp=xpos[nn]->inputForce;
524 505163 : auto & yp=ypos[nn]->inputForce;
525 505163 : auto & zp=zpos[nn]->inputForce;
526 5797849 : for(const auto & kk : a.second) {
527 5292686 : xp[kk] += forcesToApply[ind];
528 5292686 : ind++;
529 5292686 : yp[kk] += forcesToApply[ind];
530 5292686 : ind++;
531 5292686 : zp[kk] += forcesToApply[ind];
532 5292686 : ind++;
533 : }
534 : }
535 :
536 149521 : setForcesOnCell( forcesToApply, ind );
537 : }
538 :
539 149609 : void ActionAtomistic::setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind) {
540 149609 : setForcesOnCell(forcesToApply.data(),forcesToApply.size(),ind);
541 149609 : }
542 :
543 172554 : void ActionAtomistic::setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind) {
544 1725540 : for(unsigned i=0; i<9; ++i) {
545 : plumed_dbg_massert( ind<size, "problem setting forces in " + getLabel() );
546 1552986 : boxValue->addForce( i, forcesToApply[ind] );
547 1552986 : ind++;
548 : }
549 172554 : }
550 :
551 6 : Tensor ActionAtomistic::getVirial() const {
552 : Tensor vir;
553 24 : for(unsigned i=0; i<3; ++i)
554 72 : for(unsigned j=0; j<3; ++j) {
555 54 : vir[i][j] = boxValue->getForce(3*i+j);
556 : }
557 6 : return vir;
558 : }
559 :
560 0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
561 :
562 0 : for(unsigned j=0; j<indexes.size(); j++) {
563 0 : if( indexes[j].index()>pdb.size() ) {
564 0 : error("there are not enough atoms in the input pdb file");
565 : }
566 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) {
567 0 : error("there are atoms missing in the pdb file");
568 : }
569 0 : actionPositions[j]=pdb.getPositions()[indexes[j].index()];
570 : }
571 0 : for(unsigned j=0; j<indexes.size(); j++) {
572 0 : charges[j]=pdb.getBeta()[indexes[j].index()];
573 : }
574 0 : for(unsigned j=0; j<indexes.size(); j++) {
575 0 : masses[j]=pdb.getOccupancy()[indexes[j].index()];
576 : }
577 0 : }
578 :
579 14444 : unsigned ActionAtomistic::getTotAtoms()const {
580 : unsigned natoms = 0;
581 6633808 : for(unsigned i=0; i<xpos.size(); ++i ) {
582 6619364 : natoms += xpos[i]->getNumberOfValues();
583 : }
584 14444 : return natoms;
585 : }
586 :
587 154004 : void ActionAtomistic::makeWhole() {
588 154004 : if(actionMoldat && actionMoldat->isWhole()) {
589 : // make sure the tree has been constructed
590 2 : if(!tree) {
591 2 : tree=std::make_unique<Tree>(actionMoldat);
592 : }
593 2 : const auto & tree_indexes=tree->getTreeIndexes();
594 2 : const auto & root_indexes=tree->getRootIndexes();
595 2 : for(unsigned j=0; j<root_indexes.size(); j++) {
596 0 : const Vector & first (actionPositions[root_indexes[j]]);
597 0 : Vector & second (actionPositions[tree_indexes[j]]);
598 0 : second=first+pbcDistance(first,second);
599 : }
600 : } else {
601 1387990 : for(unsigned j=0; j<actionPositions.size()-1; ++j) {
602 : const Vector & first (actionPositions[j]);
603 1233988 : Vector & second (actionPositions[j+1]);
604 1233988 : second=first+pbcDistance(first,second);
605 : }
606 : }
607 154004 : }
608 :
609 8208 : void ActionAtomistic::getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const {
610 8208 : std::size_t nn = atom_value_ind[ind].first;
611 8208 : if( nn==0 ) {
612 8169 : gradients[indexes[ind]] += deriv;
613 8169 : return;
614 : }
615 39 : xpos[nn]->passGradients( deriv[0], gradients );
616 39 : ypos[nn]->passGradients( deriv[1], gradients );
617 39 : zpos[nn]->passGradients( deriv[2], gradients );
618 : }
619 :
620 240174 : void ActionAtomistic::updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l ) {
621 240174 : if( useunique ) {
622 231020 : unique_local=unique;
623 231020 : return;
624 : }
625 : // Update unique local if it needs an update
626 9154 : unique_local_needs_update=false;
627 9154 : unique_local.clear();
628 75612 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
629 66458 : if(g2l[pp->index()]>=0) {
630 26916 : unique_local.push_back(*pp); // already sorted
631 : }
632 : }
633 : }
634 :
635 : }
|
