Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73577509
colvartypes.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Jul 22, 20:37
Size
41 KB
Mime Type
text/x-c++
Expires
Wed, Jul 24, 20:37 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19226808
Attached To
rLAMMPS lammps
colvartypes.h
View Options
// -*- c++ -*-
#ifndef COLVARTYPES_H
#define COLVARTYPES_H
#include <cmath>
#include <vector>
#include "colvarmodule.h"
#ifndef PI
#define PI 3.14159265358979323846
#endif
// ----------------------------------------------------------------------
/// Linear algebra functions and data types used in the collective
/// variables implemented so far
// ----------------------------------------------------------------------
/// \brief Arbitrary size array (one dimensions) suitable for linear
/// algebra operations (i.e. for floating point numbers it can be used
/// with library functions)
template
<
class
T
>
class
colvarmodule
::
vector1d
{
protected:
std
::
vector
<
T
>
data
;
public:
/// Default constructor
inline
vector1d
(
size_t
const
n
=
0
)
{
data
.
resize
(
n
);
reset
();
}
/// Constructor from C array
inline
vector1d
(
size_t
const
n
,
T
const
*
t
)
{
data
.
resize
(
n
);
reset
();
size_t
i
;
for
(
i
=
0
;
i
<
size
();
i
++
)
{
data
[
i
]
=
t
[
i
];
}
}
/// Return a pointer to the data location
inline
T
*
c_array
()
{
if
(
data
.
size
()
>
0
)
{
return
&
(
data
[
0
]);
}
else
{
return
NULL
;
}
}
inline
~
vector1d
()
{
data
.
clear
();
}
/// Set all elements to zero
inline
void
reset
()
{
data
.
assign
(
data
.
size
(),
T
(
0.0
));
}
inline
size_t
size
()
const
{
return
data
.
size
();
}
inline
void
resize
(
size_t
const
n
)
{
data
.
resize
(
n
);
}
inline
T
&
operator
[]
(
size_t
const
i
)
{
return
data
[
i
];
}
inline
T
const
&
operator
[]
(
size_t
const
i
)
const
{
return
data
[
i
];
}
inline
static
void
check_sizes
(
vector1d
<
T
>
const
&
v1
,
vector1d
<
T
>
const
&
v2
)
{
if
(
v1
.
size
()
!=
v2
.
size
())
{
cvm
::
error
(
"Error: trying to perform an operation between vectors of different sizes, "
+
cvm
::
to_str
(
v1
.
size
())
+
" and "
+
cvm
::
to_str
(
v2
.
size
())
+
".
\n
"
);
}
}
inline
void
operator
+=
(
vector1d
<
T
>
const
&
v
)
{
check_sizes
(
*
this
,
v
);
size_t
i
;
for
(
i
=
0
;
i
<
this
->
size
();
i
++
)
{
(
*
this
)[
i
]
+=
v
[
i
];
}
}
inline
void
operator
-=
(
vector1d
<
T
>
const
&
v
)
{
check_sizes
(
*
this
,
v
);
size_t
i
;
for
(
i
=
0
;
i
<
this
->
size
();
i
++
)
{
(
*
this
)[
i
]
-=
v
[
i
];
}
}
inline
void
operator
*=
(
cvm
::
real
const
&
a
)
{
size_t
i
;
for
(
i
=
0
;
i
<
this
->
size
();
i
++
)
{
(
*
this
)[
i
]
*=
a
;
}
}
inline
void
operator
/=
(
cvm
::
real
const
&
a
)
{
size_t
i
;
for
(
i
=
0
;
i
<
this
->
size
();
i
++
)
{
(
*
this
)[
i
]
/=
a
;
}
}
inline
friend
vector1d
<
T
>
operator
+
(
vector1d
<
T
>
const
&
v1
,
vector1d
<
T
>
const
&
v2
)
{
check_sizes
(
v1
.
size
(),
v2
.
size
());
vector1d
<
T
>
result
(
v1
.
size
());
size_t
i
;
for
(
i
=
0
;
i
<
v1
.
size
();
i
++
)
{
result
[
i
]
=
v1
[
i
]
+
v2
[
i
];
}
return
result
;
}
inline
friend
vector1d
<
T
>
operator
-
(
vector1d
<
T
>
const
&
v1
,
vector1d
<
T
>
const
&
v2
)
{
check_sizes
(
v1
.
size
(),
v2
.
size
());
vector1d
<
T
>
result
(
v1
.
size
());
size_t
i
;
for
(
i
=
0
;
i
<
v1
.
size
();
i
++
)
{
result
[
i
]
=
v1
[
i
]
-
v2
[
i
];
}
return
result
;
}
inline
friend
vector1d
<
T
>
operator
*
(
vector1d
<
T
>
const
&
v
,
cvm
::
real
const
&
a
)
{
vector1d
<
T
>
result
(
v
.
size
());
size_t
i
;
for
(
i
=
0
;
i
<
v
.
size
();
i
++
)
{
result
[
i
]
=
v
[
i
]
*
a
;
}
return
result
;
}
inline
friend
vector1d
<
T
>
operator
*
(
cvm
::
real
const
&
a
,
vector1d
<
T
>
const
&
v
)
{
return
v
*
a
;
}
inline
friend
vector1d
<
T
>
operator
/
(
vector1d
<
T
>
const
&
v
,
cvm
::
real
const
&
a
)
{
vector1d
<
T
>
result
(
v
.
size
());
size_t
i
;
for
(
i
=
0
;
i
<
v
.
size
();
i
++
)
{
result
[
i
]
=
v
[
i
]
/
a
;
}
return
result
;
}
/// Inner product
inline
friend
T
operator
*
(
vector1d
<
T
>
const
&
v1
,
vector1d
<
T
>
const
&
v2
)
{
check_sizes
(
v1
.
size
(),
v2
.
size
());
T
prod
(
0.0
);
size_t
i
;
for
(
i
=
0
;
i
<
v1
.
size
();
i
++
)
{
prod
+=
v1
[
i
]
*
v2
[
i
];
}
return
prod
;
}
/// Squared norm
inline
cvm
::
real
norm2
()
const
{
cvm
::
real
result
=
0.0
;
size_t
i
;
for
(
i
=
0
;
i
<
this
->
size
();
i
++
)
{
result
+=
(
*
this
)[
i
]
*
(
*
this
)[
i
];
}
return
result
;
}
inline
cvm
::
real
norm
()
const
{
return
std
::
sqrt
(
this
->
norm2
());
}
/// Slicing
inline
vector1d
<
T
>
const
slice
(
size_t
const
i1
,
size_t
const
i2
)
const
{
if
((
i2
<
i1
)
||
(
i2
>=
this
->
size
()))
{
cvm
::
error
(
"Error: trying to slice a vector using incorrect boundaries.
\n
"
);
}
vector1d
<
T
>
result
(
i2
-
i1
);
size_t
i
;
for
(
i
=
0
;
i
<
(
i2
-
i1
);
i
++
)
{
result
[
i
]
=
(
*
this
)[
i1
+
i
];
}
return
result
;
}
/// Assign a vector to a slice of this vector
inline
void
sliceassign
(
size_t
const
i1
,
size_t
const
i2
,
vector1d
<
T
>
const
&
v
)
{
if
((
i2
<
i1
)
||
(
i2
>=
this
->
size
()))
{
cvm
::
error
(
"Error: trying to slice a vector using incorrect boundaries.
\n
"
);
}
size_t
i
;
for
(
i
=
0
;
i
<
(
i2
-
i1
);
i
++
)
{
(
*
this
)[
i1
+
i
]
=
v
[
i
];
}
}
/// Formatted output
inline
size_t
output_width
(
size_t
const
&
real_width
)
const
{
return
real_width
*
(
this
->
size
())
+
3
*
(
this
->
size
()
-
1
)
+
4
;
}
inline
friend
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
cvm
::
vector1d
<
T
>
&
v
)
{
if
(
v
.
size
()
==
0
)
return
is
;
size_t
const
start_pos
=
is
.
tellg
();
char
sep
;
if
(
!
(
is
>>
sep
)
||
!
(
sep
==
'('
)
)
{
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
return
is
;
}
size_t
count
=
0
;
while
(
(
is
>>
v
[
count
])
&&
(
count
<
(
v
.
size
()
-
1
)
?
((
is
>>
sep
)
&&
(
sep
==
','
))
:
true
)
)
{
if
(
++
count
==
v
.
size
())
break
;
}
if
(
count
<
v
.
size
())
{
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
}
return
is
;
}
inline
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
cvm
::
vector1d
<
T
>
const
&
v
)
{
std
::
streamsize
const
w
=
os
.
width
();
std
::
streamsize
const
p
=
os
.
precision
();
os
.
width
(
2
);
os
<<
"( "
;
size_t
i
;
for
(
i
=
0
;
i
<
v
.
size
()
-
1
;
i
++
)
{
os
.
width
(
w
);
os
.
precision
(
p
);
os
<<
v
[
i
]
<<
" , "
;
}
os
.
width
(
w
);
os
.
precision
(
p
);
os
<<
v
[
v
.
size
()
-
1
]
<<
" )"
;
return
os
;
}
inline
std
::
string
to_simple_string
()
const
{
if
(
this
->
size
()
==
0
)
return
std
::
string
(
""
);
std
::
ostringstream
os
;
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
os
.
precision
(
cvm
::
cv_prec
);
os
<<
(
*
this
)[
0
];
size_t
i
;
for
(
i
=
1
;
i
<
this
->
size
();
i
++
)
{
os
<<
" "
<<
(
*
this
)[
i
];
}
return
os
.
str
();
}
inline
int
from_simple_string
(
std
::
string
const
&
s
)
{
std
::
stringstream
stream
(
s
);
size_t
i
=
0
;
while
((
stream
>>
(
*
this
)[
i
])
&&
(
i
<
this
->
size
()))
{
i
++
;
}
if
(
i
<
this
->
size
())
{
return
COLVARS_ERROR
;
}
return
COLVARS_OK
;
}
};
/// \brief Arbitrary size array (two dimensions) suitable for linear
/// algebra operations (i.e. for floating point numbers it can be used
/// with library functions)
template
<
class
T
>
class
colvarmodule
::
matrix2d
{
public:
friend
class
row
;
size_t
outer_length
;
size_t
inner_length
;
protected:
class
row
{
public:
T
*
data
;
size_t
length
;
inline
row
(
T
*
const
row_data
,
size_t
const
inner_length
)
:
data
(
row_data
),
length
(
inner_length
)
{}
inline
T
&
operator
[]
(
size_t
const
j
)
{
return
*
(
data
+
j
);
}
inline
T
const
&
operator
[]
(
size_t
const
j
)
const
{
return
*
(
data
+
j
);
}
inline
operator
vector1d
<
T
>
()
const
{
return
vector1d
<
T
>
(
length
,
data
);
}
};
std
::
vector
<
T
>
data
;
std
::
vector
<
row
>
rows
;
std
::
vector
<
T
*>
pointers
;
public:
/// Allocation routine, used by all constructors
inline
void
resize
(
size_t
const
ol
,
size_t
const
il
)
{
if
((
ol
>
0
)
&&
(
il
>
0
))
{
if
(
data
.
size
()
>
0
)
{
// copy previous data
size_t
i
,
j
;
std
::
vector
<
T
>
new_data
(
ol
*
il
);
for
(
i
=
0
;
i
<
outer_length
;
i
++
)
{
for
(
j
=
0
;
j
<
inner_length
;
j
++
)
{
new_data
[
il
*
i
+
j
]
=
data
[
inner_length
*
i
+
j
];
}
}
data
.
resize
(
ol
*
il
);
// copy them back
data
=
new_data
;
}
else
{
data
.
resize
(
ol
*
il
);
}
outer_length
=
ol
;
inner_length
=
il
;
if
(
data
.
size
()
>
0
)
{
// rebuild rows
size_t
i
;
rows
.
clear
();
rows
.
reserve
(
outer_length
);
pointers
.
clear
();
pointers
.
reserve
(
outer_length
);
for
(
i
=
0
;
i
<
outer_length
;
i
++
)
{
rows
.
push_back
(
row
(
&
(
data
[
0
])
+
inner_length
*
i
,
inner_length
));
pointers
.
push_back
(
&
(
data
[
0
])
+
inner_length
*
i
);
}
}
}
else
{
// zero size
data
.
clear
();
rows
.
clear
();
}
}
/// Deallocation routine
inline
void
clear
()
{
rows
.
clear
();
data
.
clear
();
}
/// Set all elements to zero
inline
void
reset
()
{
data
.
assign
(
data
.
size
(),
T
(
0.0
));
}
inline
size_t
size
()
const
{
return
data
.
size
();
}
/// Default constructor
inline
matrix2d
()
:
outer_length
(
0
),
inner_length
(
0
)
{
this
->
resize
(
0
,
0
);
}
inline
matrix2d
(
size_t
const
ol
,
size_t
const
il
)
:
outer_length
(
ol
),
inner_length
(
il
)
{
this
->
resize
(
outer_length
,
inner_length
);
reset
();
}
/// Copy constructor
inline
matrix2d
(
matrix2d
<
T
>
const
&
m
)
:
outer_length
(
m
.
outer_length
),
inner_length
(
m
.
inner_length
)
{
// reinitialize data and rows arrays
this
->
resize
(
outer_length
,
inner_length
);
// copy data
data
=
m
.
data
;
}
/// Destructor
inline
~
matrix2d
()
{
this
->
clear
();
}
inline
row
&
operator
[]
(
size_t
const
i
)
{
return
rows
[
i
];
}
inline
row
const
&
operator
[]
(
size_t
const
i
)
const
{
return
rows
[
i
];
}
/// Assignment
inline
matrix2d
<
T
>
&
operator
=
(
matrix2d
<
T
>
const
&
m
)
{
if
((
outer_length
!=
m
.
outer_length
)
||
(
inner_length
!=
m
.
inner_length
)){
this
->
clear
();
outer_length
=
m
.
outer_length
;
inner_length
=
m
.
inner_length
;
this
->
resize
(
outer_length
,
inner_length
);
}
data
=
m
.
data
;
return
*
this
;
}
/// Return the 2-d C array
inline
T
**
c_array
()
{
if
(
rows
.
size
()
>
0
)
{
return
&
(
pointers
[
0
]);
}
else
{
return
NULL
;
}
}
inline
static
void
check_sizes
(
matrix2d
<
T
>
const
&
m1
,
matrix2d
<
T
>
const
&
m2
)
{
if
((
m1
.
outer_length
!=
m2
.
outer_length
)
||
(
m1
.
inner_length
!=
m2
.
inner_length
))
{
cvm
::
error
(
"Error: trying to perform an operation between matrices of different sizes, "
+
cvm
::
to_str
(
m1
.
outer_length
)
+
"x"
+
cvm
::
to_str
(
m1
.
inner_length
)
+
" and "
+
cvm
::
to_str
(
m2
.
outer_length
)
+
"x"
+
cvm
::
to_str
(
m2
.
inner_length
)
+
".
\n
"
);
}
}
inline
void
operator
+=
(
matrix2d
<
T
>
const
&
m
)
{
check_sizes
(
*
this
,
m
);
size_t
i
;
for
(
i
=
0
;
i
<
data
.
size
();
i
++
)
{
data
[
i
]
+=
m
.
data
[
i
];
}
}
inline
void
operator
-=
(
matrix2d
<
T
>
const
&
m
)
{
check_sizes
(
*
this
,
m
);
size_t
i
;
for
(
i
=
0
;
i
<
data
.
size
();
i
++
)
{
data
[
i
]
-=
m
.
data
[
i
];
}
}
inline
void
operator
*=
(
cvm
::
real
const
&
a
)
{
size_t
i
;
for
(
i
=
0
;
i
<
data
.
size
();
i
++
)
{
data
[
i
]
*=
a
;
}
}
inline
void
operator
/=
(
cvm
::
real
const
&
a
)
{
size_t
i
;
for
(
i
=
0
;
i
<
data
.
size
();
i
++
)
{
data
[
i
]
/=
a
;
}
}
inline
friend
matrix2d
<
T
>
operator
+
(
matrix2d
<
T
>
const
&
m1
,
matrix2d
<
T
>
const
&
m2
)
{
check_sizes
(
m1
,
m2
);
matrix2d
<
T
>
result
(
m1
.
outer_length
,
m1
.
inner_length
);
size_t
i
;
for
(
i
=
0
;
i
<
m1
.
data
.
size
();
i
++
)
{
result
.
data
[
i
]
=
m1
.
data
[
i
]
+
m2
.
data
[
i
];
}
return
result
;
}
inline
friend
matrix2d
<
T
>
operator
-
(
matrix2d
<
T
>
const
&
m1
,
matrix2d
<
T
>
const
&
m2
)
{
check_sizes
(
m1
,
m2
);
matrix2d
<
T
>
result
(
m1
.
outer_length
,
m1
.
inner_length
);
size_t
i
;
for
(
i
=
0
;
i
<
m1
.
data
.
size
();
i
++
)
{
result
.
data
[
i
]
=
m1
.
data
[
i
]
-
m1
.
data
[
i
];
}
return
result
;
}
inline
friend
matrix2d
<
T
>
operator
*
(
matrix2d
<
T
>
const
&
m
,
cvm
::
real
const
&
a
)
{
matrix2d
<
T
>
result
(
m
.
outer_length
,
m
.
inner_length
);
size_t
i
;
for
(
i
=
0
;
i
<
m
.
data
.
size
();
i
++
)
{
result
.
data
[
i
]
=
m
.
data
[
i
]
*
a
;
}
return
result
;
}
inline
friend
matrix2d
<
T
>
operator
*
(
cvm
::
real
const
&
a
,
matrix2d
<
T
>
const
&
m
)
{
return
m
*
a
;
}
inline
friend
matrix2d
<
T
>
operator
/
(
matrix2d
<
T
>
const
&
m
,
cvm
::
real
const
&
a
)
{
matrix2d
<
T
>
result
(
m
.
outer_length
,
m
.
inner_length
);
size_t
i
;
for
(
i
=
0
;
i
<
m
.
data
.
size
();
i
++
)
{
result
.
data
[
i
]
=
m
.
data
[
i
]
*
a
;
}
return
result
;
}
/// Matrix multiplication
// inline friend matrix2d<T> const & operator * (matrix2d<T> const &m1, matrix2d<T> const &m2)
// {
// matrix2d<T> result(m1.outer_length, m2.inner_length);
// if (m1.inner_length != m2.outer_length) {
// cvm::error("Error: trying to multiply two matrices of incompatible sizes, "+
// cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
// cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
// } else {
// size_t i, j, k;
// for (i = 0; i < m1.outer_length; i++) {
// for (j = 0; j < m2.inner_length; j++) {
// for (k = 0; k < m1.inner_length; k++) {
// result[i][j] += m1[i][k] * m2[k][j];
// }
// }
// }
// }
// return result;
// }
/// vector-matrix multiplication
inline
friend
vector1d
<
T
>
operator
*
(
vector1d
<
T
>
const
&
v
,
matrix2d
<
T
>
const
&
m
)
{
vector1d
<
T
>
result
(
m
.
inner_length
);
if
(
m
.
outer_length
!=
v
.
size
())
{
cvm
::
error
(
"Error: trying to multiply a vector and a matrix of incompatible sizes, "
+
cvm
::
to_str
(
v
.
size
())
+
" and "
+
cvm
::
to_str
(
m
.
outer_length
)
+
"x"
+
cvm
::
to_str
(
m
.
inner_length
)
+
".
\n
"
);
}
else
{
size_t
i
,
k
;
for
(
i
=
0
;
i
<
m
.
inner_length
;
i
++
)
{
for
(
k
=
0
;
k
<
m
.
outer_length
;
k
++
)
{
result
[
i
]
+=
m
[
k
][
i
]
*
v
[
k
];
}
}
}
return
result
;
}
// /// matrix-vector multiplication (unused for now)
// inline friend vector1d<T> const & operator * (matrix2d<T> const &m, vector1d<T> const &v)
// {
// vector1d<T> result(m.outer_length);
// if (m.inner_length != v.size()) {
// cvm::error("Error: trying to multiply a matrix and a vector of incompatible sizes, "+
// cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length)
// + " and " + cvm::to_str(v.length) + ".\n");
// } else {
// size_t i, k;
// for (i = 0; i < m.outer_length; i++) {
// for (k = 0; k < m.inner_length; k++) {
// result[i] += m[i][k] * v[k];
// }
// }
// }
// return result;
// }
/// Formatted output
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
matrix2d
<
T
>
const
&
m
)
{
std
::
streamsize
const
w
=
os
.
width
();
std
::
streamsize
const
p
=
os
.
precision
();
os
.
width
(
2
);
os
<<
"( "
;
size_t
i
;
for
(
i
=
0
;
i
<
m
.
outer_length
;
i
++
)
{
os
<<
" ( "
;
size_t
j
;
for
(
j
=
0
;
j
<
m
.
inner_length
-
1
;
j
++
)
{
os
.
width
(
w
);
os
.
precision
(
p
);
os
<<
m
[
i
][
j
]
<<
" , "
;
}
os
.
width
(
w
);
os
.
precision
(
p
);
os
<<
m
[
i
][
m
.
inner_length
-
1
]
<<
" )"
;
}
os
<<
" )"
;
return
os
;
}
inline
std
::
string
to_simple_string
()
const
{
if
(
this
->
size
()
==
0
)
return
std
::
string
(
""
);
std
::
ostringstream
os
;
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
os
.
precision
(
cvm
::
cv_prec
);
os
<<
(
*
this
)[
0
];
size_t
i
;
for
(
i
=
1
;
i
<
data
.
size
();
i
++
)
{
os
<<
" "
<<
data
[
i
];
}
return
os
.
str
();
}
inline
int
from_simple_string
(
std
::
string
const
&
s
)
{
std
::
stringstream
stream
(
s
);
size_t
i
=
0
;
while
((
stream
>>
data
[
i
])
&&
(
i
<
data
.
size
()))
{
i
++
;
}
if
(
i
<
data
.
size
())
{
return
COLVARS_ERROR
;
}
return
COLVARS_OK
;
}
};
/// vector of real numbers with three components
class
colvarmodule
::
rvector
{
public:
cvm
::
real
x
,
y
,
z
;
inline
rvector
()
:
x
(
0.0
),
y
(
0.0
),
z
(
0.0
)
{}
inline
rvector
(
cvm
::
real
const
&
x_i
,
cvm
::
real
const
&
y_i
,
cvm
::
real
const
&
z_i
)
:
x
(
x_i
),
y
(
y_i
),
z
(
z_i
)
{}
inline
rvector
(
cvm
::
vector1d
<
cvm
::
real
>
const
&
v
)
:
x
(
v
[
0
]),
y
(
v
[
1
]),
z
(
v
[
2
])
{}
inline
rvector
(
cvm
::
real
t
)
:
x
(
t
),
y
(
t
),
z
(
t
)
{}
/// \brief Set all components to a scalar value
inline
void
set
(
cvm
::
real
const
&
value
)
{
x
=
y
=
z
=
value
;
}
/// \brief Assign all components
inline
void
set
(
cvm
::
real
const
&
x_i
,
cvm
::
real
const
&
y_i
,
cvm
::
real
const
&
z_i
)
{
x
=
x_i
;
y
=
y_i
;
z
=
z_i
;
}
/// \brief Set all components to zero
inline
void
reset
()
{
x
=
y
=
z
=
0.0
;
}
/// \brief Access cartesian components by index
inline
cvm
::
real
&
operator
[]
(
int
const
&
i
)
{
return
(
i
==
0
)
?
x
:
(
i
==
1
)
?
y
:
(
i
==
2
)
?
z
:
x
;
}
/// \brief Access cartesian components by index
inline
cvm
::
real
const
&
operator
[]
(
int
const
&
i
)
const
{
return
(
i
==
0
)
?
x
:
(
i
==
1
)
?
y
:
(
i
==
2
)
?
z
:
x
;
}
inline
cvm
::
vector1d
<
cvm
::
real
>
const
as_vector
()
const
{
cvm
::
vector1d
<
cvm
::
real
>
result
(
3
);
result
[
0
]
=
this
->
x
;
result
[
1
]
=
this
->
y
;
result
[
2
]
=
this
->
z
;
return
result
;
}
inline
cvm
::
rvector
&
operator
=
(
cvm
::
real
const
&
v
)
{
x
=
v
;
y
=
v
;
z
=
v
;
return
*
this
;
}
inline
void
operator
+=
(
cvm
::
rvector
const
&
v
)
{
x
+=
v
.
x
;
y
+=
v
.
y
;
z
+=
v
.
z
;
}
inline
void
operator
-=
(
cvm
::
rvector
const
&
v
)
{
x
-=
v
.
x
;
y
-=
v
.
y
;
z
-=
v
.
z
;
}
inline
void
operator
*=
(
cvm
::
real
const
&
v
)
{
x
*=
v
;
y
*=
v
;
z
*=
v
;
}
inline
void
operator
/=
(
cvm
::
real
const
&
v
)
{
x
/=
v
;
y
/=
v
;
z
/=
v
;
}
inline
cvm
::
real
norm2
()
const
{
return
(
x
*
x
+
y
*
y
+
z
*
z
);
}
inline
cvm
::
real
norm
()
const
{
return
std
::
sqrt
(
this
->
norm2
());
}
inline
cvm
::
rvector
unit
()
const
{
const
cvm
::
real
n
=
this
->
norm
();
return
(
n
>
0.
?
cvm
::
rvector
(
x
,
y
,
z
)
/
n
:
cvm
::
rvector
(
1.
,
0.
,
0.
));
}
static
inline
size_t
output_width
(
size_t
const
&
real_width
)
{
return
3
*
real_width
+
10
;
}
static
inline
cvm
::
rvector
outer
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
cvm
::
rvector
(
v1
.
y
*
v2
.
z
-
v2
.
y
*
v1
.
z
,
-
v1
.
x
*
v2
.
z
+
v2
.
x
*
v1
.
z
,
v1
.
x
*
v2
.
y
-
v2
.
x
*
v1
.
y
);
}
friend
inline
cvm
::
rvector
operator
-
(
cvm
::
rvector
const
&
v
)
{
return
cvm
::
rvector
(
-
v
.
x
,
-
v
.
y
,
-
v
.
z
);
}
friend
inline
int
operator
==
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
(
v1
.
x
==
v2
.
x
)
&&
(
v1
.
y
==
v2
.
y
)
&&
(
v1
.
z
==
v2
.
z
);
}
friend
inline
int
operator
!=
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
(
v1
.
x
!=
v2
.
x
)
||
(
v1
.
y
!=
v2
.
y
)
||
(
v1
.
z
!=
v2
.
z
);
}
friend
inline
cvm
::
rvector
operator
+
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
cvm
::
rvector
(
v1
.
x
+
v2
.
x
,
v1
.
y
+
v2
.
y
,
v1
.
z
+
v2
.
z
);
}
friend
inline
cvm
::
rvector
operator
-
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
cvm
::
rvector
(
v1
.
x
-
v2
.
x
,
v1
.
y
-
v2
.
y
,
v1
.
z
-
v2
.
z
);
}
friend
inline
cvm
::
real
operator
*
(
cvm
::
rvector
const
&
v1
,
cvm
::
rvector
const
&
v2
)
{
return
v1
.
x
*
v2
.
x
+
v1
.
y
*
v2
.
y
+
v1
.
z
*
v2
.
z
;
}
friend
inline
cvm
::
rvector
operator
*
(
cvm
::
real
const
&
a
,
cvm
::
rvector
const
&
v
)
{
return
cvm
::
rvector
(
a
*
v
.
x
,
a
*
v
.
y
,
a
*
v
.
z
);
}
friend
inline
cvm
::
rvector
operator
*
(
cvm
::
rvector
const
&
v
,
cvm
::
real
const
&
a
)
{
return
cvm
::
rvector
(
a
*
v
.
x
,
a
*
v
.
y
,
a
*
v
.
z
);
}
friend
inline
cvm
::
rvector
operator
/
(
cvm
::
rvector
const
&
v
,
cvm
::
real
const
&
a
)
{
return
cvm
::
rvector
(
v
.
x
/
a
,
v
.
y
/
a
,
v
.
z
/
a
);
}
std
::
string
to_simple_string
()
const
;
int
from_simple_string
(
std
::
string
const
&
s
);
};
/// \brief 2-dimensional array of real numbers with three components
/// along each dimension (works with colvarmodule::rvector)
class
colvarmodule
::
rmatrix
:
public
colvarmodule
::
matrix2d
<
colvarmodule
::
real
>
{
private:
public:
/// Return the xx element
inline
cvm
::
real
&
xx
()
{
return
(
*
this
)[
0
][
0
];
}
/// Return the xy element
inline
cvm
::
real
&
xy
()
{
return
(
*
this
)[
0
][
1
];
}
/// Return the xz element
inline
cvm
::
real
&
xz
()
{
return
(
*
this
)[
0
][
2
];
}
/// Return the yx element
inline
cvm
::
real
&
yx
()
{
return
(
*
this
)[
1
][
0
];
}
/// Return the yy element
inline
cvm
::
real
&
yy
()
{
return
(
*
this
)[
1
][
1
];
}
/// Return the yz element
inline
cvm
::
real
&
yz
()
{
return
(
*
this
)[
1
][
2
];
}
/// Return the zx element
inline
cvm
::
real
&
zx
()
{
return
(
*
this
)[
2
][
0
];
}
/// Return the zy element
inline
cvm
::
real
&
zy
()
{
return
(
*
this
)[
2
][
1
];
}
/// Return the zz element
inline
cvm
::
real
&
zz
()
{
return
(
*
this
)[
2
][
2
];
}
/// Return the xx element
inline
cvm
::
real
xx
()
const
{
return
(
*
this
)[
0
][
0
];
}
/// Return the xy element
inline
cvm
::
real
xy
()
const
{
return
(
*
this
)[
0
][
1
];
}
/// Return the xz element
inline
cvm
::
real
xz
()
const
{
return
(
*
this
)[
0
][
2
];
}
/// Return the yx element
inline
cvm
::
real
yx
()
const
{
return
(
*
this
)[
1
][
0
];
}
/// Return the yy element
inline
cvm
::
real
yy
()
const
{
return
(
*
this
)[
1
][
1
];
}
/// Return the yz element
inline
cvm
::
real
yz
()
const
{
return
(
*
this
)[
1
][
2
];
}
/// Return the zx element
inline
cvm
::
real
zx
()
const
{
return
(
*
this
)[
2
][
0
];
}
/// Return the zy element
inline
cvm
::
real
zy
()
const
{
return
(
*
this
)[
2
][
1
];
}
/// Return the zz element
inline
cvm
::
real
zz
()
const
{
return
(
*
this
)[
2
][
2
];
}
/// Default constructor
inline
rmatrix
()
:
cvm
::
matrix2d
<
cvm
::
real
>
(
3
,
3
)
{}
/// Constructor component by component
inline
rmatrix
(
cvm
::
real
const
&
xxi
,
cvm
::
real
const
&
xyi
,
cvm
::
real
const
&
xzi
,
cvm
::
real
const
&
yxi
,
cvm
::
real
const
&
yyi
,
cvm
::
real
const
&
yzi
,
cvm
::
real
const
&
zxi
,
cvm
::
real
const
&
zyi
,
cvm
::
real
const
&
zzi
)
:
cvm
::
matrix2d
<
cvm
::
real
>
(
3
,
3
)
{
this
->
xx
()
=
xxi
;
this
->
xy
()
=
xyi
;
this
->
xz
()
=
xzi
;
this
->
yx
()
=
yxi
;
this
->
yy
()
=
yyi
;
this
->
yz
()
=
yzi
;
this
->
zx
()
=
zxi
;
this
->
zy
()
=
zyi
;
this
->
zz
()
=
zzi
;
}
/// Destructor
inline
~
rmatrix
()
{}
/// Return the determinant
inline
cvm
::
real
determinant
()
const
{
return
(
xx
()
*
(
yy
()
*
zz
()
-
zy
()
*
yz
()))
-
(
yx
()
*
(
xy
()
*
zz
()
-
zy
()
*
xz
()))
+
(
zx
()
*
(
xy
()
*
yz
()
-
yy
()
*
xz
()));
}
inline
cvm
::
rmatrix
transpose
()
const
{
return
cvm
::
rmatrix
(
this
->
xx
(),
this
->
yx
(),
this
->
zx
(),
this
->
xy
(),
this
->
yy
(),
this
->
zy
(),
this
->
xz
(),
this
->
yz
(),
this
->
zz
());
}
friend
cvm
::
rvector
operator
*
(
cvm
::
rmatrix
const
&
m
,
cvm
::
rvector
const
&
r
);
// friend inline cvm::rmatrix const operator * (cvm::rmatrix const &m1, cvm::rmatrix const &m2) {
// return cvm::rmatrix (m1.xx()*m2.xx() + m1.xy()*m2.yx() + m1.xz()*m2.yz(),
// m1.xx()*m2.xy() + m1.xy()*m2.yy() + m1.xz()*m2.zy(),
// m1.xx()*m2.xz() + m1.xy()*m2.yz() + m1.xz()*m2.zz(),
// m1.yx()*m2.xx() + m1.yy()*m2.yx() + m1.yz()*m2.yz(),
// m1.yx()*m2.xy() + m1.yy()*m2.yy() + m1.yz()*m2.yy(),
// m1.yx()*m2.xz() + m1.yy()*m2.yz() + m1.yz()*m2.yz(),
// m1.zx()*m2.xx() + m1.zy()*m2.yx() + m1.zz()*m2.yz(),
// m1.zx()*m2.xy() + m1.zy()*m2.yy() + m1.zz()*m2.yy(),
// m1.zx()*m2.xz() + m1.zy()*m2.yz() + m1.zz()*m2.yz());
// }
};
inline
cvm
::
rvector
operator
*
(
cvm
::
rmatrix
const
&
m
,
cvm
::
rvector
const
&
r
)
{
return
cvm
::
rvector
(
m
.
xx
()
*
r
.
x
+
m
.
xy
()
*
r
.
y
+
m
.
xz
()
*
r
.
z
,
m
.
yx
()
*
r
.
x
+
m
.
yy
()
*
r
.
y
+
m
.
yz
()
*
r
.
z
,
m
.
zx
()
*
r
.
x
+
m
.
zy
()
*
r
.
y
+
m
.
zz
()
*
r
.
z
);
}
/// Numerical recipes diagonalization
void
jacobi
(
cvm
::
real
**
a
,
cvm
::
real
*
d
,
cvm
::
real
**
v
,
int
*
nrot
);
/// Eigenvector sort
void
eigsrt
(
cvm
::
real
*
d
,
cvm
::
real
**
v
);
/// Transpose the matrix
void
transpose
(
cvm
::
real
**
v
);
/// \brief 1-dimensional vector of real numbers with four components and
/// a quaternion algebra
class
colvarmodule
::
quaternion
{
public:
cvm
::
real
q0
,
q1
,
q2
,
q3
;
/// Constructor from a 3-d vector
inline
quaternion
(
cvm
::
real
const
&
x
,
cvm
::
real
const
&
y
,
cvm
::
real
const
&
z
)
:
q0
(
0.0
),
q1
(
x
),
q2
(
y
),
q3
(
z
)
{}
/// Constructor component by component
inline
quaternion
(
cvm
::
real
const
qv
[
4
])
:
q0
(
qv
[
0
]),
q1
(
qv
[
1
]),
q2
(
qv
[
2
]),
q3
(
qv
[
3
])
{}
/// Constructor component by component
inline
quaternion
(
cvm
::
real
const
&
q0i
,
cvm
::
real
const
&
q1i
,
cvm
::
real
const
&
q2i
,
cvm
::
real
const
&
q3i
)
:
q0
(
q0i
),
q1
(
q1i
),
q2
(
q2i
),
q3
(
q3i
)
{}
inline
quaternion
(
cvm
::
vector1d
<
cvm
::
real
>
const
&
v
)
:
q0
(
v
[
0
]),
q1
(
v
[
1
]),
q2
(
v
[
2
]),
q3
(
v
[
3
])
{}
/// "Constructor" after Euler angles (in radians)
///
/// http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
inline
void
set_from_euler_angles
(
cvm
::
real
const
&
phi_in
,
cvm
::
real
const
&
theta_in
,
cvm
::
real
const
&
psi_in
)
{
q0
=
(
(
std
::
cos
(
phi_in
/
2.0
))
*
(
std
::
cos
(
theta_in
/
2.0
))
*
(
std
::
cos
(
psi_in
/
2.0
))
+
(
std
::
sin
(
phi_in
/
2.0
))
*
(
std
::
sin
(
theta_in
/
2.0
))
*
(
std
::
sin
(
psi_in
/
2.0
))
);
q1
=
(
(
std
::
sin
(
phi_in
/
2.0
))
*
(
std
::
cos
(
theta_in
/
2.0
))
*
(
std
::
cos
(
psi_in
/
2.0
))
-
(
std
::
cos
(
phi_in
/
2.0
))
*
(
std
::
sin
(
theta_in
/
2.0
))
*
(
std
::
sin
(
psi_in
/
2.0
))
);
q2
=
(
(
std
::
cos
(
phi_in
/
2.0
))
*
(
std
::
sin
(
theta_in
/
2.0
))
*
(
std
::
cos
(
psi_in
/
2.0
))
+
(
std
::
sin
(
phi_in
/
2.0
))
*
(
std
::
cos
(
theta_in
/
2.0
))
*
(
std
::
sin
(
psi_in
/
2.0
))
);
q3
=
(
(
std
::
cos
(
phi_in
/
2.0
))
*
(
std
::
cos
(
theta_in
/
2.0
))
*
(
std
::
sin
(
psi_in
/
2.0
))
-
(
std
::
sin
(
phi_in
/
2.0
))
*
(
std
::
sin
(
theta_in
/
2.0
))
*
(
std
::
cos
(
psi_in
/
2.0
))
);
}
/// \brief Default constructor
inline
quaternion
()
{
reset
();
}
/// \brief Set all components to a scalar
inline
void
set
(
cvm
::
real
const
&
value
=
0.0
)
{
q0
=
q1
=
q2
=
q3
=
value
;
}
/// \brief Set all components to zero (null quaternion)
inline
void
reset
()
{
q0
=
q1
=
q2
=
q3
=
0.0
;
}
/// \brief Set the q0 component to 1 and the others to 0 (quaternion
/// representing no rotation)
inline
void
reset_rotation
()
{
q0
=
1.0
;
q1
=
q2
=
q3
=
0.0
;
}
/// Tell the number of characters required to print a quaternion, given that of a real number
static
inline
size_t
output_width
(
size_t
const
&
real_width
)
{
return
4
*
real_width
+
13
;
}
std
::
string
to_simple_string
()
const
;
int
from_simple_string
(
std
::
string
const
&
s
);
/// \brief Formatted output operator
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
cvm
::
quaternion
const
&
q
);
/// \brief Formatted input operator
friend
std
::
istream
&
operator
>>
(
std
::
istream
&
is
,
cvm
::
quaternion
&
q
);
/// Access the quaternion as a 4-d array (return a reference)
inline
cvm
::
real
&
operator
[]
(
int
const
&
i
)
{
switch
(
i
)
{
case
0
:
return
this
->
q0
;
case
1
:
return
this
->
q1
;
case
2
:
return
this
->
q2
;
case
3
:
return
this
->
q3
;
default:
cvm
::
error
(
"Error: incorrect quaternion component.
\n
"
);
return
q0
;
}
}
/// Access the quaternion as a 4-d array (return a value)
inline
cvm
::
real
operator
[]
(
int
const
&
i
)
const
{
switch
(
i
)
{
case
0
:
return
this
->
q0
;
case
1
:
return
this
->
q1
;
case
2
:
return
this
->
q2
;
case
3
:
return
this
->
q3
;
default:
cvm
::
error
(
"Error: trying to access a quaternion "
"component which is not between 0 and 3.
\n
"
);
return
0.0
;
}
}
inline
cvm
::
vector1d
<
cvm
::
real
>
const
as_vector
()
const
{
cvm
::
vector1d
<
cvm
::
real
>
result
(
4
);
result
[
0
]
=
q0
;
result
[
1
]
=
q1
;
result
[
2
]
=
q2
;
result
[
3
]
=
q3
;
return
result
;
}
/// Square norm of the quaternion
inline
cvm
::
real
norm2
()
const
{
return
q0
*
q0
+
q1
*
q1
+
q2
*
q2
+
q3
*
q3
;
}
/// Norm of the quaternion
inline
cvm
::
real
norm
()
const
{
return
std
::
sqrt
(
this
->
norm2
());
}
/// Return the conjugate quaternion
inline
cvm
::
quaternion
conjugate
()
const
{
return
cvm
::
quaternion
(
q0
,
-
q1
,
-
q2
,
-
q3
);
}
inline
void
operator
*=
(
cvm
::
real
const
&
a
)
{
q0
*=
a
;
q1
*=
a
;
q2
*=
a
;
q3
*=
a
;
}
inline
void
operator
/=
(
cvm
::
real
const
&
a
)
{
q0
/=
a
;
q1
/=
a
;
q2
/=
a
;
q3
/=
a
;
}
inline
void
set_positive
()
{
if
(
q0
>
0.0
)
return
;
q0
=
-
q0
;
q1
=
-
q1
;
q2
=
-
q2
;
q3
=
-
q3
;
}
inline
void
operator
+=
(
cvm
::
quaternion
const
&
h
)
{
q0
+=
h
.
q0
;
q1
+=
h
.
q1
;
q2
+=
h
.
q2
;
q3
+=
h
.
q3
;
}
inline
void
operator
-=
(
cvm
::
quaternion
const
&
h
)
{
q0
-=
h
.
q0
;
q1
-=
h
.
q1
;
q2
-=
h
.
q2
;
q3
-=
h
.
q3
;
}
/// Promote a 3-vector to a quaternion
static
inline
cvm
::
quaternion
promote
(
cvm
::
rvector
const
&
v
)
{
return
cvm
::
quaternion
(
0.0
,
v
.
x
,
v
.
y
,
v
.
z
);
}
/// Return the vector component
inline
cvm
::
rvector
get_vector
()
const
{
return
cvm
::
rvector
(
q1
,
q2
,
q3
);
}
friend
inline
cvm
::
quaternion
operator
+
(
cvm
::
quaternion
const
&
h
,
cvm
::
quaternion
const
&
q
)
{
return
cvm
::
quaternion
(
h
.
q0
+
q
.
q0
,
h
.
q1
+
q
.
q1
,
h
.
q2
+
q
.
q2
,
h
.
q3
+
q
.
q3
);
}
friend
inline
cvm
::
quaternion
operator
-
(
cvm
::
quaternion
const
&
h
,
cvm
::
quaternion
const
&
q
)
{
return
cvm
::
quaternion
(
h
.
q0
-
q
.
q0
,
h
.
q1
-
q
.
q1
,
h
.
q2
-
q
.
q2
,
h
.
q3
-
q
.
q3
);
}
/// \brief Provides the quaternion product. \b NOTE: for the inner
/// product use: \code h.inner (q); \endcode
friend
inline
cvm
::
quaternion
operator
*
(
cvm
::
quaternion
const
&
h
,
cvm
::
quaternion
const
&
q
)
{
return
cvm
::
quaternion
(
h
.
q0
*
q
.
q0
-
h
.
q1
*
q
.
q1
-
h
.
q2
*
q
.
q2
-
h
.
q3
*
q
.
q3
,
h
.
q0
*
q
.
q1
+
h
.
q1
*
q
.
q0
+
h
.
q2
*
q
.
q3
-
h
.
q3
*
q
.
q2
,
h
.
q0
*
q
.
q2
+
h
.
q2
*
q
.
q0
+
h
.
q3
*
q
.
q1
-
h
.
q1
*
q
.
q3
,
h
.
q0
*
q
.
q3
+
h
.
q3
*
q
.
q0
+
h
.
q1
*
q
.
q2
-
h
.
q2
*
q
.
q1
);
}
friend
inline
cvm
::
quaternion
operator
*
(
cvm
::
real
const
&
c
,
cvm
::
quaternion
const
&
q
)
{
return
cvm
::
quaternion
(
c
*
q
.
q0
,
c
*
q
.
q1
,
c
*
q
.
q2
,
c
*
q
.
q3
);
}
friend
inline
cvm
::
quaternion
operator
*
(
cvm
::
quaternion
const
&
q
,
cvm
::
real
const
&
c
)
{
return
cvm
::
quaternion
(
q
.
q0
*
c
,
q
.
q1
*
c
,
q
.
q2
*
c
,
q
.
q3
*
c
);
}
friend
inline
cvm
::
quaternion
operator
/
(
cvm
::
quaternion
const
&
q
,
cvm
::
real
const
&
c
)
{
return
cvm
::
quaternion
(
q
.
q0
/
c
,
q
.
q1
/
c
,
q
.
q2
/
c
,
q
.
q3
/
c
);
}
/// \brief Rotate v through this quaternion (put it in the rotated
/// reference frame)
inline
cvm
::
rvector
rotate
(
cvm
::
rvector
const
&
v
)
const
{
return
((
*
this
)
*
promote
(
v
)
*
((
*
this
).
conjugate
())).
get_vector
();
}
/// \brief Rotate Q2 through this quaternion (put it in the rotated
/// reference frame)
inline
cvm
::
quaternion
rotate
(
cvm
::
quaternion
const
&
Q2
)
const
{
cvm
::
rvector
const
vq_rot
=
this
->
rotate
(
Q2
.
get_vector
());
return
cvm
::
quaternion
(
Q2
.
q0
,
vq_rot
.
x
,
vq_rot
.
y
,
vq_rot
.
z
);
}
/// Return the 3x3 matrix associated to this quaternion
inline
cvm
::
rmatrix
rotation_matrix
()
const
{
cvm
::
rmatrix
R
;
R
.
xx
()
=
q0
*
q0
+
q1
*
q1
-
q2
*
q2
-
q3
*
q3
;
R
.
yy
()
=
q0
*
q0
-
q1
*
q1
+
q2
*
q2
-
q3
*
q3
;
R
.
zz
()
=
q0
*
q0
-
q1
*
q1
-
q2
*
q2
+
q3
*
q3
;
R
.
xy
()
=
2.0
*
(
q1
*
q2
-
q0
*
q3
);
R
.
xz
()
=
2.0
*
(
q0
*
q2
+
q1
*
q3
);
R
.
yx
()
=
2.0
*
(
q0
*
q3
+
q1
*
q2
);
R
.
yz
()
=
2.0
*
(
q2
*
q3
-
q0
*
q1
);
R
.
zx
()
=
2.0
*
(
q1
*
q3
-
q0
*
q2
);
R
.
zy
()
=
2.0
*
(
q0
*
q1
+
q2
*
q3
);
return
R
;
}
/// \brief Multiply the given vector by the derivative of the given
/// (rotated) position with respect to the quaternion
cvm
::
quaternion
position_derivative_inner
(
cvm
::
rvector
const
&
pos
,
cvm
::
rvector
const
&
vec
)
const
;
/// \brief Return the cosine between the orientation frame
/// associated to this quaternion and another
inline
cvm
::
real
cosine
(
cvm
::
quaternion
const
&
q
)
const
{
cvm
::
real
const
iprod
=
this
->
inner
(
q
);
return
2.0
*
iprod
*
iprod
-
1.0
;
}
/// \brief Square distance from another quaternion on the
/// 4-dimensional unit sphere: returns the square of the angle along
/// the shorter of the two geodesics
inline
cvm
::
real
dist2
(
cvm
::
quaternion
const
&
Q2
)
const
{
cvm
::
real
const
cos_omega
=
this
->
q0
*
Q2
.
q0
+
this
->
q1
*
Q2
.
q1
+
this
->
q2
*
Q2
.
q2
+
this
->
q3
*
Q2
.
q3
;
cvm
::
real
const
omega
=
std
::
acos
(
(
cos_omega
>
1.0
)
?
1.0
:
(
(
cos_omega
<
-
1.0
)
?
-
1.0
:
cos_omega
)
);
// get the minimum distance: x and -x are the same quaternion
if
(
cos_omega
>
0.0
)
return
omega
*
omega
;
else
return
(
PI
-
omega
)
*
(
PI
-
omega
);
}
/// Gradient of the square distance: returns a 4-vector equivalent
/// to that provided by slerp
inline
cvm
::
quaternion
dist2_grad
(
cvm
::
quaternion
const
&
Q2
)
const
{
cvm
::
real
const
cos_omega
=
this
->
q0
*
Q2
.
q0
+
this
->
q1
*
Q2
.
q1
+
this
->
q2
*
Q2
.
q2
+
this
->
q3
*
Q2
.
q3
;
cvm
::
real
const
omega
=
std
::
acos
(
(
cos_omega
>
1.0
)
?
1.0
:
(
(
cos_omega
<
-
1.0
)
?
-
1.0
:
cos_omega
)
);
cvm
::
real
const
sin_omega
=
std
::
sin
(
omega
);
if
(
std
::
fabs
(
sin_omega
)
<
1.0E-14
)
{
// return a null 4d vector
return
cvm
::
quaternion
(
0.0
,
0.0
,
0.0
,
0.0
);
}
cvm
::
quaternion
const
grad1
((
-
1.0
)
*
sin_omega
*
Q2
.
q0
+
cos_omega
*
(
this
->
q0
-
cos_omega
*
Q2
.
q0
)
/
sin_omega
,
(
-
1.0
)
*
sin_omega
*
Q2
.
q1
+
cos_omega
*
(
this
->
q1
-
cos_omega
*
Q2
.
q1
)
/
sin_omega
,
(
-
1.0
)
*
sin_omega
*
Q2
.
q2
+
cos_omega
*
(
this
->
q2
-
cos_omega
*
Q2
.
q2
)
/
sin_omega
,
(
-
1.0
)
*
sin_omega
*
Q2
.
q3
+
cos_omega
*
(
this
->
q3
-
cos_omega
*
Q2
.
q3
)
/
sin_omega
);
if
(
cos_omega
>
0.0
)
{
return
2.0
*
omega
*
grad1
;
}
else
{
return
-
2.0
*
(
PI
-
omega
)
*
grad1
;
}
}
/// \brief Choose the closest between Q2 and -Q2 and save it back.
/// Not required for dist2() and dist2_grad()
inline
void
match
(
cvm
::
quaternion
&
Q2
)
const
{
cvm
::
real
const
cos_omega
=
this
->
q0
*
Q2
.
q0
+
this
->
q1
*
Q2
.
q1
+
this
->
q2
*
Q2
.
q2
+
this
->
q3
*
Q2
.
q3
;
if
(
cos_omega
<
0.0
)
Q2
*=
-
1.0
;
}
/// \brief Inner product (as a 4-d vector) with Q2; requires match()
/// if the largest overlap is looked for
inline
cvm
::
real
inner
(
cvm
::
quaternion
const
&
Q2
)
const
{
cvm
::
real
const
prod
=
this
->
q0
*
Q2
.
q0
+
this
->
q1
*
Q2
.
q1
+
this
->
q2
*
Q2
.
q2
+
this
->
q3
*
Q2
.
q3
;
return
prod
;
}
};
/// \brief A rotation between two sets of coordinates (for the moment
/// a wrapper for colvarmodule::quaternion)
class
colvarmodule
::
rotation
{
public:
/// \brief The rotation itself (implemented as a quaternion)
cvm
::
quaternion
q
;
/// \brief Eigenvalue corresponding to the optimal rotation
cvm
::
real
lambda
;
/// \brief Perform gradient tests
bool
b_debug_gradients
;
/// \brief Positions to superimpose: the rotation should brings pos1
/// into pos2
std
::
vector
<
cvm
::
atom_pos
>
pos1
,
pos2
;
cvm
::
rmatrix
C
;
cvm
::
matrix2d
<
cvm
::
real
>
S
;
cvm
::
vector1d
<
cvm
::
real
>
S_eigval
;
cvm
::
matrix2d
<
cvm
::
real
>
S_eigvec
;
/// Used for debugging gradients
cvm
::
matrix2d
<
cvm
::
real
>
S_backup
;
/// Derivatives of S
std
::
vector
<
cvm
::
matrix2d
<
cvm
::
rvector
>
>
dS_1
,
dS_2
;
/// Derivatives of leading eigenvalue
std
::
vector
<
cvm
::
rvector
>
dL0_1
,
dL0_2
;
/// Derivatives of leading eigenvector
std
::
vector
<
cvm
::
vector1d
<
cvm
::
rvector
>
>
dQ0_1
,
dQ0_2
;
/// Allocate space for the derivatives of the rotation
inline
void
request_group1_gradients
(
size_t
const
&
n
)
{
dS_1
.
resize
(
n
,
cvm
::
matrix2d
<
cvm
::
rvector
>
(
4
,
4
));
dL0_1
.
resize
(
n
,
cvm
::
rvector
(
0.0
,
0.0
,
0.0
));
dQ0_1
.
resize
(
n
,
cvm
::
vector1d
<
cvm
::
rvector
>
(
4
));
}
/// Allocate space for the derivatives of the rotation
inline
void
request_group2_gradients
(
size_t
const
&
n
)
{
dS_2
.
resize
(
n
,
cvm
::
matrix2d
<
cvm
::
rvector
>
(
4
,
4
));
dL0_2
.
resize
(
n
,
cvm
::
rvector
(
0.0
,
0.0
,
0.0
));
dQ0_2
.
resize
(
n
,
cvm
::
vector1d
<
cvm
::
rvector
>
(
4
));
}
/// \brief Calculate the optimal rotation and store the
/// corresponding eigenvalue and eigenvector in the arguments l0 and
/// q0; if the gradients have been previously requested, calculate
/// them as well
///
/// The method to derive the optimal rotation is defined in:
/// Coutsias EA, Seok C, Dill KA.
/// Using quaternions to calculate RMSD.
/// J Comput Chem. 25(15):1849-57 (2004)
/// DOI: 10.1002/jcc.20110 PubMed: 15376254
void
calc_optimal_rotation
(
std
::
vector
<
atom_pos
>
const
&
pos1
,
std
::
vector
<
atom_pos
>
const
&
pos2
);
/// Default constructor
inline
rotation
()
:
b_debug_gradients
(
false
)
{}
/// Constructor after a quaternion
inline
rotation
(
cvm
::
quaternion
const
&
qi
)
:
q
(
qi
),
b_debug_gradients
(
false
)
{
}
/// Constructor after an axis of rotation and an angle (in radians)
inline
rotation
(
cvm
::
real
const
&
angle
,
cvm
::
rvector
const
&
axis
)
:
b_debug_gradients
(
false
)
{
cvm
::
rvector
const
axis_n
=
axis
.
unit
();
cvm
::
real
const
sina
=
std
::
sin
(
angle
/
2.0
);
q
=
cvm
::
quaternion
(
std
::
cos
(
angle
/
2.0
),
sina
*
axis_n
.
x
,
sina
*
axis_n
.
y
,
sina
*
axis_n
.
z
);
}
/// Destructor
inline
~
rotation
()
{}
/// Return the rotated vector
inline
cvm
::
rvector
rotate
(
cvm
::
rvector
const
&
v
)
const
{
return
q
.
rotate
(
v
);
}
/// Return the inverse of this rotation
inline
cvm
::
rotation
inverse
()
const
{
return
cvm
::
rotation
(
this
->
q
.
conjugate
());
}
/// Return the associated 3x3 matrix
inline
cvm
::
rmatrix
matrix
()
const
{
return
q
.
rotation_matrix
();
}
/// \brief Return the spin angle (in degrees) with respect to the
/// provided axis (which MUST be normalized)
inline
cvm
::
real
spin_angle
(
cvm
::
rvector
const
&
axis
)
const
{
cvm
::
rvector
const
q_vec
=
q
.
get_vector
();
cvm
::
real
alpha
=
(
180.0
/
PI
)
*
2.0
*
std
::
atan2
(
axis
*
q_vec
,
q
.
q0
);
while
(
alpha
>
180.0
)
alpha
-=
360
;
while
(
alpha
<
-
180.0
)
alpha
+=
360
;
return
alpha
;
}
/// \brief Return the derivative of the spin angle with respect to
/// the quaternion
inline
cvm
::
quaternion
dspin_angle_dq
(
cvm
::
rvector
const
&
axis
)
const
{
cvm
::
rvector
const
q_vec
=
q
.
get_vector
();
cvm
::
real
const
iprod
=
axis
*
q_vec
;
if
(
q
.
q0
!=
0.0
)
{
// cvm::real const x = iprod/q.q0;
cvm
::
real
const
dspindx
=
(
180.0
/
PI
)
*
2.0
*
(
1.0
/
(
1.0
+
(
iprod
*
iprod
)
/
(
q
.
q0
*
q
.
q0
)));
return
cvm
::
quaternion
(
dspindx
*
(
iprod
*
(
-
1.0
)
/
(
q
.
q0
*
q
.
q0
)),
dspindx
*
((
1.0
/
q
.
q0
)
*
axis
.
x
),
dspindx
*
((
1.0
/
q
.
q0
)
*
axis
.
y
),
dspindx
*
((
1.0
/
q
.
q0
)
*
axis
.
z
));
}
else
{
// (1/(1+x^2)) ~ (1/x)^2
return
cvm
::
quaternion
((
180.0
/
PI
)
*
2.0
*
((
-
1.0
)
/
iprod
),
0.0
,
0.0
,
0.0
);
// XX TODO: What if iprod == 0? XX
}
}
/// \brief Return the projection of the orientation vector onto a
/// predefined axis
inline
cvm
::
real
cos_theta
(
cvm
::
rvector
const
&
axis
)
const
{
cvm
::
rvector
const
q_vec
=
q
.
get_vector
();
cvm
::
real
const
alpha
=
(
180.0
/
PI
)
*
2.0
*
std
::
atan2
(
axis
*
q_vec
,
q
.
q0
);
cvm
::
real
const
cos_spin_2
=
std
::
cos
(
alpha
*
(
PI
/
180.0
)
*
0.5
);
cvm
::
real
const
cos_theta_2
=
(
(
cos_spin_2
!=
0.0
)
?
(
q
.
q0
/
cos_spin_2
)
:
(
0.0
)
);
// cos(2t) = 2*cos(t)^2 - 1
return
2.0
*
(
cos_theta_2
*
cos_theta_2
)
-
1.0
;
}
/// Return the derivative of the tilt wrt the quaternion
inline
cvm
::
quaternion
dcos_theta_dq
(
cvm
::
rvector
const
&
axis
)
const
{
cvm
::
rvector
const
q_vec
=
q
.
get_vector
();
cvm
::
real
const
iprod
=
axis
*
q_vec
;
cvm
::
real
const
cos_spin_2
=
std
::
cos
(
std
::
atan2
(
iprod
,
q
.
q0
));
if
(
q
.
q0
!=
0.0
)
{
cvm
::
real
const
d_cos_theta_dq0
=
(
4.0
*
q
.
q0
/
(
cos_spin_2
*
cos_spin_2
))
*
(
1.0
-
(
iprod
*
iprod
)
/
(
q
.
q0
*
q
.
q0
)
/
(
1.0
+
(
iprod
*
iprod
)
/
(
q
.
q0
*
q
.
q0
)));
cvm
::
real
const
d_cos_theta_dqn
=
(
4.0
*
q
.
q0
/
(
cos_spin_2
*
cos_spin_2
)
*
(
iprod
/
q
.
q0
)
/
(
1.0
+
(
iprod
*
iprod
)
/
(
q
.
q0
*
q
.
q0
)));
return
cvm
::
quaternion
(
d_cos_theta_dq0
,
d_cos_theta_dqn
*
axis
.
x
,
d_cos_theta_dqn
*
axis
.
y
,
d_cos_theta_dqn
*
axis
.
z
);
}
else
{
cvm
::
real
const
d_cos_theta_dqn
=
(
4.0
/
(
cos_spin_2
*
cos_spin_2
*
iprod
));
return
cvm
::
quaternion
(
0.0
,
d_cos_theta_dqn
*
axis
.
x
,
d_cos_theta_dqn
*
axis
.
y
,
d_cos_theta_dqn
*
axis
.
z
);
}
}
/// \brief Whether to test for eigenvalue crossing
static
bool
monitor_crossings
;
/// \brief Threshold for the eigenvalue crossing test
static
cvm
::
real
crossing_threshold
;
protected:
/// \brief Previous value of the rotation (used to warn the user
/// when the structure changes too much, and there may be an
/// eigenvalue crossing)
cvm
::
quaternion
q_old
;
/// Build the overlap matrix S (used by calc_optimal_rotation())
void
build_matrix
(
std
::
vector
<
cvm
::
atom_pos
>
const
&
pos1
,
std
::
vector
<
cvm
::
atom_pos
>
const
&
pos2
,
cvm
::
matrix2d
<
cvm
::
real
>
&
S
);
/// Diagonalize the overlap matrix S (used by calc_optimal_rotation())
void
diagonalize_matrix
(
cvm
::
matrix2d
<
cvm
::
real
>
&
S
,
cvm
::
vector1d
<
cvm
::
real
>
&
S_eigval
,
cvm
::
matrix2d
<
cvm
::
real
>
&
S_eigvec
);
};
#endif
Event Timeline
Log In to Comment