Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    r112b09 rccf826  
    55 *      Author: heber
    66 */
    7 
    8 #include "Helpers/MemDebug.hpp"
    97
    108#include "atom.hpp"
    … …  
    1917#include "World.hpp"
    2018#include "Plane.hpp"
    21 #include <boost/foreach.hpp>
    22 
    2319
    2420/************************************* Functions for class molecule *********************************/
    … …  
    3228  bool status = true;
    3329  const Vector *Center = DetermineCenterOfAll();
    34   const Vector *CenterBox = DetermineCenterOfBox();
    3530  double * const cell_size = World::getInstance().getDomain();
    3631  double *M = ReturnFullMatrixforSymmetric(cell_size);
    … …  
    3934  // go through all atoms
    4035  ActOnAllVectors( &Vector::SubtractVector, *Center);
    41   ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
    4236  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    4337
    44   delete[](M);
    45   delete[](Minv);
     38  Free(&M);
     39  Free(&Minv);
    4640  delete(Center);
    4741  return status;
    … …  
    6256  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    6357
    64   delete[](M);
    65   delete[](Minv);
     58  Free(&M);
     59  Free(&Minv);
    6660  return status;
    6761};
    … …  
    7670
    7771//  Log() << Verbose(3) << "Begin of CenterEdge." << endl;
    78   molecule::const_iterator iter = begin();  // start at first in list
    79   if (iter != end()) { //list not empty?
     72  atom *ptr = start->next;  // start at first in list
     73  if (ptr != end) {  //list not empty?
    8074    for (int i=NDIM;i--;) {
    81       max->at(i) = (*iter)->x[i];
    82       min->at(i) = (*iter)->x[i];
    83     }
    84     for (; iter != end(); ++iter) {// continue with second if present
    85       //(*iter)->Output(1,1,out);
     75      max->at(i) = ptr->x[i];
     76      min->at(i) = ptr->x[i];
     77    }
     78    while (ptr->next != end) {  // continue with second if present
     79      ptr = ptr->next;
     80      //ptr->Output(1,1,out);
    8681      for (int i=NDIM;i--;) {
    87         max->at(i) = (max->at(i) < (*iter)->x[i]) ? (*iter)->x[i] : max->at(i);
    88         min->at(i) = (min->at(i) > (*iter)->x[i]) ? (*iter)->x[i] : min->at(i);
     82        max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
     83        min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(i);
    8984      }
    9085    }
    … …  
    110105{
    111106  int Num = 0;
    112   molecule::const_iterator iter = begin();  // start at first in list
     107  atom *ptr = start;  // start at first in list
    113108
    114109  Center.Zero();
    115110
    116   if (iter != end()) {   //list not empty?
    117     for (; iter != end(); ++iter) {  // continue with second if present
     111  if (ptr->next != end) {   //list not empty?
     112    while (ptr->next != end) {  // continue with second if present
     113      ptr = ptr->next;
    118114      Num++;
    119       Center += (*iter)->x;
     115      Center += ptr->x;
    120116    }
    121117    Center.Scale(-1./Num); // divide through total number (and sign for direction)
    … …  
    130126Vector * molecule::DetermineCenterOfAll() const
    131127{
    132   molecule::const_iterator iter = begin();  // start at first in list
    133   Vector *a = new Vector();
    134   double Num = 0;
    135 
    136   a->Zero();
    137 
    138   if (iter != end()) {   //list not empty?
    139     for (; iter != end(); ++iter) {  // continue with second if present
    140       Num++;
    141       (*a) += (*iter)->x;
    142     }
    143     a->Scale(1./Num); // divide through total mass (and sign for direction)
    144   }
    145   return a;
    146 };
    147 
    148 /** Returns vector pointing to center of the domain.
    149  * \return pointer to center of the domain
    150  */
    151 Vector * molecule::DetermineCenterOfBox() const
    152 {
    153   Vector *a = new Vector(0.5,0.5,0.5);
    154 
    155   const double *cell_size = World::getInstance().getDomain();
    156   double *M = ReturnFullMatrixforSymmetric(cell_size);
    157   a->MatrixMultiplication(M);
    158   delete[](M);
    159 
    160   return a;
    161 };
    162 
    163 /** Returns vector pointing to center of gravity.
    164  * \param *out output stream for debugging
    165  * \return pointer to center of gravity vector
    166  */
    167 Vector * molecule::DetermineCenterOfGravity()
    168 {
    169   molecule::const_iterator iter = begin();  // start at first in list
     128  atom *ptr = start->next;  // start at first in list
    170129  Vector *a = new Vector();
    171130  Vector tmp;
    … …  
    174133  a->Zero();
    175134
    176   if (iter != end()) {   //list not empty?
    177     for (; iter != end(); ++iter) {  // continue with second if present
    178       Num += (*iter)->type->mass;
    179       tmp = (*iter)->type->mass * (*iter)->x;
     135  if (ptr != end) {   //list not empty?
     136    while (ptr->next != end) {  // continue with second if present
     137      ptr = ptr->next;
     138      Num += 1.;
     139      tmp = ptr->x;
    180140      (*a) += tmp;
    181141    }
    182142    a->Scale(1./Num); // divide through total mass (and sign for direction)
     143  }
     144  return a;
     145};
     146
     147/** Returns vector pointing to center of gravity.
     148 * \param *out output stream for debugging
     149 * \return pointer to center of gravity vector
     150 */
     151Vector * molecule::DetermineCenterOfGravity()
     152{
     153  atom *ptr = start->next;  // start at first in list
     154  Vector *a = new Vector();
     155  Vector tmp;
     156  double Num = 0;
     157
     158  a->Zero();
     159
     160  if (ptr != end) {   //list not empty?
     161    while (ptr->next != end) {  // continue with second if present
     162      ptr = ptr->next;
     163      Num += ptr->type->mass;
     164      tmp = ptr->type->mass * ptr->x;
     165      (*a) += tmp;
     166    }
     167    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    183168  }
    184169//  Log() << Verbose(1) << "Resulting center of gravity: ";
    … …  
    218203void molecule::Scale(const double ** const factor)
    219204{
    220   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     205  atom *ptr = start;
     206
     207  while (ptr->next != end) {
     208    ptr = ptr->next;
    221209    for (int j=0;j<MDSteps;j++)
    222       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
    223     (*iter)->x.ScaleAll(*factor);
     210      ptr->Trajectory.R.at(j).ScaleAll(*factor);
     211    ptr->x.ScaleAll(*factor);
    224212  }
    225213};
    … …  
    230218void molecule::Translate(const Vector *trans)
    231219{
    232   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     220  atom *ptr = start;
     221
     222  while (ptr->next != end) {
     223    ptr = ptr->next;
    233224    for (int j=0;j<MDSteps;j++)
    234       (*iter)->Trajectory.R.at(j) += (*trans);
    235     (*iter)->x += (*trans);
     225      ptr->Trajectory.R.at(j) += (*trans);
     226    ptr->x += (*trans);
    236227  }
    237228};
    … …  
    248239
    249240  // go through all atoms
    250   ActOnAllVectors( &Vector::AddVector, *trans);
     241  ActOnAllVectors( &Vector::SubtractVector, *trans);
    251242  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    252243
    253   delete[](M);
    254   delete[](Minv);
     244  Free(&M);
     245  Free(&Minv);
    255246};
    256247
    … …  
    261252void molecule::Mirror(const Vector *n)
    262253{
    263   OBSERVE;
    264254  Plane p(*n,0);
    265   BOOST_FOREACH( atom* iter, atoms ){
    266     (*iter->node) = p.mirrorVector(*iter->node);
     255  // TODO: replace with simpler construct (e.g. Boost::foreach)
     256  // once the structure of the atom list is fully reworked
     257  atom *Walker = start;
     258  while (Walker->next != end) {
     259    Walker = Walker->next;
     260    (*Walker->node) = p.mirrorVector(*Walker->node);
    267261  }
    268262};
    … …  
    273267void molecule::DeterminePeriodicCenter(Vector &center)
    274268{
     269  atom *Walker = start;
    275270  double * const cell_size = World::getInstance().getDomain();
    276271  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    277   double *inversematrix = InverseMatrix(matrix);
     272  double *inversematrix = InverseMatrix(cell_size);
    278273  double tmp;
    279274  bool flag;
    … …  
    283278    Center.Zero();
    284279    flag = true;
    285     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     280    while (Walker->next != end) {
     281      Walker = Walker->next;
    286282#ifdef ADDHYDROGEN
    287       if ((*iter)->type->Z != 1) {
     283      if (Walker->type->Z != 1) {
    288284#endif
    289         Testvector = (*iter)->x;
     285        Testvector = Walker->x;
    290286        Testvector.MatrixMultiplication(inversematrix);
    291287        Translationvector.Zero();
    292         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    293          if ((*iter)->nr < (*Runner)->GetOtherAtom((*iter))->nr) // otherwise we shift one to, the other fro and gain nothing
     288        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
     289         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    294290            for (int j=0;j<NDIM;j++) {
    295               tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
     291              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
    296292              if ((fabs(tmp)) > BondDistance) {
    297293                flag = false;
    298                 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << (*iter)->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
     294                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    299295                if (tmp > 0)
    300296                  Translationvector[j] -= 1.;
    … …  
    310306#ifdef ADDHYDROGEN
    311307        // now also change all hydrogens
    312         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    313           if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
    314             Testvector = (*Runner)->GetOtherAtom((*iter))->x;
     308        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
     309          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
     310            Testvector = (*Runner)->GetOtherAtom(Walker)->x;
    315311            Testvector.MatrixMultiplication(inversematrix);
    316312            Testvector += Translationvector;
    … …  
    324320    }
    325321  } while (!flag);
    326   delete[](matrix);
    327   delete[](inversematrix);
    328 
    329   Center.Scale(1./static_cast<double>(getAtomCount()));
     322  Free(&matrix);
     323  Free(&inversematrix);
     324
     325  Center.Scale(1./(double)AtomCount);
    330326};
    331327
    … …  
    337333void molecule::PrincipalAxisSystem(bool DoRotate)
    338334{
     335  atom *ptr = start;  // start at first in list
    339336  double InertiaTensor[NDIM*NDIM];
    340337  Vector *CenterOfGravity = DetermineCenterOfGravity();
    … …  
    347344
    348345  // sum up inertia tensor
    349   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    350     Vector x = (*iter)->x;
     346  while (ptr->next != end) {
     347    ptr = ptr->next;
     348    Vector x = ptr->x;
    351349    //x.SubtractVector(CenterOfGravity);
    352     InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    353     InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    354     InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    355     InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    356     InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    357     InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    358     InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    359     InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    360     InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     350    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     351    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     352    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     353    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     354    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     355    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     356    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     357    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     358    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    361359  }
    362360  // print InertiaTensor for debugging
    … …  
    396394
    397395    // sum up inertia tensor
    398     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    399       Vector x = (*iter)->x;
    400       InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    401       InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    402       InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    403       InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    404       InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    405       InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    406       InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    407       InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    408       InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     396    ptr = start;
     397    while (ptr->next != end) {
     398      ptr = ptr->next;
     399      Vector x = ptr->x;
     400      //x.SubtractVector(CenterOfGravity);
     401      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     402      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     403      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     404      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     405      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     406      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     407      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     408      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     409      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    409410    }
    410411    // print InertiaTensor for debugging
    … …  
    430431void molecule::Align(Vector *n)
    431432{
     433  atom *ptr = start;
    432434  double alpha, tmp;
    433435  Vector z_axis;
    … …  
    440442  alpha = atan(-n->at(0)/n->at(2));
    441443  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    442   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    443     tmp = (*iter)->x[0];
    444     (*iter)->x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
    445     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     444  while (ptr->next != end) {
     445    ptr = ptr->next;
     446    tmp = ptr->x[0];
     447    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     448    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    446449    for (int j=0;j<MDSteps;j++) {
    447       tmp = (*iter)->Trajectory.R.at(j)[0];
    448       (*iter)->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
    449       (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     450      tmp = ptr->Trajectory.R.at(j)[0];
     451      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     452      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    450453    }
    451454  }
    … …  
    457460
    458461  // rotate on z-y plane
     462  ptr = start;
    459463  alpha = atan(-n->at(1)/n->at(2));
    460464  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    461   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    462     tmp = (*iter)->x[1];
    463     (*iter)->x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
    464     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     465  while (ptr->next != end) {
     466    ptr = ptr->next;
     467    tmp = ptr->x[1];
     468    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     469    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    465470    for (int j=0;j<MDSteps;j++) {
    466       tmp = (*iter)->Trajectory.R.at(j)[1];
    467       (*iter)->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
    468       (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     471      tmp = ptr->Trajectory.R.at(j)[1];
     472      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     473      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    469474    }
    470475  }
    … …  
    490495  Vector a,b,c,d;
    491496  struct lsq_params *par = (struct lsq_params *)params;
     497  atom *ptr = par->mol->start;
    492498
    493499  // initialize vectors
    … …  
    499505  b[2] = gsl_vector_get(x,5);
    500506  // go through all atoms
    501   for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
    502     if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
    503       c = (*iter)->x - a;
     507  while (ptr != par->mol->end) {
     508    ptr = ptr->next;
     509    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
     510      c = ptr->x - a;
    504511      t = c.ScalarProduct(b);           // get direction parameter
    505512      d = t*b;       // and create vector
Note: See TracChangeset for help on using the changeset viewer.