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 "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 : namespace PLMD {
33 :
34 : /// We assume that charges and masses are constant along the simulation
35 : /// Set this to false if you want to revert to the original (expensive) behavior
36 : static const bool shareMassAndChargeOnlyAtFirstStep=true;
37 :
38 : class PlumedMain;
39 :
40 404395 : Atoms::Atoms(PlumedMain&plumed):
41 404395 : natoms(0),
42 404395 : md_energy(0.0),
43 404395 : energy(0.0),
44 404395 : dataCanBeSet(false),
45 404395 : collectEnergy(false),
46 404395 : energyHasBeenSet(false),
47 404395 : positionsHaveBeenSet(0),
48 404395 : massesHaveBeenSet(false),
49 404395 : chargesHaveBeenSet(false),
50 404395 : boxHasBeenSet(false),
51 404395 : forcesHaveBeenSet(0),
52 404395 : virialHasBeenSet(false),
53 404395 : massAndChargeOK(false),
54 404395 : shuffledAtoms(0),
55 404395 : mdatoms(MDAtomsBase::create(sizeof(double))),
56 404395 : plumed(plumed),
57 404395 : naturalUnits(false),
58 404395 : MDnaturalUnits(false),
59 404395 : timestep(0.0),
60 404395 : forceOnEnergy(0.0),
61 404395 : zeroallforces(false),
62 404395 : kbT(0.0),
63 404395 : asyncSent(false),
64 404395 : atomsNeeded(false),
65 808790 : ddStep(0)
66 : {
67 404395 : }
68 :
69 404395 : Atoms::~Atoms() {
70 404395 : if(actions.size()>0) {
71 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";
72 : }
73 808790 : }
74 :
75 247830 : void Atoms::startStep() {
76 247830 : collectEnergy=false; energyHasBeenSet=false; positionsHaveBeenSet=0;
77 247830 : massesHaveBeenSet=false; chargesHaveBeenSet=false; boxHasBeenSet=false;
78 247830 : forcesHaveBeenSet=0; virialHasBeenSet=false; dataCanBeSet=true;
79 247830 : }
80 :
81 43333 : void Atoms::setBox(const TypesafePtr & p) {
82 43333 : mdatoms->setBox(p);
83 43333 : Tensor b; mdatoms->getBox(b); boxHasBeenSet=true;
84 43333 : }
85 :
86 47474 : void Atoms::setPositions(const TypesafePtr & p) {
87 47474 : plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
88 47474 : plumed_massert( p || gatindex.size()==0, "NULL position pointer with non-zero local atoms");
89 47474 : mdatoms->setp(p); positionsHaveBeenSet=3;
90 47474 : }
91 :
92 47478 : void Atoms::setMasses(const TypesafePtr & p) {
93 47481 : plumed_massert( dataCanBeSet,"setMasses must be called after setStep in MD code interface");
94 47475 : plumed_massert( p || gatindex.size()==0, "NULL mass pointer with non-zero local atoms");
95 47475 : mdatoms->setm(p); massesHaveBeenSet=true;
96 :
97 47475 : }
98 :
99 37570 : void Atoms::setCharges(const TypesafePtr & p) {
100 37570 : plumed_massert( dataCanBeSet, "setCharges must be called after setStep in MD code interface");
101 37570 : plumed_massert( p || gatindex.size()==0, "NULL charges pointer with non-zero local atoms");
102 37570 : mdatoms->setc(p); chargesHaveBeenSet=true;
103 37570 : }
104 :
105 37675 : void Atoms::setVirial(const TypesafePtr & p) {
106 37675 : plumed_massert( dataCanBeSet,"setVirial must be called after setStep in MD code interface");
107 37675 : mdatoms->setVirial(p); virialHasBeenSet=true;
108 37673 : }
109 :
110 9792 : void Atoms::setEnergy(const TypesafePtr & p) {
111 9792 : plumed_massert( dataCanBeSet,"setEnergy must be called after setStep in MD code interface");
112 9792 : MD2double(p,md_energy);
113 9792 : md_energy*=MDUnits.getEnergy()/units.getEnergy();
114 9792 : energyHasBeenSet=true;
115 9792 : }
116 :
117 47474 : void Atoms::setForces(const TypesafePtr & p) {
118 47474 : plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
119 47474 : plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
120 47474 : forcesHaveBeenSet=3;
121 47474 : mdatoms->setf(p);
122 47474 : }
123 :
124 3 : void Atoms::setPositions(const TypesafePtr & p,int i) {
125 3 : plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
126 3 : plumed_massert( p || gatindex.size()==0, "NULL positions pointer with non-zero local atoms");
127 3 : mdatoms->setp(p,i); positionsHaveBeenSet++;
128 3 : }
129 :
130 3 : void Atoms::setForces(const TypesafePtr & p,int i) {
131 3 : plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
132 3 : plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
133 3 : mdatoms->setf(p,i); forcesHaveBeenSet++;
134 3 : }
135 :
136 46042 : void Atoms::share() {
137 : // At first step I scatter all the atoms so as to store their mass and charge
138 : // Notice that this works with the assumption that charges and masses are
139 : // not changing during the simulation!
140 46042 : if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
141 784 : shareAll();
142 784 : return;
143 : }
144 :
145 45258 : if(!(int(gatindex.size())==natoms && shuffledAtoms==0)) {
146 19665 : for(unsigned i=0; i<actions.size(); i++) {
147 17639 : if(actions[i]->isActive()) {
148 12107 : if(!actions[i]->getUnique().empty()) {
149 10629 : atomsNeeded=true;
150 : // unique are the local atoms
151 10629 : unique.insert(actions[i]->getUniqueLocal().begin(),actions[i]->getUniqueLocal().end());
152 : }
153 : }
154 : }
155 : } else {
156 220789 : for(unsigned i=0; i<actions.size(); i++) {
157 177557 : if(actions[i]->isActive()) {
158 133775 : if(!actions[i]->getUnique().empty()) {
159 115117 : atomsNeeded=true;
160 : }
161 : }
162 : }
163 :
164 : }
165 :
166 45258 : share(unique);
167 : }
168 :
169 898 : void Atoms::shareAll() {
170 : unique.clear();
171 : // keep in unique only those atoms that are local
172 898 : if(dd && shuffledAtoms>0) {
173 17616 : for(int i=0; i<natoms; i++) if(g2l[i]>=0) unique.insert(AtomNumber::index(i));
174 : } else {
175 477612 : for(int i=0; i<natoms; i++) unique.insert(AtomNumber::index(i));
176 : }
177 898 : atomsNeeded=true;
178 898 : share(unique);
179 898 : }
180 :
181 46156 : void Atoms::share(const std::set<AtomNumber>& unique) {
182 46156 : plumed_assert( positionsHaveBeenSet==3 && massesHaveBeenSet );
183 :
184 46156 : virial.zero();
185 46156 : if(zeroallforces || int(gatindex.size())==natoms) {
186 4633356 : for(int i=0; i<natoms; i++) forces[i].zero();
187 : } else {
188 31126 : for(const auto & p : unique) forces[p.index()].zero();
189 : }
190 61703 : for(unsigned i=getNatoms(); i<positions.size(); i++) forces[i].zero(); // virtual atoms
191 46156 : forceOnEnergy=0.0;
192 46156 : mdatoms->getBox(box);
193 :
194 46156 : if(!atomsNeeded) return;
195 44914 : atomsNeeded=false;
196 :
197 44914 : if(int(gatindex.size())==natoms && shuffledAtoms==0) {
198 : // faster version, which retrieves all atoms
199 42702 : mdatoms->getPositions(0,natoms,positions);
200 : } else {
201 2212 : uniq_index.clear();
202 2212 : uniq_index.reserve(unique.size());
203 2212 : if(shuffledAtoms>0) {
204 33181 : for(const auto & p : unique) uniq_index.push_back(g2l[p.index()]);
205 : }
206 2212 : mdatoms->getPositions(unique,uniq_index,positions);
207 : }
208 :
209 :
210 : // how many double per atom should be scattered:
211 : int ndata=3;
212 44914 : if(!massAndChargeOK) {
213 : ndata=5;
214 784 : masses.assign(masses.size(),std::numeric_limits<double>::quiet_NaN());
215 784 : charges.assign(charges.size(),std::numeric_limits<double>::quiet_NaN());
216 784 : mdatoms->getCharges(gatindex,charges);
217 784 : mdatoms->getMasses(gatindex,masses);
218 : }
219 :
220 44914 : if(dd && shuffledAtoms>0) {
221 2142 : if(dd.async) {
222 9510 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) dd.mpi_request_positions[i].wait();
223 9510 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) dd.mpi_request_index[i].wait();
224 : }
225 2142 : int count=0;
226 29699 : for(const auto & p : unique) {
227 27557 : dd.indexToBeSent[count]=p.index();
228 27557 : dd.positionsToBeSent[ndata*count+0]=positions[p.index()][0];
229 27557 : dd.positionsToBeSent[ndata*count+1]=positions[p.index()][1];
230 27557 : dd.positionsToBeSent[ndata*count+2]=positions[p.index()][2];
231 27557 : if(!massAndChargeOK) {
232 2477 : dd.positionsToBeSent[ndata*count+3]=masses[p.index()];
233 2477 : dd.positionsToBeSent[ndata*count+4]=charges[p.index()];
234 : }
235 27557 : count++;
236 : }
237 2142 : if(dd.async) {
238 2122 : asyncSent=true;
239 2122 : dd.mpi_request_positions.resize(dd.Get_size());
240 2122 : dd.mpi_request_index.resize(dd.Get_size());
241 9710 : for(int i=0; i<dd.Get_size(); i++) {
242 7588 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
243 7588 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
244 : }
245 : } else {
246 20 : const int n=(dd.Get_size());
247 20 : std::vector<int> counts(n);
248 20 : std::vector<int> displ(n);
249 20 : std::vector<int> counts5(n);
250 20 : std::vector<int> displ5(n);
251 20 : dd.Allgather(count,counts);
252 20 : displ[0]=0;
253 80 : for(int i=1; i<n; ++i) displ[i]=displ[i-1]+counts[i-1];
254 100 : for(int i=0; i<n; ++i) counts5[i]=counts[i]*ndata;
255 100 : for(int i=0; i<n; ++i) displ5[i]=displ[i]*ndata;
256 20 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
257 20 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
258 20 : int tot=displ[n-1]+counts[n-1];
259 1620 : for(int i=0; i<tot; i++) {
260 1600 : positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
261 1600 : positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
262 1600 : positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
263 1600 : if(!massAndChargeOK) {
264 432 : masses[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+3];
265 432 : charges[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+4];
266 : }
267 : }
268 : }
269 : }
270 : }
271 :
272 46156 : void Atoms::wait() {
273 46156 : dataCanBeSet=false; // Everything should be set by this stage
274 : // How many double per atom should be scattered
275 : std::size_t ndata=3;
276 46156 : if(!massAndChargeOK)ndata=5;
277 :
278 46156 : if(dd) {
279 21591 : dd.Bcast(box,0);
280 : }
281 46156 : pbc.setBox(box);
282 :
283 46156 : if(collectEnergy) energy=md_energy;
284 :
285 46156 : if(dd && shuffledAtoms>0) {
286 : // receive toBeReceived
287 2142 : if(asyncSent) {
288 : Communicator::Status status;
289 : std::size_t count=0;
290 9710 : for(int i=0; i<dd.Get_size(); i++) {
291 7588 : dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
292 7588 : int c=status.Get_count<int>();
293 7588 : dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
294 7588 : count+=c;
295 : }
296 65772 : for(int i=0; i<count; i++) {
297 63650 : positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
298 63650 : positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
299 63650 : positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
300 63650 : if(!massAndChargeOK) {
301 5946 : masses[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+3];
302 5946 : charges[dd.indexToBeReceived[i]] =dd.positionsToBeReceived[ndata*i+4];
303 : }
304 : }
305 2122 : asyncSent=false;
306 : }
307 2142 : if(collectEnergy) dd.Sum(energy);
308 : }
309 : // I take note that masses and charges have been set once for all
310 : // at the beginning of the simulation.
311 46156 : if(shareMassAndChargeOnlyAtFirstStep) massAndChargeOK=true;
312 46156 : }
313 :
314 46032 : void Atoms::updateForces() {
315 46032 : plumed_assert( forcesHaveBeenSet==3 );
316 46032 : if(forceOnEnergy*forceOnEnergy>epsilon) {
317 50 : double alpha=1.0-forceOnEnergy;
318 50 : mdatoms->rescaleForces(gatindex,alpha);
319 50 : mdatoms->updateForces(gatindex,forces);
320 : } else {
321 45982 : if(int(gatindex.size())==natoms && shuffledAtoms==0) mdatoms->updateForces(gatindex,forces);
322 2098 : else mdatoms->updateForces(unique,uniq_index,forces);
323 : }
324 46032 : if( !plumed.novirial && dd.Get_rank()==0 ) {
325 31465 : plumed_assert( virialHasBeenSet );
326 31465 : mdatoms->updateVirial(virial);
327 : }
328 46032 : }
329 :
330 893 : void Atoms::setNatoms(int n) {
331 893 : natoms=n;
332 893 : positions.resize(n);
333 893 : forces.resize(n);
334 893 : masses.resize(n);
335 893 : charges.resize(n);
336 893 : gatindex.resize(n);
337 485359 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=i;
338 893 : }
339 :
340 :
341 10109 : void Atoms::add(ActionAtomistic*a) {
342 10109 : actions.push_back(a);
343 10109 : }
344 :
345 10109 : void Atoms::remove(ActionAtomistic*a) {
346 10109 : auto f=find(actions.begin(),actions.end(),a);
347 10109 : plumed_massert(f!=actions.end(),"cannot remove an action registered to atoms");
348 10109 : actions.erase(f);
349 10109 : }
350 :
351 :
352 337 : void Atoms::DomainDecomposition::enable(Communicator& c) {
353 337 : on=true;
354 337 : Set_comm(c.Get_comm());
355 337 : async=Get_size()<10;
356 337 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
357 4 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
358 4 : if(s=="yes") async=true;
359 4 : else if(s=="no") async=false;
360 0 : else plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
361 : }
362 337 : }
363 :
364 1477 : void Atoms::setAtomsNlocal(int n) {
365 1477 : gatindex.resize(n);
366 1477 : g2l.resize(natoms,-1);
367 1477 : if(dd) {
368 : // Since these vectors are sent with MPI by using e.g.
369 : // &dd.positionsToBeSent[0]
370 : // we make sure they are non-zero-sized so as to
371 : // avoid errors when doing boundary check
372 1445 : if(n==0) n++;
373 1445 : dd.positionsToBeSent.resize(n*5,0.0);
374 1445 : dd.positionsToBeReceived.resize(natoms*5,0.0);
375 1445 : dd.indexToBeSent.resize(n,0);
376 1445 : dd.indexToBeReceived.resize(natoms,0);
377 : }
378 1477 : }
379 :
380 988 : void Atoms::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
381 988 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
382 : auto gg=g.get<const int*>(gatindex.size());
383 988 : ddStep=plumed.getStep();
384 988 : if(fortran) {
385 0 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=gg[i]-1;
386 : } else {
387 22140 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=gg[i];
388 : }
389 57530 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
390 988 : if( gatindex.size()==natoms ) {
391 9 : shuffledAtoms=0;
392 1005 : for(unsigned i=0; i<gatindex.size(); i++) {
393 996 : if( gatindex[i]!=i ) { shuffledAtoms=1; break; }
394 : }
395 : } else {
396 979 : shuffledAtoms=1;
397 : }
398 988 : if(dd) {
399 956 : dd.Sum(shuffledAtoms);
400 : }
401 22140 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
402 :
403 10390 : for(unsigned i=0; i<actions.size(); i++) {
404 : // keep in unique only those atoms that are local
405 9402 : actions[i]->updateUniqueLocal();
406 : }
407 : unique.clear();
408 988 : }
409 :
410 489 : void Atoms::setAtomsContiguous(int start) {
411 489 : ddStep=plumed.getStep();
412 184532 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=start+i;
413 196844 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
414 184532 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
415 489 : if(gatindex.size()<natoms) shuffledAtoms=1;
416 737 : for(unsigned i=0; i<actions.size(); i++) {
417 : // keep in unique only those atoms that are local
418 248 : actions[i]->updateUniqueLocal();
419 : }
420 : unique.clear();
421 489 : }
422 :
423 780 : void Atoms::setRealPrecision(int p) {
424 1560 : mdatoms=MDAtomsBase::create(p);
425 780 : }
426 :
427 1688 : int Atoms::getRealPrecision()const {
428 1688 : return mdatoms->getRealPrecision();
429 : }
430 :
431 12830 : void Atoms::MD2double(const TypesafePtr & m,double&d)const {
432 12830 : plumed_assert(mdatoms); mdatoms->MD2double(m,d);
433 12830 : }
434 2365 : void Atoms::double2MD(const double&d,const TypesafePtr & m)const {
435 2365 : plumed_assert(mdatoms); mdatoms->double2MD(d,m);
436 2365 : }
437 :
438 933 : void Atoms::updateUnits() {
439 933 : mdatoms->setUnits(units,MDUnits);
440 933 : }
441 :
442 774 : void Atoms::setTimeStep(const TypesafePtr & p) {
443 774 : MD2double(p,timestep);
444 : // The following is to avoid extra digits in case the MD code uses floats
445 : // e.g.: float f=0.002 when converted to double becomes 0.002000000094995
446 : // To avoid this, we keep only up to 6 significant digits after first one
447 774 : if(getRealPrecision()<=4) {
448 0 : double magnitude=std::pow(10,std::floor(std::log10(timestep)));
449 0 : timestep=std::round(timestep/magnitude*1e6)/1e6*magnitude;
450 : }
451 774 : }
452 :
453 3158863 : double Atoms::getTimeStep()const {
454 3158863 : return timestep/units.getTime()*MDUnits.getTime();
455 : }
456 :
457 59 : void Atoms::setKbT(const TypesafePtr & p) {
458 59 : MD2double(p,kbT);
459 59 : }
460 :
461 1186 : double Atoms::getKbT()const {
462 1186 : return kbT/units.getEnergy()*MDUnits.getEnergy();
463 : }
464 :
465 :
466 116 : void Atoms::createFullList(const TypesafePtr & n) {
467 116 : if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
468 7 : n.set(int(natoms));
469 7 : fullList.resize(natoms);
470 803 : for(unsigned i=0; i<natoms; i++) fullList[i]=i;
471 : } else {
472 : // We update here the unique list defined at Atoms::unique.
473 : // This is not very clear, and probably should be coded differently.
474 : // Hopefully this fix the longstanding issue with NAMD.
475 : unique.clear();
476 892 : for(unsigned i=0; i<actions.size(); i++) {
477 783 : if(actions[i]->isActive()) {
478 625 : if(!actions[i]->getUnique().empty()) {
479 546 : atomsNeeded=true;
480 : // unique are the local atoms
481 546 : unique.insert(actions[i]->getUnique().begin(),actions[i]->getUnique().end());
482 : }
483 : }
484 : }
485 109 : fullList.resize(0);
486 109 : fullList.reserve(unique.size());
487 5012 : for(const auto & p : unique) fullList.push_back(p.index());
488 109 : n.set(int(fullList.size()));
489 : }
490 116 : }
491 :
492 116 : void Atoms::getFullList(const TypesafePtr & x) {
493 : auto xx=x.template get<const int**>();
494 116 : if(!fullList.empty()) *xx=&fullList[0];
495 6 : else *xx=NULL;
496 116 : }
497 :
498 116 : void Atoms::clearFullList() {
499 116 : fullList.resize(0);
500 116 : }
501 :
502 914 : void Atoms::init() {
503 : // Default: set domain decomposition to NO-decomposition, waiting for
504 : // further instruction
505 914 : if(dd) {
506 337 : setAtomsNlocal(natoms);
507 337 : setAtomsContiguous(0);
508 : }
509 914 : }
510 :
511 337 : void Atoms::setDomainDecomposition(Communicator& comm) {
512 337 : dd.enable(comm);
513 337 : }
514 :
515 14352 : void Atoms::resizeVectors(unsigned n) {
516 14352 : positions.resize(n);
517 14352 : forces.resize(n);
518 14352 : masses.resize(n);
519 14352 : charges.resize(n);
520 14352 : }
521 :
522 7176 : AtomNumber Atoms::addVirtualAtom(ActionWithVirtualAtom*a) {
523 7176 : unsigned n=positions.size();
524 7176 : resizeVectors(n+1);
525 7176 : virtualAtomsActions.push_back(a);
526 7176 : return AtomNumber::index(n);
527 : }
528 :
529 7176 : void Atoms::removeVirtualAtom(ActionWithVirtualAtom*a) {
530 7176 : unsigned n=positions.size();
531 7176 : plumed_massert(a==virtualAtomsActions[virtualAtomsActions.size()-1],"virtual atoms should be destroyed in reverse creation order");
532 7176 : resizeVectors(n-1);
533 : virtualAtomsActions.pop_back();
534 7176 : }
535 :
536 119 : void Atoms::insertGroup(const std::string&name,const std::vector<AtomNumber>&a) {
537 0 : plumed_massert(groups.count(name)==0,"group named "+name+" already exists");
538 119 : groups[name]=a;
539 119 : }
540 :
541 119 : void Atoms::removeGroup(const std::string&name) {
542 0 : plumed_massert(groups.count(name)==1,"cannot remove group named "+name);
543 : groups.erase(name);
544 119 : }
545 :
546 114 : void Atoms::writeBinary(std::ostream&o)const {
547 114 : o.write(reinterpret_cast<const char*>(&positions[0][0]),natoms*3*sizeof(double));
548 114 : o.write(reinterpret_cast<const char*>(&box(0,0)),9*sizeof(double));
549 114 : o.write(reinterpret_cast<const char*>(&energy),sizeof(double));
550 114 : }
551 :
552 114 : void Atoms::readBinary(std::istream&i) {
553 114 : i.read(reinterpret_cast<char*>(&positions[0][0]),natoms*3*sizeof(double));
554 114 : i.read(reinterpret_cast<char*>(&box(0,0)),9*sizeof(double));
555 114 : i.read(reinterpret_cast<char*>(&energy),sizeof(double));
556 114 : pbc.setBox(box);
557 114 : }
558 :
559 591 : double Atoms::getKBoltzmann()const {
560 591 : if(naturalUnits || MDnaturalUnits) return 1.0;
561 585 : else return kBoltzmann/units.getEnergy();
562 : }
563 :
564 0 : double Atoms::getMDKBoltzmann()const {
565 0 : if(naturalUnits || MDnaturalUnits) return 1.0;
566 0 : else return kBoltzmann/MDUnits.getEnergy();
567 : }
568 :
569 0 : void Atoms::getLocalMasses(std::vector<double>& localMasses) {
570 0 : localMasses.resize(gatindex.size());
571 0 : for(unsigned i=0; i<gatindex.size(); i++) localMasses[i] = masses[gatindex[i]];
572 0 : }
573 :
574 450 : void Atoms::getLocalPositions(std::vector<Vector>& localPositions) {
575 450 : localPositions.resize(gatindex.size());
576 450 : mdatoms->getLocalPositions(localPositions);
577 450 : }
578 :
579 450 : void Atoms::getLocalForces(std::vector<Vector>& localForces) {
580 450 : localForces.resize(gatindex.size());
581 16650 : for(unsigned i=0; i<gatindex.size(); i++) localForces[i] = forces[gatindex[i]];
582 450 : }
583 :
584 0 : void Atoms::getLocalMDForces(std::vector<Vector>& localForces) {
585 0 : localForces.resize(gatindex.size());
586 0 : for(unsigned i=0; i<gatindex.size(); i++) {
587 0 : localForces[i] = mdatoms->getMDforces(i);
588 : }
589 0 : }
590 :
591 20 : void Atoms::setExtraCV(const std::string &name,const TypesafePtr & p) {
592 20 : mdatoms->setExtraCV(name,p);
593 20 : }
594 :
595 20 : void Atoms::setExtraCVForce(const std::string &name,const TypesafePtr & p) {
596 20 : mdatoms->setExtraCVForce(name,p);
597 20 : }
598 :
599 60 : double Atoms::getExtraCV(const std::string &name) {
600 60 : return mdatoms->getExtraCV(name);
601 : }
602 :
603 40 : void Atoms::updateExtraCVForce(const std::string &name,double f) {
604 40 : mdatoms->updateExtraCVForce(name,f);
605 40 : }
606 :
607 : }
|
