10736148419

Committed 06 Sep 2024 09:40AM UTC coverage: 70.934%. First build

Build # 10736148419

Build Type

Pull #1555

github

Committed by

web-flow

Commit Message

Merge 8fa6b31f7 into 208d0516f

Pull Request Pull Request #1555: Merge main

Run Details

74 of 85 new or added lines in 9 files covered. (87.06%)

10689 of 15069 relevant lines covered (70.93%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

25.69

/pyiron_atomistics/atomistics/master/quasi.py

# coding: utf-8
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
# Distributed under the terms of "New BSD License", see the LICENSE file.

import matplotlib
import numpy as np

from pyiron_atomistics.atomistics.master.murnaghan import MurnaghanJobGenerator
from pyiron_atomistics.atomistics.master.parallel import AtomisticParallelMaster

__author__ = "Jan Janssen"
__copyright__ = (
    "Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH - "
    "Computational Materials Design (CM) Department"
)
__version__ = "0.0.1"
__maintainer__ = "Jan Janssen"
__email__ = "janssen@mpie.de"
__status__ = "development"
__date__ = "Oct 29, 2020"


def calc_v0_from_fit_funct(fit_funct, x, save_range=0.0, return_ind=False):
    fit_funct_der = fit_funct.deriv().r
    fit_funct_der_r = fit_funct_der[fit_funct_der.imag == 0].real
    fit_funct_der_val = fit_funct.deriv(2)(fit_funct_der_r)
    select = (
        (fit_funct_der_val > 0)
        & (fit_funct_der_r > np.min(x) * (1 - save_range))
        & (fit_funct_der_r < np.max(x) * (1 + save_range))
    )
    v0_lst = fit_funct_der_r[select]
    if len(v0_lst) == 1:
        if return_ind:
            return v0_lst[0], (np.abs(x - v0_lst[0])).argmin()
        else:
            return v0_lst[0]
    else:
        select = fit_funct_der_val > 0
        v0_lst = fit_funct_der_r[select]
        if len(v0_lst) == 1:
            if return_ind:
                return v0_lst[0], (np.abs(x - v0_lst[0])).argmin()
            else:
                return v0_lst[0]
        else:
            if return_ind:
                return None, None
            else:
                return None


class QuasiHarmonicJob(AtomisticParallelMaster):
    """
    Obtain finite temperature properties in the framework of quasi harmonic approximation. For the
    theoretical understanding take a look at the Wikipedia page:
    https://en.wikipedia.org/wiki/Quasi-harmonic_approximation

    Example:

    >>> pr = Project('my_project')
    >>> lmp = pr.create_job('Lammps', 'lmp')
    >>> lmp.structure = structure_of_your_choice
    >>> phono = lmp.create_job('PhonopyJob', 'phono')
    >>> qha = phono.create_job('QuasiHarmonicJob', 'qha')
    >>> qha.run()

    The final results can be obtained through `qha.optimise_volume()`.

    The temperature range defined in the input can be modified afterwards. For this, follow these
    lines:

    >>> qha.input['temperature_end'] = temperature_end
    >>> qha.input['temperature_steps'] = temperature_steps
    >>> qha.input['temperature_start'] = temperature_start
    >>> qha.collect_output()
    """

    def __init__(self, project, job_name="murnaghan"):
        """

        Args:
            project:
            job_name:
        """
        super(QuasiHarmonicJob, self).__init__(project, job_name)

        self.__version__ = "0.0.1"

        # define default input
        self.input["num_points"] = (11, "number of sample points")
        self.input["vol_range"] = (
            0.1,
            "relative volume variation around volume defined by ref_ham",
        )
        self.input["temperature_start"] = 0
        self.input["temperature_end"] = 500
        self.input["temperature_steps"] = 10
        self.input["polynomial_degree"] = 3
        self.input["strains"] = (
            None,
            "List of strains that should be calculated.  If given vol_range and num_points take no effect.",
        )
        self.input["axes"] = (
            ["x", "y", "z"],
            "Axes along which the strain will be applied",
        )
        self._job_generator = MurnaghanJobGenerator(self)

    def collect_output(self):
        free_energy_lst, entropy_lst, cv_lst, volume_lst = [], [], [], []
        for job_id in self.child_ids:
            job = self.project_hdf5.load(job_id)
            # the underlying phonopy job might create a supercell; if its
            # reference job is interactive this is reflected in its
            # job.structure and the output quantities, so we need to rescale
            # all the output here; if the structure did not change the
            # conversion_factor will simply be 1.
            conversion_factor = len(self.structure) / len(job.structure)
            thermal_properties = job.get_thermal_properties(
                temperatures=np.linspace(
                    self.input["temperature_start"],
                    self.input["temperature_end"],
                    int(self.input["temperature_steps"]),
                )
            )
            free_energy_lst.append(thermal_properties.free_energies * conversion_factor)
            entropy_lst.append(thermal_properties.entropy * conversion_factor)
            cv_lst.append(thermal_properties.cv * conversion_factor)
            volume_lst.append(job.structure.get_volume() * conversion_factor)

        arg_lst = np.argsort(volume_lst)

        self._output["free_energy"] = np.array(free_energy_lst)[arg_lst]
        self._output["entropy"] = np.array(entropy_lst)[arg_lst]
        self._output["cv"] = np.array(cv_lst)[arg_lst]

        temperature_mesh, volume_mesh = np.meshgrid(
            np.linspace(
                self.input["temperature_start"],
                self.input["temperature_end"],
                int(self.input["temperature_steps"]),
            ),
            np.array(volume_lst)[arg_lst],
        )

        self._output["volumes"] = volume_mesh
        self._output["temperatures"] = temperature_mesh
        with self.project_hdf5.open("output") as hdf5_out:
            for key, val in self._output.items():
                hdf5_out[key] = val

    def optimise_volume(self, bulk_eng):
        """
        Get finite temperature properties.

        Args:
            bulk_eng (numpy.ndarray): array of bulk energies corresponding to the box sizes given
                in the quasi harmonic calculations. For the sake of compatibility, it is strongly
                recommended to use the pyiron Murnaghan class (and make sure that you use the
                same values for `num_points` and `vol_range`).

        Returns:
            volume, free energy, entropy, heat capacity

        The corresponding temperature values can be obtained from `job['output/temperatures'][0]`
        """
        v0_lst, free_eng_lst, entropy_lst, cv_lst = [], [], [], []
        for i, [t, free_energy, cv, entropy, v] in enumerate(
            zip(
                self["output/temperatures"].T,
                self["output/free_energy"].T,
                self["output/cv"].T,
                self["output/entropy"].T,
                self["output/volumes"].T,
            )
        ):
            fit = np.poly1d(
                np.polyfit(
                    v, free_energy + bulk_eng, int(self.input["polynomial_degree"])
                )
            )
            v0, ind = calc_v0_from_fit_funct(
                fit_funct=fit, x=v, save_range=0.0, return_ind=True
            )

            if v0 is not None:
                v0_lst.append(v0)
                free_eng_lst.append(fit(v0))
                entropy_lst.append(entropy[ind])
                cv_lst.append(cv[ind])
            else:
                v0_lst.append(np.nan)
                free_eng_lst.append(np.nan)
                entropy_lst.append(np.nan)
                cv_lst.append(np.nan)

        return v0_lst, free_eng_lst, entropy_lst, cv_lst

    def plot_free_energy_volume_temperature(
        self,
        murnaghan_job=None,
        temperature_start=None,
        temperature_end=None,
        temperature_steps=None,
        color_map="coolwarm",
        axis=None,
        *args,
        **kwargs,
    ):
        """
        Plot volume vs free energy curves for defined temperatures. If no Murnaghan job is assigned, plots free energy without total electronic energy at T=0.

        Args:
            murnaghan_job: job_name or job_id of the Murnaghan job. if None, total electronic energy at T=0 is set to zero.
            temperature_start: if None, value will be used from job.input['temperature_start']
            temperature_end: if None, value will be used from job.input['temperature_end']
            temperature_steps: if None, value will be used from job.input['temperature_steps']
            color_map: colormaps options accessible via matplotlib.cm.get_cmap. Default is 'coolwarm'.
            axis (matplotlib axis, optional): plot to this axis, if not given a new one is created.
            *args: passed through to matplotlib.pyplot.plot when plotting free energies.
            **kwargs: passed through to matplotlib.pyplot.plot when plotting free energies.

        Returns:
            matplib axis: the axis the figure has been drawn to, if axis is given the same object is returned.
        """
        energy_zero = 0
        if murnaghan_job != None:
            job_murn = self.project.load(murnaghan_job)
            energy_zero = job_murn["output/energy"]

        if not self.status.finished:
            raise ValueError(
                "QHA Job must be successfully run, before calling this method."
            )

        if not job_murn.status.finished:
            raise ValueError(
                "Murnaghan Job must be successfully run, before calling this method."
            )

        if axis is None:
            _, axis = matplotlib.pyplot.subplots(1, 1)

        cmap = matplotlib.cm.get_cmap(color_map)

        if temperature_start != None:
            self.input["temperature_start"] = temperature_start
        if temperature_end != None:
            self.input["temperature_end"] = temperature_end
        if temperature_steps != None:
            self.input["temperature_steps"] = temperature_steps
        self.collect_output()

        for i, [t, free_energy, v] in enumerate(
            zip(
                self["output/temperatures"].T,
                self["output/free_energy"].T,
                self["output/volumes"].T,
            )
        ):
            color = cmap(i / len(self["output/temperatures"].T))
            axis.plot(v, free_energy + energy_zero, color=color)

        axis.set_xlabel("Volume")
        axis.set_ylabel("Free Energy")

        temperatures = self["output/temperatures"]
        normalize = matplotlib.colors.Normalize(
            vmin=temperatures.min(), vmax=temperatures.max()
        )
        scalarmappaple = matplotlib.cm.ScalarMappable(norm=normalize, cmap=cmap)
        scalarmappaple.set_array(temperatures)
        cbar = matplotlib.pyplot.colorbar(scalarmappaple, ax=axis)
        cbar.set_label("Temperature")
        return axis

1	# coding: utf-8
2	# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
3	# Distributed under the terms of "New BSD License", see the LICENSE file.
4
5	import matplotlib	1✔
6	import numpy as np	1✔
7
8	from pyiron_atomistics.atomistics.master.murnaghan import MurnaghanJobGenerator	1✔
9	from pyiron_atomistics.atomistics.master.parallel import AtomisticParallelMaster	1✔
10
11	__author__ = "Jan Janssen"	1✔
12	__copyright__ = (	1✔
13	"Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH - "
14	"Computational Materials Design (CM) Department"
15	)
16	__version__ = "0.0.1"	1✔
17	__maintainer__ = "Jan Janssen"	1✔
18	__email__ = "janssen@mpie.de"	1✔
19	__status__ = "development"	1✔
20	__date__ = "Oct 29, 2020"	1✔
21
22
23	def calc_v0_from_fit_funct(fit_funct, x, save_range=0.0, return_ind=False):	1✔
24	fit_funct_der = fit_funct.deriv().r	×
25	fit_funct_der_r = fit_funct_der[fit_funct_der.imag == 0].real	×
26	fit_funct_der_val = fit_funct.deriv(2)(fit_funct_der_r)	×
27	select = (	×
28	(fit_funct_der_val > 0)
29	& (fit_funct_der_r > np.min(x) * (1 - save_range))
30	& (fit_funct_der_r < np.max(x) * (1 + save_range))
31	)
32	v0_lst = fit_funct_der_r[select]	×
33	if len(v0_lst) == 1:	×
34	if return_ind:	×
35	return v0_lst[0], (np.abs(x - v0_lst[0])).argmin()	×
36	else:
37	return v0_lst[0]	×
38	else:
39	select = fit_funct_der_val > 0	×
40	v0_lst = fit_funct_der_r[select]	×
41	if len(v0_lst) == 1:	×
42	if return_ind:	×
43	return v0_lst[0], (np.abs(x - v0_lst[0])).argmin()	×
44	else:
45	return v0_lst[0]	×
46	else:
47	if return_ind:	×
48	return None, None	×
49	else:
50	return None	×
51
52
53	class QuasiHarmonicJob(AtomisticParallelMaster):	1✔
54	"""
55	Obtain finite temperature properties in the framework of quasi harmonic approximation. For the
56	theoretical understanding take a look at the Wikipedia page:
57	https://en.wikipedia.org/wiki/Quasi-harmonic_approximation
58
59	Example:
60
61	>>> pr = Project('my_project')
62	>>> lmp = pr.create_job('Lammps', 'lmp')
63	>>> lmp.structure = structure_of_your_choice
64	>>> phono = lmp.create_job('PhonopyJob', 'phono')
65	>>> qha = phono.create_job('QuasiHarmonicJob', 'qha')
66	>>> qha.run()
67
68	The final results can be obtained through `qha.optimise_volume()`.
69
70	The temperature range defined in the input can be modified afterwards. For this, follow these
71	lines:
72
73	>>> qha.input['temperature_end'] = temperature_end
74	>>> qha.input['temperature_steps'] = temperature_steps
75	>>> qha.input['temperature_start'] = temperature_start
76	>>> qha.collect_output()
77	"""
78
79	def __init__(self, project, job_name="murnaghan"):	1✔
80	"""
81
82	Args:
83	project:
84	job_name:
85	"""
86	super(QuasiHarmonicJob, self).__init__(project, job_name)	1✔
87
88	self.__version__ = "0.0.1"	1✔
89
90	# define default input
91	self.input["num_points"] = (11, "number of sample points")	1✔
92	self.input["vol_range"] = (	1✔
93	0.1,
94	"relative volume variation around volume defined by ref_ham",
95	)
96	self.input["temperature_start"] = 0	1✔
97	self.input["temperature_end"] = 500	1✔
98	self.input["temperature_steps"] = 10	1✔
99	self.input["polynomial_degree"] = 3	1✔
100	self.input["strains"] = (	1✔
101	None,
102	"List of strains that should be calculated. If given vol_range and num_points take no effect.",
103	)
104	self.input["axes"] = (	1✔
105	["x", "y", "z"],
106	"Axes along which the strain will be applied",
107	)
108	self._job_generator = MurnaghanJobGenerator(self)	1✔
109
110	def collect_output(self):	1✔
111	free_energy_lst, entropy_lst, cv_lst, volume_lst = [], [], [], []	×
112	for job_id in self.child_ids:	×
113	job = self.project_hdf5.load(job_id)	×
114	# the underlying phonopy job might create a supercell; if its
115	# reference job is interactive this is reflected in its
116	# job.structure and the output quantities, so we need to rescale
117	# all the output here; if the structure did not change the
118	# conversion_factor will simply be 1.
119	conversion_factor = len(self.structure) / len(job.structure)	×
120	thermal_properties = job.get_thermal_properties(	×
121	temperatures=np.linspace(
122	self.input["temperature_start"],
123	self.input["temperature_end"],
124	int(self.input["temperature_steps"]),
125	)
126	)
127	free_energy_lst.append(thermal_properties.free_energies * conversion_factor)	×
128	entropy_lst.append(thermal_properties.entropy * conversion_factor)	×
129	cv_lst.append(thermal_properties.cv * conversion_factor)	×
130	volume_lst.append(job.structure.get_volume() * conversion_factor)	×
131
132	arg_lst = np.argsort(volume_lst)	×
133
134	self._output["free_energy"] = np.array(free_energy_lst)[arg_lst]	×
135	self._output["entropy"] = np.array(entropy_lst)[arg_lst]	×
136	self._output["cv"] = np.array(cv_lst)[arg_lst]	×
137
138	temperature_mesh, volume_mesh = np.meshgrid(	×
139	np.linspace(
140	self.input["temperature_start"],
141	self.input["temperature_end"],
142	int(self.input["temperature_steps"]),
143	),
144	np.array(volume_lst)[arg_lst],
145	)
146
147	self._output["volumes"] = volume_mesh	×
148	self._output["temperatures"] = temperature_mesh	×
149	with self.project_hdf5.open("output") as hdf5_out:	×
150	for key, val in self._output.items():	×
151	hdf5_out[key] = val	×
152
153	def optimise_volume(self, bulk_eng):	1✔
154	"""
155	Get finite temperature properties.
156
157	Args:
158	bulk_eng (numpy.ndarray): array of bulk energies corresponding to the box sizes given
159	in the quasi harmonic calculations. For the sake of compatibility, it is strongly
160	recommended to use the pyiron Murnaghan class (and make sure that you use the
161	same values for `num_points` and `vol_range`).
162
163	Returns:
164	volume, free energy, entropy, heat capacity
165
166	The corresponding temperature values can be obtained from `job['output/temperatures'][0]`
167	"""
168	v0_lst, free_eng_lst, entropy_lst, cv_lst = [], [], [], []	×
169	for i, [t, free_energy, cv, entropy, v] in enumerate(	×
170	zip(
171	self["output/temperatures"].T,
172	self["output/free_energy"].T,
173	self["output/cv"].T,
174	self["output/entropy"].T,
175	self["output/volumes"].T,
176	)
177	):
178	fit = np.poly1d(	×
179	np.polyfit(
180	v, free_energy + bulk_eng, int(self.input["polynomial_degree"])
181	)
182	)
183	v0, ind = calc_v0_from_fit_funct(	×
184	fit_funct=fit, x=v, save_range=0.0, return_ind=True
185	)
186
187	if v0 is not None:	×
188	v0_lst.append(v0)	×
189	free_eng_lst.append(fit(v0))	×
190	entropy_lst.append(entropy[ind])	×
191	cv_lst.append(cv[ind])	×
192	else:
193	v0_lst.append(np.nan)	×
194	free_eng_lst.append(np.nan)	×
195	entropy_lst.append(np.nan)	×
196	cv_lst.append(np.nan)	×
197
198	return v0_lst, free_eng_lst, entropy_lst, cv_lst	×
199
200	def plot_free_energy_volume_temperature(	1✔
201	self,
202	murnaghan_job=None,
203	temperature_start=None,
204	temperature_end=None,
205	temperature_steps=None,
206	color_map="coolwarm",
207	axis=None,
208	*args,
209	**kwargs,
210	):
211	"""
212	Plot volume vs free energy curves for defined temperatures. If no Murnaghan job is assigned, plots free energy without total electronic energy at T=0.
213
214	Args:
215	murnaghan_job: job_name or job_id of the Murnaghan job. if None, total electronic energy at T=0 is set to zero.
216	temperature_start: if None, value will be used from job.input['temperature_start']
217	temperature_end: if None, value will be used from job.input['temperature_end']
218	temperature_steps: if None, value will be used from job.input['temperature_steps']
219	color_map: colormaps options accessible via matplotlib.cm.get_cmap. Default is 'coolwarm'.
220	axis (matplotlib axis, optional): plot to this axis, if not given a new one is created.
221	*args: passed through to matplotlib.pyplot.plot when plotting free energies.
222	**kwargs: passed through to matplotlib.pyplot.plot when plotting free energies.
223
224	Returns:
225	matplib axis: the axis the figure has been drawn to, if axis is given the same object is returned.
226	"""
227	energy_zero = 0	×
228	if murnaghan_job != None:	×
229	job_murn = self.project.load(murnaghan_job)	×
230	energy_zero = job_murn["output/energy"]	×
231
232	if not self.status.finished:	×
233	raise ValueError(	×
234	"QHA Job must be successfully run, before calling this method."
235	)
236
237	if not job_murn.status.finished:	×
238	raise ValueError(	×
239	"Murnaghan Job must be successfully run, before calling this method."
240	)
241
242	if axis is None:	×
243	_, axis = matplotlib.pyplot.subplots(1, 1)	×
244
245	cmap = matplotlib.cm.get_cmap(color_map)	×
246
247	if temperature_start != None:	×
248	self.input["temperature_start"] = temperature_start	×
249	if temperature_end != None:	×
250	self.input["temperature_end"] = temperature_end	×
251	if temperature_steps != None:	×
252	self.input["temperature_steps"] = temperature_steps	×
253	self.collect_output()	×
254
255	for i, [t, free_energy, v] in enumerate(	×
256	zip(
257	self["output/temperatures"].T,
258	self["output/free_energy"].T,
259	self["output/volumes"].T,
260	)
261	):
262	color = cmap(i / len(self["output/temperatures"].T))	×
263	axis.plot(v, free_energy + energy_zero, color=color)	×
264
265	axis.set_xlabel("Volume")	×
266	axis.set_ylabel("Free Energy")	×
267
268	temperatures = self["output/temperatures"]	×
269	normalize = matplotlib.colors.Normalize(	×
270	vmin=temperatures.min(), vmax=temperatures.max()
271	)
272	scalarmappaple = matplotlib.cm.ScalarMappable(norm=normalize, cmap=cmap)	×
273	scalarmappaple.set_array(temperatures)	×
NEW 274	cbar = matplotlib.pyplot.colorbar(scalarmappaple, ax=axis)	×
275	cbar.set_label("Temperature")	×
276	return axis	×

pyiron / pyiron_atomistics / 10736148419

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous