Skip to content

API documentation

Module to parse and store data from electronic structure calculations. Contains the parent class Calculation to store data from a variety of sources. Contains the child class AimsCalculation to read and store data from a FHI-aims calculation.

AimsCalculation

Bases: Calculation

Class for parsing and storing data from a FHI-AIMS total energy calculation.

Example:

BaS_calc = AimsCalculation("./aims_output/output.aims")

Attributes:

volume (float): volume of the periodic unit cell in Angstrom^3
filepath (str): path to the calculation output files
energy (float): DFT total energy in eV
xc (str): XC functional used to calculate the total energy
NAtoms (int): number of atoms in the periodic unit cell
Source code in thermopot/calculations.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
class AimsCalculation(Calculation):
    """Class for parsing and storing data from a FHI-AIMS total energy calculation.

    Example:

       BaS_calc = AimsCalculation("./aims_output/output.aims")

    Attributes:

        volume (float): volume of the periodic unit cell in Angstrom^3
        filepath (str): path to the calculation output files
        energy (float): DFT total energy in eV
        xc (str): XC functional used to calculate the total energy
        NAtoms (int): number of atoms in the periodic unit cell
    """

    def __init__(self, filepath="./calculation.out", gas=False):
        """
        Args:

            filepath (str): path to the calculation output files
            gas (bool): True if gas species, False otherwise

        Note:

            If gas is True then volume is None
        """
        super().__init__()
        self.filepath = filepath
        if not gas:
            self.volume = self.get_volume()
        self.energy = self.get_energy()
        self.xc = self.get_xc()
        self.NAtoms = self.get_NAtoms()

    def get_volume(self):
        """
        Returns:

            (float): volume of the periodic unit cell in Angstrom^3
        """
        with open(self.filepath) as contents:
            return float(
                re.findall("Unit cell volume\s+:\s*(.*)\sA", contents.read())[-1]
            )

    def get_energy(self):
        """
        Returns:

            (float): DFT total energy in eV
        """
        with open(self.filepath) as contents:
            return float(
                re.findall(
                    "Total energy of the DFT[^0-9]+(-\d*\.?\d*) eV", contents.read()
                )[-1]
            )

    def get_xc(self):
        """
        Returns:

            (str): XC functional used to calculate the total energy
        """
        with open(self.filepath) as contents:
            return re.findall("xc               (\S+)", contents.read())[0]

    def get_NAtoms(self):
        """
        Returns:

            (int): number of atoms in the periodic unit cell
        """
        with open(self.filepath) as contents:
            return int(
                re.findall("Number of atoms\s +:\s + (\S+)", contents.read())[-1]
            )

__init__(filepath='./calculation.out', gas=False)

Args:

filepath (str): path to the calculation output files
gas (bool): True if gas species, False otherwise

Note:

If gas is True then volume is None
Source code in thermopot/calculations.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def __init__(self, filepath="./calculation.out", gas=False):
    """
    Args:

        filepath (str): path to the calculation output files
        gas (bool): True if gas species, False otherwise

    Note:

        If gas is True then volume is None
    """
    super().__init__()
    self.filepath = filepath
    if not gas:
        self.volume = self.get_volume()
    self.energy = self.get_energy()
    self.xc = self.get_xc()
    self.NAtoms = self.get_NAtoms()

get_NAtoms()

Returns:

(int): number of atoms in the periodic unit cell
Source code in thermopot/calculations.py 
145
146
147
148
149
150
151
152
153
154
def get_NAtoms(self):
    """
    Returns:

        (int): number of atoms in the periodic unit cell
    """
    with open(self.filepath) as contents:
        return int(
            re.findall("Number of atoms\s +:\s + (\S+)", contents.read())[-1]
        )

get_energy()

Returns:

(float): DFT total energy in eV
Source code in thermopot/calculations.py 
123
124
125
126
127
128
129
130
131
132
133
134
def get_energy(self):
    """
    Returns:

        (float): DFT total energy in eV
    """
    with open(self.filepath) as contents:
        return float(
            re.findall(
                "Total energy of the DFT[^0-9]+(-\d*\.?\d*) eV", contents.read()
            )[-1]
        )

get_volume()

Returns:

(float): volume of the periodic unit cell in Angstrom^3
Source code in thermopot/calculations.py 
112
113
114
115
116
117
118
119
120
121
def get_volume(self):
    """
    Returns:

        (float): volume of the periodic unit cell in Angstrom^3
    """
    with open(self.filepath) as contents:
        return float(
            re.findall("Unit cell volume\s+:\s*(.*)\sA", contents.read())[-1]
        )

get_xc()

Returns:

(str): XC functional used to calculate the total energy
Source code in thermopot/calculations.py 
136
137
138
139
140
141
142
143
def get_xc(self):
    """
    Returns:

        (str): XC functional used to calculate the total energy
    """
    with open(self.filepath) as contents:
        return re.findall("xc               (\S+)", contents.read())[0]

Calculation

Parent class for parsing and storing data from electronic structure calculations.

Example:

BaS_calc = Calculation(volume=63.2552, energy=-235926.586148547, xc='pbesol', NAtoms=2)

Attributes:

volume (float): volume of the periodic unit cell in Angstrom^3
filepath (str): path to the calculation output files
energy (float): DFT total energy in eV
xc (str): XC functional used to calculate the total energy
NAtoms (int): number of atoms in the periodic unit cell

Note:

If gas is True then no volume attribute is required.
Source code in thermopot/calculations.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
class Calculation:
    """
    Parent class for parsing and storing data from electronic structure calculations.

    Example:

        BaS_calc = Calculation(volume=63.2552, energy=-235926.586148547, xc='pbesol', NAtoms=2)

    Attributes:

        volume (float): volume of the periodic unit cell in Angstrom^3
        filepath (str): path to the calculation output files
        energy (float): DFT total energy in eV
        xc (str): XC functional used to calculate the total energy
        NAtoms (int): number of atoms in the periodic unit cell

    Note:

        If gas is True then no volume attribute is required.
    """

    def __init__(
        self, energy=None, xc=None, NAtoms=None, volume=None, filepath=None, gas=False
    ):
        """
        Note:

            All attributes are None until set by derived classes or specified by user.

        Args:

            volume (float): volume of the periodic unit cell in Angstrom^3
            filepath (str, optional): path to the calculation output files
            energy (float): DFT total energy in eV
            xc (str): XC functional used to calculate the total energy
            NAtoms (int): number of atoms in the periodic unit cell
            gas (bool): True if gas species, False otherwise

        Note:

            If gas is True then volume is None
        """
        if not gas:
            self.volume = volume
        self.filepath = filepath
        self.energy = energy
        self.xc = xc
        self.NAtoms = NAtoms

        # self.check_attributes()

    def check_attributes(self):
        """Check that the Calculation class attributes make basic sense."""

        assert (
            type(self.filepath) == str or self.filepath is None
        ), "filepath must be a string"
        assert type(self.energy) == float, "energy must be a float"
        assert type(self.xc) == str, "xc must be a string"
        assert (
            type(self.NAtoms) == int
        ) and self.NAtoms >= 1, "NAtoms must be an integer >= 1"
        assert (
            type(self.volume) == float
        ) and self.volume > 0, "volume must be a float > 0"

__init__(energy=None, xc=None, NAtoms=None, volume=None, filepath=None, gas=False)

Note:

All attributes are None until set by derived classes or specified by user.

Args:

volume (float): volume of the periodic unit cell in Angstrom^3
filepath (str, optional): path to the calculation output files
energy (float): DFT total energy in eV
xc (str): XC functional used to calculate the total energy
NAtoms (int): number of atoms in the periodic unit cell
gas (bool): True if gas species, False otherwise

Note:

If gas is True then volume is None
Source code in thermopot/calculations.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def __init__(
    self, energy=None, xc=None, NAtoms=None, volume=None, filepath=None, gas=False
):
    """
    Note:

        All attributes are None until set by derived classes or specified by user.

    Args:

        volume (float): volume of the periodic unit cell in Angstrom^3
        filepath (str, optional): path to the calculation output files
        energy (float): DFT total energy in eV
        xc (str): XC functional used to calculate the total energy
        NAtoms (int): number of atoms in the periodic unit cell
        gas (bool): True if gas species, False otherwise

    Note:

        If gas is True then volume is None
    """
    if not gas:
        self.volume = volume
    self.filepath = filepath
    self.energy = energy
    self.xc = xc
    self.NAtoms = NAtoms

check_attributes()

Check that the Calculation class attributes make basic sense.

Source code in thermopot/calculations.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def check_attributes(self):
    """Check that the Calculation class attributes make basic sense."""

    assert (
        type(self.filepath) == str or self.filepath is None
    ), "filepath must be a string"
    assert type(self.energy) == float, "energy must be a float"
    assert type(self.xc) == str, "xc must be a string"
    assert (
        type(self.NAtoms) == int
    ) and self.NAtoms >= 1, "NAtoms must be an integer >= 1"
    assert (
        type(self.volume) == float
    ) and self.volume > 0, "volume must be a float > 0"

Module contains the classes Sulfur_model, Solid and IdealGas to store basic material data. Each class provides methods for calculating various thermodynamic properties.

IdealGas

Bases: Material

Class for ideal gas properties.

Example:

S2_gas = materials.IdealGas("S2", {"S":2}, "./thermo_data/S2",calculation=S2_gas_calc)

Attributes:

name (str): Identifying string stoichiometry (dict): relates element to the number of atoms in a single formula unit thermo_dile (str): path to the thermodynamics data calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class energies (dict, optional): relates xc functional to DFT total energy in eV zpe_pbesol (float, optional): zero point energy calculated using the pbesol XC-functional zpe_hse06 (float, optional): zero point energy calculated using the hse06 XC-functional zpe_lit (float, optional): zero point energy calculated using literature values

Note:

Ideal gas law PV=nRT is applied: specifically (dH/dP) at const. T = 0 and int(mu)^P2_P1 dP = kTln(P2/P1).
Enthalpy has no P dependence as volume is not restricted / expansion step is defined as isothermal
Source code in thermopot/materials.py 
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
class IdealGas(Material):
    """
    Class for ideal gas properties.

    Example:

        S2_gas = materials.IdealGas("S2", {"S":2}, "./thermo_data/S2",calculation=S2_gas_calc)

    Attributes:

       name (str): Identifying string
       stoichiometry (dict): relates element to the number of atoms in a single formula unit
       thermo_dile (str): path to the thermodynamics data
       calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class
       energies (dict, optional): relates xc functional to DFT total energy in eV
       zpe_pbesol (float, optional): zero point energy calculated using the pbesol XC-functional
       zpe_hse06 (float, optional): zero point energy calculated using the hse06 XC-functional
       zpe_lit (float, optional): zero point energy calculated using literature values

    Note:

        Ideal gas law PV=nRT is applied: specifically (dH/dP) at const. T = 0 and int(mu)^P2_P1 dP = kTln(P2/P1).
        Enthalpy has no P dependence as volume is not restricted / expansion step is defined as isothermal
    """

    def __init__(
        self,
        name,
        stoichiometry,
        thermo_file,
        calculation=False,
        energies=False,
        zpe_pbesol=0,
        zpe_hse06=0,
        zpe_lit=0,
    ):
        """
        Args:

        name (str): Identifying string
        stoichiometry (dict): relates element to the number of atoms in a single formula unit
        thermo_dile (str): path to the thermodynamics data
        calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class
        energies (dict, optional): relates xc functional to DFT total energy in eV
        zpe_pbesol (float, optional): zero point energy calculated using the pbesol XC-functional
        zpe_hse06 (float, optional): zero point energy calculated using the hse06 XC-functional
        zpe_lit (float, optional): zero point energy calculated using literature values
        """
        if calculation is not False:
            Material.__init__(
                self, name, stoichiometry, {calculation.xc: calculation.energy}
            )
        else:
            Material.__init__(self, name, stoichiometry, energies)
        self.thermo_file = materials_directory + thermo_file

        if zpe_hse06 > 0:
            self.zpe = zpe_pbesol
        elif zpe_pbesol > 0:
            self.zpe = zpe_pbesol
        elif zpe_lit > 0:
            self.zpe = zpe_lit
        else:
            self.zpe = 0

    def U(self, T, xc="pbesol", units="eV"):
        """
        Calculates the internal energy of one formula unit of ideal gas.

        Example:

            U = S2_gas.U(300,xc="pbesol",units="eV")

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Returns:

            U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of gas, or a single internal energy float when a single temperature is passed as an argument.

        """
        U_func = interpolate.get_potential_nist_table(self.thermo_file, "U")
        E_dft = self.energies[xc]
        U_eV = (
            E_dft
            + self.zpe
            + U_func(T)
            * constants.physical_constants["joule-electron volt relationship"][0]
            / constants.N_A
        )

        if units == "eV":
            return U_eV

        elif units == "J":
            return (
                U_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                U_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

    def H(self, T, xc="pbesol", units="eV"):
        """
        Calculates the Enthalpy of one formula unit of ideal gas.

        Examples:

            H = S2_gas.H(300,xc="pbesol",units="eV")
            H = S2_gas.H(np.linspace(100,700,1000))

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Returns:

            H (float/ndarray):  1D Numpy array (with the same dimensions as T) containing the enthalpy of one formula unit of gas, or a single enthalpy float when a single temperature is passed as an argument.
        """
        H_func = interpolate.get_potential_nist_table(self.thermo_file, "H")
        E_dft = self.energies[xc]
        H_eV = (
            E_dft
            + self.zpe
            + H_func(T)
            * constants.physical_constants["joule-electron volt relationship"][0]
            / constants.N_A
        )

        if units == "eV":
            return H_eV

        elif units == "J":
            return (
                H_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                H_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

    def mu(self, T, P, xc="pbesol", units="eV"):
        """
        Calculates the Gibbs Free Energy of one formula unit of ideal gas.

        Examples:

            mu = S2_gas.mu(300,xc="pbesol",units="eV")
            mu = S2_gas.mu(np.linspace(100,700,1000))

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Note:

            T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
            Other T, P arrays will result in undefined behaviour.

        Returns:

            mu (float/ndarray): Gibbs Free Energy of one formula unit of ideal gas expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
        """

        S_func = interpolate.get_potential_nist_table(self.thermo_file, "S")
        S = (
            S_func(T)
            * constants.physical_constants["joule-electron volt relationship"][0]
            / constants.N_A
        )
        H = self.H(T, xc=xc)
        mu_eV = (
            H
            - T * S
            + constants.physical_constants["Boltzmann constant in eV/K"][0]
            * T
            * np.log(P / 1e5)
        )

        if units == "eV":
            return mu_eV

        elif units == "J":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

H(T, xc='pbesol', units='eV')

Calculates the Enthalpy of one formula unit of ideal gas.

Examples:

H = S2_gas.H(300,xc="pbesol",units="eV")
H = S2_gas.H(np.linspace(100,700,1000))

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Returns:

H (float/ndarray):  1D Numpy array (with the same dimensions as T) containing the enthalpy of one formula unit of gas, or a single enthalpy float when a single temperature is passed as an argument.
Source code in thermopot/materials.py 
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
def H(self, T, xc="pbesol", units="eV"):
    """
    Calculates the Enthalpy of one formula unit of ideal gas.

    Examples:

        H = S2_gas.H(300,xc="pbesol",units="eV")
        H = S2_gas.H(np.linspace(100,700,1000))

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Returns:

        H (float/ndarray):  1D Numpy array (with the same dimensions as T) containing the enthalpy of one formula unit of gas, or a single enthalpy float when a single temperature is passed as an argument.
    """
    H_func = interpolate.get_potential_nist_table(self.thermo_file, "H")
    E_dft = self.energies[xc]
    H_eV = (
        E_dft
        + self.zpe
        + H_func(T)
        * constants.physical_constants["joule-electron volt relationship"][0]
        / constants.N_A
    )

    if units == "eV":
        return H_eV

    elif units == "J":
        return (
            H_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            H_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

U(T, xc='pbesol', units='eV')

Calculates the internal energy of one formula unit of ideal gas.

Example:

U = S2_gas.U(300,xc="pbesol",units="eV")

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Returns:

U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of gas, or a single internal energy float when a single temperature is passed as an argument.
Source code in thermopot/materials.py 
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
def U(self, T, xc="pbesol", units="eV"):
    """
    Calculates the internal energy of one formula unit of ideal gas.

    Example:

        U = S2_gas.U(300,xc="pbesol",units="eV")

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Returns:

        U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of gas, or a single internal energy float when a single temperature is passed as an argument.

    """
    U_func = interpolate.get_potential_nist_table(self.thermo_file, "U")
    E_dft = self.energies[xc]
    U_eV = (
        E_dft
        + self.zpe
        + U_func(T)
        * constants.physical_constants["joule-electron volt relationship"][0]
        / constants.N_A
    )

    if units == "eV":
        return U_eV

    elif units == "J":
        return (
            U_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            U_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

__init__(name, stoichiometry, thermo_file, calculation=False, energies=False, zpe_pbesol=0, zpe_hse06=0, zpe_lit=0)

Args:

name (str): Identifying string stoichiometry (dict): relates element to the number of atoms in a single formula unit thermo_dile (str): path to the thermodynamics data calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class energies (dict, optional): relates xc functional to DFT total energy in eV zpe_pbesol (float, optional): zero point energy calculated using the pbesol XC-functional zpe_hse06 (float, optional): zero point energy calculated using the hse06 XC-functional zpe_lit (float, optional): zero point energy calculated using literature values

Source code in thermopot/materials.py 
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
def __init__(
    self,
    name,
    stoichiometry,
    thermo_file,
    calculation=False,
    energies=False,
    zpe_pbesol=0,
    zpe_hse06=0,
    zpe_lit=0,
):
    """
    Args:

    name (str): Identifying string
    stoichiometry (dict): relates element to the number of atoms in a single formula unit
    thermo_dile (str): path to the thermodynamics data
    calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class
    energies (dict, optional): relates xc functional to DFT total energy in eV
    zpe_pbesol (float, optional): zero point energy calculated using the pbesol XC-functional
    zpe_hse06 (float, optional): zero point energy calculated using the hse06 XC-functional
    zpe_lit (float, optional): zero point energy calculated using literature values
    """
    if calculation is not False:
        Material.__init__(
            self, name, stoichiometry, {calculation.xc: calculation.energy}
        )
    else:
        Material.__init__(self, name, stoichiometry, energies)
    self.thermo_file = materials_directory + thermo_file

    if zpe_hse06 > 0:
        self.zpe = zpe_pbesol
    elif zpe_pbesol > 0:
        self.zpe = zpe_pbesol
    elif zpe_lit > 0:
        self.zpe = zpe_lit
    else:
        self.zpe = 0

mu(T, P, xc='pbesol', units='eV')

Calculates the Gibbs Free Energy of one formula unit of ideal gas.

Examples:

mu = S2_gas.mu(300,xc="pbesol",units="eV")
mu = S2_gas.mu(np.linspace(100,700,1000))

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Note:

T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
Other T, P arrays will result in undefined behaviour.

Returns:

mu (float/ndarray): Gibbs Free Energy of one formula unit of ideal gas expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
Source code in thermopot/materials.py 
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
def mu(self, T, P, xc="pbesol", units="eV"):
    """
    Calculates the Gibbs Free Energy of one formula unit of ideal gas.

    Examples:

        mu = S2_gas.mu(300,xc="pbesol",units="eV")
        mu = S2_gas.mu(np.linspace(100,700,1000))

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Note:

        T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
        Other T, P arrays will result in undefined behaviour.

    Returns:

        mu (float/ndarray): Gibbs Free Energy of one formula unit of ideal gas expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
    """

    S_func = interpolate.get_potential_nist_table(self.thermo_file, "S")
    S = (
        S_func(T)
        * constants.physical_constants["joule-electron volt relationship"][0]
        / constants.N_A
    )
    H = self.H(T, xc=xc)
    mu_eV = (
        H
        - T * S
        + constants.physical_constants["Boltzmann constant in eV/K"][0]
        * T
        * np.log(P / 1e5)
    )

    if units == "eV":
        return mu_eV

    elif units == "J":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

Material

Bases: object

Parent class for storing materials properties.

Attributes:

name (str): Identifying string
stoichiometry (dict): relates element to the number of atoms in a single formula unit
energies (dict): relates xc functional to DFT total energy in eV
N (int): number of atoms per formula unit
Source code in thermopot/materials.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
class Material(object):
    """
    Parent class for storing materials properties.

    Attributes:

        name (str): Identifying string
        stoichiometry (dict): relates element to the number of atoms in a single formula unit
        energies (dict): relates xc functional to DFT total energy in eV
        N (int): number of atoms per formula unit
    """

    def __init__(self, name, stoichiometry, energies):
        self.name = name
        self.stoichiometry = stoichiometry
        self.energies = energies
        self.N = sum(self.stoichiometry.values())

Solid

Bases: Material

Class for solid material data.

Note:

The material is assumed to be incompressible and without thermal expansion.

Example:

BaS = Solid('BaS',{'Ba':1,'S':1},"./phonon_data/Ba_S",calculation=BaS_calc)

Attributes:

name (str): Identifying string
stoichiometry (dict): relates element to the number of atoms in a single formula unit
energies (dict): relates xc functional to DFT total energy in eV
N (int): number of atoms per formula unit
fu_cell (int): number of formula units in periodic unit cell
volume (float): volume of unit cell in Angstroms^3
phonon_filepath (str): path to the phonon output data
NAtoms (int): number of atoms in periodic unit cell
Source code in thermopot/materials.py 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
class Solid(Material):
    """
    Class for solid material data.

    Note:

        The material is assumed to be incompressible and without thermal expansion.

    Example:

        BaS = Solid('BaS',{'Ba':1,'S':1},"./phonon_data/Ba_S",calculation=BaS_calc)

    Attributes:

        name (str): Identifying string
        stoichiometry (dict): relates element to the number of atoms in a single formula unit
        energies (dict): relates xc functional to DFT total energy in eV
        N (int): number of atoms per formula unit
        fu_cell (int): number of formula units in periodic unit cell
        volume (float): volume of unit cell in Angstroms^3
        phonon_filepath (str): path to the phonon output data
        NAtoms (int): number of atoms in periodic unit cell
    """

    def __init__(
        self,
        name,
        stoichiometry,
        phonon_filepath,
        calculation=False,
        volume=False,
        energies=False,
        NAtoms=1,
    ):
        """
        Args:

           name (str): Identifying string
           stoichiometry (dict): relates element to the number of atoms in a single formula unit
           phonon_filepath (str): path to the phonon output data
           calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class
           volume (float, optional): volume of unit cell in Angstroms^3
           energies (dict, optional): relates xc functional to DFT total energy in eV
           NAtoms (int): number of atoms in periodic unit cell
        """

        if calculation is not False:
            if type(calculation) is not list:
                Material.__init__(
                    self, name, stoichiometry, {calculation.xc: calculation.energy}
                )
                self.volume = calculation.volume
                self.NAtoms = calculation.NAtoms
            else:
                pass

        else:
            Material.__init__(self, name, stoichiometry, energies)

            self.NAtoms = NAtoms
            self.volume = volume

        self.fu_cell = self.NAtoms / self.N
        self.phonon_filepath = materials_directory + phonon_filepath

    def U(self, T, xc="pbesol", units="eV"):
        """
        Calculates the internal energy of one formula unit of solid.

        Example:

            U = BaS.U(300,xc="pbesol",units="eV")

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Returns:

            U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of solid, or a single internal energy float when a single temperature is passed as an argument.

        """
        U_func = interpolate.get_potential_aims(self.phonon_filepath, "U")
        E_dft = self.energies[xc]

        U_eV = (E_dft + U_func(T)) / self.fu_cell

        if units == "eV":
            return U_eV

        elif units == "J":
            return (
                U_eV(T, xc=xc)
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                U_eV(T, xc=xc)
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

    def H(self, T, P, xc="pbesol", units="eV"):
        """
        Calculates the Enthalpy (H = U + PV) of one formula unit of solid.

        Examples:

            H = BaS.H(300,1E3,xc="pbesol",units="eV")
            H = BaS.H(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Note:

            T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case H is an m x n matrix.
            Other T, P arrays will result in undefined behaviour.

        Returns:

            H (float/ndarray): Enthalpy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
        """
        U_func = interpolate.get_potential_aims(self.phonon_filepath, "U")
        PV = (
            P
            * self.volume
            * 1e-30
            * constants.physical_constants["joule-electron volt relationship"][0]
            / constants.N_A
        )
        E_dft = self.energies[xc]
        H_eV = ((E_dft + U_func(T)) + PV) / self.fu_cell

        if units == "eV":
            return H_eV

        elif units == "J":
            return (
                H_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                H_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

    def mu(self, T, P, xc="pbesol", units="eV"):
        """
        Calculates the Gibbs Free Energy (mu = U + PV - TS) of one formula unit of solid.

        Examples:

            mu = BaS.mu(300,1E3,xc="pbesol",units="eV")
            mu = BaS.mu(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Note:

            T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
            Other T, P arrays will result in undefined behaviour.

        Returns:

            mu (float/ndarray): Gibbs Free Energy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
        """
        TS_func = interpolate.get_potential_aims(self.phonon_filepath, "TS")
        H = self.H(T, P, xc=xc)
        mu_eV = H - (TS_func(T)) / self.fu_cell

        if units == "eV":
            return mu_eV

        elif units == "J":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

    def Cv(self, T, units="kB"):
        """
        Calculates the Constant-volume heat capacity of one formula unit of solid.

        Examples:

            Cv = BaS.mu(300,xc="pbesol",units="eV")
            Cv = BaS.mu(np.linspace(100,700,1000),xc="pbesol",units="kJ")

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            xc (str, optional): DFT XC functional used to calculate the ground state energy
            units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

        Returns:

            Cv (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the Constant-volume heat capacity of one formula unit of solid, or a single heat capacity float when a single temperature is passed as an argument.
        """
        Cv_func = interpolate.get_potential_aims(self.phonon_filepath, "Cv")
        Cv_kB = Cv_func(T) / self.fu_cell

        if units == "kB":
            return Cv_kB

        elif units == "eV":
            return Cv_kB * constants.physical_constants["Boltzmann constant in eV/K"][0]

        elif units == "J":
            return (
                Cv_kB
                * constants.physical_constants["Boltzmann constant"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                Cv_kB
                * constants.physical_constants["Boltzmann constant"][0]
                * constants.N_A
                * 0.001
            )

Cv(T, units='kB')

Calculates the Constant-volume heat capacity of one formula unit of solid.

Examples:

Cv = BaS.mu(300,xc="pbesol",units="eV")
Cv = BaS.mu(np.linspace(100,700,1000),xc="pbesol",units="kJ")

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Returns:

Cv (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the Constant-volume heat capacity of one formula unit of solid, or a single heat capacity float when a single temperature is passed as an argument.
Source code in thermopot/materials.py 
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def Cv(self, T, units="kB"):
    """
    Calculates the Constant-volume heat capacity of one formula unit of solid.

    Examples:

        Cv = BaS.mu(300,xc="pbesol",units="eV")
        Cv = BaS.mu(np.linspace(100,700,1000),xc="pbesol",units="kJ")

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Returns:

        Cv (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the Constant-volume heat capacity of one formula unit of solid, or a single heat capacity float when a single temperature is passed as an argument.
    """
    Cv_func = interpolate.get_potential_aims(self.phonon_filepath, "Cv")
    Cv_kB = Cv_func(T) / self.fu_cell

    if units == "kB":
        return Cv_kB

    elif units == "eV":
        return Cv_kB * constants.physical_constants["Boltzmann constant in eV/K"][0]

    elif units == "J":
        return (
            Cv_kB
            * constants.physical_constants["Boltzmann constant"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            Cv_kB
            * constants.physical_constants["Boltzmann constant"][0]
            * constants.N_A
            * 0.001
        )

H(T, P, xc='pbesol', units='eV')

Calculates the Enthalpy (H = U + PV) of one formula unit of solid.

Examples:

H = BaS.H(300,1E3,xc="pbesol",units="eV")
H = BaS.H(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Note:

T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case H is an m x n matrix.
Other T, P arrays will result in undefined behaviour.

Returns:

H (float/ndarray): Enthalpy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
Source code in thermopot/materials.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
def H(self, T, P, xc="pbesol", units="eV"):
    """
    Calculates the Enthalpy (H = U + PV) of one formula unit of solid.

    Examples:

        H = BaS.H(300,1E3,xc="pbesol",units="eV")
        H = BaS.H(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Note:

        T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case H is an m x n matrix.
        Other T, P arrays will result in undefined behaviour.

    Returns:

        H (float/ndarray): Enthalpy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
    """
    U_func = interpolate.get_potential_aims(self.phonon_filepath, "U")
    PV = (
        P
        * self.volume
        * 1e-30
        * constants.physical_constants["joule-electron volt relationship"][0]
        / constants.N_A
    )
    E_dft = self.energies[xc]
    H_eV = ((E_dft + U_func(T)) + PV) / self.fu_cell

    if units == "eV":
        return H_eV

    elif units == "J":
        return (
            H_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            H_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

U(T, xc='pbesol', units='eV')

Calculates the internal energy of one formula unit of solid.

Example:

U = BaS.U(300,xc="pbesol",units="eV")

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Returns:

U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of solid, or a single internal energy float when a single temperature is passed as an argument.
Source code in thermopot/materials.py 
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
def U(self, T, xc="pbesol", units="eV"):
    """
    Calculates the internal energy of one formula unit of solid.

    Example:

        U = BaS.U(300,xc="pbesol",units="eV")

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Returns:

        U (float/ndarray): 1D Numpy array (with the same dimensions as T) containing the internal energies of one formula unit of solid, or a single internal energy float when a single temperature is passed as an argument.

    """
    U_func = interpolate.get_potential_aims(self.phonon_filepath, "U")
    E_dft = self.energies[xc]

    U_eV = (E_dft + U_func(T)) / self.fu_cell

    if units == "eV":
        return U_eV

    elif units == "J":
        return (
            U_eV(T, xc=xc)
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            U_eV(T, xc=xc)
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

__init__(name, stoichiometry, phonon_filepath, calculation=False, volume=False, energies=False, NAtoms=1)

Args:

name (str): Identifying string stoichiometry (dict): relates element to the number of atoms in a single formula unit phonon_filepath (str): path to the phonon output data calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class volume (float, optional): volume of unit cell in Angstroms^3 energies (dict, optional): relates xc functional to DFT total energy in eV NAtoms (int): number of atoms in periodic unit cell

Source code in thermopot/materials.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def __init__(
    self,
    name,
    stoichiometry,
    phonon_filepath,
    calculation=False,
    volume=False,
    energies=False,
    NAtoms=1,
):
    """
    Args:

       name (str): Identifying string
       stoichiometry (dict): relates element to the number of atoms in a single formula unit
       phonon_filepath (str): path to the phonon output data
       calculation (thermopot.calculation.Calculation, optional): instance of the thermopot.calculation.Calculation class
       volume (float, optional): volume of unit cell in Angstroms^3
       energies (dict, optional): relates xc functional to DFT total energy in eV
       NAtoms (int): number of atoms in periodic unit cell
    """

    if calculation is not False:
        if type(calculation) is not list:
            Material.__init__(
                self, name, stoichiometry, {calculation.xc: calculation.energy}
            )
            self.volume = calculation.volume
            self.NAtoms = calculation.NAtoms
        else:
            pass

    else:
        Material.__init__(self, name, stoichiometry, energies)

        self.NAtoms = NAtoms
        self.volume = volume

    self.fu_cell = self.NAtoms / self.N
    self.phonon_filepath = materials_directory + phonon_filepath

mu(T, P, xc='pbesol', units='eV')

Calculates the Gibbs Free Energy (mu = U + PV - TS) of one formula unit of solid.

Examples:

mu = BaS.mu(300,1E3,xc="pbesol",units="eV")
mu = BaS.mu(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
xc (str, optional): DFT XC functional used to calculate the ground state energy
units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

Note:

T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
Other T, P arrays will result in undefined behaviour.

Returns:

mu (float/ndarray): Gibbs Free Energy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
Source code in thermopot/materials.py 
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
def mu(self, T, P, xc="pbesol", units="eV"):
    """
    Calculates the Gibbs Free Energy (mu = U + PV - TS) of one formula unit of solid.

    Examples:

        mu = BaS.mu(300,1E3,xc="pbesol",units="eV")
        mu = BaS.mu(np.linspace(100,700,1000),np.array(np.logspace(1, 7, 100),ndmin=2).transpose())

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.
        xc (str, optional): DFT XC functional used to calculate the ground state energy
        units (str, optional):  specifies the units as "eV", "J" (J/mol) or "kJ" (kJ/mol)

    Note:

        T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
        Other T, P arrays will result in undefined behaviour.

    Returns:

        mu (float/ndarray): Gibbs Free Energy of one formula unit of solid expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
    """
    TS_func = interpolate.get_potential_aims(self.phonon_filepath, "TS")
    H = self.H(T, P, xc=xc)
    mu_eV = H - (TS_func(T)) / self.fu_cell

    if units == "eV":
        return mu_eV

    elif units == "J":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

Sulfur_model

Bases: object

Class with parameterised model for sulfur chemical potential. From work of Jackson et al, https://doi.org/10.1039/C5SC03088A. Region of validity is 400 - 1500 K, 10^0 - 10^7 Pa. You must provide a reference energy from e.g. a DFT total energy calculation.

Source code in thermopot/materials.py 
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
class Sulfur_model(object):
    """
    Class with parameterised model for sulfur chemical potential.
    From work of Jackson et al, https://doi.org/10.1039/C5SC03088A.
    Region of validity is 400 - 1500 K, 10^0 - 10^7 Pa.
    You must provide a reference energy from e.g. a DFT total energy calculation.
    """

    def __init__(self, reference_energy):
        self.reference_energy = reference_energy

    def mu(self, T, P, units="eV", xc=None):
        """
        Returns the chemical potential of one atom of Sulfur

        Args:

            T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
            P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.

        Note:

            T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
            Other T, P arrays will result in undefined behaviour.

        Returns:

            mu (float/ndarray): Chemical potential of one sulfur atom expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
        """

        Kb = scipy.constants.physical_constants["Boltzmann constant in eV/K"][0]

        if (
            np.any(T > 1500)
            or np.any(T < 400)
            or np.any(P > 10**7)
            or np.any(P < 10**0)
        ):
            print(
                """WARNING!: You are using the sulfur model beyond the temperature and/or pressure range it was fitted to.
                 Region of validity is 400 - 1500 K, 10^0 - 10^7 Pa. """
            )

        def T_tr(P):
            return (
                5.077e2
                + 7.272e1 * np.log10(P)
                - 8.295 * np.log10(P) ** 2
                + 1.828 * np.log10(P) ** 3
            )

        def mu_S_2(T, P):
            return (
                1.207
                - 1.848e-3 * T
                - 8.566e-7 * T**2
                + 4.001e-10 * T**3
                - 8.654e-14 * T**4
                + Kb * T * np.log(P / 1e5)
            )

        def mu_S_8(T, P):
            return (
                7.62e-1
                - 2.457e-3 * T
                - 4.012e-6 * T**2
                + 1.808e-9 * T**3
                - 3.810e-13 * T**4
                + Kb * T * np.log(P / 1e5)
            )

        def a_p(P):
            return 1.465e-02 - 2.115e-03 * np.log10(P) + 6.905e-04 * np.log10(P) ** 2

        b = 10
        c = 80
        w = 100

        mu_eV = (
            0.5 * (special.erfc((T - T_tr(P)) / w) * mu_S_8(T, P) / 8)
            + 0.5 * ((special.erf((T - T_tr(P)) / w) + 1) * mu_S_2(T, P) / 2)
            - a_p(P) * np.exp(-((T - T_tr(P) + b) ** 2) / (2 * c**2))
            + self.reference_energy
        )

        if units == "eV":
            return mu_eV

        elif units == "J":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
            )

        elif units == "kJ":
            return (
                mu_eV
                * constants.physical_constants["electron volt-joule relationship"][0]
                * constants.N_A
                * 0.001
            )

mu(T, P, units='eV', xc=None)

Returns the chemical potential of one atom of Sulfur

Args:

T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.

Note:

T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
Other T, P arrays will result in undefined behaviour.

Returns:

mu (float/ndarray): Chemical potential of one sulfur atom expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
Source code in thermopot/materials.py 
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def mu(self, T, P, units="eV", xc=None):
    """
    Returns the chemical potential of one atom of Sulfur

    Args:

        T (float/ndarray): 1D Numpy array containing temperature data (in Kelvin) as floats, or a single temperature as a float.
        P (float/ndarray): 2D Numpy array with a single row containing pressure data (in Pa) as floats, or a single pressure as a float.

    Note:

        T, P are orthogonal 2D arrays of length m and n, populated in one row/column: in this case mu is an m x n matrix.
        Other T, P arrays will result in undefined behaviour.

    Returns:

        mu (float/ndarray): Chemical potential of one sulfur atom expressed as floats in a m x n Numpy array where T, P are orthogonal 2D arrays of length m and n
    """

    Kb = scipy.constants.physical_constants["Boltzmann constant in eV/K"][0]

    if (
        np.any(T > 1500)
        or np.any(T < 400)
        or np.any(P > 10**7)
        or np.any(P < 10**0)
    ):
        print(
            """WARNING!: You are using the sulfur model beyond the temperature and/or pressure range it was fitted to.
             Region of validity is 400 - 1500 K, 10^0 - 10^7 Pa. """
        )

    def T_tr(P):
        return (
            5.077e2
            + 7.272e1 * np.log10(P)
            - 8.295 * np.log10(P) ** 2
            + 1.828 * np.log10(P) ** 3
        )

    def mu_S_2(T, P):
        return (
            1.207
            - 1.848e-3 * T
            - 8.566e-7 * T**2
            + 4.001e-10 * T**3
            - 8.654e-14 * T**4
            + Kb * T * np.log(P / 1e5)
        )

    def mu_S_8(T, P):
        return (
            7.62e-1
            - 2.457e-3 * T
            - 4.012e-6 * T**2
            + 1.808e-9 * T**3
            - 3.810e-13 * T**4
            + Kb * T * np.log(P / 1e5)
        )

    def a_p(P):
        return 1.465e-02 - 2.115e-03 * np.log10(P) + 6.905e-04 * np.log10(P) ** 2

    b = 10
    c = 80
    w = 100

    mu_eV = (
        0.5 * (special.erfc((T - T_tr(P)) / w) * mu_S_8(T, P) / 8)
        + 0.5 * ((special.erf((T - T_tr(P)) / w) + 1) * mu_S_2(T, P) / 2)
        - a_p(P) * np.exp(-((T - T_tr(P) + b) ** 2) / (2 * c**2))
        + self.reference_energy
    )

    if units == "eV":
        return mu_eV

    elif units == "J":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
        )

    elif units == "kJ":
        return (
            mu_eV
            * constants.physical_constants["electron volt-joule relationship"][0]
            * constants.N_A
            * 0.001
        )

Reaction

Class for reaction data

Sets properties:

reaction.reactants (Dict relating reactant materials to a number of formula units) reaction.products (Dict relating product materials to a number of formular units)

Sets methods:

reaction.DH(T,P) : Enthalpy of formation reaction.DU(T,P) : Internal energy change reaction.Dmu(T,P) : Gibbs free energy of formation

Source code in thermopot/reactions.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
class Reaction:
    """
    Class for reaction data

    Sets properties:
    -------------------
    reaction.reactants     (Dict relating reactant materials to a number of formula units)
    reaction.products      (Dict relating product materials to a number of formular units)

    Sets methods:
    -------------------
    reaction.DH(T,P) : Enthalpy of formation
    reaction.DU(T,P) : Internal energy change
    reaction.Dmu(T,P) : Gibbs free energy of formation
    """

    def __init__(
        self,
        reactants_dictionary,
        products_dictionary,
        temperature=298.15,
        pressure=1e5,
        fu=1,
    ):
        """
        reactants_dictionary and products dictionary takes the form { class_instance : formula units }
        and can have an arbitrary number of key-value pairs. `Class instance` is an instance of the `materials.solid`
        or `materials.ideal_gas` classes.

        temperature is provided in kelvin, pressure is provided in Pa.

        fu is the number of formula units of the final reactant(s). It is used to scale the calculated changes in energy/enthalpy.

        """
        self.reactants = reactants_dictionary
        self.products = products_dictionary
        self.T = temperature
        self.P = pressure
        self.fu_scaling = fu

    def DH(self, T=None, P=None, xc="pbesol", units="eV"):
        T = self.T if T is None else T
        P = self.P if P is None else P

        reactants_enthalpy, products_enthalpy = 0, 0
        for material, fu in self.reactants.items():
            reactants_enthalpy += material.H(T, P, xc=xc, units=units) * fu
        for material, fu in self.products.items():
            products_enthalpy += material.H(T, P, xc=xc, units=units) * fu

        return potential.Potential(
            (products_enthalpy - reactants_enthalpy) / self.fu_scaling, T, P
        )

    def DU(self, T=None, P=None, xc="pbesol", units="eV"):
        T = self.T if T is None else T
        P = self.P if P is None else P

        reactants_energy, products_energy = 0, 0
        for material, fu in self.reactants.items():
            reactants_energy += material.U(T, xc=xc, units=units) * fu
        for material, fu in self.products.items():
            products_energy += material.U(T, xc=xc, units=units) * fu

        return potential.Potential(
            (products_energy - reactants_energy) / self.fu_scaling, T, P
        )

    def Dmu(self, T=None, P=None, xc="pbesol", units="eV"):
        T = self.T if T is None else T
        P = self.P if P is None else P

        reactants_energy, products_energy = 0, 0
        for material, fu in self.reactants.items():
            reactants_energy += material.mu(T, P, xc=xc, units=units) * fu
        for material, fu in self.products.items():
            products_energy += material.mu(T, P, xc=xc, units=units) * fu

        return potential.Potential(
            (products_energy - reactants_energy) / self.fu_scaling, T, P
        )

    def DE(self, xc="pbesol"):
        # units here are whatever is specified when creating material...

        reactants_energy, products_energy = 0, 0

        for material, fu in self.reactants.items():
            reactants_energy += (material.energies[xc] / material.fu_cell) * fu
        for material, fu in self.products.items():
            products_energy += (material.energies[xc] / material.fu_cell) * fu

        return (products_energy - reactants_energy) / self.fu_scaling

__init__(reactants_dictionary, products_dictionary, temperature=298.15, pressure=100000.0, fu=1)

reactants_dictionary and products dictionary takes the form { class_instance : formula units } and can have an arbitrary number of key-value pairs. Class instance is an instance of the materials.solid or materials.ideal_gas classes.

temperature is provided in kelvin, pressure is provided in Pa.

fu is the number of formula units of the final reactant(s). It is used to scale the calculated changes in energy/enthalpy.

Source code in thermopot/reactions.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def __init__(
    self,
    reactants_dictionary,
    products_dictionary,
    temperature=298.15,
    pressure=1e5,
    fu=1,
):
    """
    reactants_dictionary and products dictionary takes the form { class_instance : formula units }
    and can have an arbitrary number of key-value pairs. `Class instance` is an instance of the `materials.solid`
    or `materials.ideal_gas` classes.

    temperature is provided in kelvin, pressure is provided in Pa.

    fu is the number of formula units of the final reactant(s). It is used to scale the calculated changes in energy/enthalpy.

    """
    self.reactants = reactants_dictionary
    self.products = products_dictionary
    self.T = temperature
    self.P = pressure
    self.fu_scaling = fu

Potential

Source code in thermopot/potential.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
class Potential:
    def __init__(self, potential, T, P):
        self.potential = potential
        self.T = T
        self.P = P

    def plot_TvsP(
        self,
        potential_label="$\Delta G_f$ / kJ mol$^{-1}$",
        scale_range=[-600, 0],
        precision="%d",
        T_units="K",
        P_units="Pa",
        log_scale=True,
    ):
        """
        T is an array e.g. np.linspace(100, 1500, 100)  # K
        P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose()  # Pa
        potential is returned from a reactions.reaction method called for an instance with attributes T,P.
        If T has length m and P has length n, P will be a 2D array with dimensions m x n.
        e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol()
        potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$'
        scale_range is the scale of the colorbar e.g. [-380, -240]
        logscale determines if the y-axis Pressure is logarithmic
        """

        mpl.rcParams["font.family"] = "serif"
        mpl.rcParams["font.serif"] = "Times New Roman"
        mpl.rcParams["font.size"] = 16

        # Unit conversions (all calculations are in SI units, conversion needed for plots)
        if T_units == "K":
            x_values = self.T
            x_unitlabel = "K"
        elif T_units == "C":
            x_values = self.T - 273.15
            x_unitlabel = "$^\circ$ C"
        else:
            raise ValueError("Invalid temperature unit: {0}".format(T_units))

        if P_units == "Pa":
            y_values = self.P.flatten()
        elif P_units == "Bar" or P_units == "bar":
            y_values = self.P.flatten() * 1e-5
        elif P_units == "mbar":
            y_values = self.P.flatten() * 1e-5 * 1e3
        elif P_units == "kPa":
            y_values = self.P.flatten() * 1e-3
        elif P_units == "mmHg" or P_units == "torr":
            y_values = self.P.flatten() * 760 / (1.01325e5)
        else:
            raise ValueError("Invalid pressure unit: {0}.".format(T_units))

        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        colormap = plt.get_cmap("summer")

        a = plt.contour(
            x_values, y_values, self.potential, 10, linewidths=1, colors="k"
        )
        plt.pcolormesh(
            x_values,
            y_values,
            self.potential,
            cmap=colormap,
            vmin=scale_range[0],
            vmax=scale_range[1],
            shading="auto",
        )
        colours = plt.colorbar()
        colours.set_label(potential_label, labelpad=20)

        if log_scale:
            ax.set_yscale("log")
        plt.clabel(a, fmt=precision)

        plt.xlabel("Temperature / {0}".format(x_unitlabel))
        plt.ylabel("Pressure / {0}".format(P_units))

        return plt

plot_TvsP(potential_label='$\\Delta G_f$ / kJ mol$^{-1}$', scale_range=[-600, 0], precision='%d', T_units='K', P_units='Pa', log_scale=True)

T is an array e.g. np.linspace(100, 1500, 100) # K P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose() # Pa potential is returned from a reactions.reaction method called for an instance with attributes T,P. If T has length m and P has length n, P will be a 2D array with dimensions m x n. e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol() potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$' scale_range is the scale of the colorbar e.g. [-380, -240] logscale determines if the y-axis Pressure is logarithmic

Source code in thermopot/potential.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def plot_TvsP(
    self,
    potential_label="$\Delta G_f$ / kJ mol$^{-1}$",
    scale_range=[-600, 0],
    precision="%d",
    T_units="K",
    P_units="Pa",
    log_scale=True,
):
    """
    T is an array e.g. np.linspace(100, 1500, 100)  # K
    P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose()  # Pa
    potential is returned from a reactions.reaction method called for an instance with attributes T,P.
    If T has length m and P has length n, P will be a 2D array with dimensions m x n.
    e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol()
    potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$'
    scale_range is the scale of the colorbar e.g. [-380, -240]
    logscale determines if the y-axis Pressure is logarithmic
    """

    mpl.rcParams["font.family"] = "serif"
    mpl.rcParams["font.serif"] = "Times New Roman"
    mpl.rcParams["font.size"] = 16

    # Unit conversions (all calculations are in SI units, conversion needed for plots)
    if T_units == "K":
        x_values = self.T
        x_unitlabel = "K"
    elif T_units == "C":
        x_values = self.T - 273.15
        x_unitlabel = "$^\circ$ C"
    else:
        raise ValueError("Invalid temperature unit: {0}".format(T_units))

    if P_units == "Pa":
        y_values = self.P.flatten()
    elif P_units == "Bar" or P_units == "bar":
        y_values = self.P.flatten() * 1e-5
    elif P_units == "mbar":
        y_values = self.P.flatten() * 1e-5 * 1e3
    elif P_units == "kPa":
        y_values = self.P.flatten() * 1e-3
    elif P_units == "mmHg" or P_units == "torr":
        y_values = self.P.flatten() * 760 / (1.01325e5)
    else:
        raise ValueError("Invalid pressure unit: {0}.".format(T_units))

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    colormap = plt.get_cmap("summer")

    a = plt.contour(
        x_values, y_values, self.potential, 10, linewidths=1, colors="k"
    )
    plt.pcolormesh(
        x_values,
        y_values,
        self.potential,
        cmap=colormap,
        vmin=scale_range[0],
        vmax=scale_range[1],
        shading="auto",
    )
    colours = plt.colorbar()
    colours.set_label(potential_label, labelpad=20)

    if log_scale:
        ax.set_yscale("log")
    plt.clabel(a, fmt=precision)

    plt.xlabel("Temperature / {0}".format(x_unitlabel))
    plt.ylabel("Pressure / {0}".format(P_units))

    return plt

Potentials

Source code in thermopot/potentials.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
class Potentials:
    def __init__(self, *potentials):
        self.potentials = potentials
        self.minimum_potential = self.find_potential_minimum()
        self.T = self.potentials[0].T
        self.P = self.potentials[0].P

    # TODO: check that all the potentials are plottable over T/P range

    def plot_TvsP(
        self,
        material_labels=None,
        T_units="K",
        P_units="Pa",
        log_scale=True,
    ):
        """
        T is an array e.g. np.linspace(100, 1500, 100)  # K
        P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose()  # Pa
        potential is returned from a reactions.reaction method called for an instance with attributes T,P.
        If T has length m and P has length n, P will be a 2D array with dimensions m x n.
        e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol()
        potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$'
        scale_range is the scale of the colorbar e.g. [-380, -240]
        log_scale determines if the Pressure y-axis is logarithmic
        """

        mpl.rcParams["font.family"] = "serif"
        mpl.rcParams["font.serif"] = "Times New Roman"
        mpl.rcParams["font.size"] = 16

        # Unit conversions (all calculations are in SI units, conversion needed for plots)

        if T_units == "K":
            x_values = self.T
            x_unitlabel = "K"
        elif T_units == "C":
            x_values = self.T - 273.15
            x_unitlabel = "$^\circ$ C"
        else:
            raise ValueError("Invalid temperature unit: {0}".format(T_units))

        if P_units == "Pa":
            y_values = self.P.flatten()
        elif P_units == "Bar" or P_units == "bar":
            y_values = self.P.flatten() * 1e-5
        elif P_units == "mbar":
            y_values = self.P.flatten() * 1e-5 * 1e3
        elif P_units == "kPa":
            y_values = self.P.flatten() * 1e-3
        elif P_units == "mmHg" or P_units == "torr":
            y_values = self.P.flatten() * 760 / (1.01325e5)
        else:
            raise ValueError("Invalid pressure unit: {0}.".format(T_units))

        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        colormap = plt.get_cmap("summer")

        potential = self.find_potential_minimum()
        plt.pcolormesh(
            x_values,
            y_values,
            potential / (len(material_labels) - 1),
            cmap=colormap,
            shading="auto",
        )
        # TODO: sort the colour map out so consistent with grid. Now ranges from 0 to 1

        # Set borders in the interval [0, 1]
        bound = np.linspace(0, 1, len(material_labels))

        # AT THE MOMENT THIS IS BROKEN!!!!
        # plt.legend(
        #    [mpatches.Patch(color=colormap(i)) for i in bound],
        #    ["{:s}".format(material_labels[i]) for i in range(len(material_labels))],
        # )

        plt.xlabel("Temperature / {0}".format(x_unitlabel))
        plt.ylabel("Pressure / {0}".format(P_units))
        if log_scale:
            ax.set_yscale("log")

        return plt

    def find_potential_minimum(self):
        assert (
            len(set([potential.potential.shape for potential in self.potentials])) == 1
        ), "potential arrays must have the same dimension"

        minimum_potential = self.potentials[0].potential
        for i, potential in enumerate(self.potentials):
            minimum_potential = np.minimum(
                minimum_potential, self.potentials[i + 1].potential
            )
            if i + 2 == len(self.potentials):
                break

        for i, potential in enumerate(self.potentials):
            minimum_potential[potential.potential == minimum_potential] = i

        return minimum_potential

plot_TvsP(material_labels=None, T_units='K', P_units='Pa', log_scale=True)

T is an array e.g. np.linspace(100, 1500, 100) # K P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose() # Pa potential is returned from a reactions.reaction method called for an instance with attributes T,P. If T has length m and P has length n, P will be a 2D array with dimensions m x n. e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol() potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$' scale_range is the scale of the colorbar e.g. [-380, -240] log_scale determines if the Pressure y-axis is logarithmic

Source code in thermopot/potentials.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def plot_TvsP(
    self,
    material_labels=None,
    T_units="K",
    P_units="Pa",
    log_scale=True,
):
    """
    T is an array e.g. np.linspace(100, 1500, 100)  # K
    P is an array orthogonal to T. e.g. np.array(np.logspace(1, 7, 100), ndmin=2).transpose()  # Pa
    potential is returned from a reactions.reaction method called for an instance with attributes T,P.
    If T has length m and P has length n, P will be a 2D array with dimensions m x n.
    e.g. reactions.reaction({Ba:1,S:2}, {BaS2:1}},temperature=T,pressure=P).Dmu_eV_pbesol()
    potential_label is the label of the contour colorbar e.g. '$\Delta G_f$ / kJ mol$^{-1}$'
    scale_range is the scale of the colorbar e.g. [-380, -240]
    log_scale determines if the Pressure y-axis is logarithmic
    """

    mpl.rcParams["font.family"] = "serif"
    mpl.rcParams["font.serif"] = "Times New Roman"
    mpl.rcParams["font.size"] = 16

    # Unit conversions (all calculations are in SI units, conversion needed for plots)

    if T_units == "K":
        x_values = self.T
        x_unitlabel = "K"
    elif T_units == "C":
        x_values = self.T - 273.15
        x_unitlabel = "$^\circ$ C"
    else:
        raise ValueError("Invalid temperature unit: {0}".format(T_units))

    if P_units == "Pa":
        y_values = self.P.flatten()
    elif P_units == "Bar" or P_units == "bar":
        y_values = self.P.flatten() * 1e-5
    elif P_units == "mbar":
        y_values = self.P.flatten() * 1e-5 * 1e3
    elif P_units == "kPa":
        y_values = self.P.flatten() * 1e-3
    elif P_units == "mmHg" or P_units == "torr":
        y_values = self.P.flatten() * 760 / (1.01325e5)
    else:
        raise ValueError("Invalid pressure unit: {0}.".format(T_units))

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    colormap = plt.get_cmap("summer")

    potential = self.find_potential_minimum()
    plt.pcolormesh(
        x_values,
        y_values,
        potential / (len(material_labels) - 1),
        cmap=colormap,
        shading="auto",
    )
    # TODO: sort the colour map out so consistent with grid. Now ranges from 0 to 1

    # Set borders in the interval [0, 1]
    bound = np.linspace(0, 1, len(material_labels))

    # AT THE MOMENT THIS IS BROKEN!!!!
    # plt.legend(
    #    [mpatches.Patch(color=colormap(i)) for i in bound],
    #    ["{:s}".format(material_labels[i]) for i in range(len(material_labels))],
    # )

    plt.xlabel("Temperature / {0}".format(x_unitlabel))
    plt.ylabel("Pressure / {0}".format(P_units))
    if log_scale:
        ax.set_yscale("log")

    return plt

get_potential_aims(file, property)

Thermodynamic property interpolation function. Requires phonopy-FHI-aims output file. Reads data for S and Cv expressed in J/K/mol, F and U in kJ/mol, TS in J/mol. Outputs data for S and Cv in kB/cell, U, F and TS in eV/cell.

Source code in thermopot/interpolate.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def get_potential_aims(file, property):
    """Thermodynamic property interpolation function. Requires phonopy-FHI-aims output file.
    Reads data for S and Cv expressed in J/K/mol, F and U in kJ/mol,
    TS in J/mol.
    Outputs data for S and Cv in kB/cell, U, F and TS in eV/cell.
    """
    data = genfromtxt(file)
    T = data[:, 0]
    if property in ("Cv", "Cp", "heat_capacity", "C"):
        potential = data[:, 3] / kB2JKmol
    elif property in ("U", "internal_energy"):
        potential = data[:, 4] / eV2kJmol
    elif property in ("F", "A", "Helmholtz", "free_energy"):
        potential = data[:, 1] / eV2kJmol
    elif property in ("S", "Entropy", "entropy"):
        potential = data[:, 2] / kB2JKmol
    elif property in ("TS"):
        potential = (T * data[:, 2]) / eV2Jmol
    else:
        raise RuntimeError("Property not found")
    thefunction = interp1d(T, potential, kind="linear")

    return thefunction

get_potential_nist_table(file, property)

Thermodynamic property interpolation function. Requires NIST-JANAF table. All properties in J, mol and K

Source code in thermopot/interpolate.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def get_potential_nist_table(file, property):
    """Thermodynamic property interpolation function. Requires NIST-JANAF table. All properties in J, mol and K"""
    data = genfromtxt(file, skip_header=2)
    T = data[:, 0]
    if property in ("Cp", "C", "heat_capacity"):
        potential = data[:, 1]
    elif property in ("S", "entropy"):
        potential = data[:, 2]
    elif property in ("H", "enthalpy"):
        potential = (data[:, 4] - data[0, 4]) * 1e3
    elif property in ("U", "internal_energy"):
        # U = H - PV; for ideal gas molar PV = RT so U = H - RT
        from scipy.constants import R as R

        potential = (data[:, 4] - data[0, 4]) * 1e3 - R * data[:, 0]
    elif property in ("DH", "Delta_H", "standard_enthalpy_change"):
        potential = data[:, 4] * 1e3
    else:
        raise RuntimeError("Property not found")

    thefunction = interp1d(T, potential, kind="cubic")
    return thefunction

get_potential_sulfur_table(filename)

Read thermodynamic property as function of T, P from datafile.

Datafile should be generated by the code at http://github.com/WMD-bath/sulfur-model or follow the same format

Source code in thermopot/interpolate.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def get_potential_sulfur_table(filename):
    """
    Read thermodynamic property as function of T, P from datafile.

    Datafile should be generated by the code at http://github.com/WMD-bath/sulfur-model
    or follow the same format

    """
    # Import chemical potential in J mol-1 vs T, P from file
    data = genfromtxt(filename, comments="#", delimiter=",")
    T = data[:, 0].flatten()
    with open(filename, "r") as f:
        header = f.readline()
    P = [float(p) for p in re.findall(r"\d+.\d+", header)]
    thefunction = interp2d(T, np.log(P), data[:, 1:].transpose(), kind="cubic")

    def lin_P_function(T, P):
        return thefunction(T, np.log(P))

    return lin_P_function