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_colvar_Angle_h
      23             : #define __PLUMED_colvar_Angle_h
      24             : #include "Colvar.h"
      25             : #include "ColvarInput.h"
      26             : #include "tools/Angle.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace colvar {
      30             : 
      31             : template <typename T>
      32             : class Angle : public Colvar {
      33             :   bool pbc;
      34             :   std::vector<double> value;
      35             :   std::vector<double> derivs;
      36             : public:
      37             :   using precision=T;
      38             :   explicit Angle(const ActionOptions&);
      39             : // active methods:
      40             :   void calculate() override;
      41             :   static void registerKeywords( Keywords& keys );
      42             :   static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
      43             :   static unsigned getModeAndSetupValues( ActionWithValue* av );
      44             :   static void calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout );
      45             : };
      46             : 
      47             : template <typename T>
      48         155 : void Angle<T>::registerKeywords( Keywords& keys ) {
      49         155 :   Colvar::registerKeywords(keys);
      50         310 :   keys.setDisplayName("ANGLE");
      51         155 :   keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
      52         155 :   keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      53         310 :   keys.setValueDescription("scalar/vector","the ANGLE involving these atoms");
      54         310 :   keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
      55         155 : }
      56             : 
      57             : template <typename T>
      58         146 : void Angle<T>::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
      59         292 :   aa->parseAtomList("ATOMS",num,atoms);
      60         146 :   if(atoms.size()==3) {
      61         128 :     aa->log.printf("  between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
      62         128 :     atoms.resize(4);
      63         128 :     atoms[3]=atoms[2];
      64         128 :     atoms[2]=atoms[1];
      65          18 :   } else if(atoms.size()==4) {
      66           7 :     aa->log.printf("  between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
      67          11 :   } else if( num<0 || atoms.size()>0 ) {
      68           1 :     aa->error("Number of specified atoms should be either 3 or 4");
      69             :   }
      70         145 : }
      71             : 
      72             : template <typename T>
      73          10 : unsigned Angle<T>::getModeAndSetupValues( ActionWithValue* av ) {
      74          10 :   av->addValueWithDerivatives();
      75          10 :   av->setNotPeriodic();
      76          10 :   return 0;
      77             : }
      78             : 
      79             : template <typename T>
      80          34 : Angle<T>::Angle(const ActionOptions&ao):
      81             :   PLUMED_COLVAR_INIT(ao),
      82          34 :   pbc(true),
      83          34 :   value(1) {
      84          34 :   std::vector<AtomNumber> atoms;
      85          34 :   parseAtomList( -1, atoms, this );
      86          33 :   bool nopbc=!pbc;
      87          33 :   parseFlag("NOPBC",nopbc);
      88          33 :   pbc=!nopbc;
      89             : 
      90          33 :   if(pbc) {
      91          33 :     log.printf("  using periodic boundary conditions\n");
      92             :   } else {
      93           0 :     log.printf("  without periodic boundary conditions\n");
      94             :   }
      95             : 
      96          33 :   addValueWithDerivatives();
      97          33 :   setNotPeriodic();
      98          33 :   requestAtoms(atoms);
      99          33 :   checkRead();
     100          36 : }
     101             : 
     102             : // calculator
     103             : template <typename T>
     104         502 : void Angle<T>::calculate() {
     105             : 
     106         502 :   if(pbc) {
     107         502 :     makeWhole();
     108             :   }
     109         502 :   ColvarOutput<double> cvout = ColvarOutput<double>::createColvarOutput(value,derivs,this);
     110         502 :   Angle<double>::calculateCV( ColvarInput<double>::createColvarInput( 0, getPositions(), this ), cvout );
     111         502 :   setValue( value[0] );
     112        2510 :   for(unsigned i=0; i<getPositions().size(); ++i) {
     113        2008 :     setAtomsDerivatives( i, cvout.getAtomDerivatives(0,i) );
     114             :   }
     115         502 :   setBoxDerivatives( cvout.virial[0] );
     116         502 : }
     117             : 
     118             : template <typename T>
     119        1337 : void Angle<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
     120        1337 :   auto dij=delta(cvin.pos[2],cvin.pos[3]);
     121        1337 :   auto dik=delta(cvin.pos[1],cvin.pos[0]);
     122        1337 :   VectorTyped<T,3> ddij,ddik;
     123             : 
     124        1337 :   cvout.values[0]=PLMD::Angle::compute(dij,dik,ddij,ddik);
     125             :   cvout.derivs[0][0]=ddik;
     126             :   cvout.derivs[0][1]=-ddik;
     127             :   cvout.derivs[0][2]=-ddij;
     128             :   cvout.derivs[0][3]=ddij;
     129        1337 :   ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
     130        1337 : }
     131             : 
     132             : } //namespace colvar
     133             : } //namespace PLMD
     134             : 
     135             : #endif //__PLUMED_colvar_Angle_h