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             : #ifndef __PLUMED_core_Atoms_h
      23             : #define __PLUMED_core_Atoms_h
      24             : 
      25             : #include "tools/TypesafePtr.h"
      26             : #include "tools/Communicator.h"
      27             : #include "tools/Tensor.h"
      28             : #include "tools/Units.h"
      29             : #include "tools/Exception.h"
      30             : #include "tools/AtomNumber.h"
      31             : #include "tools/ForwardDecl.h"
      32             : #include <vector>
      33             : #include <map>
      34             : #include <string>
      35             : #include <memory>
      36             : 
      37             : namespace PLMD {
      38             : 
      39             : class MDAtomsBase;
      40             : class PlumedMain;
      41             : class ActionAtomistic;
      42             : class ActionWithVirtualAtom;
      43             : class Pbc;
      44             : 
      45             : /// Class containing atom related quantities from the MD code.
      46             : /// IT IS STILL UNDOCUMENTED. IT PROBABLY NEEDS A STRONG CLEANUP
      47             : class Atoms {
      48             :   friend class ActionAtomistic;
      49             :   friend class ActionWithVirtualAtom;
      50             :   int natoms;
      51             :   bool unique_serial=false; // use unique in serial mode
      52             :   std::vector<AtomNumber> unique;
      53             :   std::vector<unsigned> uniq_index;
      54             : /// Map global indexes to local indexes
      55             : /// E.g. g2l[i] is the position of atom i in the array passed from the MD engine.
      56             : /// Called "global to local" since originally it was used to map global indexes to local
      57             : /// ones used in domain decomposition. However, it is now also used for the NAMD-like
      58             : /// interface, where only a small number of atoms is passed to plumed.
      59             :   std::vector<int> g2l;
      60             :   std::vector<Vector> positions;
      61             :   std::vector<Vector> forces;
      62             :   std::vector<double> masses;
      63             :   std::vector<double> charges;
      64             :   std::vector<ActionWithVirtualAtom*> virtualAtomsActions;
      65             :   Tensor box;
      66             :   ForwardDecl<Pbc> pbc_fwd;
      67             :   Pbc&   pbc=*pbc_fwd;
      68             :   Tensor virial;
      69             : // this is the energy set by each processor:
      70             :   double md_energy;
      71             : // this is the summed energy:
      72             :   double energy;
      73             : 
      74             :   bool   dataCanBeSet;
      75             :   bool   collectEnergy;
      76             :   bool   energyHasBeenSet;
      77             :   unsigned positionsHaveBeenSet;
      78             :   bool massesHaveBeenSet;
      79             :   bool chargesHaveBeenSet;
      80             :   bool boxHasBeenSet;
      81             :   unsigned forcesHaveBeenSet;
      82             :   bool virialHasBeenSet;
      83             :   bool massAndChargeOK;
      84             :   unsigned shuffledAtoms;
      85             : 
      86             :   std::map<std::string,std::vector<AtomNumber> > groups;
      87             : 
      88             :   void resizeVectors(unsigned);
      89             : 
      90             :   std::vector<int> fullList;
      91             : 
      92             :   std::unique_ptr<MDAtomsBase> mdatoms;
      93             : 
      94             :   PlumedMain & plumed;
      95             : 
      96             :   Units MDUnits;
      97             :   Units units;
      98             : 
      99             :   bool naturalUnits;
     100             :   bool MDnaturalUnits;
     101             : 
     102             :   double timestep;
     103             :   double forceOnEnergy;
     104             : 
     105             : /// if set to true, all the forces in the global array are zeroes
     106             : /// at every step. It should not be necessary in general, but it is
     107             : /// for actions accessing to modifyGlobalForce() (e.g. FIT_TO_TEMPLATE).
     108             :   bool zeroallforces;
     109             : 
     110             :   double kbT;
     111             : 
     112             :   std::vector<ActionAtomistic*> actions;
     113             :   std::vector<int>    gatindex;
     114             : 
     115             :   bool asyncSent;
     116             :   bool atomsNeeded;
     117             : 
     118             :   class DomainDecomposition:
     119             :     public Communicator {
     120             :   public:
     121             :     bool on;
     122             :     bool async;
     123             : 
     124             :     std::vector<Communicator::Request> mpi_request_positions;
     125             :     std::vector<Communicator::Request> mpi_request_index;
     126             : 
     127             :     std::vector<double> positionsToBeSent;
     128             :     std::vector<double> positionsToBeReceived;
     129             :     std::vector<int>    indexToBeSent;
     130             :     std::vector<int>    indexToBeReceived;
     131             :     operator bool() const {
     132      193259 :       return on;
     133             :     }
     134      805723 :     DomainDecomposition():
     135      805723 :       on(false), async(false)
     136      805723 :     {}
     137             :     void enable(Communicator& c);
     138             :   };
     139             : 
     140             :   DomainDecomposition dd;
     141             :   long long int ddStep;  //last step in which dd happened
     142             : 
     143             :   void share(const std::vector<AtomNumber>&);
     144             : 
     145             : public:
     146             : 
     147             :   explicit Atoms(PlumedMain&plumed);
     148             :   ~Atoms();
     149             : 
     150             :   void init();
     151             : 
     152             :   void share();
     153             :   void shareAll();
     154             :   void wait();
     155             :   void updateForces();
     156             : 
     157             :   void setRealPrecision(int);
     158             :   int  getRealPrecision()const;
     159             : 
     160             :   void setTimeStep(const TypesafePtr &);
     161             :   double getTimeStep()const;
     162             : 
     163             :   void setKbT(const TypesafePtr &);
     164             :   double getKbT()const;
     165             : 
     166             :   void setNatoms(int);
     167             :   int getNatoms()const;
     168             :   int getNVirtualAtoms()const;
     169             : 
     170             :   const long long int& getDdStep()const;
     171             :   const std::vector<int>& getGatindex()const;
     172             :   const Pbc& getPbc()const;
     173             :   void getLocalMasses(std::vector<double>&);
     174             :   void getLocalPositions(std::vector<Vector>&);
     175             :   void getLocalForces(std::vector<Vector>&);
     176             :   void getLocalMDForces(std::vector<Vector>&);
     177             :   const Tensor& getVirial()const;
     178             : 
     179             :   void setCollectEnergy(bool b) {
     180        3989 :     collectEnergy=b;
     181             :   }
     182             : 
     183             :   void setDomainDecomposition(Communicator&);
     184             :   void setAtomsGatindex(const TypesafePtr &,bool);
     185             :   void setAtomsContiguous(int);
     186             :   void setAtomsNlocal(int);
     187             : 
     188             :   void startStep();
     189             :   void setEnergy(const TypesafePtr &);
     190             :   void setBox(const TypesafePtr &);
     191             :   void setVirial(const TypesafePtr &);
     192             :   void setPositions(const TypesafePtr &);
     193             :   void setPositions(const TypesafePtr &,int);
     194             :   void setForces(const TypesafePtr &);
     195             :   void setForces(const TypesafePtr &,int);
     196             :   void setMasses(const TypesafePtr &);
     197             :   void setCharges(const TypesafePtr &);
     198             :   bool chargesWereSet() const ;
     199             :   bool boxWasSet() const ;
     200             : 
     201             :   void MD2double(const TypesafePtr & m,double&d)const;
     202             :   void double2MD(const double&d,const TypesafePtr & m)const;
     203             : 
     204             :   void createFullList(const TypesafePtr &);
     205             :   void getFullList(const TypesafePtr &);
     206             :   void clearFullList();
     207             : 
     208             :   void add(ActionAtomistic*);
     209             :   void remove(ActionAtomistic*);
     210             : 
     211        3989 :   double getEnergy()const {
     212        3989 :     plumed_assert(collectEnergy && energyHasBeenSet);
     213        3989 :     return energy;
     214             :   }
     215             : 
     216             :   bool isEnergyNeeded()const {
     217          96 :     return collectEnergy;
     218             :   }
     219             : 
     220             :   void setMDEnergyUnits(double d) {
     221          45 :     MDUnits.setEnergy(d);
     222          45 :   }
     223             :   void setMDLengthUnits(double d) {
     224         837 :     MDUnits.setLength(d);
     225         837 :   }
     226             :   void setMDTimeUnits(double d) {
     227           6 :     MDUnits.setTime(d);
     228           6 :   }
     229             :   void setMDChargeUnits(double d) {
     230         837 :     MDUnits.setCharge(d);
     231         837 :   }
     232             :   void setMDMassUnits(double d) {
     233         837 :     MDUnits.setMass(d);
     234         837 :   }
     235             :   const Units& getMDUnits() {
     236             :     return MDUnits;
     237             :   }
     238             :   void setUnits(const Units&u) {
     239          20 :     units=u;
     240          20 :   }
     241             :   const Units& getUnits() {
     242             :     return units;
     243             :   }
     244             :   void updateUnits();
     245             : 
     246             :   AtomNumber addVirtualAtom(ActionWithVirtualAtom*);
     247             :   void removeVirtualAtom(ActionWithVirtualAtom*);
     248             :   ActionWithVirtualAtom* getVirtualAtomsAction(AtomNumber)const;
     249             :   bool isVirtualAtom(AtomNumber)const;
     250             :   void insertGroup(const std::string&name,const std::vector<AtomNumber>&a);
     251             :   void removeGroup(const std::string&name);
     252             :   void writeBinary(std::ostream&)const;
     253             :   void readBinary(std::istream&);
     254             :   double getKBoltzmann()const;
     255             :   double getMDKBoltzmann()const;
     256             :   bool usingNaturalUnits()const;
     257             :   void setNaturalUnits(bool n) {
     258          20 :     naturalUnits=n;
     259             :   }
     260             :   void setMDNaturalUnits(bool n) {
     261           0 :     MDnaturalUnits=n;
     262           0 :   }
     263             : 
     264             :   void setExtraCV(const std::string &name,const TypesafePtr & p);
     265             :   void setExtraCVForce(const std::string &name,const TypesafePtr & p);
     266             :   double getExtraCV(const std::string &name);
     267             :   void updateExtraCVForce(const std::string &name,double f);
     268             :   void setExtraCVNeeded(const std::string &name,bool needed=true);
     269             :   bool isExtraCVNeeded(const std::string &name) const;
     270             :   void resetExtraCVNeeded();
     271             : };
     272             : 
     273             : inline
     274             : int Atoms::getNatoms()const {
     275      808204 :   return natoms;
     276             : }
     277             : 
     278             : inline
     279             : int Atoms::getNVirtualAtoms()const {
     280          10 :   return virtualAtomsActions.size();
     281             : }
     282             : 
     283             : inline
     284             : const long long int& Atoms::getDdStep()const {
     285             :   return ddStep;
     286             : }
     287             : 
     288             : inline
     289             : const std::vector<int>& Atoms::getGatindex()const {
     290           9 :   return gatindex;
     291             : }
     292             : 
     293             : inline
     294             : const Pbc& Atoms::getPbc()const {
     295         900 :   return pbc;
     296             : }
     297             : 
     298             : inline
     299             : bool Atoms::isVirtualAtom(AtomNumber i)const {
     300     1406137 :   return i.index()>=(unsigned) getNatoms();
     301             : }
     302             : 
     303             : inline
     304             : ActionWithVirtualAtom* Atoms::getVirtualAtomsAction(AtomNumber i)const {
     305        7234 :   return virtualAtomsActions[i.index()-getNatoms()];
     306             : }
     307             : 
     308             : inline
     309             : bool Atoms::usingNaturalUnits() const {
     310        2514 :   return naturalUnits || MDnaturalUnits;
     311             : }
     312             : 
     313             : inline
     314             : bool Atoms::chargesWereSet() const {
     315      182561 :   return chargesHaveBeenSet;
     316             : }
     317             : 
     318             : inline
     319             : bool Atoms::boxWasSet() const {
     320             :   return boxHasBeenSet;
     321             : }
     322             : 
     323             : inline
     324             : const Tensor& Atoms::getVirial()const {
     325         450 :   return virial;
     326             : }
     327             : 
     328             : 
     329             : }
     330             : #endif