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_ActionAtomistic_h
      23             : #define __PLUMED_core_ActionAtomistic_h
      24             : 
      25             : #include "Action.h"
      26             : #include "tools/Tensor.h"
      27             : #include "Atoms.h"
      28             : #include "tools/Pbc.h"
      29             : #include "tools/ForwardDecl.h"
      30             : #include <vector>
      31             : #include <map>
      32             : 
      33             : namespace PLMD {
      34             : 
      35             : class Pbc;
      36             : class PDB;
      37             : 
      38             : /// \ingroup MULTIINHERIT
      39             : /// Action used to create objects that access the positions of the atoms from the MD code
      40             : class ActionAtomistic :
      41             :   virtual public Action {
      42             : 
      43             :   std::vector<AtomNumber> indexes;         // the set of needed atoms
      44             : /// unique should be an ordered set since we later create a vector containing the corresponding indexes
      45             :   std::vector<AtomNumber>  unique;
      46             : /// unique_local should be an ordered set since we later create a vector containing the corresponding indexes
      47             :   std::vector<AtomNumber>  unique_local;
      48             :   std::vector<Vector>   positions;       // positions of the needed atoms
      49             :   double                energy;
      50             :   ForwardDecl<Pbc>      pbc_fwd;
      51             :   Pbc&                  pbc=*pbc_fwd;
      52             :   Tensor                virial;
      53             :   std::vector<double>   masses;
      54             :   bool                  chargesWereSet;
      55             :   std::vector<double>   charges;
      56             : 
      57             :   std::vector<Vector>   forces;          // forces on the needed atoms
      58             :   double                forceOnEnergy;
      59             : 
      60             :   double                forceOnExtraCV;
      61             : 
      62             :   std::string           extraCV;
      63             : 
      64             :   bool                  lockRequestAtoms; // forbid changes to request atoms
      65             : 
      66             :   bool                  donotretrieve;
      67             :   bool                  donotforce;
      68             : 
      69             : protected:
      70             :   Atoms&                atoms;
      71             : 
      72             :   void setExtraCV(const std::string &name);
      73             : 
      74             : public:
      75             : /// Request an array of atoms.
      76             : /// This method is used to ask for a list of atoms. Atoms
      77             : /// should be asked for by number. If this routine is called
      78             : /// during the simulation, atoms will be available at the next step
      79             : /// MAYBE WE HAVE TO FIND SOMETHING MORE CLEAR FOR DYNAMIC
      80             : /// LISTS OF ATOMS
      81             :   void requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep=true);
      82             : /// Get position of i-th atom (access by relative index)
      83             :   const Vector & getPosition(int)const;
      84             : /// Get position of i-th atom (access by absolute AtomNumber).
      85             : /// With direct access to the global atom array.
      86             : /// \warning Should be only used by actions that need to read the shared position array.
      87             : ///          This array is insensitive to local changes such as makeWhole(), numerical derivatives, etc.
      88             :   const Vector & getGlobalPosition(AtomNumber)const;
      89             : /// Get modifiable position of i-th atom (access by absolute AtomNumber).
      90             : /// \warning Should be only used by actions that need to modify the shared position array.
      91             : ///          This array is insensitive to local changes such as makeWhole(), numerical derivatives, etc.
      92             :   Vector & modifyGlobalPosition(AtomNumber);
      93             : /// Get total number of atoms, including virtual ones.
      94             : /// Can be used to make a loop on modifyGlobalPosition or getGlobalPosition.
      95             :   unsigned getTotAtoms()const;
      96             : /// Get modifiable force of i-th atom (access by absolute AtomNumber).
      97             : /// \warning  Should be used by action that need to modify the stored atomic forces.
      98             : ///           This should be used with great care since the plumed core does
      99             : ///           not usually keep all these forces up to date. In particular,
     100             : ///           if an action require this, one should during constructor
     101             : ///           call allowToAccessGlobalForces().
     102             : ///           Notice that for efficiency reason plumed does not check if this is done!
     103             :   Vector & modifyGlobalForce(AtomNumber);
     104             : /// Get modifiable virial
     105             : /// Should be used by action that need to modify the stored virial
     106             :   Tensor & modifyGlobalVirial();
     107             : /// Get modifiable PBC
     108             : /// Should be used by action that need to modify the stored box
     109             :   Pbc & modifyGlobalPbc();
     110             : /// Get box shape
     111             :   const Tensor & getBox()const;
     112             : /// Get the array of all positions
     113             :   const std::vector<Vector> & getPositions()const;
     114             : /// Get energy
     115             :   const double & getEnergy()const;
     116             : /// Get mass of i-th atom
     117             :   double getMass(int i)const;
     118             : /// Get charge of i-th atom
     119             :   double getCharge(int i)const;
     120             : /// Get a reference to forces array
     121             :   std::vector<Vector> & modifyForces();
     122             : /// Get a reference to virial array
     123             :   Tensor & modifyVirial();
     124             : /// Get a reference to force on energy
     125             :   double & modifyForceOnEnergy();
     126             : /// Get a reference to force on extraCV
     127             :   double & modifyForceOnExtraCV();
     128             : /// Get number of available atoms
     129             :   unsigned getNumberOfAtoms()const {
     130  4194753890 :     return indexes.size();
     131             :   }
     132             : /// Compute the pbc distance between two positions
     133             :   Vector pbcDistance(const Vector&,const Vector&)const;
     134             : /// Applies  PBCs to a seriens of positions or distances
     135             :   void pbcApply(std::vector<Vector>& dlist, unsigned max_index=0) const;
     136             : /// Get the vector of absolute indexes
     137             :   virtual const std::vector<AtomNumber> & getAbsoluteIndexes()const;
     138             : /// Get the absolute index of an atom
     139             :   AtomNumber getAbsoluteIndex(int i)const;
     140             : /// Parse a list of atoms without a numbered keyword
     141             :   void parseAtomList(const std::string&key,std::vector<AtomNumber> &t);
     142             : /// Parse an list of atom with a numbred keyword
     143             :   void parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t);
     144             : /// Convert a set of read in strings into an atom list (this is used in parseAtomList)
     145             :   void interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t);
     146             : /// Change the box shape
     147             :   void changeBox( const Tensor& newbox );
     148             : /// Get reference to Pbc
     149             :   const Pbc & getPbc() const;
     150             : /// Add the forces to the atoms
     151             :   void setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned ind=0 );
     152             : /// Skip atom retrieval - use with care.
     153             : /// If this function is called during initialization, then atoms are
     154             : /// not going to be retrieved. Can be used for optimization. Notice that
     155             : /// calling getPosition(int) in an Action where DoNotRetrieve() was called might
     156             : /// lead to undefined behavior.
     157             :   void doNotRetrieve() {
     158          54 :     donotretrieve=true;
     159             :   }
     160             : /// Skip atom forces - use with care.
     161             : /// If this function is called during initialization, then forces are
     162             : /// not going to be propagated. Can be used for optimization.
     163             :   void doNotForce() {
     164          54 :     donotforce=true;
     165             :   }
     166             : /// Make atoms whole, assuming they are in the proper order
     167             :   void makeWhole();
     168             : /// Allow calls to modifyGlobalForce()
     169             :   void allowToAccessGlobalForces() {
     170           9 :     atoms.zeroallforces=true;
     171             :   }
     172             : /// updates local unique atoms
     173             :   void updateUniqueLocal();
     174             : public:
     175             : 
     176             : // virtual functions:
     177             : 
     178             :   explicit ActionAtomistic(const ActionOptions&ao);
     179             :   ~ActionAtomistic();
     180             : 
     181             :   static void registerKeywords( Keywords& keys );
     182             : 
     183             :   void clearOutputForces();
     184             : 
     185             : /// N.B. only pass an ActionWithValue to this routine if you know exactly what you
     186             : /// are doing.  The default will be correct for the vast majority of cases
     187             :   void   calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
     188             : /// Numerical derivative routine to use when using Actions that inherit from BOTH
     189             : /// ActionWithArguments and ActionAtomistic
     190             :   void calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum );
     191             : 
     192             :   virtual void retrieveAtoms();
     193             :   void applyForces();
     194             :   void lockRequests() override;
     195             :   void unlockRequests() override;
     196             :   const std::vector<AtomNumber> & getUnique()const;
     197             :   const std::vector<AtomNumber> & getUniqueLocal()const;
     198             : /// Read in an input file containing atom positions and calculate the action for the atomic
     199             : /// configuration therin
     200             :   void readAtomsFromPDB( const PDB& pdb ) override;
     201             : };
     202             : 
     203             : inline
     204             : const Vector & ActionAtomistic::getPosition(int i)const {
     205   558541035 :   return positions[i];
     206             : }
     207             : 
     208             : inline
     209             : const Vector & ActionAtomistic::getGlobalPosition(AtomNumber i)const {
     210      559990 :   return atoms.positions[i.index()];
     211             : }
     212             : 
     213             : inline
     214             : Vector & ActionAtomistic::modifyGlobalPosition(AtomNumber i) {
     215      620268 :   return atoms.positions[i.index()];
     216             : }
     217             : 
     218             : inline
     219             : Vector & ActionAtomistic::modifyGlobalForce(AtomNumber i) {
     220       16318 :   return atoms.forces[i.index()];
     221             : }
     222             : 
     223             : inline
     224             : Tensor & ActionAtomistic::modifyGlobalVirial() {
     225         125 :   return atoms.virial;
     226             : }
     227             : 
     228             : inline
     229             : double ActionAtomistic::getMass(int i)const {
     230      499772 :   return masses[i];
     231             : }
     232             : 
     233             : inline
     234       63884 : double ActionAtomistic::getCharge(int i) const {
     235       63884 :   if( !chargesWereSet ) {
     236           0 :     error("charges were not passed to plumed");
     237             :   }
     238       63884 :   return charges[i];
     239             : }
     240             : 
     241             : inline
     242         160 : const std::vector<AtomNumber> & ActionAtomistic::getAbsoluteIndexes()const {
     243         161 :   return indexes;
     244             : }
     245             : 
     246             : inline
     247             : AtomNumber ActionAtomistic::getAbsoluteIndex(int i)const {
     248   127986465 :   return indexes[i];
     249             : }
     250             : 
     251             : inline
     252             : const std::vector<Vector> & ActionAtomistic::getPositions()const {
     253      561704 :   return positions;
     254             : }
     255             : 
     256             : inline
     257             : const double & ActionAtomistic::getEnergy()const {
     258        3989 :   return energy;
     259             : }
     260             : 
     261             : inline
     262             : const Tensor & ActionAtomistic::getBox()const {
     263       25194 :   return pbc.getBox();
     264             : }
     265             : 
     266             : inline
     267             : std::vector<Vector> & ActionAtomistic::modifyForces() {
     268      119080 :   return forces;
     269             : }
     270             : 
     271             : inline
     272             : Tensor & ActionAtomistic::modifyVirial() {
     273      133545 :   return virial;
     274             : }
     275             : 
     276             : inline
     277             : double & ActionAtomistic::modifyForceOnEnergy() {
     278             :   return forceOnEnergy;
     279             : }
     280             : 
     281             : inline
     282             : double & ActionAtomistic::modifyForceOnExtraCV() {
     283             :   return forceOnExtraCV;
     284             : }
     285             : 
     286             : inline
     287             : const Pbc & ActionAtomistic::getPbc() const {
     288      629579 :   return pbc;
     289             : }
     290             : 
     291             : inline
     292      160485 : void ActionAtomistic::lockRequests() {
     293      175292 :   lockRequestAtoms=true;
     294      160485 : }
     295             : 
     296             : inline
     297      160485 : void ActionAtomistic::unlockRequests() {
     298      175292 :   lockRequestAtoms=false;
     299      160485 : }
     300             : 
     301             : inline
     302             : const std::vector<AtomNumber> & ActionAtomistic::getUnique()const {
     303         546 :   return unique;
     304             : }
     305             : 
     306             : inline
     307             : const std::vector<AtomNumber> & ActionAtomistic::getUniqueLocal()const {
     308       56782 :   return unique_local;
     309             : }
     310             : 
     311             : inline
     312             : unsigned ActionAtomistic::getTotAtoms()const {
     313       32095 :   return atoms.positions.size();
     314             : }
     315             : 
     316             : inline
     317             : Pbc & ActionAtomistic::modifyGlobalPbc() {
     318          77 :   return atoms.pbc;
     319             : }
     320             : 
     321             : inline
     322             : void ActionAtomistic::setExtraCV(const std::string &name) {
     323           3 :   extraCV=name;
     324           3 : }
     325             : 
     326             : 
     327             : 
     328             : }
     329             : 
     330             : #endif