Merge pull request #159 from ICAMS/update_cm_correction

Update cm correction

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.3.10
+current_version = 1.3.11
 commit = True
 tag = True
 

calphy/__init__.py

@@ -4,7 +4,7 @@
 from calphy.alchemy import Alchemy
 from calphy.routines import MeltingTemp
 
-__version__ = "1.3.10"
+__version__ = "1.3.11"
 
 def addtest(a,b):
     return a+b

calphy/input.py

@@ -40,7 +40,7 @@
 from ase.io import read, write
 import shutil
 
-__version__ = "1.3.10"
+__version__ = "1.3.11"
 
 def _check_equal(val):
     if not (val[0]==val[1]==val[2]):

calphy/integrators.py

@@ -38,12 +38,13 @@
 
 #Constants
 h = const.physical_constants["Planck constant in eV/Hz"][0]
+hJ = const.physical_constants["Planck constant"][0]
 hbar = h/(2*np.pi)
 kb = const.physical_constants["Boltzmann constant in eV/K"][0]
 kbJ = const.physical_constants["Boltzmann constant"][0]
 Na = const.physical_constants["Avogadro constant"][0]
 eV2J = const.eV
-
+J2eV = 6.242E18
 
 #--------------------------------------------------------------------
 #             TI PATH INTEGRATION ROUTINES
@@ -469,64 +470,96 @@ def get_einstein_crystal_fe(
     calc,
     vol, 
     k, 
-    cm_correction=True):
+    cm_correction=True,
+    return_contributions=False):
     """
     Get the free energy of einstein crystal
 
     Parameters
     ----------
-    temp : temperature, float
-        units - K
-
-    natoms : int
-        no of atoms in the system
-
-    mass : float
-        units - g/mol
+    calc : Calculation object
+        contains all input parameters
 
-    a : lattice constant, float
-        units - Angstrom
+    vol : float
+        converged volume per atom
 
     k : spring constant, float
         units - eV/Angstrom^2
 
     cm_correction : bool, optional, default - True
         add the centre of mass correction to free energy
 
+    return_contributions: bool, optional, default - True
+        If True, return individual contributions to the reference free energy.
+
     Returns
     -------
-    fe : float
-        free energy of Einstein crystal
+    F_tot : float
+        total free energy of reference crystal
+
+    F_e : float
+        Free energy of Einstein crystal without centre of mass correction. Only if `return_contributions` is True.
+
+    F_cm : float
+        centre of mass correction. Only if `return_contributions` is True.
+
+    Notes
+    -----
+    The equations for free energy of Einstein crystal and centre of mass correction are from https://doi.org/10.1063/5.0044833.
 
     """
-    #convert mass first for single particle in kg
-    mass = np.array([calc._element_dict[x]['mass'] for x in calc.element])
-    mass = (mass/Na)*1E-3
-    natoms = np.sum([calc._element_dict[x]['count'] for x in calc.element])
-    concentration = np.array([calc._element_dict[x]['composition'] for x in calc.element])
+    #temperature
+    temp = calc._temperature
 
-    #convert k from ev/A2 to J/m2
-    k = np.array(k)*(eV2J/1E-20)
-    omega = np.sqrt(k/mass)
+    #natoms
+    natoms = np.sum([calc._element_dict[x]['count'] for x in calc.element])
 
     #convert a to m3
     vol = vol*1E-30
 
-    F_harm = 0
-    F_cm = 0
+    #whats the beta
+    beta = (1/(kbJ*temp))   
 
-    for count, om in enumerate(omega):
-        if concentration[count] > 0:
-            F_harm += concentration[count]*np.log((hbar*om)/(kb*calc._temperature))
-            if cm_correction:
-                F_cm += np.log((natoms*concentration[count]/vol)*(2*np.pi*kbJ*calc._temperature/(natoms*concentration[count]*k[count]))**1.5)
-            #F_cm = 0
-    F_harm = 3*kb*calc._temperature*F_harm
-    F_cm = (kb*calc._temperature/natoms)*F_cm
-
-    F_harm = F_harm + F_cm
+    #create an array of mass
+    mass = []
+    for x in calc.element:
+        for count in range(calc._element_dict[x]['count']):
+            mass.append(calc._element_dict[x]['mass'])
+    mass = np.array(mass)
+
+    #convert mass to kg
+    mass = (mass/Na)*1E-3
 
-    return F_harm
+    #create an array of k as well
+    karr = []
+    for c, x in enumerate(calc.element):
+        for count in range(calc._element_dict[x]['count']):
+            karr.append(k[c])
+    k = np.array(karr)
+    #convert k from ev/A2 to J/m2
+    k = k*(eV2J/1E-20)
+
+    #fe of Einstein crystal
+    Z_e = ((beta**2*k*hJ**2)/(4*np.pi**2*mass))**1.5
+    F_e = np.log(Z_e)
+    F_e = kb*temp*np.sum(F_e)/natoms #*J2eV #convert back to eV
+
+    #now get the cm correction
+    if cm_correction:
+        mass_sum = np.sum(mass)
+        mu = mass/mass_sum
+        mu2_over_k = mu**2/k
+        mu2_over_k_sum = np.sum(mu2_over_k) 
+        prefactor = vol
+        F_cm = np.log(prefactor*(beta/(2*np.pi*mu2_over_k_sum))**1.5)
+        F_cm = kb*temp*F_cm/natoms #convert to eV
+    else:
+        F_cm = 0
+
+    F_tot = F_e - F_cm
+    if return_contributions:
+        return F_e, -F_cm
+    return F_tot
 
 #--------------------------------------------------------------------
 #             REF. STATE ROUTINES: LIQUID

calphy/phase.py

@@ -163,6 +163,8 @@ def __init__(self, calculation=None, simfolder=None, log_to_screen=False):
 
         self.ferr = 0
         self.fref = 0
+        self.feinstein = 0
+        self.fcm = 0
         self.fideal = 0
 
         self.w = 0
@@ -682,6 +684,8 @@ def submit_report(self, extra_dict=None):
         report["results"]["free_energy"] = float(self.fe)
         report["results"]["error"] = float(self.ferr)
         report["results"]["reference_system"] = float(self.fref)
+        report["results"]["einstein_crystal"] = float(self.feinstein)
+        report["results"]["com_correction"] = float(self.fcm)
         report["results"]["work"] = float(self.w)
         report["results"]["pv"] = float(self.pv)
         report["results"]["unit"] = "eV/atom"

calphy/solid.py

@@ -516,17 +516,20 @@ def thermodynamic_integration(self):
         Calculates the final work, energy dissipation and free energy by
         matching with Einstein crystal
         """
-        f1 = get_einstein_crystal_fe(
+        fe, fcm = get_einstein_crystal_fe(
             self.calc,  
             self.vol, 
-            self.k)
+            self.k,
+            return_contributions=True)
 
         w, q, qerr = find_w(self.simfolder, 
             self.calc,
             full=True, 
             solid=True)
 
-        self.fref = f1
+        self.fref = fe + fcm
+        self.feinstein = fe
+        self.fcm = fcm
         self.w = w
         self.ferr = qerr
 

docs/source/prologue/acknowledgements.md

@@ -24,6 +24,11 @@ We acknowledge the following people for their contribution to calphy:
 
 - Abril Azócar Guzmán for the design of calphy logo
 
+- Marvin Poul for numerous fixes, improvements, and for the help in fixing centre of mass corrections for multiple species
+
+- Sebastian Havens for the singularity recipe
+
+
 The development of this module was started at the [Interdisciplinary Centre for Advanced
 Materials Simulation](http://www.icams.de/content), at the [Ruhr
 University Bochum](https://www.ruhr-uni-bochum.de/en), Germany. Current development is carried out at the [Max-Planck-Institut für Eisenforschung GmbH](https://www.mpie.de/).

setup.py

@@ -53,7 +53,7 @@
     packages=find_packages(include=['calphy', 'calphy.*']),
     test_suite='tests',
     url='https://github.com/ICAMS/calphy',
-    version='1.3.10',
+    version='1.3.11',
     zip_safe=False,
     entry_points={
         'console_scripts': [