Changes in / [d17f3d:51c3e4]


Ignore:
Location:
src
Files:
14 added
15 edited

Legend:

Unmodified
Added
Removed

  • src/AtomSet.hpp

    rd17f3d r51c3e4  
    1111
    1212#include <functional>
     13#include <numeric>
    1314#include <algorithm>
    1415#include <boost/foreach.hpp>
    … …  
    2223
    2324#include "atom.hpp"
    24 #include <list>
    2525
    2626// this tests, whether we actually have a Vector
    … …  
    5959  template<class Function>
    6060  void transformNodes(Function f);
     61  double totalTemperatureAtStep(unsigned int i) const;
    6162
    6263private:
    … …  
    6869    }
    6970    Function &f;
     71  };
     72
     73  struct temperatureSum {
     74    temperatureSum(unsigned int _step) : step(_step),value(0) {}
     75    double operator+(atom *atom){
     76      return value + atom->getKineticEnergy(step);
     77    }
     78    double operator=(double _value){
     79      value = _value;
     80      return value;
     81    }
     82    unsigned int step;
     83    double value;
    7084  };
    7185};
    … …  
    86100}
    87101
     102template<class Set>
     103inline double AtomSetMixin<Set>::totalTemperatureAtStep(unsigned int step) const{
     104  return accumulate(this->begin(),this->end(),temperatureSum(step)).value;
     105}
     106
    88107// allows simpler definition of AtomSets
    89108#define ATOMSET(container_type) AtomSetMixin<container_type<atom*> >

  • src/LinearAlgebra/VectorContent.hpp

    rd17f3d r51c3e4  
    1414 * Because of this the pointer to the gsl_vector struct is wrapped inside another
    1515 * (dumb) object that allows for forward definitions.
     16 *
     17 * DO NOT USE OUTSIDE OF VECTOR UNLESS YOU ABSOLUTELY HAVE TO AND KNOW WHAT YOU ARE DOING.
    1618 */
    1719

  • src/Makefile.am

    rd17f3d r51c3e4  
    100100  Descriptors/MoleculePtrDescriptor.hpp \
    101101  Descriptors/MoleculeSelectionDescriptor.cpp
    102 
     102 
     103THERMOSTATSOURCE = \
     104  Thermostats/Berendsen.cpp \
     105  Thermostats/GaussianThermostat.cpp \
     106  Thermostats/Langevin.cpp \
     107  Thermostats/NoseHoover.cpp \
     108  Thermostats/NoThermostat.cpp \
     109  Thermostats/Thermostat.cpp \
     110  Thermostats/Woodcock.cpp
     111
     112THERMOSTATHEADER = \
     113  Thermostats/Berendsen.hpp \
     114  Thermostats/GaussianThermostat.hpp \
     115  Thermostats/Langevin.hpp \
     116  Thermostats/NoseHoover.hpp \
     117  Thermostats/NoThermostat.hpp \
     118  Thermostats/Thermostat.hpp \
     119  Thermostats/Woodcock.hpp
    103120
    104121QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
    … …  
    153170  ${SHAPESOURCE} \
    154171  ${DESCRIPTORSOURCE} \
     172  ${THERMOSTATSOURCE} \
    155173  ${TESSELATIONSOURCE} \
    156174  bond.cpp \
    … …  
    190208  ${SHAPEHEADER} \
    191209  ${DESCRIPTORHEADER} \
     210  ${THERMOSTATHEADER} \
    192211  ${TESSELATIONHEADER} \
    193212  bond.hpp \

  • src/Parser/PcpParser.cpp

    rd17f3d r51c3e4  
    409409    *file << "DoConstrainedMD\t" << Switches.DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    410410    ASSERT(Thermostats != NULL, "PcpParser::save() - Thermostats not initialized!");
    411     ASSERT(Thermostats->ThermostatNames != NULL, "PcpParser::save() - Thermostats not initialized!");
    412     *file << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
    413     switch(Thermostats->Thermostat) {
    414       default:
    415       case None:
    416         break;
    417       case Woodcock:
    418         *file << Thermostats->ScaleTempStep;
    419         break;
    420       case Gaussian:
    421         *file << Thermostats->ScaleTempStep;
    422         break;
    423       case Langevin:
    424         *file << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
    425         break;
    426       case Berendsen:
    427         *file << Thermostats->TempFrequency;
    428         break;
    429       case NoseHoover:
    430         *file << Thermostats->HooverMass;
    431         break;
    432     };
     411    *file << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
     412    *file << Thermostats->activeThermostat->writeParams();
    433413    *file << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    434414    *file << "CommonWannier\t" << LocalizedOrbitals.CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    … …  
    595575  // read desired Thermostat from file along with needed additional parameters
    596576  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    597     if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
    598       if (Thermostats->ThermostatImplemented[0] == 1) {
    599         Thermostats->Thermostat = None;
    600       } else {
    601         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    602         Thermostats->Thermostat = None;
    603       }
    604     } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
    605       if (Thermostats->ThermostatImplemented[1] == 1) {
    606         Thermostats->Thermostat = Woodcock;
    607         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
    608       } else {
    609         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    610         Thermostats->Thermostat = None;
    611       }
    612     } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
    613       if (Thermostats->ThermostatImplemented[2] == 1) {
    614         Thermostats->Thermostat = Gaussian;
    615         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
    616       } else {
    617         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    618         Thermostats->Thermostat = None;
    619       }
    620     } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
    621       if (Thermostats->ThermostatImplemented[3] == 1) {
    622         Thermostats->Thermostat = Langevin;
    623         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
    624         if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
    625           DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
    626         } else {
    627           Thermostats->alpha = 1.;
    628         }
    629       } else {
    630         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    631         Thermostats->Thermostat = None;
    632       }
    633     } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
    634       if (Thermostats->ThermostatImplemented[4] == 1) {
    635         Thermostats->Thermostat = Berendsen;
    636         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
    637       } else {
    638         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    639         Thermostats->Thermostat = None;
    640       }
    641     } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
    642       if (Thermostats->ThermostatImplemented[5] == 1) {
    643         Thermostats->Thermostat = NoseHoover;
    644         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
    645         Thermostats->alpha = 0.;
    646       } else {
    647         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    648         Thermostats->Thermostat = None;
    649       }
    650     } else {
    651       DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    652       Thermostats->Thermostat = None;
    653     }
     577    Thermostats->makeActive(thermo,fb);
    654578  } else {
    655579    if ((Thermostats->TargetTemp != 0))
    656580      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    657     Thermostats->Thermostat = None;
     581    Thermostats->chooseNone();
    658582  }
    659583  delete[](thermo);

  • src/ThermoStatContainer.cpp

    rd17f3d r51c3e4  
    2828#include "Helpers/Verbose.hpp"
    2929
     30#include <Thermostats/Berendsen.hpp>
     31#include <Thermostats/GaussianThermostat.hpp>
     32#include <Thermostats/Langevin.hpp>
     33#include <Thermostats/NoseHoover.hpp>
     34#include <Thermostats/NoThermostat.hpp>
     35#include <Thermostats/Woodcock.hpp>
     36
    3037/** Constructor for class ThermoStatContainer.
    3138 *
    3239 */
    3340ThermoStatContainer::ThermoStatContainer() :
    34   Thermostat(4),
    35   ThermostatImplemented(NULL),
    36   ThermostatNames(NULL),
    37   TempFrequency(2.5),
    38   alpha(0.),
    39   HooverMass(0.),
    40   TargetTemp(0.00095004455),
    41   ScaleTempStep(25)
     41  activeThermostat(0),
     42  TargetTemp(0.00095004455)
    4243{
    43   ThermostatImplemented = new int[MaxThermostats];
    44   ThermostatNames = new char *[MaxThermostats];
    45   for (int j=0;j<MaxThermostats;j++)
    46     ThermostatNames[j] = new char[12];
     44  ThermostatTraits<Thermostat> *BerendsenTrait = new ThermostatTraits<Berendsen>();
     45  availThermostats[BerendsenTrait->getName()] = BerendsenTrait;
     46  ThermostatTraits<Thermostat> *GaussianTrait = new ThermostatTraits<GaussianThermostat>();
     47  availThermostats[GaussianTrait->getName()] = GaussianTrait;
     48  ThermostatTraits<Thermostat> *LangevinTrait = new ThermostatTraits<Langevin>();
     49  availThermostats[LangevinTrait->getName()] = LangevinTrait;
     50  ThermostatTraits<Thermostat> *NoseHooverTrait = new ThermostatTraits<NoseHoover>();
     51  availThermostats[NoseHooverTrait->getName()] = NoseHooverTrait;
     52  ThermostatTraits<Thermostat> *NoThermostatTrait = new ThermostatTraits<NoThermostat>();
     53  availThermostats[NoThermostatTrait->getName()] = NoThermostatTrait;
     54  ThermostatTraits<Thermostat> *WoodcockTrait = new ThermostatTraits<Woodcock>();
     55  availThermostats[WoodcockTrait->getName()] = WoodcockTrait;
     56  for(traitsMap::iterator iter = availThermostats.begin();iter!=availThermostats.end();++iter){
     57    cout << iter->first << endl;
     58  }
     59  ASSERT(availThermostats.size()==6,"Not all implemented thermostats referenced!\nDid you check the names in the traits?");
     60  activeThermostat=new Berendsen();
     61}
    4762
    48   strcpy(ThermostatNames[0],"None");
    49   ThermostatImplemented[0] = 1;
    50   strcpy(ThermostatNames[1],"Woodcock");
    51   ThermostatImplemented[1] = 1;
    52   strcpy(ThermostatNames[2],"Gaussian");
    53   ThermostatImplemented[2] = 1;
    54   strcpy(ThermostatNames[3],"Langevin");
    55   ThermostatImplemented[3] = 1;
    56   strcpy(ThermostatNames[4],"Berendsen");
    57   ThermostatImplemented[4] = 1;
    58   strcpy(ThermostatNames[5],"NoseHoover");
    59   ThermostatImplemented[5] = 1;
     63ThermostatTraits<Thermostat> *ThermoStatContainer::getTraitByName(const std::string name){
     64  if(!availThermostats.count(name))
     65    return 0;
     66  return availThermostats[name];
     67}
     68
     69Thermostat *ThermoStatContainer::makeByName(const std::string name,class ConfigFileBuffer * const fb){
     70  ThermostatTraits<Thermostat>* trait = getTraitByName(name);
     71  if(trait){
     72    Thermostat *res = trait->make(fb);
     73    res->addToContainer(this);
     74    return res;
     75  }
     76  else{
     77    return 0;
     78  }
     79}
     80
     81void ThermoStatContainer::makeActive(const std::string name,class ConfigFileBuffer * const fb){
     82  Thermostat* newThermostat = makeByName(name,fb);
     83  if(newThermostat){
     84    if(activeThermostat){
     85      delete activeThermostat;
     86    }
     87    activeThermostat = newThermostat;
     88  }
     89}
     90
     91void ThermoStatContainer::chooseNone(){
     92  if(activeThermostat){
     93    delete activeThermostat;
     94  }
     95  activeThermostat = new NoThermostat();
    6096}
    6197
    … …  
    65101ThermoStatContainer::~ThermoStatContainer()
    66102{
    67   delete[](ThermostatImplemented);
    68   for (int j=0;j<MaxThermostats;j++)
    69     delete[](ThermostatNames[j]);
    70   delete[](ThermostatNames);
     103  if(activeThermostat){
     104    delete activeThermostat;
     105  }
     106  for(traitsMap::iterator iter= availThermostats.begin();iter!=availThermostats.end();++iter){
     107    delete (iter->second);
     108  }
    71109}
    72110

  • src/ThermoStatContainer.hpp

    rd17f3d r51c3e4  
    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   void ParseThermostats(class ConfigFileBuffer * const fb);
    27 
    28   int Thermostat;
    29   int *ThermostatImplemented;
    30   char **ThermostatNames;
    31   double TempFrequency;
    32   double alpha;
    33   double HooverMass;
     29  Thermostat* activeThermostat;
     30  typedef std::map<std::string,ThermostatTraits<Thermostat>* > traitsMap;
     31  traitsMap availThermostats;
    3432  double TargetTemp;
    35   int ScaleTempStep;
    3633};
    3734

  • src/atom_trajectoryparticle.cpp

    rd17f3d r51c3e4  
    4242};
    4343
     44/**
     45 *  returns the kinetic energy of this atom at a given time step
     46 */
    4447
    45 /** Adds kinetic energy of this atom to given temperature value.
    46  * \param *temperature add on this value
    47  * \param step given step of trajectory to add
    48  */
    49 void TrajectoryParticle::AddKineticToTemperature(double *temperature, int step) const
    50 {
    51   for (int i=NDIM;i--;)
    52     *temperature += getType()->getMass() * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i];
    53 };
     48double TrajectoryParticle::getKineticEnergy(unsigned int step) const{
     49  return getType()->getMass() * Trajectory.U.at(step).NormSquared();
     50}
    5451
    5552/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
    … …  
    155152};
    156153
    157 /** Scales velocity of atom according to Woodcock thermostat.
    158  * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
    159  * \param Step MD step to scale
    160  * \param *ekin sum of kinetic energy
    161  */
    162 void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
    163 {
    164   Vector &U = Trajectory.U.at(Step);
    165   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    166     for (int d=0; d<NDIM; d++) {
    167       U[d] *= ScaleTempFactor;
    168       *ekin += 0.5*getType()->getMass() * U[d]*U[d];
    169     }
    170 };
    171 
    172 /** Scales velocity of atom according to Gaussian thermostat.
    173  * \param Step MD step to scale
    174  * \param *G
    175  * \param *E
    176  */
    177 void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
    178 {
    179   Vector &U = Trajectory.U.at(Step);
    180   Vector &F = Trajectory.F.at(Step);
    181   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    182     for (int d=0; d<NDIM; d++) {
    183       *G += U[d] * F[d];
    184       *E += U[d]*U[d]*getType()->getMass();
    185     }
    186 };
    187 
    188 /** Determines scale factors according to Gaussian thermostat.
    189  * \param Step MD step to scale
    190  * \param GE G over E ratio
    191  * \param *ekin sum of kinetic energy
    192  * \param *configuration configuration class with TempFrequency and TargetTemp
    193  */
    194 void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
    195 {
    196   Vector &U = Trajectory.U.at(Step);
    197   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
    198     for (int d=0; d<NDIM; d++) {
    199       U[d] += configuration->Deltat/getType()->getMass() * ( (G_over_E) * (U[d]*getType()->getMass()) );
    200       *ekin += getType()->getMass() * U[d]*U[d];
    201     }
    202 };
    203 
    204 /** Scales velocity of atom according to Langevin thermostat.
    205  * \param Step MD step to scale
    206  * \param *r random number generator
    207  * \param *ekin sum of kinetic energy
    208  * \param *configuration configuration class with TempFrequency and TargetTemp
    209  */
    210 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
    211 {
    212   double sigma  = sqrt(configuration->Thermostats->TargetTemp/getType()->getMass()); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    213   Vector &U = Trajectory.U.at(Step);
    214   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    215     // throw a dice to determine whether it gets hit by a heat bath particle
    216     if (((((rand()/(double)RAND_MAX))*configuration->Thermostats->TempFrequency) < 1.)) {
    217       DoLog(3) && (Log() << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ");
    218       // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
    219       for (int d=0; d<NDIM; d++) {
    220         U[d] = gsl_ran_gaussian (r, sigma);
    221       }
    222       DoLog(2) && (Log() << Verbose(2) << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl);
    223     }
    224     for (int d=0; d<NDIM; d++)
    225       *ekin += 0.5*getType()->getMass() * U[d]*U[d];
    226   }
    227 };
    228 
    229 /** Scales velocity of atom according to Berendsen thermostat.
    230  * \param Step MD step to scale
    231  * \param ScaleTempFactor factor to scale energy (not velocity!) with
    232  * \param *ekin sum of kinetic energy
    233  * \param *configuration configuration class with TempFrequency and Deltat
    234  */
    235 void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
    236 {
    237   Vector &U = Trajectory.U.at(Step);
    238   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    239     for (int d=0; d<NDIM; d++) {
    240       U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1));
    241       *ekin += 0.5*getType()->getMass() * U[d]*U[d];
    242     }
    243   }
    244 };
    245 
    246 /** Initializes current run of NoseHoover thermostat.
    247  * \param Step MD step to scale
    248  * \param *delta_alpha additional sum of kinetic energy on return
    249  */
    250 void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
    251 {
    252   Vector &U = Trajectory.U.at(Step);
    253   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    254     for (int d=0; d<NDIM; d++) {
    255       *delta_alpha += U[d]*U[d]*getType()->getMass();
    256     }
    257   }
    258 };
    259 
    260 /** Initializes current run of NoseHoover thermostat.
    261  * \param Step MD step to scale
    262  * \param *ekin sum of kinetic energy
    263  * \param *configuration configuration class with TempFrequency and Deltat
    264  */
    265 void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
    266 {
    267   Vector &U = Trajectory.U.at(Step);
    268   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    269     for (int d=0; d<NDIM; d++) {
    270         U[d] += configuration->Deltat/getType()->getMass() * (configuration->Thermostats->alpha * (U[d] * getType()->getMass()));
    271         *ekin += (0.5*getType()->getMass()) * U[d]*U[d];
    272       }
    273   }
    274 };
    275 
    276 
    277154std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const
    278155{

  • src/atom_trajectoryparticle.hpp

    rd17f3d r51c3e4  
    4343
    4444  // constraint potential and dynamics stuff
    45   void AddKineticToTemperature(double *temperature, int step) const;
    4645  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force) const;
    4746  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
    … …  
    5251  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
    5352  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const;
    54 
    55   // thermostats
    56   void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
    57   void Thermostat_Gaussian_init(int Step, double *G, double *E);
    58   void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
    59   void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
    60   void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
    61   void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
    62   void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
     53  double getKineticEnergy(unsigned int step) const;
    6354
    6455  std::ostream & operator << (std::ostream &ost) const;

  • src/config.cpp

    rd17f3d r51c3e4  
    629629};
    630630
    631 
    632 /** Initializes config file structure by loading elements from a give file.
    633  * \param *file input file stream being the opened config file
    634  * \param BondGraphFileName file name of the bond length table file, if string is left blank, no table is parsed.
    635  * \param *periode pointer to a periodentafel class with all elements
    636  * \param *&MolList pointer to MoleculeListClass, on return containing all parsed molecules in system
    637  */
    638 void config::Load(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList)
    639 {
    640   molecule *mol = World::getInstance().createMolecule();
    641   ifstream *file = new ifstream(filename);
    642   if (file == NULL) {
    643     DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    644     return;
    645   }
    646   file->close();
    647   delete(file);
    648 
    649   // ParseParameterFile
    650   class ConfigFileBuffer *FileBuffer = new ConfigFileBuffer(filename);
    651 
    652   /* Oeffne Hauptparameterdatei */
    653   int di = 0;
    654   double BoxLength[9];
    655   string zeile;
    656   string dummy;
    657   int verbose = 0;
    658  
    659   //TODO: This is actually sensible?: if (MaxOuterStep > 0)
    660   ParseThermostats(FileBuffer);
    661  
    662   /* Namen einlesen */
    663 
    664   // 1. parse in options
    665   ParseForParameter(verbose,FileBuffer, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
    666   ParseForParameter(verbose,FileBuffer, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
    667   ParseForParameter(verbose,FileBuffer, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
    668   ParseForParameter(verbose,FileBuffer,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
    669   ParseForParameter(verbose,FileBuffer,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);
    670 
    671   if (!ParseForParameter(verbose,FileBuffer,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))
    672     config::Seed = 1;
    673 
    674   if(!ParseForParameter(verbose,FileBuffer,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {
    675     config::DoOutOrbitals = 0;
    676   } else {
    677     if (config::DoOutOrbitals < 0) config::DoOutOrbitals = 0;
    678     if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1;
    679   }
    680   ParseForParameter(verbose,FileBuffer,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
    681   if (config::DoOutVis < 0) config::DoOutVis = 0;
    682   if (config::DoOutVis > 1) config::DoOutVis = 1;
    683   if (!ParseForParameter(verbose,FileBuffer,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))
    684     config::VectorPlane = -1;
    685   if (!ParseForParameter(verbose,FileBuffer,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))
    686     config::VectorCut = 0.;
    687   ParseForParameter(verbose,FileBuffer,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
    688   if (config::DoOutMes < 0) config::DoOutMes = 0;
    689   if (config::DoOutMes > 1) config::DoOutMes = 1;
    690   if (!ParseForParameter(verbose,FileBuffer,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))
    691     config::DoOutCurrent = 0;
    692   if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
    693   if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
    694   ParseForParameter(verbose,FileBuffer,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
    695   if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
    696   if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
    697   if(!ParseForParameter(verbose,FileBuffer,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {
    698     config::DoWannier = 0;
    699   } else {
    700     if (config::DoWannier < 0) config::DoWannier = 0;
    701     if (config::DoWannier > 1) config::DoWannier = 1;
    702   }
    703   if(!ParseForParameter(verbose,FileBuffer,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {
    704     config::CommonWannier = 0;
    705   } else {
    706     if (config::CommonWannier < 0) config::CommonWannier = 0;
    707     if (config::CommonWannier > 4) config::CommonWannier = 4;
    708   }
    709   if(!ParseForParameter(verbose,FileBuffer,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {
    710     config::SawtoothStart = 0.01;
    711   } else {
    712     if (config::SawtoothStart < 0.) config::SawtoothStart = 0.;
    713     if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
    714   }
    715  
    716   if (ParseForParameter(verbose,FileBuffer,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))
    717     if (config::DoConstrainedMD < 0)
    718       config::DoConstrainedMD = 0;
    719   ParseForParameter(verbose,FileBuffer,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
    720   if (!ParseForParameter(verbose,FileBuffer,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
    721     config::Deltat = 1;
    722   ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    723   ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    724   ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
    725   //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
    726   if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
    727     config::EpsWannier = 1e-8;
    728 
    729   // stop conditions
    730   //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    731   ParseForParameter(verbose,FileBuffer,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    732   if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    733 
    734   ParseForParameter(verbose,FileBuffer,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    735   ParseForParameter(verbose,FileBuffer,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    736   ParseForParameter(verbose,FileBuffer,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
    737   ParseForParameter(verbose,FileBuffer,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);
    738   ParseForParameter(verbose,FileBuffer,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);
    739   if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
    740   if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    741   if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
    742 
    743   ParseForParameter(verbose,FileBuffer,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    744   ParseForParameter(verbose,FileBuffer,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    745   ParseForParameter(verbose,FileBuffer,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
    746   ParseForParameter(verbose,FileBuffer,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
    747   ParseForParameter(verbose,FileBuffer,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);
    748   if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
    749   if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
    750   if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
    751 
    752   // Unit cell and magnetic field
    753   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    754   double * cell_size = new double[6];
    755   cell_size[0] = BoxLength[0];
    756   cell_size[1] = BoxLength[3];
    757   cell_size[2] = BoxLength[4];
    758   cell_size[3] = BoxLength[6];
    759   cell_size[4] = BoxLength[7];
    760   cell_size[5] = BoxLength[8];
    761   World::getInstance().setDomain(cell_size);
    762   delete cell_size;
    763   //if (1) fprintf(stderr,"\n");
    764 
    765   ParseForParameter(verbose,FileBuffer,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);
    766   ParseForParameter(verbose,FileBuffer,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);
    767   if (!ParseForParameter(verbose,FileBuffer,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))
    768     config::DoFullCurrent = 0;
    769   if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
    770   if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
    771   if (config::DoOutNICS < 0) config::DoOutNICS = 0;
    772   if (config::DoOutNICS > 2) config::DoOutNICS = 2;
    773   if (config::DoPerturbation == 0) {
    774     config::DoFullCurrent = 0;
    775     config::DoOutNICS = 0;
    776   }
    777 
    778   ParseForParameter(verbose,FileBuffer,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
    779   ParseForParameter(verbose,FileBuffer,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
    780   ParseForParameter(verbose,FileBuffer,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
    781   if (config::Lev0Factor < 2) {
    782     config::Lev0Factor = 2;
    783   }
    784   ParseForParameter(verbose,FileBuffer,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
    785   if (di >= 0 && di < 2) {
    786     config::RiemannTensor = di;
    787   } else {
    788     fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    789     exit(1);
    790   }
    791   switch (config::RiemannTensor) {
    792     case 0: //UseNoRT
    793       if (config::MaxLevel < 2) {
    794         config::MaxLevel = 2;
    795       }
    796       config::LevRFactor = 2;
    797       config::RTActualUse = 0;
    798       break;
    799     case 1: // UseRT
    800       if (config::MaxLevel < 3) {
    801         config::MaxLevel = 3;
    802       }
    803       ParseForParameter(verbose,FileBuffer,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
    804       if (config::RiemannLevel < 2) {
    805         config::RiemannLevel = 2;
    806       }
    807       if (config::RiemannLevel > config::MaxLevel-1) {
    808         config::RiemannLevel = config::MaxLevel-1;
    809       }
    810       ParseForParameter(verbose,FileBuffer,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
    811       if (config::LevRFactor < 2) {
    812         config::LevRFactor = 2;
    813       }
    814       config::Lev0Factor = 2;
    815       config::RTActualUse = 2;
    816       break;
    817   }
    818   ParseForParameter(verbose,FileBuffer,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
    819   if (di >= 0 && di < 2) {
    820     config::PsiType = di;
    821   } else {
    822     fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    823     exit(1);
    824   }
    825   switch (config::PsiType) {
    826   case 0: // SpinDouble
    827     ParseForParameter(verbose,FileBuffer,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    828     ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    829     break;
    830   case 1: // SpinUpDown
    831     if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    832     ParseForParameter(verbose,FileBuffer,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    833     ParseForParameter(verbose,FileBuffer,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    834     ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    835     break;
    836   }
    837 
    838   // IonsInitRead
    839 
    840   ParseForParameter(verbose,FileBuffer,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    841   ParseForParameter(verbose,FileBuffer,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    842   ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
    843   if (!ParseForParameter(verbose,FileBuffer,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))
    844     config::RelativeCoord = 0;
    845   if (!ParseForParameter(verbose,FileBuffer,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
    846     config::StructOpt = 0;
    847 
    848   // 2. parse the bond graph file if given
    849   if (BG == NULL) {
    850     BG = new BondGraph(IsAngstroem);
    851     if (BG->LoadBondLengthTable(BondGraphFileName)) {
    852       DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
    853     } else {
    854       DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
    855     }
    856   }
    857 
    858   // 3. parse the molecule in
    859   LoadMolecule(mol, FileBuffer, periode, FastParsing);
    860   mol->SetNameFromFilename(filename);
    861   mol->ActiveFlag = true;
    862   MolList->insert(mol);
    863 
    864   // 4. dissect the molecule into connected subgraphs
    865   // don't do this here ...
    866   //FragmentationSubgraphDissection();
    867   //delete(mol);
    868 
    869   delete(FileBuffer);
    870 };
    871 
    872 /** Initializes config file structure by loading elements from a give file with old pcp syntax.
    873  * \param *file input file stream being the opened config file with old pcp syntax
    874  * \param BondGraphFileName file name of the bond length table file, if string is left blank, no table is parsed.
    875  * \param *periode pointer to a periodentafel class with all elements
    876  * \param *&MolList pointer to MoleculeListClass, on return containing all parsed molecules in system
    877  */
    878 void config::LoadOld(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList)
    879 {
    880   molecule *mol = World::getInstance().createMolecule();
    881   ifstream *file = new ifstream(filename);
    882   if (file == NULL) {
    883     DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    884     return;
    885   }
    886   // ParseParameters
    887 
    888   /* Oeffne Hauptparameterdatei */
    889   int l = 0;
    890   int i = 0;
    891   int di = 0;
    892   double a = 0.;
    893   double b = 0.;
    894   double BoxLength[9];
    895   string zeile;
    896   string dummy;
    897   const element *elementhash[128];
    898   int Z = -1;
    899   int No = -1;
    900   int AtomNo = -1;
    901   int found = 0;
    902   int verbose = 0;
    903 
    904   mol->ActiveFlag = true;
    905   MolList->insert(mol);
    906   /* Namen einlesen */
    907 
    908   ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
    909   ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
    910   ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
    911   ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
    912   ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 2, 1, int_type, &(config::ProcPEPsi), 1, critical);
    913   config::Seed = 1;
    914   config::DoOutOrbitals = 0;
    915   ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
    916   if (config::DoOutVis < 0) config::DoOutVis = 0;
    917   if (config::DoOutVis > 1) config::DoOutVis = 1;
    918   config::VectorPlane = -1;
    919   config::VectorCut = 0.;
    920   ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
    921   if (config::DoOutMes < 0) config::DoOutMes = 0;
    922   if (config::DoOutMes > 1) config::DoOutMes = 1;
    923   config::DoOutCurrent = 0;
    924   ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
    925   if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
    926   if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
    927   config::CommonWannier = 0;
    928   config::SawtoothStart = 0.01;
    929 
    930   ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, double_type, &(config::MaxOuterStep), 1, critical);
    931   ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
    932   ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    933   ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    934   ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
    935   ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(Thermostats->ScaleTempStep), 1, optional);
    936   config::EpsWannier = 1e-8;
    937 
    938   // stop conditions
    939   //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    940   ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    941   if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    942 
    943   ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    944   ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    945   ParseForParameter(verbose,file,"MaxMinStep", 0, 3, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
    946   ParseForParameter(verbose,file,"MaxMinStep", 0, 4, 1, int_type, &(config::MaxMinStopStep), 1, critical);
    947   if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
    948   if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    949   config::MaxMinGapStopStep = 1;
    950 
    951   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    952   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    953   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 3, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
    954   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 4, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
    955   if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
    956   if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
    957   config::InitMaxMinGapStopStep = 1;
    958 
    959   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    960   double * cell_size = new double[6];
    961   cell_size[0] = BoxLength[0];
    962   cell_size[1] = BoxLength[3];
    963   cell_size[2] = BoxLength[4];
    964   cell_size[3] = BoxLength[6];
    965   cell_size[4] = BoxLength[7];
    966   cell_size[5] = BoxLength[8];
    967   World::getInstance().setDomain(cell_size);
    968   delete[] cell_size;
    969   if (1) fprintf(stderr,"\n");
    970   config::DoPerturbation = 0;
    971   config::DoFullCurrent = 0;
    972 
    973   ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
    974   ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
    975   ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
    976   if (config::Lev0Factor < 2) {
    977     config::Lev0Factor = 2;
    978   }
    979   ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
    980   if (di >= 0 && di < 2) {
    981     config::RiemannTensor = di;
    982   } else {
    983     fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    984     exit(1);
    985   }
    986   switch (config::RiemannTensor) {
    987     case 0: //UseNoRT
    988       if (config::MaxLevel < 2) {
    989         config::MaxLevel = 2;
    990       }
    991       config::LevRFactor = 2;
    992       config::RTActualUse = 0;
    993       break;
    994     case 1: // UseRT
    995       if (config::MaxLevel < 3) {
    996         config::MaxLevel = 3;
    997       }
    998       ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
    999       if (config::RiemannLevel < 2) {
    1000         config::RiemannLevel = 2;
    1001       }
    1002       if (config::RiemannLevel > config::MaxLevel-1) {
    1003         config::RiemannLevel = config::MaxLevel-1;
    1004       }
    1005       ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
    1006       if (config::LevRFactor < 2) {
    1007         config::LevRFactor = 2;
    1008       }
    1009       config::Lev0Factor = 2;
    1010       config::RTActualUse = 2;
    1011       break;
    1012   }
    1013   ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
    1014   if (di >= 0 && di < 2) {
    1015     config::PsiType = di;
    1016   } else {
    1017     fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    1018     exit(1);
    1019   }
    1020   switch (config::PsiType) {
    1021   case 0: // SpinDouble
    1022     ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    1023     config::AddPsis = 0;
    1024     break;
    1025   case 1: // SpinUpDown
    1026     if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    1027     ParseForParameter(verbose,file,"MaxPsiUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    1028     ParseForParameter(verbose,file,"MaxPsiDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    1029     config::AddPsis = 0;
    1030     break;
    1031   }
    1032 
    1033   // IonsInitRead
    1034 
    1035   ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    1036   ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    1037   config::RelativeCoord = 0;
    1038   config::StructOpt = 0;
    1039 
    1040 
    1041   // 2. parse the bond graph file if given
    1042   BG = new BondGraph(IsAngstroem);
    1043   if (BG->LoadBondLengthTable(BondGraphFileName)) {
    1044     DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
    1045   } else {
    1046     DoLog(0) && (Log() << Verbose(0) << "Bond length table loading failed." << endl);
    1047   }
    1048 
    1049   // Routine from builder.cpp
    1050 
    1051   for (i=MAX_ELEMENTS;i--;)
    1052     elementhash[i] = NULL;
    1053   DoLog(0) && (Log() << Verbose(0) << "Parsing Ions ..." << endl);
    1054   No=0;
    1055   found = 0;
    1056   while (getline(*file,zeile,'\n')) {
    1057     if (zeile.find("Ions_Data") == 0) {
    1058       DoLog(1) && (Log() << Verbose(1) << "found Ions_Data...begin parsing" << endl);
    1059       found ++;
    1060     }
    1061     if (found > 0) {
    1062       if (zeile.find("Ions_Data") == 0)
    1063         getline(*file,zeile,'\n'); // read next line and parse this one
    1064       istringstream input(zeile);
    1065       input >> AtomNo;  // number of atoms
    1066       input >> Z;       // atomic number
    1067       input >> a;
    1068       input >> l;
    1069       input >> l;
    1070       input >> b;     // element mass
    1071       elementhash[No] = periode->FindElement(Z);
    1072       DoLog(1) && (Log() << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"  << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl);
    1073       for(i=0;i<AtomNo;i++) {
    1074         if (!getline(*file,zeile,'\n')) {// parse on and on
    1075           DoLog(2) && (Log() << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl);
    1076           // return 1;
    1077         } else {
    1078           //Log() << Verbose(2) << "Reading line: " << zeile << endl;
    1079         }
    1080         istringstream input2(zeile);
    1081         atom *neues = World::getInstance().createAtom();
    1082         double tmp;
    1083         for (int j=0;j<NDIM;j++) {
    1084           input2 >> tmp;
    1085           neues->set(j,tmp);
    1086         }
    1087         input2 >> l;
    1088         neues->setType(elementhash[No]); // find element type
    1089         mol->AddAtom(neues);
    1090       }
    1091       No++;
    1092     }
    1093   }
    1094   file->close();
    1095   delete(file);
    1096 };
    1097 
    1098631/** Stores all elements of config structure from which they can be re-read.
    1099632 * \param *filename name of file
    … …  
    1123656    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    1124657    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    1125     *output << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
    1126     switch(Thermostats->Thermostat) {
    1127       default:
    1128       case None:
    1129         break;
    1130       case Woodcock:
    1131         *output << Thermostats->ScaleTempStep;
    1132         break;
    1133       case Gaussian:
    1134         *output << Thermostats->ScaleTempStep;
    1135         break;
    1136       case Langevin:
    1137         *output << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
    1138         break;
    1139       case Berendsen:
    1140         *output << Thermostats->TempFrequency;
    1141         break;
    1142       case NoseHoover:
    1143         *output << Thermostats->HooverMass;
    1144         break;
    1145     };
     658    *output << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
     659    *output << Thermostats->activeThermostat->writeParams();
    1146660    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    1147661    *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    … …  
    21281642  // read desired Thermostat from file along with needed additional parameters
    21291643  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    2130     if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
    2131       if (Thermostats->ThermostatImplemented[0] == 1) {
    2132         Thermostats->Thermostat = None;
    2133       } else {
    2134         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2135         Thermostats->Thermostat = None;
    2136       }
    2137     } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
    2138       if (Thermostats->ThermostatImplemented[1] == 1) {
    2139         Thermostats->Thermostat = Woodcock;
    2140         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
    2141       } else {
    2142         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2143         Thermostats->Thermostat = None;
    2144       }
    2145     } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
    2146       if (Thermostats->ThermostatImplemented[2] == 1) {
    2147         Thermostats->Thermostat = Gaussian;
    2148         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
    2149       } else {
    2150         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2151         Thermostats->Thermostat = None;
    2152       }
    2153     } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
    2154       if (Thermostats->ThermostatImplemented[3] == 1) {
    2155         Thermostats->Thermostat = Langevin;
    2156         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
    2157         if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
    2158           DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
    2159         } else {
    2160           Thermostats->alpha = 1.;
    2161         }
    2162       } else {
    2163         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2164         Thermostats->Thermostat = None;
    2165       }
    2166     } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
    2167       if (Thermostats->ThermostatImplemented[4] == 1) {
    2168         Thermostats->Thermostat = Berendsen;
    2169         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
    2170       } else {
    2171         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2172         Thermostats->Thermostat = None;
    2173       }
    2174     } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
    2175       if (Thermostats->ThermostatImplemented[5] == 1) {
    2176         Thermostats->Thermostat = NoseHoover;
    2177         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
    2178         Thermostats->alpha = 0.;
    2179       } else {
    2180         DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    2181         Thermostats->Thermostat = None;
    2182       }
    2183     } else {
    2184       DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    2185       Thermostats->Thermostat = None;
    2186     }
     1644    Thermostats->makeActive(thermo,fb);
    21871645  } else {
    21881646    if ((Thermostats->TargetTemp != 0))
    21891647      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    2190     Thermostats->Thermostat = None;
     1648    Thermostats->chooseNone();
    21911649  }
    21921650  delete[](thermo);

  • src/config.hpp

    rd17f3d r51c3e4  
    114114
    115115  int TestSyntax(const char * const filename, const periodentafel * const periode) const;
    116   void Load(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList);
    117   void LoadOld(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList);
    118116  bool Save(const char * const filename, const periodentafel * const periode, molecule * const mol) const;
    119117  bool SaveMPQC(const char * const filename, const molecule * const mol) const;

  • src/molecule.cpp

    rd17f3d r51c3e4  
    10551055    *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
    10561056  for (int step=startstep;step < endstep; step++) { // loop over all time steps
    1057     temperature = 0.;
    1058     ActOnAllAtoms( &TrajectoryParticle::AddKineticToTemperature, &temperature, step);
     1057    temperature = atoms.totalTemperatureAtStep(step);
    10591058    *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
    10601059  }

  • src/molecule.hpp

    rd17f3d r51c3e4  
    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

    rd17f3d r51c3e4  
    3131#include "LinearAlgebra/Plane.hpp"
    3232#include "ThermoStatContainer.hpp"
     33#include "Thermostats/Berendsen.hpp"
    3334
    3435#include <gsl/gsl_matrix.h>
    … …  
    511512
    512513  // check whether we have sufficient space in Trajectories for each atom
    513   ActOnAllAtoms( &atom::ResizeTrajectory, MaxSteps );
     514  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps));
    514515  // push endstep to last one
    515   ActOnAllAtoms( &atom::CopyStepOnStep, MaxSteps, endstep );
     516  for_each(atoms.begin(),atoms.end(),boost::bind(&atom::CopyStepOnStep,_1,MaxSteps,endstep));
    516517  endstep = MaxSteps;
    517518
    … …  
    622623  ActualTemp = 0.;
    623624  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
    624   Thermostats(configuration, ActualTemp, Berendsen);
     625  Berendsen berendsen = Berendsen();
     626  berendsen.addToContainer(configuration.Thermostats);
     627  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
     628  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
    625629  MDSteps++;
    626630
    … …  
    628632  return true;
    629633};
    630 
    631 /** Implementation of various thermostats.
    632  * All these thermostats apply an additional force which has the following forms:
    633  * -# Woodcock
    634  *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
    635  * -# Gaussian
    636  *  \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$
    637  * -# Langevin
    638  *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$
    639  * -# Berendsen
    640  *  \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$
    641  * -# Nose-Hoover
    642  *  \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$
    643  * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
    644  * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
    645  * belatedly and the constraint force set.
    646  * \param *P Problem at hand
    647  * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
    648  * \sa InitThermostat()
    649  */
    650 void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
    651 {
    652   double ekin = 0.;
    653   double E = 0., G = 0.;
    654   double delta_alpha = 0.;
    655   double ScaleTempFactor;
    656   gsl_rng * r;
    657   const gsl_rng_type * T;
    658 
    659   // calculate scale configuration
    660   ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
    661 
    662   // differentating between the various thermostats
    663   switch(Thermostat) {
    664      case None:
    665       DoLog(2) && (Log() << Verbose(2) <<  "Applying no thermostat..." << endl);
    666       break;
    667      case Woodcock:
    668       if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
    669         DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
    670         ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
    671       }
    672       break;
    673      case Gaussian:
    674       DoLog(2) && (Log() << Verbose(2) <<  "Applying Gaussian thermostat..." << endl);
    675       ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
    676 
    677       DoLog(1) && (Log() << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl);
    678       ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
    679 
    680       break;
    681      case Langevin:
    682       DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
    683       // init random number generator
    684       gsl_rng_env_setup();
    685       T = gsl_rng_default;
    686       r = gsl_rng_alloc (T);
    687       // Go through each ion
    688       ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
    689       break;
    690 
    691      case Berendsen:
    692       DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
    693       ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
    694       break;
    695 
    696      case NoseHoover:
    697       DoLog(2) && (Log() << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl);
    698       // dynamically evolve alpha (the additional degree of freedom)
    699       delta_alpha = 0.;
    700       ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    701       delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
    702       configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
    703       DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
    704       // apply updated alpha as additional force
    705       ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
    706       break;
    707   }
    708   DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
    709 };

  • src/molecule_fragmentation.cpp

    rd17f3d r51c3e4  
    810810  DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
    811811  if (file.good()) {
    812     ActOnAllAtoms( &atom::OutputOrder, &file );
     812    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputOrder), &file));
    813813    file.close();
    814814    DoLog(1) && (Log() << Verbose(1) << "done." << endl);

  • src/molecule_graph.cpp

    rd17f3d r51c3e4  
    220220
    221221    // output bonds for debugging (if bond chain list was correctly installed)
    222     ActOnAllAtoms( &atom::OutputBondOfAtom );
     222    for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
    223223  } else
    224224    DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << getAtomCount() << ", thus no bonds, no connections!." << endl);
    … …  
    638638{
    639639  DoLog(1) && (Log() << Verbose(1) << "Final graph info for each atom is:" << endl);
    640   ActOnAllAtoms( &atom::OutputGraphInfo );
     640  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputGraphInfo));
    641641}
    642642;
Note: See TracChangeset for help on using the changeset viewer.