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_tools_Torsion_h
      23             : #define __PLUMED_tools_Torsion_h
      24             : 
      25             : #include "Tensor.h"
      26             : 
      27             : namespace PLMD {
      28             : 
      29             : /// \ingroup TOOLBOX
      30             : /// Class to compute torsional angles.
      31             : /// I define it as a class even if it does not contain anything. The reason
      32             : /// is that in the future I would like to extend it to contain options about
      33             : /// how the calculation should be done. So, for now use it as
      34             : /// Torsion t;
      35             : /// T angle=t.compute(v1,v2,v3);
      36             : /// I know it is a bit misleading. If we really do not need to store "options"
      37             : /// inside the Torsion class, we can remove it later and write compute as
      38             : /// a static function.
      39             : class Torsion {
      40             : // still empty, but may accommodate some options in the future
      41             : public:
      42             : /// Compute the angle between the projections of v1 and v3 on the plane orthogonal
      43             : /// to v2. To have a "normal" definition (= plumed1), use it as
      44             : /// compute(r01,r12,r23);
      45             : /// See ColvarTorsion for a practical usage...
      46             :   template <typename T>
      47             :   static T compute(const VectorTyped<T,3>& v1,
      48             :                    const VectorTyped<T,3>& v2,
      49             :                    const VectorTyped<T,3>& v3) {
      50             :     using V3=PLMD::VectorTyped<T,3>;
      51             : 
      52             :     const V3 nv2(v2*(1.0/v2.modulo()));
      53             :     const V3 a(crossProduct(nv2,v1));
      54             :     const V3 b(crossProduct(v3,nv2));
      55             :     const T cosangle=dotProduct(a,b);
      56             :     const T sinangle=dotProduct(crossProduct(a,b),nv2);
      57             :     return std::atan2(-sinangle,cosangle);
      58             :   }
      59             : 
      60             : /// This is the version which also computes the derivatives wrt the arguments.
      61             :   template <typename T>
      62     1112800 :   static T compute(const VectorTyped<T,3>& v1,
      63             :                    const VectorTyped<T,3>& v2,
      64             :                    const VectorTyped<T,3>& v3,
      65             :                    VectorTyped<T,3>& d1,
      66             :                    VectorTyped<T,3>& d2,
      67             :                    VectorTyped<T,3>& d3) {
      68             :     using V3=PLMD::VectorTyped<T,3>;
      69             :     using T33=PLMD::TensorTyped<T,3,3>;
      70             : 
      71     1112800 :     const T modv2(T(1.0)/v2.modulo());
      72             :     const V3 nv2(v2*modv2);
      73     1112800 :     const T33 dnv2_v2((T33::identity()-extProduct(nv2,nv2))*modv2);
      74             : 
      75     1112800 :     const V3 a(crossProduct(v2,v1));
      76             :     const T33 da_dv2(dcrossDv1(v2,v1));
      77             :     const T33 da_dv1(dcrossDv2(v2,v1));
      78     1112800 :     const V3 b(crossProduct(v3,v2));
      79             :     const T33 db_dv3(dcrossDv1(v3,v2));
      80             :     const T33 db_dv2(dcrossDv2(v3,v2));
      81             :     const T cosangle=dotProduct(a,b);
      82     1112800 :     const V3 dcosangle_dv1=matmul(b,da_dv1);
      83     1112800 :     const V3 dcosangle_dv2=matmul(b,da_dv2) + matmul(a,db_dv2);
      84     1112800 :     const V3 dcosangle_dv3=matmul(a,db_dv3);
      85             : 
      86     1112800 :     const V3 cab(crossProduct(a,b));
      87     1112800 :     const T33 dcab_dv1(matmul(dcrossDv1(a,b),da_dv1));
      88     1112800 :     const T33 dcab_dv2(matmul(dcrossDv1(a,b),da_dv2) + matmul(dcrossDv2(a,b),db_dv2));
      89     1112800 :     const T33 dcab_dv3(matmul(dcrossDv2(a,b),db_dv3));
      90             : 
      91             :     const T sinangle=dotProduct(cab,nv2);
      92     1112800 :     const V3 dsinangle_dv1=matmul(nv2,dcab_dv1);
      93     1112800 :     const V3 dsinangle_dv2=matmul(nv2,dcab_dv2)+matmul(cab,dnv2_v2);
      94     1112800 :     const V3 dsinangle_dv3=matmul(nv2,dcab_dv3);
      95             : 
      96     1112800 :     const T torsion=std::atan2(-sinangle,cosangle);
      97             : // this is required since v1 and v3 are not normalized:
      98     1112800 :     const T invR2=T(1.0)/(cosangle*cosangle+sinangle*sinangle);
      99             : 
     100     1112800 :     d1= ( -dsinangle_dv1*cosangle + sinangle * dcosangle_dv1 ) *invR2;
     101     1112800 :     d2= ( -dsinangle_dv2*cosangle + sinangle * dcosangle_dv2 ) *invR2;
     102     1112800 :     d3= ( -dsinangle_dv3*cosangle + sinangle * dcosangle_dv3 ) *invR2;
     103             : 
     104     1112800 :     return torsion;
     105             :   }
     106             : 
     107             : };
     108             : 
     109             : }
     110             : 
     111             : #endif