Changeset f915e1

Timestamp:

Apr 18, 2008, 1:30:19 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

995ff3

Parents:

719746

Message:

function Thermostats() added

File:

• pcp/src/run.c (modified) (6 diffs)

pcp/src/run.c ¶

@@ -33 +33 @@
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_min.h>
+#include <gsl/gsl_randist.h>
 #include "mpi.h"
 #include "data.h"
@@ -1198 +1199 @@
  * \param *params additional arguments, here Problem at hand
  */
-double my_f(const gsl_vector *v, void *params)
+double StructOpt_f(const gsl_vector *v, void *params)
 {
   struct Problem *P = (struct Problem *)params;
@@ -1258 +1259 @@
 }
 
-void my_df(const gsl_vector *v, void *params, gsl_vector *df) 
+void StructOpt_df(const gsl_vector *v, void *params, gsl_vector *df) 
 {
   struct Problem *P = (struct Problem *)params;
@@ -1273 +1274 @@
 }
 
-void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df)
+void StructOpt_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df)
 {
-  *f = my_f(x, params);
-  my_df(x, params, df);
+  *f = StructOpt_f(x, params);
+  StructOpt_df(x, params, df);
 }
 
@@ -1312 +1313 @@
 
   /* Initialize method and iterate */
-  minex_func.f = &my_f;
-  minex_func.df = &my_df;
-  minex_func.fdf = &my_fdf;
+  minex_func.f = &StructOpt_f;
+  minex_func.df = &StructOpt_df;
+  minex_func.fdf = &StructOpt_fdf;
   minex_func.n = np;
   minex_func.params = (void *)P;
@@ -1341 +1342 @@
   gsl_vector_free(x);
   gsl_multimin_fdfminimizer_free (s);
+}
+
+/** 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 Thermostats() should be called after UpdateIonsU(), 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 Thermostats(struct Problem *P, enum thermostats i)
+{
+  struct FileData *Files = &P->Files;
+  struct Ions *I = &P->Ion;
+  int is, ia, d;
+  double *U;
+  double a, ekin = 0.;
+  double E = 0., F = 0.;
+  double delta_alpha = 0.;
+  const int delta_t = P->R.delta_t;
+  double ScaleTempFactor;
+  double sigma;
+  gsl_rng * r;
+  const gsl_rng_type * T;
+  
+  // calculate current temperature
+  CalculateEnergyIonsU(P); // Temperature now in I->ActualTemp
+  ScaleTempFactor = P->R.TargetTemp/I->ActualTemp;
+  //if ((P->Par.me == 0) && (I->ActualTemp < MYEPSILON)) fprintf(stderr,"Thermostat: (1) I->ActualTemp = %lg",I->ActualTemp);
+  if (Files->MeOutMes) fprintf(Files->TemperatureFile, "%d\t%lg",P->R.OuterStep, I->ActualTemp);
+    
+  // differentating between the various thermostats
+  switch(i) {
+     case None:
+      debug(P, "Applying no thermostat...");
+      break;
+     case Woodcock:
+      if ((P->R.ScaleTempStep > 0) && ((P->R.OuterStep-1) % P->R.ScaleTempStep == 0)) {
+        debug(P, "Applying Woodcock thermostat...");
+        for (is=0; is < I->Max_Types; is++) {
+          a = 0.5*I->I[is].IonMass;
+          for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+            U = &I->I[is].U[NDIM*ia];
+            if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+              for (d=0; d<NDIM; d++) {
+                U[d] *= sqrt(ScaleTempFactor);
+                ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
+              }
+          }
+        }
+      }
+      break;
+     case Gaussian:
+      debug(P, "Applying Gaussian thermostat...");
+      for (is=0; is < I->Max_Types; is++) { // sum up constraint constant
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              F += U[d] * I->I[is].FIon[d+NDIM*ia]; 
+              E += U[d]*U[d]*I->I[is].IonMass;
+            }
+        }
+      }
+      if (P->Call.out[ValueOut]) fprintf(stderr, "(%i) Gaussian Least Constraint constant is %lg\n", P->Par.me, F/E);
+      for (is=0; is < I->Max_Types; is++) { // apply constraint constant on constraint force and on velocities
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              I->I[is].FConstraint[d+NDIM*ia] = (F/E) * (U[d]*I->I[is].IonMass);
+              U[d] += delta_t/I->I[is].IonMass * (I->I[is].FConstraint[d+NDIM*ia]);
+              ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
+            }
+        }
+      }
+      break;
+     case Langevin:
+      debug(P, "Applying Langevin thermostat...");
+      // init random number generator
+      gsl_rng_env_setup();
+      T = gsl_rng_default;
+      r = gsl_rng_alloc (T);
+      // Go through each ion
+      for (is=0; is < I->Max_Types; is++) {
+        sigma  = sqrt(P->R.TargetTemp/I->I[is].IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          // throw a dice to determine whether it gets hit by a heat bath particle
+          if (((((rand()/(double)RAND_MAX))*P->R.TempFrequency) < 1.)) {  // (I->I[is].IMT[ia] == MoveIon) &&  even FixedIon moves, only not by other's forces
+            if (P->Par.me == 0) fprintf(stderr,"(%i) Particle %i,%i was hit (sigma %lg): %lg -> ", P->Par.me, is, ia, 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);
+            }
+            if (P->Par.me == 0) fprintf(stderr,"%lg\n", sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]));
+          }
+          for (d=0; d<NDIM; d++)
+            ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
+        }
+      }
+      break;
+     case Berendsen:
+      debug(P, "Applying Berendsen-VanGunsteren thermostat...");
+      for (is=0; is < I->Max_Types; is++) {
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              U[d] *= sqrt(1+(P->R.delta_t/P->R.TempFrequency)*(ScaleTempFactor-1));
+              ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
+            }
+        }
+      }
+      break;
+     case NoseHoover:
+      debug(P, "Applying Nose-Hoover thermostat...");
+      // dynamically evolve alpha (the additional degree of freedom)
+      delta_alpha = 0.;
+      for (is=0; is < I->Max_Types; is++) { // sum up constraint constant
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              delta_alpha += U[d]*U[d]*I->I[is].IonMass;
+            }
+        }
+      }
+      delta_alpha = (delta_alpha - (3.*I->Max_TotalIons+1.) * P->R.TargetTemp)/(P->R.HooverMass*Units2Electronmass);
+      P->R.alpha += delta_alpha*delta_t;
+      if (P->Par.me == 0) fprintf(stderr,"(%i) alpha = %lg * %i = %lg\n", P->Par.me, delta_alpha, delta_t, P->R.alpha);
+      // apply updated alpha as additional force
+      for (is=0; is < I->Max_Types; is++) { // apply constraint constant on constraint force and on velocities
+        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
+          U = &I->I[is].U[NDIM*ia];
+          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              I->I[is].FConstraint[d+NDIM*ia] = - P->R.alpha * (U[d] * I->I[is].IonMass);
+              U[d] += delta_t/I->I[is].IonMass * (I->I[is].FConstraint[d+NDIM*ia]);
+              ekin += (0.5*I->I[is].IonMass) * U[d]*U[d];
+            }
+        }
+      }
+      break;
+  }    
+  I->EKin = ekin;
+  I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);    
+  //if ((P->Par.me == 0) && (I->ActualTemp < MYEPSILON)) fprintf(stderr,"Thermostat: (2) I->ActualTemp = %lg",I->ActualTemp);
+  if (Files->MeOutMes) { fprintf(Files->TemperatureFile, "\t%lg\n", I->ActualTemp); fflush(Files->TemperatureFile); }
 }