Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

colvarmodule.C File Reference

#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias.h"
#include "colvarbias_meta.h"
#include "colvarbias_abf.h"

Go to the source code of this file.

Defines

#define ROTATE(a, i, j, k, l)

Functions

std::ostream & operator<< (std::ostream &os, colvarmodule::rvector const &v)
std::istream & operator>> (std::istream &is, colvarmodule::rvector &v)
std::ostream & operator<< (std::ostream &os, colvarmodule::quaternion const &q)
std::istream & operator>> (std::istream &is, colvarmodule::quaternion &q)
void jacobi (cvm::real **a, int n, cvm::real d[], cvm::real **v, int *nrot)
 Numerical recipes diagonalization.
void eigsrt (cvm::real d[], cvm::real **v, int n)
 Eigenvector sort.
void transpose (cvm::real **v, int n)
 Transpose the matrix.


Define Documentation

#define ROTATE a,
i,
j,
k,
 ) 
 

Value:

g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);    \
  a[k][l]=h+s*(g-h*tau);

Definition at line 1161 of file colvarmodule.C.

Referenced by jacobi().


Function Documentation

void eigsrt cvm::real  d[],
cvm::real **  v,
int  n
 

Eigenvector sort.

Definition at line 1244 of file colvarmodule.C.

References j.

01245 {
01246   int k,j,i;
01247   cvm::real p;
01248 
01249   for (i=0;i<n;i++) {
01250     p=d[k=i];
01251     for (j=i+1;j<n;j++)
01252       if (d[j] >= p) p=d[k=j];
01253     if (k != i) {
01254       d[k]=d[i];
01255       d[i]=p;
01256       for (j=0;j<n;j++) {
01257         p=v[j][i];
01258         v[j][i]=v[j][k];
01259         v[j][k]=p;
01260       }
01261     }
01262   }
01263 }

void jacobi cvm::real **  a,
int  n,
cvm::real  d[],
cvm::real **  v,
int *  nrot
 

Numerical recipes diagonalization.

Definition at line 1163 of file colvarmodule.C.

References colvarmodule::fatal_error(), j, and ROTATE.

01164 {
01165   int j,iq,ip,i;
01166   cvm::real tresh,theta,tau,t,sm,s,h,g,c;
01167 
01168   std::vector<cvm::real> b (n, 0.0);
01169   std::vector<cvm::real> z (n, 0.0);
01170 
01171   for (ip=0;ip<n;ip++) {
01172     for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
01173     v[ip][ip]=1.0;
01174   }
01175   for (ip=0;ip<n;ip++) {
01176     b[ip]=d[ip]=a[ip][ip];
01177     z[ip]=0.0;
01178   }
01179   *nrot=0;
01180   for (i=0;i<=50;i++) {
01181     sm=0.0;
01182     for (ip=0;ip<n-1;ip++) {
01183       for (iq=ip+1;iq<n;iq++)
01184         sm += fabs(a[ip][iq]);
01185     }
01186     if (sm == 0.0) {
01187       return;
01188     }
01189     if (i < 4)
01190       tresh=0.2*sm/(n*n);
01191     else
01192       tresh=0.0;
01193     for (ip=0;ip<n-1;ip++) {
01194       for (iq=ip+1;iq<n;iq++) {
01195         g=100.0*fabs(a[ip][iq]);
01196         if (i > 4 && (cvm::real)(fabs(d[ip])+g) == (cvm::real)fabs(d[ip])
01197             && (cvm::real)(fabs(d[iq])+g) == (cvm::real)fabs(d[iq]))
01198           a[ip][iq]=0.0;
01199         else if (fabs(a[ip][iq]) > tresh) {
01200           h=d[iq]-d[ip];
01201           if ((cvm::real)(fabs(h)+g) == (cvm::real)fabs(h))
01202             t=(a[ip][iq])/h;
01203           else {
01204             theta=0.5*h/(a[ip][iq]);
01205             t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
01206             if (theta < 0.0) t = -t;
01207           }
01208           c=1.0/sqrt(1+t*t);
01209           s=t*c;
01210           tau=s/(1.0+c);
01211           h=t*a[ip][iq];
01212           z[ip] -= h;
01213           z[iq] += h;
01214           d[ip] -= h;
01215           d[iq] += h;
01216           a[ip][iq]=0.0;
01217           for (j=0;j<=ip-1;j++) {
01218             ROTATE(a,j,ip,j,iq)
01219               }
01220           for (j=ip+1;j<=iq-1;j++) {
01221             ROTATE(a,ip,j,j,iq)
01222               }
01223           for (j=iq+1;j<n;j++) {
01224             ROTATE(a,ip,j,iq,j)
01225               }
01226           for (j=0;j<n;j++) {
01227             ROTATE(v,j,ip,j,iq)
01228               }
01229           ++(*nrot);
01230         }
01231       }
01232     }
01233     for (ip=0;ip<n;ip++) {
01234       b[ip] += z[ip];
01235       d[ip]=b[ip];
01236       z[ip]=0.0;
01237     }
01238   }
01239   cvm::fatal_error ("Too many iterations in routine jacobi.\n");
01240 }

std::ostream& operator<< std::ostream &  os,
colvarmodule::quaternion const &  q
 

Definition at line 773 of file colvarmodule.C.

References colvarmodule::quaternion::q0, colvarmodule::quaternion::q1, colvarmodule::quaternion::q2, and colvarmodule::quaternion::q3.

00774 {
00775   std::streamsize const w = os.width();
00776   std::streamsize const p = os.precision();
00777 
00778   os.width (2);
00779   os << "( ";
00780   os.width (w); os.precision (p);
00781   os << q.q0 << " , ";
00782   os.width (w); os.precision (p);
00783   os << q.q1 << " , ";
00784   os.width (w); os.precision (p);
00785   os << q.q2 << " , ";
00786   os.width (w); os.precision (p);
00787   os << q.q3 << " )";
00788   return os;
00789 }

std::ostream& operator<< std::ostream &  os,
colvarmodule::rvector const &  v
 

Definition at line 738 of file colvarmodule.C.

References colvarmodule::rvector::x, colvarmodule::rvector::y, and colvarmodule::rvector::z.

00739 {
00740   std::streamsize const w = os.width();
00741   std::streamsize const p = os.precision();
00742 
00743   os.width (2);
00744   os << "( ";
00745   os.width (w); os.precision (p);
00746   os << v.x << " , ";
00747   os.width (w); os.precision (p);
00748   os << v.y << " , ";
00749   os.width (w); os.precision (p);
00750   os << v.z << " )";
00751   return os;
00752 }

std::istream& operator>> std::istream &  is,
colvarmodule::quaternion q
 

Definition at line 792 of file colvarmodule.C.

References colvarmodule::quaternion::q0, colvarmodule::quaternion::q1, colvarmodule::quaternion::q2, colvarmodule::quaternion::q3, and colvarparse::to_lower_cppstr().

00793 {
00794   size_t const start_pos = is.tellg();
00795 
00796   std::string euler ("");
00797 
00798   if ( (is >> euler) && (colvarparse::to_lower_cppstr (euler) ==
00799                          std::string ("euler")) ) {
00800 
00801     // parse the Euler angles
00802     
00803     char sep;
00804     cvm::real phi, theta, psi;
00805     if ( !(is >> sep)   || !(sep == '(') ||
00806          !(is >> phi)   || !(is >> sep)  || !(sep == ',') ||
00807          !(is >> theta) || !(is >> sep)  || !(sep == ',') ||
00808          !(is >> psi)   || !(is >> sep)  || !(sep == ')') ) {
00809       is.clear();
00810       is.seekg (start_pos, std::ios::beg);
00811       is.setstate (std::ios::failbit);
00812       return is;
00813     }
00814 
00815     q = colvarmodule::quaternion (phi, theta, psi);
00816 
00817   } else {
00818 
00819     // parse the quaternion components
00820 
00821     is.seekg (start_pos, std::ios::beg);      
00822     char sep;
00823     if ( !(is >> sep)  || !(sep == '(') ||
00824          !(is >> q.q0) || !(is >> sep)  || !(sep == ',') ||
00825          !(is >> q.q1) || !(is >> sep)  || !(sep == ',') ||
00826          !(is >> q.q2) || !(is >> sep)  || !(sep == ',') ||
00827          !(is >> q.q3) || !(is >> sep)  || !(sep == ')') ) {
00828       is.clear();
00829       is.seekg (start_pos, std::ios::beg);
00830       is.setstate (std::ios::failbit);
00831       return is;
00832     }
00833   }
00834 
00835   return is;
00836 }

std::istream& operator>> std::istream &  is,
colvarmodule::rvector v
 

Definition at line 755 of file colvarmodule.C.

References colvarmodule::rvector::x, colvarmodule::rvector::y, and colvarmodule::rvector::z.

00756 {
00757   size_t const start_pos = is.tellg();
00758   char sep;
00759   if ( !(is >> sep) || !(sep == '(') ||
00760        !(is >> v.x) || !(is >> sep)  || !(sep == ',') ||
00761        !(is >> v.y) || !(is >> sep)  || !(sep == ',') ||
00762        !(is >> v.z) || !(is >> sep)  || !(sep == ')') ) {
00763     is.clear();
00764     is.seekg (start_pos, std::ios::beg);
00765     is.setstate (std::ios::failbit);
00766     return is;
00767   }
00768   return is;
00769 }

void transpose cvm::real **  v,
int  n
 

Transpose the matrix.

Definition at line 1266 of file colvarmodule.C.

References j.

01267 {
01268   cvm::real p;
01269   for (int i=0;i<n;i++) {
01270     for (int j=i+1;j<n;j++) {
01271       p=v[i][j];
01272       v[i][j]=v[j][i];
01273       v[j][i]=p;
01274     }
01275   }
01276 }


Generated on Mon Nov 23 04:59:26 2009 for NAMD by  doxygen 1.3.9.1