Merge pull request #161 from ICAMS/main

update branch

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.3.10
+current_version = 1.3.12
 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.12"
 
 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.12"
 
 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/postprocessing.py

@@ -1,6 +1,8 @@
 import os
 import numpy as np
 import yaml
+import matplotlib.pyplot as plt
+import warnings
 
 def read_report(folder):
     """
@@ -145,4 +147,83 @@ def gather_results(mainfolder):
                 datadict['error_code'][-1] = _extract_error(errfile)
 
     df = pd.DataFrame(data=datadict)
-    return df
+    return df
+
+def find_transition_temperature(folder1, folder2, fit_order=4, plot=True):
+    """
+    Find transition temperature where free energy of two phases are equal.
+
+    Parameters
+    ----------
+    folder1: string
+        directory with temperature scale calculation
+
+    folder2: string
+        directory with temperature scale calculation
+
+    fit_order: int, optional
+        default 4. Order for polynomial fit of temperature vs free energy
+    
+    plot: bool, optional
+        default True. Plot the results.
+    """
+    file1 = os.path.join(folder1, 'temperature_sweep.dat')
+    file2 = os.path.join(folder2, 'temperature_sweep.dat')
+    if not os.path.exists(file1):
+        raise FileNotFoundError(f'{file1} does not exist')
+    if not os.path.exists(file2):
+        raise FileNotFoundError(f'{file2} does not exist')
+
+    t1, f1 = np.loadtxt(file1, unpack=True, usecols=(0,1))
+    t2, f2 = np.loadtxt(file2, unpack=True, usecols=(0,1))
+
+    #do some fitting to determine temps
+    t1min = np.min(t1)
+    t2min = np.min(t2)
+    t1max = np.max(t1)
+    t2max = np.max(t2)
+
+    tmin = np.min([t1min, t2min])
+    tmax = np.max([t1max, t2max])
+
+    #warn about extrapolation
+    if not t1min == t2min:
+        warnings.warn(f'free energy is being extrapolated!')
+    if not t1max == t2max:
+        warnings.warn(f'free energy is being extrapolated!')
+
+    #now fit
+    f1fit = np.polyfit(t1, f1, fit_order)
+    f2fit = np.polyfit(t2, f2, fit_order)
+
+    #reevaluate over the new range
+    fit_t = np.arange(tmin, tmax+1, 1)
+    fit_f1 = np.polyval(f1fit, fit_t)
+    fit_f2 = np.polyval(f2fit, fit_t)
+
+    #now evaluate the intersection temp
+    arg = np.argsort(np.abs(fit_f1-fit_f2))[0]
+    transition_temp = fit_t[arg]
+
+    #warn if the temperature is shady
+    if np.abs(transition_temp-tmin) < 1E-3:
+        warnings.warn('It is likely there is no intersection of free energies')
+    elif np.abs(transition_temp-tmax) < 1E-3:
+        warnings.warn('It is likely there is no intersection of free energies')
+
+    #plot
+    if plot:
+        c1lo = '#ef9a9a'
+        c1hi = '#b71c1c'
+        c2lo = '#90caf9'
+        c2hi = '#0d47a1'
+
+        plt.plot(fit_t, fit_f1, color=c1lo, label=f'{folder1} fit')
+        plt.plot(fit_t, fit_f2, color=c2lo, label=f'{folder2} fit')
+        plt.plot(t1, f1, color=c1hi, label=folder1, ls='dashed')
+        plt.plot(t2, f2, color=c2hi, label=folder2, ls='dashed')
+        plt.axvline(transition_temp, ls='dashed', c='#37474f')
+        plt.ylabel('Free energy (eV/atom)')
+        plt.xlabel('Temperature (K)')
+        plt.legend(frameon=False)
+    return transition_temp

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/index.rst

@@ -47,6 +47,7 @@ publication <https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMateri
    Using calphy <running_calphy/running_calphy>
    inputfile
    examples
+   research
    calphy
    faqs
    prologue/extending

docs/source/inputfile.md

@@ -546,6 +546,59 @@ folder_prefix: set1
 
 Prefix string to be added to folder names for calculation. Folders for calculations in calphy are named as `mode-lattice-temperature-pressure`. Therefore, if more than one calculation is run with the same parameters, they will be overwritten. To prevent this, `folder_prefix` can be used. If `folder_prefix` is provided, the folders will be named as `folder_prefix-mode-lattice-temperature-pressure`.
 
+
+---
+
+(script_mode)=
+#### `script_mode`        
+
+_type_: bool \
+_default_: False \
+_example_:
+```
+script_mode: False
+```  
+
+If True, a LAMMPS executable script is written and executed instead of the library interface of LAMMPS.
+Works only with `reference_phase: solid`, and `mode: fe`.
+Needs specification of [`lammps_executable`](lammps_executable) and [`mpi_executable`](mpi_executable).
+
+
+---
+
+(lammps_executable)=
+#### `lammps_executable`        
+
+_type_: string \
+_default_: None \
+_example_:
+```
+lammps_executable: lmp_mpi
+```  
+
+LAMMPS executable to run the calculations with.
+Works only with `reference_phase: solid`, and `mode: fe`.
+Works only if [`script_mode`](script_mode) is `True`.
+
+
+
+---
+
+(mpi_executable)=
+#### `mpi_executable`        
+
+_type_: string \
+_default_: None \
+_example_:
+```
+mpi_executable: mpiexec
+```  
+
+MPI executable to run the LAMMPS with.
+Works only with `reference_phase: solid`, and `mode: fe`.
+Works only if [`script_mode`](script_mode) is `True`.
+
+
 ---
 ---
 

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/).