Changeset 4a7776a for src/atom.cpp

Timestamp:

Oct 12, 2009, 10:30:02 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:

5034e1

Parents:

ccd9f5

Message:

Complete refactoring of molecule_dynamics.cpp

• new small functions: SumDistanceOfTrajectories(), PenalizeEqualTargets(), (), CreateInitialLists(), TryNextNearestNeighbourForInjectivePermutation(), MakeInjectivePermutation()
• Thermostat routines (case blocks) have been shifted to small functions in atom class (and soon to be extracted from there into its own class: ThermostatedParticle ?)
• old functions are a lot smaller now: MinimiseConstrainedPotential(), molecule::EvaluateConstrainedForces(), VerletForceIntegration(), also using ActOnAllAtoms
• atom class has a lot of new small functions:
  • trajectory: ResizeTrajectory(), CopyStepOnStep(), VelocityVerletUpdate(), VelocityVerletUpdate()
  • constraint potential: CorrectVelocity(), Thermostat_...()

File:

• src/atom.cpp (modified) (3 diffs)

src/atom.cpp

@@ -7 +7 @@
 #include "atom.hpp"
 #include "bond.hpp"
+#include "config.hpp"
 #include "element.hpp"
 #include "memoryallocator.hpp"
@@ -265 +266 @@
 };
 
+/** Extends the trajectory STL vector to the new size.
+ * Does nothing if \a MaxSteps is smaller than current size.
+ * \param MaxSteps
+ */
+void atom::ResizeTrajectory(int MaxSteps)
+{
+  if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
+    //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
+    Trajectory.R.resize(MaxSteps+1);
+    Trajectory.U.resize(MaxSteps+1);
+    Trajectory.F.resize(MaxSteps+1);
+  }
+};
+
+/** Copies a given trajectory step \a src onto another \a dest
+ * \param dest index of destination step
+ * \param src index of source step
+ */
+void atom::CopyStepOnStep(int dest, int src)
+{
+  if (dest == src)  // self assignment check
+    return;
+
+  for (int n=NDIM;n--;) {
+    Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
+    Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
+    Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[n];
+  }
+};
+
+/** Performs a velocity verlet update of the trajectory.
+ * Parameters are according to those in configuration class.
+ * \param NextStep index of sequential step to set
+ * \param *configuration pointer to configuration with parameters
+ * \param *Force matrix with forces
+ */
+void atom::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
+{
+  //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
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
+    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
+    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
+  }
+  // Update U
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
+    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
+  }
+  // Update R (and F)
+//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectory.U.at(NextStep).x[d] << " ";          // next step
+//      out << ")" << endl;
+};
+
+/** Sums up mass and kinetics.
+ * \param Step step to sum for
+ * \param *TotalMass pointer to total mass sum
+ * \param *TotalVelocity pointer to tota velocity sum
+ */
+void atom::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
+{
+  *TotalMass += type->mass;  // sum up total mass
+  for(int d=0;d<NDIM;d++) {
+    TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
+  }
+};
+
 /** Returns squared distance to a given vector.
  * \param origin vector to calculate distance to
@@ -312 +386 @@
     Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
 };
+
+/** Correct velocity against the summed \a CoGVelocity for \a step.
+ * \param *ActualTemp sum up actual temperature meanwhile
+ * \param Step MD step in atom::Tracjetory
+ * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
+ */
+void atom::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
+{
+  for(int d=0;d<NDIM;d++) {
+    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
+    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
+  }
+};
+
+/** Scales velocity of atom according to Woodcock thermostat.
+ * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ */
+void atom::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= ScaleTempFactor;
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param *G
+ * \param *E
+ */
+void atom::Thermostat_Gaussian_init(int Step, double *G, double *E)
+{
+  double *U = Trajectory.U.at(Step).x;
+  double *F = Trajectory.F.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *G += U[d] * F[d];
+      *E += U[d]*U[d]*type->mass;
+    }
+};
+
+/** Determines scale factors according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param GE G over E ratio
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void atom::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
+      *ekin += type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Langevin thermostat.
+ * \param Step MD step to scale
+ * \param *r random number generator
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void atom::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
+{
+  double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    // throw a dice to determine whether it gets hit by a heat bath particle
+    if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
+      cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
+      // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
+      for (int d=0; d<NDIM; d++) {
+        U[d] = gsl_ran_gaussian (r, sigma);
+      }
+      cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+    }
+    for (int d=0; d<NDIM; d++)
+      *ekin += 0.5*type->mass * U[d]*U[d];
+  }
+};
+
+/** Scales velocity of atom according to Berendsen thermostat.
+ * \param Step MD step to scale
+ * \param ScaleTempFactor factor to scale energy (not velocity!) with
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void atom::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *delta_alpha additional sum of kinetic energy on return
+ */
+void atom::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *delta_alpha += U[d]*U[d]*type->mass;
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void atom::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+        U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
+        *ekin += (0.5*type->mass) * U[d]*U[d];
+      }
+  }
+};