diff --git a/.bumpversion.cfg b/.bumpversion.cfg index a22973e..5553454 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.3.10 +current_version = 1.3.12 commit = True tag = True diff --git a/calphy/__init__.py b/calphy/__init__.py index d8bc8ea..7e57954 100644 --- a/calphy/__init__.py +++ b/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 diff --git a/calphy/input.py b/calphy/input.py index 58020d6..749e4ef 100644 --- a/calphy/input.py +++ 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]): diff --git a/calphy/integrators.py b/calphy/integrators.py index e89751e..4b34324 100644 --- a/calphy/integrators.py +++ b/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,23 +470,18 @@ 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 @@ -493,40 +489,77 @@ def get_einstein_crystal_fe( 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 diff --git a/calphy/phase.py b/calphy/phase.py index a510c05..569d4db 100644 --- a/calphy/phase.py +++ b/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" diff --git a/calphy/postprocessing.py b/calphy/postprocessing.py index 6d82e3d..25936ca 100644 --- a/calphy/postprocessing.py +++ b/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 \ No newline at end of file + 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 \ No newline at end of file diff --git a/calphy/solid.py b/calphy/solid.py index 053a143..3f5b442 100644 --- a/calphy/solid.py +++ b/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 diff --git a/docs/source/index.rst b/docs/source/index.rst index 2dce4b3..ee0601d 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -47,6 +47,7 @@ publication inputfile examples + research calphy faqs prologue/extending diff --git a/docs/source/inputfile.md b/docs/source/inputfile.md index aa5d01d..af80be9 100644 --- a/docs/source/inputfile.md +++ b/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`. + + --- --- diff --git a/docs/source/prologue/acknowledgements.md b/docs/source/prologue/acknowledgements.md index 90c35e1..67435cc 100644 --- a/docs/source/prologue/acknowledgements.md +++ b/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/). diff --git a/docs/source/research.md b/docs/source/research.md new file mode 100644 index 0000000..f54b75a --- /dev/null +++ b/docs/source/research.md @@ -0,0 +1,16 @@ +# Research using `calphy` + +The following research works employed `calphy`: + +| Year | Material system | Calculation | Interatomic Potential | Publication | +|------ |----------------- |--------------------- |----------------------- |------------- | +| 2024 | MgAl, MgCa, AlCa | phase diagram | MTP | [Poul et. al.](https://www.researchsquare.com/article/rs-4732459/v1) | +| 2024 | Mo, Si, Mo$_3$Si, Mo$_5$Si$_3$, MoSi$_2$ | melting temperature | ACE | [Lenchuk et. al.](https://doi.org/10.1002/adem.202302043) | +| 2024 | AlLi | phase diagram | EAM, HDNNP, ACE | [Menon et. al.](https://arxiv.org/abs/2403.05724) | +| 2024 | CuZr | phase diagram | ACE | [Leimeroth et. al.](https://doi.org/10.1103/PhysRevMaterials.8.043602) | +| 2024 | SiO | phase diagram | ACE | [Erhard et. al.](https://doi.org/10.1038/s41467-024-45840-9) | +| 2023 | Mg | phase diagram | ACE | [Ibrahim et. al.](https://doi.org/10.1103/PhysRevMaterials.7.113801) | +| 2022 | ZnO | free energies | ML | [Goniakowski et. al.](https://doi.org/10.1021/acs.jpcc.2c06341) | +| 2019 | Ti, Si | phase diagram | EAM, SW | [Menon et. al.](https://doi.org/10.1103/PhysRevMaterials.5.103801) | + +Missing works? Reach out to [us](mailto:s.menon@mpie.de). \ No newline at end of file diff --git a/setup.py b/setup.py index 25c8b93..bdc761a 100644 --- a/setup.py +++ b/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.12', zip_safe=False, entry_points={ 'console_scripts': [