Changeset 6c438f for src/molecule.cpp


Ignore:
Timestamp:
Aug 28, 2010, 3:21:11 AM (15 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, Candidate_v1.7.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
a6e6b5, f8982c
Parents:
2ad482 (diff), fd4905 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Frederik Heber <heber@…> (08/28/10 03:17:48)
git-committer:
Frederik Heber <heber@…> (08/28/10 03:21:11)
Message:

Merge branch 'StructureRefactoring' into Shapes

Conflicts:

src/Box.cpp
src/Box.hpp
src/Descriptors/AtomShapeDescriptor.cpp
src/Descriptors/AtomShapeDescriptor.hpp
src/Descriptors/AtomShapeDescriptor_impl.hpp
src/LinearAlgebra/Line.cpp
src/LinearAlgebra/Line.hpp
src/LinearAlgebra/Matrix.cpp
src/LinearAlgebra/Matrix.hpp
src/Makefile.am
src/Shapes/BaseShapes.cpp
src/Shapes/BaseShapes_impl.hpp
src/Shapes/Shape.cpp
src/Shapes/Shape.hpp
src/Shapes/ShapeOps_impl.hpp
src/Shapes/Shape_impl.hpp
src/unittests/ShapeUnittest.cpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule.cpp

    r2ad482 r6c438f  
     1/*
     2 * Project: MoleCuilder
     3 * Description: creates and alters molecular systems
     4 * Copyright (C)  2010 University of Bonn. All rights reserved.
     5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
     6 */
     7
    18/** \file molecules.cpp
    29 *
    … …  
    512 */
    613
     14// include config.h
    715#ifdef HAVE_CONFIG_H
    816#include <config.h>
    … …  
    2432#include "element.hpp"
    2533#include "graph.hpp"
    26 #include "helpers.hpp"
     34#include "Helpers/helpers.hpp"
    2735#include "leastsquaremin.hpp"
    2836#include "linkedcell.hpp"
    2937#include "lists.hpp"
    30 #include "log.hpp"
     38#include "Helpers/Log.hpp"
    3139#include "molecule.hpp"
    3240
    … …  
    3442#include "stackclass.hpp"
    3543#include "tesselation.hpp"
    36 #include "vector.hpp"
    37 #include "Matrix.hpp"
     44#include "LinearAlgebra/Vector.hpp"
     45#include "LinearAlgebra/Matrix.hpp"
    3846#include "World.hpp"
    3947#include "Box.hpp"
    40 #include "Plane.hpp"
     48#include "LinearAlgebra/Plane.hpp"
    4149#include "Exceptions/LinearDependenceException.hpp"
    4250
    … …  
    7583
    7684// getter and setter
    77 const std::string molecule::getName(){
     85const std::string molecule::getName() const{
    7886  return std::string(name);
    7987}
    … …  
    8997}
    9098
    91 moleculeId_t molecule::getId(){
     99bool molecule::changeId(moleculeId_t newId){
     100  // first we move ourselves in the world
     101  // the world lets us know if that succeeded
     102  if(World::getInstance().changeMoleculeId(id,newId,this)){
     103    id = newId;
     104    return true;
     105  }
     106  else{
     107    return false;
     108  }
     109}
     110
     111
     112moleculeId_t molecule::getId() const {
    92113  return id;
    93114}
    … …  
    97118}
    98119
    99 const Formula &molecule::getFormula(){
     120const Formula &molecule::getFormula() const {
    100121  return formula;
    101122}
    102123
    103 unsigned int molecule::getElementCount(){
     124unsigned int molecule::getElementCount() const{
    104125  return formula.getElementCount();
    105126}
    … …  
    151172molecule::const_iterator molecule::erase( const_iterator loc )
    152173{
     174  OBSERVE;
    153175  molecule::const_iterator iter = loc;
    154176  iter--;
    … …  
    156178  atomIds.erase( atom->getId() );
    157179  atoms.remove( atom );
     180  formula-=atom->getType();
    158181  atom->removeFromMolecule();
    159182  return iter;
    … …  
    162185molecule::const_iterator molecule::erase( atom * key )
    163186{
     187  OBSERVE;
    164188  molecule::const_iterator iter = find(key);
    165189  if (iter != end()){
    166190    atomIds.erase( key->getId() );
    167191    atoms.remove( key );
     192    formula-=key->getType();
    168193    key->removeFromMolecule();
    169194  }
    … …  
    183208pair<molecule::iterator,bool> molecule::insert ( atom * const key )
    184209{
     210  OBSERVE;
    185211  pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
    186212  if (res.second) { // push atom if went well
    187213    atoms.push_back(key);
     214    formula+=key->getType();
    188215    return pair<iterator,bool>(molecule::iterator(--end()),res.second);
    189216  } else {
    … …  
    206233  if (pointer != NULL) {
    207234    pointer->sort = &pointer->nr;
    208     if (pointer->type != NULL) {
    209       formula += pointer->type;
    210       if (pointer->type->Z != 1)
     235    if (pointer->getType() != NULL) {
     236      if (pointer->getType()->getAtomicNumber() != 1)
    211237        NoNonHydrogen++;
    212238      if(pointer->getName() == "Unknown"){
    213239        stringstream sstr;
    214         sstr << pointer->type->symbol << pointer->nr+1;
     240        sstr << pointer->getType()->getSymbol() << pointer->nr+1;
    215241        pointer->setName(sstr.str());
    216242      }
    … …  
    236262    walker->nr = last_atom++;  // increase number within molecule
    237263    insert(walker);
    238     if ((pointer->type != NULL) && (pointer->type->Z != 1))
     264    if ((pointer->getType() != NULL) && (pointer->getType()->getAtomicNumber() != 1))
    239265      NoNonHydrogen++;
     266    walker->setMolecule(this);
    240267    retval=walker;
    241268  }
    … …  
    294321//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    295322  // create vector in direction of bond
    296   InBondvector = TopReplacement->x - TopOrigin->x;
     323  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
    297324  bondlength = InBondvector.Norm();
    298325
    … …  
    306333    Orthovector1.Zero();
    307334    for (int i=NDIM;i--;) {
    308       l = TopReplacement->x[i] - TopOrigin->x[i];
     335      l = TopReplacement->at(i) - TopOrigin->at(i);
    309336      if (fabs(l) > BondDistance) { // is component greater than bond distance
    310337        Orthovector1[i] = (l < 0) ? -1. : +1.;
    … …  
    321348  InBondvector.Normalize();
    322349  // get typical bond length and store as scale factor for later
    323   ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
    324   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
     350  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
     351  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->BondDegree-1);
    325352  if (BondRescale == -1) {
    326353    DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
    … …  
    336363    case 1:
    337364      FirstOtherAtom = World::getInstance().createAtom();    // new atom
    338       FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
    339       FirstOtherAtom->v = TopReplacement->v; // copy velocity
     365      FirstOtherAtom->setType(1);  // element is Hydrogen
     366      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    340367      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
    341       if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
     368      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
    342369        FirstOtherAtom->father = TopReplacement;
    343370        BondRescale = bondlength;
    … …  
    346373      }
    347374      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
    348       FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
    349       FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
     375      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
    350376      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    351377//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
    … …  
    380406        // determine the plane of these two with the *origin
    381407        try {
    382           Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
     408          Orthovector1 =Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
    383409        }
    384410        catch(LinearDependenceException &excp){
    … …  
    401427      FirstOtherAtom = World::getInstance().createAtom();
    402428      SecondOtherAtom = World::getInstance().createAtom();
    403       FirstOtherAtom->type = elemente->FindElement(1);
    404       SecondOtherAtom->type = elemente->FindElement(1);
    405       FirstOtherAtom->v = TopReplacement->v; // copy velocity
     429      FirstOtherAtom->setType(1);
     430      SecondOtherAtom->setType(1);
     431      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    406432      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
    407       SecondOtherAtom->v = TopReplacement->v; // copy velocity
     433      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    408434      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
    409435      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
    410436      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    411       bondangle = TopOrigin->type->HBondAngle[1];
     437      bondangle = TopOrigin->getType()->getHBondAngle(1);
    412438      if (bondangle == -1) {
    413439        DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
    … …  
    423449//      Log() << Verbose(0) << endl;
    424450//      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
    425       FirstOtherAtom->x.Zero();
    426       SecondOtherAtom->x.Zero();
     451      FirstOtherAtom->Zero();
     452      SecondOtherAtom->Zero();
    427453      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
    428         FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
    429         SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
    430       }
    431       FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
    432       SecondOtherAtom->x *= BondRescale;
     454        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
     455        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
     456      }
     457      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
     458      SecondOtherAtom->Scale(BondRescale);
    433459      //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
    434       for(int i=NDIM;i--;) { // and make relative to origin atom
    435         FirstOtherAtom->x[i] += TopOrigin->x[i];
    436         SecondOtherAtom->x[i] += TopOrigin->x[i];
    437       }
     460      *FirstOtherAtom += TopOrigin->getPosition();
     461      *SecondOtherAtom += TopOrigin->getPosition();
    438462      // ... and add to molecule
    439463      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    … …  
    457481      SecondOtherAtom = World::getInstance().createAtom();
    458482      ThirdOtherAtom = World::getInstance().createAtom();
    459       FirstOtherAtom->type = elemente->FindElement(1);
    460       SecondOtherAtom->type = elemente->FindElement(1);
    461       ThirdOtherAtom->type = elemente->FindElement(1);
    462       FirstOtherAtom->v = TopReplacement->v; // copy velocity
     483      FirstOtherAtom->setType(1);
     484      SecondOtherAtom->setType(1);
     485      ThirdOtherAtom->setType(1);
     486      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    463487      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
    464       SecondOtherAtom->v = TopReplacement->v; // copy velocity
     488      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    465489      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
    466       ThirdOtherAtom->v = TopReplacement->v; // copy velocity
     490      ThirdOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
    467491      ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
    468492      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    … …  
    487511
    488512      // create correct coordination for the three atoms
    489       alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
     513      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
    490514      l = BondRescale;        // desired bond length
    491515      b = 2.*l*sin(alpha);    // base length of isosceles triangle
    … …  
    498522      factors[1] = f;
    499523      factors[2] = 0.;
    500       FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     524      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    501525      factors[1] = -0.5*f;
    502526      factors[2] = g;
    503       SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     527      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    504528      factors[2] = -g;
    505       ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     529      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    506530
    507531      // rescale each to correct BondDistance
    … …  
    511535
    512536      // and relative to *origin atom
    513       FirstOtherAtom->x += TopOrigin->x;
    514       SecondOtherAtom->x += TopOrigin->x;
    515       ThirdOtherAtom->x += TopOrigin->x;
     537      *FirstOtherAtom += TopOrigin->getPosition();
     538      *SecondOtherAtom += TopOrigin->getPosition();
     539      *ThirdOtherAtom += TopOrigin->getPosition();
    516540
    517541      // ... and add to molecule
    … …  
    589613    *item >> x[1];
    590614    *item >> x[2];
    591     Walker->type = elemente->FindElement(shorthand);
    592     if (Walker->type == NULL) {
     615    Walker->setType(elemente->FindElement(shorthand));
     616    if (Walker->getType() == NULL) {
    593617      DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
    594       Walker->type = elemente->FindElement(1);
     618      Walker->setType(1);
    595619    }
    596620    if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
    … …  
    599623      Walker->Trajectory.F.resize(MDSteps+10);
    600624    }
     625    Walker->setPosition(Vector(x));
    601626    for(j=NDIM;j--;) {
    602       Walker->x[j] = x[j];
    603627      Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
    604628      Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
    … …  
    619643{
    620644  molecule *copy = World::getInstance().createMolecule();
    621   atom *LeftAtom = NULL, *RightAtom = NULL;
    622645
    623646  // copy all atoms
    624   ActOnCopyWithEachAtom ( &molecule::AddCopyAtom, copy );
     647  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
    625648
    626649  // copy all bonds
    627   bond *Binder = NULL;
    628   bond *NewBond = NULL;
    629650  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
    630     for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
     651    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
    631652      if ((*BondRunner)->leftatom == *AtomRunner) {
    632         Binder = (*BondRunner);
    633 
     653        bond *Binder = (*BondRunner);
    634654        // get the pendant atoms of current bond in the copy molecule
    635         copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
    636         copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );
    637 
    638         NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
     655        atomSet::iterator leftiter=find_if(copy->atoms.begin(),copy->atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->leftatom));
     656        atomSet::iterator rightiter=find_if(copy->atoms.begin(),copy->atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->rightatom));
     657        ASSERT(leftiter!=copy->atoms.end(),"No copy of original left atom for bond copy found");
     658        ASSERT(leftiter!=copy->atoms.end(),"No copy of original right atom for bond copy found");
     659        atom *LeftAtom = *leftiter;
     660        atom *RightAtom = *rightiter;
     661
     662        bond *NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
    639663        NewBond->Cyclic = Binder->Cyclic;
    640664        if (Binder->Cyclic)
    … …  
    643667      }
    644668  // correct fathers
    645   ActOnAllAtoms( &atom::CorrectFather );
     669  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
    646670
    647671  // copy values
    … …  
    689713  ASSERT(atom2, "Second atom in bond-creation was an invalid pointer");
    690714  ASSERT(FindAtom(atom1->nr),"First atom in bond-creation was not part of molecule");
    691   ASSERT(FindAtom(atom2->nr),"Second atom in bond-creation was not parto of molecule");
     715  ASSERT(FindAtom(atom2->nr),"Second atom in bond-creation was not part of molecule");
    692716
    693717  Binder = new bond(atom1, atom2, degree, BondCount++);
    694718  atom1->RegisterBond(Binder);
    695719  atom2->RegisterBond(Binder);
    696   if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
     720  if ((atom1->getType() != NULL) && (atom1->getType()->getAtomicNumber() != 1) && (atom2->getType() != NULL) && (atom2->getType()->getAtomicNumber() != 1))
    697721    NoNonBonds++;
    698722
    … …  
    760784};
    761785
    762 /** Removes atom from molecule list and deletes it.
     786/** Removes atom from molecule list and removes all of its bonds.
    763787 * \param *pointer atom to be removed
    764788 * \return true - succeeded, false - atom not found in list
    … …  
    768792  ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
    769793  OBSERVE;
    770   formula-=pointer->type;
    771794  RemoveBonds(pointer);
    772795  erase(pointer);
    … …  
    782805  if (pointer == NULL)
    783806    return false;
    784   formula-=pointer->type;
    785807  erase(pointer);
    786808  return true;
    … …  
    793815{
    794816  for (molecule::iterator iter = begin(); !empty(); iter = begin())
    795       erase(iter);
     817    erase(*iter);
    796818  return empty();
    797819};
    … …  
    852874 * \param *out output stream
    853875 */
    854 bool molecule::Output(ofstream * const output)
    855 {
    856   int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
    857 
    858   for (int i=0;i<MAX_ELEMENTS;++i) {
    859     AtomNo[i] = 0;
    860     ElementNo[i] = 0;
    861   }
     876bool molecule::Output(ostream * const output)
     877{
    862878  if (output == NULL) {
    863879    return false;
    864880  } else {
     881    int AtomNo[MAX_ELEMENTS];
     882    memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
     883    enumeration<const element*> elementLookup = formula.enumerateElements();
    865884    *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
    866     SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
    867     int current=1;
    868     for (int i=0;i<MAX_ELEMENTS;++i) {
    869       if (ElementNo[i] == 1)
    870         ElementNo[i] = current++;
    871     }
    872     ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     885    for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
    873886    return true;
    874887  }
    … …  
    880893bool molecule::OutputTrajectories(ofstream * const output)
    881894{
    882   int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
    883 
    884895  if (output == NULL) {
    885896    return false;
    … …  
    891902        *output << "# ====== MD step " << step << " =========" << endl;
    892903      }
    893       for (int i=0;i<MAX_ELEMENTS;++i) {
    894         AtomNo[i] = 0;
    895         ElementNo[i] = 0;
    896       }
    897       SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
    898       int current=1;
    899       for (int i=0;i<MAX_ELEMENTS;++i) {
    900         if (ElementNo[i] == 1)
    901           ElementNo[i] = current++;
    902       }
    903       ActOnAllAtoms( &atom::OutputTrajectory, output, (const int *)ElementNo, AtomNo, (const int)step );
     904      int AtomNo[MAX_ELEMENTS];
     905      memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
     906      enumeration<const element*> elementLookup = formula.enumerateElements();
     907      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectory,_1,output,elementLookup, AtomNo, (const int)step));
    904908    }
    905909    return true;
    … …  
    913917{
    914918  DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
    915   ActOnAllAtoms (&atom::OutputBondOfAtom );
     919  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
    916920  DoLog(0) && (Log() << Verbose(0) << endl);
    917921};
    … …  
    936940    for (int step=0;step<MDSteps;step++) {
    937941      *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
    938       ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     942      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
    939943    }
    940944    return true;
    … …  
    953957    now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    954958    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
    955     ActOnAllAtoms( &atom::OutputXYZLine, output );
     959    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
    956960    return true;
    957961  } else
    … …  
    969973  for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    970974    (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
    971     if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
     975    if ((*iter)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it
    972976      NoNonHydrogen++;
    973977    stringstream sstr;
    974     sstr << (*iter)->type->symbol << (*iter)->nr+1;
     978    sstr << (*iter)->getType()->getSymbol() << (*iter)->nr+1;
    975979    (*iter)->setName(sstr.str());
    976980    DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
    … …  
    978982  }
    979983  return res;
    980 };
    981 
    982 /** Determines whether two molecules actually contain the same atoms and coordination.
    983  * \param *out output stream for debugging
    984  * \param *OtherMolecule the molecule to compare this one to
    985  * \param threshold upper limit of difference when comparing the coordination.
    986  * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
    987  */
    988 int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)
    989 {
    990   int flag;
    991   double *Distances = NULL, *OtherDistances = NULL;
    992   Vector CenterOfGravity, OtherCenterOfGravity;
    993   size_t *PermMap = NULL, *OtherPermMap = NULL;
    994   int *PermutationMap = NULL;
    995   bool result = true; // status of comparison
    996 
    997   DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
    998   /// first count both their atoms and elements and update lists thereby ...
    999   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
    1000 
    1001   /// ... and compare:
    1002   /// -# AtomCount
    1003   if (result) {
    1004     if (getAtomCount() != OtherMolecule->getAtomCount()) {
    1005       DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
    1006       result = false;
    1007     } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
    1008   }
    1009   /// -# Formula
    1010   if (result) {
    1011     if (formula != OtherMolecule->formula) {
    1012       DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);
    1013       result = false;
    1014     } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;
    1015   }
    1016   /// then determine and compare center of gravity for each molecule ...
    1017   if (result) {
    1018     DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
    1019     DeterminePeriodicCenter(CenterOfGravity);
    1020     OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    1021     DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);
    1022     DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);
    1023     if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {
    1024       DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
    1025       result = false;
    1026     }
    1027   }
    1028 
    1029   /// ... then make a list with the euclidian distance to this center for each atom of both molecules
    1030   if (result) {
    1031     DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
    1032     Distances = new double[getAtomCount()];
    1033     OtherDistances = new double[getAtomCount()];
    1034     SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
    1035     SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
    1036     for(int i=0;i<getAtomCount();i++) {
    1037       Distances[i] = 0.;
    1038       OtherDistances[i] = 0.;
    1039     }
    1040 
    1041     /// ... sort each list (using heapsort (o(N log N)) from GSL)
    1042     DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
    1043     PermMap = new size_t[getAtomCount()];
    1044     OtherPermMap = new size_t[getAtomCount()];
    1045     for(int i=0;i<getAtomCount();i++) {
    1046       PermMap[i] = 0;
    1047       OtherPermMap[i] = 0;
    1048     }
    1049     gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
    1050     gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
    1051     PermutationMap = new int[getAtomCount()];
    1052     for(int i=0;i<getAtomCount();i++)
    1053       PermutationMap[i] = 0;
    1054     DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
    1055     for(int i=getAtomCount();i--;)
    1056       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    1057 
    1058     /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
    1059     DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
    1060     flag = 0;
    1061     for (int i=0;i<getAtomCount();i++) {
    1062       DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl);
    1063       if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
    1064         flag = 1;
    1065     }
    1066 
    1067     // free memory
    1068     delete[](PermMap);
    1069     delete[](OtherPermMap);
    1070     delete[](Distances);
    1071     delete[](OtherDistances);
    1072     if (flag) { // if not equal
    1073       delete[](PermutationMap);
    1074       result = false;
    1075     }
    1076   }
    1077   /// return pointer to map if all distances were below \a threshold
    1078   DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
    1079   if (result) {
    1080     DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
    1081     return PermutationMap;
    1082   } else {
    1083     DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
    1084     return NULL;
    1085   }
    1086984};
    1087985
    … …  
    11421040    *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
    11431041  for (int step=startstep;step < endstep; step++) { // loop over all time steps
    1144     temperature = 0.;
    1145     ActOnAllAtoms( &TrajectoryParticle::AddKineticToTemperature, &temperature, step);
     1042    temperature = atoms.totalTemperatureAtStep(step);
    11461043    *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
    11471044  }
    11481045  return true;
    1149 };
    1150 
    1151 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
    1152 {
    1153   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    1154     array[((*iter)->*index)] = (*iter);
    1155   }
    11561046};
    11571047
Note: See TracChangeset for help on using the changeset viewer.