Changeset 14c57a for src/molecule_dynamics.cpp

Timestamp:

Aug 25, 2010, 12:04:11 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

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:

51c3e4

Parents:

8d1dd4

Message:

Made all parts of the programm use the new thermostat structure

File:

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

src/molecule_dynamics.cpp

@@ -24 +24 @@
 #include "LinearAlgebra/Plane.hpp"
 #include "ThermoStatContainer.hpp"
+#include "Thermostats/Berendsen.hpp"
 
 #include <gsl/gsl_matrix.h>
@@ -615 +616 @@
   ActualTemp = 0.;
   ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
-  Thermostats(configuration, ActualTemp, Berendsen);
+  Berendsen berendsen = Berendsen();
+  berendsen.addToContainer(configuration.Thermostats);
+  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
+  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
   MDSteps++;
 
@@ -621 +625 @@
   return true;
 };
-
-/** Implementation of various thermostats.
- * All these thermostats apply an additional force which has the following forms:
- * -# Woodcock
- *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
- * -# Gaussian
- *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
- * -# Langevin
- *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$
- * -# Berendsen
- *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
- * -# Nose-Hoover
- *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
- * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
- * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
- * belatedly and the constraint force set.
- * \param *P Problem at hand
- * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
- * \sa InitThermostat()
- */
-void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
-{
-  double ekin = 0.;
-  double E = 0., G = 0.;
-  double delta_alpha = 0.;
-  double ScaleTempFactor;
-  gsl_rng * r;
-  const gsl_rng_type * T;
-
-  // calculate scale configuration
-  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
-
-  // differentating between the various thermostats
-  switch(Thermostat) {
-     case None:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying no thermostat..." << endl);
-      break;
-     case Woodcock:
-      if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
-        DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
-        ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
-      }
-      break;
-     case Gaussian:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Gaussian thermostat..." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
-
-      DoLog(1) && (Log() << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
-
-      break;
-     case Langevin:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
-      // init random number generator
-      gsl_rng_env_setup();
-      T = gsl_rng_default;
-      r = gsl_rng_alloc (T);
-      // Go through each ion
-      ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
-      break;
-
-     case Berendsen:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
-      break;
-
-     case NoseHoover:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl);
-      // dynamically evolve alpha (the additional degree of freedom)
-      delta_alpha = 0.;
-      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
-      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
-      configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
-      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
-      // apply updated alpha as additional force
-      ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
-      break;
-  }
-  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
-};