diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..8609886 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,3 @@ +*.1relax linguist-detectable=false +*.2relax linguist-detectable=false +*.3static linguist-detectable=false diff --git a/LICENSE b/LICENSE deleted file mode 100644 index 5533d21..0000000 --- a/LICENSE +++ /dev/null @@ -1,20 +0,0 @@ -The MIT License (MIT) -Copyright (c) 2011-2012 MIT & The Regents of the University of California, -through Lawrence Berkeley National Laboratory - -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the "Software"), to deal in -the Software without restriction, including without limitation the rights to -use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of -the Software, and to permit persons to whom the Software is furnished to do so, -subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS -FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR -COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER -IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN -CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..59ac984 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Nigel Lee En Hew + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/examples/codespace/Fe3Pt.ipynb b/examples/codespace/Fe3Pt.ipynb index d60bf9b..a11dead 100644 --- a/examples/codespace/Fe3Pt.ipynb +++ b/examples/codespace/Fe3Pt.ipynb @@ -19,6 +19,24 @@ "4. You are now ready to run the cells below!" ] }, + { + "cell_type": "markdown", + "id": "9ae53c39", + "metadata": {}, + "source": [ + "# Install pyzentropy in editable mode" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "124181fd", + "metadata": {}, + "outputs": [], + "source": [ + "!cd ../../ && pip install -e ." + ] + }, { "cell_type": "markdown", "id": "8de285bd", @@ -85,20 +103,7 @@ "source": [ "# PyZentropy Configuration Objects\n", "\n", - "This section constructs PyZentropy `Configuration` objects for each Fe₃Pt configuration. Properties at 0 K are excluded to avoid division by zero errors. For each configuration, internal energies and partition functions are calculated." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "eba55270", - "metadata": {}, - "outputs": [], - "source": [ - "# Extract basic properties from the reference configuration\n", - "number_of_atoms = input[\"FM\"]['number_of_atoms']\n", - "volumes = input[\"FM\"]['volumes']\n", - "temperatures = input[\"FM\"]['temperatures'] " + "This section constructs PyZentropy `Configuration` objects for each Fe₃Pt configuration. For each configuration, the internal energies are calculated automatically." ] }, { @@ -118,38 +123,19 @@ " config_objects[config_name] = Configuration(\n", " name=config_name,\n", " multiplicity=config['multiplicity'],\n", - " number_of_atoms=number_of_atoms,\n", - " volumes=volumes,\n", - " temperatures=temperatures,\n", + " number_of_atoms=config['number_of_atoms'],\n", + " volumes=config['volumes'],\n", + " temperatures=config['temperatures'],\n", " # Extract Helmholtz energies and their derivatives from DFTTK results\n", - " helmholtz_energies=config['helmholtz_energy'],\n", - " helmholtz_energies_dV=config['helmholtz_energy_dV'],\n", - " helmholtz_energies_d2V2=config['helmholtz_energy_d2V2'],\n", + " helmholtz_energies=config['helmholtz_energies'],\n", + " helmholtz_energies_dV=config['helmholtz_energies_dV'],\n", + " helmholtz_energies_d2V2=config['helmholtz_energies_d2V2'],\n", " # Extract entropy and heat capacity data\n", - " entropies=config['entropy'],\n", - " heat_capacities=config['heat_capacity']\n", + " entropies=config['entropies'],\n", + " heat_capacities=config['heat_capacities']\n", " )" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "3e1b202a", - "metadata": {}, - "outputs": [], - "source": [ - "# Calculate internal energies for each configuration\n", - "for config in config_objects.values():\n", - " config.calculate_internal_energies()\n", - "\n", - "# Use the ground state configuration as reference for partition function calculations\n", - "reference_helmholtz_energy = config_objects[\"FM\"].helmholtz_energies\n", - "\n", - "# Calculate partition functions for each configuration using the reference energy\n", - "for config in config_objects.values():\n", - " config.calculate_partition_functions(reference_helmholtz_energy)" - ] - }, { "cell_type": "markdown", "id": "2979bf50", @@ -177,7 +163,7 @@ "metadata": {}, "outputs": [], "source": [ - "config_objects[\"FM\"].plot(\"helmholtz_energy_vs_volume\")" + "config_objects[\"FM\"].plot_vt(\"helmholtz_energy_vs_volume\")" ] }, { @@ -190,7 +176,7 @@ "# You can also plot for selected temperatures by passing an array of temperatures.\n", "# The code will automatically find and use the closest available temperatures.\n", "selected_temperatures = np.array([100, 300, 600, 900])\n", - "config_objects[\"FM\"].plot(\"helmholtz_energy_vs_volume\", selected_temperatures=selected_temperatures)" + "config_objects[\"FM\"].plot_vt(\"helmholtz_energy_vs_volume\", selected_temperatures=selected_temperatures)" ] }, { @@ -202,8 +188,8 @@ "source": [ "# You can also plot for selected volumes by passing an array of volumes.\n", "# The code will automatically find and use the closest available volumes.\n", - "selected_volumes = np.array([100, 120, 140, 160])\n", - "config_objects[\"FM\"].plot(\"helmholtz_energy_vs_temperature\", selected_volumes=selected_volumes)" + "selected_volumes = np.array([136, 140, 160])\n", + "config_objects[\"FM\"].plot_vt(\"helmholtz_energy_vs_temperature\", selected_volumes=selected_volumes)" ] }, { @@ -224,17 +210,7 @@ "outputs": [], "source": [ "# Create a System object from the configuration objects\n", - "system = System(config_objects)\n", - "\n", - "# Perform thermodynamic calculations for the system\n", - "system.calculate_partition_functions()\n", - "system.calculate_probabilities()\n", - "system.calculate_helmholtz_energies(reference_helmholtz_energy)\n", - "system.calculate_entropies()\n", - "system.calculate_helmholtz_energies_dV()\n", - "system.calculate_bulk_moduli()\n", - "system.calculate_helmholtz_energies_d2V2()\n", - "system.calculate_heat_capacities()" + "system = System(config_objects, ground_state=\"FM\")" ] }, { @@ -246,7 +222,7 @@ "source": [ "# Calculate phase diagrams with specified ground state configuration. \n", "# This also calculates all the properties at constant pressure.\n", - "system.calculate_phase_diagrams(ground_state=\"FM\")" + "system.calculate_phase_diagrams()" ] }, { @@ -254,7 +230,7 @@ "id": "dc6c34ee", "metadata": {}, "source": [ - "The following are example plots vs. volume or temperature. For all of these plots, you can plot for `selected_temperatures` and `selected_volumes`. The available plots are:\n", + "The following are example plots vs. volume or temperature. For all of these plots (apart from vt_phase_diagram), you can plot for `selected_temperatures` and `selected_volumes`. The available plots are:\n", "\n", "**Helmholtz Energy**\n", "- \"helmholtz_energy_vs_volume\"\n", @@ -292,36 +268,6 @@ "system.plot_vt(\"helmholtz_energy_vs_volume\")" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "4f4b7ab4", - "metadata": {}, - "outputs": [], - "source": [ - "system.plot_vt(\"entropy_vs_temperature\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "987276f9", - "metadata": {}, - "outputs": [], - "source": [ - "system.plot_vt(\"heat_capacity_vs_temperature\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2f736a5a", - "metadata": {}, - "outputs": [], - "source": [ - "system.plot_vt(\"bulk_modulus_vs_temperature\")" - ] - }, { "cell_type": "code", "execution_count": null, @@ -337,7 +283,7 @@ "id": "678a6f4f", "metadata": {}, "source": [ - "The following are example plots for properties at constant pressure. The available plots are:\n", + "The following are example plots for properties at constant pressure. You can plot for `selected_temperatures` for helmholtz_energy_pv_vs_volume. The available plots are:\n", "\n", "- \"probability_vs_temperature\"\n", "- \"helmholtz_energy_pv_vs_volume\"\n", @@ -369,30 +315,7 @@ "metadata": {}, "outputs": [], "source": [ - "system.plot_pt(\"probability_vs_temperature\", ground_state=\"FM\", P=0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2a171cbd", - "metadata": {}, - "outputs": [], - "source": [ - "system.plot_pt(\"helmholtz_energy_pv_vs_volume\", P=0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1ca0418c", - "metadata": {}, - "outputs": [], - "source": [ - "# You can also plot for selected temperatures by passing an array of temperatures.\n", - "# The code will automatically find and use the closest available temperatures.\n", - "selected_temperatures = np.array([100, 300, 600, 900])\n", - "system.plot_pt(\"helmholtz_energy_pv_vs_volume\", P=0, selected_temperatures=selected_temperatures)" + "system.plot_pt(\"probability_vs_temperature\", P=0)" ] }, { @@ -434,16 +357,6 @@ "source": [ "system.plot_pt(\"heat_capacity_vs_temperature\", P=0)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bb30e4e6", - "metadata": {}, - "outputs": [], - "source": [ - "system.plot_pt(\"gibbs_energy_vs_temperature\", P=0)" - ] } ], "metadata": { diff --git a/examples/codespace/Fe3Pt_input.pkl b/examples/codespace/Fe3Pt_input.pkl index 62e196c..19772a4 100644 Binary files a/examples/codespace/Fe3Pt_input.pkl and b/examples/codespace/Fe3Pt_input.pkl differ diff --git a/pyzentropy/configuration.py b/pyzentropy/configuration.py index bb84a73..b1c0397 100644 --- a/pyzentropy/configuration.py +++ b/pyzentropy/configuration.py @@ -99,7 +99,6 @@ def __init__( helmholtz_energies: np.ndarray, helmholtz_energies_dV: np.ndarray, helmholtz_energies_d2V2: np.ndarray, - reference_helmholtz_energies: np.ndarray, entropies: np.ndarray = None, heat_capacities: np.ndarray = None, ): @@ -139,39 +138,12 @@ def __init__( if self.entropies is not None: self.calculate_internal_energies() - # Always calculate partition functions using the provided reference Helmholtz energies - self.calculate_partition_functions(reference_helmholtz_energies) - def calculate_internal_energies(self) -> None: """ Calculate internal energies using the formula: E = F + T*S. """ self.internal_energies = self.helmholtz_energies + self.temperatures[:, np.newaxis] * self.entropies - def calculate_partition_functions(self, reference_helmholtz_energies: np.ndarray) -> None: - """ - Calculate the partition function for each state, using a reference Helmholtz energy using the formula: Z = exp(-(F - F_ref) / (k_B * T)). - - Args: - reference_helmholtz_energies (np.ndarray): Reference Helmholtz energies to shift by (shape: [n_temperatures, n_volumes]) - - Raises: - ValueError: If reference_helmholtz_energies does not match the expected shape. - """ - - if reference_helmholtz_energies.shape != (len(self.temperatures), len(self.volumes)): - raise ValueError( - f"reference_helmholtz_energies must have shape {(len(self.temperatures), len(self.volumes))} " - f"Received array with shape {reference_helmholtz_energies.shape}." - ) - - helmholtz_energy_shifted = self.helmholtz_energies - reference_helmholtz_energies - exponent = -helmholtz_energy_shifted / (BOLTZMANN_CONSTANT * self.temperatures[:, np.newaxis]) - with np.errstate(over="ignore"): - partition_functions = np.exp(exponent) - partition_functions[np.isinf(partition_functions)] = np.nan - self.partition_functions = partition_functions - def plot_vt( self, type: str, diff --git a/pyzentropy/system.py b/pyzentropy/system.py index 12077de..2917928 100644 --- a/pyzentropy/system.py +++ b/pyzentropy/system.py @@ -21,6 +21,7 @@ class System: Attributes: configurations (dict[str, Configuration]): Dictionary mapping configuration names to Configuration objects. + ground_state (str): Name of the ground state configuration. number_of_atoms (int): Number of atoms in the system. temperatures (np.ndarray): Array of temperatures considered (shape: [n_temperatures]). volumes (np.ndarray): Array of volumes considered (shape: [n_volumes]). @@ -28,6 +29,7 @@ class System: helmholtz_energies (np.ndarray): Helmholtz energies for the system (shape: [n_temperatures, n_volumes]). helmholtz_energies_dV (np.ndarray): First derivative of Helmholtz energies with respect to volume (shape: [n_temperatures, n_volumes]). helmholtz_energies_d2V2 (np.ndarray): Second derivative of Helmholtz energies with respect to volume (shape: [n_temperatures, n_volumes]). + ground_state_helmholtz_energies (np.ndarray): Helmholtz energies of the ground state configuration (shape: [n_temperatures, n_volumes]). entropies (np.ndarray): Total entropies for the system (shape: [n_temperatures, n_volumes]). configurational_entropies (np.ndarray): Configurational entropies for the system (shape: [n_temperatures, n_volumes]). bulk_moduli (np.ndarray): Bulk moduli for the system (shape: [n_temperatures, n_volumes]). @@ -37,25 +39,28 @@ class System: vt_phase_diagram (dict): Volume-temperature phase diagram data. """ - def __init__(self, configurations: dict[str, Configuration], reference_helmholtz_energies: np.ndarray): + def __init__(self, configurations: dict[str, Configuration], ground_state: str) -> None: """ Initialize the System with a dictionary of Configuration objects. Args: configurations (dict[str, Configuration]): Dictionary mapping configuration names to Configuration objects. - reference_helmholtz_energies (np.ndarray): Reference Helmholtz energies to shift by (shape: [n_temperatures, n_volumes]) - when calculating system Helmholtz energies. + ground_state (str): Name of the ground state configuration. Raises: ValueError: If configurations have inconsistent number of atoms, volumes, or temperatures. """ self.configurations = configurations + self.ground_state = ground_state + + # Get reference values from the ground state configuration + ground_state = next(config for name, config in configurations.items() if config.name == self.ground_state) - # Get reference values from the first configuration - first_config = next(iter(configurations.values())) - ref_num_atoms = first_config.number_of_atoms - ref_volumes = first_config.volumes - ref_temperatures = first_config.temperatures + ref_num_atoms = ground_state.number_of_atoms + ref_volumes = ground_state.volumes + ref_temperatures = ground_state.temperatures + self.ground_state = ground_state.name + self.ground_state_helmholtz_energies = ground_state.helmholtz_energies # Ensure all configurations are consistent in number of atoms, temperatures, and volumes for name, config in configurations.items(): @@ -66,10 +71,6 @@ def __init__(self, configurations: dict[str, Configuration], reference_helmholtz if not np.array_equal(config.temperatures, ref_temperatures): raise ValueError("Temperatures for configurations are not the same.") - # Ensure that reference_helmholtz_energies has the correct shape - if reference_helmholtz_energies.shape != (len(ref_temperatures), len(ref_volumes)): - raise ValueError("reference_helmholtz_energies must have shape (n_temperatures, n_volumes).") - self.number_of_atoms = ref_num_atoms self.temperatures = ref_temperatures self.volumes = ref_volumes @@ -104,28 +105,54 @@ def __init__(self, configurations: dict[str, Configuration], reference_helmholtz # Perform initial calculations self.calculate_partition_functions() self.calculate_probabilities() - self.calculate_helmholtz_energies(reference_helmholtz_energies) + self.calculate_helmholtz_energies(self.ground_state_helmholtz_energies) self.calculate_bulk_moduli() self.calculate_entropies() self.calculate_heat_capacities() def calculate_partition_functions(self) -> None: """ + Calculate the partition function for each configuration, using a reference ground state + Helmholtz energy using the formula: Zk = exp(-(Fk - Fk_ref) / (k_B * T)). Calculate the total partition function for the system by summing over all configurations. - Each configuration's partition function is weighted by its multiplicity. - - Raises: - ValueError: If any configuration is missing its partition function. + Each configuration's partition function is weighted by its multiplicity, Z = Σk mk * Zk. """ + + # Initialize config partition functions + for config in self.configurations.values(): + config.partition_functions = np.zeros((self._n_temps, self._n_vols)) + + # Handle T = 0 K case + # Stack all F(0K) arrays: shape (n_configs, n_vols) + F_0K_all = np.array([config.helmholtz_energies[0, :] for config in self.configurations.values()]) + min_F_0K = np.min(F_0K_all, axis=0) # shape: (n_vols,) + + for config_idx, config in enumerate(self.configurations.values()): + # For each volume, set Z=1 if F is minimal, else 0 + is_min = np.isclose(config.helmholtz_energies[0, :], min_F_0K) + config.partition_functions[0, :] = is_min.astype(float) + + # Handle T > 0 K case + for config in self.configurations.values(): + # Calculate partition function for T > 0 K + with np.errstate(over="ignore", divide="ignore", invalid="ignore"): + for t_idx, temperature in enumerate(self.temperatures): + if temperature == 0: + continue # Skip T = 0 K, already handled + exponent = -( + config.helmholtz_energies[t_idx, :] - self.ground_state_helmholtz_energies[t_idx, :] + ) / (BOLTZMANN_CONSTANT * temperature) + with np.errstate(over="ignore"): + config.partition_functions[t_idx, :] = np.exp(exponent) + # Replace inf values with nan + config.partition_functions[np.isinf(config.partition_functions)] = np.nan + # Initialize partition functions array partition_functions = np.zeros((self._n_temps, self._n_vols)) # Check that the partition functions are calculated for each configuration - with np.errstate(invalid="ignore", over="ignore"): - for config in self.configurations.values(): - if config.partition_functions is None: - raise ValueError(f"partition_functions not set for configuration '{config.name}'.") - partition_functions += config.partition_functions * config.multiplicity + for config in self.configurations.values(): + partition_functions += config.partition_functions * config.multiplicity # Replace inf values with nan partition_functions[np.isinf(partition_functions)] = np.nan @@ -149,12 +176,12 @@ def calculate_probabilities(self) -> None: probabilities = config.multiplicity * config.partition_functions / self.partition_functions config.probabilities = probabilities - def calculate_helmholtz_energies(self, reference_helmholtz_energies: np.ndarray) -> None: + def calculate_helmholtz_energies(self, ground_state_helmholtz_energies: np.ndarray) -> None: """ Calculate the Helmholtz energy for the system at each temperature and volume. Args: - reference_helmholtz_energies (np.ndarray): Reference Helmholtz energies to shift the calculated values. + ground_state_helmholtz_energies (np.ndarray): Reference Helmholtz energies to shift the calculated values. Raises: ValueError: If the system partition function is not calculated. @@ -163,11 +190,26 @@ def calculate_helmholtz_energies(self, reference_helmholtz_energies: np.ndarray) if self.partition_functions is None: raise ValueError("Partition functions not calculated. Call calculate_partition_functions() first.") + # Initialize Helmholtz energies array + self.helmholtz_energies = np.zeros((self._n_temps, self._n_vols)) + # Calculate Helmholtz energies - self.helmholtz_energies = ( - -BOLTZMANN_CONSTANT * self.temperatures[:, np.newaxis] * np.log(self.partition_functions) - + reference_helmholtz_energies - ) + # If T = 0 K + if self.temperatures[0] == 0: + for config in self.configurations.values(): + self.helmholtz_energies[0, :] += config.probabilities[0, :] * config.helmholtz_energies[0, :] + + self.helmholtz_energies[1:, :] = ( + -BOLTZMANN_CONSTANT * self.temperatures[1:, np.newaxis] * np.log(self.partition_functions[1:, :]) + + ground_state_helmholtz_energies[1:, :] + ) + + # For T > 0 K and no T = 0 K + else: + self.helmholtz_energies = ( + -BOLTZMANN_CONSTANT * self.temperatures[:, np.newaxis] * np.log(self.partition_functions) + + ground_state_helmholtz_energies + ) # Calculate derivatives automatically self.calculate_helmholtz_energies_dV() @@ -203,29 +245,38 @@ def calculate_helmholtz_energies_d2V2(self) -> None: Raises: ValueError: If probabilities, dF/dV, or d2F/dV2 are missing for any configuration. """ + # Check that the probabilities, dF/dV, and d2F/dV2 are calculated for each configuration + for config in self.configurations.values(): + if config.probabilities is None: + raise ValueError( + f"Probabilities not set for configuration '{config.name}'. Call calculate_probabilities() first." + ) + if config.helmholtz_energies_dV is None: + raise ValueError(f"helmholtz_energies_dV not set for configuration '{config.name}'.") + if config.helmholtz_energies_d2V2 is None: + raise ValueError(f"helmholtz_energies_d2V2 not set for configuration '{config.name}'.") # Initialize array self.helmholtz_energies_d2V2 = np.zeros((self._n_temps, self._n_vols)) # Loop over temperatures - for i in range(self._n_temps): + for i, temperature in enumerate(self.temperatures): + # At T = 0 K, set to ground state d2F/dV2 + if temperature == 0: + for config in self.configurations.values(): + self.helmholtz_energies_d2V2[i, :] += ( + config.probabilities[i, :] * config.helmholtz_energies_d2V2[i, :] + ) + continue + # For T > 0 K # Initialize averages d2F_dV2_avg = np.zeros(self._n_vols) dF_dV_avg = np.zeros(self._n_vols) dF_dV_sq_avg = np.zeros(self._n_vols) - # Check that the probabilities, dF/dV, and d2F/dV2 are calculated for each configuration + # Accumulate contributions from each configuration for config in self.configurations.values(): - if config.probabilities is None: - raise ValueError( - f"Probabilities not set for configuration '{config.name}'. Call calculate_probabilities() first." - ) - if config.helmholtz_energies_dV is None: - raise ValueError(f"helmholtz_energies_dV not set for configuration '{config.name}'.") - if config.helmholtz_energies_d2V2 is None: - raise ValueError(f"helmholtz_energies_d2V2 not set for configuration '{config.name}'.") - pk = config.probabilities[i] dF_dV = config.helmholtz_energies_dV[i] dF_dV_avg += pk * dF_dV @@ -256,11 +307,7 @@ def calculate_entropies(self) -> None: Raises: ValueError: If probabilities or entropies are missing for any configuration. """ - # Initialize arrays - self.configurational_entropies = np.zeros((self._n_temps, self._n_vols)) - intra_configurational_entropies = np.zeros((self._n_temps, self._n_vols)) - self.entropies = np.zeros((self._n_temps, self._n_vols)) - + # Check that probabilities and entropies are calculated for each configuration for config in self.configurations.values(): if config.probabilities is None: raise ValueError( @@ -269,13 +316,29 @@ def calculate_entropies(self) -> None: if config.entropies is None: raise ValueError(f"Entropies not set for configuration '{config.name}'.") - multiplicity = config.multiplicity - pk = config.probabilities - self.configurational_entropies += ( - -BOLTZMANN_CONSTANT * multiplicity * (pk / multiplicity) * np.log(pk / multiplicity) - ) - intra_configurational_entropies += pk * config.entropies - self.entropies = self.configurational_entropies + intra_configurational_entropies + # Initialize arrays + self.configurational_entropies = np.zeros((self._n_temps, self._n_vols)) + intra_configurational_entropies = np.zeros((self._n_temps, self._n_vols)) + self.entropies = np.zeros((self._n_temps, self._n_vols)) + + for i, temperature in enumerate(self.temperatures): + # For T = 0 K, set entropy to the ground state + if temperature == 0: + for config in self.configurations.values(): + self.entropies[i] += config.probabilities[i, :] * config.entropies[i, :] + continue + + # For T > 0 K + for config in self.configurations.values(): + multiplicity = config.multiplicity + pk = config.probabilities[i] + entropies = config.entropies[i] + + self.configurational_entropies[i] += ( + -BOLTZMANN_CONSTANT * multiplicity * (pk / multiplicity) * np.log(pk / multiplicity) + ) + intra_configurational_entropies[i] += pk * entropies + self.entropies[i] = self.configurational_entropies[i] + intra_configurational_entropies[i] def calculate_heat_capacities(self) -> None: """ @@ -284,12 +347,6 @@ def calculate_heat_capacities(self) -> None: Raises: ValueError: If probabilities, heat capacities, or internal energies are missing for any configuration. """ - # Initialize terms - Cv_avg = np.zeros((self._n_temps, self._n_vols)) - E_sq_avg = np.zeros((self._n_temps, self._n_vols)) - E_avg = np.zeros((self._n_temps, self._n_vols)) - factor = 1 / (BOLTZMANN_CONSTANT * self.temperatures[:, np.newaxis] ** 2) - # Check that the probabilities, heat capacities, and internal energies are calculated for each configuration for config in self.configurations.values(): if config.probabilities is None: @@ -301,13 +358,34 @@ def calculate_heat_capacities(self) -> None: if config.internal_energies is None: raise ValueError(f"Internal energies not set for configuration '{config.name}'.") - # Accumulate contributions from each configuration - Cv_avg += config.probabilities * config.heat_capacities - E_sq_avg += config.probabilities * config.internal_energies**2 - E_avg += config.probabilities * config.internal_energies + # Initialize terms + self.heat_capacities = np.zeros((self._n_temps, self._n_vols)) + Cv_avg = np.zeros((self._n_temps, self._n_vols)) + E_sq_avg = np.zeros((self._n_temps, self._n_vols)) + E_avg = np.zeros((self._n_temps, self._n_vols)) + factor = np.zeros((self._n_temps, 1)) + if self.temperatures[0] == 0: + factor[0] = 0 + factor[1:] = 1 / (BOLTZMANN_CONSTANT * self.temperatures[1:, np.newaxis] ** 2) + else: + factor = 1 / (BOLTZMANN_CONSTANT * self.temperatures[:, np.newaxis] ** 2) - # Final heat capacity calculation - self.heat_capacities = Cv_avg + factor * (E_sq_avg - E_avg**2) + for i, temperature in enumerate(self.temperatures): + # For T = 0 K, set entropy to the ground state + if temperature == 0: + for config in self.configurations.values(): + self.heat_capacities[i, :] += config.probabilities[i, :] * config.heat_capacities[i, :] + continue + + # For T > 0 K + for config in self.configurations.values(): + # Accumulate contributions from each configuration + Cv_avg[i] += config.probabilities[i] * config.heat_capacities[i] + E_sq_avg[i] += config.probabilities[i] * config.internal_energies[i] ** 2 + E_avg[i] += config.probabilities[i] * config.internal_energies[i] + + # Final heat capacity calculation + self.heat_capacities[i] = Cv_avg[i] + factor[i] * (E_sq_avg[i] - E_avg[i] ** 2) def calculate_pressure_properties(self, P: float) -> None: """ @@ -510,27 +588,20 @@ def calculate_pressure_properties(self, P: float) -> None: dS_dT = np.diff(self.pt_properties[f"{P:.2f}_GPa"]["S0"]) / np.diff(self.temperatures) self.pt_properties[f"{P:.2f}_GPa"]["Cp"] = self.temperatures[:-1] * dS_dT - def calculate_phase_diagrams( - self, ground_state: str, dP: float = 0.2, volume_step_size: float = 1e-4, atol: float = 1e-6 - ) -> None: + def calculate_phase_diagrams(self, dP: float = 0.2, volume_step_size: float = 1e-4, atol: float = 1e-6) -> None: """ Calculate pressure-temperature and volume-temperature phase diagrams for the system. Args: - ground_state (str): Name of the ground state configuration. dP (float): Pressure increment in GPa. Default is 0.2 GPa. volume_step_size (float): Step size for volume when searching for the miscibility gap. Default is 1e-4. atol (float): Absolute tolerance for convergence in the common tangent search. Default is 1e-6. Raises: - ValueError: If any of the Helmholtz energies or their derivatives are not calculated, + ValueError: If any of the Helmholtz energies or their first derivatives are not calculated, or if probabilities for any configuration are not calculated. """ - - if self.helmholtz_energies_d2V2 is None: - raise ValueError( - "Helmholtz energies second derivative not calculated. Call calculate_helmholtz_energies_d2V2() first." - ) + # Raise errors if required attributes are missing if self.helmholtz_energies_dV is None: raise ValueError( "Helmholtz energies first derivative not calculated. Call calculate_helmholtz_energies_dV() first." @@ -558,7 +629,7 @@ def calculate_phase_diagrams( while True: try: self.calculate_pressure_properties(P) - gs_probabilities = self.configurations[ground_state].probabilities_at_P[ + gs_probabilities = self.configurations[self.ground_state].probabilities_at_P[ f"{P:.2f}_GPa" ] # probabilities vs T V0 = self.pt_properties[f"{P:.2f}_GPa"]["V0"] # V0 vs T @@ -612,45 +683,48 @@ def calculate_phase_diagrams( # Only consider 0 GPa for miscibility gap for T-V diagram # Current method can only handle one common tangent + helmholtz_energies = self.helmholtz_energies[index, :] helmholtz_energies_dV = self.helmholtz_energies_dV[index, :] - helmholtz_energies_d2V2 = self.helmholtz_energies_d2V2[index, :] - # Interpolators for dF/dV and d^2F/dV^2 + # Interpolators for F and dF/dV + F_interpolator = PchipInterpolator(self.volumes, helmholtz_energies) dV_interpolator = PchipInterpolator(self.volumes, helmholtz_energies_dV) - d2V2_interpolator = PchipInterpolator(self.volumes, helmholtz_energies_d2V2) - # Find sign changes in d^2F/dV^2 (roots) - roots = [] - for i in range(len(self.volumes) - 1): - if ( - helmholtz_energies_d2V2[i] * helmholtz_energies_d2V2[i + 1] - ) < 0: # True if there is a change of sign - res = root_scalar( - d2V2_interpolator, - bracket=(self.volumes[i], self.volumes[i + 1]), - method="brentq", - xtol=1e-10, - rtol=1e-12, - ) - if res.converged: - roots.append(res.root) + # Find when dF/dV changes sign between consecutive volumes + dy = np.diff(helmholtz_energies_dV) + sign_change = np.diff(np.sign(dy)) + idx = np.where(sign_change != 0)[0] + 1 # +1 to correct index after diff # Move left from the first root and right from the second root until dF/dV values are approximately equal (miscibility gap volumes) - if roots: - left_root = roots[0] - right_root = roots[-1] + if len(idx) > 0: + left_root = idx[0] + right_root = idx[-1] - left_volume = left_root - volume_step_size - right_volume = right_root + volume_step_size + left_volume = self.volumes[left_root] + right_volume = self.volumes[right_root] left_helmholtz_energy_dV = dV_interpolator(left_volume) right_helmholtz_energy_dV = dV_interpolator(right_volume) - while ~np.isclose(left_helmholtz_energy_dV, right_helmholtz_energy_dV, atol=atol): + while not (np.isclose(left_helmholtz_energy_dV, right_helmholtz_energy_dV, atol=atol)): left_volume -= volume_step_size right_volume += volume_step_size left_helmholtz_energy_dV = dV_interpolator(left_volume) right_helmholtz_energy_dV = dV_interpolator(right_volume) + left_helmholtz_energy = F_interpolator(left_volume) + right_helmholtz_energy = F_interpolator(right_volume) + secant_slope = (right_helmholtz_energy - left_helmholtz_energy) / (right_volume - left_volume) + + # Check that all three values are within a tolerance of 0.001 + if not ( + np.isclose(secant_slope, left_helmholtz_energy_dV, atol=0.001) + and np.isclose(secant_slope, right_helmholtz_energy_dV, atol=0.001) + and np.isclose(left_helmholtz_energy_dV, right_helmholtz_energy_dV, atol=0.001) + ): + print("Warning: Secant slope and derivatives are not all within a tolerance of 0.001.") + print(secant_slope, left_helmholtz_energy_dV, right_helmholtz_energy_dV) + break + self.pt_phase_diagram["first_order"]["P"] = np.append( self.pt_phase_diagram["first_order"]["P"], np.round(-left_helmholtz_energy_dV * EV_PER_CUBIC_ANGSTROM_TO_GPA, 2), @@ -929,7 +1003,6 @@ def plot_pt( type: str, P: float = 0.00, selected_temperatures: np.ndarray = None, - ground_state: str = None, width: int = 650, height: int = 600, ): @@ -940,7 +1013,6 @@ def plot_pt( type (str): The type of plot to generate (e.g., "helmholtz_energy_pv_vs_volume"). P (float): Pressure in GPa for the plot. Default is 0.00 GPa. selected_temperatures (np.ndarray, optional): Temperatures to highlight in the helmholtz_energy_pv_vs_volume plot. - ground_state (str, optional): Name of the ground state configuration for probability plots. width (int, optional): Width of the plot in pixels. height (int, optional): Height of the plot in pixels. @@ -1072,63 +1144,34 @@ def plot_pt( if type == "probability_vs_temperature": traces = [] - # Add ground state first if present - if ground_state is not None and ground_state in y_data: - traces.append( - go.Scatter( - x=x_data, - y=y_data[ground_state], - mode="lines", - name=ground_state, - ) - ) - gs_probs = np.array(y_data[ground_state]) - # Interpolate to find temperature where probability crosses 0.5 - try: - interp = PchipInterpolator(x_data, gs_probs) - roots = [] - for i in range(len(x_data) - 1): - # Check if the probability crosses 0.5 between x_data[i] and x_data[i+1] - if (gs_probs[i] - 0.5) * (gs_probs[i + 1] - 0.5) < 0: - res = root_scalar( - lambda T: interp(T) - 0.5, - bracket=(x_data[i], x_data[i + 1]), - method="brentq", - xtol=1e-10, - rtol=1e-12, - ) - if res.converged: - roots.append(res.root) - except Exception: - pass - # Sum all probabilities except the ground state - excited_probs = np.zeros_like(x_data, dtype=float) - for name, probs in y_data.items(): - if name != ground_state: - excited_probs += probs - traces.append( - go.Scatter( - x=x_data, - y=excited_probs, - mode="lines", - name="Sum (Non-GS)", - line=dict(color="black"), - ) + + # Add ground state first + traces.append( + go.Scatter( + x=x_data, + y=y_data[self.ground_state], + mode="lines", + name=self.ground_state, ) - # Add the rest (non-ground states) - for name, probs in y_data.items(): - if name != ground_state: - traces.append( - go.Scatter( - x=x_data, - y=probs, - mode="lines", - name=name, - ) - ) + ) - else: - for name, probs in y_data.items(): + # Sum all probabilities except the ground state + excited_probs = np.zeros_like(x_data, dtype=float) + for name, probs in y_data.items(): + if name != self.ground_state: + excited_probs += probs + traces.append( + go.Scatter( + x=x_data, + y=excited_probs, + mode="lines", + name="Sum (Non-GS)", + line=dict(color="black"), + ) + ) + # Add the rest (non-ground states) + for name, probs in y_data.items(): + if name != self.ground_state: traces.append( go.Scatter( x=x_data, @@ -1144,7 +1187,7 @@ def plot_pt( y_label = "Probability" fig.update_layout( title=dict( - text=f"P = {P:.2f} GPa" if P is not None else "P = None", + text=f"P = {P:.2f} GPa", font=dict(size=22, color="rgb(0,0,0)"), ), yaxis=dict(range=[0, None]), @@ -1226,7 +1269,7 @@ def plot_pt( ) ) if fixed_by == "pressure" or type == "helmholtz_energy_pv_vs_volume": - title_text = f"P = {P:.2f} GPa" if P is not None else "P = None" + title_text = f"P = {P:.2f} GPa" fig.update_layout(title=dict(text=title_text, font=dict(size=22, color="rgb(0,0,0)"))) x_label = "Temperature (K)" if "temperature" in type else "Volume (ų)" diff --git a/tests/test_configuration.py b/tests/test_configuration.py index 5ac9486..375b438 100644 --- a/tests/test_configuration.py +++ b/tests/test_configuration.py @@ -35,7 +35,6 @@ def test_initialization(): helmholtz_energies=arr_wrong, # wrong shape here helmholtz_energies_dV=arr, helmholtz_energies_d2V2=arr, - reference_helmholtz_energies=arr, entropies=arr, heat_capacities=arr, ) @@ -45,7 +44,6 @@ def test_initialization(): def test_configuration_param(config_key): # Test Configuration internal energy and partition function calculations. config = config_data[config_key] - reference_helmholtz_energies = config_data["FM"].helmholtz_energies instance = Configuration( name=config.name, multiplicity=config.multiplicity, @@ -55,20 +53,10 @@ def test_configuration_param(config_key): helmholtz_energies=config.helmholtz_energies, helmholtz_energies_dV=config.helmholtz_energies_dV, helmholtz_energies_d2V2=config.helmholtz_energies_d2V2, - reference_helmholtz_energies=reference_helmholtz_energies, entropies=config.entropies, heat_capacities=config.heat_capacities, ) assert np.allclose(instance.internal_energies, config.internal_energies) - assert np.allclose(instance.partition_functions, config.partition_functions, equal_nan=True) - - # Test invalid reference_helmholtz_energies shape raises ValueError - wrong_shape_ref = np.zeros((len(instance.temperatures) - 1, len(instance.volumes))) - with pytest.raises( - ValueError, - match=rf"reference_helmholtz_energies must have shape \({len(instance.temperatures)}, {len(instance.volumes)}\) Received array with shape \({wrong_shape_ref.shape[0]}, {wrong_shape_ref.shape[1]}\)\.", - ): - instance.calculate_partition_functions(wrong_shape_ref) def test_plot_vt(): diff --git a/tests/test_data/Fe3Pt_system.pkl b/tests/test_data/Fe3Pt_system.pkl index dc46168..1efd3ee 100644 Binary files a/tests/test_data/Fe3Pt_system.pkl and b/tests/test_data/Fe3Pt_system.pkl differ diff --git a/tests/test_data/Fe3Pt_three_configs.pkl b/tests/test_data/Fe3Pt_three_configs.pkl index 9264062..05135bd 100644 Binary files a/tests/test_data/Fe3Pt_three_configs.pkl and b/tests/test_data/Fe3Pt_three_configs.pkl differ diff --git a/tests/test_data/expected_results.ipynb b/tests/test_data/expected_results.ipynb index 8380084..18eb694 100644 --- a/tests/test_data/expected_results.ipynb +++ b/tests/test_data/expected_results.ipynb @@ -199,28 +199,24 @@ "source": [ "# Initialize dictionary to hold PyZentropy Configuration objects\n", "new_config_objects = {}\n", - "reference_helmholtz_energies = config_objects[\"config_0\"].qha.methods['debye']['helmholtz_energy']['values'][1:, :]\n", "\n", "for i, name in enumerate(config_names):\n", " config = config_objects[name]\n", "\n", - " # Remove 0 K results from each configuration.\n", - " # TODO: we want to include this in the near future\n", " # Create a Configuration object for each config and store in the dictionary\n", " new_config_objects[name] = PyZentropyConfiguration(\n", " name=name,\n", " multiplicity=config.multiplicity,\n", " number_of_atoms=config.qha.number_of_atoms,\n", " volumes=config.qha.volumes,\n", - " temperatures=config.qha.temperatures[1:],\n", + " temperatures=config.qha.temperatures,\n", " # Extract Helmholtz energies and their derivatives from DFTTK results\n", - " helmholtz_energies=config.qha.methods['debye']['helmholtz_energy']['values'][1:, :],\n", - " helmholtz_energies_dV=config.qha.methods['debye']['helmholtz_energy']['derivatives'][1:, :],\n", - " helmholtz_energies_d2V2=config.qha.methods['debye']['helmholtz_energy']['derivatives2'][1:, :],\n", - " reference_helmholtz_energies=reference_helmholtz_energies,\n", + " helmholtz_energies=config.qha.methods['debye']['helmholtz_energy']['values'][:, :],\n", + " helmholtz_energies_dV=config.qha.methods['debye']['helmholtz_energy']['derivatives'][:, :],\n", + " helmholtz_energies_d2V2=config.qha.methods['debye']['helmholtz_energy']['derivatives2'][:, :],\n", " # Extract entropy and heat capacity data\n", - " entropies=config.qha.methods['debye']['entropy']['values'][1:, :],\n", - " heat_capacities=config.qha.methods['debye']['heat_capacity']['values'][1:, :]\n", + " entropies=config.qha.methods['debye']['entropy']['values'][:, :],\n", + " heat_capacities=config.qha.methods['debye']['heat_capacity']['values'][:, :]\n", " )" ] }, @@ -255,6 +251,38 @@ " pickle.dump(new_config_objects, f)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "95b56472", + "metadata": {}, + "outputs": [], + "source": [ + "# Save Fe3Pt input for tutorial\n", + "Fe3Pt_input = {\n", + " \"FM\": {},\n", + " \"SF22\": {},\n", + " \"SF28\": {}\n", + "}\n", + "\n", + "for config in new_config_objects.values():\n", + " Fe3Pt_input[config.name] = {\n", + " \"multiplicity\": config.multiplicity,\n", + " \"number_of_atoms\": config.number_of_atoms,\n", + " \"volumes\": config.volumes,\n", + " \"temperatures\": config.temperatures,\n", + " \"helmholtz_energies\": config.helmholtz_energies,\n", + " \"helmholtz_energies_dV\": config.helmholtz_energies_dV,\n", + " \"helmholtz_energies_d2V2\": config.helmholtz_energies_d2V2,\n", + " \"entropies\": config.entropies,\n", + " \"heat_capacities\": config.heat_capacities\n", + " }\n", + "\n", + "# Save the input dictionary to a pickle file\n", + "with open(\"Fe3Pt_input.pkl\", \"wb\") as f:\n", + " pickle.dump(Fe3Pt_input, f)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -329,7 +357,7 @@ "source": [ "# Create a System object from the configuration objects\n", "# Thermodynamic calculations for the system will be performed automatically during initialization\n", - "system = System(new_config_objects, reference_helmholtz_energies)" + "system = System(new_config_objects, ground_state=\"FM\")" ] }, { @@ -340,7 +368,7 @@ "outputs": [], "source": [ "# Calculate phase diagrams\n", - "system.calculate_phase_diagrams(ground_state=\"FM\", dP=0.5)" + "system.calculate_phase_diagrams(dP=0.5)" ] }, { @@ -370,10 +398,10 @@ "# heat_capacity_vs_volume, heat_capacity_vs_temperature\n", "# bulk_modulus_vs_volume, bulk_modulus_vs_temperature\n", "# vt_phase_diagram\n", - "selected_temperatures = np.array([300, 400, 500, 600, 700, 800, 900, 1000])\n", + "selected_temperatures = np.array([0, 10, 20, 100, 300, 400, 500, 600, 700, 800, 900, 1000])\n", "selected_volumes = np.array([100, 150, 200])\n", - "system.plot_vt(\"entropy_vs_volume\")\n", - "#system.plot_vt(\"helmholtz_energy_vs_temperature\", selected_temperatures=selected_temperatures, selected_volumes=selected_volumes)" + "system.plot_vt(\"helmholtz_energy_dV_vs_volume\", selected_temperatures=selected_temperatures, selected_volumes=selected_volumes)\n", + "#system.plot_vt(\"bulk_modulus_vs_temperature\", selected_temperatures=selected_temperatures, selected_volumes=selected_volumes)" ] }, { @@ -385,7 +413,7 @@ "source": [ "# Options - helmholtz_energy_pt_vs_volume\n", "selected_temperatures = np.array([300, 400, 500, 600, 700, 800, 900, 1000])\n", - "system.plot_pt(\"volume_vs_temperature\", P=0, selected_temperatures=selected_temperatures, ground_state=\"FM\")\n", + "system.plot_pt(\"helmholtz_energy_pv_vs_volume\", P=0, selected_temperatures=selected_temperatures)\n", "\n", "# volume_vs_temperature, CTE_vs_temperature, LCTE_vs_temperature\n", "# entropy_vs_temperature, configurational_entropy_vs_temperature\n", diff --git a/tests/test_system.py b/tests/test_system.py index f4cf66c..5985cbf 100644 --- a/tests/test_system.py +++ b/tests/test_system.py @@ -18,7 +18,6 @@ test_data_path = os.path.join(os.path.dirname(__file__), "test_data", "Fe3Pt_three_configs.pkl") with open(test_data_path, "rb") as f: config_data = pickle.load(f) -reference_helmholtz_energies = config_data["FM"].helmholtz_energies expected_results_path = os.path.join(os.path.dirname(__file__), "test_data", "Fe3Pt_system.pkl") with open(expected_results_path, "rb") as f: @@ -41,15 +40,14 @@ def make_config(name, number_of_atoms=1, volumes=None, temperatures=None): helmholtz_energies=arr, helmholtz_energies_dV=arr, helmholtz_energies_d2V2=arr, - reference_helmholtz_energies=arr, entropies=arr, heat_capacities=arr, ) -def make_system(configs, reference_helmholtz_energies): +def make_system(configs): """Create a System and run all calculations up to heat capacities.""" - system = System(configs, reference_helmholtz_energies) + system = System(configs, ground_state="FM") return system @@ -58,46 +56,42 @@ def test_system_inconsistent_number_of_atoms(): """Test error if configs have different number of atoms.""" config1 = make_config("A", number_of_atoms=1, volumes=np.arange(2), temperatures=np.arange(3)) config2 = make_config("B", number_of_atoms=2, volumes=np.arange(2), temperatures=np.arange(3)) - reference_helmholtz_energies = np.zeros((3, 2)) with pytest.raises(ValueError, match="Number of atoms for configurations are not the same"): - System({"A": config1, "B": config2}, reference_helmholtz_energies=reference_helmholtz_energies) + System({"A": config1, "B": config2}, ground_state="A") def test_system_inconsistent_volumes(): """Test error if configs have different volumes.""" config1 = make_config("A", volumes=np.arange(2), temperatures=np.arange(3)) config2 = make_config("B", volumes=np.arange(3), temperatures=np.arange(3)) - reference_helmholtz_energies = np.zeros((3, 2)) with pytest.raises(ValueError, match="Volumes for configurations are not the same"): - System({"A": config1, "B": config2}, reference_helmholtz_energies=reference_helmholtz_energies) + System({"A": config1, "B": config2}, ground_state="A") def test_system_inconsistent_temperatures(): """Test error if configs have different temperatures.""" config1 = make_config("A", volumes=np.arange(2), temperatures=np.arange(3)) config2 = make_config("B", volumes=np.arange(2), temperatures=np.arange(4)) - reference_helmholtz_energies = np.zeros((3, 2)) with pytest.raises(ValueError, match="Temperatures for configurations are not the same"): - System({"A": config1, "B": config2}, reference_helmholtz_energies=reference_helmholtz_energies) + System({"A": config1, "B": config2}, ground_state="A") # Calculation Tests def test_partition_functions(): - """Test partition functions and error on missing config partition functions.""" + """Test partition functions.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) - assert np.allclose(system.partition_functions, expected_results.partition_functions, equal_nan=True) - # Test error if config partition_functions is None + system = make_system(local_config_data) for config in system.configurations.values(): - config.partition_functions = None - with pytest.raises(ValueError, match=f"partition_functions not set for configuration 'FM'"): - system.calculate_partition_functions() + assert np.allclose( + config.partition_functions, expected_results.configurations[config.name].partition_functions, equal_nan=True + ) + assert np.allclose(system.partition_functions, expected_results.partition_functions, equal_nan=True) def test_probabilities(): """Test probabilities and error on missing system partition functions.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) for config in system.configurations.values(): assert np.allclose( config.probabilities, expected_results.configurations[config.name].probabilities, equal_nan=True @@ -119,19 +113,20 @@ def test_probabilities(): def test_helmholtz_energies(): """Test Helmholtz energies and error on missing system partition functions.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.helmholtz_energies, expected_results.helmholtz_energies, equal_nan=True) system.partition_functions = None + ground_state_helmholtz_energies = expected_results.configurations["FM"].helmholtz_energies with pytest.raises( ValueError, match=re.escape("Partition functions not calculated. Call calculate_partition_functions() first.") ): - system.calculate_helmholtz_energies(reference_helmholtz_energies) + system.calculate_helmholtz_energies(ground_state_helmholtz_energies) def test_helmholtz_energies_dV(): """Test dF/dV and error on missing config dF/dV or probabilities.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.helmholtz_energies_dV, expected_results.helmholtz_energies_dV, equal_nan=True) # Error if config dF/dV is None for config in system.configurations.values(): @@ -151,7 +146,7 @@ def test_helmholtz_energies_dV(): def test_helmholtz_energies_d2V2(): """Test d2F/dV2 and errors for missing config data.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.helmholtz_energies_d2V2, expected_results.helmholtz_energies_d2V2, equal_nan=True) # Error if config d2F/dV2 is None for config in system.configurations.values(): @@ -176,7 +171,7 @@ def test_helmholtz_energies_d2V2(): def test_bulk_moduli(): """Test bulk moduli and error on missing system d^2F/dV2.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.bulk_moduli, expected_results.bulk_moduli, equal_nan=True) # Error if config d2F/dV2 is None system.helmholtz_energies_d2V2 = None @@ -190,7 +185,7 @@ def test_bulk_moduli(): def test_entropies(): """Test entropies and error on missing config data.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.entropies, expected_results.entropies, equal_nan=True) assert np.allclose(system.configurational_entropies, expected_results.configurational_entropies, equal_nan=True) # Error if config entropies is None @@ -211,7 +206,7 @@ def test_entropies(): def test_heat_capacities(): """Test heat capacities and error on missing config data.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) assert np.allclose(system.heat_capacities, expected_results.heat_capacities, equal_nan=True) # Error if config internal_energies is None for config in system.configurations.values(): @@ -236,7 +231,7 @@ def test_heat_capacities(): def test_calculate_pressure_properties(): """Test pressure properties and edge cases (all entropies None, bulk_moduli None, helmholtz_energies None, helmholtz_energies_dV None).""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.calculate_pressure_properties(P=0) # Test values against expected results assert np.allclose( @@ -271,7 +266,7 @@ def test_calculate_pressure_properties(): ) # Test with all entropies set to None - system2 = make_system(local_config_data, reference_helmholtz_energies) + system2 = make_system(local_config_data) system2.configurational_entropies = None system2.entropies = None system2.calculate_pressure_properties(P=0) @@ -301,7 +296,7 @@ def test_calculate_pressure_properties(): ) # Test with helmholtz_energies set to None - system3 = make_system(local_config_data, reference_helmholtz_energies) + system3 = make_system(local_config_data) system3.helmholtz_energies = None with pytest.raises( ValueError, match=re.escape("Helmholtz energies not calculated. Call calculate_helmholtz_energies() first.") @@ -309,7 +304,7 @@ def test_calculate_pressure_properties(): system3.calculate_pressure_properties(P=0) # Test with helmholtz_energies_dV set to None - system4 = make_system(local_config_data, reference_helmholtz_energies) + system4 = make_system(local_config_data) system4.helmholtz_energies_dV = None with pytest.raises( ValueError, @@ -318,7 +313,7 @@ def test_calculate_pressure_properties(): system4.calculate_pressure_properties(P=0) # Test with config probabilities set to None - system5 = make_system(local_config_data, reference_helmholtz_energies) + system5 = make_system(local_config_data) for config in system5.configurations.values(): config.probabilities = None with pytest.raises( @@ -330,8 +325,8 @@ def test_calculate_pressure_properties(): def test_calculate_phase_diagrams(): local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) - system.calculate_phase_diagrams(ground_state="FM", dP=0.5) + system = make_system(local_config_data) + system.calculate_phase_diagrams(dP=0.5) # Test values against expected results assert np.allclose( system.pt_phase_diagram["first_order"]["P"], expected_results.pt_phase_diagram["first_order"]["P"] @@ -362,18 +357,8 @@ def test_calculate_phase_diagrams(): system.vt_phase_diagram["second_order"]["T"], expected_results.vt_phase_diagram["second_order"]["T"] ) - # Test with helmholtz_energies_d2V2 to None - system.helmholtz_energies_d2V2 = None - with pytest.raises( - ValueError, - match=re.escape( - "Helmholtz energies second derivative not calculated. Call calculate_helmholtz_energies_d2V2() first." - ), - ): - system.calculate_phase_diagrams(ground_state="FM") - # Test with helmholtz_energies_dV set to None - system2 = make_system(local_config_data, reference_helmholtz_energies) + system2 = make_system(local_config_data) system2.helmholtz_energies_dV = None with pytest.raises( ValueError, @@ -381,25 +366,25 @@ def test_calculate_phase_diagrams(): "Helmholtz energies first derivative not calculated. Call calculate_helmholtz_energies_dV() first." ), ): - system2.calculate_phase_diagrams(ground_state="FM") + system2.calculate_phase_diagrams() # Test with helmholtz_energies set to None - system3 = make_system(local_config_data, reference_helmholtz_energies) + system3 = make_system(local_config_data) system3.helmholtz_energies = None with pytest.raises( ValueError, match=re.escape("Helmholtz energies not calculated. Call calculate_helmholtz_energies() first.") ): - system3.calculate_phase_diagrams(ground_state="FM") + system3.calculate_phase_diagrams() # Test with config probabilities set to None - system4 = make_system(local_config_data, reference_helmholtz_energies) + system4 = make_system(local_config_data) for config in system4.configurations.values(): config.probabilities = None with pytest.raises( ValueError, match=re.escape("Probabilities not set for configuration 'FM'. Call calculate_probabilities() first."), ): - system4.calculate_phase_diagrams(ground_state="FM") + system4.calculate_phase_diagrams() # Plotting Tests for plot_vt @@ -426,16 +411,16 @@ def test_calculate_phase_diagrams(): def test_plot_vt_smoke(plot_type): """Test that plot_vt runs without error for all supported plot types.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) if plot_type == "vt_phase_diagram": - system.calculate_phase_diagrams(ground_state="FM") + system.calculate_phase_diagrams() system.plot_vt(plot_type) def test_plot_vt_invalid_type(): """Test that an invalid plot type raises ValueError.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) with pytest.raises(ValueError, match="Invalid plot type"): system.plot_vt("not_a_real_plot_type") @@ -463,7 +448,7 @@ def test_plot_vt_invalid_type(): def test_plot_vt_missing_data(plot_type, attr): """Test that missing required data for each plot type raises ValueError.""" local_config_data = copy.deepcopy(config_data) - system = System(local_config_data, reference_helmholtz_energies=reference_helmholtz_energies) + system = System(local_config_data, ground_state="FM") if plot_type != "vt_phase_diagram": setattr(system, attr, None) with pytest.raises( @@ -487,7 +472,7 @@ def test_plot_vt_missing_data(plot_type, attr): def test_plot_vt_selected_temperatures(plot_type): """Test that plot_vt works with selected_temperatures argument for relevant plot types.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.plot_vt(plot_type, selected_temperatures=np.array([300, 600, 900])) @@ -506,7 +491,7 @@ def test_plot_vt_selected_temperatures(plot_type): def test_plot_vt_selected_volumes(plot_type): """Test that plot_vt works with selected_volumes argument.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.plot_vt(plot_type, selected_volumes=np.array([100, 150, 200])) @@ -530,17 +515,17 @@ def test_plot_vt_selected_volumes(plot_type): def test_plot_pt_smoke(plot_type): """Test that plot_pt runs without error for all supported plot types.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.calculate_pressure_properties(P=0) if plot_type == "pt_phase_diagram": - system.calculate_phase_diagrams(ground_state="FM") + system.calculate_phase_diagrams() system.plot_pt(plot_type) def test_plot_pt_properties_at_P_missing(): """Test that plot_pt raises ValueError if pressure properties are not calculated.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) with pytest.raises( ValueError, match="Properties at 0.00 GPa not calculated. Run calculate_pressure_properties() first." ): @@ -550,7 +535,7 @@ def test_plot_pt_properties_at_P_missing(): def test_plot_pt_invalid_type(): """Test that an invalid plot type raises ValueError.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.calculate_pressure_properties(P=0) with pytest.raises(ValueError, match="Invalid plot type"): system.plot_pt("not_a_real_plot_type") @@ -563,7 +548,7 @@ def test_plot_pt_invalid_type(): def test_plot_pt_missing_data(plot_type): """Test that missing required data for each plot type raises ValueError.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.calculate_pressure_properties(P=0) expected_msg = re.escape(f"{plot_type} data not calculated. Run the appropriate calculation method first.") with pytest.raises(ValueError, match=expected_msg): @@ -577,14 +562,6 @@ def test_plot_pt_missing_data(plot_type): def test_plot_pt_selected_temperatures(plot_type): """Test that plot_pt works with selected_temperatures argument for relevant plot types.""" local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) + system = make_system(local_config_data) system.calculate_pressure_properties(P=0) system.plot_pt(plot_type, selected_temperatures=np.array([300, 600, 900])) - - -def test_plot_pt_probabilities_ground_state(): - """Test that plot_pt works with ground_state argument.""" - local_config_data = copy.deepcopy(config_data) - system = make_system(local_config_data, reference_helmholtz_energies) - system.calculate_pressure_properties(P=0) - system.plot_pt("probability_vs_temperature", ground_state="FM")