ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12811 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2014-12-02 02:09:53 +00:00
parent 9c0f849594
commit c4af1cdbbd
34 changed files with 5927 additions and 4786 deletions
Download Patch File Download Diff File Expand all files Collapse all files

301
lib/colvars/colvartypes.cpp
View File

@ -1,22 +1,47 @@
/// -*- c++ -*-
// -*- c++ -*-
#include <stdlib.h>
#include <string.h>
#include "colvarmodule.h"
#include "colvartypes.h"
#include "colvarparse.h"
std::string cvm::rvector::to_simple_string() const
{
std::ostringstream os;
os.setf(std::ios::scientific, std::ios::floatfield);
os.precision(cvm::cv_prec);
os << x << " " << y << " " << z;
return os.str();
}
int cvm::rvector::from_simple_string(std::string const &s)
{
std::stringstream stream(s);
if ( !(stream >> x) ||
!(stream >> y) ||
!(stream >> z) ) {
return COLVARS_ERROR;
}
return COLVARS_OK;
}
std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os.width(2);
os << "( ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << v.x << " , ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << v.y << " , ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << v.z << " )";
return os;
}
@ -31,29 +56,48 @@ std::istream & operator >> (std::istream &is, colvarmodule::rvector &v)
!(is >> v.y) || !(is >> sep) || !(sep == ',') ||
!(is >> v.z) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
return is;
}
std::string cvm::quaternion::to_simple_string() const
{
std::ostringstream os;
os.setf(std::ios::scientific, std::ios::floatfield);
os.precision(cvm::cv_prec);
os << q0 << " " << q1 << " " << q2 << " " << q3;
return os.str();
}
int cvm::quaternion::from_simple_string(std::string const &s)
{
std::stringstream stream(s);
if ( !(stream >> q0) ||
!(stream >> q1) ||
!(stream >> q2) ||
!(stream >> q3) ) {
return COLVARS_ERROR;
}
return COLVARS_OK;
}
std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os.width(2);
os << "( ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << q.q0 << " , ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << q.q1 << " , ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << q.q2 << " , ";
os.width (w); os.precision (p);
os.width(w); os.precision(p);
os << q.q3 << " )";
return os;
}
@ -63,10 +107,10 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
{
size_t const start_pos = is.tellg();
std::string euler ("");
std::string euler("");
if ( (is >> euler) && (colvarparse::to_lower_cppstr (euler) ==
std::string ("euler")) ) {
if ( (is >> euler) && (colvarparse::to_lower_cppstr(euler) ==
std::string("euler")) ) {
// parse the Euler angles
@ -77,18 +121,18 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
!(is >> theta) || !(is >> sep) || !(sep == ',') ||
!(is >> psi) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
q = colvarmodule::quaternion (phi, theta, psi);
q = colvarmodule::quaternion(phi, theta, psi);
} else {
// parse the quaternion components
is.seekg (start_pos, std::ios::beg);
is.seekg(start_pos, std::ios::beg);
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> q.q0) || !(is >> sep) || !(sep == ',') ||
@ -96,8 +140,8 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
!(is >> q.q2) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q3) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
}
@ -107,10 +151,10 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
cvm::quaternion
cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
cvm::quaternion::position_derivative_inner(cvm::rvector const &pos,
cvm::rvector const &vec) const
{
cvm::quaternion result (0.0, 0.0, 0.0, 0.0);
cvm::quaternion result(0.0, 0.0, 0.0, 0.0);
result.q0 = 2.0 * pos.x * q0 * vec.x
@ -174,22 +218,20 @@ cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
// Calculate the optimal rotation between two groups, and implement it
// as a quaternion. Uses the method documented 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 colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2,
matrix2d<cvm::real, 4, 4> &S)
void colvarmodule::rotation::build_matrix(std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2,
cvm::matrix2d<cvm::real> &S)
{
cvm::rmatrix C;
// build the correlation matrix
C.resize(3, 3);
C.reset();
for (size_t i = 0; i < pos1.size(); i++) {
size_t i;
for (i = 0; i < pos1.size(); i++) {
C.xx() += pos1[i].x * pos2[i].x;
C.xy() += pos1[i].x * pos2[i].y;
C.xz() += pos1[i].x * pos2[i].z;
@ -219,120 +261,101 @@ void colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos
S[3][2] = C.yz() + C.zy();
S[2][3] = S[3][2];
S[3][3] = - C.xx() - C.yy() + C.zz();
// if (cvm::debug()) {
// for (size_t i = 0; i < 4; i++) {
// std::string line ("");
// for (size_t j = 0; j < 4; j++) {
// line += std::string (" S["+cvm::to_str (i)+
// "]["+cvm::to_str (j)+"] ="+cvm::to_str (S[i][j]));
// }
// cvm::log (line+"\n");
// }
// }
}
void colvarmodule::rotation::diagonalize_matrix (matrix2d<cvm::real, 4, 4> &S,
cvm::real S_eigval[4],
matrix2d<cvm::real, 4, 4> &S_eigvec)
void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
cvm::vector1d<cvm::real> &S_eigval,
cvm::matrix2d<cvm::real> &S_eigvec)
{
S_eigval.resize(4);
S_eigval.reset();
S_eigvec.resize(4,4);
S_eigvec.reset();
// diagonalize
int jac_nrot = 0;
jacobi (S, S_eigval, S_eigvec, &jac_nrot);
eigsrt (S_eigval, S_eigvec);
jacobi(S.c_array(), S_eigval.c_array(), S_eigvec.c_array(), &jac_nrot);
eigsrt(S_eigval.c_array(), S_eigvec.c_array());
// jacobi saves eigenvectors by columns
transpose (S_eigvec);
transpose(S_eigvec.c_array());
// normalize eigenvectors
for (size_t ie = 0; ie < 4; ie++) {
cvm::real norm2 = 0.0;
size_t i;
for (i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2));
cvm::real const norm = std::sqrt (norm2);
for (i = 0; i < 4; i++) S_eigvec[ie][i] /= norm;
for (i = 0; i < 4; i++) {
norm2 += std::pow(S_eigvec[ie][i], int(2));
}
cvm::real const norm = std::sqrt(norm2);
for (i = 0; i < 4; i++) {
S_eigvec[ie][i] /= norm;
}
}
}
// Calculate the rotation, plus its derivatives
void colvarmodule::rotation::calc_optimal_rotation
(std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2)
void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2)
{
matrix2d<cvm::real, 4, 4> S;
matrix2d<cvm::real, 4, 4> S_backup;
cvm::real S_eigval[4];
matrix2d<cvm::real, 4, 4> S_eigvec;
S.resize(4,4);
S.reset();
// if (cvm::debug()) {
// cvm::atom_pos cog1 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos1.size(); i++) {
// cog1 += pos1[i];
// }
// cog1 /= cvm::real (pos1.size());
// cvm::atom_pos cog2 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos2.size(); i++) {
// cog2 += pos2[i];
// }
// cog2 /= cvm::real (pos1.size());
// cvm::log ("calc_optimal_rotation: centers of geometry are: "+
// cvm::to_str (cog1, cvm::cv_width, cvm::cv_prec)+
// " and "+cvm::to_str (cog2, cvm::cv_width, cvm::cv_prec)+".\n");
// }
build_matrix(pos1, pos2, S);
build_matrix (pos1, pos2, S);
S_backup.resize(4,4);
S_backup = S;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("S = "+cvm::to_str (cvm::to_str (S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
cvm::log("S = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
}
}
diagonalize_matrix (S, S_eigval, S_eigvec);
diagonalize_matrix(S, S_eigval, S_eigvec);
// eigenvalues and eigenvectors
cvm::real const &L0 = S_eigval[0];
cvm::real const &L1 = S_eigval[1];
cvm::real const &L2 = S_eigval[2];
cvm::real const &L3 = S_eigval[3];
cvm::real const *Q0 = S_eigvec[0];
cvm::real const *Q1 = S_eigvec[1];
cvm::real const *Q2 = S_eigvec[2];
cvm::real const *Q3 = S_eigvec[3];
cvm::real const L0 = S_eigval[0];
cvm::real const L1 = S_eigval[1];
cvm::real const L2 = S_eigval[2];
cvm::real const L3 = S_eigval[3];
cvm::quaternion const Q0(S_eigvec[0]);
cvm::quaternion const Q1(S_eigvec[1]);
cvm::quaternion const Q2(S_eigvec[2]);
cvm::quaternion const Q3(S_eigvec[3]);
lambda = L0;
q = cvm::quaternion (Q0);
q = Q0;
if (q_old.norm2() > 0.0) {
q.match (q_old);
if (q_old.inner (q) < (1.0 - crossing_threshold)) {
cvm::log ("Warning: one molecular orientation has changed by more than "+
cvm::to_str (crossing_threshold)+": discontinuous rotation ?\n");
q.match(q_old);
if (q_old.inner(q) < (1.0 - crossing_threshold)) {
cvm::log("Warning: one molecular orientation has changed by more than "+
cvm::to_str(crossing_threshold)+": discontinuous rotation ?\n");
}
}
q_old = q;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str (cvm::quaternion (Q0), cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q0)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L1 = "+cvm::to_str (L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str (cvm::quaternion (Q1), cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q1)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L2 = "+cvm::to_str (L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str (cvm::quaternion (Q2), cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q2)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L3 = "+cvm::to_str (L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str (cvm::quaternion (Q3), cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q3)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
@ -346,7 +369,7 @@ void colvarmodule::rotation::calc_optimal_rotation
cvm::real const &a2y = pos2[ia].y;
cvm::real const &a2z = pos2[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_1 = dS_1[ia];
cvm::matrix2d<cvm::rvector> &ds_1 = dS_1[ia];
// derivative of the S matrix
ds_1.reset();
@ -367,8 +390,8 @@ void colvarmodule::rotation::calc_optimal_rotation
ds_1[2][3] = ds_1[3][2];
ds_1[3][3].set(-a2x, -a2y, a2z);
cvm::rvector &dl0_1 = dL0_1[ia];
vector1d<cvm::rvector, 4> &dq0_1 = dQ0_1[ia];
cvm::rvector &dl0_1 = dL0_1[ia];
cvm::vector1d<cvm::rvector> &dq0_1 = dQ0_1[ia];
// matrix multiplications; derivatives of L_0 and Q_0 are
// calculated using Hellmann-Feynman theorem (i.e. exploiting the
@ -401,7 +424,7 @@ void colvarmodule::rotation::calc_optimal_rotation
cvm::real const &a1y = pos1[ia].y;
cvm::real const &a1z = pos1[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_2 = dS_2[ia];
cvm::matrix2d<cvm::rvector> &ds_2 = dS_2[ia];
ds_2.reset();
ds_2[0][0].set( a1x, a1y, a1z);
@ -421,8 +444,8 @@ void colvarmodule::rotation::calc_optimal_rotation
ds_2[2][3] = ds_2[3][2];
ds_2[3][3].set(-a1x, -a1y, a1z);
cvm::rvector &dl0_2 = dL0_2[ia];
vector1d<cvm::rvector, 4> &dq0_2 = dQ0_2[ia];
cvm::rvector &dl0_2 = dL0_2[ia];
cvm::vector1d<cvm::rvector> &dq0_2 = dQ0_2[ia];
dl0_2.reset();
for (size_t i = 0; i < 4; i++) {
@ -447,9 +470,9 @@ void colvarmodule::rotation::calc_optimal_rotation
if (b_debug_gradients) {
matrix2d<cvm::real, 4, 4> S_new;
cvm::real S_new_eigval[4];
matrix2d<cvm::real, 4, 4> S_new_eigvec;
cvm::matrix2d<cvm::real> S_new(4, 4);
cvm::vector1d<cvm::real> S_new_eigval(4);
cvm::matrix2d<cvm::real> S_new_eigvec(4, 4);
// make an infitesimal move along each cartesian coordinate of
// this atom, and solve again the eigenvector problem
@ -464,25 +487,25 @@ void colvarmodule::rotation::calc_optimal_rotation
}
}
// cvm::log ("S_new = "+cvm::to_str (cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
// cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix (S_new, S_new_eigval, S_new_eigvec);
diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::real const *Q0_new = S_new_eigvec[0];
cvm::quaternion const Q0_new(S_new_eigvec[0]);
cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
cvm::quaternion const q0 (Q0);
cvm::quaternion const DQ0 (dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::quaternion const q0(Q0);
cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::log ( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str (std::fabs (L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str ((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
}
@ -494,21 +517,24 @@ void colvarmodule::rotation::calc_optimal_rotation
// Numerical Recipes routine for diagonalization
#define ROTATE(a,i,j,k,l) g=a[i][j]; \
h=a[k][l]; \
a[i][j]=g-s*(h+g*tau); \
h=a[k][l]; \
a[i][j]=g-s*(h+g*tau); \
a[k][l]=h+s*(g-h*tau);
#define n 4
void jacobi(cvm::real **a, cvm::real d[], cvm::real **v, int *nrot)
void jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
{
int j,iq,ip,i;
cvm::real tresh,theta,tau,t,sm,s,h,g,c;
std::vector<cvm::real> b (n, 0.0);
std::vector<cvm::real> z (n, 0.0);
cvm::vector1d<cvm::real> b(n);
cvm::vector1d<cvm::real> z(n);
for (ip=0;ip<n;ip++) {
for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
for (iq=0;iq<n;iq++) {
v[ip][iq]=0.0;
}
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++) {
@ -575,10 +601,10 @@ void jacobi(cvm::real **a, cvm::real d[], cvm::real **v, int *nrot)
z[ip]=0.0;
}
}
cvm::error ("Too many iterations in routine jacobi.\n");
cvm::error("Too many iterations in routine jacobi.\n");
}
void eigsrt(cvm::real d[], cvm::real **v)
void eigsrt(cvm::real *d, cvm::real **v)
{
int k,j,i;
cvm::real p;
@ -602,8 +628,9 @@ void eigsrt(cvm::real d[], cvm::real **v)
void transpose(cvm::real **v)
{
cvm::real p;
for (int i=0;i<n;i++) {
for (int j=i+1;j<n;j++) {
int i,j;
for (i=0;i<n;i++) {
for (j=i+1;j<n;j++) {
p=v[i][j];
v[i][j]=v[j][i];
v[j][i]=p;