Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "Atoms.h"
23 : #include "ActionAtomistic.h"
24 : #include "MDAtoms.h"
25 : #include "PlumedMain.h"
26 : #include "tools/Pbc.h"
27 : #include <algorithm>
28 : #include <iostream>
29 : #include <string>
30 : #include <cmath>
31 :
32 : using namespace std;
33 :
34 : namespace PLMD {
35 :
36 : /// We assume that charges and masses are constant along the simulation
37 : /// Set this to false if you want to revert to the original (expensive) behavior
38 : static const bool shareMassAndChargeOnlyAtFirstStep=true;
39 :
40 : class PlumedMain;
41 :
42 2726 : Atoms::Atoms(PlumedMain&plumed):
43 : natoms(0),
44 : md_energy(0.0),
45 : energy(0.0),
46 : dataCanBeSet(false),
47 : collectEnergy(false),
48 : energyHasBeenSet(false),
49 : positionsHaveBeenSet(0),
50 : massesHaveBeenSet(false),
51 : chargesHaveBeenSet(false),
52 : boxHasBeenSet(false),
53 : forcesHaveBeenSet(0),
54 : virialHasBeenSet(false),
55 : massAndChargeOK(false),
56 : shuffledAtoms(0),
57 : mdatoms(MDAtomsBase::create(sizeof(double))),
58 : plumed(plumed),
59 : naturalUnits(false),
60 : MDnaturalUnits(false),
61 : timestep(0.0),
62 : forceOnEnergy(0.0),
63 : zeroallforces(false),
64 : kbT(0.0),
65 : asyncSent(false),
66 : atomsNeeded(false),
67 13630 : ddStep(0)
68 : {
69 2726 : }
70 :
71 8178 : Atoms::~Atoms() {
72 2726 : if(actions.size()>0) {
73 0 : std::cerr<<"WARNING: there is some inconsistency in action added to atoms, as some of them were not properly destroyed. This might indicate an internal bug!!\n";
74 : }
75 2726 : }
76 :
77 85685 : void Atoms::startStep() {
78 85685 : collectEnergy=false; energyHasBeenSet=false; positionsHaveBeenSet=0;
79 85685 : massesHaveBeenSet=false; chargesHaveBeenSet=false; boxHasBeenSet=false;
80 85685 : forcesHaveBeenSet=0; virialHasBeenSet=false; dataCanBeSet=true;
81 85685 : }
82 :
83 38493 : void Atoms::setBox(void*p) {
84 38493 : mdatoms->setBox(p);
85 76986 : Tensor b; mdatoms->getBox(b); boxHasBeenSet=true;
86 38493 : }
87 :
88 42333 : void Atoms::setPositions(void*p) {
89 42333 : plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
90 42350 : plumed_massert( p || gatindex.size()==0, "NULL position pointer with non-zero local atoms");
91 42333 : mdatoms->setp(p); positionsHaveBeenSet=3;
92 42333 : }
93 :
94 42336 : void Atoms::setMasses(void*p) {
95 42342 : plumed_massert( dataCanBeSet,"setMasses must be called after setStep in MD code interface");
96 42350 : plumed_massert( p || gatindex.size()==0, "NULL mass pointer with non-zero local atoms");
97 42333 : mdatoms->setm(p); massesHaveBeenSet=true;
98 :
99 42333 : }
100 :
101 32740 : void Atoms::setCharges(void*p) {
102 32740 : plumed_massert( dataCanBeSet, "setCharges must be called after setStep in MD code interface");
103 32757 : plumed_massert( p || gatindex.size()==0, "NULL charges pointer with non-zero local atoms");
104 32740 : mdatoms->setc(p); chargesHaveBeenSet=true;
105 32740 : }
106 :
107 32835 : void Atoms::setVirial(void*p) {
108 32835 : plumed_massert( dataCanBeSet,"setVirial must be called after setStep in MD code interface");
109 32835 : mdatoms->setVirial(p); virialHasBeenSet=true;
110 32835 : }
111 :
112 9488 : void Atoms::setEnergy(void*p) {
113 9488 : plumed_massert( dataCanBeSet,"setEnergy must be called after setStep in MD code interface");
114 9488 : MD2double(p,md_energy);
115 9488 : md_energy*=MDUnits.getEnergy()/units.getEnergy();
116 9488 : energyHasBeenSet=true;
117 9488 : }
118 :
119 42333 : void Atoms::setForces(void*p) {
120 42333 : plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
121 42350 : plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
122 42333 : forcesHaveBeenSet=3;
123 42333 : mdatoms->setf(p);
124 42333 : }
125 :
126 0 : void Atoms::setPositions(void*p,int i) {
127 0 : plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
128 0 : plumed_massert( p || gatindex.size()==0, "NULL positions pointer with non-zero local atoms");
129 0 : mdatoms->setp(p,i); positionsHaveBeenSet++;
130 0 : }
131 :
132 0 : void Atoms::setForces(void*p,int i) {
133 0 : plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
134 0 : plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
135 0 : mdatoms->setf(p,i); forcesHaveBeenSet++;
136 0 : }
137 :
138 40890 : void Atoms::share() {
139 : // At first step I scatter all the atoms so as to store their mass and charge
140 : // Notice that this works with the assumption that charges and masses are
141 : // not changing during the simulation!
142 40890 : if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
143 675 : shareAll();
144 41565 : return;
145 : }
146 :
147 40215 : if(!(int(gatindex.size())==natoms && shuffledAtoms==0)) {
148 36831 : for(unsigned i=0; i<actions.size(); i++) {
149 34838 : if(actions[i]->isActive()) {
150 11851 : if(!actions[i]->getUnique().empty()) {
151 10406 : atomsNeeded=true;
152 : // unique are the local atoms
153 10406 : unique.insert(actions[i]->getUniqueLocal().begin(),actions[i]->getUniqueLocal().end());
154 : }
155 : }
156 : }
157 : } else {
158 370616 : for(unsigned i=0; i<actions.size(); i++) {
159 332394 : if(actions[i]->isActive()) {
160 123128 : if(!actions[i]->getUnique().empty()) {
161 105349 : atomsNeeded=true;
162 : }
163 : }
164 : }
165 :
166 : }
167 :
168 40215 : share(unique);
169 : }
170 :
171 789 : void Atoms::shareAll() {
172 : unique.clear();
173 : // keep in unique only those atoms that are local
174 789 : if(dd && shuffledAtoms>0) {
175 42203 : for(int i=0; i<natoms; i++) if(g2l[i]>=0) unique.insert(AtomNumber::index(i));
176 : } else {
177 795576 : for(int i=0; i<natoms; i++) unique.insert(AtomNumber::index(i));
178 : }
179 789 : atomsNeeded=true;
180 789 : share(unique);
181 789 : }
182 :
183 41004 : void Atoms::share(const std::set<AtomNumber>& unique) {
184 41004 : plumed_assert( positionsHaveBeenSet==3 && massesHaveBeenSet );
185 :
186 41004 : virial.zero();
187 81912 : if(zeroallforces || int(gatindex.size())==natoms) {
188 7451584 : for(int i=0; i<natoms; i++) forces[i].zero();
189 : } else {
190 86434 : for(const auto & p : unique) forces[p.index()].zero();
191 : }
192 98702 : for(unsigned i=getNatoms(); i<positions.size(); i++) forces[i].zero(); // virtual atoms
193 41004 : forceOnEnergy=0.0;
194 41004 : mdatoms->getBox(box);
195 :
196 82008 : if(!atomsNeeded) return;
197 39772 : atomsNeeded=false;
198 :
199 39772 : if(int(gatindex.size())==natoms && shuffledAtoms==0) {
200 : // faster version, which retrieves all atoms
201 37595 : mdatoms->getPositions(0,natoms,positions);
202 : } else {
203 : uniq_index.clear();
204 2177 : uniq_index.reserve(unique.size());
205 2177 : if(shuffledAtoms>0) {
206 121950 : for(const auto & p : unique) uniq_index.push_back(g2l[p.index()]);
207 : }
208 2177 : mdatoms->getPositions(unique,uniq_index,positions);
209 : }
210 :
211 :
212 : // how many double per atom should be scattered:
213 : int ndata=3;
214 39772 : if(!massAndChargeOK) {
215 : ndata=5;
216 2025 : masses.assign(masses.size(),std::numeric_limits<double>::quiet_NaN());
217 2025 : charges.assign(charges.size(),std::numeric_limits<double>::quiet_NaN());
218 675 : mdatoms->getCharges(gatindex,charges);
219 675 : mdatoms->getMasses(gatindex,masses);
220 : }
221 :
222 39772 : if(dd && shuffledAtoms>0) {
223 2118 : if(dd.async) {
224 16786 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) dd.mpi_request_positions[i].wait();
225 16786 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) dd.mpi_request_index[i].wait();
226 : }
227 2118 : int count=0;
228 30245 : for(const auto & p : unique) {
229 52018 : dd.indexToBeSent[count]=p.index();
230 78027 : dd.positionsToBeSent[ndata*count+0]=positions[p.index()][0];
231 78027 : dd.positionsToBeSent[ndata*count+1]=positions[p.index()][1];
232 78027 : dd.positionsToBeSent[ndata*count+2]=positions[p.index()][2];
233 26009 : if(!massAndChargeOK) {
234 7035 : dd.positionsToBeSent[ndata*count+3]=masses[p.index()];
235 7035 : dd.positionsToBeSent[ndata*count+4]=charges[p.index()];
236 : }
237 26009 : count++;
238 : }
239 2118 : if(dd.async) {
240 2098 : asyncSent=true;
241 2098 : dd.mpi_request_positions.resize(dd.Get_size());
242 2098 : dd.mpi_request_index.resize(dd.Get_size());
243 7540 : for(int i=0; i<dd.Get_size(); i++) {
244 22620 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
245 15080 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
246 : }
247 : } else {
248 20 : const int n=(dd.Get_size());
249 20 : vector<int> counts(n);
250 20 : vector<int> displ(n);
251 20 : vector<int> counts5(n);
252 20 : vector<int> displ5(n);
253 20 : dd.Allgather(count,counts);
254 20 : displ[0]=0;
255 200 : for(int i=1; i<n; ++i) displ[i]=displ[i-1]+counts[i-1];
256 160 : for(int i=0; i<n; ++i) counts5[i]=counts[i]*ndata;
257 160 : for(int i=0; i<n; ++i) displ5[i]=displ[i]*ndata;
258 40 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
259 40 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
260 60 : int tot=displ[n-1]+counts[n-1];
261 1620 : for(int i=0; i<tot; i++) {
262 6400 : positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
263 4800 : positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
264 4800 : positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
265 1600 : if(!massAndChargeOK) {
266 1296 : masses[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+3];
267 1296 : charges[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+4];
268 : }
269 : }
270 : }
271 : }
272 : }
273 :
274 41004 : void Atoms::wait() {
275 41004 : dataCanBeSet=false; // Everything should be set by this stage
276 : // How many double per atom should be scattered
277 : int ndata=3;
278 41004 : if(!massAndChargeOK)ndata=5;
279 :
280 41004 : if(dd) {
281 20039 : dd.Bcast(box,0);
282 : }
283 41004 : pbc.setBox(box);
284 :
285 41004 : if(collectEnergy) energy=md_energy;
286 :
287 41004 : if(dd && shuffledAtoms>0) {
288 : // receive toBeReceived
289 2118 : if(asyncSent) {
290 : Communicator::Status status;
291 : int count=0;
292 7540 : for(int i=0; i<dd.Get_size(); i++) {
293 15080 : dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
294 : int c=status.Get_count<int>();
295 15080 : dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
296 7540 : count+=c;
297 : }
298 60554 : for(int i=0; i<count; i++) {
299 242216 : positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
300 181662 : positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
301 181662 : positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
302 60554 : if(!massAndChargeOK) {
303 17046 : masses[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+3];
304 17046 : charges[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+4];
305 : }
306 : }
307 2098 : asyncSent=false;
308 : }
309 2118 : if(collectEnergy) dd.Sum(energy);
310 : }
311 : // I take note that masses and charges have been set once for all
312 : // at the beginning of the simulation.
313 41004 : if(shareMassAndChargeOnlyAtFirstStep) massAndChargeOK=true;
314 41004 : }
315 :
316 40890 : void Atoms::updateForces() {
317 40890 : plumed_assert( forcesHaveBeenSet==3 );
318 40890 : if(forceOnEnergy*forceOnEnergy>epsilon) {
319 50 : double alpha=1.0-forceOnEnergy;
320 50 : mdatoms->rescaleForces(gatindex,alpha);
321 50 : mdatoms->updateForces(gatindex,forces);
322 : } else {
323 79617 : if(int(gatindex.size())==natoms && shuffledAtoms==0) mdatoms->updateForces(gatindex,forces);
324 2063 : else mdatoms->updateForces(unique,uniq_index,forces);
325 : }
326 40890 : if( !plumed.novirial && dd.Get_rank()==0 ) {
327 26967 : plumed_assert( virialHasBeenSet );
328 26967 : mdatoms->updateVirial(virial);
329 : }
330 40890 : }
331 :
332 771 : void Atoms::setNatoms(int n) {
333 771 : natoms=n;
334 771 : positions.resize(n);
335 771 : forces.resize(n);
336 771 : masses.resize(n);
337 771 : charges.resize(n);
338 771 : gatindex.resize(n);
339 809650 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=i;
340 771 : }
341 :
342 :
343 2819 : void Atoms::add(ActionAtomistic*a) {
344 2819 : actions.push_back(a);
345 2819 : }
346 :
347 2819 : void Atoms::remove(ActionAtomistic*a) {
348 : auto f=find(actions.begin(),actions.end(),a);
349 2819 : plumed_massert(f!=actions.end(),"cannot remove an action registered to atoms");
350 2819 : actions.erase(f);
351 2819 : }
352 :
353 :
354 293 : void Atoms::DomainDecomposition::enable(Communicator& c) {
355 293 : on=true;
356 293 : Set_comm(c.Get_comm());
357 293 : async=Get_size()<10;
358 293 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
359 4 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
360 4 : if(s=="yes") async=true;
361 4 : else if(s=="no") async=false;
362 0 : else plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
363 : }
364 293 : }
365 :
366 1421 : void Atoms::setAtomsNlocal(int n) {
367 1421 : gatindex.resize(n);
368 1421 : g2l.resize(natoms,-1);
369 1421 : if(dd) {
370 : // Since these vectors are sent with MPI by using e.g.
371 : // &dd.positionsToBeSent[0]
372 : // we make sure they are non-zero-sized so as to
373 : // avoid errors when doing boundary check
374 1391 : if(n==0) n++;
375 1391 : dd.positionsToBeSent.resize(n*5,0.0);
376 1391 : dd.positionsToBeReceived.resize(natoms*5,0.0);
377 1391 : dd.indexToBeSent.resize(n,0);
378 1391 : dd.indexToBeReceived.resize(natoms,0);
379 : }
380 1421 : }
381 :
382 976 : void Atoms::setAtomsGatindex(int*g,bool fortran) {
383 981 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
384 1952 : ddStep=plumed.getStep();
385 976 : if(fortran) {
386 0 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=g[i]-1;
387 : } else {
388 41692 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=g[i];
389 : }
390 110892 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
391 976 : if( gatindex.size()==natoms ) {
392 8 : shuffledAtoms=0;
393 872 : for(unsigned i=0; i<gatindex.size(); i++) {
394 864 : if( gatindex[i]!=i ) { shuffledAtoms=1; break; }
395 : }
396 : } else {
397 968 : shuffledAtoms=1;
398 : }
399 976 : if(dd) {
400 946 : dd.Sum(shuffledAtoms);
401 : }
402 62050 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
403 :
404 19592 : for(unsigned i=0; i<actions.size(); i++) {
405 : // keep in unique only those atoms that are local
406 9308 : actions[i]->updateUniqueLocal();
407 : }
408 : unique.clear();
409 976 : }
410 :
411 445 : void Atoms::setAtomsContiguous(int start) {
412 890 : ddStep=plumed.getStep();
413 246636 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=start+i;
414 270815 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
415 369064 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
416 445 : if(gatindex.size()<natoms) shuffledAtoms=1;
417 941 : for(unsigned i=0; i<actions.size(); i++) {
418 : // keep in unique only those atoms that are local
419 248 : actions[i]->updateUniqueLocal();
420 : }
421 : unique.clear();
422 445 : }
423 :
424 664 : void Atoms::setRealPrecision(int p) {
425 1327 : mdatoms=MDAtomsBase::create(p);
426 663 : }
427 :
428 771 : int Atoms::getRealPrecision()const {
429 771 : return mdatoms->getRealPrecision();
430 : }
431 :
432 12119 : void Atoms::MD2double(const void*m,double&d)const {
433 12119 : plumed_assert(mdatoms); mdatoms->MD2double(m,d);
434 12119 : }
435 2303 : void Atoms::double2MD(const double&d,void*m)const {
436 2303 : plumed_assert(mdatoms); mdatoms->double2MD(d,m);
437 2303 : }
438 :
439 771 : void Atoms::updateUnits() {
440 771 : mdatoms->setUnits(units,MDUnits);
441 771 : }
442 :
443 671 : void Atoms::setTimeStep(void*p) {
444 671 : MD2double(p,timestep);
445 671 : }
446 :
447 1431898 : double Atoms::getTimeStep()const {
448 1431898 : return timestep/units.getTime()*MDUnits.getTime();
449 : }
450 :
451 53 : void Atoms::setKbT(void*p) {
452 53 : MD2double(p,kbT);
453 53 : }
454 :
455 977 : double Atoms::getKbT()const {
456 977 : return kbT/units.getEnergy()*MDUnits.getEnergy();
457 : }
458 :
459 :
460 104 : void Atoms::createFullList(int*n) {
461 104 : if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
462 6 : *n=natoms;
463 6 : fullList.resize(natoms);
464 670 : for(unsigned i=0; i<natoms; i++) fullList[i]=i;
465 : } else {
466 : // We update here the unique list defined at Atoms::unique.
467 : // This is not very clear, and probably should be coded differently.
468 : // Hopefully this fix the longstanding issue with NAMD.
469 : unique.clear();
470 1694 : for(unsigned i=0; i<actions.size(); i++) {
471 1498 : if(actions[i]->isActive()) {
472 591 : if(!actions[i]->getUnique().empty()) {
473 523 : atomsNeeded=true;
474 : // unique are the local atoms
475 523 : unique.insert(actions[i]->getUnique().begin(),actions[i]->getUnique().end());
476 : }
477 : }
478 : }
479 98 : fullList.resize(0);
480 98 : fullList.reserve(unique.size());
481 14773 : for(const auto & p : unique) fullList.push_back(p.index());
482 98 : *n=fullList.size();
483 : }
484 104 : }
485 :
486 104 : void Atoms::getFullList(int**x) {
487 104 : if(!fullList.empty()) *x=&fullList[0];
488 6 : else *x=NULL;
489 104 : }
490 :
491 104 : void Atoms::clearFullList() {
492 104 : fullList.resize(0);
493 104 : }
494 :
495 771 : void Atoms::init() {
496 : // Default: set domain decomposition to NO-decomposition, waiting for
497 : // further instruction
498 771 : if(dd) {
499 293 : setAtomsNlocal(natoms);
500 293 : setAtomsContiguous(0);
501 : }
502 771 : }
503 :
504 293 : void Atoms::setDomainDecomposition(Communicator& comm) {
505 293 : dd.enable(comm);
506 293 : }
507 :
508 386 : void Atoms::resizeVectors(unsigned n) {
509 386 : positions.resize(n);
510 386 : forces.resize(n);
511 386 : masses.resize(n);
512 386 : charges.resize(n);
513 386 : }
514 :
515 193 : AtomNumber Atoms::addVirtualAtom(ActionWithVirtualAtom*a) {
516 193 : unsigned n=positions.size();
517 193 : resizeVectors(n+1);
518 193 : virtualAtomsActions.push_back(a);
519 193 : return AtomNumber::index(n);
520 : }
521 :
522 193 : void Atoms::removeVirtualAtom(ActionWithVirtualAtom*a) {
523 193 : unsigned n=positions.size();
524 386 : plumed_massert(a==virtualAtomsActions[virtualAtomsActions.size()-1],"virtual atoms should be destroyed in reverse creation order");
525 193 : resizeVectors(n-1);
526 : virtualAtomsActions.pop_back();
527 193 : }
528 :
529 104 : void Atoms::insertGroup(const std::string&name,const std::vector<AtomNumber>&a) {
530 104 : plumed_massert(groups.count(name)==0,"group named "+name+" already exists");
531 104 : groups[name]=a;
532 104 : }
533 :
534 104 : void Atoms::removeGroup(const std::string&name) {
535 104 : plumed_massert(groups.count(name)==1,"cannot remove group named "+name);
536 : groups.erase(name);
537 104 : }
538 :
539 114 : void Atoms::writeBinary(std::ostream&o)const {
540 228 : o.write(reinterpret_cast<const char*>(&positions[0][0]),natoms*3*sizeof(double));
541 114 : o.write(reinterpret_cast<const char*>(&box(0,0)),9*sizeof(double));
542 114 : o.write(reinterpret_cast<const char*>(&energy),sizeof(double));
543 114 : }
544 :
545 114 : void Atoms::readBinary(std::istream&i) {
546 228 : i.read(reinterpret_cast<char*>(&positions[0][0]),natoms*3*sizeof(double));
547 114 : i.read(reinterpret_cast<char*>(&box(0,0)),9*sizeof(double));
548 114 : i.read(reinterpret_cast<char*>(&energy),sizeof(double));
549 114 : pbc.setBox(box);
550 114 : }
551 :
552 467 : double Atoms::getKBoltzmann()const {
553 467 : if(naturalUnits || MDnaturalUnits) return 1.0;
554 461 : else return kBoltzmann/units.getEnergy();
555 : }
556 :
557 0 : double Atoms::getMDKBoltzmann()const {
558 0 : if(naturalUnits || MDnaturalUnits) return 1.0;
559 0 : else return kBoltzmann/MDUnits.getEnergy();
560 : }
561 :
562 0 : void Atoms::getLocalMasses(std::vector<double>& localMasses) {
563 0 : localMasses.resize(gatindex.size());
564 0 : for(unsigned i=0; i<gatindex.size(); i++) localMasses[i] = masses[gatindex[i]];
565 0 : }
566 :
567 450 : void Atoms::getLocalPositions(std::vector<Vector>& localPositions) {
568 450 : localPositions.resize(gatindex.size());
569 450 : mdatoms->getLocalPositions(localPositions);
570 450 : }
571 :
572 450 : void Atoms::getLocalForces(std::vector<Vector>& localForces) {
573 450 : localForces.resize(gatindex.size());
574 49500 : for(unsigned i=0; i<gatindex.size(); i++) localForces[i] = forces[gatindex[i]];
575 450 : }
576 :
577 0 : void Atoms::getLocalMDForces(std::vector<Vector>& localForces) {
578 0 : localForces.resize(gatindex.size());
579 0 : for(unsigned i=0; i<gatindex.size(); i++) {
580 0 : localForces[i] = mdatoms->getMDforces(i);
581 : }
582 0 : }
583 :
584 20 : void Atoms::setExtraCV(const std::string &name,void*p) {
585 20 : mdatoms->setExtraCV(name,p);
586 20 : }
587 :
588 20 : void Atoms::setExtraCVForce(const std::string &name,void*p) {
589 20 : mdatoms->setExtraCVForce(name,p);
590 20 : }
591 :
592 40 : double Atoms::getExtraCV(const std::string &name) {
593 40 : return mdatoms->getExtraCV(name);
594 : }
595 :
596 40 : void Atoms::updateExtraCVForce(const std::string &name,double f) {
597 40 : mdatoms->updateExtraCVForce(name,f);
598 40 : }
599 :
600 5874 : }
|
