Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-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_Position_h
      23             : #define __PLUMED_colvar_Position_h
      24             : #include "Colvar.h"
      25             : #include "tools/Pbc.h"
      26             : #include "ColvarInput.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace colvar {
      30             : 
      31             : template <typename T>
      32             : class Position : public Colvar {
      33             :   enum components:unsigned {scaled=1,standard=0} mode{standard};
      34             :   bool pbc;
      35             :   std::vector<double> value;
      36             :   std::vector<double> derivs;
      37             : public:
      38             :   using precision=T;
      39             :   static void registerKeywords( Keywords& keys );
      40             :   explicit Position(const ActionOptions&);
      41             :   static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
      42             :   static unsigned getModeAndSetupValues( ActionWithValue* av );
      43             : // active methods:
      44             :   void calculate() override;
      45             :   static void calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout );
      46             : };
      47             : 
      48             : template <typename T>
      49         505 : void Position<T>::registerKeywords( Keywords& keys ) {
      50         505 :   Colvar::registerKeywords( keys );
      51        1010 :   keys.setDisplayName("POSITION");
      52         505 :   keys.add("atoms","ATOM","the atom number");
      53         505 :   keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of");
      54         505 :   keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before");
      55         505 :   keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c");
      56        1010 :   keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position");
      57        1010 :   keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position");
      58        1010 :   keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position");
      59        1010 :   keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position");
      60        1010 :   keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position");
      61        1010 :   keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position");
      62         505 :   keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      63        1010 :   keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
      64         505 : }
      65             : 
      66             : template <typename T>
      67          95 : Position<T>::Position(const ActionOptions&ao):
      68             :   PLUMED_COLVAR_INIT(ao),
      69          95 :   pbc(true),
      70          95 :   value(3) {
      71          95 :   std::vector<AtomNumber> atoms;
      72          95 :   parseAtomList(-1,atoms,this);
      73          94 :   mode=static_cast<components>(getModeAndSetupValues(this));
      74             : 
      75          94 :   bool nopbc=!pbc;
      76          94 :   parseFlag("NOPBC",nopbc);
      77          94 :   pbc=!nopbc;
      78          94 :   checkRead();
      79             : 
      80          94 :   if(pbc) {
      81          88 :     log.printf("  using periodic boundary conditions\n");
      82             :   } else {
      83           6 :     log.printf("  without periodic boundary conditions\n");
      84             :   }
      85             : 
      86          94 :   requestAtoms(atoms);
      87          97 : }
      88             : 
      89             : template <typename T>
      90         103 : void Position<T>::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
      91         206 :   aa->parseAtomList("ATOM",num,t);
      92         103 :   if( t.size()==1 ) {
      93         100 :     aa->log.printf("  for atom %d\n",t[0].serial());
      94           3 :   } else if( num<0 || t.size()!=0 ) {
      95           1 :     aa->error("Number of specified atoms should be 1");
      96             :   }
      97         102 : }
      98             : 
      99             : template <typename T>
     100         142 : unsigned Position<T>::getModeAndSetupValues( ActionWithValue* av ) {
     101             :   bool sc;
     102         142 :   av->parseFlag("SCALED_COMPONENTS",sc);
     103         142 :   if(sc) {
     104          30 :     av->addComponentWithDerivatives("a");
     105          30 :     av->componentIsPeriodic("a","-0.5","+0.5");
     106          30 :     av->addComponentWithDerivatives("b");
     107          30 :     av->componentIsPeriodic("b","-0.5","+0.5");
     108          30 :     av->addComponentWithDerivatives("c");
     109          30 :     av->componentIsPeriodic("c","-0.5","+0.5");
     110          15 :     return components::scaled;
     111             :   }
     112         254 :   av->addComponentWithDerivatives("x");
     113         127 :   av->componentIsNotPeriodic("x");
     114         254 :   av->addComponentWithDerivatives("y");
     115         127 :   av->componentIsNotPeriodic("y");
     116         254 :   av->addComponentWithDerivatives("z");
     117         127 :   av->componentIsNotPeriodic("z");
     118         127 :   av->log<<"  WARNING: components will not have the proper periodicity - see manual\n";
     119             :   return components::standard;
     120             : }
     121             : 
     122             : // calculator
     123             : template <typename T>
     124        8078 : void Position<T>::calculate() {
     125             : 
     126        8078 :   std::vector<Vector> distance(1);
     127        8078 :   if(pbc) {
     128        8022 :     distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
     129             :   } else {
     130          56 :     distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
     131             :   }
     132             : 
     133        8078 :   ColvarOutput<double> cvout = ColvarOutput<double>::createColvarOutput(value,derivs,this);
     134        8078 :   Position<double>::calculateCV( ColvarInput<double>::createColvarInput( mode, distance, this ), cvout );
     135        8078 :   if(mode==components::scaled) {
     136          56 :     Value* valuea=getPntrToComponent("a");
     137          56 :     Value* valueb=getPntrToComponent("b");
     138         112 :     Value* valuec=getPntrToComponent("c");
     139          56 :     setAtomsDerivatives (valuea,0,cvout.getAtomDerivatives(0,0));
     140          56 :     valuea->set(value[0]);
     141          56 :     setAtomsDerivatives (valueb,0,cvout.getAtomDerivatives(1,0));
     142          56 :     valueb->set(value[1]);
     143          56 :     setAtomsDerivatives (valuec,0,cvout.getAtomDerivatives(2,0));
     144          56 :     valuec->set(value[2]);
     145             :   } else {
     146        8022 :     Value* valuex=getPntrToComponent("x");
     147        8022 :     Value* valuey=getPntrToComponent("y");
     148       16044 :     Value* valuez=getPntrToComponent("z");
     149             : 
     150        8022 :     setAtomsDerivatives (valuex,0,cvout.getAtomDerivatives(0,0));
     151        8022 :     setBoxDerivatives   (valuex,cvout.virial[0]);
     152        8022 :     valuex->set(value[0]);
     153             : 
     154        8022 :     setAtomsDerivatives (valuey,0,cvout.getAtomDerivatives(1,0));
     155        8022 :     setBoxDerivatives   (valuey,cvout.virial[1]);
     156        8022 :     valuey->set(value[1]);
     157             : 
     158        8022 :     setAtomsDerivatives (valuez,0,cvout.getAtomDerivatives(2,0));
     159        8022 :     setBoxDerivatives   (valuez,cvout.virial[2]);
     160        8022 :     valuez->set(value[2]);
     161             :   }
     162        8078 : }
     163             : 
     164             : template <typename T>
     165      195645 : void Position<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
     166      195645 :   if( cvin.mode==components::scaled ) {
     167       11865 :     auto d=cvin.pbc.realToScaled(VectorTyped<T,3>(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]));
     168       11865 :     cvout.values[0]=Tools::pbc(d[0]);
     169       11865 :     cvout.values[1]=Tools::pbc(d[1]);
     170       11865 :     cvout.values[2]=Tools::pbc(d[2]);
     171       11865 :     cvout.derivs[0][0]=matmul(cvin.pbc.getInvBox(),Versors::xp<T>);
     172       11865 :     cvout.derivs[1][0]=matmul(cvin.pbc.getInvBox(),Versors::yp<T>);
     173       11865 :     cvout.derivs[2][0]=matmul(cvin.pbc.getInvBox(),Versors::zp<T>);
     174             :   } else {
     175      735120 :     for(unsigned i=0; i<3; ++i) {
     176      551340 :       cvout.values[i]=cvin.pos[0][i];
     177             :     }
     178             :     cvout.derivs[0][0]=Versors::xp<T>;
     179             :     cvout.derivs[1][0]=Versors::yp<T>;
     180             :     cvout.derivs[2][0]=Versors::zp<T>;
     181      183780 :     cvout.virial.set(0, TensorTyped<T,3,3>(VectorTyped<T,3>(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),
     182             :                                            Versors::xm<T>));
     183      183780 :     cvout.virial.set(1, TensorTyped<T,3,3>(VectorTyped<T,3>(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),
     184             :                                            Versors::ym<T>));
     185      183780 :     cvout.virial.set(2, TensorTyped<T,3,3>(VectorTyped<T,3>(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),
     186             :                                            Versors::zm<T>));
     187             :   }
     188      195645 : }
     189             : 
     190             : }
     191             : }
     192             : #endif //__PLUMED_colvar_Position_h