LCOV - code coverage report
Current view: top level - core - Atoms.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 313 349 89.7 %
Date: 2018-12-19 07:49:13 Functions: 46 53 86.8 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2018 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        1159 : Atoms::Atoms(PlumedMain&plumed):
      43             :   natoms(0),
      44        1159 :   pbc(*new Pbc),
      45             :   md_energy(0.0),
      46             :   energy(0.0),
      47             :   dataCanBeSet(false),
      48             :   collectEnergy(false),
      49             :   energyHasBeenSet(false),
      50             :   positionsHaveBeenSet(0),
      51             :   massesHaveBeenSet(false),
      52             :   chargesHaveBeenSet(false),
      53             :   boxHasBeenSet(false),
      54             :   forcesHaveBeenSet(0),
      55             :   virialHasBeenSet(false),
      56             :   massAndChargeOK(false),
      57             :   shuffledAtoms(0),
      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        2318 :   ddStep(0)
      68             : {
      69        1159 :   mdatoms=MDAtomsBase::create(sizeof(double));
      70        1159 : }
      71             : 
      72        2318 : Atoms::~Atoms() {
      73        1159 :   if(actions.size()>0) {
      74           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";
      75             :   }
      76        1159 :   delete mdatoms;
      77        1159 :   delete &pbc;
      78        1159 : }
      79             : 
      80       21537 : void Atoms::startStep() {
      81       21537 :   collectEnergy=false; energyHasBeenSet=false; positionsHaveBeenSet=0;
      82       21537 :   massesHaveBeenSet=false; chargesHaveBeenSet=false; boxHasBeenSet=false;
      83       21537 :   forcesHaveBeenSet=0; virialHasBeenSet=false; dataCanBeSet=true;
      84       21537 : }
      85             : 
      86       20963 : void Atoms::setBox(void*p) {
      87       20963 :   mdatoms->setBox(p);
      88       20963 :   Tensor b; mdatoms->getBox(b); boxHasBeenSet=true;
      89       20963 : }
      90             : 
      91       20991 : void Atoms::setPositions(void*p) {
      92       20991 :   plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
      93       20991 :   plumed_massert( p || gatindex.size()==0, "NULL position pointer with non-zero local atoms");
      94       20991 :   mdatoms->setp(p); positionsHaveBeenSet=3;
      95       20991 : }
      96             : 
      97       20991 : void Atoms::setMasses(void*p) {
      98       20991 :   plumed_massert( dataCanBeSet,"setMasses must be called after setStep in MD code interface");
      99       20991 :   plumed_massert( p || gatindex.size()==0, "NULL mass pointer with non-zero local atoms");
     100       20991 :   mdatoms->setm(p); massesHaveBeenSet=true;
     101             : 
     102       20991 : }
     103             : 
     104       18773 : void Atoms::setCharges(void*p) {
     105       18773 :   plumed_massert( dataCanBeSet, "setCharges must be called after setStep in MD code interface");
     106       18773 :   plumed_massert( p || gatindex.size()==0, "NULL charges pointer with non-zero local atoms");
     107       18773 :   mdatoms->setc(p); chargesHaveBeenSet=true;
     108       18773 : }
     109             : 
     110       18841 : void Atoms::setVirial(void*p) {
     111       18841 :   plumed_massert( dataCanBeSet,"setVirial must be called after setStep in MD code interface");
     112       18841 :   mdatoms->setVirial(p); virialHasBeenSet=true;
     113             : 
     114       18841 : }
     115             : 
     116        2150 : void Atoms::setEnergy(void*p) {
     117        2150 :   plumed_massert( dataCanBeSet,"setEnergy must be called after setStep in MD code interface");
     118        2150 :   MD2double(p,md_energy);
     119        2150 :   md_energy*=MDUnits.getEnergy()/units.getEnergy();
     120        2150 :   energyHasBeenSet=true;
     121        2150 : }
     122             : 
     123       20991 : void Atoms::setForces(void*p) {
     124       20991 :   plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
     125       20991 :   plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
     126       20991 :   forcesHaveBeenSet=3;
     127       20991 :   mdatoms->setf(p);
     128       20991 : }
     129             : 
     130           0 : void Atoms::setPositions(void*p,int i) {
     131           0 :   plumed_massert( dataCanBeSet,"setPositions must be called after setStep in MD code interface");
     132           0 :   plumed_massert( p || gatindex.size()==0, "NULL positions pointer with non-zero local atoms");
     133           0 :   mdatoms->setp(p,i); positionsHaveBeenSet++;
     134           0 : }
     135             : 
     136           0 : void Atoms::setForces(void*p,int i) {
     137           0 :   plumed_massert( dataCanBeSet,"setForces must be called after setStep in MD code interface");
     138           0 :   plumed_massert( p || gatindex.size()==0, "NULL force pointer with non-zero local atoms");
     139           0 :   mdatoms->setf(p,i); forcesHaveBeenSet++;
     140           0 : }
     141             : 
     142       19518 : void Atoms::share() {
     143       19518 :   std::set<AtomNumber> unique;
     144             : // At first step I scatter all the atoms so as to store their mass and charge
     145             : // Notice that this works with the assumption that charges and masses are
     146             : // not changing during the simulation!
     147       19518 :   if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
     148         322 :     shareAll();
     149       19840 :     return;
     150             :   }
     151       79111 :   for(unsigned i=0; i<actions.size(); i++) if(actions[i]->isActive()) {
     152       54392 :       if(dd && shuffledAtoms>0) {
     153        1132 :         unique.insert(actions[i]->getUnique().begin(),actions[i]->getUnique().end());
     154             :       }
     155       54392 :       if(!actions[i]->getUnique().empty()) atomsNeeded=true;
     156             :     }
     157       19196 :   share(unique);
     158             : }
     159             : 
     160         394 : void Atoms::shareAll() {
     161         394 :   std::set<AtomNumber> unique;
     162         394 :   if(dd && shuffledAtoms>0)
     163         106 :     for(int i=0; i<natoms; i++) unique.insert(AtomNumber::index(i));
     164         394 :   atomsNeeded=true;
     165         394 :   share(unique);
     166         394 : }
     167             : 
     168       19590 : void Atoms::share(const std::set<AtomNumber>& unique) {
     169       19590 :   plumed_assert( positionsHaveBeenSet==3 && massesHaveBeenSet );
     170       19590 :   virial.zero();
     171       19590 :   if(zeroallforces || int(gatindex.size())==natoms) {
     172       19314 :     for(int i=0; i<natoms; i++) forces[i].zero();
     173             :   } else {
     174         276 :     for(unsigned i=0; i<gatindex.size(); i++) forces[gatindex[i]].zero();
     175             :   }
     176       19590 :   for(unsigned i=getNatoms(); i<positions.size(); i++) forces[i].zero(); // virtual atoms
     177       19590 :   forceOnEnergy=0.0;
     178       19590 :   mdatoms->getBox(box);
     179             : 
     180       39180 :   if(!atomsNeeded) return;
     181             : 
     182       18883 :   atomsNeeded=false;
     183             : 
     184       18883 :   if(int(gatindex.size())==natoms && shuffledAtoms==0) {
     185             : // faster version, which retrieves all atoms
     186       18559 :     mdatoms->getPositions(0,natoms,positions);
     187             :   } else {
     188             : // version that picks only atoms that are available on this proc
     189         324 :     mdatoms->getPositions(gatindex,positions);
     190             :   }
     191             : // how many double per atom should be scattered:
     192       18883 :   int ndata=3;
     193       18883 :   if(!massAndChargeOK) {
     194         322 :     ndata=5;
     195         322 :     masses.assign(masses.size(),NAN);
     196         322 :     charges.assign(charges.size(),NAN);
     197         322 :     mdatoms->getCharges(gatindex,charges);
     198         322 :     mdatoms->getMasses(gatindex,masses);
     199             :   }
     200             : 
     201       18883 :   if(dd && shuffledAtoms>0) {
     202         324 :     if(dd.async) {
     203         304 :       for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) dd.mpi_request_positions[i].wait();
     204         304 :       for(unsigned i=0; i<dd.mpi_request_index.size(); i++)     dd.mpi_request_index[i].wait();
     205             :     }
     206         324 :     int count=0;
     207       26246 :     for(std::set<AtomNumber>::const_iterator p=unique.begin(); p!=unique.end(); ++p) {
     208       25922 :       if(dd.g2l[p->index()]>=0) {
     209       11765 :         dd.indexToBeSent[count]=p->index();
     210       11765 :         dd.positionsToBeSent[ndata*count+0]=positions[p->index()][0];
     211       11765 :         dd.positionsToBeSent[ndata*count+1]=positions[p->index()][1];
     212       11765 :         dd.positionsToBeSent[ndata*count+2]=positions[p->index()][2];
     213       11765 :         if(!massAndChargeOK) {
     214        1245 :           dd.positionsToBeSent[ndata*count+3]=masses[p->index()];
     215        1245 :           dd.positionsToBeSent[ndata*count+4]=charges[p->index()];
     216             :         }
     217       11765 :         count++;
     218             :       }
     219             :     }
     220         324 :     if(dd.async) {
     221         304 :       asyncSent=true;
     222         304 :       dd.mpi_request_positions.resize(dd.Get_size());
     223         304 :       dd.mpi_request_index.resize(dd.Get_size());
     224         992 :       for(int i=0; i<dd.Get_size(); i++) {
     225         688 :         dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
     226         688 :         dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
     227             :       }
     228             :     } else {
     229          20 :       const int n=(dd.Get_size());
     230          20 :       vector<int> counts(n);
     231          40 :       vector<int> displ(n);
     232          40 :       vector<int> counts5(n);
     233          40 :       vector<int> displ5(n);
     234          20 :       dd.Allgather(count,counts);
     235          20 :       displ[0]=0;
     236          20 :       for(int i=1; i<n; ++i) displ[i]=displ[i-1]+counts[i-1];
     237          20 :       for(int i=0; i<n; ++i) counts5[i]=counts[i]*ndata;
     238          20 :       for(int i=0; i<n; ++i) displ5[i]=displ[i]*ndata;
     239          20 :       dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
     240          20 :       dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
     241          20 :       int tot=displ[n-1]+counts[n-1];
     242        1620 :       for(int i=0; i<tot; i++) {
     243        1600 :         positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
     244        1600 :         positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
     245        1600 :         positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
     246        1600 :         if(!massAndChargeOK) {
     247         432 :           masses[dd.indexToBeReceived[i]]      =dd.positionsToBeReceived[ndata*i+3];
     248         432 :           charges[dd.indexToBeReceived[i]]     =dd.positionsToBeReceived[ndata*i+4];
     249             :         }
     250          20 :       }
     251             :     }
     252             :   }
     253             : }
     254             : 
     255       19590 : void Atoms::wait() {
     256       19590 :   dataCanBeSet=false; // Everything should be set by this stage
     257             : // How many double per atom should be scattered
     258       19590 :   int ndata=3;
     259       19590 :   if(!massAndChargeOK)ndata=5;
     260             : 
     261       19590 :   if(dd) {
     262        8274 :     dd.Bcast(box,0);
     263             :   }
     264       19590 :   pbc.setBox(box);
     265             : 
     266       19590 :   if(collectEnergy) energy=md_energy;
     267             : 
     268       19590 :   if(dd && shuffledAtoms>0) {
     269             : // receive toBeReceived
     270         324 :     if(asyncSent) {
     271             :       Communicator::Status status;
     272         304 :       int count=0;
     273         992 :       for(int i=0; i<dd.Get_size(); i++) {
     274         688 :         dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
     275         688 :         int c=status.Get_count<int>();
     276         688 :         dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
     277         688 :         count+=c;
     278             :       }
     279       24626 :       for(int i=0; i<count; i++) {
     280       24322 :         positions[dd.indexToBeReceived[i]][0]=dd.positionsToBeReceived[ndata*i+0];
     281       24322 :         positions[dd.indexToBeReceived[i]][1]=dd.positionsToBeReceived[ndata*i+1];
     282       24322 :         positions[dd.indexToBeReceived[i]][2]=dd.positionsToBeReceived[ndata*i+2];
     283       24322 :         if(!massAndChargeOK) {
     284        2706 :           masses[dd.indexToBeReceived[i]]      =dd.positionsToBeReceived[ndata*i+3];
     285        2706 :           charges[dd.indexToBeReceived[i]]     =dd.positionsToBeReceived[ndata*i+4];
     286             :         }
     287             :       }
     288         304 :       asyncSent=false;
     289             :     }
     290         324 :     if(collectEnergy) dd.Sum(energy);
     291             :   }
     292             : // I take note that masses and charges have been set once for all
     293             : // at the beginning of the simulation.
     294       19590 :   if(shareMassAndChargeOnlyAtFirstStep) massAndChargeOK=true;
     295       19590 : }
     296             : 
     297       19518 : void Atoms::updateForces() {
     298       19518 :   plumed_assert( forcesHaveBeenSet==3 );
     299       19518 :   if(forceOnEnergy*forceOnEnergy>epsilon) {
     300          50 :     double alpha=1.0-forceOnEnergy;
     301          50 :     mdatoms->rescaleForces(gatindex,alpha);
     302             :   }
     303       19518 :   mdatoms->updateForces(gatindex,forces);
     304       19518 :   if( !plumed.novirial && dd.Get_rank()==0 ) {
     305       14811 :     plumed_assert( virialHasBeenSet );
     306       14811 :     mdatoms->updateVirial(virial);
     307             :   }
     308       19518 : }
     309             : 
     310         331 : void Atoms::setNatoms(int n) {
     311         331 :   natoms=n;
     312         331 :   positions.resize(n);
     313         331 :   forces.resize(n);
     314         331 :   masses.resize(n);
     315         331 :   charges.resize(n);
     316         331 :   gatindex.resize(n);
     317         331 :   for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=i;
     318         331 : }
     319             : 
     320             : 
     321        1356 : void Atoms::add(const ActionAtomistic*a) {
     322        1356 :   actions.push_back(a);
     323        1356 : }
     324             : 
     325        1356 : void Atoms::remove(const ActionAtomistic*a) {
     326        1356 :   vector<const ActionAtomistic*>::iterator f=find(actions.begin(),actions.end(),a);
     327        1356 :   plumed_massert(f!=actions.end(),"cannot remove an action registered to atoms");
     328        1356 :   actions.erase(f);
     329        1356 : }
     330             : 
     331             : 
     332         115 : void Atoms::DomainDecomposition::enable(Communicator& c) {
     333         115 :   on=true;
     334         115 :   Set_comm(c.Get_comm());
     335         115 :   async=Get_size()<10;
     336         115 :   if(std::getenv("PLUMED_ASYNC_SHARE")) {
     337           4 :     std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
     338           4 :     if(s=="yes") async=true;
     339           4 :     else if(s=="no") async=false;
     340           0 :     else plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
     341             :   }
     342         115 : }
     343             : 
     344         265 : void Atoms::setAtomsNlocal(int n) {
     345         265 :   gatindex.resize(n);
     346         265 :   if(dd) {
     347         265 :     dd.g2l.resize(natoms,-1);
     348             : // Since these vectors are sent with MPI by using e.g.
     349             : // &dd.positionsToBeSent[0]
     350             : // we make sure they are non-zero-sized so as to
     351             : // avoid errors when doing boundary check
     352         265 :     if(n==0) n++;
     353         265 :     dd.positionsToBeSent.resize(n*5,0.0);
     354         265 :     dd.positionsToBeReceived.resize(natoms*5,0.0);
     355         265 :     dd.indexToBeSent.resize(n,0);
     356         265 :     dd.indexToBeReceived.resize(natoms,0);
     357             :   };
     358         265 : }
     359             : 
     360         134 : void Atoms::setAtomsGatindex(int*g,bool fortran) {
     361         134 :   plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
     362         134 :   ddStep=plumed.getStep();
     363         134 :   if(fortran) {
     364           0 :     for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=g[i]-1;
     365             :   } else {
     366         134 :     for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=g[i];
     367             :   }
     368         134 :   for(unsigned i=0; i<dd.g2l.size(); i++) dd.g2l[i]=-1;
     369         134 :   if( gatindex.size()==natoms ) {
     370           0 :     shuffledAtoms=0;
     371           0 :     for(unsigned i=0; i<gatindex.size(); i++) {
     372           0 :       if( gatindex[i]!=i ) { shuffledAtoms=1; break; }
     373             :     }
     374             :   } else {
     375         134 :     shuffledAtoms=1;
     376             :   }
     377         134 :   if(dd) {
     378         134 :     dd.Sum(shuffledAtoms);
     379         134 :     for(unsigned i=0; i<gatindex.size(); i++) dd.g2l[gatindex[i]]=i;
     380             :   }
     381         134 : }
     382             : 
     383         131 : void Atoms::setAtomsContiguous(int start) {
     384         131 :   ddStep=plumed.getStep();
     385         131 :   for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=start+i;
     386         131 :   for(unsigned i=0; i<dd.g2l.size(); i++) dd.g2l[i]=-1;
     387         131 :   if(dd) for(unsigned i=0; i<gatindex.size(); i++) dd.g2l[gatindex[i]]=i;
     388         131 :   if(gatindex.size()<natoms) shuffledAtoms=1;
     389         131 : }
     390             : 
     391         292 : void Atoms::setRealPrecision(int p) {
     392         292 :   delete mdatoms;
     393         292 :   mdatoms=MDAtomsBase::create(p);
     394         292 : }
     395             : 
     396         331 : int Atoms::getRealPrecision()const {
     397         331 :   return mdatoms->getRealPrecision();
     398             : }
     399             : 
     400        3310 : void Atoms::MD2double(const void*m,double&d)const {
     401        3310 :   plumed_assert(mdatoms); mdatoms->MD2double(m,d);
     402        3310 : }
     403         191 : void Atoms::double2MD(const double&d,void*m)const {
     404         191 :   plumed_assert(mdatoms); mdatoms->double2MD(d,m);
     405         191 : }
     406             : 
     407         331 : void Atoms::updateUnits() {
     408         331 :   mdatoms->setUnits(units,MDUnits);
     409         331 : }
     410             : 
     411         292 : void Atoms::setTimeStep(void*p) {
     412         292 :   MD2double(p,timestep);
     413         292 : }
     414             : 
     415      676533 : double Atoms::getTimeStep()const {
     416      676533 :   return timestep/units.getTime()*MDUnits.getTime();
     417             : }
     418             : 
     419           4 : void Atoms::setKbT(void*p) {
     420           4 :   MD2double(p,kbT);
     421           4 : }
     422             : 
     423         417 : double Atoms::getKbT()const {
     424         417 :   return kbT/units.getEnergy()*MDUnits.getEnergy();
     425             : }
     426             : 
     427             : 
     428           5 : void Atoms::createFullList(int*n) {
     429           5 :   vector<AtomNumber> fullListTmp;
     430          15 :   for(unsigned i=0; i<actions.size(); i++) if(actions[i]->isActive())
     431           5 :       fullListTmp.insert(fullListTmp.end(),actions[i]->getUnique().begin(),actions[i]->getUnique().end());
     432           5 :   std::sort(fullListTmp.begin(),fullListTmp.end());
     433           5 :   int nn=std::unique(fullListTmp.begin(),fullListTmp.end())-fullListTmp.begin();
     434           5 :   fullList.resize(nn);
     435           5 :   for(int i=0; i<nn; ++i) fullList[i]=fullListTmp[i].index();
     436           5 :   *n=nn;
     437           5 : }
     438             : 
     439           5 : void Atoms::getFullList(int**x) {
     440           5 :   if(!fullList.empty()) *x=&fullList[0];
     441           1 :   else *x=NULL;
     442           5 : }
     443             : 
     444           5 : void Atoms::clearFullList() {
     445           5 :   fullList.resize(0);
     446           5 : }
     447             : 
     448         331 : void Atoms::init() {
     449             : // Default: set domain decomposition to NO-decomposition, waiting for
     450             : // further instruction
     451         331 :   if(dd) {
     452         115 :     setAtomsNlocal(natoms);
     453         115 :     setAtomsContiguous(0);
     454             :   }
     455         331 : }
     456             : 
     457         115 : void Atoms::setDomainDecomposition(Communicator& comm) {
     458         115 :   dd.enable(comm);
     459         115 : }
     460             : 
     461         234 : void Atoms::resizeVectors(unsigned n) {
     462         234 :   positions.resize(n);
     463         234 :   forces.resize(n);
     464         234 :   masses.resize(n);
     465         234 :   charges.resize(n);
     466         234 : }
     467             : 
     468         117 : AtomNumber Atoms::addVirtualAtom(ActionWithVirtualAtom*a) {
     469         117 :   unsigned n=positions.size();
     470         117 :   resizeVectors(n+1);
     471         117 :   virtualAtomsActions.push_back(a);
     472         117 :   return AtomNumber::index(n);
     473             : }
     474             : 
     475         117 : void Atoms::removeVirtualAtom(ActionWithVirtualAtom*a) {
     476         117 :   unsigned n=positions.size();
     477         117 :   plumed_massert(a==virtualAtomsActions[virtualAtomsActions.size()-1],"virtual atoms should be destroyed in reverse creation order");
     478         117 :   resizeVectors(n-1);
     479         117 :   virtualAtomsActions.pop_back();
     480         117 : }
     481             : 
     482          63 : void Atoms::insertGroup(const std::string&name,const std::vector<AtomNumber>&a) {
     483          63 :   plumed_massert(groups.count(name)==0,"group named "+name+" already exists");
     484          63 :   groups[name]=a;
     485          63 : }
     486             : 
     487          63 : void Atoms::removeGroup(const std::string&name) {
     488          63 :   plumed_massert(groups.count(name)==1,"cannot remove group named "+name);
     489          63 :   groups.erase(name);
     490          63 : }
     491             : 
     492          72 : void Atoms::writeBinary(std::ostream&o)const {
     493          72 :   o.write(reinterpret_cast<const char*>(&positions[0][0]),natoms*3*sizeof(double));
     494          72 :   o.write(reinterpret_cast<const char*>(&box(0,0)),9*sizeof(double));
     495          72 :   o.write(reinterpret_cast<const char*>(&energy),sizeof(double));
     496          72 : }
     497             : 
     498          72 : void Atoms::readBinary(std::istream&i) {
     499          72 :   i.read(reinterpret_cast<char*>(&positions[0][0]),natoms*3*sizeof(double));
     500          72 :   i.read(reinterpret_cast<char*>(&box(0,0)),9*sizeof(double));
     501          72 :   i.read(reinterpret_cast<char*>(&energy),sizeof(double));
     502          72 :   pbc.setBox(box);
     503          72 : }
     504             : 
     505          53 : double Atoms::getKBoltzmann()const {
     506          53 :   if(naturalUnits || MDnaturalUnits) return 1.0;
     507          46 :   else return kBoltzmann/units.getEnergy();
     508             : }
     509             : 
     510           0 : double Atoms::getMDKBoltzmann()const {
     511           0 :   if(naturalUnits || MDnaturalUnits) return 1.0;
     512           0 :   else return kBoltzmann/MDUnits.getEnergy();
     513             : }
     514             : 
     515           0 : void Atoms::getLocalMasses(std::vector<double>& localMasses) {
     516           0 :   localMasses.resize(gatindex.size());
     517           0 :   for(unsigned i=0; i<gatindex.size(); i++) localMasses[i] = masses[gatindex[i]];
     518           0 : }
     519             : 
     520           0 : void Atoms::getLocalPositions(std::vector<Vector>& localPositions) {
     521           0 :   localPositions.resize(gatindex.size());
     522           0 :   mdatoms->getLocalPositions(localPositions);
     523           0 : }
     524             : 
     525           0 : void Atoms::getLocalForces(std::vector<Vector>& localForces) {
     526           0 :   localForces.resize(gatindex.size());
     527           0 :   for(unsigned i=0; i<gatindex.size(); i++) localForces[i] = forces[gatindex[i]];
     528           0 : }
     529             : 
     530           0 : void Atoms::getLocalMDForces(std::vector<Vector>& localForces) {
     531           0 :   localForces.resize(gatindex.size());
     532           0 :   for(unsigned i=0; i<gatindex.size(); i++) {
     533           0 :     localForces[i] = mdatoms->getMDforces(i);
     534             :   }
     535           0 : }
     536             : 
     537        2523 : }

Generated by: LCOV version 1.13