Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    r673c7f ra67d19  
    2727  bool status = true;
    2828  const Vector *Center = DetermineCenterOfAll();
    29   double * const cell_size = World::getInstance().getDomain();
     29  double * const cell_size = World::get()->cell_size;
    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::getInstance().getDomain();
     50  double * const cell_size = World::get()->cell_size;
    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->at(i) = ptr->x[i];
    75       min->at(i) = ptr->x[i];
     74      max->x[i] = ptr->x.x[i];
     75      min->x[i] = ptr->x.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->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);
     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];
    8383      }
    8484    }
    … …  
    8989//    Log() << Verbose(0) << endl;
    9090    min->Scale(-1.);
    91     (*max) += (*min);
     91    max->AddVector(min);
    9292    Translate(min);
    9393    Center.Zero();
    … …  
    109109
    110110  if (ptr->next != end) {   //list not empty?
    111     while (ptr->next != end) {
     111    while (ptr->next != end) {  // continue with second if present
    112112      ptr = ptr->next;
    113113      Num++;
    114       Center += ptr->x;
     114      Center.AddVector(&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  */
    148 Vector * molecule::DetermineCenterOfGravity()
    149126{
    150127  atom *ptr = start->next;  // start at first in list
    … …  
    158135    while (ptr->next != end) {  // continue with second if present
    159136      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 */
     150Vector * 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;
    160162      Num += ptr->type->mass;
    161       tmp = ptr->type->mass * ptr->x;
    162       (*a) += tmp;
     163      tmp.CopyVector(&ptr->x);
     164      tmp.Scale(ptr->type->mass);  // scale by mass
     165      a->AddVector(&tmp);
    163166    }
    164167    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    … …  
    186189void molecule::CenterAtVector(Vector *newcenter)
    187190{
    188   Center = *newcenter;
     191  Center.CopyVector(newcenter);
    189192};
    190193
    … …  
    192195/** Scales all atoms by \a *factor.
    193196 * \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)
    199197 */
    200198void molecule::Scale(const double ** const factor)
    … …  
    205203    ptr = ptr->next;
    206204    for (int j=0;j<MDSteps;j++)
    207       ptr->Trajectory.R.at(j).ScaleAll(*factor);
    208     ptr->x.ScaleAll(*factor);
     205      ptr->Trajectory.R.at(j).Scale(factor);
     206    ptr->x.Scale(factor);
    209207  }
    210208};
    … …  
    220218    ptr = ptr->next;
    221219    for (int j=0;j<MDSteps;j++)
    222       ptr->Trajectory.R.at(j) += (*trans);
    223     ptr->x += (*trans);
     220      ptr->Trajectory.R.at(j).Translate(trans);
     221    ptr->x.Translate(trans);
    224222  }
    225223};
    … …  
    231229void molecule::TranslatePeriodically(const Vector *trans)
    232230{
    233   double * const cell_size = World::getInstance().getDomain();
     231  double * const cell_size = World::get()->cell_size;
    234232  double *M = ReturnFullMatrixforSymmetric(cell_size);
    235233  double *Minv = InverseMatrix(M);
    236234
    237235  // go through all atoms
    238   ActOnAllVectors( &Vector::SubtractVector, *trans);
     236  ActOnAllVectors( &Vector::SubtractVector, trans);
    239237  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    240238
    … …  
    249247void molecule::Mirror(const Vector *n)
    250248{
    251   ActOnAllVectors( &Vector::Mirror, *n );
     249  ActOnAllVectors( &Vector::Mirror, n );
    252250};
    253251
    … …  
    258256{
    259257  atom *Walker = start;
    260   double * const cell_size = World::getInstance().getDomain();
     258  double * const cell_size = World::get()->cell_size;
    261259  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    262260  double *inversematrix = InverseMatrix(cell_size);
    … …  
    273271      if (Walker->type->Z != 1) {
    274272#endif
    275         Testvector = Walker->x;
     273        Testvector.CopyVector(&Walker->x);
    276274        Testvector.MatrixMultiplication(inversematrix);
    277275        Translationvector.Zero();
    … …  
    279277         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    280278            for (int j=0;j<NDIM;j++) {
    281               tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
     279              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
    282280              if ((fabs(tmp)) > BondDistance) {
    283281                flag = false;
    284282                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    285283                if (tmp > 0)
    286                   Translationvector[j] -= 1.;
     284                  Translationvector.x[j] -= 1.;
    287285                else
    288                   Translationvector[j] += 1.;
     286                  Translationvector.x[j] += 1.;
    289287              }
    290288            }
    291289        }
    292         Testvector += Translationvector;
     290        Testvector.AddVector(&Translationvector);
    293291        Testvector.MatrixMultiplication(matrix);
    294         Center += Testvector;
    295         Log() << Verbose(1) << "vector is: " << Testvector << endl;
     292        Center.AddVector(&Testvector);
     293        DoLog(1) && (Log() << Verbose(1) << "vector is: ");
     294        Testvector.Output();
     295        DoLog(0) && (Log() << Verbose(0) << 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 = (*Runner)->GetOtherAtom(Walker)->x;
     300            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
    301301            Testvector.MatrixMultiplication(inversematrix);
    302             Testvector += Translationvector;
     302            Testvector.AddVector(&Translationvector);
    303303            Testvector.MatrixMultiplication(matrix);
    304             Center += Testvector;
    305             Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
     304            Center.AddVector(&Testvector);
     305            DoLog(1) && (Log() << Verbose(1) << "Hydrogen vector is: ");
     306            Testvector.Output();
     307            DoLog(0) && (Log() << Verbose(0) << endl);
    306308          }
    307309        }
    … …  
    336338  while (ptr->next != end) {
    337339    ptr = ptr->next;
    338     Vector x = ptr->x;
     340    Vector x;
     341    x.CopyVector(&ptr->x);
    339342    //x.SubtractVector(CenterOfGravity);
    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]);
     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]);
    349352  }
    350353  // print InertiaTensor for debugging
    … …  
    387390    while (ptr->next != end) {
    388391      ptr = ptr->next;
    389       Vector x = ptr->x;
     392      Vector x;
     393      x.CopyVector(&ptr->x);
    390394      //x.SubtractVector(CenterOfGravity);
    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]);
     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]);
    400404    }
    401405    // print InertiaTensor for debugging
    … …  
    424428  double alpha, tmp;
    425429  Vector z_axis;
    426   z_axis[0] = 0.;
    427   z_axis[1] = 0.;
    428   z_axis[2] = 1.;
     430  z_axis.x[0] = 0.;
     431  z_axis.x[1] = 0.;
     432  z_axis.x[2] = 1.;
    429433
    430434  // rotate on z-x plane
    431435  DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
    432   alpha = atan(-n->at(0)/n->at(2));
     436  alpha = atan(-n->x[0]/n->x[2]);
    433437  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    434438  while (ptr->next != end) {
    435439    ptr = ptr->next;
    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];
     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];
    439443    for (int j=0;j<MDSteps;j++) {
    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];
     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];
    443447    }
    444448  }
    445449  // rotate n vector
    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);
     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);
    450456
    451457  // rotate on z-y plane
    452458  ptr = start;
    453   alpha = atan(-n->at(1)/n->at(2));
     459  alpha = atan(-n->x[1]/n->x[2]);
    454460  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    455461  while (ptr->next != end) {
    456462    ptr = ptr->next;
    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];
     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];
    460466    for (int j=0;j<MDSteps;j++) {
    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];
     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];
    464470    }
    465471  }
    466472  // rotate n vector (for consistency check)
    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);
     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);
    473480  DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
    474481};
    … …  
    488495
    489496  // initialize vectors
    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);
     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);
    496503  // go through all atoms
    497504  while (ptr != par->mol->end) {
    498505    ptr = ptr->next;
    499506    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
    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
     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
    505514    }
    506515  }
Note: See TracChangeset for help on using the changeset viewer.