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             : #include "Torsion.h"
      23             : #include "Tensor.h"
      24             : 
      25             : #include <cmath>
      26             : #include <iostream>
      27             : 
      28             : namespace PLMD {
      29             : 
      30           0 : double Torsion::compute(const Vector& v1,const Vector& v2,const Vector& v3)const {
      31           0 :   const Vector nv2(v2*(1.0/v2.modulo()));
      32           0 :   const Vector a(crossProduct(nv2,v1));
      33           0 :   const Vector b(crossProduct(v3,nv2));
      34           0 :   const double cosangle=dotProduct(a,b);
      35           0 :   const double sinangle=dotProduct(crossProduct(a,b),nv2);
      36           0 :   return std::atan2(-sinangle,cosangle);
      37             : }
      38             : 
      39     1084670 : double Torsion::compute(const Vector& v1,const Vector& v2,const Vector& v3,Vector& d1,Vector& d2,Vector& d3)const {
      40             : 
      41     1084670 :   const double modv2(1./v2.modulo());
      42     1084670 :   const Vector nv2(v2*modv2);
      43     1084670 :   const Tensor dnv2_v2((Tensor::identity()-extProduct(nv2,nv2))*modv2);
      44             : 
      45     1084670 :   const Vector a(crossProduct(v2,v1));
      46     1084670 :   const Tensor da_dv2(dcrossDv1(v2,v1));
      47     1084670 :   const Tensor da_dv1(dcrossDv2(v2,v1));
      48     1084670 :   const Vector b(crossProduct(v3,v2));
      49     1084670 :   const Tensor db_dv3(dcrossDv1(v3,v2));
      50     1084670 :   const Tensor db_dv2(dcrossDv2(v3,v2));
      51     1084670 :   const double cosangle=dotProduct(a,b);
      52     1084670 :   const Vector dcosangle_dv1=matmul(b,da_dv1);
      53     1084670 :   const Vector dcosangle_dv2=matmul(b,da_dv2) + matmul(a,db_dv2);
      54     1084670 :   const Vector dcosangle_dv3=matmul(a,db_dv3);
      55             : 
      56     1084670 :   const Vector cab(crossProduct(a,b));
      57     1084670 :   const Tensor dcab_dv1(matmul(dcrossDv1(a,b),da_dv1));
      58     1084670 :   const Tensor dcab_dv2(matmul(dcrossDv1(a,b),da_dv2) + matmul(dcrossDv2(a,b),db_dv2));
      59     1084670 :   const Tensor dcab_dv3(matmul(dcrossDv2(a,b),db_dv3));
      60             : 
      61     1084670 :   const double sinangle=dotProduct(cab,nv2);
      62     1084670 :   const Vector dsinangle_dv1=matmul(nv2,dcab_dv1);
      63     1084670 :   const Vector dsinangle_dv2=matmul(nv2,dcab_dv2)+matmul(cab,dnv2_v2);
      64     1084670 :   const Vector dsinangle_dv3=matmul(nv2,dcab_dv3);
      65             : 
      66     1084670 :   const double torsion=std::atan2(-sinangle,cosangle);
      67             : // this is required since v1 and v3 are not normalized:
      68     1084670 :   const double invR2=1.0/(cosangle*cosangle+sinangle*sinangle);
      69             : 
      70     1084670 :   d1= ( -dsinangle_dv1*cosangle + sinangle * dcosangle_dv1 ) *invR2;
      71     1084670 :   d2= ( -dsinangle_dv2*cosangle + sinangle * dcosangle_dv2 ) *invR2;
      72     1084670 :   d3= ( -dsinangle_dv3*cosangle + sinangle * dcosangle_dv3 ) *invR2;
      73             : 
      74     1084670 :   return torsion;
      75             : }
      76             : 
      77             : }
      78             : 
      79             : 
      80             :