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