Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    ra67d19 r68f03d  
    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();
    … …  
    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)
    … …  
    136136      ptr = ptr->next;
    137137      Num += 1.;
    138       tmp.CopyVector(&ptr->x);
    139       a->AddVector(&tmp);
     138      tmp = ptr->x;
     139      (*a) += tmp;
    140140    }
    141141    a->Scale(1./Num); // divide through total mass (and sign for direction)
    … …  
    161161      ptr = ptr->next;
    162162      Num += ptr->type->mass;
    163       tmp.CopyVector(&ptr->x);
    164       tmp.Scale(ptr->type->mass);  // scale by mass
    165       a->AddVector(&tmp);
     163      tmp = ptr->type->mass * ptr->x;
     164      (*a) += tmp;
    166165    }
    167166    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    … …  
    189188void molecule::CenterAtVector(Vector *newcenter)
    190189{
    191   Center.CopyVector(newcenter);
     190  Center = *newcenter;
    192191};
    193192
    … …  
    195194/** Scales all atoms by \a *factor.
    196195 * \param *factor pointer to scaling factor
     196 *
     197 * TODO: Is this realy what is meant, i.e.
     198 * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
     199 * or rather
     200 * x=(**factor) * x (as suggested by comment)
    197201 */
    198202void molecule::Scale(const double ** const factor)
    … …  
    203207    ptr = ptr->next;
    204208    for (int j=0;j<MDSteps;j++)
    205       ptr->Trajectory.R.at(j).Scale(factor);
    206     ptr->x.Scale(factor);
     209      ptr->Trajectory.R.at(j).ScaleAll(*factor);
     210    ptr->x.ScaleAll(*factor);
    207211  }
    208212};
    … …  
    218222    ptr = ptr->next;
    219223    for (int j=0;j<MDSteps;j++)
    220       ptr->Trajectory.R.at(j).Translate(trans);
    221     ptr->x.Translate(trans);
     224      ptr->Trajectory.R.at(j) += (*trans);
     225    ptr->x += (*trans);
    222226  }
    223227};
    … …  
    229233void molecule::TranslatePeriodically(const Vector *trans)
    230234{
    231   double * const cell_size = World::get()->cell_size;
     235  double * const cell_size = World::getInstance().getDomain();
    232236  double *M = ReturnFullMatrixforSymmetric(cell_size);
    233237  double *Minv = InverseMatrix(M);
    234238
    235239  // go through all atoms
    236   ActOnAllVectors( &Vector::SubtractVector, trans);
     240  ActOnAllVectors( &Vector::SubtractVector, *trans);
    237241  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    238242
    … …  
    247251void molecule::Mirror(const Vector *n)
    248252{
    249   ActOnAllVectors( &Vector::Mirror, n );
     253  ActOnAllVectors( &Vector::Mirror, *n );
    250254};
    251255
    … …  
    256260{
    257261  atom *Walker = start;
    258   double * const cell_size = World::get()->cell_size;
     262  double * const cell_size = World::getInstance().getDomain();
    259263  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    260264  double *inversematrix = InverseMatrix(cell_size);
    … …  
    271275      if (Walker->type->Z != 1) {
    272276#endif
    273         Testvector.CopyVector(&Walker->x);
     277        Testvector = Walker->x;
    274278        Testvector.MatrixMultiplication(inversematrix);
    275279        Translationvector.Zero();
    … …  
    277281         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    278282            for (int j=0;j<NDIM;j++) {
    279               tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
     283              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
    280284              if ((fabs(tmp)) > BondDistance) {
    281285                flag = false;
    282                 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
     286                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    283287                if (tmp > 0)
    284                   Translationvector.x[j] -= 1.;
     288                  Translationvector[j] -= 1.;
    285289                else
    286                   Translationvector.x[j] += 1.;
     290                  Translationvector[j] += 1.;
    287291              }
    288292            }
    289293        }
    290         Testvector.AddVector(&Translationvector);
     294        Testvector += Translationvector;
    291295        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);
     296        Center += Testvector;
     297        Log() << Verbose(1) << "vector is: " << Testvector << endl;
    296298#ifdef ADDHYDROGEN
    297299        // now also change all hydrogens
    298300        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
    299301          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
    300             Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
     302            Testvector = (*Runner)->GetOtherAtom(Walker)->x;
    301303            Testvector.MatrixMultiplication(inversematrix);
    302             Testvector.AddVector(&Translationvector);
     304            Testvector += Translationvector;
    303305            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);
     306            Center += Testvector;
     307            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
    308308          }
    309309        }
    … …  
    338338  while (ptr->next != end) {
    339339    ptr = ptr->next;
    340     Vector x;
    341     x.CopyVector(&ptr->x);
     340    Vector x = ptr->x;
    342341    //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]);
     342    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     343    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     344    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     345    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     346    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     347    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     348    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     349    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     350    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    352351  }
    353352  // print InertiaTensor for debugging
    … …  
    390389    while (ptr->next != end) {
    391390      ptr = ptr->next;
    392       Vector x;
    393       x.CopyVector(&ptr->x);
     391      Vector x = ptr->x;
    394392      //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]);
     393      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     394      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     395      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     396      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     397      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     398      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     399      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     400      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     401      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    404402    }
    405403    // print InertiaTensor for debugging
    … …  
    428426  double alpha, tmp;
    429427  Vector z_axis;
    430   z_axis.x[0] = 0.;
    431   z_axis.x[1] = 0.;
    432   z_axis.x[2] = 1.;
     428  z_axis[0] = 0.;
     429  z_axis[1] = 0.;
     430  z_axis[2] = 1.;
    433431
    434432  // rotate on z-x plane
    435433  DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
    436   alpha = atan(-n->x[0]/n->x[2]);
     434  alpha = atan(-n->at(0)/n->at(2));
    437435  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    438436  while (ptr->next != end) {
    439437    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];
     438    tmp = ptr->x[0];
     439    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     440    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    443441    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];
     442      tmp = ptr->Trajectory.R.at(j)[0];
     443      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     444      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    447445    }
    448446  }
    449447  // 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);
     448  tmp = n->at(0);
     449  n->at(0) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
     450  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
     451  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
    456452
    457453  // rotate on z-y plane
    458454  ptr = start;
    459   alpha = atan(-n->x[1]/n->x[2]);
     455  alpha = atan(-n->at(1)/n->at(2));
    460456  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    461457  while (ptr->next != end) {
    462458    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];
     459    tmp = ptr->x[1];
     460    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     461    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    466462    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];
     463      tmp = ptr->Trajectory.R.at(j)[1];
     464      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     465      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    470466    }
    471467  }
    472468  // 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);
     469  tmp = n->at(1);
     470  n->at(1) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
     471  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
     472
     473
     474  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
    480475  DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
    481476};
    … …  
    495490
    496491  // 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);
     492  a[0] = gsl_vector_get(x,0);
     493  a[1] = gsl_vector_get(x,1);
     494  a[2] = gsl_vector_get(x,2);
     495  b[0] = gsl_vector_get(x,3);
     496  b[1] = gsl_vector_get(x,4);
     497  b[2] = gsl_vector_get(x,5);
    503498  // go through all atoms
    504499  while (ptr != par->mol->end) {
    505500    ptr = ptr->next;
    506501    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
     502      c = ptr->x - a;
     503      t = c.ScalarProduct(b);           // get direction parameter
     504      d = t*b;       // and create vector
     505      c -= d;   // ... yielding distance vector
     506      res += d.ScalarProduct(d);        // add squared distance
    514507    }
    515508  }
Note: See TracChangeset for help on using the changeset viewer.