From d3fffaa3290ce756d87bda83633c69d3a7da16fe Mon Sep 17 00:00:00 2001 From: Jiongzhi ZHENG Date: Sun, 30 Jun 2024 14:54:04 -0400 Subject: [PATCH] Update phonons.py --- src/atomate2/common/schemas/phonons.py | 262 ++++++++++++++++++++++++- 1 file changed, 260 insertions(+), 2 deletions(-) diff --git a/src/atomate2/common/schemas/phonons.py b/src/atomate2/common/schemas/phonons.py index a594afa623..8de61e0695 100644 --- a/src/atomate2/common/schemas/phonons.py +++ b/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")