Changeset 62f793 for src


Ignore:
Timestamp:
Sep 25, 2008, 5:57:19 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
aa5702
Parents:
16a52b
Message:

thermostats enumerator, necessary variables included in class config, molecule::Thermostats() and ..::Constrained
Potential()

Thermostat header definitions implemented. Can be parsed from the ESPACK config file into class config
ConstrainedPotential() calculates the transformation between two given configurations my minimsing a penalty func
tional of the distance to travel per atom. This was implemented due to Michael Griebel's talk during the MultiMat

Closing meeting in order to produce some visuals. It basically mimics a "Bain Transformation". But it is not yet
perfectly implemented.

Also, MD was corrected in VerletIntegration(). However, forces are still wrong with mpqc, as the kinetic energy increases dramatically during the MD simulation.

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r16a52b r62f793  
    10041004              SaveFlag = true;
    10051005              cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
    1006               if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration.Deltat, configuration.GetIsAngstroem(), configuration.DoConstrainedMD))
     1006              if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration))
    10071007                cout << Verbose(2) << "File not found." << endl;
    10081008              else {

  • src/config.cpp

    r16a52b r62f793  
    1818  configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    1919  configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     20  ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
     21  ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
     22  for (int j=0;j<MaxThermostats;j++)
     23    ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
     24  Thermostat = 4;
     25  alpha = 0.;
     26  ScaleTempStep = 25;
     27  TempFrequency = 2.5;
    2028  strcpy(mainname,"pcp");
    2129  strcpy(defaultpath,"not specified");
    … …  
    2331  configpath[0]='\0';
    2432  configname[0]='\0';
    25  
     33
     34  strcpy(ThermostatNames[0],"None");
     35  ThermostatImplemented[0] = 1;
     36  strcpy(ThermostatNames[1],"Woodcock");
     37  ThermostatImplemented[1] = 1;
     38  strcpy(ThermostatNames[2],"Gaussian");
     39  ThermostatImplemented[2] = 1;
     40  strcpy(ThermostatNames[3],"Langevin");
     41  ThermostatImplemented[3] = 1;
     42  strcpy(ThermostatNames[4],"Berendsen");
     43  ThermostatImplemented[4] = 1;
     44  strcpy(ThermostatNames[5],"NoseHoover");
     45  ThermostatImplemented[5] = 1;
     46
    2647  FastParsing = false;
    2748  ProcPEGamma=8;
    … …  
    96117  Free((void **)&configname, "config::~config: *configname");
    97118};
     119
     120/** Readin of Thermostat related values from parameter file.
     121 * \param *source parameter file
     122 */
     123void config::InitThermostats(ifstream *source)
     124{
     125  char *thermo = MallocString(12, "IonsInitRead: thermo");
     126  int verbose = 0;
     127
     128  // read desired Thermostat from file along with needed additional parameters
     129  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     130    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
     131      if (ThermostatImplemented[0] == 1) {
     132        Thermostat = None;
     133      } else {
     134        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     135        Thermostat = None;
     136      }
     137    } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
     138      if (ThermostatImplemented[1] == 1) {
     139        Thermostat = Woodcock;
     140        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
     141      } else {
     142        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     143        Thermostat = None;
     144      }
     145    } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
     146      if (ThermostatImplemented[2] == 1) {
     147        Thermostat = Gaussian;
     148        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
     149      } else {
     150        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     151        Thermostat = None;
     152      }
     153    } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
     154      if (ThermostatImplemented[3] == 1) {
     155        Thermostat = Langevin;
     156        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
     157        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
     158          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
     159        } else {
     160          alpha = 1.;
     161        }
     162      } else {
     163        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     164        Thermostat = None;
     165      }
     166    } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
     167      if (ThermostatImplemented[4] == 1) {
     168        Thermostat = Berendsen;
     169        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
     170      } else {
     171        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     172        Thermostat = None;
     173      }
     174    } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
     175      if (ThermostatImplemented[5] == 1) {
     176        Thermostat = NoseHoover;
     177        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
     178        alpha = 0.;
     179      } else {
     180        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     181        Thermostat = None;
     182      }
     183    } else {
     184      cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
     185      Thermostat = None;
     186    }
     187  } else {
     188    if ((MaxOuterStep > 0) && (TargetTemp != 0))
     189      cout << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
     190    Thermostat = None;
     191  }
     192  Free((void **)&thermo, "InitThermostats: thermo");
     193};
     194
    98195
    99196/** Displays menu for editing each entry of the config file.
    … …  
    465562  double value[3];
    466563 
     564  InitThermostats(file);
     565 
    467566  /* Namen einlesen */
    468567
    … …  
    10011100    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    10021101    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
     1102    *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
     1103    switch(Thermostat) {
     1104      default:
     1105      case None:
     1106        break;
     1107      case Woodcock:
     1108        *output << ScaleTempStep;
     1109        break;
     1110      case Gaussian:
     1111        *output << ScaleTempStep;
     1112        break;
     1113      case Langevin:
     1114        *output << TempFrequency << "\t" << alpha;
     1115        break;
     1116      case Berendsen:
     1117        *output << TempFrequency;
     1118        break;
     1119      case NoseHoover:
     1120        *output << HooverMass;
     1121        break;
     1122    };
     1123    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    10031124    *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    10041125    *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;

  • src/molecules.cpp

    r16a52b r62f793  
    10171017 * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
    10181018 * \param *out output stream for debugging
    1019  * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$V^{con}(x)\f$)
     1019 * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
    10201020 * \param startstep start configuration (MDStep in molecule::trajectories)
    10211021 * \param endstep end configuration (MDStep in molecule::trajectories)
    … …  
    14791479 * \param *out output stream for debugging
    14801480 * \param *file filename
     1481 * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
    14811482 * \param delta_t time step width in atomic units
    14821483 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    … …  
    14851486 * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
    14861487 */
    1487 bool molecule::VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained)
    1488 {
    1489   element *runner = elemente->start;
     1488bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
     1489{
    14901490  atom *walker = NULL;
    14911491  int AtomNo;
    … …  
    14931493  string token;
    14941494  stringstream item;
    1495   double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy;
     1495  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp = 0.;
    14961496  ForceMatrix Force;
    14971497
    … …  
    15231523      }
    15241524    // solve a constrained potential if we are meant to
    1525     if (DoConstrained) {
     1525    if (configuration.DoConstrainedMD) {
    15261526      // calculate forces and potential
    15271527      atom **PermutationMap = NULL;
    1528       ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap, DoConstrained, 0, IsAngstroem);
    1529       EvaluateConstrainedForces(out, DoConstrained, 0, PermutationMap, &Force);
     1528      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
     1529      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
    15301530      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
    15311531    }
    15321532   
    15331533    // and perform Verlet integration for each atom with position, velocity and force vector
    1534     runner = elemente->start;
    1535     while (runner->next != elemente->end) { // go through every element
    1536       runner = runner->next;
    1537       IonMass = runner->mass;
    1538       a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    1539       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
    1540         AtomNo = 0;
     1534    walker = start;
     1535    while (walker->next != end) { // go through every atom of this element
     1536      walker = walker->next;
     1537      a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
     1538      // check size of vectors
     1539      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
     1540        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
     1541        Trajectories[walker].R.resize(MDSteps+10);
     1542        Trajectories[walker].U.resize(MDSteps+10);
     1543        Trajectories[walker].F.resize(MDSteps+10);
     1544      }
     1545     
     1546      // Update R (and F)
     1547      for (int d=0; d<NDIM; d++) {
     1548        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     1549        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
     1550        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);
     1551        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
     1552      }
     1553      // Update U
     1554      for (int d=0; d<NDIM; d++) {
     1555        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
     1556        Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
     1557      }
     1558//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
     1559//      for (int d=0;d<NDIM;d++)
     1560//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
     1561//      out << ")\t(";
     1562//      for (int d=0;d<NDIM;d++)
     1563//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
     1564//      out << ")" << endl;
     1565            // next atom
     1566    }
     1567  }
     1568  // correct velocities (rather momenta) so that center of mass remains motionless
     1569  for(int d=0;d<NDIM;d++)
     1570    Vector[d] = 0.;
     1571  IonMass = 0.;
     1572  walker = start;
     1573  while (walker->next != end) { // go through every atom
     1574    walker = walker->next;
     1575    IonMass += walker->type->mass;  // sum up total mass
     1576    for(int d=0;d<NDIM;d++) {
     1577      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
     1578    }
     1579  }
     1580  walker = start;
     1581  while (walker->next != end) { // go through every atom of this element
     1582    walker = walker->next;
     1583    for(int d=0;d<NDIM;d++) {
     1584      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
     1585      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
     1586    }
     1587  }
     1588  Thermostats(configuration, ActualTemp, Berendsen);
     1589  MDSteps++;
     1590
     1591 
     1592  // exit
     1593  return true;
     1594};
     1595
     1596
     1597/** Implementation of various thermostats.
     1598 * All these thermostats apply an additional force which has the following forms:
     1599 * -# Woodcock
     1600 *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
     1601 * -# Gaussian
     1602 *  \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$
     1603 * -# Langevin
     1604 *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$ 
     1605 * -# Berendsen
     1606 *  \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$
     1607 * -# Nose-Hoover
     1608 *  \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$
     1609 * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
     1610 * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
     1611 * belatedly and the constraint force set.
     1612 * \param *P Problem at hand
     1613 * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
     1614 * \sa InitThermostat()
     1615 */
     1616void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
     1617{
     1618  double ekin = 0.;
     1619  double E = 0., G = 0.;
     1620  double delta_alpha = 0.;
     1621  double ScaleTempFactor;
     1622  double sigma;
     1623  double IonMass;
     1624  int d;
     1625  gsl_rng * r;
     1626  const gsl_rng_type * T;
     1627  double *U = NULL, *F = NULL, FConstraint[NDIM];
     1628  atom *walker = NULL;
     1629 
     1630  // calculate scale configuration
     1631  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
     1632   
     1633  // differentating between the various thermostats
     1634  switch(Thermostat) {
     1635     case None:
     1636      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
     1637      break;
     1638     case Woodcock:
     1639      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
     1640        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
    15411641        walker = start;
    15421642        while (walker->next != end) { // go through every atom of this element
    15431643          walker = walker->next;
    1544           if (walker->type == runner) { // if this atom fits to element
    1545             // check size of vectors
    1546             if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
    1547               //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    1548               Trajectories[walker].R.resize(MDSteps+10);
    1549               Trajectories[walker].U.resize(MDSteps+10);
    1550               Trajectories[walker].F.resize(MDSteps+10);
     1644          IonMass = walker->type->mass;
     1645          U = Trajectories[walker].U.at(MDSteps).x;
     1646          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1647            for (d=0; d<NDIM; d++) {
     1648              U[d] *= sqrt(ScaleTempFactor);
     1649              ekin += 0.5*IonMass * U[d]*U[d];
    15511650            }
    1552            
    1553             // Update R (and F)
    1554             for (int d=0; d<NDIM; d++) {
    1555               Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
    1556               Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
    1557               Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
    1558               Trajectories[walker].R.at(MDSteps).x[d] += delta_t*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
     1651        }
     1652      }
     1653      break;
     1654     case Gaussian:
     1655      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
     1656      walker = start;
     1657      while (walker->next != end) { // go through every atom of this element
     1658        walker = walker->next;
     1659        IonMass = walker->type->mass;
     1660        U = Trajectories[walker].U.at(MDSteps).x;
     1661        F = Trajectories[walker].F.at(MDSteps).x;
     1662        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1663          for (d=0; d<NDIM; d++) {
     1664            G += U[d] * F[d];
     1665            E += U[d]*U[d]*IonMass;
     1666          }
     1667      }
     1668      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
     1669      walker = start;
     1670      while (walker->next != end) { // go through every atom of this element
     1671        walker = walker->next;
     1672        IonMass = walker->type->mass;
     1673        U = Trajectories[walker].U.at(MDSteps).x;
     1674        F = Trajectories[walker].F.at(MDSteps).x;
     1675        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1676          for (d=0; d<NDIM; d++) {
     1677            FConstraint[d] = (G/E) * (U[d]*IonMass);
     1678            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1679            ekin += IonMass * U[d]*U[d];
     1680          }
     1681      }
     1682      break;
     1683     case Langevin:
     1684      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
     1685      // init random number generator
     1686      gsl_rng_env_setup();
     1687      T = gsl_rng_default;
     1688      r = gsl_rng_alloc (T);
     1689      // Go through each ion
     1690      walker = start;
     1691      while (walker->next != end) { // go through every atom of this element
     1692        walker = walker->next;
     1693        IonMass = walker->type->mass;
     1694        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
     1695        U = Trajectories[walker].U.at(MDSteps).x;
     1696        F = Trajectories[walker].F.at(MDSteps).x;
     1697        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1698          // throw a dice to determine whether it gets hit by a heat bath particle
     1699          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) { 
     1700            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
     1701            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
     1702            for (d=0; d<NDIM; d++) {
     1703              U[d] = gsl_ran_gaussian (r, sigma);
    15591704            }
    1560             // Update U
    1561             for (int d=0; d<NDIM; d++) {
    1562               Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
    1563               Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
    1564             }
    1565 //            out << "Integrated position&velocity of step " << (MDSteps) << ": (";
    1566 //            for (int d=0;d<NDIM;d++)
    1567 //              out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
    1568 //            out << ")\t(";
    1569 //            for (int d=0;d<NDIM;d++)
    1570 //              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
    1571 //            out << ")" << endl;
    1572             // next atom
    1573             AtomNo++;
     1705            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
     1706          }
     1707          for (d=0; d<NDIM; d++)
     1708            ekin += 0.5*IonMass * U[d]*U[d];
     1709        }
     1710      }
     1711      break;
     1712     case Berendsen:
     1713      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
     1714      walker = start;
     1715      while (walker->next != end) { // go through every atom of this element
     1716        walker = walker->next;
     1717        IonMass = walker->type->mass;
     1718        U = Trajectories[walker].U.at(MDSteps).x;
     1719        F = Trajectories[walker].F.at(MDSteps).x;
     1720        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1721          for (d=0; d<NDIM; d++) {
     1722            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
     1723            ekin += 0.5*IonMass * U[d]*U[d];
    15741724          }
    15751725        }
    15761726      }
    1577     }
    1578   }
    1579 //  // correct velocities (rather momenta) so that center of mass remains motionless
    1580 //  for(int d=0;d<NDIM;d++)
    1581 //    Vector[d] = 0.;
    1582 //  IonMass = 0.;
    1583 //  walker = start;
    1584 //  while (walker->next != end) { // go through every atom
    1585 //    walker = walker->next;
    1586 //    IonMass += walker->type->mass;  // sum up total mass
    1587 //    for(int d=0;d<NDIM;d++) {
    1588 //      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
    1589 //    }
    1590 //  }
    1591 //  walker = start;
    1592 //  while (walker->next != end) { // go through every atom of this element
    1593 //    walker = walker->next;
    1594 //    for(int d=0;d<NDIM;d++) {
    1595 //      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
    1596 //    }
    1597 //  }
    1598   MDSteps++;
    1599 
    1600  
    1601   // exit
    1602   return true;
     1727      break;
     1728     case NoseHoover:
     1729      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
     1730      // dynamically evolve alpha (the additional degree of freedom)
     1731      delta_alpha = 0.;
     1732      walker = start;
     1733      while (walker->next != end) { // go through every atom of this element
     1734        walker = walker->next;
     1735        IonMass = walker->type->mass;
     1736        U = Trajectories[walker].U.at(MDSteps).x;
     1737        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1738          for (d=0; d<NDIM; d++) {
     1739            delta_alpha += U[d]*U[d]*IonMass;
     1740          }
     1741        }
     1742      }
     1743      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     1744      configuration.alpha += delta_alpha*configuration.Deltat;
     1745      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
     1746      // apply updated alpha as additional force
     1747      walker = start;
     1748      while (walker->next != end) { // go through every atom of this element
     1749        walker = walker->next;
     1750        IonMass = walker->type->mass;
     1751        U = Trajectories[walker].U.at(MDSteps).x;
     1752        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1753          for (d=0; d<NDIM; d++) {
     1754              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
     1755              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1756              ekin += (0.5*IonMass) * U[d]*U[d];
     1757            }
     1758        }
     1759      }
     1760      break;
     1761  }   
     1762  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
    16031763};
    16041764

  • src/molecules.hpp

    r16a52b r62f793  
    1616#include <gsl/gsl_multimin.h>
    1717#include <gsl/gsl_vector.h>
     18#include <gsl/gsl_randist.h>
    1819
    1920// STL headers
    … …  
    192193};
    193194
     195
    194196ostream & operator << (ostream &ost, bond &b);
    195197
    196198class MoleculeLeafClass;
     199
     200
     201#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
     202enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
     203
    197204
    198205/** The complete molecule.
    … …  
    260267        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
    261268        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
    262         bool VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained);
    263         double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
     269       
     270        bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
     271  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
     272
     273  double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
    264274        double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
    265275        void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    … …  
    392402    int DoConstrainedMD;
    393403    int MaxOuterStep;
     404    int Thermostat;
     405    int *ThermostatImplemented;
     406    char **ThermostatNames;
     407    double TempFrequency;
     408    double alpha;
     409    double HooverMass;
     410    double TargetTemp;
     411    int ScaleTempStep;
    394412   
    395413    private:
    … …  
    415433    int OutVisStep;
    416434    int OutSrcStep;
    417     double TargetTemp;
    418     int ScaleTempStep;
    419435    int MaxPsiStep;
    420436    double EpsWannier;
    … …  
    463479  char *GetDefaultPath() const;
    464480  void SetDefaultPath(const char *path);
     481  void InitThermostats(ifstream *source);
    465482};
    466483
Note: See TracChangeset for help on using the changeset viewer.