Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_geometry.cpp

    re138de ra67d19  
    1515#include "memoryallocator.hpp"
    1616#include "molecule.hpp"
     17#include "World.hpp"
    1718
    1819/************************************* Functions for class molecule *********************************/
    … …  
    2627  bool status = true;
    2728  const Vector *Center = DetermineCenterOfAll();
     29  double * const cell_size = World::get()->cell_size;
    2830  double *M = ReturnFullMatrixforSymmetric(cell_size);
    2931  double *Minv = InverseMatrix(M);
    … …  
    3335  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    3436
    35   delete(M);
    36   delete(Minv);
     37  Free(&M);
     38  Free(&Minv);
    3739  delete(Center);
    3840  return status;
    … …  
    4648{
    4749  bool status = true;
     50  double * const cell_size = World::get()->cell_size;
    4851  double *M = ReturnFullMatrixforSymmetric(cell_size);
    4952  double *Minv = InverseMatrix(M);
    … …  
    5255  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    5356
    54   delete(M);
    55   delete(Minv);
     57  Free(&M);
     58  Free(&Minv);
    5659  return status;
    5760};
    … …  
    101104{
    102105  int Num = 0;
    103   atom *ptr = start->next;  // start at first in list
     106  atom *ptr = start;  // start at first in list
    104107
    105108  Center.Zero();
    106109
    107   if (ptr != end) {   //list not empty?
     110  if (ptr->next != end) {   //list not empty?
    108111    while (ptr->next != end) {  // continue with second if present
    109112      ptr = ptr->next;
    … …  
    226229void molecule::TranslatePeriodically(const Vector *trans)
    227230{
     231  double * const cell_size = World::get()->cell_size;
    228232  double *M = ReturnFullMatrixforSymmetric(cell_size);
    229233  double *Minv = InverseMatrix(M);
    … …  
    233237  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    234238
    235   delete(M);
    236   delete(Minv);
     239  Free(&M);
     240  Free(&Minv);
    237241};
    238242
    … …  
    252256{
    253257  atom *Walker = start;
     258  double * const cell_size = World::get()->cell_size;
    254259  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     260  double *inversematrix = InverseMatrix(cell_size);
    255261  double tmp;
    256262  bool flag;
    … …  
    266272#endif
    267273        Testvector.CopyVector(&Walker->x);
    268         Testvector.InverseMatrixMultiplication(matrix);
     274        Testvector.MatrixMultiplication(inversematrix);
    269275        Translationvector.Zero();
    270276        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
    … …  
    274280              if ((fabs(tmp)) > BondDistance) {
    275281                flag = false;
    276                 Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
     282                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    277283                if (tmp > 0)
    278284                  Translationvector.x[j] -= 1.;
    … …  
    285291        Testvector.MatrixMultiplication(matrix);
    286292        Center.AddVector(&Testvector);
    287         Log() << Verbose(1) << "vector is: ";
     293        DoLog(1) && (Log() << Verbose(1) << "vector is: ");
    288294        Testvector.Output();
    289         Log() << Verbose(0) << endl;
     295        DoLog(0) && (Log() << Verbose(0) << endl);
    290296#ifdef ADDHYDROGEN
    291297        // now also change all hydrogens
    … …  
    293299          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
    294300            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
    295             Testvector.InverseMatrixMultiplication(matrix);
     301            Testvector.MatrixMultiplication(inversematrix);
    296302            Testvector.AddVector(&Translationvector);
    297303            Testvector.MatrixMultiplication(matrix);
    298304            Center.AddVector(&Testvector);
    299             Log() << Verbose(1) << "Hydrogen vector is: ";
     305            DoLog(1) && (Log() << Verbose(1) << "Hydrogen vector is: ");
    300306            Testvector.Output();
    301             Log() << Verbose(0) << endl;
     307            DoLog(0) && (Log() << Verbose(0) << endl);
    302308          }
    303309        }
    … …  
    307313  } while (!flag);
    308314  Free(&matrix);
     315  Free(&inversematrix);
     316
    309317  Center.Scale(1./(double)AtomCount);
    310318};
    … …  
    344352  }
    345353  // print InertiaTensor for debugging
    346   Log() << Verbose(0) << "The inertia tensor is:" << endl;
     354  DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
    347355  for(int i=0;i<NDIM;i++) {
    348356    for(int j=0;j<NDIM;j++)
    349       Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ";
    350     Log() << Verbose(0) << endl;
    351   }
    352   Log() << Verbose(0) << endl;
     357      DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
     358    DoLog(0) && (Log() << Verbose(0) << endl);
     359  }
     360  DoLog(0) && (Log() << Verbose(0) << endl);
    353361
    354362  // diagonalize to determine principal axis system
    … …  
    362370
    363371  for(int i=0;i<NDIM;i++) {
    364     Log() << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
    365     Log() << Verbose(0) << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
     372    DoLog(1) && (Log() << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i));
     373    DoLog(0) && (Log() << Verbose(0) << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl);
    366374  }
    367375
    368376  // check whether we rotate or not
    369377  if (DoRotate) {
    370     Log() << Verbose(1) << "Transforming molecule into PAS ... ";
     378    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    371379    // the eigenvectors specify the transformation matrix
    372380    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
    373     Log() << Verbose(0) << "done." << endl;
     381    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    374382
    375383    // summing anew for debugging (resulting matrix has to be diagonal!)
    … …  
    396404    }
    397405    // print InertiaTensor for debugging
    398     Log() << Verbose(0) << "The inertia tensor is:" << endl;
     406    DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
    399407    for(int i=0;i<NDIM;i++) {
    400408      for(int j=0;j<NDIM;j++)
    401         Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ";
    402       Log() << Verbose(0) << endl;
    403     }
    404     Log() << Verbose(0) << endl;
     409        DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
     410      DoLog(0) && (Log() << Verbose(0) << endl);
     411    }
     412    DoLog(0) && (Log() << Verbose(0) << endl);
    405413  }
    406414
    … …  
    425433
    426434  // rotate on z-x plane
    427   Log() << Verbose(0) << "Begin of Aligning all atoms." << endl;
     435  DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
    428436  alpha = atan(-n->x[0]/n->x[2]);
    429   Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ";
     437  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    430438  while (ptr->next != end) {
    431439    ptr = ptr->next;
    … …  
    443451  n->x[0] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    444452  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    445   Log() << Verbose(1) << "alignment vector after first rotation: ";
     453  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: ");
    446454  n->Output();
    447   Log() << Verbose(0) << endl;
     455  DoLog(0) && (Log() << Verbose(0) << endl);
    448456
    449457  // rotate on z-y plane
    450458  ptr = start;
    451459  alpha = atan(-n->x[1]/n->x[2]);
    452   Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ";
     460  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    453461  while (ptr->next != end) {
    454462    ptr = ptr->next;
    … …  
    467475  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    468476
    469   Log() << Verbose(1) << "alignment vector after second rotation: ";
     477  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: ");
    470478  n->Output();
    471   Log() << Verbose(1) << endl;
    472   Log() << Verbose(0) << "End of Aligning all atoms." << endl;
     479  DoLog(1) && (Log() << Verbose(1) << endl);
     480  DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
    473481};
    474482
Note: See TracChangeset for help on using the changeset viewer.