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 18779 : ActionAtomistic::~ActionAtomistic() {
41 : // empty destructor to delete unique_ptr
42 56337 : }
43 :
44 18779 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
45 : Action(ao),
46 18779 : unique_local_needs_update(true),
47 18779 : boxValue(NULL),
48 18779 : lockRequestAtoms(false),
49 18779 : donotretrieve(false),
50 18779 : donotforce(false),
51 18779 : massesWereSet(false),
52 18779 : chargesWereSet(false) {
53 18779 : ActionWithValue* bv = plumed.getActionSet().selectWithLabel<ActionWithValue*>("Box");
54 18779 : actionMoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
55 18779 : if( bv ) {
56 18615 : boxValue=bv->copyOutput(0);
57 : }
58 : // We now get all the information about atoms that are lying about
59 18779 : getAtomValuesFromPlumedObject( plumed, xpos, ypos, zpos, masv, chargev );
60 18779 : 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 18779 : }
64 :
65 28325 : 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 : // the lists are kept by the ActionSet, so that they are not rebuilt by scanning
67 : // all the actions every time a new atomistic action is created
68 28325 : plumed.getActionSet().getAtomValues( xpos, ypos, zpos, masv, chargev );
69 28325 : }
70 :
71 34911 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
72 : (void) keys; // avoid warning
73 34911 : }
74 :
75 :
76 15549 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
77 15549 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
78 : // this makes the EMST invalid so that it will be regenerated at next request
79 : tree.reset();
80 15549 : int nat=a.size();
81 15549 : indexes=a;
82 15549 : actionPositions.resize(nat);
83 15549 : forces.resize(nat);
84 15549 : masses.resize(nat);
85 15549 : charges.resize(nat);
86 15549 : atom_value_ind.resize( a.size() );
87 15549 : unsigned n=getTotAtoms();
88 15549 : if(clearDep) {
89 14715 : clearDependencies();
90 : }
91 15549 : unique.clear();
92 15549 : std::vector<bool> requirements( xpos.size(), false );
93 15549 : if( boxValue ) {
94 15549 : addDependency( boxValue->getPntrToAction() );
95 : }
96 1508592 : for(unsigned i=0; i<indexes.size(); i++) {
97 1493044 : if(indexes[i].index()>=n) {
98 : std::string num;
99 1 : Tools::convert( indexes[i].serial(),num );
100 3 : error("atom " + num + " out of range");
101 : }
102 1493043 : atom_value_ind[i] = getValueIndices( indexes[i] );
103 1493043 : requirements[atom_value_ind[i].first] = true;
104 1493043 : if( atom_value_ind[i].first==0 ) {
105 1477761 : unique.push_back(indexes[i]);
106 15282 : } else if( atom_value_ind[i].second>0 ) {
107 0 : error("action atomistic is not set up to deal with multiple vectors in position input");
108 : }
109 : }
110 :
111 15548 : atom_value_ind_grouped.clear();
112 :
113 15548 : if(atom_value_ind.size()>0) {
114 15494 : auto prev_nn=atom_value_ind[0].first;
115 : atom_value_ind_grouped.push_back(
116 15494 : std::pair<std::size_t,std::vector<std::size_t>>(
117 : prev_nn, {}));
118 15494 : atom_value_ind_grouped.back().second.push_back(atom_value_ind[0].second);
119 1493036 : for(unsigned i=1; i<atom_value_ind.size(); i++) {
120 1477542 : auto nn = atom_value_ind[i].first;
121 1477542 : auto kk = atom_value_ind[i].second;
122 1477542 : if(nn!=prev_nn)
123 16966 : atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
124 1477542 : atom_value_ind_grouped.back().second.push_back(kk);
125 : prev_nn=nn;
126 : }
127 : }
128 :
129 : // Add the dependencies to the actions that we require
130 15548 : Tools::removeDuplicates(unique);
131 15548 : value_depends.resize(0);
132 6636057 : for(unsigned i=0; i<requirements.size(); ++i ) {
133 6620509 : if( !requirements[i] ) {
134 6591039 : continue;
135 : }
136 29471 : value_depends.push_back( i );
137 29470 : addDependency( xpos[i]->getPntrToAction() );
138 29470 : addDependency( ypos[i]->getPntrToAction() );
139 29470 : addDependency( zpos[i]->getPntrToAction() );
140 29470 : addDependency( masv[i]->getPntrToAction() );
141 29470 : addDependency( chargev[i]->getPntrToAction() );
142 : }
143 15548 : unique_local_needs_update=true;
144 15548 : }
145 :
146 0 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
147 0 : actionPbc.apply(dlist, max_index);
148 0 : }
149 :
150 186 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
151 186 : calculateAtomicNumericalDerivatives( a, 0 );
152 186 : }
153 :
154 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
155 0 : actionPbc.setBox( newbox );
156 0 : }
157 :
158 261 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
159 261 : if(!a) {
160 261 : a=castToActionWithValue();
161 261 : plumed_massert(a,"only Actions with a value can be differentiated");
162 : }
163 :
164 : const size_t nval=a->getNumberOfComponents();
165 261 : const size_t natoms=getNumberOfAtoms();
166 261 : std::vector<Vector> value(nval*natoms);
167 261 : std::vector<Tensor> valuebox(nval);
168 261 : std::vector<Vector> savedPositions(natoms);
169 : const double delta=std::sqrt(epsilon);
170 :
171 15445 : for(unsigned i=0; i<natoms; i++)
172 60736 : for(int k=0; k<3; k++) {
173 45552 : savedPositions[i][k]=actionPositions[i][k];
174 45552 : actionPositions[i][k]=actionPositions[i][k]+delta;
175 45552 : a->calculate();
176 45552 : actionPositions[i][k]=savedPositions[i][k];
177 118764 : for(unsigned j=0; j<nval; j++) {
178 73212 : value[j*natoms+i][k]=a->getOutputQuantity(j);
179 : }
180 : }
181 261 : Tensor box(actionPbc.getBox());
182 1044 : for(int i=0; i<3; i++)
183 3132 : for(int k=0; k<3; k++) {
184 2349 : double arg0=box(i,k);
185 139005 : for(unsigned j=0; j<natoms; j++) {
186 136656 : actionPositions[j]=actionPbc.realToScaled(actionPositions[j]);
187 : }
188 2349 : box(i,k)=box(i,k)+delta;
189 2349 : actionPbc.setBox(box);
190 139005 : for(unsigned j=0; j<natoms; j++) {
191 136656 : actionPositions[j]=actionPbc.scaledToReal(actionPositions[j]);
192 : }
193 2349 : a->calculate();
194 2349 : box(i,k)=arg0;
195 2349 : actionPbc.setBox(box);
196 139005 : for(unsigned j=0; j<natoms; j++) {
197 136656 : actionPositions[j]=savedPositions[j];
198 : }
199 15165 : for(unsigned j=0; j<nval; j++) {
200 12816 : valuebox[j](i,k)=a->getOutputQuantity(j);
201 : }
202 : }
203 :
204 261 : a->calculate();
205 261 : a->clearDerivatives();
206 1685 : for(unsigned j=0; j<nval; j++) {
207 1424 : Value* v=a->copyOutput(j);
208 1424 : double ref=v->get();
209 1424 : if(v->hasDerivatives()) {
210 15874 : for(unsigned i=0; i<natoms; i++)
211 62224 : for(int k=0; k<3; k++) {
212 46668 : double d=(value[j*natoms+i][k]-ref)/delta;
213 46668 : v->addDerivative(startnum+3*i+k,d);
214 : }
215 318 : Tensor virial;
216 1272 : for(int i=0; i<3; i++)
217 3816 : for(int k=0; k<3; k++) {
218 2862 : virial(i,k)= (valuebox[j](i,k)-ref)/delta;
219 : }
220 : // BE CAREFUL WITH NON ORTHORHOMBIC CELL
221 318 : virial=-matmul(box.transpose(),virial);
222 1272 : for(int i=0; i<3; i++)
223 3816 : for(int k=0; k<3; k++) {
224 2862 : v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
225 : }
226 : }
227 : }
228 261 : }
229 :
230 383545 : bool ActionAtomistic::actionHasForces() {
231 383545 : ActionWithValue* av = castToActionWithValue();
232 383545 : if( av ) {
233 383545 : return !av->doNotCalculateDerivatives();
234 : }
235 0 : if( indexes.size()>0 ) {
236 0 : plumed_merror("you have to overwrite the function actionHasForce to tell plumed if you method applies forces");
237 : }
238 : return true;
239 : }
240 :
241 14305 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
242 14305 : parseAtomList(key,-1,t);
243 14305 : }
244 :
245 255239 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
246 258008 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden") || keywords.style(key,"deprecated"), "keyword " + key + " should be registered as atoms");
247 255239 : std::vector<std::string> strings;
248 255239 : if( num<0 ) {
249 19308 : parseVector(key,strings);
250 19308 : if(strings.empty()) {
251 : return;
252 : }
253 : } else {
254 235931 : if ( !parseNumberedVector(key,num,strings) ) {
255 : return;
256 : }
257 : }
258 249985 : t.resize(0);
259 249985 : interpretAtomList( strings, xpos, this, t );
260 255239 : }
261 :
262 40 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
263 40 : interpretAtomList( strings, xpos, this, t );
264 40 : }
265 :
266 250661 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t) {
267 250661 : Tools::interpretRanges(strings);
268 1879664 : for(unsigned i=0; i<strings.size(); ++i) {
269 : AtomNumber atom;
270 1629003 : bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
271 1629003 : if(ok) {
272 1603240 : t.push_back(atom);
273 : }
274 : // here we check if this is a special symbol for MOLINFO
275 25763 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
276 859 : std::string symbol=strings[i].substr(1);
277 859 : if(symbol=="allatoms") {
278 : unsigned n=0;
279 24 : for(unsigned ii=0; ii<xpos.size(); ++ii) {
280 14 : n += xpos[ii]->getNumberOfValues();
281 : }
282 10 : t.reserve(n);
283 365 : for(unsigned ii=0; ii<n; ii++) {
284 355 : t.push_back(AtomNumber::index(ii));
285 : }
286 : ok=true;
287 849 : } else if(symbol=="mdatoms") {
288 16 : const auto n=xpos[0]->getNumberOfValues();
289 16 : t.reserve(t.size()+n);
290 1246 : for(unsigned ii=0; ii<n; ii++) {
291 1230 : t.push_back(AtomNumber::index(ii));
292 : }
293 : ok=true;
294 833 : } else if(Tools::startWith(symbol,"ndx:")) {
295 64 : auto words=Tools::getWords(symbol.substr(4));
296 : std::string ndxfile,ndxgroup;
297 32 : if(words.size()==1) {
298 : ndxfile=words[0];
299 30 : } else if(words.size()==2) {
300 : ndxfile=words[0];
301 : ndxgroup=words[1];
302 : } else {
303 0 : plumed_error()<<"Cannot intepret selection "<<symbol;
304 : }
305 :
306 32 : if(ndxgroup.size()>0) {
307 90 : action->log<<" importing group '"+ndxgroup+"'";
308 : } else {
309 2 : action->log<<" importing first group";
310 : }
311 32 : action->log<<" from index file "<<ndxfile<<"\n";
312 :
313 32 : IFile ifile;
314 32 : ifile.open(ndxfile);
315 : std::string line;
316 : std::string groupname;
317 : bool firstgroup=true;
318 : bool groupfound=false;
319 5990 : while(ifile.getline(line)) {
320 5958 : std::vector<std::string> groupWords=Tools::getWords(line);
321 5958 : if(groupWords.size()>=3 && groupWords[0]=="[" && groupWords[2]=="]") {
322 288 : if(groupname.length()>0) {
323 : firstgroup=false;
324 : }
325 : groupname=groupWords[1];
326 288 : if(groupname==ndxgroup || ndxgroup.length()==0) {
327 : groupfound=true;
328 : }
329 5670 : } else if(groupname==ndxgroup || (firstgroup && ndxgroup.length()==0)) {
330 9532 : for(unsigned ii=0; ii<groupWords.size(); ii++) {
331 : AtomNumber at;
332 8912 : Tools::convert(groupWords[ii],at);
333 8912 : t.push_back(at);
334 : }
335 : }
336 5958 : }
337 32 : if(!groupfound) {
338 0 : plumed_error()<<"group has not been found in index file";
339 : }
340 : ok=true;
341 64 : } else {
342 801 : auto* moldat=action->plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
343 801 : if( moldat ) {
344 801 : std::vector<AtomNumber> atom_list;
345 801 : moldat->interpretSymbol( symbol, atom_list );
346 801 : if( atom_list.size()>0 ) {
347 : ok=true;
348 801 : t.insert(t.end(),atom_list.begin(),atom_list.end());
349 : } else {
350 0 : action->error(strings[i] + " is not a label plumed knows");
351 : }
352 : } else {
353 0 : action->error("atoms specified using @ symbol but no MOLINFO was available");
354 : }
355 : }
356 : }
357 : // here we check if the atom name is the name of a group
358 1628144 : if(!ok) {
359 24904 : Group* mygrp=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]);
360 24904 : if(mygrp) {
361 509 : std::vector<std::string> grp_str( mygrp->getGroupAtoms() );
362 509 : interpretAtomList( grp_str, xpos, action, t );
363 : ok=true;
364 509 : } else {
365 24395 : Group* mygrp2=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]+"_grp");
366 24395 : if(mygrp2) {
367 118 : std::vector<std::string> grp_str( mygrp2->getGroupAtoms() );
368 118 : interpretAtomList( grp_str, xpos, action, t );
369 : ok=true;
370 118 : }
371 : }
372 : }
373 : // here we check if the atom name is the name of an added virtual atom
374 : if(!ok) {
375 : unsigned ind = 0;
376 13305227 : for(unsigned j=0; j<xpos.size(); ++j) {
377 13305227 : if( xpos[j]->getPntrToAction()->getLabel()==strings[i] ) {
378 24277 : t.push_back( AtomNumber::index(ind) );
379 : ok=true;
380 : break;
381 : }
382 13280950 : ind = ind + xpos[j]->getNumberOfValues();
383 : }
384 : }
385 : if(!ok) {
386 0 : action->error("it was not possible to interpret atom name " + strings[i]);
387 : }
388 : }
389 250661 : }
390 :
391 2034198 : std::pair<std::size_t, std::size_t> ActionAtomistic::getValueIndices( const AtomNumber& i ) const {
392 2034198 : std::size_t valno=0, k = i.index();
393 15479728 : for(unsigned j=0; j<xpos.size(); ++j) {
394 15479728 : if( k<xpos[j]->getNumberOfValues() ) {
395 : valno=j;
396 : break;
397 : }
398 13445530 : k = k - xpos[j]->getNumberOfValues();
399 : }
400 2034198 : return std::pair<std::size_t, std::size_t>( valno, k );
401 : }
402 :
403 606519 : void ActionAtomistic::retrieveAtoms( const bool& force ) {
404 606519 : if( boxValue ) {
405 576218 : auto* ptr=boxValue->getPntrToAction();
406 576218 : plumed_assert(ptr); // needed for following calls, see #1046
407 576218 : PbcAction* pbca = ptr->castToPbcAction();
408 576218 : plumed_assert( pbca );
409 576218 : actionPbc=pbca->pbc;
410 : }
411 606519 : if( donotretrieve || indexes.size()==0 ) {
412 : return;
413 : }
414 407739 : auto * mtr=masv[0]->getPntrToAction();
415 407739 : plumed_assert(mtr); // needed for following calls, see #1046
416 407739 : ActionToPutData* mv = mtr->castToActionToPutData();
417 407739 : if(mv) {
418 407737 : massesWereSet=mv->hasBeenSet();
419 2 : } else if( (masv[0]->getPntrToAction())->getName()=="CONSTANT" ) {
420 2 : massesWereSet=true; // Read masses from PDB file
421 : }
422 407739 : auto * ptr=chargev[0]->getPntrToAction();
423 407739 : plumed_assert(ptr); // needed for following calls, see #1046
424 407739 : ActionToPutData* cv = ptr->castToActionToPutData();
425 407739 : if(cv) {
426 407737 : chargesWereSet=cv->hasBeenSet();
427 2 : } else if( (chargev[0]->getPntrToAction())->getName()=="CONSTANT" ) {
428 2 : chargesWereSet=true; // Read masses from PDB file
429 : }
430 : unsigned j = 0;
431 :
432 :
433 1218514 : for(const auto & a : atom_value_ind_grouped) {
434 810775 : const auto nn=a.first;
435 810775 : auto & xp=xpos[nn]->data;
436 810775 : auto & yp=ypos[nn]->data;
437 810775 : auto & zp=zpos[nn]->data;
438 810775 : auto & ch=chargev[nn]->data;
439 810775 : auto & ma=masv[nn]->data;
440 11147898 : for(const auto & kk : a.second) {
441 10337123 : actionPositions[j][0] = xp[kk];
442 10337123 : actionPositions[j][1] = yp[kk];
443 10337123 : actionPositions[j][2] = zp[kk];
444 10337123 : charges[j] = ch[kk];
445 10337123 : masses[j] = ma[kk];
446 10337123 : j++;
447 : }
448 : }
449 :
450 : }
451 :
452 428635 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned& ind) {
453 428635 : if( donotforce || (indexes.size()==0 && getName()!="FIXEDATOM") ) {
454 : return;
455 : }
456 676968 : for(unsigned i=0; i<value_depends.size(); ++i) {
457 361139 : xpos[value_depends[i]]->hasForce = true;
458 361139 : ypos[value_depends[i]]->hasForce = true;
459 361139 : zpos[value_depends[i]]->hasForce = true;
460 : }
461 :
462 :
463 987586 : for(const auto & a : atom_value_ind_grouped) {
464 671757 : const auto nn=a.first;
465 : plumed_dbg_assert(ind<forcesToApply.size()) << "problem setting forces in " << getLabel();
466 671757 : auto & xp=xpos[nn]->inputForce;
467 671757 : auto & yp=ypos[nn]->inputForce;
468 671757 : auto & zp=zpos[nn]->inputForce;
469 7155870 : for(const auto & kk : a.second) {
470 6484113 : xp[kk] += forcesToApply[ind];
471 6484113 : ind++;
472 6484113 : yp[kk] += forcesToApply[ind];
473 6484113 : ind++;
474 6484113 : zp[kk] += forcesToApply[ind];
475 6484113 : ind++;
476 : }
477 : }
478 :
479 315829 : setForcesOnCell( forcesToApply, ind );
480 : }
481 :
482 315917 : void ActionAtomistic::setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind) {
483 315917 : setForcesOnCell(forcesToApply.data(),forcesToApply.size(),ind);
484 315917 : }
485 :
486 338924 : void ActionAtomistic::setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind) {
487 3389240 : for(unsigned i=0; i<9; ++i) {
488 : plumed_dbg_massert( ind<size, "problem setting forces in " + getLabel() );
489 3050316 : boxValue->addForce( i, forcesToApply[ind] );
490 3050316 : ind++;
491 : }
492 338924 : }
493 :
494 6 : Tensor ActionAtomistic::getVirial() const {
495 6 : Tensor vir;
496 24 : for(unsigned i=0; i<3; ++i)
497 72 : for(unsigned j=0; j<3; ++j) {
498 54 : vir[i][j] = boxValue->getForce(3*i+j);
499 : }
500 6 : return vir;
501 : }
502 :
503 0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
504 :
505 0 : for(unsigned j=0; j<indexes.size(); j++) {
506 0 : if( indexes[j].index()>pdb.size() ) {
507 0 : error("there are not enough atoms in the input pdb file");
508 : }
509 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) {
510 0 : error("there are atoms missing in the pdb file");
511 : }
512 0 : actionPositions[j]=pdb.getPositions()[indexes[j].index()];
513 : }
514 0 : for(unsigned j=0; j<indexes.size(); j++) {
515 0 : charges[j]=pdb.getBeta()[indexes[j].index()];
516 : }
517 0 : for(unsigned j=0; j<indexes.size(); j++) {
518 0 : masses[j]=pdb.getOccupancy()[indexes[j].index()];
519 : }
520 0 : }
521 :
522 15571 : unsigned ActionAtomistic::getTotAtoms()const {
523 : unsigned natoms = 0;
524 6636103 : for(unsigned i=0; i<xpos.size(); ++i ) {
525 6620532 : natoms += xpos[i]->getNumberOfValues();
526 : }
527 15571 : return natoms;
528 : }
529 :
530 : namespace {
531 : /// matmul(rotation,Vector(x,y,z)) written out for the loops on all the atoms below, where it is several times faster.
532 : /// The sums start from zero and are done in the same order as in matmul, so that the result is identical.
533 19100 : inline void rotateVector( const Tensor& r, const double x, const double y, const double z,
534 : double& rx, double& ry, double& rz ) {
535 19100 : rx=((0.0+r(0,0)*x)+r(0,1)*y)+r(0,2)*z;
536 19100 : ry=((0.0+r(1,0)*x)+r(1,1)*y)+r(1,2)*z;
537 19100 : rz=((0.0+r(2,0)*x)+r(2,1)*y)+r(2,2)*z;
538 19100 : }
539 : }
540 :
541 48 : void ActionAtomistic::translateAllGlobalPositions(const Vector& shift) {
542 96 : for(unsigned j=0; j<xpos.size(); ++j) {
543 48 : auto & xp=xpos[j]->data;
544 48 : auto & yp=ypos[j]->data;
545 48 : auto & zp=zpos[j]->data;
546 48 : const std::size_t n=xpos[j]->getNumberOfValues();
547 6384 : for(std::size_t i=0; i<n; ++i) {
548 6336 : xp[i]+=shift[0];
549 6336 : yp[i]+=shift[1];
550 6336 : zp[i]+=shift[2];
551 : }
552 : }
553 48 : }
554 :
555 17 : void ActionAtomistic::rotateAllGlobalPositions(const Tensor& rotation) {
556 46 : for(unsigned j=0; j<xpos.size(); ++j) {
557 29 : auto & xp=xpos[j]->data;
558 29 : auto & yp=ypos[j]->data;
559 29 : auto & zp=zpos[j]->data;
560 29 : const std::size_t n=xpos[j]->getNumberOfValues();
561 1635 : for(std::size_t i=0; i<n; ++i) {
562 : double rx, ry, rz;
563 1606 : rotateVector( rotation, xp[i], yp[i], zp[i], rx, ry, rz );
564 1606 : xp[i]=rx;
565 1606 : yp[i]=ry;
566 1606 : zp[i]=rz;
567 : }
568 : }
569 17 : }
570 :
571 60 : void ActionAtomistic::rotateAllGlobalPositions(const Tensor& rotation, const Vector& shiftBefore, const Vector& shiftAfter) {
572 144 : for(unsigned j=0; j<xpos.size(); ++j) {
573 84 : auto & xp=xpos[j]->data;
574 84 : auto & yp=ypos[j]->data;
575 84 : auto & zp=zpos[j]->data;
576 84 : const std::size_t n=xpos[j]->getNumberOfValues();
577 8028 : for(std::size_t i=0; i<n; ++i) {
578 : double rx, ry, rz;
579 7944 : rotateVector( rotation, xp[i]-shiftBefore[0], yp[i]-shiftBefore[1], zp[i]-shiftBefore[2], rx, ry, rz );
580 7944 : xp[i]=rx+shiftAfter[0];
581 7944 : yp[i]=ry+shiftAfter[1];
582 7944 : zp[i]=rz+shiftAfter[2];
583 : }
584 : }
585 60 : }
586 :
587 48 : Vector ActionAtomistic::getTotalForceOnAllAtoms() const {
588 48 : Vector totForce;
589 96 : for(unsigned j=0; j<xpos.size(); ++j) {
590 48 : const auto & fx=xpos[j]->inputForce;
591 48 : const auto & fy=ypos[j]->inputForce;
592 48 : const auto & fz=zpos[j]->inputForce;
593 48 : const std::size_t n=xpos[j]->getNumberOfValues();
594 6384 : for(std::size_t i=0; i<n; ++i) {
595 6336 : totForce+=Vector(fx[i],fy[i],fz[i]);
596 : }
597 : }
598 48 : return totForce;
599 : }
600 :
601 77 : Vector ActionAtomistic::rotateForcesOnAllAtoms(const Tensor& rotation) {
602 77 : Vector totForce;
603 190 : for(unsigned j=0; j<xpos.size(); ++j) {
604 113 : auto & fx=xpos[j]->inputForce;
605 113 : auto & fy=ypos[j]->inputForce;
606 113 : auto & fz=zpos[j]->inputForce;
607 113 : const std::size_t n=xpos[j]->getNumberOfValues();
608 113 : if( n>0 ) {
609 113 : xpos[j]->hasForce=true;
610 113 : ypos[j]->hasForce=true;
611 113 : zpos[j]->hasForce=true;
612 : }
613 113 : double tx=totForce[0], ty=totForce[1], tz=totForce[2];
614 9663 : for(std::size_t i=0; i<n; ++i) {
615 : double nx, ny, nz;
616 9550 : rotateVector( rotation, fx[i], fy[i], fz[i], nx, ny, nz );
617 : // the difference is added (as addForce would do) so that the result is identical
618 9550 : fx[i]+=nx-fx[i];
619 9550 : fy[i]+=ny-fy[i];
620 9550 : fz[i]+=nz-fz[i];
621 9550 : tx+=nx;
622 9550 : ty+=ny;
623 9550 : tz+=nz;
624 : }
625 113 : totForce=Vector(tx,ty,tz);
626 : }
627 77 : return totForce;
628 : }
629 :
630 180703 : void ActionAtomistic::makeWhole() {
631 180703 : if(actionMoldat && actionMoldat->isWhole()) {
632 : // make sure the tree has been constructed
633 17 : if(!tree) {
634 12 : tree=std::make_unique<Tree>(actionMoldat);
635 6 : tree->buildTree(getAbsoluteIndexes());
636 : }
637 17 : const auto & tree_indexes=tree->getTreeIndexes();
638 17 : const auto & root_indexes=tree->getRootIndexes();
639 300 : for(unsigned j=0; j<root_indexes.size(); j++) {
640 283 : const Vector & first (actionPositions[root_indexes[j]]);
641 283 : Vector & second (actionPositions[tree_indexes[j+1]]);
642 283 : second=first+pbcDistance(first,second);
643 : }
644 180686 : } else if(actionPositions.size()>1) {
645 1554827 : for(unsigned j=0; j<actionPositions.size()-1; ++j) {
646 : const Vector & first (actionPositions[j]);
647 1374728 : Vector & second (actionPositions[j+1]);
648 1374728 : second=first+pbcDistance(first,second);
649 : }
650 : }
651 180703 : }
652 :
653 15002 : void ActionAtomistic::getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const {
654 15002 : std::size_t nn = atom_value_ind[ind].first;
655 15002 : if( nn==0 ) {
656 14963 : gradients[indexes[ind]] += deriv;
657 14963 : return;
658 : }
659 39 : xpos[nn]->passGradients( deriv[0], gradients );
660 39 : ypos[nn]->passGradients( deriv[1], gradients );
661 39 : zpos[nn]->passGradients( deriv[2], gradients );
662 : }
663 :
664 12048 : void ActionAtomistic::updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l ) {
665 : // The flag is set again when the atoms are requested again or when the local atoms change
666 12048 : unique_local_needs_update=false;
667 12048 : if( useunique ) {
668 6311 : unique_local=unique;
669 6311 : return;
670 : }
671 5737 : unique_local.clear();
672 60739 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
673 55002 : if(g2l[pp->index()]>=0) {
674 23045 : unique_local.push_back(*pp); // already sorted
675 : }
676 : }
677 : }
678 :
679 : }
|
