Changeset fcd7b6 for src


Ignore:
Timestamp:
Oct 7, 2009, 1:11:28 PM (16 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, 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:
681a8a
Parents:
e9f8f9
git-author:
Frederik Heber <heber@…> (10/07/09 12:14:15)
git-committer:
Frederik Heber <heber@…> (10/07/09 13:11:28)
Message:

In molecule::OutputTrajectories() ActOnAllAtoms() with new function atom::OutputTrajectory() is used.

For this to work, I had to change the Trajectory struct that was so far included in molecule.hpp to be incorporated directly into the class atom.
NOTE: This incorporation is incomplete and a ticket (#34) has been filed to remind of this issue.
However, the trajectory is better suited to reside in atom anyway and was probably just put in molecule due to memory prejudices against STL vector<>.
Functions in molecule.cpp, config.cpp, molecule_geometry.cpp and molecule_dynamics.cpp were adapted (changed from Trajectories[atom *] to atom *->Trajectory).
And the atom pointer in the Trajectory structure was removed.

Location:
src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • src/atom.cpp

    re9f8f9 rfcd7b6  
    6464  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
    6565  Free<char>(&Name, "atom::~atom: *Name");
     66  Trajectory.R.clear();
     67  Trajectory.U.clear();
     68  Trajectory.F.clear();
    6669};
    6770
    … …  
    118121 * \param *comment commentary after '#' sign
    119122 */
    120 bool atom::Output(int ElementNo, int AtomNo, ofstream *out, const char *comment) const
     123bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
    121124{
    122125  if (out != NULL) {
    … …  
    134137    return false;
    135138};
    136 bool atom::Output(int *ElementNo, int *AtomNo, ofstream *out, const char *comment)
     139bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
    137140{
    138141  AtomNo[type->Z]++;  // increment number
    … …  
    164167};
    165168
     169/** Output of a single atom as one lin in xyz file.
     170 * \param *out stream to output to
     171 */
     172bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
     173{
     174  AtomNo[type->Z]++;
     175  if (out != NULL) {
     176    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
     177    *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
     178    *out << "\t" << FixedIon;
     179    if (Trajectory.U.at(step).Norm() > MYEPSILON)
     180      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
     181    if (Trajectory.F.at(step).Norm() > MYEPSILON)
     182      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t";
     183    *out << "\t# Number in molecule " << nr << endl;
     184    return true;
     185  } else
     186    return false;
     187};
     188
    166189ostream & operator << (ostream &ost, const atom &a)
    167190{

  • src/atom.hpp

    re9f8f9 rfcd7b6  
    1818
    1919#include <iostream>
     20#include <vector>
    2021
    2122#include "element.hpp"
    … …  
    2829class atom : public TesselPoint {
    2930  public:
    30     Vector x;       //!< coordinate array of atom, giving position within cell
    31     Vector v;       //!< velocity array of atom
     31    struct
     32    {
     33      vector<Vector> R;  //!< position vector
     34      vector<Vector> U;  //!< velocity vector
     35      vector<Vector> F;  //!< last force vector
     36    } Trajectory;
     37
     38    Vector x;       //!< coordinate vector of atom, giving last position within cell
     39    Vector v;       //!< velocity vector of atom, giving last velocity within cell
     40    Vector F;       //!< Force vector of atom, giving last force within cell
    3241    element *type;  //!< pointing to element
    3342    atom *previous; //!< previous atom in molecule list
    … …  
    5160  virtual ~atom();
    5261
    53   bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
    54   bool Output(int *ElementNo, int *AtomNo, ofstream *out, const char *comment = NULL);
     62  bool Output(ofstream *out, int ElementNo, int AtomNo, const char *comment = NULL) const;
     63  bool Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment = NULL);
    5564  bool OutputXYZLine(ofstream *out) const;
     65  bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
    5666  void EqualsFather ( atom *ptr, atom **res );
    5767  void CorrectFather();
    … …  
    6676};
    6777
     78
    6879ostream & operator << (ostream &ost, const atom &a);
    6980

  • src/boundary.cpp

    re9f8f9 rfcd7b6  
    10291029
    10301030  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
    1031   mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
     1031  //mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
    10321032//  mol->GoToFirst();
    10331033//  class TesselPoint *Runner = NULL;

  • src/boundary.hpp

    re9f8f9 rfcd7b6  
    1616
    1717#define DEBUG 1
    18 #define DoSingleStepOutput 1
     18#define DoSingleStepOutput 0
    1919#define SingleStepWidth 1
    2020

  • src/config.cpp

    re9f8f9 rfcd7b6  
    964964
    965965            // check size of vectors
    966             if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
     966            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
    967967              //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
    968               mol->Trajectories[neues].R.resize(repetition+10);
    969               mol->Trajectories[neues].U.resize(repetition+10);
    970               mol->Trajectories[neues].F.resize(repetition+10);
     968              neues->Trajectory.R.resize(repetition+10);
     969              neues->Trajectory.U.resize(repetition+10);
     970              neues->Trajectory.F.resize(repetition+10);
    971971            }
    972972
    973973            // put into trajectories list
    974974            for (int d=0;d<NDIM;d++)
    975               mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
     975              neues->Trajectory.R.at(repetition).x[d] = neues->x.x[d];
    976976
    977977            // parse velocities if present
    … …  
    983983              neues->v.x[2] = 0.;
    984984            for (int d=0;d<NDIM;d++)
    985               mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
     985              neues->Trajectory.U.at(repetition).x[d] = neues->v.x[d];
    986986
    987987            // parse forces if present
    … …  
    993993              value[2] = 0.;
    994994            for (int d=0;d<NDIM;d++)
    995               mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
     995              neues->Trajectory.F.at(repetition).x[d] = value[d];
    996996
    997997  //            cout << "Parsed position of step " << (repetition) << ": (";
    998998  //            for (int d=0;d<NDIM;d++)
    999   //              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
     999  //              cout << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
    10001000  //            cout << ")\t(";
    10011001  //            for (int d=0;d<NDIM;d++)
    1002   //              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
     1002  //              cout << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
    10031003  //            cout << ")\t(";
    10041004  //            for (int d=0;d<NDIM;d++)
    1005   //              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
     1005  //              cout << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
    10061006  //            cout << ")" << endl;
    10071007          }

  • src/molecule.cpp

    re9f8f9 rfcd7b6  
    458458      Walker->type = elemente->FindElement(1);
    459459    }
    460     if (Trajectories[Walker].R.size() <= (unsigned int)MDSteps) {
    461       Trajectories[Walker].R.resize(MDSteps+10);
    462       Trajectories[Walker].U.resize(MDSteps+10);
    463       Trajectories[Walker].F.resize(MDSteps+10);
     460    if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
     461      Walker->Trajectory.R.resize(MDSteps+10);
     462      Walker->Trajectory.U.resize(MDSteps+10);
     463      Walker->Trajectory.F.resize(MDSteps+10);
    464464    }
    465465    for(j=NDIM;j--;) {
    466466      Walker->x.x[j] = x[j];
    467       Trajectories[Walker].R.at(MDSteps-1).x[j] = x[j];
    468       Trajectories[Walker].U.at(MDSteps-1).x[j] = 0;
    469       Trajectories[Walker].F.at(MDSteps-1).x[j] = 0;
     467      Walker->Trajectory.R.at(MDSteps-1).x[j] = x[j];
     468      Walker->Trajectory.U.at(MDSteps-1).x[j] = 0;
     469      Walker->Trajectory.F.at(MDSteps-1).x[j] = 0;
    470470    }
    471471    AddAtom(Walker);  // add to molecule
    … …  
    624624  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    625625    ElementCount--;
    626   Trajectories.erase(pointer);
    627626  return remove(pointer, start, end);
    628627};
    … …  
    642641  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    643642    ElementCount--;
    644   Trajectories.erase(pointer);
    645643  unlink(pointer);
    646644  return true;
    … …  
    726724        ElementNo[i] = current++;
    727725    }
    728     ActOnAllAtoms( &atom::Output, ElementNo, AtomNo, out, (const char *) NULL ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
     726    ActOnAllAtoms( &atom::Output, out, ElementNo, AtomNo, (const char *) NULL ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
    729727    return true;
    730728  }
    … …  
    759757          ElementNo[i] = current++;
    760758      }
    761       walker = start;
    762       while (walker->next != end) { // go through every atom of this element
    763         walker = walker->next;
    764         AtomNo[walker->type->Z]++;
    765         *out << "Ion_Type" << ElementNo[walker->type->Z] << "_" << AtomNo[walker->type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
    766         *out << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2];
    767         *out << "\t" << walker->FixedIon;
    768         if (Trajectories[walker].U.at(step).Norm() > MYEPSILON)
    769           *out << "\t" << scientific << setprecision(6) << Trajectories[walker].U.at(step).x[0] << "\t" << Trajectories[walker].U.at(step).x[1] << "\t" << Trajectories[walker].U.at(step).x[2] << "\t";
    770         if (Trajectories[walker].F.at(step).Norm() > MYEPSILON)
    771           *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t";
    772         *out << "\t# Number in molecule " << walker->nr << endl;
    773       }
     759      ActOnAllAtoms( &atom::OutputTrajectory, out, ElementNo, AtomNo, step ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
    774760    }
    775761    return true;
    … …  
    829815      while (walker->next != end) { // go through every atom of this element
    830816        walker = walker->next;
    831         *out << walker->type->symbol << "\t" << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2] << endl;
     817        *out << walker->type->symbol << "\t" << walker->Trajectory.R.at(step).x[0] << "\t" << walker->Trajectory.R.at(step).x[1] << "\t" << walker->Trajectory.R.at(step).x[2] << endl;
    832818      }
    833819    }
    … …  
    12371223      Walker = Walker->next;
    12381224      for (int i=NDIM;i--;)
    1239         temperature += Walker->type->mass * Trajectories[Walker].U.at(step).x[i]* Trajectories[Walker].U.at(step).x[i];
     1225        temperature += Walker->type->mass * Walker->Trajectory.U.at(step).x[i]* Walker->Trajectory.U.at(step).x[i];
    12401226    }
    12411227    *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;

  • src/molecule.hpp

    re9f8f9 rfcd7b6  
    6565{
    6666  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
    67 };
    68 
    69 struct Trajectory
    70 {
    71   vector<Vector> R;  //!< position vector
    72   vector<Vector> U;  //!< velocity vector
    73   vector<Vector> F;  //!< last force vector
    74   atom *ptr;         //!< pointer to atom whose trajectory we contain
    7567};
    7668
    … …  
    10193    bond *last;         //!< end of bond list
    10294    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
    103     map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
    10495    int MDSteps;        //!< The number of MD steps in Trajectories
    10596    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has

  • src/molecule_dynamics.cpp

    re9f8f9 rfcd7b6  
    4646    // first term: distance to target
    4747    Runner = PermutationMap[Walker->nr];   // find target point
    48     tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
     48    tmp = (Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)));
    4949    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    5050    result += constants[0] * tmp;
    … …  
    5959      // determine normalized trajectories direction vector (n1, n2)
    6060      Sprinter = PermutationMap[Walker->nr];   // find first target point
    61       trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
    62       trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
     61      trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));
     62      trajectory1.SubtractVector(&Walker->Trajectory.R.at(startstep));
    6363      trajectory1.Normalize();
    6464      Norm1 = trajectory1.Norm();
    6565      Sprinter = PermutationMap[Runner->nr];   // find second target point
    66       trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
    67       trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
     66      trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));
     67      trajectory2.SubtractVector(&Runner->Trajectory.R.at(startstep));
    6868      trajectory2.Normalize();
    6969      Norm2 = trajectory1.Norm();
    7070      // check whether either is zero()
    7171      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    72         tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     72        tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
    7373      } else if (Norm1 < MYEPSILON) {
    7474        Sprinter = PermutationMap[Walker->nr];   // find first target point
    75         trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
    76         trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
     75        trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy first offset
     76        trajectory1.SubtractVector(&Runner->Trajectory.R.at(startstep));  // subtract second offset
    7777        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    7878        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
    … …  
    8080      } else if (Norm2 < MYEPSILON) {
    8181        Sprinter = PermutationMap[Runner->nr];   // find second target point
    82         trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
    83         trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
     82        trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy second offset
     83        trajectory2.SubtractVector(&Walker->Trajectory.R.at(startstep));  // subtract first offset
    8484        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    8585        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
    … …  
    9090//        *out << " and ";
    9191//        *out << trajectory2;
    92         tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     92        tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
    9393//        *out << " with distance " << tmp << "." << endl;
    9494      } else { // determine distance by finding minimum distance
    … …  
    109109          gsl_matrix_set(A, 1, i, trajectory2.x[i]);
    110110          gsl_matrix_set(A, 2, i, normal.x[i]);
    111           gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
     111          gsl_vector_set(x,i, (Walker->Trajectory.R.at(startstep).x[i] - Runner->Trajectory.R.at(startstep).x[i]));
    112112        }
    113113        // solve the linear system by Householder transformations
    … …  
    117117//        *out << " with distance " << tmp << "." << endl;
    118118        // test whether we really have the intersection (by checking on c_1 and c_2)
    119         TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
     119        TestVector.CopyVector(&Runner->Trajectory.R.at(startstep));
    120120        trajectory2.Scale(gsl_vector_get(x,1));
    121121        TestVector.AddVector(&trajectory2);
    122122        normal.Scale(gsl_vector_get(x,2));
    123123        TestVector.AddVector(&normal);
    124         TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
     124        TestVector.SubtractVector(&Walker->Trajectory.R.at(startstep));
    125125        trajectory1.Scale(gsl_vector_get(x,0));
    126126        TestVector.SubtractVector(&trajectory1);
    … …  
    148148        Sprinter = PermutationMap[Walker->nr];
    149149//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
    150 //        *out << Trajectories[Sprinter].R.at(endstep);
     150//        *out << Sprinter->Trajectory.R.at(endstep);
    151151//        *out << ", penalting." << endl;
    152152        result += constants[2];
    … …  
    240240    while(Runner->next != end) {
    241241      Runner = Runner->next;
    242       DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
     242      DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)), Runner) );
    243243    }
    244244  }
    … …  
    398398    // set forces
    399399    for (int i=NDIM;i++;)
    400       Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
     400      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Walker->Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
    401401  }
    402402  *out << Verbose(1) << "done." << endl;
    … …  
    436436  while (Walker->next != end) {
    437437    Walker = Walker->next;
    438     if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
     438    if (Walker->Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
    439439      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
    440       Trajectories[Walker].R.resize(MaxSteps+1);
    441       Trajectories[Walker].U.resize(MaxSteps+1);
    442       Trajectories[Walker].F.resize(MaxSteps+1);
     440      Walker->Trajectory.R.resize(MaxSteps+1);
     441      Walker->Trajectory.U.resize(MaxSteps+1);
     442      Walker->Trajectory.F.resize(MaxSteps+1);
    443443    }
    444444  }
    … …  
    448448    Walker = Walker->next;
    449449    for (int n=NDIM;n--;) {
    450       Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
    451       Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
    452       Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
     450      Walker->Trajectory.R.at(MaxSteps).x[n] = Walker->Trajectory.R.at(endstep).x[n];
     451      Walker->Trajectory.U.at(MaxSteps).x[n] = Walker->Trajectory.U.at(endstep).x[n];
     452      Walker->Trajectory.F.at(MaxSteps).x[n] = Walker->Trajectory.F.at(endstep).x[n];
    453453    }
    454454  }
    … …  
    466466      Sprinter = mol->AddCopyAtom(Walker);
    467467      for (int n=NDIM;n--;) {
    468         Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
     468        Sprinter->x.x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
    469469        // add to Trajectories
    470470        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    471471        if (step < MaxSteps) {
    472           Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
    473           Trajectories[Walker].U.at(step).x[n] = 0.;
    474           Trajectories[Walker].F.at(step).x[n] = 0.;
     472          Walker->Trajectory.R.at(step).x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
     473          Walker->Trajectory.U.at(step).x[n] = 0.;
     474          Walker->Trajectory.F.at(step).x[n] = 0.;
    475475        }
    476476      }
    … …  
    554554      //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    555555      // check size of vectors
    556       if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
     556      if (walker->Trajectory.R.size() <= (unsigned int)(MDSteps)) {
    557557        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    558         Trajectories[walker].R.resize(MDSteps+10);
    559         Trajectories[walker].U.resize(MDSteps+10);
    560         Trajectories[walker].F.resize(MDSteps+10);
     558        walker->Trajectory.R.resize(MDSteps+10);
     559        walker->Trajectory.U.resize(MDSteps+10);
     560        walker->Trajectory.F.resize(MDSteps+10);
    561561      }
    562562
    563563      // Update R (and F)
    564564      for (int d=0; d<NDIM; d++) {
    565         Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
    566         Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
    567         Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
    568         Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
     565        walker->Trajectory.F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     566        walker->Trajectory.R.at(MDSteps).x[d] = walker->Trajectory.R.at(MDSteps-1).x[d];
     567        walker->Trajectory.R.at(MDSteps).x[d] += configuration.Deltat*(walker->Trajectory.U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
     568        walker->Trajectory.R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(walker->Trajectory.F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
    569569      }
    570570      // Update U
    571571      for (int d=0; d<NDIM; d++) {
    572         Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
    573         Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
     572        walker->Trajectory.U.at(MDSteps).x[d] = walker->Trajectory.U.at(MDSteps-1).x[d];
     573        walker->Trajectory.U.at(MDSteps).x[d] += configuration.Deltat * (walker->Trajectory.F.at(MDSteps).x[d]+walker->Trajectory.F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
    574574      }
    575575//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
    576576//      for (int d=0;d<NDIM;d++)
    577 //        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
     577//        out << walker->Trajectory.R.at(MDSteps).x[d] << " ";          // next step
    578578//      out << ")\t(";
    579579//      for (int d=0;d<NDIM;d++)
    580 //        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
     580//        cout << walker->Trajectory.U.at(MDSteps).x[d] << " ";          // next step
    581581//      out << ")" << endl;
    582582            // next atom
    … …  
    592592    IonMass += walker->type->mass;  // sum up total mass
    593593    for(int d=0;d<NDIM;d++) {
    594       Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
     594      Vector[d] += walker->Trajectory.U.at(MDSteps).x[d]*walker->type->mass;
    595595    }
    596596  }
    … …  
    603603    walker = walker->next;
    604604    for(int d=0;d<NDIM;d++) {
    605       Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
    606       ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
     605      walker->Trajectory.U.at(MDSteps).x[d] -= Vector[d];
     606      ActualTemp += 0.5 * walker->type->mass * walker->Trajectory.U.at(MDSteps).x[d] * walker->Trajectory.U.at(MDSteps).x[d];
    607607    }
    608608  }
    … …  
    663663          walker = walker->next;
    664664          IonMass = walker->type->mass;
    665           U = Trajectories[walker].U.at(MDSteps).x;
     665          U = walker->Trajectory.U.at(MDSteps).x;
    666666          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    667667            for (d=0; d<NDIM; d++) {
    … …  
    678678        walker = walker->next;
    679679        IonMass = walker->type->mass;
    680         U = Trajectories[walker].U.at(MDSteps).x;
    681         F = Trajectories[walker].F.at(MDSteps).x;
     680        U = walker->Trajectory.U.at(MDSteps).x;
     681        F = walker->Trajectory.F.at(MDSteps).x;
    682682        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    683683          for (d=0; d<NDIM; d++) {
    … …  
    691691        walker = walker->next;
    692692        IonMass = walker->type->mass;
    693         U = Trajectories[walker].U.at(MDSteps).x;
    694         F = Trajectories[walker].F.at(MDSteps).x;
     693        U = walker->Trajectory.U.at(MDSteps).x;
     694        F = walker->Trajectory.F.at(MDSteps).x;
    695695        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    696696          for (d=0; d<NDIM; d++) {
    … …  
    713713        IonMass = walker->type->mass;
    714714        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    715         U = Trajectories[walker].U.at(MDSteps).x;
    716         F = Trajectories[walker].F.at(MDSteps).x;
     715        U = walker->Trajectory.U.at(MDSteps).x;
     716        F = walker->Trajectory.F.at(MDSteps).x;
    717717        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    718718          // throw a dice to determine whether it gets hit by a heat bath particle
    … …  
    736736        walker = walker->next;
    737737        IonMass = walker->type->mass;
    738         U = Trajectories[walker].U.at(MDSteps).x;
    739         F = Trajectories[walker].F.at(MDSteps).x;
     738        U = walker->Trajectory.U.at(MDSteps).x;
     739        F = walker->Trajectory.F.at(MDSteps).x;
    740740        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    741741          for (d=0; d<NDIM; d++) {
    … …  
    754754        walker = walker->next;
    755755        IonMass = walker->type->mass;
    756         U = Trajectories[walker].U.at(MDSteps).x;
     756        U = walker->Trajectory.U.at(MDSteps).x;
    757757        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    758758          for (d=0; d<NDIM; d++) {
    … …  
    769769        walker = walker->next;
    770770        IonMass = walker->type->mass;
    771         U = Trajectories[walker].U.at(MDSteps).x;
     771        U = walker->Trajectory.U.at(MDSteps).x;
    772772        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    773773          for (d=0; d<NDIM; d++) {

  • src/molecule_geometry.cpp

    re9f8f9 rfcd7b6  
    200200    ptr = ptr->next;
    201201    for (int j=0;j<MDSteps;j++)
    202       Trajectories[ptr].R.at(j).Scale(factor);
     202      ptr->Trajectory.R.at(j).Scale(factor);
    203203    ptr->x.Scale(factor);
    204204  }
    … …  
    215215    ptr = ptr->next;
    216216    for (int j=0;j<MDSteps;j++)
    217       Trajectories[ptr].R.at(j).Translate(trans);
     217      ptr->Trajectory.R.at(j).Translate(trans);
    218218    ptr->x.Translate(trans);
    219219  }
    … …  
    438438    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    439439    for (int j=0;j<MDSteps;j++) {
    440       tmp = Trajectories[ptr].R.at(j).x[0];
    441       Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    442       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     440      tmp = ptr->Trajectory.R.at(j).x[0];
     441      ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
     442      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
    443443    }
    444444  }
    … …  
    461461    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    462462    for (int j=0;j<MDSteps;j++) {
    463       tmp = Trajectories[ptr].R.at(j).x[1];
    464       Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    465       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     463      tmp = ptr->Trajectory.R.at(j).x[1];
     464      ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
     465      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
    466466    }
    467467  }
Note: See TracChangeset for help on using the changeset viewer.