From bb4efaa0142b703cc19006e540b5089f807d4163 Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Wed, 23 Oct 2024 16:58:57 +0200 Subject: [PATCH 1/9] update einstein --- calphy/integrators.py | 67 ++++++++++++++++++++++++++++--------------- 1 file changed, 44 insertions(+), 23 deletions(-) diff --git a/calphy/integrators.py b/calphy/integrators.py index e89751e..95e4df0 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -43,7 +43,7 @@ kbJ = const.physical_constants["Boltzmann constant"][0] Na = const.physical_constants["Avogadro constant"][0] eV2J = const.eV - +J2eV = 6.242E18 #-------------------------------------------------------------------- # TI PATH INTEGRATION ROUTINES @@ -499,34 +499,55 @@ def get_einstein_crystal_fe( free energy of Einstein crystal """ - #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 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]) - - #convert k from ev/A2 to J/m2 - k = np.array(k)*(eV2J/1E-20) - omega = np.sqrt(k/mass) #convert a to m3 vol = vol*1E-30 - F_harm = 0 - F_cm = 0 - - 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 + #whats the beta + beta = (1/(kb*temp)) + beta = beta*eV2J + + #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) + + #evaluate DeBroglie wavelength of each atom + gamma_sqrd = (beta*h**2)/(2*np.pi*mass) + + #fe of Einstein crystal + F_e = np.log(((beta*k*gamma_sqrd)/(2*np.pi))**1.5) + F_e = np.sum(F_e)*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/natoms + F_cm = np.log(prefactor*(beta/(2*np.pi*mu2_over_k_sum))**1.5) + F_cm = F_cm*J2eV #convert to eV + else: + F_cm = 0 + + F_tot = F_e - F_cm + return F_tot #-------------------------------------------------------------------- # REF. STATE ROUTINES: LIQUID From a6f1f43e3db1a115ef099aa8df743b809a66ddb6 Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 00:02:22 +0200 Subject: [PATCH 2/9] update fe --- calphy/integrators.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/calphy/integrators.py b/calphy/integrators.py index 95e4df0..ae9cb66 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -38,6 +38,7 @@ #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] @@ -506,8 +507,7 @@ def get_einstein_crystal_fe( vol = vol*1E-30 #whats the beta - beta = (1/(kb*temp)) - beta = beta*eV2J + beta = (1/(kbJ*temp)) #create an array of mass mass = [] @@ -515,6 +515,7 @@ def get_einstein_crystal_fe( 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 @@ -531,8 +532,9 @@ def get_einstein_crystal_fe( gamma_sqrd = (beta*h**2)/(2*np.pi*mass) #fe of Einstein crystal - F_e = np.log(((beta*k*gamma_sqrd)/(2*np.pi))**1.5) - F_e = np.sum(F_e)*J2eV #convert back to eV + Z_e = ((beta**2*k*h**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: @@ -540,9 +542,9 @@ def get_einstein_crystal_fe( mu = mass/mass_sum mu2_over_k = mu**2/k mu2_over_k_sum = np.sum(mu2_over_k) - prefactor = vol/natoms + prefactor = vol F_cm = np.log(prefactor*(beta/(2*np.pi*mu2_over_k_sum))**1.5) - F_cm = F_cm*J2eV #convert to eV + F_cm = kb*temp*F_cm/natoms #convert to eV else: F_cm = 0 From 2c9a3739a0aadaba4b964cd8a1c821d7d50b7e1f Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 00:05:12 +0200 Subject: [PATCH 3/9] add individual contributions to report --- calphy/integrators.py | 5 ++++- calphy/solid.py | 9 ++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/calphy/integrators.py b/calphy/integrators.py index ae9cb66..86730a8 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -470,7 +470,8 @@ def get_einstein_crystal_fe( calc, vol, k, - cm_correction=True): + cm_correction=True, + return_contributions=False): """ Get the free energy of einstein crystal @@ -549,6 +550,8 @@ def get_einstein_crystal_fe( F_cm = 0 F_tot = F_e - F_cm + if return_contributions: + return F_e, -F_cm return F_tot #-------------------------------------------------------------------- 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 From 3c4ad96bbfe78313ea841d928dcb26fa80b60c59 Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 00:06:51 +0200 Subject: [PATCH 4/9] update reports --- calphy/phase.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/calphy/phase.py b/calphy/phase.py index 426fda9..294d9c0 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" From 91c992529ac6ec4c2a295e15517720d81e67313f Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 09:08:28 +0200 Subject: [PATCH 5/9] update temp --- calphy/integrators.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/calphy/integrators.py b/calphy/integrators.py index 86730a8..b8c86d8 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -501,6 +501,9 @@ def get_einstein_crystal_fe( free energy of Einstein crystal """ + #temperature + temp = calc._temperature + #natoms natoms = np.sum([calc._element_dict[x]['count'] for x in calc.element]) From d3c3e7d476989e3f0dc815699af4370275868156 Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 09:31:43 +0200 Subject: [PATCH 6/9] update units for h --- calphy/integrators.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/calphy/integrators.py b/calphy/integrators.py index b8c86d8..68b258e 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -532,11 +532,8 @@ def get_einstein_crystal_fe( #convert k from ev/A2 to J/m2 k = k*(eV2J/1E-20) - #evaluate DeBroglie wavelength of each atom - gamma_sqrd = (beta*h**2)/(2*np.pi*mass) - #fe of Einstein crystal - Z_e = ((beta**2*k*h**2)/(4*np.pi**2*mass))**1.5 + 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 From eccda27e1b7042a0335f52ec8131cf0531df6dc3 Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 10:38:14 +0200 Subject: [PATCH 7/9] update docs of Einstein free energy fn --- calphy/integrators.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/calphy/integrators.py b/calphy/integrators.py index 68b258e..4b34324 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -477,17 +477,11 @@ def get_einstein_crystal_fe( Parameters ---------- - temp : temperature, float - units - K + calc : Calculation object + contains all input parameters - natoms : int - no of atoms in the system - - mass : float - units - g/mol - - a : lattice constant, float - units - Angstrom + vol : float + converged volume per atom k : spring constant, float units - eV/Angstrom^2 @@ -495,10 +489,23 @@ 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. """ #temperature From 1af12a96b89148f8424b6d40adc2b5884c21c19a Mon Sep 17 00:00:00 2001 From: Sarath Menon Date: Thu, 24 Oct 2024 10:41:01 +0200 Subject: [PATCH 8/9] update docs --- docs/source/prologue/acknowledgements.md | 5 +++++ 1 file changed, 5 insertions(+) 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/). From 5593b904c5b277df2ff35264e1fab8787b8ec303 Mon Sep 17 00:00:00 2001 From: Sarath Date: Thu, 24 Oct 2024 10:43:21 +0200 Subject: [PATCH 9/9] =?UTF-8?q?Bump=20version:=201.3.10=20=E2=86=92=201.3.?= =?UTF-8?q?11?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- calphy/__init__.py | 2 +- calphy/input.py | 2 +- setup.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index a22973e..ef99446 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.3.10 +current_version = 1.3.11 commit = True tag = True diff --git a/calphy/__init__.py b/calphy/__init__.py index d8bc8ea..3806d22 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.11" def addtest(a,b): return a+b diff --git a/calphy/input.py b/calphy/input.py index 2d0f6ba..c0d8804 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.11" def _check_equal(val): if not (val[0]==val[1]==val[2]): diff --git a/setup.py b/setup.py index 25c8b93..27ff91c 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.11', zip_safe=False, entry_points={ 'console_scripts': [