Changeset e0c5b1 for molecuilder/src/vector.cpp

Timestamp:

Dec 3, 2008, 2:12:05 PM (17 years ago)

Author:

Christian Neuen <neuen@…>

Children:

3e20fe

Parents:

f5b58e

Message:

Multiple changes to boundary, currently not fully operational.
Molecules has a new routine to create adjacency lists, reading bonds from a dbond file instead of looking for the distances by itself.
Vector function Project onto plane has been updated.

File:

• molecuilder/src/vector.cpp (modified) (34 diffs)

molecuilder/src/vector.cpp

@@ -1 +1 @@
 /** \file vector.cpp
- * 
+ *
  * Function implementations for the class vector.
- * 
- */
- 
+ *
+ */
+
 #include "molecules.hpp"
- 
+
 
 /************************************ Functions for class vector ************************************/
@@ -31 +31 @@
   for (int i=NDIM;i--;)
     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (res);  
+  return (res);
 };
 
@@ -69 +69 @@
         if (tmp < res) res = tmp;
       }
-  return (res);  
+  return (res);
 };
 
@@ -112 +112 @@
   for (int i=NDIM;i--;)
     res += x[i]*y->x[i];
-  return (res);  
+  return (res);
 };
 
@@ -121 +121 @@
  *  \param *y array to vector with which to calculate crossproduct
  *  \return \f$ x \times y \f&
- */ 
+ */
 void Vector::VectorProduct(const Vector *y)
 {
   Vector tmp;
-  tmp[0] = x[1]*y->x[2] - x[2]*y->x[1];
-  tmp[1] = x[2]*y->x[0] - x[0]*y->x[2];
-  tmp[2] = x[0]*y->x[1] - x[1]*Y->x[0];
+  tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
+  tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
+  tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
   this->CopyVector(&tmp);
 
@@ -134 +134 @@
 
 /** projects this vector onto plane defined by \a *y.
- * \param *y array to normal vector of plane
+ * \param *y normal vector of plane
  * \return \f$\langle x, y \rangle\f$
  */
@@ -141 +141 @@
   Vector tmp;
   tmp.CopyVector(y);
-  tmp.Scale(Projection(y));
+  tmp.Normalize();
+  tmp.Scale(ScalarProduct(&tmp));
   this->SubtractVector(&tmp);
 };
@@ -162 +163 @@
   for (int i=NDIM;i--;)
     res += this->x[i]*this->x[i];
-  return (sqrt(res));  
+  return (sqrt(res));
 };
 
@@ -196 +197 @@
  */
 void Vector::Init(double x1, double x2, double x3)
-{ 
+{
   x[0] = x1;
   x[1] = x2;
@@ -202 +203 @@
 };
 
-/** Calculates the angle between this and another vector. 
+/** Calculates the angle between this and another vector.
  * \param *y array to second vector
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
@@ -208 +209 @@
 double Vector::Angle(Vector *y) const
 {
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
 };
 
@@ -256 +257 @@
 
 /** Sums two vectors \a  and \b component-wise.
- * \param a first vector 
+ * \param a first vector
  * \param b second vector
  * \return a + b
@@ -269 +270 @@
 
 /** Factors given vector \a a times \a m.
- * \param a vector 
+ * \param a vector
  * \param m factor
  * \return a + b
@@ -327 +328 @@
 };
 
-/** Translate atom by given vector. 
+/** Translate atom by given vector.
  * \param trans[] translation vector.
  */
@@ -334 +335 @@
   for (int i=NDIM;i--;)
     x[i] += trans->x[i];
-}; 
+};
 
 /** Do a matrix multiplication.
@@ -373 +374 @@
     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-  
+
     // do the matrix multiplication
     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
@@ -397 +398 @@
 {
   for(int i=NDIM;i--;)
-    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i]; 
-};
-
-/** Mirrors atom against a given plane. 
+    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
+};
+
+/** Mirrors atom against a given plane.
  * \param n[] normal vector of mirror plane.
  */
@@ -416 +417 @@
   Output((ofstream *)&cout);
   cout << endl;
-}; 
+};
 
 /** Calculates normal vector for three given vectors (being three points in space).
@@ -448 +449 @@
   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
   Normalize();
-  
+
   return true;
 };
@@ -506 +507 @@
 /** Creates this vector as one of the possible orthonormal ones to the given one.
  * Just scan how many components of given *vector are unequal to zero and
- * try to get the skp of both to be zero accordingly.  
+ * try to get the skp of both to be zero accordingly.
  * \param *vector given vector
  * \return true - success, false - failure (null vector given)
@@ -527 +528 @@
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
-        
+
   switch(Last) {
     case 3:  // threecomponent system
@@ -540 +541 @@
     case 1: // one component system
       // set sole non-zero component to 0, and one of the other zero component pendants to 1
-      x[(Components[0]+2)%NDIM] = 0.;   
-      x[(Components[0]+1)%NDIM] = 1.;   
-      x[Components[0]] = 0.;   
+      x[(Components[0]+2)%NDIM] = 0.;
+      x[(Components[0]+1)%NDIM] = 1.;
+      x[Components[0]] = 0.;
       return true;
       break;
@@ -559 +560 @@
 {
 //  cout << Verbose(3) << "For comparison: ";
-//  cout << "A " << A->Projection(this) << "\t"; 
-//  cout << "B " << B->Projection(this) << "\t"; 
-//  cout << "C " << C->Projection(this) << "\t"; 
+//  cout << "A " << A->Projection(this) << "\t";
+//  cout << "B " << B->Projection(this) << "\t";
+//  cout << "C " << C->Projection(this) << "\t";
 //  cout << endl;
   return A->Projection(this);
@@ -571 +572 @@
  * \return true if success, false if failed due to linear dependency
  */
-bool Vector::LSQdistance(Vector **vectors, int num) 
+bool Vector::LSQdistance(Vector **vectors, int num)
 {
         int j;
-                                
+
         for (j=0;j<num;j++) {
                 cout << Verbose(1) << j << "th atom's vector: ";
@@ -583 +584 @@
   int np = 3;
         struct LSQ_params par;
-    
+
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -589 +590 @@
    gsl_vector *ss, *y;
    gsl_multimin_function minex_func;
- 
+
    size_t iter = 0, i;
    int status;
    double size;
- 
+
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
    y = gsl_vector_alloc (np);
- 
+
    /* Set all step sizes to 1 */
    gsl_vector_set_all (ss, 1.0);
- 
+
    /* Starting point */
    par.vectors = vectors;
          par.num = num;
-        
+
          for (i=NDIM;i--;)
-                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
-         
+                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
+
    /* Initialize method and iterate */
    minex_func.f = &LSQ;
    minex_func.n = np;
    minex_func.params = (void *)&par;
- 
+
    s = gsl_multimin_fminimizer_alloc (T, np);
    gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
- 
+
    do
      {
        iter++;
        status = gsl_multimin_fminimizer_iterate(s);
- 
+
        if (status)
          break;
- 
+
        size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-2);
- 
+
        if (status == GSL_SUCCESS)
          {
            printf ("converged to minimum at\n");
          }
- 
+
        printf ("%5d ", (int)iter);
        for (i = 0; i < (size_t)np; i++)
@@ -640 +641 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
- 
+
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
@@ -706 +707 @@
  * \param alpha first angle
  * \param beta second angle
- * \param c norm of final vector 
+ * \param c norm of final vector
  * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
- * \bug this is not yet working properly 
+ * \bug this is not yet working properly
  */
 bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
@@ -724 +725 @@
   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     if (fabs(x1->x[1]) > MYEPSILON) {
-      flag = 1;    
+      flag = 1;
     } else if (fabs(x1->x[2]) > MYEPSILON) {
        flag = 2;
@@ -757 +758 @@
   // now comes the case system
   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
-  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 
+  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
   D3 = y->x[0]/x1->x[0]*A-B1;
   cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
   if (fabs(D1) < MYEPSILON) {
-    cout << Verbose(2) << "D1 == 0!\n"; 
+    cout << Verbose(2) << "D1 == 0!\n";
     if (fabs(D2) > MYEPSILON) {
-      cout << Verbose(3) << "D2 != 0!\n"; 
+      cout << Verbose(3) << "D2 != 0!\n";
       x[2] = -D3/D2;
       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
@@ -773 +774 @@
       cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
       if (fabs(F1) < MYEPSILON) {
-        cout << Verbose(4) << "F1 == 0!\n"; 
+        cout << Verbose(4) << "F1 == 0!\n";
         cout << Verbose(4) << "Gleichungssystem linear\n";
-        x[1] = F3/(2.*F2); 
+        x[1] = F3/(2.*F2);
       } else {
         p = F2/F1;
         q = p*p - F3/F1;
-        cout << Verbose(4) << "p " << p << "\tq " << q << endl;  
+        cout << Verbose(4) << "p " << p << "\tq " << q << endl;
         if (q < 0) {
           cout << Verbose(4) << "q < 0" << endl;
@@ -800 +801 @@
     cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     if (fabs(F1) < MYEPSILON) {
-      cout << Verbose(3) << "F1 == 0!\n"; 
+      cout << Verbose(3) << "F1 == 0!\n";
       cout << Verbose(3) << "Gleichungssystem linear\n";
-      x[2] = F3/(2.*F2);     
+      x[2] = F3/(2.*F2);
     } else {
       p = F2/F1;
       q = p*p - F3/F1;
-      cout << Verbose(3) << "p " << p << "\tq " << q << endl;  
+      cout << Verbose(3) << "p " << p << "\tq " << q << endl;
       if (q < 0) {
         cout << Verbose(3) << "q < 0" << endl;
@@ -850 +851 @@
     }
     cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
-    // apply sign matrix 
+    // apply sign matrix
     for (j=NDIM;j--;)
       x[j] *= sign[j];
@@ -856 +857 @@
     ang = x2->Angle (this);
     cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
-    if (fabs(ang - cos(beta)) < MYEPSILON) {  
+    if (fabs(ang - cos(beta)) < MYEPSILON) {
       break;
     }