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 "Angle.h"
23 #include "Tools.h"
24 #include <cmath>
25
26 namespace PLMD {
27
compute(const Vector & v1,const Vector & v2) const28 double Angle::compute(const Vector& v1,const Vector& v2)const {
29 return std::acos(dotProduct(v1,v2)/(v1.modulo()*v2.modulo()));
30 }
31
compute(const Vector & v1,const Vector & v2,Vector & d1,Vector & d2) const32 double Angle::compute(const Vector& v1,const Vector& v2,Vector& d1,Vector& d2)const {
33 const double dp(dotProduct(v1,v2));
34 const Vector& dp_dv1(v2);
35 const Vector& dp_dv2(v1);
36 const double sv1(v1.modulo2());
37 const double sv2(v2.modulo2());
38 const Vector dsv1_dv1(2*v1);
39 const Vector dsv2_dv2(2*v2);
40 const double nn(1.0/sqrt(sv1*sv2));
41 const Vector dnn_dv1(-0.5*nn/sv1*dsv1_dv1);
42 const Vector dnn_dv2(-0.5*nn/sv2*dsv2_dv2);
43
44 const double dpnn(dp*nn);
45
46 if(dpnn>=1.0-epsilon) {
47 d1=Vector(0.0,0.0,0.0);
48 d2=Vector(0.0,0.0,0.0);
49 return 0.0;
50 }
51 if(dpnn<=-1.0+epsilon) {
52 d1=Vector(0.0,0.0,0.0);
53 d2=Vector(0.0,0.0,0.0);
54 return pi;
55 }
56 const Vector ddpnn_dv1(dp*dnn_dv1+dp_dv1*nn);
57 const Vector ddpnn_dv2(dp*dnn_dv2+dp_dv2*nn);
58
59 const double x(-1.0/sqrt(1-dpnn*dpnn));
60
61 d1=x*ddpnn_dv1;
62 d2=x*ddpnn_dv2;
63
64
65 return std::acos(dpnn);
66 }
67
68 }
69