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 18306 : ActionAtomistic::~ActionAtomistic() {
40 : // empty destructor to delete unique_ptr
41 18306 : }
42 :
43 18306 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
44 : Action(ao),
45 18306 : unique_local_needs_update(true),
46 18306 : boxValue(NULL),
47 18306 : lockRequestAtoms(false),
48 18306 : donotretrieve(false),
49 18306 : donotforce(false),
50 18306 : massesWereSet(false),
51 18306 : chargesWereSet(false)
52 : {
53 18306 : ActionWithValue* bv = plumed.getActionSet().selectWithLabel<ActionWithValue*>("Box");
54 18306 : if( bv ) boxValue=bv->copyOutput(0);
55 : // We now get all the information about atoms that are lying about
56 18306 : getAtomValuesFromPlumedObject( plumed, xpos, ypos, zpos, masv, chargev );
57 18306 : if( xpos.size()!=ypos.size() || xpos.size()!=zpos.size() || xpos.size()!=masv.size() || xpos.size()!=chargev.size() )
58 0 : error("mismatch between value arrays");
59 18306 : }
60 :
61 18314 : 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 ) {
62 18314 : std::vector<ActionShortcut*> shortcuts = plumed.getActionSet().select<ActionShortcut*>(); bool foundpdb=false;
63 6737322 : for(const auto & ss : shortcuts ) {
64 6719008 : if( ss->getName()=="READMASSCHARGE" ) {
65 : foundpdb=true;
66 2 : ActionWithValue* mv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_mass");
67 2 : plumed_assert( mv ); masv.push_back( mv->copyOutput(0) );
68 2 : ActionWithValue* qv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_charges");
69 2 : plumed_assert( qv ); chargev.push_back( qv->copyOutput(0) );
70 : }
71 : }
72 18314 : std::vector<ActionWithValue*> vatoms = plumed.getActionSet().select<ActionWithValue*>();
73 7539490 : for(const auto & vv : vatoms ) {
74 7521176 : plumed_assert(vv); // needed for following calls, see #1046
75 7521176 : ActionToPutData* ap = vv->castToActionToPutData();
76 7521176 : if( ap ) {
77 253584 : if( ap->getRole()=="x" ) xpos.push_back( ap->copyOutput(0) );
78 253584 : if( ap->getRole()=="y" ) ypos.push_back( ap->copyOutput(0) );
79 253584 : if( ap->getRole()=="z" ) zpos.push_back( ap->copyOutput(0) );
80 380348 : if( !foundpdb && ap->getRole()=="m" ) masv.push_back( ap->copyOutput(0) );
81 380348 : if( !foundpdb && ap->getRole()=="q" ) chargev.push_back( ap->copyOutput(0) );
82 : }
83 7521176 : ActionWithVirtualAtom* av = vv->castToActionWithVirtualAtom();
84 7521176 : if( av || vv->getName()=="ARGS2VATOM" ) {
85 6711834 : xpos.push_back( vv->copyOutput( vv->getLabel() + ".x") );
86 6711834 : ypos.push_back( vv->copyOutput( vv->getLabel() + ".y") );
87 6711834 : zpos.push_back( vv->copyOutput( vv->getLabel() + ".z") );
88 6711834 : masv.push_back( vv->copyOutput( vv->getLabel() + ".mass") );
89 6711834 : chargev.push_back( vv->copyOutput( vv->getLabel() + ".charge") );
90 : }
91 : }
92 18314 : }
93 :
94 30874 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
95 : (void) keys; // avoid warning
96 30874 : }
97 :
98 :
99 14077 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
100 14077 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
101 14077 : int nat=a.size();
102 14077 : indexes=a;
103 14077 : positions.resize(nat);
104 14077 : forces.resize(nat);
105 14077 : masses.resize(nat);
106 14077 : charges.resize(nat);
107 14077 : atom_value_ind.resize( a.size() );
108 14077 : int n=getTotAtoms();
109 14077 : if(clearDep) clearDependencies();
110 14606 : unique.clear(); std::vector<bool> requirements( xpos.size(), false );
111 14077 : if( boxValue ) addDependency( boxValue->getPntrToAction() );
112 1337488 : for(unsigned i=0; i<indexes.size(); i++) {
113 1323414 : if(indexes[i].index()>=n) { std::string num; Tools::convert( indexes[i].serial(),num ); error("atom " + num + " out of range"); }
114 1323411 : atom_value_ind[i] = getValueIndices( indexes[i] ); requirements[atom_value_ind[i].first] = true;
115 1323411 : if( atom_value_ind[i].first==0 ) unique.push_back(indexes[i]);
116 18024 : else if( atom_value_ind[i].second>0 ) error("action atomistic is not set up to deal with multiple vectors in position input");
117 : }
118 :
119 14076 : atom_value_ind_grouped.clear();
120 :
121 14076 : if(atom_value_ind.size()>0) {
122 14035 : auto nn = atom_value_ind[0].first;
123 14035 : auto kk = atom_value_ind[0].second;
124 14035 : atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
125 14035 : atom_value_ind_grouped.back().second.push_back(kk);
126 : auto prev_nn=nn;
127 1323404 : for(unsigned i=1; i<atom_value_ind.size(); i++) {
128 1309369 : auto nn = atom_value_ind[i].first;
129 1309369 : auto kk = atom_value_ind[i].second;
130 1329052 : if(nn!=prev_nn) atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
131 1309369 : atom_value_ind_grouped.back().second.push_back(kk);
132 : prev_nn=nn;
133 : }
134 : }
135 :
136 : // Add the dependencies to the actions that we require
137 14076 : Tools::removeDuplicates(unique); value_depends.resize(0);
138 6661371 : for(unsigned i=0; i<requirements.size(); ++i ) {
139 6647295 : if( !requirements[i] ) continue;
140 30783 : value_depends.push_back( i );
141 30782 : addDependency( xpos[i]->getPntrToAction() );
142 30782 : addDependency( ypos[i]->getPntrToAction() );
143 30782 : addDependency( zpos[i]->getPntrToAction() );
144 30782 : addDependency( masv[i]->getPntrToAction() );
145 30782 : addDependency( chargev[i]->getPntrToAction() );
146 : }
147 14076 : unique_local_needs_update=true;
148 14076 : }
149 :
150 405876 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
151 405876 : pbc.apply(dlist, max_index);
152 405876 : }
153 :
154 156 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
155 156 : calculateAtomicNumericalDerivatives( a, 0 );
156 156 : }
157 :
158 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
159 0 : pbc.setBox( newbox );
160 0 : }
161 :
162 240 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
163 240 : if(!a) {
164 240 : a=castToActionWithValue();
165 240 : plumed_massert(a,"only Actions with a value can be differentiated");
166 : }
167 :
168 240 : const size_t nval=a->getNumberOfComponents();
169 240 : const size_t natoms=getNumberOfAtoms();
170 240 : std::vector<Vector> value(nval*natoms);
171 240 : std::vector<Tensor> valuebox(nval);
172 240 : std::vector<Vector> savedPositions(natoms);
173 : const double delta=std::sqrt(epsilon);
174 :
175 82068 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
176 61371 : savedPositions[i][k]=positions[i][k];
177 61371 : positions[i][k]=positions[i][k]+delta;
178 61371 : a->calculate();
179 61371 : positions[i][k]=savedPositions[i][k];
180 188208 : for(int j=0; j<nval; j++) {
181 126837 : value[j*natoms+i][k]=a->getOutputQuantity(j);
182 : }
183 : }
184 240 : Tensor box(pbc.getBox());
185 3120 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
186 2160 : double arg0=box(i,k);
187 186273 : for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
188 2160 : box(i,k)=box(i,k)+delta;
189 2160 : pbc.setBox(box);
190 186273 : for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
191 2160 : a->calculate();
192 2160 : box(i,k)=arg0;
193 2160 : pbc.setBox(box);
194 186273 : for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
195 14886 : for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
196 : }
197 :
198 240 : a->calculate();
199 240 : a->clearDerivatives();
200 1654 : for(int j=0; j<nval; j++) {
201 1414 : Value* v=a->copyOutput(j);
202 1414 : double ref=v->get();
203 1414 : if(v->hasDerivatives()) {
204 89779 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
205 66975 : double d=(value[j*natoms+i][k]-ref)/delta;
206 66975 : v->addDerivative(startnum+3*i+k,d);
207 : }
208 479 : Tensor virial;
209 6227 : for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
210 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
211 479 : virial=-matmul(box.transpose(),virial);
212 6227 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
213 : }
214 : }
215 240 : }
216 :
217 105422 : bool ActionAtomistic::actionHasForces() {
218 105422 : ActionWithValue* av = castToActionWithValue();
219 105422 : if( av ) return !av->doNotCalculateDerivatives();
220 0 : if( indexes.size()>0 ) plumed_merror("you have to overwrite the function actionHasForce to tell plumed if you method applies forces");
221 : return true;
222 : }
223 :
224 13922 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
225 13922 : parseAtomList(key,-1,t);
226 13921 : }
227 :
228 51540 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
229 51540 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
230 : std::vector<std::string> strings;
231 51540 : if( num<0 ) {
232 17481 : parseVector(key,strings);
233 17480 : if(strings.empty()) return;
234 : } else {
235 34059 : if ( !parseNumberedVector(key,num,strings) ) return;
236 : }
237 47363 : t.resize(0); interpretAtomList( strings, xpos, this, t );
238 51540 : }
239 :
240 28 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
241 28 : interpretAtomList( strings, xpos, this, t );
242 28 : }
243 :
244 47955 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t) {
245 47955 : Tools::interpretRanges(strings);
246 1159863 : for(unsigned i=0; i<strings.size(); ++i) {
247 : AtomNumber atom;
248 1111908 : bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
249 1111908 : if(ok) t.push_back(atom);
250 : // here we check if this is a special symbol for MOLINFO
251 1111908 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
252 705 : std::string symbol=strings[i].substr(1);
253 705 : if(symbol=="allatoms") {
254 24 : auto n=0; for(unsigned i=0; i<xpos.size(); ++i) n += xpos[i]->getNumberOfValues();
255 365 : t.reserve(n); for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
256 : ok=true;
257 695 : } else if(symbol=="mdatoms") {
258 14 : const auto n=xpos[0]->getNumberOfValues();
259 14 : t.reserve(t.size()+n);
260 1166 : for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
261 : ok=true;
262 1362 : } else if(Tools::startWith(symbol,"ndx:")) {
263 40 : auto words=Tools::getWords(symbol.substr(4));
264 : std::string ndxfile,ndxgroup;
265 20 : if(words.size()==1) {
266 : ndxfile=words[0];
267 18 : } else if(words.size()==2) {
268 : ndxfile=words[0];
269 : ndxgroup=words[1];
270 0 : } else plumed_error()<<"Cannot intepret selection "<<symbol;
271 :
272 38 : if(ndxgroup.size()>0) action->log<<" importing group '"+ndxgroup+"'";
273 2 : else action->log<<" importing first group";
274 20 : action->log<<" from index file "<<ndxfile<<"\n";
275 :
276 20 : IFile ifile;
277 20 : ifile.open(ndxfile);
278 : std::string line;
279 : std::string groupname;
280 : bool firstgroup=true;
281 : bool groupfound=false;
282 3849 : while(ifile.getline(line)) {
283 3829 : std::vector<std::string> words=Tools::getWords(line);
284 7765 : if(words.size()>=3 && words[0]=="[" && words[2]=="]") {
285 168 : if(groupname.length()>0) firstgroup=false;
286 : groupname=words[1];
287 168 : if(groupname==ndxgroup || ndxgroup.length()==0) groupfound=true;
288 3661 : } else if(groupname==ndxgroup || (firstgroup && ndxgroup.length()==0)) {
289 6186 : for(unsigned i=0; i<words.size(); i++) {
290 5782 : AtomNumber at; Tools::convert(words[i],at);
291 5782 : t.push_back(at);
292 : }
293 : }
294 3829 : }
295 20 : if(!groupfound) plumed_error()<<"group has not been found in index file";
296 : ok=true;
297 40 : } else {
298 661 : auto* moldat=action->plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
299 661 : if( moldat ) {
300 661 : std::vector<AtomNumber> atom_list; moldat->interpretSymbol( symbol, atom_list );
301 661 : if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
302 0 : else { action->error(strings[i] + " is not a label plumed knows"); }
303 : } else {
304 0 : action->error("atoms specified using @ symbol but no MOLINFO was available");
305 : }
306 : }
307 : }
308 : // here we check if the atom name is the name of a group
309 1111203 : if(!ok) {
310 24226 : Group* mygrp=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]);
311 24226 : if(mygrp) {
312 445 : std::vector<std::string> grp_str( mygrp->getGroupAtoms() );
313 445 : interpretAtomList( grp_str, xpos, action, t ); ok=true;
314 445 : } else {
315 23781 : Group* mygrp2=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]+"_grp");
316 23781 : if(mygrp2) {
317 111 : std::vector<std::string> grp_str( mygrp2->getGroupAtoms() );
318 111 : interpretAtomList( grp_str, xpos, action, t ); ok=true;
319 111 : }
320 : }
321 : }
322 : // here we check if the atom name is the name of an added virtual atom
323 1111352 : if(!ok) {
324 : unsigned ind = 0;
325 13269158 : for(unsigned j=0; j<xpos.size(); ++j) {
326 13269158 : if( xpos[j]->getPntrToAction()->getLabel()==strings[i] ) {
327 23670 : t.push_back( AtomNumber::index(ind) ); ok=true; break;
328 : }
329 13245488 : ind = ind + xpos[j]->getNumberOfValues();
330 : }
331 : }
332 1088238 : if(!ok) action->error("it was not possible to interpret atom name " + strings[i]);
333 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
334 : }
335 47955 : }
336 :
337 1761306 : std::pair<std::size_t, std::size_t> ActionAtomistic::getValueIndices( const AtomNumber& i ) const {
338 1761306 : std::size_t valno=0, k = i.index();
339 15486520 : for(unsigned j=0; j<xpos.size(); ++j) {
340 15486520 : if( k<xpos[j]->getNumberOfValues() ) { valno=j; break; }
341 13725214 : k = k - xpos[j]->getNumberOfValues();
342 : }
343 1761306 : return std::pair<std::size_t, std::size_t>( valno, k );
344 : }
345 :
346 542621 : void ActionAtomistic::retrieveAtoms( const bool& force ) {
347 542621 : if( boxValue ) {
348 526106 : auto* ptr=boxValue->getPntrToAction();
349 526106 : plumed_assert(ptr); // needed for following calls, see #1046
350 526106 : PbcAction* pbca = ptr->castToPbcAction();
351 526106 : plumed_assert( pbca ); pbc=pbca->pbc;
352 : }
353 542621 : if( donotretrieve || indexes.size()==0 ) return;
354 225004 : auto * mtr=masv[0]->getPntrToAction();
355 225004 : plumed_assert(mtr); // needed for following calls, see #1046
356 225004 : ActionToPutData* mv = mtr->castToActionToPutData();
357 225004 : if(mv) massesWereSet=mv->hasBeenSet();
358 2 : else if( (masv[0]->getPntrToAction())->getName()=="CONSTANT" ) massesWereSet=true; // Read masses from PDB file
359 225004 : auto * ptr=chargev[0]->getPntrToAction();
360 225004 : plumed_assert(ptr); // needed for following calls, see #1046
361 225004 : ActionToPutData* cv = ptr->castToActionToPutData();
362 225004 : if(cv) chargesWereSet=cv->hasBeenSet();
363 2 : else if( (chargev[0]->getPntrToAction())->getName()=="CONSTANT" ) chargesWereSet=true; // Read masses from PDB file
364 : unsigned j = 0;
365 :
366 : // for(const auto & a : atom_value_ind) {
367 : // std::size_t nn = a.first, kk = a.second;
368 : // positions[j][0] = xpos[nn]->data[kk];
369 : // positions[j][1] = ypos[nn]->data[kk];
370 : // positions[j][2] = zpos[nn]->data[kk];
371 : // charges[j] = chargev[nn]->data[kk];
372 : // masses[j] = masv[nn]->data[kk];
373 : // j++;
374 : // }
375 :
376 840649 : for(const auto & a : atom_value_ind_grouped) {
377 615645 : const auto nn=a.first;
378 615645 : auto & xp=xpos[nn]->data;
379 615645 : auto & yp=ypos[nn]->data;
380 615645 : auto & zp=zpos[nn]->data;
381 615645 : auto & ch=chargev[nn]->data;
382 615645 : auto & ma=masv[nn]->data;
383 6608291 : for(const auto & kk : a.second) {
384 5992646 : positions[j][0] = xp[kk];
385 5992646 : positions[j][1] = yp[kk];
386 5992646 : positions[j][2] = zp[kk];
387 5992646 : charges[j] = ch[kk];
388 5992646 : masses[j] = ma[kk];
389 5992646 : j++;
390 : }
391 : }
392 :
393 : }
394 :
395 247404 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned& ind) {
396 247404 : if( donotforce || (indexes.size()==0 && getName()!="FIXEDATOM") ) return;
397 486697 : for(unsigned i=0; i<value_depends.size(); ++i) {
398 343121 : xpos[value_depends[i]]->hasForce = true;
399 343121 : ypos[value_depends[i]]->hasForce = true;
400 343121 : zpos[value_depends[i]]->hasForce = true;
401 : }
402 :
403 : // for(const auto & a : atom_value_ind) {
404 : // plumed_dbg_massert( ind<forcesToApply.size(), "problem setting forces in " + getLabel() );
405 : // std::size_t nn = a.first, kk = a.second;
406 : // xpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
407 : // ypos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
408 : // zpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
409 : // }
410 :
411 797085 : for(const auto & a : atom_value_ind_grouped) {
412 653509 : const auto nn=a.first;
413 : plumed_dbg_assert(ind<forcesToApply.size()) << "problem setting forces in " << getLabel();
414 653509 : auto & xp=xpos[nn]->inputForce;
415 653509 : auto & yp=ypos[nn]->inputForce;
416 653509 : auto & zp=zpos[nn]->inputForce;
417 3980595 : for(const auto & kk : a.second) {
418 3327086 : xp[kk] += forcesToApply[ind]; ind++;
419 3327086 : yp[kk] += forcesToApply[ind]; ind++;
420 3327086 : zp[kk] += forcesToApply[ind]; ind++;
421 : }
422 : }
423 :
424 143576 : setForcesOnCell( forcesToApply, ind );
425 : }
426 :
427 143664 : void ActionAtomistic::setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind) {
428 143664 : setForcesOnCell(forcesToApply.data(),forcesToApply.size(),ind);
429 143664 : }
430 :
431 166565 : void ActionAtomistic::setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind) {
432 1665650 : for(unsigned i=0; i<9; ++i) {
433 : plumed_dbg_massert( ind<size, "problem setting forces in " + getLabel() );
434 1499085 : boxValue->addForce( i, forcesToApply[ind] ); ind++;
435 : }
436 166565 : }
437 :
438 6 : Tensor ActionAtomistic::getVirial() const {
439 78 : Tensor vir; for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) vir[i][j] = boxValue->getForce(3*i+j);
440 6 : return vir;
441 : }
442 :
443 0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
444 :
445 0 : for(unsigned j=0; j<indexes.size(); j++) {
446 0 : if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
447 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
448 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
449 : }
450 0 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
451 0 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
452 0 : }
453 :
454 14340 : unsigned ActionAtomistic::getTotAtoms()const {
455 : unsigned natoms = 0;
456 6661971 : for(unsigned i=0; i<xpos.size(); ++i ) natoms += xpos[i]->getNumberOfValues();
457 14340 : return natoms;
458 : }
459 :
460 149867 : void ActionAtomistic::makeWhole() {
461 1375530 : for(unsigned j=0; j<positions.size()-1; ++j) {
462 : const Vector & first (positions[j]);
463 1225663 : Vector & second (positions[j+1]);
464 1225663 : second=first+pbcDistance(first,second);
465 : }
466 149867 : }
467 :
468 8208 : void ActionAtomistic::getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const {
469 8208 : std::size_t nn = atom_value_ind[ind].first;
470 8208 : if( nn==0 ) { gradients[indexes[ind]] += deriv; return; }
471 39 : xpos[nn]->passGradients( deriv[0], gradients );
472 39 : ypos[nn]->passGradients( deriv[1], gradients );
473 39 : zpos[nn]->passGradients( deriv[2], gradients );
474 : }
475 :
476 256390 : void ActionAtomistic::updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l ) {
477 256390 : if( useunique ) { unique_local=unique; return; }
478 : // Update unique local if it needs an update
479 9154 : unique_local_needs_update=false; unique_local.clear();
480 75612 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
481 66458 : if(g2l[pp->index()]>=0) unique_local.push_back(*pp); // already sorted
482 : }
483 : }
484 :
485 : }
|