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