Changeset 62f793 for src/molecules.cpp

Timestamp:

Sep 25, 2008, 5:57:19 PM (17 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.

File:

• src/molecules.cpp (modified) (5 diffs)

src/molecules.cpp

@@ -1017 +1017 @@
  * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
  * \param *out output stream for debugging
- * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$V^{con}(x)\f$)
+ * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
  * \param startstep start configuration (MDStep in molecule::trajectories)
  * \param endstep end configuration (MDStep in molecule::trajectories)
@@ -1479 +1479 @@
  * \param *out output stream for debugging
  * \param *file filename
+ * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
  * \param delta_t time step width in atomic units
  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
@@ -1485 +1486 @@
  * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
  */
-bool molecule::VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained)
-{
-  element *runner = elemente->start;
+bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
+{
   atom *walker = NULL;
   int AtomNo;
@@ -1493 +1493 @@
   string token;
   stringstream item;
-  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy;
+  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp = 0.;
   ForceMatrix Force;
 
@@ -1523 +1523 @@
       }
     // solve a constrained potential if we are meant to
-    if (DoConstrained) {
+    if (configuration.DoConstrainedMD) {
       // calculate forces and potential
       atom **PermutationMap = NULL;
-      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap, DoConstrained, 0, IsAngstroem);
-      EvaluateConstrainedForces(out, DoConstrained, 0, PermutationMap, &Force);
+      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
+      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
       Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     }
     
     // and perform Verlet integration for each atom with position, velocity and force vector
-    runner = elemente->start;
-    while (runner->next != elemente->end) { // go through every element
-      runner = runner->next;
-      IonMass = runner->mass;
-      a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
-        AtomNo = 0;
+    walker = start;
+    while (walker->next != end) { // go through every atom of this element
+      walker = walker->next;
+      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
+      // check size of vectors
+      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
+        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
+        Trajectories[walker].R.resize(MDSteps+10);
+        Trajectories[walker].U.resize(MDSteps+10);
+        Trajectories[walker].F.resize(MDSteps+10);
+      }
+      
+      // Update R (and F)
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
+        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);
+        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
+      }
+      // Update U
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
+        Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
+      }
+//      out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
+//      out << ")" << endl;
+            // next atom
+    }
+  }
+  // correct velocities (rather momenta) so that center of mass remains motionless
+  for(int d=0;d<NDIM;d++)
+    Vector[d] = 0.;
+  IonMass = 0.;
+  walker = start;
+  while (walker->next != end) { // go through every atom
+    walker = walker->next;
+    IonMass += walker->type->mass;  // sum up total mass
+    for(int d=0;d<NDIM;d++) {
+      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
+    }
+  }
+  walker = start;
+  while (walker->next != end) { // go through every atom of this element
+    walker = walker->next;
+    for(int d=0;d<NDIM;d++) {
+      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
+      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
+    }
+  }
+  Thermostats(configuration, ActualTemp, Berendsen);
+  MDSteps++;
+
+  
+  // exit
+  return true;
+};
+
+
+/** Implementation of various thermostats.
+ * All these thermostats apply an additional force which has the following forms:
+ * -# Woodcock
+ *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
+ * -# Gaussian
+ *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
+ * -# Langevin
+ *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$  
+ * -# Berendsen
+ *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
+ * -# Nose-Hoover
+ *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
+ * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
+ * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
+ * belatedly and the constraint force set.
+ * \param *P Problem at hand
+ * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
+ * \sa InitThermostat()
+ */
+void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
+{
+  double ekin = 0.;
+  double E = 0., G = 0.;
+  double delta_alpha = 0.;
+  double ScaleTempFactor;
+  double sigma;
+  double IonMass;
+  int d;
+  gsl_rng * r;
+  const gsl_rng_type * T;
+  double *U = NULL, *F = NULL, FConstraint[NDIM];
+  atom *walker = NULL;
+  
+  // calculate scale configuration
+  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
+    
+  // differentating between the various thermostats
+  switch(Thermostat) {
+     case None:
+      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
+      break;
+     case Woodcock:
+      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
+        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
         walker = start;
         while (walker->next != end) { // go through every atom of this element
           walker = walker->next;
-          if (walker->type == runner) { // if this atom fits to element
-            // check size of vectors
-            if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
-              //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-              Trajectories[walker].R.resize(MDSteps+10);
-              Trajectories[walker].U.resize(MDSteps+10);
-              Trajectories[walker].F.resize(MDSteps+10);
+          IonMass = walker->type->mass;
+          U = Trajectories[walker].U.at(MDSteps).x;
+          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              U[d] *= sqrt(ScaleTempFactor);
+              ekin += 0.5*IonMass * U[d]*U[d];
             }
-            
-            // Update R (and F)
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
-              Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
-              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
+        }
+      }
+      break;
+     case Gaussian:
+      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            G += U[d] * F[d]; 
+            E += U[d]*U[d]*IonMass;
+          }
+      }
+      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            FConstraint[d] = (G/E) * (U[d]*IonMass);
+            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+            ekin += IonMass * U[d]*U[d];
+          }
+      }
+      break;
+     case Langevin:
+      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
+      // init random number generator
+      gsl_rng_env_setup();
+      T = gsl_rng_default;
+      r = gsl_rng_alloc (T);
+      // Go through each ion
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          // throw a dice to determine whether it gets hit by a heat bath particle
+          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) {  
+            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
+            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
+            for (d=0; d<NDIM; d++) {
+              U[d] = gsl_ran_gaussian (r, sigma);
             }
-            // Update U
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
-              Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
-            }
-//            out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
-//            for (int d=0;d<NDIM;d++)
-//              out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
-//            out << ")\t(";
-//            for (int d=0;d<NDIM;d++)
-//              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
-//            out << ")" << endl;
-            // next atom
-            AtomNo++;
+            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+          }
+          for (d=0; d<NDIM; d++)
+            ekin += 0.5*IonMass * U[d]*U[d];
+        }
+      }
+      break;
+     case Berendsen:
+      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
+            ekin += 0.5*IonMass * U[d]*U[d];
           }
         }
       }
-    }
-  }
-//  // correct velocities (rather momenta) so that center of mass remains motionless
-//  for(int d=0;d<NDIM;d++)
-//    Vector[d] = 0.;
-//  IonMass = 0.;
-//  walker = start;
-//  while (walker->next != end) { // go through every atom
-//    walker = walker->next;
-//    IonMass += walker->type->mass;  // sum up total mass
-//    for(int d=0;d<NDIM;d++) {
-//      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
-//    }
-//  }
-//  walker = start;
-//  while (walker->next != end) { // go through every atom of this element
-//    walker = walker->next;
-//    for(int d=0;d<NDIM;d++) {
-//      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
-//    }
-//  }
-  MDSteps++;
-
-  
-  // exit
-  return true;
+      break;
+     case NoseHoover:
+      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
+      // dynamically evolve alpha (the additional degree of freedom)
+      delta_alpha = 0.;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            delta_alpha += U[d]*U[d]*IonMass;
+          }
+        }
+      }
+      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
+      configuration.alpha += delta_alpha*configuration.Deltat;
+      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
+      // apply updated alpha as additional force
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
+              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+              ekin += (0.5*IonMass) * U[d]*U[d];
+            }
+        }
+      }
+      break;
+  }    
+  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
 };