Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    ra67d19 r673c7f  
    2727  bool status = true;
    2828  const Vector *Center = DetermineCenterOfAll();
    29   double * const cell_size = World::get()->cell_size;
     29  double * const cell_size = World::getInstance().getDomain();
    3030  double *M = ReturnFullMatrixforSymmetric(cell_size);
    3131  double *Minv = InverseMatrix(M);
    3232
    3333  // go through all atoms
    34   ActOnAllVectors( &Vector::SubtractVector, Center);
     34  ActOnAllVectors( &Vector::SubtractVector, *Center);
    3535  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    3636
    … …  
    4848{
    4949  bool status = true;
    50   double * const cell_size = World::get()->cell_size;
     50  double * const cell_size = World::getInstance().getDomain();
    5151  double *M = ReturnFullMatrixforSymmetric(cell_size);
    5252  double *Minv = InverseMatrix(M);
    … …  
    7272  if (ptr != end) {   //list not empty?
    7373    for (int i=NDIM;i--;) {
    74       max->x[i] = ptr->x.x[i];
    75       min->x[i] = ptr->x.x[i];
     74      max->at(i) = ptr->x[i];
     75      min->at(i) = ptr->x[i];
    7676    }
    7777    while (ptr->next != end) {  // continue with second if present
    … …  
    7979      //ptr->Output(1,1,out);
    8080      for (int i=NDIM;i--;) {
    81         max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
    82         min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
     81        max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
     82        min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(i);
    8383      }
    8484    }
    … …  
    8989//    Log() << Verbose(0) << endl;
    9090    min->Scale(-1.);
    91     max->AddVector(min);
     91    (*max) += (*min);
    9292    Translate(min);
    9393    Center.Zero();
    … …  
    109109
    110110  if (ptr->next != end) {   //list not empty?
    111     while (ptr->next != end) {  // continue with second if present
     111    while (ptr->next != end) {
    112112      ptr = ptr->next;
    113113      Num++;
    114       Center.AddVector(&ptr->x);
     114      Center += ptr->x;
    115115    }
    116116    Center.Scale(-1./Num); // divide through total number (and sign for direction)
    … …  
    124124 */
    125125Vector * molecule::DetermineCenterOfAll() const
     126{
     127  atom *ptr = start;  // start at first in list
     128  Vector *a = new Vector();
     129  double Num = 0;
     130
     131  a->Zero();
     132
     133  if (ptr->next != end) {   //list not empty?
     134    while (ptr->next != end) {
     135      ptr = ptr->next;
     136      Num += 1.;
     137      (*a) += ptr->x;
     138    }
     139    a->Scale(1./Num); // divide through total mass (and sign for direction)
     140  }
     141  return a;
     142};
     143
     144/** Returns vector pointing to center of gravity.
     145 * \param *out output stream for debugging
     146 * \return pointer to center of gravity vector
     147 */
     148Vector * molecule::DetermineCenterOfGravity()
    126149{
    127150  atom *ptr = start->next;  // start at first in list
    … …  
    135158    while (ptr->next != end) {  // continue with second if present
    136159      ptr = ptr->next;
    137       Num += 1.;
    138       tmp.CopyVector(&ptr->x);
    139       a->AddVector(&tmp);
    140     }
    141     a->Scale(1./Num); // divide through total mass (and sign for direction)
    142   }
    143   return a;
    144 };
    145 
    146 /** Returns vector pointing to center of gravity.
    147  * \param *out output stream for debugging
    148  * \return pointer to center of gravity vector
    149  */
    150 Vector * molecule::DetermineCenterOfGravity()
    151 {
    152   atom *ptr = start->next;  // start at first in list
    153   Vector *a = new Vector();
    154   Vector tmp;
    155   double Num = 0;
    156 
    157   a->Zero();
    158 
    159   if (ptr != end) {   //list not empty?
    160     while (ptr->next != end) {  // continue with second if present
    161       ptr = ptr->next;
    162160      Num += ptr->type->mass;
    163       tmp.CopyVector(&ptr->x);
    164       tmp.Scale(ptr->type->mass);  // scale by mass
    165       a->AddVector(&tmp);
     161      tmp = ptr->type->mass * ptr->x;
     162      (*a) += tmp;
    166163    }
    167164    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    … …  
    189186void molecule::CenterAtVector(Vector *newcenter)
    190187{
    191   Center.CopyVector(newcenter);
     188  Center = *newcenter;
    192189};
    193190
    … …  
    195192/** Scales all atoms by \a *factor.
    196193 * \param *factor pointer to scaling factor
     194 *
     195 * TODO: Is this realy what is meant, i.e.
     196 * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
     197 * or rather
     198 * x=(**factor) * x (as suggested by comment)
    197199 */
    198200void molecule::Scale(const double ** const factor)
    … …  
    203205    ptr = ptr->next;
    204206    for (int j=0;j<MDSteps;j++)
    205       ptr->Trajectory.R.at(j).Scale(factor);
    206     ptr->x.Scale(factor);
     207      ptr->Trajectory.R.at(j).ScaleAll(*factor);
     208    ptr->x.ScaleAll(*factor);
    207209  }
    208210};
    … …  
    218220    ptr = ptr->next;
    219221    for (int j=0;j<MDSteps;j++)
    220       ptr->Trajectory.R.at(j).Translate(trans);
    221     ptr->x.Translate(trans);
     222      ptr->Trajectory.R.at(j) += (*trans);
     223    ptr->x += (*trans);
    222224  }
    223225};
    … …  
    229231void molecule::TranslatePeriodically(const Vector *trans)
    230232{
    231   double * const cell_size = World::get()->cell_size;
     233  double * const cell_size = World::getInstance().getDomain();
    232234  double *M = ReturnFullMatrixforSymmetric(cell_size);
    233235  double *Minv = InverseMatrix(M);
    234236
    235237  // go through all atoms
    236   ActOnAllVectors( &Vector::SubtractVector, trans);
     238  ActOnAllVectors( &Vector::SubtractVector, *trans);
    237239  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    238240
    … …  
    247249void molecule::Mirror(const Vector *n)
    248250{
    249   ActOnAllVectors( &Vector::Mirror, n );
     251  ActOnAllVectors( &Vector::Mirror, *n );
    250252};
    251253
    … …  
    256258{
    257259  atom *Walker = start;
    258   double * const cell_size = World::get()->cell_size;
     260  double * const cell_size = World::getInstance().getDomain();
    259261  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    260262  double *inversematrix = InverseMatrix(cell_size);
    … …  
    271273      if (Walker->type->Z != 1) {
    272274#endif
    273         Testvector.CopyVector(&Walker->x);
     275        Testvector = Walker->x;
    274276        Testvector.MatrixMultiplication(inversematrix);
    275277        Translationvector.Zero();
    … …  
    277279         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    278280            for (int j=0;j<NDIM;j++) {
    279               tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
     281              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
    280282              if ((fabs(tmp)) > BondDistance) {
    281283                flag = false;
    282284                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    283285                if (tmp > 0)
    284                   Translationvector.x[j] -= 1.;
     286                  Translationvector[j] -= 1.;
    285287                else
    286                   Translationvector.x[j] += 1.;
     288                  Translationvector[j] += 1.;
    287289              }
    288290            }
    289291        }
    290         Testvector.AddVector(&Translationvector);
     292        Testvector += Translationvector;
    291293        Testvector.MatrixMultiplication(matrix);
    292         Center.AddVector(&Testvector);
    293         DoLog(1) && (Log() << Verbose(1) << "vector is: ");
    294         Testvector.Output();
    295         DoLog(0) && (Log() << Verbose(0) << endl);
     294        Center += Testvector;
     295        Log() << Verbose(1) << "vector is: " << Testvector << endl;
    296296#ifdef ADDHYDROGEN
    297297        // now also change all hydrogens
    298298        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
    299299          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
    300             Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
     300            Testvector = (*Runner)->GetOtherAtom(Walker)->x;
    301301            Testvector.MatrixMultiplication(inversematrix);
    302             Testvector.AddVector(&Translationvector);
     302            Testvector += Translationvector;
    303303            Testvector.MatrixMultiplication(matrix);
    304             Center.AddVector(&Testvector);
    305             DoLog(1) && (Log() << Verbose(1) << "Hydrogen vector is: ");
    306             Testvector.Output();
    307             DoLog(0) && (Log() << Verbose(0) << endl);
     304            Center += Testvector;
     305            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
    308306          }
    309307        }
    … …  
    338336  while (ptr->next != end) {
    339337    ptr = ptr->next;
    340     Vector x;
    341     x.CopyVector(&ptr->x);
     338    Vector x = ptr->x;
    342339    //x.SubtractVector(CenterOfGravity);
    343     InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
    344     InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
    345     InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
    346     InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
    347     InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
    348     InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
    349     InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
    350     InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
    351     InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
     340    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     341    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     342    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     343    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     344    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     345    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     346    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     347    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     348    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    352349  }
    353350  // print InertiaTensor for debugging
    … …  
    390387    while (ptr->next != end) {
    391388      ptr = ptr->next;
    392       Vector x;
    393       x.CopyVector(&ptr->x);
     389      Vector x = ptr->x;
    394390      //x.SubtractVector(CenterOfGravity);
    395       InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
    396       InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
    397       InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
    398       InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
    399       InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
    400       InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
    401       InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
    402       InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
    403       InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
     391      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     392      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     393      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     394      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     395      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     396      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     397      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     398      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     399      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    404400    }
    405401    // print InertiaTensor for debugging
    … …  
    428424  double alpha, tmp;
    429425  Vector z_axis;
    430   z_axis.x[0] = 0.;
    431   z_axis.x[1] = 0.;
    432   z_axis.x[2] = 1.;
     426  z_axis[0] = 0.;
     427  z_axis[1] = 0.;
     428  z_axis[2] = 1.;
    433429
    434430  // rotate on z-x plane
    435431  DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
    436   alpha = atan(-n->x[0]/n->x[2]);
     432  alpha = atan(-n->at(0)/n->at(2));
    437433  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    438434  while (ptr->next != end) {
    439435    ptr = ptr->next;
    440     tmp = ptr->x.x[0];
    441     ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    442     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     436    tmp = ptr->x[0];
     437    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     438    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    443439    for (int j=0;j<MDSteps;j++) {
    444       tmp = ptr->Trajectory.R.at(j).x[0];
    445       ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
    446       ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
     440      tmp = ptr->Trajectory.R.at(j)[0];
     441      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     442      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    447443    }
    448444  }
    449445  // rotate n vector
    450   tmp = n->x[0];
    451   n->x[0] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    452   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    453   DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: ");
    454   n->Output();
    455   DoLog(0) && (Log() << Verbose(0) << endl);
     446  tmp = n->at(0);
     447  n->at(0) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
     448  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
     449  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
    456450
    457451  // rotate on z-y plane
    458452  ptr = start;
    459   alpha = atan(-n->x[1]/n->x[2]);
     453  alpha = atan(-n->at(1)/n->at(2));
    460454  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    461455  while (ptr->next != end) {
    462456    ptr = ptr->next;
    463     tmp = ptr->x.x[1];
    464     ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    465     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     457    tmp = ptr->x[1];
     458    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     459    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    466460    for (int j=0;j<MDSteps;j++) {
    467       tmp = ptr->Trajectory.R.at(j).x[1];
    468       ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
    469       ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
     461      tmp = ptr->Trajectory.R.at(j)[1];
     462      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     463      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    470464    }
    471465  }
    472466  // rotate n vector (for consistency check)
    473   tmp = n->x[1];
    474   n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    475   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    476 
    477   DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: ");
    478   n->Output();
    479   DoLog(1) && (Log() << Verbose(1) << endl);
     467  tmp = n->at(1);
     468  n->at(1) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
     469  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
     470
     471
     472  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
    480473  DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
    481474};
    … …  
    495488
    496489  // initialize vectors
    497   a.x[0] = gsl_vector_get(x,0);
    498   a.x[1] = gsl_vector_get(x,1);
    499   a.x[2] = gsl_vector_get(x,2);
    500   b.x[0] = gsl_vector_get(x,3);
    501   b.x[1] = gsl_vector_get(x,4);
    502   b.x[2] = gsl_vector_get(x,5);
     490  a[0] = gsl_vector_get(x,0);
     491  a[1] = gsl_vector_get(x,1);
     492  a[2] = gsl_vector_get(x,2);
     493  b[0] = gsl_vector_get(x,3);
     494  b[1] = gsl_vector_get(x,4);
     495  b[2] = gsl_vector_get(x,5);
    503496  // go through all atoms
    504497  while (ptr != par->mol->end) {
    505498    ptr = ptr->next;
    506499    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
    507       c.CopyVector(&ptr->x);  // copy vector to temporary one
    508       c.SubtractVector(&a);   // subtract offset vector
    509       t = c.ScalarProduct(&b);           // get direction parameter
    510       d.CopyVector(&b);       // and create vector
    511       d.Scale(&t);
    512       c.SubtractVector(&d);   // ... yielding distance vector
    513       res += d.ScalarProduct((const Vector *)&d);        // add squared distance
     500      c = ptr->x - a;
     501      t = c.ScalarProduct(b);           // get direction parameter
     502      d = t*b;       // and create vector
     503      c -= d;   // ... yielding distance vector
     504      res += d.ScalarProduct(d);        // add squared distance
    514505    }
    515506  }
Note: See TracChangeset for help on using the changeset viewer.