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