Update phonons.py

src/atomate2/common/schemas/phonons.py

@@ -29,6 +29,17 @@
 from pymatgen.symmetry.kpath import KPathSeek
 from typing_extensions import Self
 
+# imported lib by jiongzhi zheng 
+from ase.io import read
+from hiphive import ClusterSpace, ForceConstantPotential, ForceConstants, enforce_rotational_sum_rules
+from hiphive.cutoffs import estimate_maximum_cutoff
+from hiphive.utilities import extract_parameters
+import subprocess
+from phonopy.file_IO import write_FORCE_CONSTANTS, parse_FORCE_CONSTANTS
+from phonopy.interface.vasp import write_vasp
+from phonopy.interface.vasp import read_vasp
+
+
 from atomate2.aims.utils.units import omegaToTHz
 
 logger = logging.getLogger(__name__)
@@ -297,9 +308,115 @@ def from_forces_born(
             symprec=symprec,
             is_symmetry=sym_reduce,
         )
+        #phonon.generate_displacements(distance=displacement)
+        #set_of_forces = [np.array(forces) for forces in displacement_data["forces"]]
+
+
+        # modified by jiongzhi zheng start from here 
+        supercell = phonon.get_supercell()
+        write_vasp("POSCAR", cell)
+        write_vasp("SPOSCAR", supercell)
+
         phonon.generate_displacements(distance=displacement)
-        set_of_forces = [np.array(forces) for forces in displacement_data["forces"]]
+        disps_j = phonon.displacements
+
+        f_disp_n1 = len(disps_j)
+
+        # jiongzhi zheng
+        # Print keys and their corresponding values
+        for key, value in displacement_data.items():
+            print(f"{key}: {value}")
+
+
+        if f_disp_n1 < 1:
+            set_of_forces = [np.array(forces) for forces in displacement_data["forces"]]
+            set_of_forces_a = np.array(set_of_forces)
+            set_of_disps = [np.array(disps.cart_coords) for disps in displacement_data["displaced_structures"]]
+            set_of_disps_m = np.array(set_of_disps)
+            print(set_of_disps_m)
+
+            print(np.shape(set_of_disps_m))
+            n_shape = set_of_disps_m.shape[0]
+            # set_of_disps_d = {}
+            # set_of_disps_d['displacements'] = set_of_disps_m
+            set_of_disps_d = {'displacements': set_of_disps_m, 'dtype': 'double', 'order': 'C'}
+        else:
+            set_of_forces = [np.array(forces) for forces in displacement_data["forces"]]
+            set_of_forces_a_o = np.array(set_of_forces)
+            set_of_forces_a_t = set_of_forces_a_o - set_of_forces_a_o[-1, :, :]
+            set_of_forces_a = set_of_forces_a_t[:-1, :, :]
+            set_of_disps = [np.array(disps.cart_coords) for disps in displacement_data["displaced_structures"]]
+            set_of_disps_m_o = np.round((set_of_disps - supercell.get_positions()), decimals=16).astype('double')
+            set_of_disps_m = set_of_disps_m_o[:-1, :, :]
+            print(set_of_disps_m)
+
+            print(np.shape(set_of_disps_m))
+            n_shape = set_of_disps_m.shape[0]
+            # set_of_disps_d = {}
+            # set_of_disps_d['displacements'] = set_of_disps_m
+            set_of_disps_d = {'displacements': set_of_disps_m, 'dtype': 'double', 'order': 'C'}
+            # phonon.set_displacement_dataset(set_of_disps_d)
+            # print (set_of_disps)
+
+        import pickle
+        with open("disp_matrix.pkl","wb") as file: 
+             pickle.dump(set_of_disps_m,file)
+        with open("force_matrix.pkl","wb") as file:
+             pickle.dump(set_of_forces_a,file)
+
+        """"put a if else condition here to determine whether we need to calculate the higher order force constants """
+        Calc_anharmonic_FCs = False
+        if Calc_anharmonic_FCs:
+            from alm import ALM
+            supercell_z = phonon.supercell
+            lattice = supercell_z.cell
+            positions = supercell_z.scaled_positions
+            numbers = supercell_z.numbers
+            natom = len(numbers)
+            with ALM(lattice, positions, numbers) as alm:
+                alm.define(1)
+                alm.suggest()
+                n_fp = alm._get_number_of_irred_fc_elements(1)
+
+            num = int(np.ceil(n_fp / (3.0 * natom)))
+            num_round = int(np.round((n_fp / (3.0 * natom))))
+
+            if num > num_round:
+                num_d = num
+                displacement_t = 0.01
+                phonon.generate_displacements(displacement_t)
+                num_disp_t = len(phonon.displacements)
+                int_num = int(num_disp_t / num_d)
+                if int_num > 3:
+                    num_d = int(np.ceil(int_num / 3.0))
+                else:
+                    num_d = int(np.ceil(int_num / 3.0) + 1)
+            else:
+                num_d = num
+                displacement_t = 0.01
+                phonon.generate_displacements(displacement_t)
+                num_disp_t = len(phonon.displacements)
+                int_num = int(num_disp_t / num_d)
+                if int_num >= 3:
+                    num_d = int(np.ceil(int_num / 3.0))
+                else:
+                    num_d = int(num + 1)
+
+            displacement_f = 0.01
+            phonon.generate_displacements(distance=displacement_f)
+            disps = phonon.displacements
+
+            f_disp_n = int(len(disps))
+            if f_disp_n > 2:
+                num_har = num_d
+            else:
+                num_har = f_disp_n
+        else:
+            num_har = n_shape
 
+
+        #set_of_forces = [np.array(forces) for forces in displacement_data["forces"]]
+
         if born is not None and epsilon_static is not None:
             if len(structure) == len(born):
                 borns, epsilon = symmetrize_borns_and_epsilon(
@@ -325,9 +442,47 @@ def from_forces_born(
         else:
             borns = None
             epsilon = None
+
+        # Produces all force constants
+        #phonon.produce_force_constants(forces=set_of_forces)
+        #phonon.produce_force_constants(forces=set_of_forces,  fc_calculator="alm") 
+        #phonon.symmetrize_force_constants()
+        #fcs = phonon.force_constants
+        #write_FORCE_CONSTANTS(fcs)
+        prim = read('POSCAR')
+        supercell = read('SPOSCAR')
+
+
+        pheasy_cmd_1 = 'pheasy --dim "{0}" "{1}" "{2}" -s -w 2 --symprec 1e-3 --nbody 2'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]))
+        pheasy_cmd_2 = 'pheasy --dim "{0}" "{1}" "{2}" -c --symprec 1e-3 -w 2'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]))
+        pheasy_cmd_3 = 'pheasy --dim "{0}" "{1}" "{2}" -w 2 -d --symprec 1e-3 --ndata "{3}" --disp_file'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+
+        displacement_f = 0.01
+        phonon.generate_displacements(distance=displacement_f)
+        disps = phonon.displacements
+        num_judge = len(disps)
+
+        if num_judge > 3:
+           pheasy_cmd_4 = 'pheasy --dim "{0}" "{1}" "{2}" -f --full_ifc -w 2 --symprec 1e-3 -l LASSO --std --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+        else:
+            pheasy_cmd_4 = 'pheasy --dim "{0}" "{1}" "{2}" -f --full_ifc -w 2 --symprec 1e-3 --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+            #pheasy_cmd_4_1 = 'pheasy --dim "{0}" "{1}" "{2}" -f --full_ifc -w 2 --hdf5 --symprec 1e-3 --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+
+        print ("Start running pheasy in andes8: Dartmouth College Clusters or NERSC Perlmutter")
+        subprocess.call(pheasy_cmd_1, shell=True)
+        subprocess.call(pheasy_cmd_2, shell=True)
+        subprocess.call(pheasy_cmd_3, shell=True)
+        subprocess.call(pheasy_cmd_4, shell=True)
+        #subprocess.call(pheasy_cmd_4_1, shell=True)
+
+
 
         # Produces all force constants
-        phonon.produce_force_constants(forces=set_of_forces)
+        # phonon.produce_force_constants(forces=set_of_forces)
+
+        force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS")
+        phonon.force_constants = force_constants
+        phonon.symmetrize_force_constants()
 
         # with phonopy.load("phonopy.yaml") the phonopy API can be used
         phonon.save("phonopy.yaml")
@@ -369,6 +524,109 @@ def from_forces_born(
             tol=kwargs.get("tol_imaginary_modes", 1e-5)
         )
 
+         # jiongzhi zheng
+        if imaginary_modes:
+            # Define a cluster space using the largest cutoff you can
+            max_cutoff = estimate_maximum_cutoff(supercell) - 0.01
+            cutoffs = [max_cutoff]  # only second order needed
+            cs = ClusterSpace(prim, cutoffs)
+
+            # import the phonopy force constants using the correct supercell also provided by phonopy
+            fcs = ForceConstants.read_phonopy(supercell, 'FORCE_CONSTANTS')
+            # Find the parameters that best fits the force constants given you cluster space
+            parameters = extract_parameters(fcs, cs)
+            # Enforce the rotational sum rules
+            parameters_rot = enforce_rotational_sum_rules(cs, parameters, ['Huang','Born-Huang'], alpha=1e-6)
+            # use the new parameters to make a fcp and then create the force constants and write to a phonopy file
+
+            fcp = ForceConstantPotential(cs, parameters_rot)
+            fcs = fcp.get_force_constants(supercell)
+            fcs.write_to_phonopy('FORCE_CONSTANTS_new', format='text')
+
+            force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS_new")
+            phonon.force_constants = force_constants
+            phonon.symmetrize_force_constants()
+
+            phonon.run_band_structure(qpoints, path_connections=connections, with_eigenvectors=True)
+            phonon.write_yaml_band_structure(filename=filename_band_yaml)
+            bs_symm_line = get_ph_bs_symm_line(filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None)
+
+            new_plotter = PhononBSPlotter(bs=bs_symm_line)
+
+            new_plotter.save_plot("phonon_band_structure.eps",img_format=kwargs.get("img_format", "eps"),units=kwargs.get("units", "THz"),)
+
+            imaginary_modes_hiphive = bs_symm_line.has_imaginary_freq(
+            tol=kwargs.get("tol_imaginary_modes", 1e-5)
+        )
+
+        else:
+           imaginary_modes_hiphive = False
+
+        if imaginary_modes_hiphive:
+            pheasy_cmd_11 = 'pheasy --dim "{0}" "{1}" "{2}" -s -w 2 --c2 10.0 --symprec 1e-3 --nbody 2'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]))
+            pheasy_cmd_12 = 'pheasy --dim "{0}" "{1}" "{2}" -c --symprec 1e-3 --c2 10.0 -w 2'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]))
+            pheasy_cmd_13 = 'pheasy --dim "{0}" "{1}" "{2}" -w 2 -d --symprec 1e-3 --c2 10.0 --ndata "{3}" --disp_file'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+
+            displacement_f = 0.01
+            phonon.generate_displacements(distance=displacement_f)
+            disps = phonon.displacements
+            num_judge = len(disps)
+
+            if num_judge > 3:
+                pheasy_cmd_14 = 'pheasy --dim "{0}" "{1}" "{2}" -f --c2 10.0 --full_ifc -w 2 --symprec 1e-3 -l LASSO --std --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+            else:
+                pheasy_cmd_14 = 'pheasy --dim "{0}" "{1}" "{2}" -f --full_ifc --c2 10.0 -w 2 --symprec 1e-3 --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+            #pheasy_cmd_4_1 = 'pheasy --dim "{0}" "{1}" "{2}" -f --full_ifc -w 2 --hdf5 --symprec 1e-3 --rasr BHH --ndata "{3}"'.format(int(supercell_matrix[0][0]), int(supercell_matrix[1][1]), int(supercell_matrix[2][2]), int(num_har))
+
+            print ("Start running pheasy in discovery")
+            subprocess.call(pheasy_cmd_11, shell=True)
+            subprocess.call(pheasy_cmd_12, shell=True)
+            subprocess.call(pheasy_cmd_13, shell=True)
+            subprocess.call(pheasy_cmd_14, shell=True)
+
+            force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS")
+            phonon.force_constants = force_constants
+            phonon.symmetrize_force_constants()
+
+
+            #phonon.force_constants = fcs
+
+            # with phonon.load("phonopy.yaml") the phonopy API can be used
+            phonon.save("phonopy.yaml")
+
+            # get phonon band structure
+            kpath_dict, kpath_concrete = cls.get_kpath(
+                structure=get_pmg_structure(phonon.primitive),
+                kpath_scheme=kpath_scheme,
+                symprec=symprec,
+            )
+
+            npoints_band = kwargs.get("npoints_band", 101)
+            qpoints, connections = get_band_qpoints_and_path_connections(
+                kpath_concrete, npoints=kwargs.get("npoints_band", 101)
+            )
+
+            # phonon band structures will always be cmouted
+            filename_band_yaml = "phonon_band_structure.yaml"
+            phonon.run_band_structure(
+                qpoints, path_connections=connections, with_eigenvectors=True
+            )
+            phonon.write_yaml_band_structure(filename=filename_band_yaml)
+            bs_symm_line = get_ph_bs_symm_line(
+                filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None
+            )
+            new_plotter = PhononBSPlotter(bs=bs_symm_line)
+
+            new_plotter.save_plot(
+                "phonon_band_structure.eps",
+                img_format=kwargs.get("img_format", "eps"),
+                units=kwargs.get("units", "THz"),
+            )
+
+        else:
+            pass
+
+
         # gets data for visualization on website - yaml is also enough
         if kwargs.get("band_structure_eigenvectors"):
             bs_symm_line.write_phononwebsite("phonon_website.json")