forked from plumed/plumed2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Torsion.cpp
86 lines (65 loc) · 3.09 KB
/
Torsion.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Copyright (c) 2013 The plumed team
(see the PEOPLE file at the root of the distribution for a list of names)
See http://www.plumed-code.org for more information.
This file is part of plumed, version 2.0.
plumed is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
plumed is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "Torsion.h"
#include "Tensor.h"
#include <cmath>
#include <iostream>
namespace PLMD{
double Torsion::compute(const Vector& v1,const Vector& v2,const Vector& v3)const{
const Vector nv2(v2*(1.0/v2.modulo()));
const Vector a(crossProduct(nv2,v1));
const Vector b(crossProduct(v3,nv2));
const double cosangle=dotProduct(a,b);
const double sinangle=dotProduct(crossProduct(a,b),nv2);
return std::atan2(-sinangle,cosangle);
}
double Torsion::compute(const Vector& v1,const Vector& v2,const Vector& v3,Vector& d1,Vector& d2,Vector& d3)const{
d1.zero();
d2.zero();
d3.zero();
const double modv2(v2.modulo());
const Vector nv2(v2/modv2);
const Tensor dnv2_v2((Tensor::identity()-extProduct(nv2,nv2))/modv2);
const Vector a(crossProduct(v2,v1));
const Tensor da_dv2(dcrossDv1(v2,v1));
const Tensor da_dv1(dcrossDv2(v2,v1));
const Vector b(crossProduct(v3,v2));
const Tensor db_dv3(dcrossDv1(v3,v2));
const Tensor db_dv2(dcrossDv2(v3,v2));
const double cosangle=dotProduct(a,b);
const Vector dcosangle_dv1=matmul(b,da_dv1);
const Vector dcosangle_dv2=matmul(b,da_dv2) + matmul(a,db_dv2);
const Vector dcosangle_dv3=matmul(a,db_dv3);
const Vector cab(crossProduct(a,b));
const Tensor dcab_dv1(matmul(dcrossDv1(a,b),da_dv1));
const Tensor dcab_dv2(matmul(dcrossDv1(a,b),da_dv2) + matmul(dcrossDv2(a,b),db_dv2));
const Tensor dcab_dv3(matmul(dcrossDv2(a,b),db_dv3));
const double sinangle=dotProduct(cab,nv2);
const Vector dsinangle_dv1=matmul(nv2,dcab_dv1);
const Vector dsinangle_dv2=matmul(nv2,dcab_dv2)+matmul(cab,dnv2_v2);
const Vector dsinangle_dv3=matmul(nv2,dcab_dv3);
const double x=-sinangle/cosangle;
const Vector dx_dv1=-dsinangle_dv1/cosangle + sinangle/(cosangle*cosangle) * dcosangle_dv1;
const Vector dx_dv2=-dsinangle_dv2/cosangle + sinangle/(cosangle*cosangle) * dcosangle_dv2;
const Vector dx_dv3=-dsinangle_dv3/cosangle + sinangle/(cosangle*cosangle) * dcosangle_dv3;
const double torsion=std::atan2(-sinangle,cosangle);
d1=1.0/(1.0+x*x) * dx_dv1;
d2=1.0/(1.0+x*x) * dx_dv2;
d3=1.0/(1.0+x*x) * dx_dv3;
return torsion;
}
}