(CW) AbInitioThermodynamics: Added option for diatomic molecules

AIMSTools/aimstools/AbInitioThermodynamics.py

@@ -268,6 +268,7 @@ class PhaseDiagram:
         colors=None,
         vibrational_data=None,
         pressure_unit="Pa",
+        diatomic=False,
     ):
         """The energy is assumed to be given per area."""
         # adsorption energy and coverage are reshape into column vectors
@@ -307,6 +308,8 @@ class PhaseDiagram:
             assert self.n_phases == len(colors), "Number of colors does not match number of phases"
             self.colors = colors
 
+        self.diatomic = diatomic
+
         self.grid_data = {
             "t_min": None,
             "t_max": None,
@@ -328,7 +331,11 @@ class PhaseDiagram:
             self.f_vib = [lambda x: 0] * self.n_phases # always return zero
 
     def getDeltaMu(self, temperature, pressure):
-        return calculateDeltaMu(temperature, pressure, *self.molecule_data)
+        if self.diatomic is False:
+            delta_mu = calculateDeltaMu(temperature, pressure, *self.molecule_data)
+        else:
+            delta_mu = calculateDeltaMu(temperature, pressure, *self.molecule_data) / 2
+        return delta_mu
 
     def getGammaVsMu(self, mu):
         """Calculate the Gibbs free energy of adsorption for each
@@ -404,7 +411,10 @@ class PhaseDiagram:
         mu_rot = calculateMuRot(temperature, *self.molecule_data)
 
         # need exception for T=0
+        if self.diatomic is False:
             pressure = c * np.exp((delta - mu_rot) / units.kB / temperature)
+        else:
+            pressure = c * np.exp((2 * delta - mu_rot) / units.kB / temperature)
 
         return pressure