Add vibrational properties

README.md

@@ -30,8 +30,8 @@ Suggestions, improvements, and edits are most welcome.
 The following snippet can be used to test the code (provided that you have some VASP data or intermediate `.mat` files.
 
 ```python
-import GEMDAT
-GEMDAT.analyse_md('<path to data>', 'Li', 'argyrodite')
+import gemdat
+gemdat.analyse_md('<path to data>', 'Li', 'argyrodite')
 
 ```
 

docs/api.md

@@ -1 +1 @@
-::: GEMDAT
+::: gemdat

pyproject.toml

@@ -65,7 +65,7 @@ publishing = [
 ]
 
 [tool.setuptools]
-package-dir = {"GEMDAT" = "src", "pymatgen.analysis.GEMDAT" = "src"}
+package-dir = {"gemdat" = "src", "pymatgen.analysis.gemdat" = "src"}
 include-package-data = true
 
 [tool.coverage.run]

src/__init__.py

@@ -4,10 +4,15 @@
     plot_displacement_per_element,
     plot_displacement_per_site,
 )
+from .project import load_project
+from .vibration import plot_frequency_vs_occurence, plot_vibrational_amplitudes
 
 __all__ = [
     'calculate_displacements',
     'plot_displacement_per_site',
     'plot_displacement_per_element',
     'plot_displacement_histogram',
+    'plot_frequency_vs_occurence',
+    'plot_vibrational_amplitudes',
+    'load_project',
 ]

src/constants.py

@@ -0,0 +1,3 @@
+e_charge = 1.60217657e-19
+k_boltzmann = 1.3806488e-23
+avogadro = 6.022140857e23

src/displacements.py

@@ -59,19 +59,25 @@ def calculate_lengths(vectors: np.ndarray,
     return np.sqrt(total_displacement)
 
 
-def calculate_displacements(traj_coords: np.ndarray) -> np.ndarray:
+def calculate_displacements(traj_coords: np.ndarray,
+                            lattice: Lattice,
+                            equilibration_steps: int = 0) -> np.ndarray:
     """Calculate displacements from first set of positions.
 
     Corrects for elements jumping to the next unit cell.
 
     Parameters
     ----------
-    traj_coords : np.array[i, j, k]
+    traj_coords : np.ndarray
         3-dimensional numpy array with dimensions i: time_steps, j: sites, k: coordinates
+    lattice : Lattice
+        Lattice parameters
+    equilibration_steps : int, optional
+        Number of steps to skip before equilibration
 
     Returns
     -------
-    displacements : np.ndarray[i, j]
+    displacements : np.ndarray
         Displacements from first set of positions.
     """
     offsets = calculate_cell_offsets_from_coords(traj_coords)
@@ -90,7 +96,7 @@ def calculate_displacements(traj_coords: np.ndarray) -> np.ndarray:
 
     displacements = np.array(displacements)
 
-    return displacements
+    return displacements.T
 
 
 def plot_displacement_per_site(displacements: np.ndarray):
@@ -103,7 +109,7 @@ def plot_displacement_per_site(displacements: np.ndarray):
     """
     fig, ax = plt.subplots()
 
-    for site_displacement in displacements.T:
+    for site_displacement in displacements:
         ax.plot(site_displacement, lw=0.3)
 
     ax.set(title='Displacement of diffusing element',
@@ -128,7 +134,7 @@ def plot_displacement_per_element(species: list[Element],
 
     grouped = defaultdict(list)
 
-    for specie, displacement in zip(species, displacements.T):
+    for specie, displacement in zip(species, displacements):
         grouped[specie.name].append(displacement)
 
     fig, ax = plt.subplots()
@@ -154,69 +160,32 @@ def plot_displacement_histogram(displacements: np.ndarray):
         Numpy array with displacements
     """
     fig, ax = plt.subplots()
-    ax.hist(displacements[-1])
+    ax.hist(displacements[:, -1])
     ax.set(title='Histogram of displacement of diffusing element',
            xlabel='Displacement (Angstrom)',
            ylabel='Nr. of atoms')
     plt.show()
 
 
 if __name__ == '__main__':
-    from pathlib import Path
-
-    import yaml
-    from pymatgen.io.vasp.outputs import Vasprun
+    from gemdat import load_project
 
     vasp_xml = '/run/media/stef/Scratch/md-analysis-matlab-example/vasprun.xml'
 
-    path_coords = Path('traj_coords.npy')
-    path_data = Path('data.yaml')
-
     # skip first timesteps
     equilibration_steps = 1250
 
-    if not (path_coords.exists() and path_data.exists()):
-        vasprun = Vasprun(
-            vasp_xml,
-            parse_dos=False,
-            parse_eigen=False,
-            parse_projected_eigen=False,
-            parse_potcar_file=False,
-        )
-
-        traj = vasprun.get_trajectory()
-
-        structure = vasprun.structures[0]
-
-        lattice = structure.lattice
-        species = structure.species
-
-        data = {
-            'species': [e.as_dict() for e in species],
-            'lattice': lattice.as_dict()
-        }
-
-        with open('data.yaml', 'w') as f:
-            yaml.dump(data, f)
-
-        traj.to_positions()
-        np.save('traj_coords.npy', traj.coords)
-
-        traj_coords = traj.coords
-    else:
-        traj_coords = np.load('traj_coords.npy')
-
-        with open('data.yaml') as f:
-            data = yaml.safe_load(f)
+    traj_coords, data = load_project(vasp_xml, diffusing_element='Li')
 
-        species = [Element.from_dict(e) for e in data['species']]
-        lattice = Lattice.from_dict(data['lattice'])
+    species = data['species']
+    lattice = data['lattice']
 
     print(species)
     print(lattice)
     print(traj_coords.shape)
 
-    displacements = calculate_displacements(traj_coords)
+    displacements = calculate_displacements(
+        traj_coords, lattice=lattice, equilibration_steps=equilibration_steps)
 
     plot_displacement_per_site(displacements)
 

src/project.py

@@ -0,0 +1,70 @@
+from pathlib import Path
+from typing import Optional
+
+import numpy as np
+import yaml
+from pymatgen.core import Element, Lattice
+from pymatgen.io.vasp.outputs import Vasprun
+
+PATH_COORDS = Path('traj_coords.npy')
+PATH_DATA = Path('data.yaml')
+VASP_XML = 'vasprun.xml'
+
+
+def load_project(directory: str | Path,
+                 diffusing_element: Optional[str] = None):
+    """Read coordinates and data from vasp_xml."""
+    vasp_xml = Path(directory) / VASP_XML
+
+    if not (PATH_COORDS.exists() and PATH_DATA.exists()):
+        assert vasp_xml.exists()
+
+        vasprun = Vasprun(
+            vasp_xml,
+            parse_dos=False,
+            parse_eigen=False,
+            parse_projected_eigen=False,
+            parse_potcar_file=False,
+        )
+
+        traj = vasprun.get_trajectory()
+
+        structure = vasprun.structures[0]
+
+        data = {
+            'species': [e.as_dict() for e in structure.species],
+            'lattice':
+            structure.lattice.as_dict(),
+            'parameters':
+            vasprun.parameters,
+            # Size of the timestep (*1E-15 = in femtoseconds)
+            'time_step':
+            vasprun.time_step['POTIM'] * 1e-15,
+            # Temperature of the MD simulation
+            'temperature':
+            vasprun.parameters['TEBEG'],
+            # Number of diffusing elements
+            'nr_diffusing':
+            sum([
+                e.name == diffusing_element
+                for e in vasprun.initial_structure.species
+            ])
+        }
+
+        with open('data.yaml', 'w') as f:
+            yaml.dump(data, f)
+
+        traj.to_positions()
+        np.save('traj_coords.npy', traj.coords)
+
+        traj_coords = traj.coords
+    else:
+        traj_coords = np.load('traj_coords.npy')
+
+        with open('data.yaml') as f:
+            data = yaml.unsafe_load(f)
+
+    data['species'] = [Element.from_dict(e) for e in data['species']]
+    data['lattice'] = Lattice.from_dict(data['lattice'])
+
+    return traj_coords, data