Changeset 14c57a


Ignore:
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

Location:
src
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • src/Parser/PcpParser.cpp

    r8d1dd4 r14c57a  
    402402    *file << "DoConstrainedMD\t" << Switches.DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    403403    ASSERT(Thermostats != NULL, "PcpParser::save() - Thermostats not initialized!");
    404     ASSERT(Thermostats->ThermostatNames != NULL, "PcpParser::save() - Thermostats not initialized!");
    405     *file << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
    406     switch(Thermostats->Thermostat) {
    407       default:
    408       case None:
    409         break;
    410       case Woodcock:
    411         *file << Thermostats->ScaleTempStep;
    412         break;
    413       case Gaussian:
    414         *file << Thermostats->ScaleTempStep;
    415         break;
    416       case Langevin:
    417         *file << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
    418         break;
    419       case Berendsen:
    420         *file << Thermostats->TempFrequency;
    421         break;
    422       case NoseHoover:
    423         *file << Thermostats->HooverMass;
    424         break;
    425     };
     404    *file << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
     405    *file << Thermostats->activeThermostat->writeParams();
    426406    *file << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    427407    *file << "CommonWannier\t" << LocalizedOrbitals.CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    … …  
    588568  // read desired Thermostat from file along with needed additional parameters
    589569  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    590     if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
    591       if (Thermostats->ThermostatImplemented[0] == 1) {
    592         Thermostats->Thermostat = None;
    593       } else {
    594         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    595         Thermostats->Thermostat = None;
    596       }
    597     } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
    598       if (Thermostats->ThermostatImplemented[1] == 1) {
    599         Thermostats->Thermostat = Woodcock;
    600         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
    601       } else {
    602         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    603         Thermostats->Thermostat = None;
    604       }
    605     } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
    606       if (Thermostats->ThermostatImplemented[2] == 1) {
    607         Thermostats->Thermostat = Gaussian;
    608         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
    609       } else {
    610         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    611         Thermostats->Thermostat = None;
    612       }
    613     } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
    614       if (Thermostats->ThermostatImplemented[3] == 1) {
    615         Thermostats->Thermostat = Langevin;
    616         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
    617         if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
    618           DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
    619         } else {
    620           Thermostats->alpha = 1.;
    621         }
    622       } else {
    623         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    624         Thermostats->Thermostat = None;
    625       }
    626     } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
    627       if (Thermostats->ThermostatImplemented[4] == 1) {
    628         Thermostats->Thermostat = Berendsen;
    629         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
    630       } else {
    631         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    632         Thermostats->Thermostat = None;
    633       }
    634     } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
    635       if (Thermostats->ThermostatImplemented[5] == 1) {
    636         Thermostats->Thermostat = NoseHoover;
    637         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
    638         Thermostats->alpha = 0.;
    639       } else {
    640         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    641         Thermostats->Thermostat = None;
    642       }
    643     } else {
    644       DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    645       Thermostats->Thermostat = None;
    646     }
     570    Thermostats->makeActive(thermo,fb);
    647571  } else {
    648572    if ((Thermostats->TargetTemp != 0))
    649573      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    650     Thermostats->Thermostat = None;
     574    Thermostats->chooseNone();
    651575  }
    652576  delete[](thermo);

  • src/ThermoStatContainer.cpp

    r8d1dd4 r14c57a  
    2121#include "Helpers/Verbose.hpp"
    2222
     23#include <Thermostats/Berendsen.hpp>
     24#include <Thermostats/GaussianThermostat.hpp>
     25#include <Thermostats/Langevin.hpp>
     26#include <Thermostats/NoseHoover.hpp>
     27#include <Thermostats/NoThermostat.hpp>
     28#include <Thermostats/Woodcock.hpp>
     29
    2330/** Constructor for class ThermoStatContainer.
    2431 *
    2532 */
    2633ThermoStatContainer::ThermoStatContainer() :
    27   Thermostat(4),
    28   ThermostatImplemented(NULL),
    29   ThermostatNames(NULL),
    30   TempFrequency(2.5),
    31   alpha(0.),
    32   HooverMass(0.),
    33   TargetTemp(0.00095004455),
    34   ScaleTempStep(25)
     34  activeThermostat(0),
     35  TargetTemp(0.00095004455)
    3536{
    36   ThermostatImplemented = new int[MaxThermostats];
    37   ThermostatNames = new char *[MaxThermostats];
    38   for (int j=0;j<MaxThermostats;j++)
    39     ThermostatNames[j] = new char[12];
     37  ThermostatTraits<Thermostat> *BerendsenTrait = new ThermostatTraits<Berendsen>();
     38  availThermostats[BerendsenTrait->getName()] = BerendsenTrait;
     39  ThermostatTraits<Thermostat> *GaussianTrait = new ThermostatTraits<GaussianThermostat>();
     40  availThermostats[GaussianTrait->getName()] = GaussianTrait;
     41  ThermostatTraits<Thermostat> *LangevinTrait = new ThermostatTraits<Langevin>();
     42  availThermostats[LangevinTrait->getName()] = LangevinTrait;
     43  ThermostatTraits<Thermostat> *NoseHooverTrait = new ThermostatTraits<NoseHoover>();
     44  availThermostats[NoseHooverTrait->getName()] = NoseHooverTrait;
     45  ThermostatTraits<Thermostat> *NoThermostatTrait = new ThermostatTraits<NoThermostat>();
     46  availThermostats[NoThermostatTrait->getName()] = NoThermostatTrait;
     47  ThermostatTraits<Thermostat> *WoodcockTrait = new ThermostatTraits<Woodcock>();
     48  availThermostats[WoodcockTrait->getName()] = WoodcockTrait;
     49  for(traitsMap::iterator iter = availThermostats.begin();iter!=availThermostats.end();++iter){
     50    cout << iter->first << endl;
     51  }
     52  ASSERT(availThermostats.size()==6,"Not all implemented thermostats referenced!\nDid you check the names in the traits?");
     53  activeThermostat=new Berendsen();
     54}
    4055
    41   strcpy(ThermostatNames[0],"None");
    42   ThermostatImplemented[0] = 1;
    43   strcpy(ThermostatNames[1],"Woodcock");
    44   ThermostatImplemented[1] = 1;
    45   strcpy(ThermostatNames[2],"Gaussian");
    46   ThermostatImplemented[2] = 1;
    47   strcpy(ThermostatNames[3],"Langevin");
    48   ThermostatImplemented[3] = 1;
    49   strcpy(ThermostatNames[4],"Berendsen");
    50   ThermostatImplemented[4] = 1;
    51   strcpy(ThermostatNames[5],"NoseHoover");
    52   ThermostatImplemented[5] = 1;
     56ThermostatTraits<Thermostat> *ThermoStatContainer::getTraitByName(const std::string name){
     57  if(!availThermostats.count(name))
     58    return 0;
     59  return availThermostats[name];
     60}
     61
     62Thermostat *ThermoStatContainer::makeByName(const std::string name,class ConfigFileBuffer * const fb){
     63  ThermostatTraits<Thermostat>* trait = getTraitByName(name);
     64  if(trait){
     65    Thermostat *res = trait->make(fb);
     66    res->addToContainer(this);
     67    return res;
     68  }
     69  else{
     70    return 0;
     71  }
     72}
     73
     74void ThermoStatContainer::makeActive(const std::string name,class ConfigFileBuffer * const fb){
     75  Thermostat* newThermostat = makeByName(name,fb);
     76  if(newThermostat){
     77    if(activeThermostat){
     78      delete activeThermostat;
     79    }
     80    activeThermostat = newThermostat;
     81  }
     82}
     83
     84void ThermoStatContainer::chooseNone(){
     85  if(activeThermostat){
     86    delete activeThermostat;
     87  }
     88  activeThermostat = new NoThermostat();
    5389}
    5490
    … …  
    5894ThermoStatContainer::~ThermoStatContainer()
    5995{
    60   delete[](ThermostatImplemented);
    61   for (int j=0;j<MaxThermostats;j++)
    62     delete[](ThermostatNames[j]);
    63   delete[](ThermostatNames);
     96  if(activeThermostat){
     97    delete activeThermostat;
     98  }
     99  for(traitsMap::iterator iter= availThermostats.begin();iter!=availThermostats.end();++iter){
     100    delete (iter->second);
     101  }
    64102}
    65103

  • src/ThermoStatContainer.hpp

    r8d1dd4 r14c57a  
    1414#endif
    1515
     16#include "Thermostats/Thermostat.hpp"
     17
    1618class ConfigFileBuffer;
    17 
    18 enum { MaxThermostats=6 };      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
    19 enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
    2019
    2120class ThermoStatContainer {
    … …  
    2322  ThermoStatContainer();
    2423  ~ThermoStatContainer();
     24  ThermostatTraits<Thermostat> *getTraitByName(const std::string name);
     25  Thermostat *makeByName(const std::string name,class ConfigFileBuffer * const fb);
     26  void makeActive(const std::string name,class ConfigFileBuffer * const fb);
     27  void chooseNone();
    2528
    26   int Thermostat;
    27   int *ThermostatImplemented;
    28   char **ThermostatNames;
    29   double TempFrequency;
    30   double alpha;
    31   double HooverMass;
     29  Thermostat* activeThermostat;
     30  typedef std::map<std::string,ThermostatTraits<Thermostat>* > traitsMap;
     31  traitsMap availThermostats;
    3232  double TargetTemp;
    33   int ScaleTempStep;
    3433};
    3534

  • src/Thermostats/Berendsen.cpp

    r8d1dd4 r14c57a  
    2525{}
    2626
    27 ThermostatTraits<class Berendsen>::ThermostatTraits() :
    28   name("Berendsen")
    29 {}
     27const char *ThermostatTraits<Berendsen>::name = "Berendsen";
    3028
    31 Thermostat *ThermostatTraits<class Berendsen>::make(class ConfigFileBuffer * const fb){
     29std::string ThermostatTraits<Berendsen>::getName(){
     30  return ThermostatTraits<Berendsen>::name;
     31}
     32
     33Thermostat *ThermostatTraits<Berendsen>::make(class ConfigFileBuffer * const fb){
    3234  double TempFrequency;
    3335  const int verbose = 0;
    3436  ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
    35   return new class Berendsen(TempFrequency);
     37  return new Berendsen(TempFrequency);
    3638}
    3739
    … …  
    5153double Berendsen::doScaleAtoms(unsigned int step,double ActualTemp,ForwardIterator begin, ForwardIterator end){
    5254  DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
    53   double ekin;
     55  double ekin=0;
    5456  double ScaleTempFactor = getContainer().TargetTemp/ActualTemp;
    5557  for(ForwardIterator iter=begin;iter!=end;++iter){
    … …  
    6466
    6567std::string Berendsen::name(){
    66   return ThermostatTraits<Berendsen>().name;
     68  return ThermostatTraits<Berendsen>::name;
    6769}
    6870

  • src/Thermostats/Berendsen.hpp

    r8d1dd4 r14c57a  
    3636struct ThermostatTraits<class Berendsen> : public ThermostatTraits<Thermostat>
    3737{
    38   ThermostatTraits();
    3938  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    40   const char* name;
     39  virtual std::string getName();
     40  static const char* name;
    4141};
    4242

  • src/Thermostats/GaussianThermostat.cpp

    r8d1dd4 r14c57a  
    3030{}
    3131
    32 ThermostatTraits<GaussianThermostat>::ThermostatTraits() :
    33   name("Gaussian")
    34 {}
     32const char *ThermostatTraits<GaussianThermostat>::name = "Gaussian";
     33
     34std::string ThermostatTraits<GaussianThermostat>::getName(){
     35  return ThermostatTraits<GaussianThermostat>::name;
     36}
    3537
    3638Thermostat *ThermostatTraits<GaussianThermostat>::make(class ConfigFileBuffer * const fb){
    … …  
    9395
    9496std::string GaussianThermostat::name(){
    95   return ThermostatTraits<GaussianThermostat>().name;
     97  return ThermostatTraits<GaussianThermostat>::name;
    9698}
    9799

  • src/Thermostats/GaussianThermostat.hpp

    r8d1dd4 r14c57a  
    4141struct ThermostatTraits<GaussianThermostat> : public ThermostatTraits<Thermostat>
    4242{
    43   ThermostatTraits();
    4443  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    45   const char* name;
     44  virtual std::string getName();
     45  static const char* name;
    4646};
    4747

  • src/Thermostats/Langevin.cpp

    r8d1dd4 r14c57a  
    3232}
    3333
    34 ThermostatTraits<class Langevin>::ThermostatTraits() :
    35   name("Langevin")
    36 {}
     34const char *ThermostatTraits<Langevin>::name = "Langevin";
    3735
    38 Thermostat *ThermostatTraits<class Langevin>::make(class ConfigFileBuffer * const fb){
     36std::string ThermostatTraits<Langevin>::getName(){
     37  return ThermostatTraits<Langevin>::name;
     38}
     39
     40Thermostat *ThermostatTraits<Langevin>::make(class ConfigFileBuffer * const fb){
    3941  double TempFrequency;
    4042  double alpha;
    … …  
    4648    alpha = 1.;
    4749  }
    48   return new class Langevin(TempFrequency,alpha);
     50  return new Langevin(TempFrequency,alpha);
    4951}
    5052
    … …  
    8587
    8688std::string Langevin::name(){
    87   return ThermostatTraits<Langevin>().name;
     89  return ThermostatTraits<Langevin>::name;
    8890}
    8991

  • src/Thermostats/Langevin.hpp

    r8d1dd4 r14c57a  
    3939struct ThermostatTraits<class Langevin> : public ThermostatTraits<Thermostat>
    4040{
    41   ThermostatTraits();
    4241  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    43   const char* name;
     42  virtual std::string getName();
     43  static const char* name;
    4444};
    4545

  • src/Thermostats/NoThermostat.cpp

    r8d1dd4 r14c57a  
    1414{}
    1515
    16 ThermostatTraits<class NoThermostat>::ThermostatTraits() :
    17   name("None")
    18 {}
     16const char *ThermostatTraits<NoThermostat>::name = "None";
    1917
    20 Thermostat *ThermostatTraits<class NoThermostat>::make(class ConfigFileBuffer * const fb){
    21   return new class NoThermostat();
     18std::string ThermostatTraits<NoThermostat>::getName(){
     19  return ThermostatTraits<NoThermostat>::name;
     20}
     21
     22Thermostat *ThermostatTraits<NoThermostat>::make(class ConfigFileBuffer * const fb){
     23  return new NoThermostat();
    2224}
    2325
    … …  
    3537
    3638std::string NoThermostat::name(){
    37   return ThermostatTraits<NoThermostat>().name;
     39  return ThermostatTraits<NoThermostat>::name;
    3840}
    3941

  • src/Thermostats/NoThermostat.hpp

    r8d1dd4 r14c57a  
    2828struct ThermostatTraits<class NoThermostat>: public ThermostatTraits<Thermostat>
    2929{
    30   ThermostatTraits();
    3130  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    32   const char* name;
     31  virtual std::string getName();
     32  static const char* name;
    3333};
    3434

  • src/Thermostats/NoseHoover.cpp

    r8d1dd4 r14c57a  
    2727{}
    2828
    29 ThermostatTraits<class NoseHoover>::ThermostatTraits() :
    30   name("NoseHoover")
    31 {}
     29const char *ThermostatTraits<NoseHoover>::name = "NoseHoover";
    3230
    33 Thermostat *ThermostatTraits<class NoseHoover>::make(class ConfigFileBuffer * const fb){
     31std::string ThermostatTraits<NoseHoover>::getName(){
     32  return ThermostatTraits<NoseHoover>::name;
     33}
     34
     35Thermostat *ThermostatTraits<NoseHoover>::make(class ConfigFileBuffer * const fb){
    3436  double HooverMass;
    3537  const int verbose = 0;
    3638  ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
    37   return new class NoseHoover(HooverMass);
     39  return new NoseHoover(HooverMass);
    3840}
    3941
    … …  
    8284
    8385std::string NoseHoover::name(){
    84   return ThermostatTraits<NoseHoover>().name;
     86  return ThermostatTraits<NoseHoover>::name;
    8587}
    8688

  • src/Thermostats/NoseHoover.hpp

    r8d1dd4 r14c57a  
    4040struct ThermostatTraits<class NoseHoover> : public ThermostatTraits<Thermostat>
    4141{
    42   ThermostatTraits();
    4342  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    44   const char* name;
     43  virtual std::string getName();
     44  static const char* name;
    4545};
    4646

  • src/Thermostats/Thermostat.cpp

    r8d1dd4 r14c57a  
    1717{}
    1818
    19 ThermostatTraits<Thermostat>::ThermostatTraits() :
    20   name("Base")
    21 {}
     19const char *ThermostatTraits<Thermostat>::name = "Base";
    2220
    23 void Thermostat::addToContainer(ThermoStatContainer &_container){
    24   container = &_container;
     21std::string ThermostatTraits<Thermostat>::getName(){
     22  return ThermostatTraits<Thermostat>::name;
     23}
     24
     25void Thermostat::addToContainer(ThermoStatContainer *_container){
     26  container = _container;
    2527}
    2628

  • src/Thermostats/Thermostat.hpp

    r8d1dd4 r14c57a  
    3333
    3434  virtual std::string writeParams()=0;
     35  void addToContainer(ThermoStatContainer *);
    3536protected:
    36   void addToContainer(ThermoStatContainer &);
    3737  ThermoStatContainer &getContainer();
    3838private:
    … …  
    4545template <>
    4646struct ThermostatTraits<Thermostat>{
    47   ThermostatTraits();
    4847  virtual Thermostat *make(class ConfigFileBuffer * const fb)=0;
    49   const char* name;
     48  virtual std::string getName();
     49  static const char* name;
    5050};
    5151

  • src/Thermostats/Woodcock.cpp

    r8d1dd4 r14c57a  
    2727{}
    2828
    29 ThermostatTraits<class Woodcock>::ThermostatTraits() :
    30   name("Woodcock")
    31 {}
     29const char *ThermostatTraits<Woodcock>::name = "Woodcock";
    3230
    33 Thermostat *ThermostatTraits<class Woodcock>::make(class ConfigFileBuffer * const fb){
     31std::string ThermostatTraits<Woodcock>::getName(){
     32  return ThermostatTraits<Woodcock>::name;
     33}
     34
     35Thermostat *ThermostatTraits<Woodcock>::make(class ConfigFileBuffer * const fb){
    3436  int ScaleTempStep;
    3537  const int verbose = 0;
    3638  ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
    37   return new class Woodcock(ScaleTempStep);
     39  return new Woodcock(ScaleTempStep);
    3840}
    3941
    … …  
    6971
    7072std::string Woodcock::name(){
    71   return ThermostatTraits<Woodcock>().name;
     73  return ThermostatTraits<Woodcock>::name;
    7274}
    7375

  • src/Thermostats/Woodcock.hpp

    r8d1dd4 r14c57a  
    3636struct ThermostatTraits<class Woodcock> : public ThermostatTraits<Thermostat>
    3737{
    38   ThermostatTraits();
    3938  virtual Thermostat *make(class ConfigFileBuffer * const fb);
    40   const char* name;
     39  virtual std::string getName();
     40  static const char* name;
    4141};
    4242

  • src/atom_trajectoryparticle.cpp

    r8d1dd4 r14c57a  
    144144};
    145145
    146 /** Scales velocity of atom according to Woodcock thermostat.
    147  * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
    148  * \param Step MD step to scale
    149  * \param *ekin sum of kinetic energy
    150  */
    151 void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
    152 {
    153   Vector &U = Trajectory.U.at(Step);
    154   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    155     for (int d=0; d<NDIM; d++) {
    156       U[d] *= ScaleTempFactor;
    157       *ekin += 0.5*getType()->mass * U[d]*U[d];
    158     }
    159 };
    160 
    161 /** Scales velocity of atom according to Gaussian thermostat.
    162  * \param Step MD step to scale
    163  * \param *G
    164  * \param *E
    165  */
    166 void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
    167 {
    168   Vector &U = Trajectory.U.at(Step);
    169   Vector &F = Trajectory.F.at(Step);
    170   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    171     for (int d=0; d<NDIM; d++) {
    172       *G += U[d] * F[d];
    173       *E += U[d]*U[d]*getType()->mass;
    174     }
    175 };
    176 
    177 /** Determines scale factors according to Gaussian thermostat.
    178  * \param Step MD step to scale
    179  * \param GE G over E ratio
    180  * \param *ekin sum of kinetic energy
    181  * \param *configuration configuration class with TempFrequency and TargetTemp
    182  */
    183 void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
    184 {
    185   Vector &U = Trajectory.U.at(Step);
    186   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    187     for (int d=0; d<NDIM; d++) {
    188       U[d] += configuration->Deltat/getType()->mass * ( (G_over_E) * (U[d]*getType()->mass) );
    189       *ekin += getType()->mass * U[d]*U[d];
    190     }
    191 };
    192 
    193 /** Scales velocity of atom according to Langevin thermostat.
    194  * \param Step MD step to scale
    195  * \param *r random number generator
    196  * \param *ekin sum of kinetic energy
    197  * \param *configuration configuration class with TempFrequency and TargetTemp
    198  */
    199 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
    200 {
    201   double sigma  = sqrt(configuration->Thermostats->TargetTemp/getType()->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    202   Vector &U = Trajectory.U.at(Step);
    203   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    204     // throw a dice to determine whether it gets hit by a heat bath particle
    205     if (((((rand()/(double)RAND_MAX))*configuration->Thermostats->TempFrequency) < 1.)) {
    206       DoLog(3) && (Log() << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ");
    207       // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
    208       for (int d=0; d<NDIM; d++) {
    209         U[d] = gsl_ran_gaussian (r, sigma);
    210       }
    211       DoLog(2) && (Log() << Verbose(2) << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl);
    212     }
    213     for (int d=0; d<NDIM; d++)
    214       *ekin += 0.5*getType()->mass * U[d]*U[d];
    215   }
    216 };
    217 
    218 /** Scales velocity of atom according to Berendsen thermostat.
    219  * \param Step MD step to scale
    220  * \param ScaleTempFactor factor to scale energy (not velocity!) with
    221  * \param *ekin sum of kinetic energy
    222  * \param *configuration configuration class with TempFrequency and Deltat
    223  */
    224 void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
    225 {
    226   Vector &U = Trajectory.U.at(Step);
    227   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    228     for (int d=0; d<NDIM; d++) {
    229       U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1));
    230       *ekin += 0.5*getType()->mass * U[d]*U[d];
    231     }
    232   }
    233 };
    234 
    235 /** Initializes current run of NoseHoover thermostat.
    236  * \param Step MD step to scale
    237  * \param *delta_alpha additional sum of kinetic energy on return
    238  */
    239 void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
    240 {
    241   Vector &U = Trajectory.U.at(Step);
    242   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    243     for (int d=0; d<NDIM; d++) {
    244       *delta_alpha += U[d]*U[d]*getType()->mass;
    245     }
    246   }
    247 };
    248 
    249 /** Initializes current run of NoseHoover thermostat.
    250  * \param Step MD step to scale
    251  * \param *ekin sum of kinetic energy
    252  * \param *configuration configuration class with TempFrequency and Deltat
    253  */
    254 void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
    255 {
    256   Vector &U = Trajectory.U.at(Step);
    257   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    258     for (int d=0; d<NDIM; d++) {
    259         U[d] += configuration->Deltat/getType()->mass * (configuration->Thermostats->alpha * (U[d] * getType()->mass));
    260         *ekin += (0.5*getType()->mass) * U[d]*U[d];
    261       }
    262   }
    263 };
    264 
    265 
    266146std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const
    267147{

  • src/atom_trajectoryparticle.hpp

    r8d1dd4 r14c57a  
    5454  double getKineticEnergy(unsigned int step) const;
    5555
    56   // thermostats
    57   void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
    58   void Thermostat_Gaussian_init(int Step, double *G, double *E);
    59   void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
    60   void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
    61   void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
    62   void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
    63   void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
    64 
    6556  std::ostream & operator << (std::ostream &ost) const;
    6657

  • src/config.cpp

    r8d1dd4 r14c57a  
    622622};
    623623
    624 
    625624/** Stores all elements of config structure from which they can be re-read.
    626625 * \param *filename name of file
    … …  
    650649    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    651650    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    652     *output << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
    653     switch(Thermostats->Thermostat) {
    654       default:
    655       case None:
    656         break;
    657       case Woodcock:
    658         *output << Thermostats->ScaleTempStep;
    659         break;
    660       case Gaussian:
    661         *output << Thermostats->ScaleTempStep;
    662         break;
    663       case Langevin:
    664         *output << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
    665         break;
    666       case Berendsen:
    667         *output << Thermostats->TempFrequency;
    668         break;
    669       case NoseHoover:
    670         *output << Thermostats->HooverMass;
    671         break;
    672     };
     651    *output << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
     652    *output << Thermostats->activeThermostat->writeParams();
    673653    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    674654    *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    … …  
    16551635  // read desired Thermostat from file along with needed additional parameters
    16561636  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    1657     if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
    1658       if (Thermostats->ThermostatImplemented[0] == 1) {
    1659         Thermostats->Thermostat = None;
    1660       } else {
    1661         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1662         Thermostats->Thermostat = None;
    1663       }
    1664     } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
    1665       if (Thermostats->ThermostatImplemented[1] == 1) {
    1666         Thermostats->Thermostat = Woodcock;
    1667         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
    1668       } else {
    1669         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1670         Thermostats->Thermostat = None;
    1671       }
    1672     } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
    1673       if (Thermostats->ThermostatImplemented[2] == 1) {
    1674         Thermostats->Thermostat = Gaussian;
    1675         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
    1676       } else {
    1677         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1678         Thermostats->Thermostat = None;
    1679       }
    1680     } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
    1681       if (Thermostats->ThermostatImplemented[3] == 1) {
    1682         Thermostats->Thermostat = Langevin;
    1683         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
    1684         if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
    1685           DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
    1686         } else {
    1687           Thermostats->alpha = 1.;
    1688         }
    1689       } else {
    1690         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1691         Thermostats->Thermostat = None;
    1692       }
    1693     } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
    1694       if (Thermostats->ThermostatImplemented[4] == 1) {
    1695         Thermostats->Thermostat = Berendsen;
    1696         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
    1697       } else {
    1698         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1699         Thermostats->Thermostat = None;
    1700       }
    1701     } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
    1702       if (Thermostats->ThermostatImplemented[5] == 1) {
    1703         Thermostats->Thermostat = NoseHoover;
    1704         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
    1705         Thermostats->alpha = 0.;
    1706       } else {
    1707         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    1708         Thermostats->Thermostat = None;
    1709       }
    1710     } else {
    1711       DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    1712       Thermostats->Thermostat = None;
    1713     }
     1637    Thermostats->makeActive(thermo,fb);
    17141638  } else {
    17151639    if ((Thermostats->TargetTemp != 0))
    17161640      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    1717     Thermostats->Thermostat = None;
     1641    Thermostats->chooseNone();
    17181642  }
    17191643  delete[](thermo);

  • src/molecule.hpp

    r8d1dd4 r14c57a  
    266266  void ScanForPeriodicCorrection();
    267267  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
    268   void Thermostats(config &configuration, double ActualTemp, int Thermostat);
    269268  double VolumeOfConvexEnvelope(bool IsAngstroem);
    270269

  • src/molecule_dynamics.cpp

    r8d1dd4 r14c57a  
    2424#include "LinearAlgebra/Plane.hpp"
    2525#include "ThermoStatContainer.hpp"
     26#include "Thermostats/Berendsen.hpp"
    2627
    2728#include <gsl/gsl_matrix.h>
    … …  
    615616  ActualTemp = 0.;
    616617  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
    617   Thermostats(configuration, ActualTemp, Berendsen);
     618  Berendsen berendsen = Berendsen();
     619  berendsen.addToContainer(configuration.Thermostats);
     620  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
     621  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
    618622  MDSteps++;
    619623
    … …  
    621625  return true;
    622626};
    623 
    624 /** Implementation of various thermostats.
    625  * All these thermostats apply an additional force which has the following forms:
    626  * -# Woodcock
    627  *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
    628  * -# Gaussian
    629  *  \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$
    630  * -# Langevin
    631  *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$
    632  * -# Berendsen
    633  *  \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$
    634  * -# Nose-Hoover
    635  *  \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$
    636  * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
    637  * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
    638  * belatedly and the constraint force set.
    639  * \param *P Problem at hand
    640  * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
    641  * \sa InitThermostat()
    642  */
    643 void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
    644 {
    645   double ekin = 0.;
    646   double E = 0., G = 0.;
    647   double delta_alpha = 0.;
    648   double ScaleTempFactor;
    649   gsl_rng * r;
    650   const gsl_rng_type * T;
    651 
    652   // calculate scale configuration
    653   ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
    654 
    655   // differentating between the various thermostats
    656   switch(Thermostat) {
    657      case None:
    658       DoLog(2) && (Log() << Verbose(2) <<  "Applying no thermostat..." << endl);
    659       break;
    660      case Woodcock:
    661       if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
    662         DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
    663         ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
    664       }
    665       break;
    666      case Gaussian:
    667       DoLog(2) && (Log() << Verbose(2) <<  "Applying Gaussian thermostat..." << endl);
    668       ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
    669 
    670       DoLog(1) && (Log() << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl);
    671       ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
    672 
    673       break;
    674      case Langevin:
    675       DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
    676       // init random number generator
    677       gsl_rng_env_setup();
    678       T = gsl_rng_default;
    679       r = gsl_rng_alloc (T);
    680       // Go through each ion
    681       ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
    682       break;
    683 
    684      case Berendsen:
    685       DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
    686       ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
    687       break;
    688 
    689      case NoseHoover:
    690       DoLog(2) && (Log() << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl);
    691       // dynamically evolve alpha (the additional degree of freedom)
    692       delta_alpha = 0.;
    693       ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    694       delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
    695       configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
    696       DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
    697       // apply updated alpha as additional force
    698       ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
    699       break;
    700   }
    701   DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
    702 };
Note: See TracChangeset for help on using the changeset viewer.