Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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_multicolvar_MultiColvarBase_h
      23             : #define __PLUMED_multicolvar_MultiColvarBase_h
      24             : 
      25             : #include "core/ActionAtomistic.h"
      26             : #include "core/ActionWithValue.h"
      27             : #include "tools/DynamicList.h"
      28             : #include "tools/LinkCells.h"
      29             : #include "vesselbase/StoreDataVessel.h"
      30             : #include "vesselbase/ActionWithVessel.h"
      31             : #include "CatomPack.h"
      32             : #include <vector>
      33             : 
      34             : namespace PLMD {
      35             : namespace multicolvar {
      36             : 
      37             : class AtomValuePack;
      38             : class BridgedMultiColvarFunction;
      39             : class ActionVolume;
      40             : 
      41             : class MultiColvarBase :
      42             :   public ActionAtomistic,
      43             :   public ActionWithValue,
      44             :   public vesselbase::ActionWithVessel
      45             : {
      46             :   friend class BridgedMultiColvarFunction;
      47             :   friend class VolumeGradientBase;
      48             :   friend class MultiColvarFilter;
      49             :   friend class AtomValuePack;
      50             : private:
      51             : /// Use periodic boundary conditions
      52             :   bool usepbc;
      53             : /// The forces we are going to apply to things
      54             :   std::vector<double> forcesToApply;
      55             : /// We use this to say that all the atoms in the third block should are in the tasks
      56             :   bool allthirdblockintasks;
      57             : /// In certain cases we can make three atom link cells faster
      58             :   bool uselinkforthree;
      59             : /// Number of atoms that are active on this step
      60             :   unsigned nactive_atoms;
      61             : /// Stuff for link cells - this is used to make coordination number like variables faster
      62             :   LinkCells linkcells;
      63             : /// Link cells for third block of atoms
      64             :   LinkCells threecells;
      65             : /// Number of atoms that are being used for central atom position
      66             :   unsigned ncentral;
      67             : /// Bool vector telling us which atoms are required to calculate central atom position
      68             :   std::vector<bool> use_for_central_atom;
      69             : /// Vector of tempory holders for central atom values
      70             :   std::vector<CatomPack> my_tmp_capacks;
      71             : /// 1/number of atoms involved in central atoms
      72             :   double numberForCentralAtom;
      73             : /// Ensures that setup is only performed once per loop
      74             :   bool setup_completed;
      75             : /// Ensures that retrieving of atoms is only done once per calculation loop
      76             :   bool atomsWereRetrieved;
      77             : /// Add derivatives of center of mass position
      78             :   void addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const ;
      79             : protected:
      80             : /// This is used to keep track of what is calculated where
      81             :   std::vector< std::pair<unsigned,unsigned> > atom_lab;
      82             : /// The vessels in these multicolvars in which the data is stored
      83             :   std::vector<vesselbase::StoreDataVessel*> mybasedata;
      84             : /// The multicolvars from which we construct these quantities
      85             :   std::vector<MultiColvarBase*> mybasemulticolvars;
      86             : /// This remembers where the boundaries are for the tasks. It makes link cells work fast
      87             :   Matrix<std::pair<unsigned,unsigned> > bookeeping;
      88             : /// Function that recursively checks if filters have been used in the input to a multicolvar
      89             : /// we need this to ensure that setupLinkCells is run in calculate with some actions
      90             :   bool filtersUsedAsInput();
      91             : /// This resizes the arrays that are used for link cell update
      92             :   void resizeBookeepingArray( const unsigned& num1, const unsigned& num2 );
      93             : /// Are we doing sums of matrix rows
      94             :   bool matsums;
      95             : /// Using the species keyword to read in atoms
      96             :   bool usespecies;
      97             : /// Number of atoms in each block
      98             :   unsigned nblock;
      99             : /// This is used when turning cvcodes into atom numbers
     100             :   std::vector<unsigned> decoder;
     101             : /// Blocks of atom numbers
     102             :   std::vector< std::vector<unsigned> > ablocks;
     103             : /// Add a task to the list of tasks
     104             :   void addTaskToList( const unsigned& taskCode );
     105             : /// Finish setting up the multicolvar base
     106             :   void setupMultiColvarBase( const std::vector<AtomNumber>& atoms );
     107             : /// This routine take the vector of input derivatives and adds all the vectors to ivalth output derivatives
     108             : /// In other words end-start sets of derivatives come in and one set of derivatives come out
     109             :   void mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end, const unsigned& jatom,
     110             :                               const std::vector<double>& der, MultiValue& myder, AtomValuePack& myatoms ) const ;
     111             : /// This routine take the ith set of input derivatives and adds it to each of the (end-start) output derivatives
     112             : /// In other words one set of derivatives comes in and end-start sets of derivatives come out
     113             :   void splitInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
     114             :                               const unsigned& jatom, const std::vector<double>& der,
     115             :                               MultiValue& myder, AtomValuePack& myatoms ) const ;
     116             : /// This is used to accumulate functions of the coordination sphere.  Ensures weights are taken into account
     117             :   void accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const ;
     118             : /// Set which atoms are to be used to calculate the central atom position
     119             :   void setAtomsForCentralAtom( const std::vector<bool>& catom_ind );
     120             : /// Set the value of the cutoff for the link cells
     121             :   void setLinkCellCutoff( const double& lcut, double tcut=-1.0 );
     122             : /// Setup the link cells and neighbour list stuff
     123             :   void setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label );
     124             : /// Setup link cells in order to make this calculation faster
     125             :   void setupLinkCells();
     126             : /// Get the cutoff for the link cells
     127             :   double getLinkCellCutoff()  const ;
     128             : /// This does setup of link cell stuff that is specific to the non-use of the usespecies keyword
     129             :   void setupNonUseSpeciesLinkCells( const unsigned& );
     130             : /// This sets up the list of atoms that are involved in this colvar
     131             :   bool setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const ;
     132             : /// Decode indices if there are 2 or 3 atoms involved
     133             :   void decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const ;
     134             : /// Read in some atoms
     135             :   bool parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t);
     136             : /// Read in ATOMS keyword
     137             :   void readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms );
     138             : /// Read in two groups of atoms and setup multicolvars to calculate
     139             :   void readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms );
     140             : /// Read in three groups of atoms
     141             :   void readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
     142             :                           const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms );
     143             : /// Read in three groups of atoms and construct CVs involving at least one
     144             :   void readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
     145             :                         const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms );
     146             : /// Build sets by taking one multicolvar from each base
     147             :   void buildSets();
     148             : public:
     149             :   explicit MultiColvarBase(const ActionOptions&);
     150        1400 :   ~MultiColvarBase() {}
     151             :   static void registerKeywords( Keywords& keys );
     152             : /// Turn on the derivatives
     153             :   void turnOnDerivatives() override;
     154             : /// Get the separation between a pair of vectors
     155             :   Vector getSeparation( const Vector& vec1, const Vector& vec2 ) const ;
     156             : /// Do we use pbc to calculate this quantity
     157             :   bool usesPbc() const ;
     158             : /// Apply PBCs over a set of distance vectors
     159             :   void applyPbc(std::vector<Vector>& dlist, unsigned max_index=0) const;
     160             : /// Is it safe to use multithreading
     161        4461 :   bool threadSafe() const override { return !(mybasemulticolvars.size()>0); }
     162             : /// Do some setup before the calculation
     163             :   void prepare() override;
     164             : /// This is overwritten here in order to make sure that we do not retrieve atoms multiple times
     165             :   void retrieveAtoms() override;
     166             : /// Do the calculation
     167             :   void calculate() override;
     168             : /// Calculate numerical derivatives
     169             :   void calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
     170             : /// Perform one of the tasks
     171             :   void performTask( const unsigned&, const unsigned&, MultiValue& ) const override;
     172             : /// Update the active atoms
     173             :   virtual void updateActiveAtoms( AtomValuePack& myatoms ) const ;
     174             : /// This gets the position of an atom for the link cell setup
     175             :   virtual Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
     176             : /// Get the absolute index of the central atom
     177             :   virtual AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& i ) const ;
     178             : /// This is replaced once we have a function to calculate the cv
     179             :   virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const=0;
     180             : /// Apply the forces from this action
     181             :   void apply() override;
     182             : /// Get the number of derivatives for this action
     183             :   unsigned getNumberOfDerivatives() override;  // N.B. This is replacing the virtual function in ActionWithValue
     184             : /// Checks if an task is being performed at the present time
     185             :   virtual bool isCurrentlyActive( const unsigned& code );
     186             : /// Add some derivatives to a particular component of a particular atom
     187             :   void addAtomDerivatives( const int&, const unsigned&, const Vector&, multicolvar::AtomValuePack& ) const ;
     188             : ///
     189             :   virtual void getCentralAtomPack( const unsigned& basn, const unsigned& curr, CatomPack& mypack);
     190             : /// Get the index where the central atom is stored
     191             :   virtual Vector getCentralAtomPos( const unsigned& curr );
     192             : /// You can use this to screen contributions that are very small so we can avoid expensive (and pointless) calculations
     193             :   virtual double calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const ;
     194             : /// Is this a density?
     195     1814585 :   virtual bool isDensity() const { return false; }
     196             : /// Is the iatom'th stored value currently active
     197             :   bool storedValueIsActive( const unsigned& iatom );
     198             : /// This is true if multicolvar is calculating a vector or if the multicolvar is the density
     199           0 :   virtual bool hasDifferentiableOrientation() const { return false; }
     200             : /// This makes sure we are not calculating the director when we do LocalAverage
     201           0 :   virtual void doNotCalculateDirector() {}
     202             : /// Ensure that derivatives are only calculated when needed
     203             :   bool doNotCalculateDerivatives() const override;
     204             : /// Get the icolv th base multicolvar
     205             :   MultiColvarBase* getBaseMultiColvar( const unsigned& icolv ) const ;
     206             : /// Get the number of base multicolvars
     207             :   unsigned getNumberOfBaseMultiColvars() const ;
     208             : /// Get an input data
     209             :   virtual void getInputData( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms, std::vector<double>& orient ) const ;
     210             : /// Retrieve the input derivatives
     211             :   virtual MultiValue& getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const ;
     212             : };
     213             : 
     214             : inline
     215      917133 : bool MultiColvarBase::isCurrentlyActive( const unsigned& code ) {
     216      917133 :   if( setup_completed && atom_lab[code].first>0 ) {
     217       17644 :     unsigned mmc=atom_lab[code].first - 1;
     218       17644 :     return mybasedata[mmc]->storedValueIsActive( atom_lab[code].second );
     219             :   }
     220             :   return true;
     221             : }
     222             : 
     223             : inline
     224         156 : AtomNumber MultiColvarBase::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
     225             :   plumed_dbg_assert( iatom<atom_lab.size() );
     226         156 :   if( atom_lab[iatom].first>0  ) {
     227           0 :     unsigned mmc=atom_lab[iatom].first - 1;
     228           0 :     return mybasemulticolvars[mmc]->getAbsoluteIndexOfCentralAtom( atom_lab[iatom].second );
     229             :   }
     230             :   plumed_dbg_assert( usespecies );
     231         156 :   return ActionAtomistic::getAbsoluteIndex( atom_lab[getTaskCode(iatom)].second );
     232             : }
     233             : 
     234             : inline
     235   408128313 : Vector MultiColvarBase::getPositionOfAtomForLinkCells( const unsigned& iatom ) const {
     236             :   plumed_dbg_assert( iatom<atom_lab.size() );
     237   408128313 :   if( atom_lab[iatom].first>0  ) {
     238     2487982 :     unsigned mmc=atom_lab[iatom].first - 1;
     239     2487982 :     return mybasemulticolvars[mmc]->getCentralAtomPos( atom_lab[iatom].second );
     240             :   }
     241   405640331 :   return ActionAtomistic::getPosition( atom_lab[iatom].second );
     242             : }
     243             : 
     244             : inline
     245     8911885 : unsigned MultiColvarBase::getNumberOfDerivatives() {
     246     8911885 :   return 3*getNumberOfAtoms()+9;
     247             : }
     248             : 
     249             : inline
     250             : bool MultiColvarBase::usesPbc() const {
     251      220567 :   return usepbc;
     252             : }
     253             : 
     254             : inline
     255    69055090 : bool MultiColvarBase::doNotCalculateDerivatives() const {
     256    69055090 :   if( !dertime ) return true;
     257    68367370 :   return ActionWithValue::doNotCalculateDerivatives();
     258             : }
     259             : 
     260             : inline
     261             : unsigned MultiColvarBase::getNumberOfBaseMultiColvars() const {
     262         178 :   return mybasemulticolvars.size();
     263             : }
     264             : 
     265             : inline
     266             : MultiColvarBase* MultiColvarBase::getBaseMultiColvar( const unsigned& icolv ) const {
     267             :   plumed_dbg_assert( icolv<mybasemulticolvars.size() );
     268       25010 :   return mybasemulticolvars[icolv];
     269             : }
     270             : 
     271             : }
     272             : }
     273             : 
     274             : #endif