Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    rccf826 r112b09  
    55 *      Author: heber
    66 */
     7
     8#include "Helpers/MemDebug.hpp"
    79
    810#include "atom.hpp"
    … …  
    1719#include "World.hpp"
    1820#include "Plane.hpp"
     21#include <boost/foreach.hpp>
     22
    1923
    2024/************************************* Functions for class molecule *********************************/
    … …  
    2832  bool status = true;
    2933  const Vector *Center = DetermineCenterOfAll();
     34  const Vector *CenterBox = DetermineCenterOfBox();
    3035  double * const cell_size = World::getInstance().getDomain();
    3136  double *M = ReturnFullMatrixforSymmetric(cell_size);
    … …  
    3439  // go through all atoms
    3540  ActOnAllVectors( &Vector::SubtractVector, *Center);
     41  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
    3642  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    3743
    38   Free(&M);
    39   Free(&Minv);
     44  delete[](M);
     45  delete[](Minv);
    4046  delete(Center);
    4147  return status;
    … …  
    5662  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    5763
    58   Free(&M);
    59   Free(&Minv);
     64  delete[](M);
     65  delete[](Minv);
    6066  return status;
    6167};
    … …  
    7076
    7177//  Log() << Verbose(3) << "Begin of CenterEdge." << endl;
    72   atom *ptr = start->next;  // start at first in list
    73   if (ptr != end) {  //list not empty?
     78  molecule::const_iterator iter = begin();  // start at first in list
     79  if (iter != end()) { //list not empty?
    7480    for (int i=NDIM;i--;) {
    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);
     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);
    8186      for (int i=NDIM;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);
     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);
    8489      }
    8590    }
    … …  
    105110{
    106111  int Num = 0;
    107   atom *ptr = start;  // start at first in list
     112  molecule::const_iterator iter = begin();  // start at first in list
    108113
    109114  Center.Zero();
    110115
    111   if (ptr->next != end) {   //list not empty?
    112     while (ptr->next != end) {  // continue with second if present
    113       ptr = ptr->next;
     116  if (iter != end()) {   //list not empty?
     117    for (; iter != end(); ++iter) {  // continue with second if present
    114118      Num++;
    115       Center += ptr->x;
     119      Center += (*iter)->x;
    116120    }
    117121    Center.Scale(-1./Num); // divide through total number (and sign for direction)
    … …  
    126130Vector * molecule::DetermineCenterOfAll() const
    127131{
    128   atom *ptr = start->next;  // start at first in list
     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 */
     151Vector * 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 */
     167Vector * molecule::DetermineCenterOfGravity()
     168{
     169  molecule::const_iterator iter = begin();  // start at first in list
    129170  Vector *a = new Vector();
    130171  Vector tmp;
    … …  
    133174  a->Zero();
    134175
    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;
     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;
    140180      (*a) += tmp;
    141181    }
    142182    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  */
    151 Vector * 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)
    168183  }
    169184//  Log() << Verbose(1) << "Resulting center of gravity: ";
    … …  
    203218void molecule::Scale(const double ** const factor)
    204219{
    205   atom *ptr = start;
    206 
    207   while (ptr->next != end) {
    208     ptr = ptr->next;
     220  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    209221    for (int j=0;j<MDSteps;j++)
    210       ptr->Trajectory.R.at(j).ScaleAll(*factor);
    211     ptr->x.ScaleAll(*factor);
     222      (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
     223    (*iter)->x.ScaleAll(*factor);
    212224  }
    213225};
    … …  
    218230void molecule::Translate(const Vector *trans)
    219231{
    220   atom *ptr = start;
    221 
    222   while (ptr->next != end) {
    223     ptr = ptr->next;
     232  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    224233    for (int j=0;j<MDSteps;j++)
    225       ptr->Trajectory.R.at(j) += (*trans);
    226     ptr->x += (*trans);
     234      (*iter)->Trajectory.R.at(j) += (*trans);
     235    (*iter)->x += (*trans);
    227236  }
    228237};
    … …  
    239248
    240249  // go through all atoms
    241   ActOnAllVectors( &Vector::SubtractVector, *trans);
     250  ActOnAllVectors( &Vector::AddVector, *trans);
    242251  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    243252
    244   Free(&M);
    245   Free(&Minv);
     253  delete[](M);
     254  delete[](Minv);
    246255};
    247256
    … …  
    252261void molecule::Mirror(const Vector *n)
    253262{
     263  OBSERVE;
    254264  Plane p(*n,0);
    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);
     265  BOOST_FOREACH( atom* iter, atoms ){
     266    (*iter->node) = p.mirrorVector(*iter->node);
    261267  }
    262268};
    … …  
    267273void molecule::DeterminePeriodicCenter(Vector &center)
    268274{
    269   atom *Walker = start;
    270275  double * const cell_size = World::getInstance().getDomain();
    271276  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    272   double *inversematrix = InverseMatrix(cell_size);
     277  double *inversematrix = InverseMatrix(matrix);
    273278  double tmp;
    274279  bool flag;
    … …  
    278283    Center.Zero();
    279284    flag = true;
    280     while (Walker->next != end) {
    281       Walker = Walker->next;
     285    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    282286#ifdef ADDHYDROGEN
    283       if (Walker->type->Z != 1) {
     287      if ((*iter)->type->Z != 1) {
    284288#endif
    285         Testvector = Walker->x;
     289        Testvector = (*iter)->x;
    286290        Testvector.MatrixMultiplication(inversematrix);
    287291        Translationvector.Zero();
    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
     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
    290294            for (int j=0;j<NDIM;j++) {
    291               tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
     295              tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
    292296              if ((fabs(tmp)) > BondDistance) {
    293297                flag = false;
    294                 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
     298                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << (*iter)->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    295299                if (tmp > 0)
    296300                  Translationvector[j] -= 1.;
    … …  
    306310#ifdef ADDHYDROGEN
    307311        // now also change all hydrogens
    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;
     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;
    311315            Testvector.MatrixMultiplication(inversematrix);
    312316            Testvector += Translationvector;
    … …  
    320324    }
    321325  } while (!flag);
    322   Free(&matrix);
    323   Free(&inversematrix);
    324 
    325   Center.Scale(1./(double)AtomCount);
     326  delete[](matrix);
     327  delete[](inversematrix);
     328
     329  Center.Scale(1./static_cast<double>(getAtomCount()));
    326330};
    327331
    … …  
    333337void molecule::PrincipalAxisSystem(bool DoRotate)
    334338{
    335   atom *ptr = start;  // start at first in list
    336339  double InertiaTensor[NDIM*NDIM];
    337340  Vector *CenterOfGravity = DetermineCenterOfGravity();
    … …  
    344347
    345348  // sum up inertia tensor
    346   while (ptr->next != end) {
    347     ptr = ptr->next;
    348     Vector x = ptr->x;
     349  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     350    Vector x = (*iter)->x;
    349351    //x.SubtractVector(CenterOfGravity);
    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]);
     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]);
    359361  }
    360362  // print InertiaTensor for debugging
    … …  
    394396
    395397    // sum up inertia tensor
    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]);
     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]);
    410409    }
    411410    // print InertiaTensor for debugging
    … …  
    431430void molecule::Align(Vector *n)
    432431{
    433   atom *ptr = start;
    434432  double alpha, tmp;
    435433  Vector z_axis;
    … …  
    442440  alpha = atan(-n->at(0)/n->at(2));
    443441  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    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];
     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];
    449446    for (int j=0;j<MDSteps;j++) {
    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];
     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];
    453450    }
    454451  }
    … …  
    460457
    461458  // rotate on z-y plane
    462   ptr = start;
    463459  alpha = atan(-n->at(1)/n->at(2));
    464460  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    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];
     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];
    470465    for (int j=0;j<MDSteps;j++) {
    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];
     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];
    474469    }
    475470  }
    … …  
    495490  Vector a,b,c,d;
    496491  struct lsq_params *par = (struct lsq_params *)params;
    497   atom *ptr = par->mol->start;
    498492
    499493  // initialize vectors
    … …  
    505499  b[2] = gsl_vector_get(x,5);
    506500  // go through all atoms
    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;
     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;
    511504      t = c.ScalarProduct(b);           // get direction parameter
    512505      d = t*b;       // and create vector
Note: See TracChangeset for help on using the changeset viewer.