Changeset 6b919f8 for src/atom.cpp


Ignore:
Timestamp:
Oct 20, 2009, 8:55:17 AM (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:
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@…>

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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 };
Note: See TracChangeset for help on using the changeset viewer.