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_Torsion_h
      23             : #define __PLUMED_colvar_Torsion_h
      24             : #include "Colvar.h"
      25             : #include "tools/Torsion.h"
      26             : #include "ColvarInput.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace colvar {
      30             : 
      31             : template <typename T=double>
      32             : class Torsion : public Colvar {
      33             :   bool pbc;
      34             :   bool do_cosine;
      35             : 
      36             :   std::vector<double> value;
      37             :   std::vector<double> derivs;
      38             : public:
      39             :   using precision=T;
      40             :   explicit Torsion(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             :   static void registerKeywords(Keywords& keys);
      47             : };
      48             : 
      49             : template <typename T>
      50        2506 : void Torsion<T>::registerKeywords(Keywords& keys) {
      51        2506 :   Colvar::registerKeywords( keys );
      52        5012 :   keys.setDisplayName("TORSION");
      53        2506 :   keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
      54        2506 :   keys.add("atoms-2","AXIS","two atoms that define an axis.  You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords.");
      55        2506 :   keys.add("atoms-2","VECTORA","two atoms that define a vector.  You can use this in combination with VECTOR2 and AXIS");
      56        2506 :   keys.add("atoms-2","VECTORB","two atoms that define a vector.  You can use this in combination with VECTOR1 and AXIS");
      57        2506 :   keys.addDeprecatedKeyword("VECTOR1","VECTORA");
      58        2506 :   keys.addDeprecatedKeyword("VECTOR2","VECTORB");
      59        2506 :   keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
      60        2506 :   keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      61        5012 :   keys.setValueDescription("scalar/vector","the TORSION involving these atoms");
      62        5012 :   keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
      63        2506 : }
      64             : 
      65             : template <typename T>
      66         806 : Torsion<T>::Torsion(const ActionOptions&ao):
      67             :   PLUMED_COLVAR_INIT(ao),
      68         806 :   pbc(true),
      69         806 :   do_cosine(false),
      70         806 :   value(1) {
      71         806 :   std::vector<AtomNumber> atoms;
      72         806 :   std::vector<AtomNumber> v1;
      73        1612 :   ActionAtomistic::parseAtomList("VECTOR1",v1);
      74         806 :   if( v1.size()>0 ) {
      75           0 :     std::vector<AtomNumber> v2;
      76           0 :     ActionAtomistic::parseAtomList("VECTOR2",v2);
      77           0 :     std::vector<AtomNumber> axis;
      78           0 :     ActionAtomistic::parseAtomList("AXIS",axis);
      79           0 :     if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
      80           0 :       error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
      81             :     }
      82           0 :     atoms.resize(6);
      83           0 :     atoms[0]=v1[1];
      84           0 :     atoms[1]=v1[0];
      85           0 :     atoms[2]=axis[0];
      86           0 :     atoms[3]=axis[1];
      87           0 :     atoms[4]=v2[0];
      88           0 :     atoms[5]=v2[1];
      89           0 :     log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
      90             :                v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
      91           0 :   } else {
      92         806 :     parseAtomList(-1,atoms,this);
      93             :   }
      94         804 :   unsigned mode=getModeAndSetupValues(this);
      95         804 :   if( mode==1 ) {
      96           3 :     do_cosine=true;
      97             :   }
      98             : 
      99         804 :   bool nopbc=!pbc;
     100         804 :   parseFlag("NOPBC",nopbc);
     101         804 :   pbc=!nopbc;
     102         804 :   checkRead();
     103             : 
     104         804 :   if(pbc) {
     105         680 :     log.printf("  using periodic boundary conditions\n");
     106             :   } else {
     107         124 :     log.printf("  without periodic boundary conditions\n");
     108             :   }
     109         804 :   requestAtoms(atoms);
     110         812 : }
     111             : 
     112             : template <typename T>
     113         950 : void Torsion<T>::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
     114         950 :   std::vector<AtomNumber> v1,v2,axis;
     115         950 :   aa->parseAtomList("ATOMS",num,t);
     116         950 :   aa->parseAtomList("VECTORA",num,v1);
     117         950 :   aa->parseAtomList("VECTORB",num,v2);
     118        1900 :   aa->parseAtomList("AXIS",num,axis);
     119             : 
     120         950 :   if(t.size()==4) {
     121         884 :     if(!(v1.empty() && v2.empty() && axis.empty())) {
     122           1 :       aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords");
     123             :     }
     124         883 :     aa->log.printf("  between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial());
     125         883 :     t.resize(6);
     126         883 :     t[5]=t[3];
     127         883 :     t[4]=t[2];
     128         883 :     t[3]=t[2];
     129         883 :     t[2]=t[1];
     130          66 :   } else if(t.empty()) {
     131          65 :     if( num>0 && v1.empty() && v2.empty() && axis.empty() ) {
     132             :       return;
     133             :     }
     134          43 :     if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
     135           0 :       aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each");
     136             :     }
     137          43 :     aa->log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
     138             :                    v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
     139          43 :     t.resize(6);
     140          43 :     t[0]=v1[1];
     141          43 :     t[1]=v1[0];
     142          43 :     t[2]=axis[0];
     143          43 :     t[3]=axis[1];
     144          43 :     t[4]=v2[0];
     145          43 :     t[5]=v2[1];
     146             :   } else if( t.size()!=4 ) {
     147           1 :     aa->error("ATOMS should specify 4 atoms");
     148             :   }
     149         954 : }
     150             : 
     151             : template <typename T>
     152         826 : unsigned Torsion<T>::getModeAndSetupValues( ActionWithValue* av ) {
     153             :   bool do_cos;
     154         826 :   av->parseFlag("COSINE",do_cos);
     155         826 :   if(do_cos) {
     156           4 :     av->log.printf("  calculating cosine instead of torsion\n");
     157             :   }
     158             : 
     159         826 :   av->addValueWithDerivatives();
     160         826 :   if(!do_cos) {
     161        1644 :     av->setPeriodic("-pi","pi");
     162         822 :     return 0;
     163             :   }
     164           4 :   av->setNotPeriodic();
     165             :   return 1;
     166             : }
     167             : 
     168             : // calculator
     169             : template <typename T>
     170       54207 : void Torsion<T>::calculate() {
     171       54207 :   if(pbc) {
     172       51639 :     makeWhole();
     173             :   }
     174       54207 :   auto cvout = ColvarOutput<double>::createColvarOutput(value,derivs,this);
     175       54207 :   if(do_cosine) {
     176          15 :     Torsion<double>::calculateCV(
     177          30 :       ColvarInput<double>::createColvarInput( 1, getPositions(), this ),
     178             :       cvout );
     179             :   } else {
     180       54192 :     Torsion<double>::calculateCV(
     181      108384 :       ColvarInput<double>::createColvarInput( 0, getPositions(), this ),
     182             :       cvout );
     183             :   }
     184      379449 :   for(unsigned i=0; i<6; ++i) {
     185      325242 :     setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
     186             :   }
     187       54207 :   setValue(value[0]);
     188       54207 :   setBoxDerivatives( cvout.virial[0] );
     189       54207 : }
     190             : 
     191             : template <typename T>
     192       55582 : void Torsion<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
     193       55582 :   auto d0=delta(cvin.pos[1],cvin.pos[0]);
     194       55582 :   auto d1=delta(cvin.pos[3],cvin.pos[2]);
     195       55582 :   auto d2=delta(cvin.pos[5],cvin.pos[4]);
     196       55582 :   VectorTyped<T,3> dd0,dd1,dd2;
     197       55582 :   cvout.values[0] = PLMD::Torsion::compute(d0,d1,d2,dd0,dd1,dd2);
     198       55582 :   if(cvin.mode==1) {
     199          30 :     dd0 *= -std::sin(cvout.values[0]);
     200             :     dd1 *= -std::sin(cvout.values[0]);
     201             :     dd2 *= -std::sin(cvout.values[0]);
     202          30 :     cvout.values[0] = std::cos(cvout.values[0]);
     203             :   }
     204             :   cvout.derivs[0][0] = dd0;
     205             :   cvout.derivs[0][1] = -dd0;
     206             :   cvout.derivs[0][2] = dd1;
     207             :   cvout.derivs[0][3] = -dd1;
     208             :   cvout.derivs[0][4] = dd2;
     209             :   cvout.derivs[0][5] = -dd2;
     210       55582 :   ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
     211       55582 : }
     212             : 
     213             : }
     214             : }
     215             : #endif // __PLUMED_colvar_Torsion.h