Extract maximum project force constant - max_pfc

📝 Note

Authors like to thank Dr. Jan Hempelmann for sharing the initial code from his private repository available at https://github.com/JanHempelmann/projectionscript

The read_force_constants function has been reused as is from the above mentioned git repository

The get_max_pfc function uses parts of generate_projected_force_constants function from the above mentioned git repository

import numpy as np
import pandas as pd
import multiprocessing as mp
from tqdm.autonotebook import tqdm
from pymatgen.core import Structure
import warnings

warnings.filterwarnings("ignore")

def read_force_constants(file_name: str = "FORCE_CONSTANTS"):
    """
    Read FORCE_CONSTANTS file.
    """
    data: dict[str, list] = {"n_atoms": [], "atomic_pairs": [], "force_matrices": []}
    mtrx = []
    with open(file_name) as f:
        for i, lin in enumerate(f):
            line_split = lin.strip().split()
            if i == 0:
                data["n_atoms"] = [int(x) for x in line_split]
            elif len(line_split) == 2:
                data["atomic_pairs"].append([int(x) for x in line_split])
            else:
                mtrx.append([float(x) for x in line_split])
                if len(mtrx) == 3:
                    data["force_matrices"].append(mtrx)
                    mtrx = []
    return data

Please set the default value of phonopy_outputs_path to path where you stored the Phonopy_fc after conversion of ddbs file

For demonstration here, only one entry is included in the path below

# Directory with phonopy files (we need only FORCECONSTANTS, primitive and supercell POSCARS)
phonopy_outputs_path = "example_phonon_db_files/phonopy_fc/"

def get_max_pfc(mpid, phonopy_outputs_path=phonopy_outputs_path, cutoff=6):
    """
    Calculate projected force constant and return max
    """
    
    force_constants = read_force_constants(f"{phonopy_outputs_path}/{mpid}/FORCE_CONSTANTS")
    phonopy_primitive = Structure.from_file(f"{phonopy_outputs_path}/{mpid}/POSCAR")
    phonopy_sc = Structure.from_file(f"{phonopy_outputs_path}/{mpid}/supercell_POSCAR")

    # create a dict with atom-pairs and associated force constant matrices
    pair_fc_dict = { "-".join([str(i) for i in pair]) : force_constants["force_matrices"][ix] 
            for ix, pair in enumerate(force_constants["atomic_pairs"])}

    # get list of neighbours for primitive and supercell cell based on cutoff
    neigh_sc = phonopy_sc.get_neighbor_list(r=cutoff)
    neigh_primitive = phonopy_primitive.get_neighbor_list(r=cutoff)

    mean_pfcs = []
    
    for _, (src, dst, img) in enumerate(zip(neigh_sc[0], neigh_sc[1], neigh_sc[2])):
        site_a = phonopy_sc[src].frac_coords
        site_b = phonopy_sc[dst].frac_coords + img
    
        bv_frac = site_b - site_a
    
        bv_cart = np.dot(bv_frac, phonopy_sc.lattice.matrix)

        
        bl_here = np.linalg.norm(bv_cart)
    
        unit_vector = bv_cart /  bl_here
    
        if pair_fc_dict.get(f"{src+1}-{dst+1}"):

            # project force constant matrix onto unit vector
            pfc = np.matmul(pair_fc_dict.get(f"{src+1}-{dst+1}"), unit_vector)
            pfc_transpose = np.transpose(pfc)
            norm_pfc = np.linalg.norm(pfc)
            norm_pfc_transpose = np.linalg.norm(pfc_transpose)
        
            mean_pfc = np.mean([norm_pfc, norm_pfc_transpose])
        
            mean_pfcs.append(mean_pfc)

    if len(mean_pfcs) != neigh_primitive[-1].shape[0]:
        warnings.warn(f"{mpid} : Different number of interactions for primitive and supercell under same cutoff of {cutoff} Angs")

    return {mpid: max(mean_pfcs)}

The code below will show for one example entry using the data obtained after conversion of the source DDB file in example_phonon_db_files/phonopy_fc/

result = get_max_pfc(mpid="mp-66")

Use the code snippet below once you have FORCECONSTANTS, primitive and supercell POSCARS extracted from the ddb files using the convert_ddb_to_phonopy.ipynb script. All the ddb source files can be downloaded from here

# load the mpids for which lobster calcs are available in our LOBSTER database and are part of matbench benchmark dataset for last phonon dos peak. PFC data for particular mpids will be extracted
with open("pfc_mpids.txt", "r", encoding="utf-8") as f:
    mpids = [line.rstrip("\n") for line in f]

all_data = {}
with (
    mp.Pool(processes=8, maxtasksperchild=1) as pool,
    tqdm(total=len(mpids), desc="Extracting max_pfc") as pbar,
    ):
        for _, result in enumerate(pool.imap_unordered(get_max_pfc, mpids, chunksize=1)):
            pbar.update()
            all_data.update(result)

# Save the complete target dataset df
pd.DataFrame.from_dict(data=all_data, orient="index", columns=["max_pfc"]).to_json("max_pfc.json")