Changeset 6b919f8


Ignore:
Timestamp:
Oct 20, 2009, 8:55:17 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, 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:
831a14
Parents:
c75de1
Message:

Refactored atom.cpp into multiple files.

After the class atom was refactored into multiple base classes that are inherited, these base classes are also all put into separate files. This is basically a preparatory step for the like-wise refactoring of class molecule into inherited base classes and splitting up (that is there done already). Finally, we will also separate the relations, i.e. not have "atom.hpp" included everywhere and use class atom, but rather the subclasses such as TrajectoryParticle and its header files only.
Signed-off-by: Frederik Heber <heber@…>

Location:
src
Files:
16 added
5 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rc75de1 r6b919f8  
    1 SOURCE = analysis_correlation.cpp atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp graph.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
    2 HEADER = analysis_correlation.hpp atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp graph.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
     1ATOMSOURCE = atom.cpp atom_atominfo.cpp atom_bondedparticle.cpp atom_bondedparticleinfo.cpp atom_graphnode.cpp atom_graphnodeinfo.cpp atom_particleinfo.cpp atom_trajectoryparticle.cpp atom_trajectoryparticleinfo.cpp
     2ATOMHEADER = atom.hpp atom_atominfo.hpp atom_bondedparticle.hpp atom_bondedparticleinfo.hpp atom_graphnode.hpp atom_graphnodeinfo.hpp atom_particleinfo.hpp atom_trajectoryparticle.hpp atom_trajectoryparticleinfo.hpp
     3
     4SOURCE = analysis_correlation.cpp ${ATOMSOURCE} bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp graph.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
     5HEADER = analysis_correlation.hpp ${ATOMHEADER} bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp graph.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
    36
    47BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)

  • src/atom.cpp

    rc75de1 r6b919f8  
    1515
    1616/************************************* Functions for class atom *************************************/
    17 
    18 /** Constructor of class AtomInfo.
    19  */
    20 AtomInfo::AtomInfo() : type(NULL) {};
    21 
    22 /** Destructor of class AtomInfo.
    23  */
    24 AtomInfo::~AtomInfo()
    25 {
    26 };
    27 
    28 /** Constructor of class TrajectoryParticleInfo.
    29  */
    30 TrajectoryParticleInfo::TrajectoryParticleInfo() : FixedIon(0) {};
    31 
    32 /** Destructor of class TrajectoryParticleInfo.
    33  */
    34 TrajectoryParticleInfo::~TrajectoryParticleInfo()
    35 {
    36   Trajectory.R.clear();
    37   Trajectory.U.clear();
    38   Trajectory.F.clear();
    39 };
    40 
    41 /** Constructor of class TrajectoryParticle.
    42  */
    43 TrajectoryParticle::TrajectoryParticle()
    44 {
    45 };
    46 
    47 /** Destructor of class TrajectoryParticle.
    48  */
    49 TrajectoryParticle::~TrajectoryParticle()
    50 {
    51 };
    52 
    53 /** Constructor of class GraphNodeInfo.
    54  */
    55 GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
    56 
    57 /** Destructor of class GraphNodeInfo.
    58  */
    59 GraphNodeInfo::~GraphNodeInfo()
    60 {
    61   Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
    62 };
    63 
    64 /** Constructor of class GraphNode.
    65  */
    66 GraphNode::GraphNode()
    67 {
    68 };
    69 
    70 /** Destructor of class GraphNode.
    71  */
    72 GraphNode::~GraphNode()
    73 {
    74 };
    75 
    76 /** Constructor of class BondedParticleInfo.
    77  */
    78 BondedParticleInfo::BondedParticleInfo() : AdaptiveOrder(0), MaxOrder(false) {};
    79 
    80 /** Destructor of class BondedParticleInfo.
    81  */
    82 BondedParticleInfo::~BondedParticleInfo()
    83 {
    84 };
    85 
    86 /** Constructor of class BondedParticle.
    87  */
    88 BondedParticle::BondedParticle(){};
    89 
    90 /** Destructor of class BondedParticle.
    91  */
    92 BondedParticle::~BondedParticle()
    93 {
    94   BondList::const_iterator Runner;
    95   while (!ListOfBonds.empty()) {
    96     Runner = ListOfBonds.begin();
    97     removewithoutcheck(*Runner);
    98   }
    99 };
    10017
    10118/** Constructor of class atom.
    … …  
    360277};
    361278
    362 /** Output graph info of this atom.
    363  * \param *out output stream
    364  */
    365 void GraphNode::OutputGraphInfo(ofstream *out) const
    366 {
    367   *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
    368   OutputComponentNumber(out);
    369   *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
    370 };
    371 
    372 /** Output a list of flags, stating whether the bond was visited or not.
    373  * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
    374  * \param *out output stream for debugging
    375  */
    376 void GraphNode::OutputComponentNumber(ofstream *out) const
    377 {
    378   if (ComponentNr != NULL) {
    379     for (int i=0; ComponentNr[i] != -1; i++)
    380       *out << ComponentNr[i] << " ";
    381   }
    382 };
    383 
    384 /** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
    385  * \param *file output stream
    386  */
    387 void BondedParticle::OutputOrder(ofstream *file)
    388 {
    389   *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
    390   //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
    391 };
    392 
    393 /** Prints all bonds of this atom with total degree.
    394  * \param *out stream to output to
    395  * \return true - \a *out present, false - \a *out is NULL
    396  */
    397 bool BondedParticle::OutputBondOfAtom(ofstream *out) const
    398 {
    399   if (out != NULL) {
    400     *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
    401     int TotalDegree = 0;
    402     for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
    403       *out << **Runner << "\t";
    404       TotalDegree += (*Runner)->BondDegree;
    405     }
    406     *out << " -- TotalDegree: " << TotalDegree << endl;
    407     return true;
    408   } else
    409     return false;
    410 };
    411 
    412 /** Output of atom::nr along with all bond partners.
    413  * \param *AdjacencyFile output stream
    414  */
    415 void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
    416 {
    417   *AdjacencyFile << nr << "\t";
    418   for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
    419     *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
    420   *AdjacencyFile << endl;
    421 };
    422 
    423 /** Puts a given bond into atom::ListOfBonds.
    424  * \param *Binder bond to insert
    425  */
    426 bool BondedParticle::RegisterBond(bond *Binder)
    427 {
    428   bool status = false;
    429   if (Binder != NULL) {
    430     if (Binder->Contains(this)) {
    431       ListOfBonds.push_back(Binder);
    432       status = true;
    433     } else {
    434       cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
    435     }
    436   } else {
    437     cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
    438   }
    439   return status;
    440 };
    441 
    442 /** Removes a given bond from atom::ListOfBonds.
    443  * \param *Binder bond to remove
    444  */
    445 bool BondedParticle::UnregisterBond(bond *Binder)
    446 {
    447   bool status = false;
    448   if (Binder != NULL) {
    449     if (Binder->Contains(this)) {
    450       ListOfBonds.remove(Binder);
    451       status = true;
    452     } else {
    453       cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
    454     }
    455   } else {
    456     cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
    457   }
    458   return status;
    459 };
    460 
    461 /** Removes all bonds from atom::ListOfBonds.
    462  * \note Does not do any memory de-allocation.
    463  */
    464 void BondedParticle::UnregisterAllBond()
    465 {
    466   ListOfBonds.clear();
    467 };
    468 
    469 /** Corrects the bond degree by one at most if necessary.
    470  * \param *out output stream for debugging
    471  */
    472 int BondedParticle::CorrectBondDegree(ofstream *out)
    473 {
    474   int NoBonds = 0;
    475   int OtherNoBonds = 0;
    476   int FalseBondDegree = 0;
    477   atom *OtherWalker = NULL;
    478   bond *CandidateBond = NULL;
    479 
    480   *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
    481   NoBonds = CountBonds();
    482   if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
    483     for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
    484       OtherWalker = (*Runner)->GetOtherAtom(this);
    485       OtherNoBonds = OtherWalker->CountBonds();
    486       *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
    487       if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
    488         if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
    489           CandidateBond = (*Runner);
    490           *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
    491         }
    492       }
    493     }
    494     if ((CandidateBond != NULL)) {
    495       CandidateBond->BondDegree++;
    496       *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
    497     } else {
    498       *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
    499       FalseBondDegree++;
    500     }
    501   }
    502   return FalseBondDegree;
    503 };
    504 
    505 /** Adds kinetic energy of this atom to given temperature value.
    506  * \param *temperature add on this value
    507  * \param step given step of trajectory to add
    508  */
    509 void TrajectoryParticle::AddKineticToTemperature(double *temperature, int step) const
    510 {
    511   for (int i=NDIM;i--;)
    512     *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
    513 };
    514 
    515 /** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
    516  * \param startstep trajectory begins at
    517  * \param endstep trajectory ends at
    518  * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
    519  * \param *Force Force matrix to store result in
    520  */
    521 void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
    522 {
    523   double constant = 10.;
    524   TrajectoryParticle *Sprinter = PermutationMap[nr];
    525   // set forces
    526   for (int i=NDIM;i++;)
    527     Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
    528 };
    529 
    530 /** Correct velocity against the summed \a CoGVelocity for \a step.
    531  * \param *ActualTemp sum up actual temperature meanwhile
    532  * \param Step MD step in atom::Tracjetory
    533  * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
    534  */
    535 void TrajectoryParticle::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
    536 {
    537   for(int d=0;d<NDIM;d++) {
    538     Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
    539     *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
    540   }
    541 };
    542 
    543 /** Extends the trajectory STL vector to the new size.
    544  * Does nothing if \a MaxSteps is smaller than current size.
    545  * \param MaxSteps
    546  */
    547 void TrajectoryParticle::ResizeTrajectory(int MaxSteps)
    548 {
    549   if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
    550     //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
    551     Trajectory.R.resize(MaxSteps+1);
    552     Trajectory.U.resize(MaxSteps+1);
    553     Trajectory.F.resize(MaxSteps+1);
    554   }
    555 };
    556 
    557 /** Copies a given trajectory step \a src onto another \a dest
    558  * \param dest index of destination step
    559  * \param src index of source step
    560  */
    561 void TrajectoryParticle::CopyStepOnStep(int dest, int src)
    562 {
    563   if (dest == src)  // self assignment check
    564     return;
    565 
    566   for (int n=NDIM;n--;) {
    567     Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
    568     Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
    569     Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[n];
    570   }
    571 };
    572 
    573 /** Performs a velocity verlet update of the trajectory.
    574  * Parameters are according to those in configuration class.
    575  * \param NextStep index of sequential step to set
    576  * \param *configuration pointer to configuration with parameters
    577  * \param *Force matrix with forces
    578  */
    579 void TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
    580 {
    581   //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
    582   for (int d=0; d<NDIM; d++) {
    583     Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
    584     Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
    585     Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
    586     Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
    587   }
    588   // Update U
    589   for (int d=0; d<NDIM; d++) {
    590     Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
    591     Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
    592   }
    593   // Update R (and F)
    594 //      out << "Integrated position&velocity of step " << (NextStep) << ": (";
    595 //      for (int d=0;d<NDIM;d++)
    596 //        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
    597 //      out << ")\t(";
    598 //      for (int d=0;d<NDIM;d++)
    599 //        cout << Trajectory.U.at(NextStep).x[d] << " ";          // next step
    600 //      out << ")" << endl;
    601 };
    602 
    603 /** Sums up mass and kinetics.
    604  * \param Step step to sum for
    605  * \param *TotalMass pointer to total mass sum
    606  * \param *TotalVelocity pointer to tota velocity sum
    607  */
    608 void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
    609 {
    610   *TotalMass += type->mass;  // sum up total mass
    611   for(int d=0;d<NDIM;d++) {
    612     TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
    613   }
    614 };
    615 
    616 /** Scales velocity of atom according to Woodcock thermostat.
    617  * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
    618  * \param Step MD step to scale
    619  * \param *ekin sum of kinetic energy
    620  */
    621 void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
    622 {
    623   double *U = Trajectory.U.at(Step).x;
    624   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    625     for (int d=0; d<NDIM; d++) {
    626       U[d] *= ScaleTempFactor;
    627       *ekin += 0.5*type->mass * U[d]*U[d];
    628     }
    629 };
    630 
    631 /** Scales velocity of atom according to Gaussian thermostat.
    632  * \param Step MD step to scale
    633  * \param *G
    634  * \param *E
    635  */
    636 void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
    637 {
    638   double *U = Trajectory.U.at(Step).x;
    639   double *F = Trajectory.F.at(Step).x;
    640   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    641     for (int d=0; d<NDIM; d++) {
    642       *G += U[d] * F[d];
    643       *E += U[d]*U[d]*type->mass;
    644     }
    645 };
    646 
    647 /** Determines scale factors according to Gaussian thermostat.
    648  * \param Step MD step to scale
    649  * \param GE G over E ratio
    650  * \param *ekin sum of kinetic energy
    651  * \param *configuration configuration class with TempFrequency and TargetTemp
    652  */
    653 void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
    654 {
    655   double *U = Trajectory.U.at(Step).x;
    656   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    657     for (int d=0; d<NDIM; d++) {
    658       U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
    659       *ekin += type->mass * U[d]*U[d];
    660     }
    661 };
    662 
    663 /** Scales velocity of atom according to Langevin thermostat.
    664  * \param Step MD step to scale
    665  * \param *r random number generator
    666  * \param *ekin sum of kinetic energy
    667  * \param *configuration configuration class with TempFrequency and TargetTemp
    668  */
    669 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
    670 {
    671   double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    672   double *U = Trajectory.U.at(Step).x;
    673   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    674     // throw a dice to determine whether it gets hit by a heat bath particle
    675     if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
    676       cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
    677       // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
    678       for (int d=0; d<NDIM; d++) {
    679         U[d] = gsl_ran_gaussian (r, sigma);
    680       }
    681       cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
    682     }
    683     for (int d=0; d<NDIM; d++)
    684       *ekin += 0.5*type->mass * U[d]*U[d];
    685   }
    686 };
    687 
    688 /** Scales velocity of atom according to Berendsen thermostat.
    689  * \param Step MD step to scale
    690  * \param ScaleTempFactor factor to scale energy (not velocity!) with
    691  * \param *ekin sum of kinetic energy
    692  * \param *configuration configuration class with TempFrequency and Deltat
    693  */
    694 void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
    695 {
    696   double *U = Trajectory.U.at(Step).x;
    697   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    698     for (int d=0; d<NDIM; d++) {
    699       U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
    700       *ekin += 0.5*type->mass * U[d]*U[d];
    701     }
    702   }
    703 };
    704 
    705 /** Initializes current run of NoseHoover thermostat.
    706  * \param Step MD step to scale
    707  * \param *delta_alpha additional sum of kinetic energy on return
    708  */
    709 void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
    710 {
    711   double *U = Trajectory.U.at(Step).x;
    712   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    713     for (int d=0; d<NDIM; d++) {
    714       *delta_alpha += U[d]*U[d]*type->mass;
    715     }
    716   }
    717 };
    718 
    719 /** Initializes current run of NoseHoover thermostat.
    720  * \param Step MD step to scale
    721  * \param *ekin sum of kinetic energy
    722  * \param *configuration configuration class with TempFrequency and Deltat
    723  */
    724 void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
    725 {
    726   double *U = Trajectory.U.at(Step).x;
    727   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    728     for (int d=0; d<NDIM; d++) {
    729         U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
    730         *ekin += (0.5*type->mass) * U[d]*U[d];
    731       }
    732   }
    733 };

  • src/atom.hpp

    rc75de1 r6b919f8  
    2222#include <vector>
    2323
    24 #include <gsl/gsl_randist.h>
    25 
     24#include "atom_atominfo.hpp"
     25#include "atom_bondedparticle.hpp"
     26#include "atom_graphnode.hpp"
     27#include "atom_particleinfo.hpp"
     28#include "atom_trajectoryparticle.hpp"
    2629#include "tesselation.hpp"
    2730
    2831/****************************************** forward declarations *****************************/
    2932
    30 class atom;
    31 class bond;
    32 class config;
    33 class element;
    34 class ForceMatrix;
    3533class Vector;
    3634
    37 #define BondList list<bond *>
    38 
    3935/********************************************** declarations *******************************/
    40 
    41 class AtomInfo {
    42 public:
    43   Vector x;       //!< coordinate vector of atom, giving last position within cell
    44   Vector v;       //!< velocity vector of atom, giving last velocity within cell
    45   Vector F;       //!< Force vector of atom, giving last force within cell
    46   element *type;  //!< pointing to element
    47 
    48   AtomInfo();
    49   ~AtomInfo();
    50 
    51 private:
    52 };
    53 
    54 class TrajectoryParticleInfo {
    55 public:
    56   struct
    57   {
    58     vector<Vector> R;  //!< position vector
    59     vector<Vector> U;  //!< velocity vector
    60     vector<Vector> F;  //!< last force vector
    61   } Trajectory;
    62   int FixedIon;   //!< config variable that states whether forces act on the ion or not
    63 
    64   TrajectoryParticleInfo();
    65   ~TrajectoryParticleInfo();
    66 
    67 private:
    68 };
    69 
    70 class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
    71 public:
    72 
    73   TrajectoryParticle();
    74   virtual ~TrajectoryParticle();
    75 
    76   // constraint potential and dynamics stuff
    77   void AddKineticToTemperature(double *temperature, int step) const;
    78   void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    79   void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
    80 
    81   // trajectory stuff
    82   void ResizeTrajectory(int MaxSteps);
    83   void CopyStepOnStep(int dest, int src);
    84   void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
    85   void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
    86 
    87   // thermostats
    88   void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
    89   void Thermostat_Gaussian_init(int Step, double *G, double *E);
    90   void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
    91   void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
    92   void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
    93   void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
    94   void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
    95 
    96 private:
    97 
    98 };
    99 
    100 class GraphNode;
    101 
    102 class GraphNodeInfo {
    103 public:
    104 
    105   int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
    106   int *ComponentNr;//!< belongs to this non-separable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
    107   int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect non-separable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
    108   bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
    109   bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
    110   atom *Ancestor; //!< "Father" in Depth-First-Search
    111 
    112   GraphNodeInfo();
    113   ~GraphNodeInfo();
    114 private:
    115 
    116 };
    117 
    118 
    119 class GraphNode : public GraphNodeInfo, public virtual ParticleInfo {
    120 public:
    121 
    122   GraphNode();
    123   virtual ~GraphNode();
    124 
    125   void OutputGraphInfo(ofstream *out) const;
    126   void OutputComponentNumber(ofstream *out) const;
    127 
    128 private:
    129 
    130 };
    131 
    132 class BondedParticleInfo {
    133 public:
    134   unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
    135   bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
    136   BondList ListOfBonds; //!< list of all bonds
    137 
    138   BondedParticleInfo();
    139   ~BondedParticleInfo();
    140 
    141 private:
    142 
    143 };
    144 
    145 class BondedParticle : public BondedParticleInfo, public virtual ParticleInfo, public virtual AtomInfo {
    146 public:
    147   BondedParticle();
    148   virtual ~BondedParticle();
    149 
    150   bool RegisterBond(bond *Binder);
    151   bool UnregisterBond(bond *Binder);
    152   void UnregisterAllBond();
    153   int CountBonds() const;
    154   int CorrectBondDegree(ofstream *out);
    155   bool OutputBondOfAtom(ofstream *out) const;
    156   void OutputAdjacency(ofstream *AdjacencyFile) const;
    157   void OutputOrder(ofstream *file);
    158 
    159 private:
    160 
    161 };
    16236
    16337/** Single atom.

  • src/tesselation.cpp

    rc75de1 r6b919f8  
    1616
    1717class molecule;
    18 
    19 // =========================================== ParticleInfo ====================================
    20 
    21 /** Constructor of ParticleInfo.
    22  */
    23 ParticleInfo::ParticleInfo() : nr(-1), Name(NULL) {};
    24 
    25 /** Destructor of ParticleInfo.
    26  */
    27 ParticleInfo::~ParticleInfo()
    28 {
    29   Free(&Name);
    30 };
    3118
    3219// ======================================== Points on Boundary =================================

  • src/tesselation.hpp

    rc75de1 r6b919f8  
    2424#include <set>
    2525
     26#include "atom_particleinfo.hpp"
     27#include "helpers.hpp"
    2628#include "vector.hpp"
    27 #include "helpers.hpp"
    2829
    2930/****************************************** forward declarations *****************************/
    … …  
    144145
    145146// =========================================================== class TESSELPOINT ===========================================
    146 
    147 class ParticleInfo {
    148 public:
    149   int nr;       // index to easierly identify
    150   char *Name;   // some name to reference to on output
    151 
    152   ParticleInfo();
    153   ~ParticleInfo();
    154 
    155   ostream & operator << (ostream &ost);
    156 
    157 private:
    158 };
    159 
    160 ostream & operator << (ostream &ost, const ParticleInfo &a);
    161147
    162148/** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
Note: See TracChangeset for help on using the changeset viewer.