From 369881a4df0f4e77429962790c537616744595be Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 4 Mar 2024 08:11:19 +0000 Subject: [PATCH 01/14] Bump pymatgen from 2024.2.23 to 2024.3.1 (#759) Bumps [pymatgen](https://github.com/materialsproject/pymatgen) from 2024.2.23 to 2024.3.1. - [Release notes](https://github.com/materialsproject/pymatgen/releases) - [Changelog](https://github.com/materialsproject/pymatgen/blob/master/docs/CHANGES.md) - [Commits](https://github.com/materialsproject/pymatgen/compare/v2024.2.23...v2024.3.1) --- updated-dependencies: - dependency-name: pymatgen dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index d9f9e0accb..596e7c50bb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -88,7 +88,7 @@ strict = [ "pydantic-settings==2.2.1", "pydantic==2.6.2", "pymatgen-analysis-defects==2024.2.24", - "pymatgen==2024.2.23", + "pymatgen==2024.3.1", "quippy-ase==0.9.14", "seekpath==2.1.0", "typing-extensions==4.10.0", From 27bde10a36df3a908ea8bdabaa6dd323c0050632 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 4 Mar 2024 08:30:11 +0000 Subject: [PATCH 02/14] Bump pydantic from 2.6.2 to 2.6.3 (#763) Bumps [pydantic](https://github.com/pydantic/pydantic) from 2.6.2 to 2.6.3. - [Release notes](https://github.com/pydantic/pydantic/releases) - [Changelog](https://github.com/pydantic/pydantic/blob/main/HISTORY.md) - [Commits](https://github.com/pydantic/pydantic/compare/v2.6.2...v2.6.3) --- updated-dependencies: - dependency-name: pydantic dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 596e7c50bb..bf42b10544 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -86,7 +86,7 @@ strict = [ "numpy", "phonopy==2.21.2", "pydantic-settings==2.2.1", - "pydantic==2.6.2", + "pydantic==2.6.3", "pymatgen-analysis-defects==2024.2.24", "pymatgen==2024.3.1", "quippy-ase==0.9.14", From 8c2072496d6b763a34c1edb4205a90f64aedff6d Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Mon, 4 Mar 2024 15:12:49 -0800 Subject: [PATCH 03/14] commenting out quality control temporarily commenting out quality control job --- src/atomate2/common/flows/hiphive.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index ff4184f98e..629d7297da 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -15,7 +15,7 @@ from atomate2.common.jobs.hiphive import ( hiphive_static_calcs, - quality_control, + # quality_control, run_fc_to_pdos, run_hiphive, run_hiphive_renormalization, From a0775604bee3b77d5c108ddc0d9d1ec923d534c2 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 8 Mar 2024 00:01:42 -0800 Subject: [PATCH 04/14] adding renormalization changes from Junsoo --- src/atomate2/common/flows/hiphive.py | 2 +- src/atomate2/common/jobs/hiphive.py | 592 +++++++++++++-------------- 2 files changed, 288 insertions(+), 306 deletions(-) diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index 629d7297da..c58421c5a6 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -131,7 +131,7 @@ class BaseHiphiveMaker(Maker, ABC): i * 100 for i in range(21) ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) T_RENORM: ClassVar[list[int]] = [ - 400 # 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000 + 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000 # 400 ] # [i*100 for i in range(0,16)] # Temp. at which renorm. is to be performed # Temperature at which lattice thermal conductivity is calculated # If renormalization is performed, diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index 0faa8d4097..aaacfcfe73 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -83,9 +83,9 @@ IMAGINARY_TOL = 0.025 # in THz FIT_METHOD = "rfe" -ev2j = sp.constants.elementary_charge -hbar = sp.constants.hbar # J-s -kb = sp.constants.Boltzmann # J/K +eV2J = sp.constants.elementary_charge +hbar = sp.constants.hbar # J-s +kB = sp.constants.Boltzmann # J/K __all__ = [ "hiphive_static_calcs", @@ -947,14 +947,13 @@ def harmonic_properties( structure: Structure, supercell_matrix: np.ndarray, fcs: ForceConstants, - temperature: list, - imaginary_tol: float | None = IMAGINARY_TOL, -) -> tuple[dict, Phonopy]: + temperature: List, + imaginary_tol: float = IMAGINARY_TOL +) -> Tuple[Dict,Phonopy]: """ - Calculate harmonic properties. - Uses the force constants to extract phonon properties. Used for comparing the accuracy of force constant fits. + Args: structure: The parent structure. supercell_matrix: The supercell transformation matrix. @@ -962,242 +961,229 @@ def harmonic_properties( imaginary_tol: Tolerance used to decide if a phonon mode is imaginary, in THz. - Returns - ------- + Returns: A tuple of the number of imaginary modes at Gamma, the minimum phonon - frequency at Gamma, and the free energy, entropy, and heat capacity. + frequency at Gamma, and the free energy, entropy, and heat capacity """ - logger.info("Evaluating harmonic properties...") + + logger.info('Evaluating harmonic properties...') fcs2 = fcs.get_fc_array(2) - fcs.get_fc_array(3) + fcs3 = fcs.get_fc_array(3) parent_phonopy = get_phonopy_structure(structure) phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) natom = phonopy.primitive.get_number_of_atoms() - mesh = supercell_matrix.diagonal() * 2 + mesh = supercell_matrix.diagonal()*2 + logger.info(f'Mesh: {mesh}') phonopy.set_force_constants(fcs2) - phonopy.set_mesh( - mesh, is_eigenvectors=True, is_mesh_symmetry=False - ) # run_mesh(is_gamma_center=True) + phonopy.set_mesh(mesh,is_eigenvectors=True,is_mesh_symmetry=False) #run_mesh(is_gamma_center=True) phonopy.run_thermal_properties(temperatures=temperature) - logger.info("Thermal properties successfully run!") + logger.info('Thermal properties successfully run!') _, free_energy, entropy, heat_capacity = phonopy.get_thermal_properties() - free_energy *= 1000 / sp.constants.Avogadro / ev2j / natom # kJ/mol to eV/atom - entropy *= 1 / sp.constants.Avogadro / ev2j / natom # J/K/mol to eV/K/atom - heat_capacity *= 1 / sp.constants.Avogadro / ev2j / natom # J/K/mol to eV/K/atom - - freq = phonopy.mesh.frequencies # in THz + free_energy *= 1000/sp.constants.Avogadro/eV2J/natom # kJ/mol to eV/atom + entropy *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom + heat_capacity *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom + logger.info(f"Heat_capacity_harmonic_property: {heat_capacity}") + freq = phonopy.mesh.frequencies # in THz + logger.info(f'Frequencies: {freq}') # find imaginary modes at gamma - # phonopy.run_qpoints([0, 0, 0]) - # gamma_eigs = phonopy.get_qpoints_dict()["frequencies"] +# phonopy.run_qpoints([0, 0, 0]) +# gamma_eigs = phonopy.get_qpoints_dict()["frequencies"] n_imaginary = int(np.sum(freq < -np.abs(imaginary_tol))) - np.min(freq) + min_freq = np.min(freq) if n_imaginary == 0: - logger.info("No imaginary modes!") - else: # do not calculate these if imaginary modes exist - logger.warning("Imaginary modes found!") + logger.info('No imaginary modes!') + else: # do not calculate these if imaginary modes exist + logger.warning('Imaginary modes found!') - # if len(temperature) == 1: + # if len(temperature)==1: # temperature = temperature[0] # free_energy = free_energy[0] # entropy = entropy[0] # heat_capacity = heat_capacity[0] - + logger.info(f"Heat_capacity_harmonic_property[0]: {heat_capacity}") + temperature = temperature[0] return { "temperature": temperature, "free_energy": free_energy, "entropy": entropy, "heat_capacity": heat_capacity, - "n_imaginary": n_imaginary, - }, phonopy + "n_imaginary": n_imaginary + }, phonopy def anharmonic_properties( phonopy: Phonopy, fcs: ForceConstants, - temperature: list, + temperature: List, heat_capacity: np.ndarray, n_imaginary: float, - bulk_modulus: float | None = None, -) -> dict: + bulk_modulus: float = None +) -> Dict: + if n_imaginary == 0: - logger.info("Evaluating anharmonic properties...") + logger.info('Evaluating anharmonic properties...') fcs2 = fcs.get_fc_array(2) fcs3 = fcs.get_fc_array(3) - grun, cte, dLfrac = gruneisen( - phonopy, fcs2, fcs3, temperature, heat_capacity, bulk_modulus=bulk_modulus - ) - else: # do not calculate these if imaginary modes exist - logger.warning( - "Gruneisen and thermal expansion cannot be calculated with " - "imaginary modes. All set to 0." - ) - grun = [np.zeros((len(temperature), 3))] - cte = [np.zeros((len(temperature), 3))] - dLfrac = np.zeros((len(temperature), 3)) + grun, cte, dLfrac = gruneisen(phonopy,fcs2,fcs3,temperature,heat_capacity,bulk_modulus=bulk_modulus) + else: # do not calculate these if imaginary modes exist + logger.warning('Gruneisen and thermal expansion cannot be calculated with imaginary modes. All set to 0.') + grun = np.zeros((len(temperature),3)) + cte = np.zeros((len(temperature),3)) + dLfrac = np.zeros((len(temperature),3)) return { "gruneisen": grun, "thermal_expansion": cte, "expansion_fraction": dLfrac, - } + } def get_total_grun( - omega: np.ndarray, grun: np.ndarray, kweight: np.ndarray, T: float + omega: np.ndarray, + grun: np.ndarray, + kweight: np.ndarray, + T: float ) -> np.ndarray: - # total = 0 - total = np.zeros((3, 3)) + total = 0 weight = 0 nptk = omega.shape[0] nbands = omega.shape[1] - omega = abs(omega) * 1e12 * 2 * np.pi - if T == 0: - total = np.zeros((3, 3)) + omega = abs(omega)*1e12*2*np.pi + if T==0: + total = np.zeros((3,3)) grun_total_diag = np.zeros(3) else: for i in range(nptk): for j in range(nbands): - x = hbar * omega[i, j] / (2.0 * kb * T) - dBE = (x / np.sinh(x)) ** 2 - weight += dBE * kweight[i] - total += dBE * kweight[i] * grun[i, j] - total = total / weight - grun_total_diag = np.array([total[0, 2], total[1, 1], total[2, 0]]) - - def percent_diff(a: float, - b: float) -> float: - return abs((a - b) / b) - + x = hbar*omega[i,j]/(2.0*kB*T) + dBE = (x/np.sinh(x))**2 + weight += dBE*kweight[i] + total += dBE*kweight[i]*grun[i,j] + total = total/weight + grun_total_diag = np.array([total[0,2],total[1,1],total[2,0]]) + + def percent_diff(a,b): + return abs((a-b)/b) # This process preserves cell symmetry upon thermal expansion, i.e., it prevents # symmetry-identical directions from inadvertently expanding by different ratios - # when the Gruneisen routine returns slightly different ratios for - # those directions - avg012 = np.mean((grun_total_diag[0], grun_total_diag[1], grun_total_diag[2])) - avg01 = np.mean((grun_total_diag[0], grun_total_diag[1])) - avg02 = np.mean((grun_total_diag[0], grun_total_diag[2])) - avg12 = np.mean((grun_total_diag[1], grun_total_diag[2])) - if percent_diff(grun_total_diag[0], avg012) < 0.1: - if percent_diff(grun_total_diag[1], avg012) < 0.1: - if percent_diff(grun_total_diag[2], avg012) < 0.1: # all siilar + # when/if the Gruneisen routine returns slightly different ratios for those directions + avg012 = np.mean((grun_total_diag[0],grun_total_diag[1],grun_total_diag[2])) + avg01 = np.mean((grun_total_diag[0],grun_total_diag[1])) + avg02 = np.mean((grun_total_diag[0],grun_total_diag[2])) + avg12 = np.mean((grun_total_diag[1],grun_total_diag[2])) + if percent_diff(grun_total_diag[0],avg012) < 0.1: + if percent_diff(grun_total_diag[1],avg012) < 0.1: + if percent_diff(grun_total_diag[2],avg012) < 0.1: # all siilar grun_total_diag[0] = avg012 grun_total_diag[1] = avg012 grun_total_diag[2] = avg012 - elif percent_diff(grun_total_diag[2], avg02) < 0.1: # 0 and 2 similar + elif percent_diff(grun_total_diag[2],avg02) < 0.1: # 0 and 2 similar grun_total_diag[0] = avg02 grun_total_diag[2] = avg02 - elif percent_diff(grun_total_diag[2], avg12) < 0.1: # 1 and 2 similar + elif percent_diff(grun_total_diag[2],avg12) < 0.1: # 1 and 2 similar grun_total_diag[1] = avg12 grun_total_diag[2] = avg12 else: pass - elif percent_diff(grun_total_diag[1], avg01) < 0.1: # 0 and 1 similar + elif percent_diff(grun_total_diag[1],avg01) < 0.1: # 0 and 1 similar grun_total_diag[0] = avg01 grun_total_diag[1] = avg01 - elif percent_diff(grun_total_diag[1], avg12) < 0.1: # 1 and 2 similar + elif percent_diff(grun_total_diag[1],avg12) < 0.1: # 1 and 2 similar grun_total_diag[1] = avg12 grun_total_diag[2] = avg12 else: pass - elif percent_diff(grun_total_diag[0], avg01) < 0.1: # 0 and 1 similar + elif percent_diff(grun_total_diag[0],avg01) < 0.1: # 0 and 1 similar grun_total_diag[0] = avg01 grun_total_diag[1] = avg01 - elif percent_diff(grun_total_diag[0], avg02) < 0.1: # 0 and 2 similar + elif percent_diff(grun_total_diag[0],avg02) < 0.1: # 0 and 2 similar grun_total_diag[0] = avg02 grun_total_diag[2] = avg02 - else: # nothing similar + else: # nothing similar pass return grun_total_diag def gruneisen( - phonopy: Phonopy, - fcs2: np.ndarray, - fcs3: np.ndarray, - temperature: list, - heat_capacity: np.ndarray, # in eV/K/atom - bulk_modulus: float | None = None, # in GPa -) -> tuple[list, list, np.ndarray]: - gruneisen = Gruneisen(fcs2, fcs3, phonopy.supercell, phonopy.primitive) - gruneisen.set_sampling_mesh(phonopy.mesh_numbers, is_gamma_center=True) + phonopy: Phonopy, + fcs2: np.ndarray, + fcs3: np.ndarray, + temperature: List, + heat_capacity: np.ndarray, # in eV/K/atom + bulk_modulus: float = None # in GPa +) -> Tuple[List,List]: + + gruneisen = Gruneisen(fcs2,fcs3,phonopy.supercell,phonopy.primitive) + gruneisen.set_sampling_mesh(phonopy.mesh_numbers,is_gamma_center=True) gruneisen.run() - grun = gruneisen.get_gruneisen_parameters() # (nptk,nmode,3,3) + grun = gruneisen.get_gruneisen_parameters() # (nptk,nmode,3,3) omega = gruneisen._frequencies + qp = gruneisen._qpoints kweight = gruneisen._weights - grun_tot = [] - # for temp in temperature: - # grun_tot.append(get_total_grun(omega, grun, kweight, temp)) - # grun_tot = np.nan_to_num(np.array(grun_tot)) - - grun_tot = [ - np.nan_to_num(np.array(get_total_grun(omega, grun, kweight, temp))) - for temp in temperature - ] + grun_tot = list() + for temp in temperature: + grun_tot.append(get_total_grun(omega,grun,kweight,temp)) + grun_tot = np.nan_to_num(np.array(grun_tot)) + + # grun_tot = [ + # np.nan_to_num(np.array(get_total_grun(omega, grun, kweight, temp))) + # for temp in temperature + # ] - # linear thermal expansion coefficient and fraction + # linear thermal expansion coefficeint and fraction if bulk_modulus is None: cte = None dLfrac = None else: - heat_capacity *= ( - ev2j * phonopy.primitive.get_number_of_atoms() - ) # eV/K/atom to J/K + heat_capacity *= eV2J*phonopy.primitive.get_number_of_atoms() # eV/K/atom to J/K + # Convert heat_capacity to an array if it's a scalar + # heat_capacity = np.array([heat_capacity]) + logger.info(f"heat capacity = {heat_capacity}") vol = phonopy.primitive.get_volume() - # cte = grun_tot*heat_capacity.repeat(3)/(vol/10**30)/(bulk_modulus*10**9)/3 - # Check the shapes of the arrays + logger.info(f"grun_tot: {grun_tot}") + logger.info(f"grun_tot shape: {grun_tot.shape}") logger.info(f"heat_capacity shape: {heat_capacity.shape}") logger.info(f"heat_capacity: {heat_capacity}") logger.info(f"vol: {vol}") logger.info(f"bulk_modulus: {bulk_modulus}") - - # cte = np.array(grun_tot)*heat_capacity.repeat(3).reshape(len(heat_capacity),3)/(vol/10**30)/(bulk_modulus[0]*10**9)/3 - # cte = grun_tot*heat_capacity.repeat(3).reshape(len(heat_capacity),3)/(vol/10**30)/(bulk_modulus[0]*10**9)/3 +# cte = grun_tot*heat_capacity.repeat(3)/(vol/10**30)/(bulk_modulus*10**9)/3 cte = grun_tot*heat_capacity.repeat(3).reshape(len(heat_capacity),3)/(vol/10**30)/(bulk_modulus*10**9)/3 - # cte = ( - # grun_tot - # * heat_capacity.repeat(3).reshape(len(heat_capacity), 3) - # / (vol / 10**30) - # / (bulk_modulus * 10**9) - # / 3 - # ) + # cte = grun_tot*heat_capacity.repeat(3).reshape(1,3)/(vol/10**30)/(bulk_modulus*10**9)/3 cte = np.nan_to_num(cte) - dLfrac = thermal_expansion(temperature, cte) - if len(temperature) == 1: - dLfrac = dLfrac[-1] - logger.info(f"Gruneisen: \n {grun_tot}") - logger.info(f"Coefficient of Thermal Expansion: \n {cte}") - logger.info(f"Linear Expansion Fraction: \n {dLfrac}") - + if len(temperature)==1: + dLfrac = cte*temperature + else: + dLfrac = thermal_expansion(temperature, cte) + logger.info('Gruneisen: \n {}'.format(grun_tot)) + logger.info('Coefficient of Thermal Expansion: \n {}'.format(cte)) + logger.info('Linear Expansion Fraction: \n {}'.format(dLfrac)) + return grun_tot, cte, dLfrac def thermal_expansion( - temperature: list, - cte: np.ndarray, + temperature: list, + cte: np.array, ) -> np.ndarray: - assert len(temperature) == len(cte) + assert len(temperature)==len(cte) if 0 not in temperature: - temperature = [0, *temperature] - cte = np.array([np.array([0, 0, 0]), *list(cte)]) - temperature_np = np.array(temperature) - # ind = np.argsort(temperature_np) - ind = np.array(np.argsort(temperature_np)) - logger.info(f"ind = {ind}") - logger.info(f"temperature_np = {temperature_np}") - logger.info(f"cte = {cte}") - # temperature_np = temperature[ind] - temperature_np = temperature - # cte = np.array(cte)[ind] - cte = np.array(cte) + temperature = [0] + temperature + cte = np.array([np.array([0,0,0])] + list(cte)) + temperature = np.array(temperature) + ind = np.argsort(temperature) + temperature = temperature[ind] + cte = np.array(cte)[ind] # linear expansion fraction dLfrac = copy(cte) - for t in range(len(temperature_np)): - dLfrac[t, :] = np.trapz(cte[: t + 1, :], temperature_np[: t + 1], axis=0) - return np.nan_to_num(dLfrac) + for t in range(len(temperature)): + dLfrac[t,:] = np.trapz(cte[:t+1,:],temperature[:t+1],axis=0) + dLfrac = np.nan_to_num(dLfrac) + return dLfrac @job def run_thermal_cond_solver( @@ -1364,7 +1350,7 @@ def run_fc_to_pdos( logger.info("Finished inserting force constants and phonon data") else: - renorm_thermal_data = loadfn("renorm_thermal_data.json") + renorm_thermal_data = loadfn("thermal_data.json") # renorm_thermal_data.json fcs = ForceConstants.read("force_constants.fcs") T = renorm_thermal_data["temperature"] @@ -1444,126 +1430,92 @@ def run_hiphive_renormalization( structure_data = loadfn(f"structure_data_{loop}.json") phonopy_orig = phpy.load("phonopy_params.yaml") - thermal_data = loadfn("thermal_data.json") - thermal_data = thermal_data["heat_capacity"] + # thermal_data = loadfn("thermal_data.json") + # thermal_data = thermal_data["heat_capacity"] cutoffs = fitting_data["cutoffs"] fit_method = fitting_data["fit_method"] - fitting_data["n_imaginary"] - fitting_data["imaginary_tol"] + + + # fitting_data["n_imaginary"] + # fitting_data["imaginary_tol"] parent_structure = structure_data["structure"] supercell_structure = structure_data["supercell_structure"] - supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) supercell_matrix = np.array(structure_data["supercell_matrix"]) + parent_atoms = AseAtomsAdaptor.get_atoms(parent_structure) + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + + # Renormalization with DFT lattice - renorm_data = run_renormalization( - parent_structure, - supercell_atoms, - supercell_matrix, - cs, - fcs, - param, - temperature, - nconfig, - max_iter, - conv_thresh, - renorm_method, - fit_method, - bulk_modulus, - phonopy_orig, - ) + TD_data = run_renormalization(parent_structure, supercell_structure, supercell_matrix, + cs, fcs, param, temperature, nconfig, renorm_method, + fit_method, bulk_modulus, phonopy_orig) + TD_structure_data = copy(structure_data) + TD_structure_data["structure"] = parent_structure + TD_structure_data["supercell_structure"] = supercell_structure # Additional renormalization with thermal expansion - # optional - just single "iteration" for now if renorm_TE_iter: n_TE_iter = 1 for i in range(n_TE_iter): - if renorm_data is None or renorm_data["n_imaginary"] < 0: + if TD_data is None or TD_data["n_imaginary"] < 0: # failed, incomplete, or still imaginary break logger.info( f"Renormalizing with thermally expanded lattice - iteration {i}" ) - dLfrac = renorm_data["expansion_fraction"] - param = renorm_data["param"] - - parent_structure_TE, supercell_atoms_TE, cs_TE, fcs_TE = setup_TE_iter( - cs, - cutoffs, - parent_structure, - param, - temperature, - dLfrac, - supercell_matrix, - ) - prim_TE_atoms = AseAtomsAdaptor.get_atoms(parent_structure_TE) - prim_TE_phonopy = PhonopyAtoms( - symbols=prim_TE_atoms.get_chemical_symbols(), - scaled_positions=prim_TE_atoms.get_scaled_positions(), - cell=prim_TE_atoms.cell, - ) - phonopy_TE = Phonopy( - prim_TE_phonopy, - supercell_matrix=supercell_matrix, - primitive_matrix=None, - ) - - renorm_data = run_renormalization( - parent_structure_TE, - supercell_atoms_TE, - supercell_matrix, - cs_TE, - fcs_TE, - param, - temperature, - nconfig, - max_iter, - conv_thresh, - renorm_method, - fit_method, - bulk_modulus, - phonopy_TE, - ) - structure_data["structure"] = parent_structure_TE - structure_data["supercell_structure"] = AseAtomsAdaptor.get_structure( - supercell_atoms_TE - ) - + dLfrac = TD_data["expansion_fraction"] + param_TD = TD_data["param"] + + a, b, c, d, e, failed = setup_TE_renorm( + cs,cutoffs,parent_atoms,supercell_atoms,param_TD,dLfrac,supercell_matrix + ) + if not failed: + parent_structure_TD, supercell_structure_TD, cs_TD, phonopy_TD, fcs_TD = a, b, c, d, e + TD_data = run_renormalization(parent_structure_TD, supercell_structure_TD, supercell_matrix, + cs_TD, fcs, param, temperature, nconfig, + renorm_method, fit_method, bulk_modulus, + phonopy_TD, param_TD, fcs_TD + ) + TD_structure_data["structure"] = parent_structure_TD + TD_structure_data["supercell_structure"] = supercell_structure_TD + + # Thermodynamic integration for anharmonic free energy + TD_data = thermodynamic_integration_ifc(TD_data, # everything TD + fcs, # original + param, # original + ) # write results logger.info("Writing renormalized results") - # renorm_thermal_data = {} - renorm_thermal_data: dict[str, Any] = {} - fcs = renorm_data["fcs"] - fcs.write("force_constants.fcs") - thermal_keys = [ - "temperature", - "free_energy", - "entropy", - "heat_capacity", - "gruneisen", - "thermal_expansion", - "expansion_fraction", - "free_energy_correction_S", - "free_energy_correction_SC", - ] - renorm_thermal_data = {key: [] for key in thermal_keys} + fcs_TD = TD_data['fcs'] + fcs_TD.write("force_constants.fcs") + if "n_imaginary" in TD_data: + thermal_keys = ["temperature","free_energy","entropy","heat_capacity", + "free_energy_correction_S","free_energy_correction_SC", + "free_energy_correction_TI"] + else: + thermal_keys = ["temperature","free_energy","entropy","heat_capacity", + "gruneisen","thermal_expansion","expansion_fraction", + "free_energy_correction_S","free_energy_correction_SC", + "free_energy_correction_TI"] + TD_thermal_data = {key: [] for key in thermal_keys} for key in thermal_keys: - renorm_thermal_data[key].append(renorm_data[key]) + TD_thermal_data[key].append(TD_data[key]) - logger.info(f"DEBUG: {renorm_data}") - if renorm_data["n_imaginary"] > 0: - logger.warning(f"Imaginary modes remain for {temperature} K!") - logger.warning("ShengBTE files not written") - logger.warning("No renormalization with thermal expansion") + logger.info("DEBUG: ",TD_data) + if TD_data["n_imaginary"] > 0: + logger.warning('Imaginary modes remain still exist') + logger.warning('ShengBTE FORCE_CONSTANTS_2ND not written') else: - logger.info("No imaginary modes! Writing ShengBTE files") - fcs.write_to_phonopy("FORCE_CONSTANTS_2ND".format(), format="text") + logger.info("No imaginary modes! Writing ShengBTE FORCE_CONSTANTS_2ND...") + fcs_TD.write_to_phonopy("FORCE_CONSTANTS_2ND".format(temperature), format="text") - dumpfn(structure_data, "structure_data.json".format()) - dumpfn(renorm_thermal_data, "renorm_thermal_data.json".format()) + dumpfn(TD_structure_data, "structure_data.json") + dumpfn(TD_thermal_data, "thermal_data.json") current_dir = os.getcwd() @@ -1572,19 +1524,19 @@ def run_hiphive_renormalization( def run_renormalization( structure: Structure, - supercell: Atoms, + supercell_structure: Structure, supercell_matrix: np.ndarray, cs: ClusterSpace, fcs: ForceConstants, param: np.ndarray, T: float, nconfig: int, - max_iter: int, - conv_tresh: float, renorm_method: str, fit_method: str, bulk_modulus: float = None, phonopy_orig: Phonopy = None, + param_TD: np.ndarray = None, + fcs_TD: ForceConstants = None, imaginary_tol: float = IMAGINARY_TOL, ) -> dict: """ @@ -1608,78 +1560,108 @@ def run_renormalization( frequency at Gamma, and the free energy, entropy, and heat capacity. """ nconfig = int(nconfig) - renorm = Renormalization(cs, supercell, fcs, param, T, renorm_method, fit_method) - fcp, fcs, param = renorm.renormalize(nconfig) # ,conv_tresh) - - renorm_data, phonopy = harmonic_properties( - structure, supercell_matrix, fcs, [T], imaginary_tol - ) + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + renorm = Renormalization(cs,supercell_atoms,param,fcs,T,renorm_method,fit_method,param_TD=param_TD,fcs_TD=fcs_TD) + fcp_TD, fcs_TD, param_TD = renorm.renormalize(nconfig)#,conv_tresh) - if renorm_data["n_imaginary"] == 0: - logger.info(f"Renormalized phonon is completely real at T = {T} K!") - anharmonic_data = anharmonic_properties( - phonopy, - fcs, - [T], - renorm_data["heat_capacity"], - renorm_data["n_imaginary"], - bulk_modulus=bulk_modulus, + TD_data, phonopy_TD = harmonic_properties( + structure, supercell_matrix, fcs_TD, [T], imaginary_tol ) - # else: - # anharmonic_data = dict() - # anharmonic_data["temperature"] = T - # anharmonic_data["gruneisen"] = np.array([0,0,0]) - # anharmonic_data["thermal_expansion"] = np.array([0,0,0]) - # anharmonic_data["expansion_fraction"] = np.array([0,0,0]) - renorm_data.update(anharmonic_data) - - phonopy_orig.run_mesh() - phonopy.run_mesh() - omega0 = phonopy_orig.mesh.frequencies # THz - omega_TD = phonopy.mesh.frequencies # THz - evec = phonopy.mesh.eigenvectors - # natom = phonopy.primitive.get_number_of_atoms() - correction_S, correction_SC = free_energy_correction( - omega0, omega_TD, evec, [T] - ) # eV/atom - - renorm_data["free_energy_correction_S"] = correction_S[0] - renorm_data["free_energy_correction_SC"] = correction_SC[0] - renorm_data["fcp"] = fcp - renorm_data["fcs"] = fcs - renorm_data["param"] = param - - return renorm_data - -def setup_TE_iter( - cs, cutoffs, parent_structure, param, temperatures, dLfrac, supercell_matrix -): - - new_atoms = AseAtomsAdaptor.get_atoms(parent_structure) - new_cell = Cell( - np.transpose( - [new_atoms.get_cell()[:, i] * (1 + dLfrac[0, i]) for i in range(3)] + logger.info(f"Heat capacity_TD_DATA: {TD_data['heat_capacity']}") + if TD_data["n_imaginary"] == 0: + logger.info(f'Renormalized phonon is completely real at T = {T} K!') + anharmonic_data = anharmonic_properties( + phonopy_TD, fcs_TD, [T], TD_data["heat_capacity"], TD_data["n_imaginary"], bulk_modulus=bulk_modulus ) - ) - new_atoms.set_cell(new_cell, scale_atoms=True) - parent_structure_TE = AseAtomsAdaptor.get_structure(new_atoms) - supercell_atoms_TE = AseAtomsAdaptor.get_atoms( - parent_structure_TE * supercell_matrix - ) - new_cutoffs = [i * (1 + np.linalg.norm(dLfrac)) for i in cutoffs] - fcs_TE = [] # not sure if this is correct - + TD_data.update(anharmonic_data) + + # phonopy_orig.run_mesh() + # phonopy_TD.run_mesh() + mesh = supercell_matrix.diagonal()*2 + phonopy_orig.set_mesh(mesh,is_eigenvectors=True,is_mesh_symmetry=False) + phonopy_TD.set_mesh(mesh,is_eigenvectors=True,is_mesh_symmetry=False) + omega_h = phonopy_orig.mesh.frequencies # THz + evec_h = phonopy_orig.mesh.eigenvectors + omega_TD = phonopy_TD.mesh.frequencies # THz + evec_TD = phonopy_TD.mesh.eigenvectors + logger.info(f'TD_data = {TD_data}') + logger.info(f'omega_h = {omega_h}') + logger.info(f'omega_TD = {omega_TD}') + logger.info(f'shape of omega_h = {omega_h.shape}') + logger.info(f'shape of omega_TD = {omega_TD.shape}') + logger.info(f'evec_h = {evec_h}') + logger.info(f'evec_TD = {evec_TD}') + logger.info(f"phonopy_orig.mesh = {phonopy_orig.mesh}") + logger.info(f"phonopy_TD.mesh = {phonopy_TD.mesh}") + correction_S, correction_SC = free_energy_correction(omega_h,omega_TD,evec_h,evec_TD,[T]) # eV/atom + + TD_data["supercell_structure"] = supercell_structure + TD_data["free_energy_correction_S"] = correction_S # S = -(dF/dT)_V quasiparticle correction + TD_data["free_energy_correction_SC"] = correction_SC # SCPH 4th-order correction (perturbation theory) + TD_data["fcp"] = fcp_TD + TD_data["fcs"] = fcs_TD + TD_data["param"] = param_TD + TD_data['cs'] = cs + + return TD_data + +def thermodynamic_integration_ifc( + TD_data: dict, + fcs: ForceConstants, + param: np.ndarray, + lambda_array: np.ndarray = np.array([0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]), + TI_nconfig=3, +) -> dict: + supercell_structure = TD_data["supercell_structure"] + cs = TD_data['cs'] + fcs_TD = TD_data["fcs"] + param_TD = TD_data["param"] + T = TD_data['temperature'] + logger.info(f"Temperature = {T}") + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + renorm = Renormalization(cs, supercell_atoms, param, fcs, T, 'least_squares', 'rfe', param_TD, fcs_TD) + matcov_TD, matcov_BO, matcov_TDBO = renorm.born_oppenheimer_qcv(TI_nconfig) + correction_TI = renorm.thermodynamic_integration(lambda_array, matcov_TD, matcov_BO, matcov_TDBO, TI_nconfig) + TD_data["free_energy_correction_TI"] = correction_TI + return TD_data + +def setup_TE_renorm(cs,cutoffs,parent_atoms,supercell_atoms,param,dLfrac,supercell_matrix): + parent_atoms_TE = copy(parent_atoms) + new_cell = Cell(np.transpose([parent_atoms_TE.get_cell()[:,i]*(1+dLfrac[0,i]) for i in range(3)])) + parent_atoms_TE.set_cell(new_cell,scale_atoms=True) + parent_structure_TE = AseAtomsAdaptor.get_structure(parent_atoms_TE) + supercell_atoms_TE = copy(supercell_atoms) + new_supercell = Cell(np.transpose([supercell_atoms_TE.get_cell()[:,i]*(1+dLfrac[0,i]) for i in range(3)])) + supercell_atoms_TE.set_cell(new_supercell,scale_atoms=True) + supercell_structure_TE = AseAtomsAdaptor.get_structure(supercell_atoms_TE) + count = 0 + failed = False + cs_TE = ClusterSpace(parent_atoms_TE,cutoffs,symprec=1e-2,acoustic_sum_rules=True) while True: - cs_TE = ClusterSpace(atoms, new_cutoffs, 1e-3, acoustic_sum_rules=True) + count += 1 if cs_TE.n_dofs == cs.n_dofs: break - if cs_TE.n_dofs > cs.n_dofs: - new_cutoffs = [i * 0.999 for i in new_cutoffs] - if cs_TE.n_dofs < cs.n_dofs: - new_cutoffs = [i * 1.001 for i in new_cutoffs] - new_fcp = ForceConstantPotential(cs_TE, param) - fcs_TE.append(new_fcp.get_force_constants(supercell_atoms_TE)) - return parent_structure_TE, supercell_atoms_TE, cs_TE, fcs_TE + elif count>10: + logger.warning("Could not find ClusterSpace for expanded cell identical to the original cluster space!") + failed = True + break + elif count==1: + cutoffs_TE = [i*(1+np.linalg.norm(dLfrac)) for i in cutoffs] + elif cs_TE.n_dofs > cs.n_dofs: + cutoffs_TE = [i*0.999 for i in cutoffs_TE] + elif cs_TE.n_dofs < cs.n_dofs: + cutoffs_TE = [i*1.001 for i in cutoffs_TE] + cs_TE = ClusterSpace(parent_atoms_TE,cutoffs_TE,symprec=1e-2,acoustic_sum_rules=True) + if failed: + return None, None, None, None, None, failed + else: + fcp_TE = ForceConstantPotential(cs_TE, param) + fcs_TE = fcp_TE.get_force_constants(supercell_atoms_TE) + prim_TE_phonopy = PhonopyAtoms(symbols=parent_atoms_TE.get_chemical_symbols(), + scaled_positions=parent_atoms_TE.get_scaled_positions(), + cell=parent_atoms_TE.cell) + phonopy_TE = Phonopy(prim_TE_phonopy, supercell_matrix=supercell_matrix, primitive_matrix=None) + return parent_structure_TE, supercell_structure_TE, cs_TE, phonopy_TE, fcs_TE, failed @job From 11c7ff16b8fe5e91556f91363dcd11d9b600cfcb Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 13:45:01 -0700 Subject: [PATCH 05/14] fixes 1. using replace instead of add in hiphive_static_calcs job 2. ensuring that the given structure. in primitive, if not then converting to prim structure. 3. updated logic to write the FC_3rd order file 4. ruff fixes 5. added support of amabte, shengbte, and phono3py in LTC calculation 6. updated logic to only use lambda=0 for TI, when imag_modes_bool = False --- src/atomate2/common/jobs/hiphive.py | 114 +++++++++++++++++----------- 1 file changed, 70 insertions(+), 44 deletions(-) diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index aaacfcfe73..75095400c8 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -48,6 +48,7 @@ # Phonopy & Phono3py from phonopy import Phonopy +from phonopy.interface.hiphive_interface import phonopy_atoms_to_ase from phonopy.structure.atoms import PhonopyAtoms from pymatgen.io.ase import AseAtomsAdaptor from pymatgen.io.phonopy import ( @@ -57,6 +58,7 @@ get_phonopy_structure, ) from pymatgen.io.shengbte import Control +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer from pymatgen.transformations.standard_transformations import SupercellTransformation from atomate2.common.jobs.phonons import get_supercell_size, run_phonon_displacements @@ -78,7 +80,8 @@ # Temperature at which lattice thermal conductivity is calculated # If renorm. is performed, T_RENORM overrides T_KLAT for lattice thermal conductivity # T_KLAT = {"t_min":100,"t_max":1500,"t_step":100} #[i*100 for i in range(0,11)] -T_KLAT = 300 # [i*100 for i in range(0,11)] +# T_KLAT = [300] # [i*100 for i in range(0,11)] +T_KLAT = {"min":100,"max":1000,"step":100} #[i*100 for i in range(0,11)] T_THERMAL_CONDUCTIVITY = [0, 100, 200, 300] # [i*100 for i in range(0,16)] IMAGINARY_TOL = 0.025 # in THz FIT_METHOD = "rfe" @@ -170,7 +173,7 @@ def hiphive_static_calcs( outputs["structure_data"] = json_saver.output[2] outputs["current_dir"] = json_saver.output[3] - return Response(addition=jobs, output=outputs) + return Response(replace=jobs, output=outputs) @job def generate_hiphive_displacements( @@ -447,6 +450,7 @@ def run_hiphive( """ logger.info(f"cutoffs = {cutoffs}") logger.info(f"disp_cut is {disp_cut}") + logger.info(f"prev_dir_json_saver is {prev_dir_json_saver} this is for def run_hiphive") copy_hiphive_outputs(prev_dir_json_saver) @@ -458,9 +462,12 @@ def run_hiphive( supercell_structure = structure_data["supercell_structure"] supercell_matrix = np.array(structure_data["supercell_matrix"]) + parent_structure = SpacegroupAnalyzer(parent_structure).find_primitive() #TODO refactor this later + if cutoffs is None: cutoffs = get_cutoffs(supercell_structure) logger.info(f"cutoffs is {cutoffs}") + # cutoffs = [cutoffs[7], cutoffs[18], cutoffs[21], cutoffs[24], cutoffs[25]] #TODO, for testing purpose. Remove this line later else: pass @@ -492,20 +499,20 @@ def run_hiphive( mean_displacements = np.linalg.norm(displacements, axis=1).mean() logger.info(f"Mean displacements while reading individual displacements: {mean_displacements}") - logger.info(f'forces in 0th structure are {structures[0].get_array("forces")}') - logger.info( - f'displacements in 0th structure are {structures[0].get_array("displacements")}' - ) + # logger.info(f'forces in 0th structure are {structures[0].get_array("forces")}') + # logger.info( + # f'displacements in 0th structure are {structures[0].get_array("displacements")}' + # ) - logger.info(f'forces in 1st structure are {structures[1].get_array("forces")}') - logger.info( - f'displacements in 1st structure are {structures[1].get_array("displacements")}' - ) + # logger.info(f'forces in 1st structure are {structures[1].get_array("forces")}') + # logger.info( + # f'displacements in 1st structure are {structures[1].get_array("displacements")}' + # ) - logger.info(f'forces in 2nd structure are {structures[2].get_array("forces")}') - logger.info( - f'displacements in 2nd structure are {structures[2].get_array("displacements")}' - ) + # logger.info(f'forces in 2nd structure are {structures[2].get_array("forces")}') + # logger.info( + # f'displacements in 2nd structure are {structures[2].get_array("displacements")}' + # ) all_cutoffs = cutoffs logger.info(f"all_cutoffs is {all_cutoffs}") @@ -560,10 +567,15 @@ def run_hiphive( if fitting_data["n_imaginary"] == 0: logger.info("No imaginary modes! Writing ShengBTE files") atoms = AseAtomsAdaptor.get_atoms(parent_structure) - fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) fcs.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") ForceConstants.write_to_phonopy(fcs, "fc2.hdf5", "hdf5") - ForceConstants.write_to_phono3py(fcs, "fc3.hdf5", 3) + ForceConstants.write_to_phono3py(fcs, "fc3.hdf5", "hdf5") + ### detour from hdf5 + supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) + else: logger.info("ShengBTE files not written due to imaginary modes.") logger.info("You may want to perform phonon renormalization.") @@ -947,9 +959,9 @@ def harmonic_properties( structure: Structure, supercell_matrix: np.ndarray, fcs: ForceConstants, - temperature: List, + temperature: list, imaginary_tol: float = IMAGINARY_TOL -) -> Tuple[Dict,Phonopy]: +) -> tuple[dict,Phonopy]: """ Uses the force constants to extract phonon properties. Used for comparing the accuracy of force constant fits. @@ -1129,11 +1141,6 @@ def gruneisen( grun_tot.append(get_total_grun(omega,grun,kweight,temp)) grun_tot = np.nan_to_num(np.array(grun_tot)) - # grun_tot = [ - # np.nan_to_num(np.array(get_total_grun(omega, grun, kweight, temp))) - # for temp in temperature - # ] - # linear thermal expansion coefficeint and fraction if bulk_modulus is None: cte = None @@ -1153,7 +1160,6 @@ def gruneisen( logger.info(f"bulk_modulus: {bulk_modulus}") # cte = grun_tot*heat_capacity.repeat(3)/(vol/10**30)/(bulk_modulus*10**9)/3 cte = grun_tot*heat_capacity.repeat(3).reshape(len(heat_capacity),3)/(vol/10**30)/(bulk_modulus*10**9)/3 - # cte = grun_tot*heat_capacity.repeat(3).reshape(1,3)/(vol/10**30)/(bulk_modulus*10**9)/3 cte = np.nan_to_num(cte) if len(temperature)==1: dLfrac = cte*temperature @@ -1214,11 +1220,11 @@ def run_thermal_cond_solver( file. """ if therm_cond_solver == "almabte": - therm_cond_solver_cmd = Atomate2Settings.ALMABTE_CMD + therm_cond_solver_cmd = Atomate2Settings().ALMABTE_CMD elif therm_cond_solver == "shengbte": - therm_cond_solver_cmd = Atomate2Settings.SHENGBTE_CMD + therm_cond_solver_cmd = Atomate2Settings().SHENGBTE_CMD elif therm_cond_solver == "phono3py": - therm_cond_solver_cmd = Atomate2Settings.PHONO3PY_CMD + therm_cond_solver_cmd = Atomate2Settings().PHONO3PY_CMD logger.info(f"therm_cond_solver_cmd = {therm_cond_solver_cmd}") therm_cond_solver_cmd = expandvars(therm_cond_solver_cmd) @@ -1233,10 +1239,14 @@ def run_thermal_cond_solver( structure_data = loadfn("structure_data.json") structure = structure_data["structure"] supercell_matrix = structure_data["supercell_matrix"] + + structure = SpacegroupAnalyzer(structure).find_primitive() #TODO refactor this later logger.info(f"Temperature = {temperature}") temperature = temperature if temperature is not None else T_KLAT + logger.info(f"Temperature = {temperature}") + logger.info(f"type of temperature = {type(temperature)}") renormalized = renormalized if renormalized is not None else False @@ -1246,9 +1256,9 @@ def run_thermal_cond_solver( if isinstance(temperature, (int, float)): pass elif isinstance(temperature, dict): - temperature["t_min"] - temperature["t_max"] - temperature["t_step"] + temperature["min"] + temperature["max"] + temperature["step"] else: raise ValueError("Unsupported temperature type, must be float or dict") @@ -1296,7 +1306,7 @@ def run_thermal_cond_solver( @job def run_fc_to_pdos( renormalized: bool | None = None, - renorm_temperature: str | None = None, + # renorm_temperature: str | None = None, mesh_density: float | None = None, prev_dir_json_saver: str | None = None, loop: int | None = None, @@ -1328,7 +1338,7 @@ def run_fc_to_pdos( logger.info(f"loop = {loop}") renormalized = renormalized if renormalized else False - renorm_temperature = renorm_temperature if renorm_temperature else None + # renorm_temperature = renorm_temperature if renorm_temperature else None mesh_density = mesh_density if mesh_density else 100.0 structure_data = loadfn(f"structure_data_{loop}.json") @@ -1366,7 +1376,7 @@ def run_fc_to_pdos( f"Finished inserting renormalized force constants and phonon data at {T} K" ) - return uniform_bs, lm_bs, dos + return uniform_bs, lm_bs, dos, prev_dir_json_saver def _get_fc_fsid(structure: Structure, @@ -1397,14 +1407,11 @@ def run_hiphive_renormalization( temperature: float, renorm_method: str, nconfig: int, - conv_thresh: float, - max_iter: int, renorm_TE_iter: bool, bulk_modulus: float, - mesh_density: float, prev_dir_hiphive: str, loop: int, -) -> list[str]: +) -> list[str, dict[str, Any]]: """ Phonon renormalization using hiPhive. @@ -1484,16 +1491,21 @@ def run_hiphive_renormalization( TD_structure_data["structure"] = parent_structure_TD TD_structure_data["supercell_structure"] = supercell_structure_TD + imag_modes_bool = TD_data["n_imaginary"] > 0 # True if imaginary modes exist + # Thermodynamic integration for anharmonic free energy - TD_data = thermodynamic_integration_ifc(TD_data, # everything TD - fcs, # original - param, # original - ) + TD_data = thermodynamic_integration_ifc( + TD_data, # everything TD + fcs, # original + param, # original + imag_modes_bool, # if False, only uses lambda=0 + ) # write results logger.info("Writing renormalized results") fcs_TD = TD_data['fcs'] fcs_TD.write("force_constants.fcs") - if "n_imaginary" in TD_data: + # if "n_imaginary" in TD_data: + if TD_data["n_imaginary"] != 0: thermal_keys = ["temperature","free_energy","entropy","heat_capacity", "free_energy_correction_S","free_energy_correction_SC", "free_energy_correction_TI"] @@ -1519,7 +1531,7 @@ def run_hiphive_renormalization( current_dir = os.getcwd() - return [current_dir] + return [current_dir, TD_thermal_data] def run_renormalization( @@ -1611,6 +1623,7 @@ def thermodynamic_integration_ifc( param: np.ndarray, lambda_array: np.ndarray = np.array([0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]), TI_nconfig=3, + imag_modes_bool: bool = True ) -> dict: supercell_structure = TD_data["supercell_structure"] cs = TD_data['cs'] @@ -1621,6 +1634,8 @@ def thermodynamic_integration_ifc( supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) renorm = Renormalization(cs, supercell_atoms, param, fcs, T, 'least_squares', 'rfe', param_TD, fcs_TD) matcov_TD, matcov_BO, matcov_TDBO = renorm.born_oppenheimer_qcv(TI_nconfig) + if not imag_modes_bool: + lambda_array = np.array([0]) correction_TI = renorm.thermodynamic_integration(lambda_array, matcov_TD, matcov_BO, matcov_TDBO, TI_nconfig) TD_data["free_energy_correction_TI"] = correction_TI return TD_data @@ -1699,6 +1714,7 @@ def run_lattice_thermal_conductivity( # files needed to run ShengBTE logger.info("We are in Lattice Thermal Conductivity... 1") + logger.info(f"previ_dir_hiphive in def run_lattice_thermal_conductivity = {prev_dir_hiphive}") if renormalized: @@ -1727,5 +1743,15 @@ def run_lattice_thermal_conductivity( loop=loop, therm_cond_solver=therm_cond_solver ) + shengbte.update_config({"manager_config": {"_fworker": "gpu_reg_fworker"}}) #change to gpu_fworker + shengbte.name += f" {temperature} {loop}" + shengbte.metadata.update( + { + "tag": [ + f"run_thermal_cond_solver_{loop}", + f"loop={loop}", + ] + } + ) - return Response(replace=shengbte) + return Response(addition=shengbte) From bad7f5e2d086771fdb67c961cfa05507cfbccff7 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 13:54:19 -0700 Subject: [PATCH 06/14] fixes 1. ruff fixes 2. updated class attributes 3. added bulk modulus to the metadata 4. #TODO: Implement Quality Control Job 5. updated fireworks config for each job 6. updated definition to extract phonon dispersion 7. updated definition to extract LTC, depending upon whether or not renormalization was performed. 8. updated name of the flow --- src/atomate2/common/flows/hiphive.py | 211 ++++++++++++++++----------- 1 file changed, 126 insertions(+), 85 deletions(-) diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index c58421c5a6..8ac63883f4 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -23,7 +23,6 @@ ) # Atomate2 packages -from atomate2.vasp.jobs.core import StaticMaker from atomate2.vasp.jobs.phonons import PhononDisplacementMaker from atomate2.vasp.sets.core import StaticSetGenerator @@ -97,9 +96,6 @@ class BaseHiphiveMaker(Maker, ABC): Temperatures for renormalization calculations, default is [1500]. T_KLAT (int): Temperature for lattice thermal conductivity calculation, default is 300. - T_THERMAL_CONDUCTIVITY (list): - Temperatures for thermal conductivity calculations, - default is [0, 100, 200, 300]. FIT_METHOD (str): Method for fitting force constants, default is "rfe". RENORM_METHOD (str): @@ -131,20 +127,13 @@ class BaseHiphiveMaker(Maker, ABC): i * 100 for i in range(21) ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) T_RENORM: ClassVar[list[int]] = [ - 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000 # 400 + 1500, 2000, 2500, 2700, 3000 # 300, 500, 600, 700, 800, 900, 1000 ] # [i*100 for i in range(0,16)] # Temp. at which renorm. is to be performed - # Temperature at which lattice thermal conductivity is calculated # If renormalization is performed, # T_RENORM overrides T_KLAT for lattice thermal conductivity - T_KLAT: ClassVar[list[int]] = [100, 200, 300, 400] # [i*100 for i in range(0,11)] - T_THERMAL_CONDUCTIVITY: ClassVar[list[int]] = [ - 0, - 100, - 200, - 300, - ] # [i*100 for i in range(0,16)] + T_KLAT: ClassVar[dict] = {"min":100,"max":1000,"step":100} #[i*100 for i in range(0,11)] FIT_METHOD = "rfe" - RENORM_METHOD = "least_squares" + RENORM_METHOD = "least_squares" # pseudoinverse RENORM_NCONFIG = 5 # Changed from 50 RENORM_CONV_THRESH = 0.1 # meV/atom RENORM_MAX_ITER = 30 # Changed from 20 @@ -239,6 +228,7 @@ def make( if self.prev_calc_dir_argname is not None: bulk_kwargs[self.prev_calc_dir_argname] = prev_dir bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) + bulk.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker jobs.append(bulk) outputs.append(bulk.output) structure = bulk.output.structure @@ -248,6 +238,7 @@ def make( { "tag": [ f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", f"relax_{loops}", f"nConfigsPerStd={n_structures}", f"fixedDispls={fixed_displs}", @@ -277,7 +268,6 @@ def make( ) jobs.append(static_calcs) - # 3. Hiphive Fitting of FCPs upto 4th order fit_force_constant = run_hiphive( fit_method=fit_method, @@ -291,12 +281,14 @@ def make( cutoffs=cutoffs ) fit_force_constant.name += f" {loops}" + fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) jobs.append(fit_force_constant) outputs.append(fit_force_constant.output) fit_force_constant.metadata.update( { "tag": [ f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", f"fit_force_constant_{loops}", f"nConfigsPerStd={n_structures}", f"fixedDispls={fixed_displs}", @@ -307,48 +299,50 @@ def make( } ) - # # 7. Quality Control Job to check if the desired Test RMSE is achieved, - # # if not, then increase the number of structures -- - # # Using "addintion" feature of jobflow - # loops += 1 - # n_structures += 1 - # logger.info(f"Number of structures increased to {n_structures}") - # logger.info(f"loop = {loops}") - # error_check_job = quality_control( - # rmse_test=fit_force_constant.output[5], - # n_structures=n_structures, - # fixedDispls=fixed_displs, - # loop=loops, - # fit_method=fit_method, - # disp_cut=disp_cut, - # bulk_modulus=bulk_modulus, - # temperature_qha=temperature_qha, - # mesh_density=mesh_density, - # imaginary_tol=imaginary_tol, - # prev_dir_json_saver=static_calcs.output["current_dir"], - # prev_dir=prev_dir, - # supercell_matrix=supercell_matrix, - # # supercell_matrix_kwargs=supercell_matrix_kwargs, - # ) - # error_check_job.name += f" {loops}" - # jobs.append(error_check_job) - # outputs.append(error_check_job.output) - # error_check_job.metadata.update( - # { - # "tag": [ - # f"error_check_job_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # # f"supercell_matrix_kwargs={supercell_matrix_kwargs}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) + #TODO: Implement Quality Control Job + # # # # # # 7. Quality Control Job to check if the desired Test RMSE is achieved, + # # # # # # if not, then increase the number of structures -- + # # # # # # Using "addintion" feature of jobflow + # # # # # loops += 1 + # # # # # n_structures += 1 + # # # # # logger.info(f"Number of structures increased to {n_structures}") + # # # # # logger.info(f"loop = {loops}") + # # # # # error_check_job = quality_control( + # # # # # rmse_test=fit_force_constant.output[5], + # # # # # n_structures=n_structures, + # # # # # fixedDispls=fixed_displs, + # # # # # loop=loops, + # # # # # fit_method=fit_method, + # # # # # disp_cut=disp_cut, + # # # # # bulk_modulus=bulk_modulus, + # # # # # temperature_qha=temperature_qha, + # # # # # mesh_density=mesh_density, + # # # # # imaginary_tol=imaginary_tol, + # # # # # prev_dir_json_saver=static_calcs.output["current_dir"], + # # # # # prev_dir=prev_dir, + # # # # # supercell_matrix=supercell_matrix, + # # # # # # supercell_matrix_kwargs=supercell_matrix_kwargs, + # # # # # ) + # # # # # error_check_job.name += f" {loops}" + # # # # # jobs.append(error_check_job) + # # # # # outputs.append(error_check_job.output) + # # # # # error_check_job.metadata.update( + # # # # # { + # # # # # "tag": [ + # # # # # f"error_check_job_{loops}", + # # # # # f"nConfigsPerStd={n_structures}", + # # # # # f"fixedDispls={fixed_displs}", + # # # # # f"dispCut={disp_cut}", + # # # # # # f"supercell_matrix_kwargs={supercell_matrix_kwargs}", + # # # # # f"supercell_matrix={supercell_matrix}", + # # # # # f"loop={loops}", + # # # # # ] + # # # # # } + # # # # # ) # 4. Perform phonon renormalization to obtain temperature-dependent # force constants using hiPhive + outputs_renorm = [] if renormalize: for temperature in renormalize_temperature: nconfig = renormalize_nconfig * (1 + temperature // 100) @@ -365,11 +359,15 @@ def make( loop=loops, ) renormalization.name += f" {temperature} {loops}" + renormalization.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) jobs.append(renormalization) + outputs_renorm.append(renormalization.output) outputs.append(renormalization.output) renormalization.metadata.update( { "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", f"run_renormalization_{loops}", f"nConfigsPerStd={n_structures}", f"fixedDispls={fixed_displs}", @@ -381,29 +379,22 @@ def make( ) # 5. Extract Phonon Band structure & DOS from FC - if renormalize: - fc_pdos_pb_to_db = run_fc_to_pdos( - renormalized=renormalize, - renorm_temperature=renormalize_temperature, - mesh_density=mesh_density, - prev_dir_json_saver=renormalization.output[0], - loop=loops, - ) - else: - fc_pdos_pb_to_db = run_fc_to_pdos( + # for 0K + fc_pdos_pb_to_db = run_fc_to_pdos( renormalized=renormalize, - renorm_temperature=renormalize_temperature, mesh_density=mesh_density, prev_dir_json_saver=fit_force_constant.output[4], loop=loops, ) - fc_pdos_pb_to_db.name += f" {loops}" + fc_pdos_pb_to_db.name += f" {loops} 0K" jobs.append(fc_pdos_pb_to_db) outputs.append(fc_pdos_pb_to_db.output) fc_pdos_pb_to_db.metadata.update( { "tag": [ f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"temperature=0K" f"fc_pdos_pb_to_db_{loops}", f"nConfigsPerStd={n_structures}", f"fixedDispls={fixed_displs}", @@ -413,7 +404,35 @@ def make( ] } ) + # for finite temperatures + if renormalize: + for i, temperature in enumerate(renormalize_temperature): + fc_pdos_pb_to_db = run_fc_to_pdos( + renormalized=renormalize, + mesh_density=mesh_density, + prev_dir_json_saver=outputs_renorm[i][0], + loop=loops, + ) + fc_pdos_pb_to_db.name += f" {loops} {temperature}K" + jobs.append(fc_pdos_pb_to_db) + outputs.append(fc_pdos_pb_to_db.output) + fc_pdos_pb_to_db.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"temperature={temperature}K" + f"fc_pdos_pb_to_db_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-03-28-05-45-16-259235-19145" # 6. Lattice thermal conductivity calculation using Sheng BTE if calculate_lattice_thermal_conductivity: if renormalize: @@ -434,36 +453,58 @@ def make( prev_dir_hiphive=fit_force_constant.output[4], therm_cond_solver= self.THERM_COND_SOLVER ) + lattice_thermal_conductivity.name += f" {temperatures} {loops}" + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker + jobs.append(lattice_thermal_conductivity) + outputs.append(lattice_thermal_conductivity.output) + lattice_thermal_conductivity.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"run_lattice_thermal_conductivity_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) else: - for _, T in enumerate(temperatures): + for t, T in enumerate(temperatures): if T == 0: continue lattice_thermal_conductivity = run_lattice_thermal_conductivity( renormalized=renormalize, temperature=T, loop=loops, - prev_dir_hiphive=fit_force_constant.output[4], + prev_dir_hiphive=outputs_renorm[t][0], therm_cond_solver= self.THERM_COND_SOLVER ) - lattice_thermal_conductivity.name += f" {loops}" - jobs.append(lattice_thermal_conductivity) - outputs.append(lattice_thermal_conductivity.output) - lattice_thermal_conductivity.metadata.update( - { - "tag": [ - f"run_lattice_thermal_conductivity_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - # f"supercell_matrix_kwargs={supercell_matrix_kwargs}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) + lattice_thermal_conductivity.name += f" {T} {loops}" + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker + jobs.append(lattice_thermal_conductivity) + outputs.append(lattice_thermal_conductivity.output) + lattice_thermal_conductivity.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"run_lattice_thermal_conductivity_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) - return Flow(jobs=jobs, output=outputs, name=f"{mpid}_ShengBTE_" + return Flow(jobs=jobs, output=outputs, name=f"{mpid}_{self.THERM_COND_SOLVER}_" f"{fixed_displs}_" f"{self.name}") From c45b616d4431d2ce67c806a87e3fb012285693cc Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 13:55:28 -0700 Subject: [PATCH 07/14] customized INCAR tags using tags in accordance with atomate LD workflow --- src/atomate2/vasp/flows/hiphive.py | 38 ++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/atomate2/vasp/flows/hiphive.py b/src/atomate2/vasp/flows/hiphive.py index 3ff759b4e9..56e5c884bd 100644 --- a/src/atomate2/vasp/flows/hiphive.py +++ b/src/atomate2/vasp/flows/hiphive.py @@ -84,17 +84,41 @@ class HiphiveMaker(BaseHiphiveMaker): input_set_generator = StaticSetGenerator( user_kpoints_settings={"reciprocal_density": 500}, user_incar_settings={ - "IBRION": 2, + # "IBRION": 2, + # "ISIF": 3, + # "ENCUT": 600, # Changed this from 600 + # "EDIFF": 1e-7, + # "LAECHG": False, + # "ALGO": "Normal", + # "NSW": 0, + # "LCHARG": False, + # "LREAL": True + + "ADDGRID": True, + "ALGO": "Normal", + "EDIFF": 1e-06, + "ENCUT": 600, + "GGA": "PS", + "IBRION": -1, "ISIF": 3, - "ENCUT": 700, # Changed this from 600 - "EDIFF": 1e-7, + "ISMEAR": 0, + "ISPIN": 2, "LAECHG": False, - "ALGO": "Normal", - "NSW": 0, + "LASPH": True, "LCHARG": False, - "LREAL": True + "LORBIT": 11, + "LREAL": "Auto", + "LVHAR": False, + "LVTOT": False, + "LWAVE": False, + # "MAGMOM": 250*0.6, + "NCORE": 6, + "NELM": 100, + "NSW": 0, + "PREC": "Accurate", + "SIGMA": 0.1, }, - auto_ispin=True, + # auto_ispin=True, ) ) ) From 7250e2a6dbb42f89ec4f5831ddde03c6f16dbbe4 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 13:58:01 -0700 Subject: [PATCH 08/14] adding LTC calculators 1. supports shengBTE, almaBTE & phono3py as of now. 2. before running any of these calculators, make sure to successfully compile them as explained on their respective websites. 3.Also, remember to put the execution command for each of these calculators in atomate2.yaml file --- src/atomate2/settings.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atomate2/settings.py b/src/atomate2/settings.py index 79b222660e..c57bdbaef4 100644 --- a/src/atomate2/settings.py +++ b/src/atomate2/settings.py @@ -194,12 +194,12 @@ class Atomate2Settings(BaseSettings): # ShengBTE settings SHENGBTE_CMD: str = Field( - "shengbte", description="Command to run 3 phonon ShengBTE." + "ShengBTE", description="Command to run 3 phonon ShengBTE." ) # AlmaBTE settings ALMABTE_CMD: str = Field( - "almabte", description="Command to run 4 phonon AlmaBTE." + "ShengBTE", description="Command to run 4 phonon AlmaBTE." ) # Phono3py settings From 82ef1c0649fefafd8e29fe85fca9185b65596c34 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 13:58:54 -0700 Subject: [PATCH 09/14] added anharmonicphonons block in pyproject.yaml --- pyproject.toml | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 2c2c3c22ee..c535cad1f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,7 +52,10 @@ forcefields = [ "quippy-ase>=0.9.14", ] anharmonicphonons = [ - "phono3py", + "phonopy>=2.21.0", + "hiphive @ git+https://gitlab.com/jsyony37/hiphive.git@personal#egg=hiphive", + "f90nml==1.4.4", + "spglib>=2.2.0", "ase>=3.22.1",] docs = [ "FireWorks==2.0.3", @@ -96,7 +99,6 @@ strict = [ "seekpath==2.1.0", "typing-extensions==4.10.0", "python-ulid==2.2.0", - "phono3py==2.5.1" ] [project.scripts] @@ -170,6 +172,7 @@ lint.ignore = [ "FBT002", "FIX002", "G004", # logging uses fstring + "ISC001", "PD011", # pandas-use-of-dot-values "PERF203", # try-except-in-loop "PLR", # pylint-refactor @@ -179,10 +182,10 @@ lint.ignore = [ "PTH", # prefer Pathlib to os.path "RUF013", # implicit-optional "S324", # use of insecure hash function + "S507", # paramiko auto trust "SLF", # private member accessed outside class "TD", # TODOs "TRY003", # long message outside exception class - "S507", # paramiko auto trust ] lint.pydocstyle.convention = "numpy" lint.isort.known-first-party = ["atomate2"] From 8889d3edfccb91080d08372762629d24fe9be14d Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Fri, 29 Mar 2024 15:03:31 -0700 Subject: [PATCH 10/14] update the logic to use lambda=0 for TI the decision is now made based on whether n_imaginary > 0 from the fitting_data.json @ 0K --- src/atomate2/common/jobs/hiphive.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index 75095400c8..31a97e90c6 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -1440,6 +1440,9 @@ def run_hiphive_renormalization( # thermal_data = loadfn("thermal_data.json") # thermal_data = thermal_data["heat_capacity"] + n_imaginary_orig = fitting_data["n_imaginary"] + imag_modes_bool = n_imaginary_orig > 0 # True if imaginary modes exist + cutoffs = fitting_data["cutoffs"] fit_method = fitting_data["fit_method"] @@ -1491,8 +1494,6 @@ def run_hiphive_renormalization( TD_structure_data["structure"] = parent_structure_TD TD_structure_data["supercell_structure"] = supercell_structure_TD - imag_modes_bool = TD_data["n_imaginary"] > 0 # True if imaginary modes exist - # Thermodynamic integration for anharmonic free energy TD_data = thermodynamic_integration_ifc( TD_data, # everything TD From 6e9da6710ee82a6d61f9ae2a998349a595899649 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Wed, 3 Apr 2024 18:15:32 -0700 Subject: [PATCH 11/14] updating the logic to generate fixed displacements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. largest displacement is based on the avg row (periodic table) number of the compound. 2. displacements are skewed towards the large displacement 3. # of displaced structures to generate depends on the n_structures, which depends upon the symmetry of the primitive structure 4. smallest displacement is fixed at 0.01 Å --- src/atomate2/common/jobs/hiphive.py | 162 +++++++++++++++++++++------- 1 file changed, 121 insertions(+), 41 deletions(-) diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index 31a97e90c6..f7dd4c3bb8 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -13,6 +13,7 @@ from itertools import product from os.path import expandvars from typing import TYPE_CHECKING, Any, Union +from dataclasses import field import numpy as np import phonopy as phpy @@ -39,7 +40,7 @@ from hiphive.utilities import get_displacements # Jobflow packages -from jobflow import Response, job +from jobflow import Flow, Response, job from joblib import Parallel, delayed # Pymatgen packages @@ -75,7 +76,7 @@ logger = initialize_logger(level=3) T_QHA = [ - i * 100 for i in range(21) + i * 100 for i in range(31) ] # Temp. for straight-up phonopy calculation of thermo. props. (free energy etc.) # Temperature at which lattice thermal conductivity is calculated # If renorm. is performed, T_RENORM overrides T_KLAT for lattice thermal conductivity @@ -94,6 +95,7 @@ "hiphive_static_calcs", "generate_hiphive_displacements", "quality_control", + "hiphive_cutoff" "run_hiphive", "run_thermal_cond_solver", "run_fc_to_pdos", @@ -111,7 +113,7 @@ def hiphive_static_calcs( min_length: float | None = None, prefer_90_degrees: bool = True, n_structures: int | None = None, - fixed_displs: list[float] = [0.01, 0.03, 0.08, 0.1], + # fixed_displs: list[float] = [0.01, 0.03, 0.08, 0.1], loops: int | None = None, prev_dir: str | None = None, phonon_displacement_maker: BaseVaspMaker | None = None, @@ -144,7 +146,7 @@ def hiphive_static_calcs( structure=structure, supercell_matrix=supercell_matrix, n_structures=n_structures, - fixed_displs=fixed_displs, + # fixed_displs=fixed_displs, loop=loops, ) jobs.append(displacements) @@ -185,34 +187,55 @@ def generate_hiphive_displacements( ) -> list[Structure]: """Job generates the perturbed structures for hiPhive fitting.""" if fixed_displs is None: - fixed_displs = [0.01, 0.03, 0.08, 0.1] + + # fixed_displs = [0.01, 0.03, 0.08, 0.1] #TODO, this list is used in + # the paper + smallest_disp = 0.01 + + # dictionary of largest displacements for each period + largest_displacement = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} + + row = int( + np.around(np.array([s.row for s in structure.species]).mean(), 0) + ) + + largest_disp = largest_displacement[row] + + ratio = int(largest_disp/smallest_disp) + factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) + amplitudes = (smallest_disp*factors).round(3) + logger.info(f"amplitudes = {amplitudes}") + fixed_displs = amplitudes.tolist() + logger.info(f"list_amplitudes = {fixed_displs}") + logger.info(f"supercell_matrix = {supercell_matrix}") supercell_structure = SupercellTransformation( scaling_matrix=supercell_matrix ).apply_transformation(structure) logger.info(f"supercell_structure = {supercell_structure}") - # Generate the rattled structures #### structures_ase_all = [] + logger.info(f"n_structures = {n_structures}") # Convert to ASE atoms for i in range(len(fixed_displs)): supercell_ase = AseAtomsAdaptor.get_atoms(supercell_structure) structures_ase = generate_displaced_structures( - supercell_ase, n_structures, fixed_displs[i], loop + atoms=supercell_ase, n_structures=1, distance=fixed_displs[i], loop=loop ) - structures_ase_all.append(structures_ase) + structures_ase_all.extend(structures_ase) logger.info(f"structures_ase_all: {structures_ase_all}") + logger.info(f"len(structures_ase_all): {len(structures_ase_all)}") # Convert back to pymatgen structure structures_pymatgen = [] for atoms_ase in range(len(structures_ase_all)): logger.info(f"atoms: {atoms_ase}") - logger.info(f"structures_ase_all[atoms]: {structures_ase_all[atoms_ase][0]}") - structure_i = AseAtomsAdaptor.get_structure(structures_ase_all[atoms_ase][0]) + logger.info(f"structures_ase_all[atoms]: {structures_ase_all[atoms_ase]}") + structure_i = AseAtomsAdaptor.get_structure(structures_ase_all[atoms_ase]) structures_pymatgen.append(structure_i) - + logger.info(f"len(structures_pymatgen): {len(structures_pymatgen)}") for i in range(len(structures_pymatgen)): structures_pymatgen[i].to(f"POSCAR_{i}", "poscar") @@ -222,8 +245,8 @@ def generate_hiphive_displacements( def collect_perturbed_structures( structure: Structure, supercell_matrix: list[list[int]] | None = None, - rattled_structures: list[Structure] | None = None, # Add type annotation for rattled_structures - forces: Any | None = None, # Add type annotation for forces + rattled_structures: list[Structure] | None = None, + forces: Any | None = None, loop: int | None = None, prev_dir_json_saver: str | None = None, ) -> list: @@ -279,7 +302,6 @@ def collect_perturbed_structures( with open(f"perturbed_forces_{loop}_new.json", "w") as f: json.dump(output, f) - # all_forces_loop = loadfn(f"perturbed_forces_{loop}_new.json") with open(f"perturbed_forces_{loop}_new.json") as file: all_forces_loop = json.load(file) @@ -414,6 +436,64 @@ def collect_perturbed_structures( # quality_control_job.name = f"quality_control_job {loop}" # return Response(addition=flow) +@job +def hiphive_cutoff( + mpid: str = None, + cutoffs: list[list] | None = None, + fit_method: str | None = None, + disp_cut: float | None = None, + bulk_modulus: float | None = None, + temperature_qha: float | None = None, + imaginary_tol: float | None = None, + prev_dir_json_saver: str | None = None, + loop: int | None = None, +) -> None: + """Run hiPhive with different cutoffs.""" + copy_hiphive_outputs(prev_dir_json_saver) + + structure_data = loadfn(f"structure_data_{loop}.json") + + supercell_structure = structure_data["supercell_structure"] + + if cutoffs is None: + cutoffs = get_cutoffs(supercell_structure) + logger.info(f"cutoffs is {cutoffs}") + else: + pass + + jobs_hiphive = [] + outputs_hiphive = [] + for _, i in enumerate(cutoffs): + logger.info(f"cutoffs is {i}") + run_hiphive_cutoffs = run_hiphive( + cutoffs=[i], + fit_method=fit_method, + disp_cut=disp_cut, + bulk_modulus=bulk_modulus, + temperature_qha=temperature_qha, + imaginary_tol=imaginary_tol, + prev_dir_json_saver=prev_dir_json_saver, + loop=loop, + ) + run_hiphive_cutoffs.name += f" {loop} {i}" + run_hiphive_cutoffs.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) + jobs_hiphive.append(run_hiphive_cutoffs) + outputs_hiphive.append(run_hiphive_cutoffs.output) + run_hiphive_cutoffs.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"cutoffs={i}", + f"fit_method={fit_method}", + f"disp_cut={disp_cut}", + f"bulk_modulus={bulk_modulus}", + f"imaginary_tol={imaginary_tol}", + f"prev_dir_json_saver={prev_dir_json_saver}", + ] + } + ) + hiphive_cutoff_flow = Flow(jobs_hiphive, outputs_hiphive) + return Response(replace=hiphive_cutoff_flow) @job @@ -467,7 +547,6 @@ def run_hiphive( if cutoffs is None: cutoffs = get_cutoffs(supercell_structure) logger.info(f"cutoffs is {cutoffs}") - # cutoffs = [cutoffs[7], cutoffs[18], cutoffs[21], cutoffs[24], cutoffs[25]] #TODO, for testing purpose. Remove this line later else: pass @@ -499,21 +578,6 @@ def run_hiphive( mean_displacements = np.linalg.norm(displacements, axis=1).mean() logger.info(f"Mean displacements while reading individual displacements: {mean_displacements}") - # logger.info(f'forces in 0th structure are {structures[0].get_array("forces")}') - # logger.info( - # f'displacements in 0th structure are {structures[0].get_array("displacements")}' - # ) - - # logger.info(f'forces in 1st structure are {structures[1].get_array("forces")}') - # logger.info( - # f'displacements in 1st structure are {structures[1].get_array("displacements")}' - # ) - - # logger.info(f'forces in 2nd structure are {structures[2].get_array("forces")}') - # logger.info( - # f'displacements in 2nd structure are {structures[2].get_array("displacements")}' - # ) - all_cutoffs = cutoffs logger.info(f"all_cutoffs is {all_cutoffs}") @@ -981,6 +1045,7 @@ def harmonic_properties( logger.info('Evaluating harmonic properties...') fcs2 = fcs.get_fc_array(2) fcs3 = fcs.get_fc_array(3) + logger.info('fcs2 & fcs3 read...') parent_phonopy = get_phonopy_structure(structure) phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) natom = phonopy.primitive.get_number_of_atoms() @@ -1437,19 +1502,12 @@ def run_hiphive_renormalization( structure_data = loadfn(f"structure_data_{loop}.json") phonopy_orig = phpy.load("phonopy_params.yaml") - # thermal_data = loadfn("thermal_data.json") - # thermal_data = thermal_data["heat_capacity"] - n_imaginary_orig = fitting_data["n_imaginary"] imag_modes_bool = n_imaginary_orig > 0 # True if imaginary modes exist cutoffs = fitting_data["cutoffs"] fit_method = fitting_data["fit_method"] - - # fitting_data["n_imaginary"] - # fitting_data["imaginary_tol"] - parent_structure = structure_data["structure"] supercell_structure = structure_data["supercell_structure"] supercell_matrix = np.array(structure_data["supercell_matrix"]) @@ -1522,10 +1580,32 @@ def run_hiphive_renormalization( logger.info("DEBUG: ",TD_data) if TD_data["n_imaginary"] > 0: logger.warning('Imaginary modes remain still exist') - logger.warning('ShengBTE FORCE_CONSTANTS_2ND not written') + logger.warning('ShengBTE FORCE_CONSTANTS_2ND & FORCE_CONSTANTS_3RD not written') else: - logger.info("No imaginary modes! Writing ShengBTE FORCE_CONSTANTS_2ND...") - fcs_TD.write_to_phonopy("FORCE_CONSTANTS_2ND".format(temperature), format="text") + logger.info("No imaginary modes! Writing ShengBTE files") + + parent_atoms_TD = copy(parent_atoms) + logger.info(f"TD_data exp frac: {TD_data['expansion_fraction']}") + logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,0]}") + logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,1]}") + logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,2]}") + new_cell = Cell(np.transpose([parent_atoms_TD.get_cell()[:,i]*(1+TD_data["expansion_fraction"][0,i]) for i in range(3)])) + parent_atoms_TD.set_cell(new_cell,scale_atoms=True) + + prim_TD_phonopy = PhonopyAtoms(symbols=parent_atoms_TD.get_chemical_symbols(), + scaled_positions=parent_atoms_TD.get_scaled_positions(), + cell=parent_atoms_TD.cell) + phonopy_TD = Phonopy(prim_TD_phonopy, supercell_matrix=supercell_matrix, primitive_matrix=None) + + atoms = AseAtomsAdaptor.get_atoms(parent_structure_TD) + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) + fcs_TD.write_to_phonopy(f"FORCE_CONSTANTS_2ND_{temperature}", format="text") + ForceConstants.write_to_phonopy(fcs_TD, "fc2.hdf5", "hdf5") + ForceConstants.write_to_phono3py(fcs_TD, "fc3.hdf5", "hdf5") + ### detour from hdf5 + supercell_atoms = phonopy_atoms_to_ase(phonopy_TD.supercell) + FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD_{temperature}", atoms, order=3, fc_tol=1e-4) dumpfn(TD_structure_data, "structure_data.json") dumpfn(TD_thermal_data, "thermal_data.json") @@ -1622,9 +1702,9 @@ def thermodynamic_integration_ifc( TD_data: dict, fcs: ForceConstants, param: np.ndarray, + imag_modes_bool: bool = True, lambda_array: np.ndarray = np.array([0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]), TI_nconfig=3, - imag_modes_bool: bool = True ) -> dict: supercell_structure = TD_data["supercell_structure"] cs = TD_data['cs'] From c18fe9342311978a51ae6f97f678bbe111388af7 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Thu, 4 Apr 2024 17:38:35 -0700 Subject: [PATCH 12/14] can now run hiphive fitting individually 1. added a job that creates multiples sub jobs for each hiphive fitting and then collects the fitting data & the fcs. fcp, cluster space and other properties corresponding to the best cutoff! 2. setting the min # of displ to 4 & max # of displ to 60 --- src/atomate2/common/flows/hiphive.py | 64 +++++++---- src/atomate2/common/jobs/hiphive.py | 157 +++++++++++++++++++++++---- 2 files changed, 174 insertions(+), 47 deletions(-) diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index 8ac63883f4..5ed1c6705f 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -18,6 +18,7 @@ # quality_control, run_fc_to_pdos, run_hiphive, + run_hiphive_individually, run_hiphive_renormalization, run_lattice_thermal_conductivity, ) @@ -35,7 +36,6 @@ from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker from atomate2.vasp.flows.core import DoubleRelaxMaker from atomate2.vasp.jobs.base import BaseVaspMaker -from emmet.core.math import Matrix3D logger = logging.getLogger(__name__) @@ -127,7 +127,7 @@ class BaseHiphiveMaker(Maker, ABC): i * 100 for i in range(21) ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) T_RENORM: ClassVar[list[int]] = [ - 1500, 2000, 2500, 2700, 3000 # 300, 500, 600, 700, 800, 900, 1000 + 10 # 300, 500, 600, 700, 800, 900, 1000, 1500, 2500, 2700, 3000 ] # [i*100 for i in range(0,16)] # Temp. at which renorm. is to be performed # If renormalization is performed, # T_RENORM overrides T_KLAT for lattice thermal conductivity @@ -228,7 +228,7 @@ def make( if self.prev_calc_dir_argname is not None: bulk_kwargs[self.prev_calc_dir_argname] = prev_dir bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) - bulk.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker + bulk.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) jobs.append(bulk) outputs.append(bulk.output) structure = bulk.output.structure @@ -260,7 +260,7 @@ def make( min_length=self.min_length, prefer_90_degrees=self.prefer_90_degrees, n_structures=n_structures, - fixed_displs=fixed_displs, + # fixed_displs=fixed_displs, loops=loops, prev_dir=prev_dir, phonon_displacement_maker=self.phonon_displacement_maker, @@ -268,18 +268,34 @@ def make( ) jobs.append(static_calcs) + # 3. Hiphive Fitting of FCPs upto 4th order - fit_force_constant = run_hiphive( - fit_method=fit_method, - disp_cut=disp_cut, - bulk_modulus=bulk_modulus, - temperature_qha=temperature_qha, - mesh_density=mesh_density, - imaginary_tol=imaginary_tol, - prev_dir_json_saver=static_calcs.output["current_dir"], - loop=loops, - cutoffs=cutoffs - ) + if n_structures >= 10: + fit_force_constant = run_hiphive_individually( + mpid = mpid, + cutoffs = cutoffs, + fit_method = fit_method, + disp_cut = disp_cut, + bulk_modulus = bulk_modulus, + temperature_qha = temperature_qha, + imaginary_tol = imaginary_tol, + prev_dir_json_saver = static_calcs.output["current_dir"], + # prev_dir_json_saver = prev_dir_json_saver, + loop = loops, + ) + else: + fit_force_constant = run_hiphive( + fit_method=fit_method, + disp_cut=disp_cut, + bulk_modulus=bulk_modulus, + temperature_qha=temperature_qha, + mesh_density=mesh_density, + imaginary_tol=imaginary_tol, + prev_dir_json_saver=static_calcs.output["current_dir"], + # prev_dir_json_saver=prev_dir_json_saver, + loop=loops, + cutoffs=cutoffs + ) fit_force_constant.name += f" {loops}" fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) jobs.append(fit_force_constant) @@ -350,12 +366,10 @@ def make( temperature=temperature, renorm_method=renormalize_method, nconfig=nconfig, - conv_thresh=renormalize_conv_thresh, - max_iter=renormalize_max_iter, renorm_TE_iter=renormalize_thermal_expansion_iter, bulk_modulus=bulk_modulus, - mesh_density=mesh_density, - prev_dir_hiphive=fit_force_constant.output[4], + prev_dir_hiphive=fit_force_constant.output["current_dir"], + # prev_dir_hiphive=prev_dir_hiphive, loop=loops, ) renormalization.name += f" {temperature} {loops}" @@ -378,12 +392,14 @@ def make( } ) + # 5. Extract Phonon Band structure & DOS from FC # for 0K fc_pdos_pb_to_db = run_fc_to_pdos( renormalized=renormalize, mesh_density=mesh_density, - prev_dir_json_saver=fit_force_constant.output[4], + prev_dir_json_saver=fit_force_constant.output["current_dir"], + # prev_dir_json_saver=prev_dir_hiphive, loop=loops, ) fc_pdos_pb_to_db.name += f" {loops} 0K" @@ -432,7 +448,7 @@ def make( } ) - # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-03-28-05-45-16-259235-19145" + # 6. Lattice thermal conductivity calculation using Sheng BTE if calculate_lattice_thermal_conductivity: if renormalize: @@ -450,11 +466,11 @@ def make( renormalized=renormalize, temperature=temperatures, loop=loops, - prev_dir_hiphive=fit_force_constant.output[4], + prev_dir_hiphive=fit_force_constant.output["current_dir"], therm_cond_solver= self.THERM_COND_SOLVER ) lattice_thermal_conductivity.name += f" {temperatures} {loops}" - lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) jobs.append(lattice_thermal_conductivity) outputs.append(lattice_thermal_conductivity.output) lattice_thermal_conductivity.metadata.update( @@ -485,7 +501,7 @@ def make( ) lattice_thermal_conductivity.name += f" {T} {loops}" - lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #change to gpu_fworker + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) jobs.append(lattice_thermal_conductivity) outputs.append(lattice_thermal_conductivity.output) lattice_thermal_conductivity.metadata.update( diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index f7dd4c3bb8..864b5cf70d 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -12,14 +12,12 @@ from copy import copy from itertools import product from os.path import expandvars -from typing import TYPE_CHECKING, Any, Union -from dataclasses import field +from typing import TYPE_CHECKING, Any import numpy as np import phonopy as phpy import psutil import scipy as sp -from ase import atoms from ase.cell import Cell # Hiphive packages @@ -95,7 +93,7 @@ "hiphive_static_calcs", "generate_hiphive_displacements", "quality_control", - "hiphive_cutoff" + "run_hiphive_individually" "run_hiphive", "run_thermal_cond_solver", "run_fc_to_pdos", @@ -202,7 +200,13 @@ def generate_hiphive_displacements( largest_disp = largest_displacement[row] ratio = int(largest_disp/smallest_disp) - factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) + if 60 >= n_structures >= 4: + logger.info(f"n_structures inside if statement >= 4 is {n_structures}") + factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) + elif n_structures < 4 : + factors = np.sqrt(np.linspace(1,ratio**2,4)) + else: + factors = np.sqrt(np.linspace(1,ratio**2,60)) amplitudes = (smallest_disp*factors).round(3) logger.info(f"amplitudes = {amplitudes}") fixed_displs = amplitudes.tolist() @@ -437,7 +441,7 @@ def collect_perturbed_structures( # return Response(addition=flow) @job -def hiphive_cutoff( +def run_hiphive_individually( mpid: str = None, cutoffs: list[list] | None = None, fit_method: str | None = None, @@ -461,7 +465,14 @@ def hiphive_cutoff( else: pass - jobs_hiphive = [] + jobs = [] + outputs: dict[str, list] = { + "thermal_data": [], + "anharmonic_data": [], + "fitting_data": [], + "param": [], + "current_dir": [] + } outputs_hiphive = [] for _, i in enumerate(cutoffs): logger.info(f"cutoffs is {i}") @@ -477,7 +488,7 @@ def hiphive_cutoff( ) run_hiphive_cutoffs.name += f" {loop} {i}" run_hiphive_cutoffs.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) - jobs_hiphive.append(run_hiphive_cutoffs) + jobs.append(run_hiphive_cutoffs) outputs_hiphive.append(run_hiphive_cutoffs.output) run_hiphive_cutoffs.metadata.update( { @@ -492,9 +503,92 @@ def hiphive_cutoff( ] } ) - hiphive_cutoff_flow = Flow(jobs_hiphive, outputs_hiphive) - return Response(replace=hiphive_cutoff_flow) + job_collect_hiphive_outputs = collect_hiphive_outputs( + fit_method=fit_method, + disp_cut=disp_cut, + imaginary_tol=imaginary_tol, + outputs=outputs_hiphive + ) + job_collect_hiphive_outputs.name += f" {loop} {i}" + job_collect_hiphive_outputs.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + jobs.append(job_collect_hiphive_outputs) + + outputs["thermal_data"] = job_collect_hiphive_outputs.output[0] + outputs["anharmonic_data"] = job_collect_hiphive_outputs.output[1] + outputs["fitting_data"] = job_collect_hiphive_outputs.output[2] + outputs["param"] = job_collect_hiphive_outputs.output[3] + outputs["current_dir"] = job_collect_hiphive_outputs.output[4] + + job_collect_hiphive_outputs.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"cutoffs={i}", + f"fit_method={fit_method}", + f"disp_cut={disp_cut}", + f"bulk_modulus={bulk_modulus}", + f"imaginary_tol={imaginary_tol}", + f"prev_dir_json_saver={prev_dir_json_saver}", + ] + } + ) + + return Response(replace=jobs, output=outputs) + +@job +def collect_hiphive_outputs( + fit_method: str | None = None, + disp_cut: float | None = None, + imaginary_tol: float | None = None, + outputs: list[dict] | None = None, +) -> list : + logger.info("We are in collect_hiphive_outputs") + + # Initialize best_fit with high initial values for comparison + fitting_data: dict[str, Any] = { + "cutoffs": [], + "rmse_test": [], + "fit_method": fit_method, + "disp_cut": disp_cut, + "imaginary_tol": imaginary_tol, + "best_cutoff": None, + "best_rmse": np.inf, + "n_imaginary": None + } + + best_fit = { + "rmse_test": np.inf, + "directory": None, + } + # Assuming outputs_hiphive is a list of dictionaries with the results + for result in outputs: + if result is None: + continue + + # Assuming result is a dictionary with keys: 'cutoffs', 'rmse_test', etc. + fitting_data["cutoffs"].append(result["fitting_data"]["cutoffs"][0]) + fitting_data["rmse_test"].append(result["fitting_data"]["rmse_test"][0]) + + # Update best_fit if the current result has a lower rmse_test + # Add additional conditions as needed + if ( + result["fitting_data"]["rmse_test"][0] < best_fit["rmse_test"] + ): + best_fit["directory"] = result["current_dir"] + fitting_data["best_cutoff"] = result["fitting_data"]["cutoffs"][0] + fitting_data["best_rmse"] = result["fitting_data"]["rmse_test"][0] + best_fit["rmse_test"] = result["fitting_data"]["rmse_test"][0] + fitting_data["n_imaginary"] = result["fitting_data"]["n_imaginary"] + + copy_hiphive_outputs(best_fit["directory"]) + thermal_data = loadfn("thermal_data.json") + dumpfn(fitting_data, "fitting_data.json") + param = np.loadtxt("parameters.txt") + + current_dir = os.getcwd() + logger.info(f"current_dir = {current_dir}") + return [thermal_data, thermal_data, fitting_data, param, current_dir] @job def run_hiphive( @@ -503,7 +597,6 @@ def run_hiphive( disp_cut: float | None = None, bulk_modulus: float | None = None, temperature_qha: float | None = None, - mesh_density: float | None = None, imaginary_tol: float | None = None, prev_dir_json_saver: str | None = None, loop: int | None = None, @@ -546,6 +639,7 @@ def run_hiphive( if cutoffs is None: cutoffs = get_cutoffs(supercell_structure) + cutoffs = [[6, 3.5, 3]] logger.info(f"cutoffs is {cutoffs}") else: pass @@ -576,7 +670,8 @@ def run_hiphive( structures.append(atoms) mean_displacements = np.linalg.norm(displacements, axis=1).mean() - logger.info(f"Mean displacements while reading individual displacements: {mean_displacements}") + logger.info(f"Mean displacements while reading individual displacements: " + f"{mean_displacements}") all_cutoffs = cutoffs logger.info(f"all_cutoffs is {all_cutoffs}") @@ -594,6 +689,7 @@ def run_hiphive( if fcs is None: raise RuntimeError("Could not find a force constant solution") + logger.info("Saving Harmonic props") thermal_data, phonopy = harmonic_properties( parent_structure, supercell_matrix, fcs, t_qha, imaginary_tol ) @@ -606,25 +702,31 @@ def run_hiphive( thermal_data["n_imaginary"], bulk_modulus, ) - + logger.info("Saving phonopy_params") phonopy.save("phonopy_params.yaml") fitting_data["n_imaginary"] = thermal_data.pop("n_imaginary") thermal_data.update(anharmonic_data) + logger.info("Saving fitting_data") dumpfn(fitting_data, "fitting_data.json") + logger.info("Saving thermal_data") dumpfn(thermal_data, "thermal_data.json") logger.info("Writing cluster space and force_constants") logger.info(f"{type(fcs)}") if isinstance(fcs, ForceConstants): + logger.info("Writing force_constants") fcs.write("force_constants.fcs") else: logger.info("fcs is not an instance of ForceConstants") + logger.info("Saving parameters") np.savetxt("parameters.txt", param) if isinstance(cs, ClusterSpace): + logger.info("Writing cluster_space") cs.write("cluster_space.cs") + logger.info("cluster_space writing is complete") else: logger.info("cs is not an instance of ClusterSpace") @@ -639,14 +741,21 @@ def run_hiphive( supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) - else: logger.info("ShengBTE files not written due to imaginary modes.") logger.info("You may want to perform phonon renormalization.") current_dir = os.getcwd() - return [thermal_data, anharmonic_data, fitting_data, param, current_dir] + outputs: dict[str, list] = { + "thermal_data": thermal_data, + "anharmonic_data": anharmonic_data, + "fitting_data": fitting_data, + "param": param, + "current_dir": current_dir + } + + return outputs def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: @@ -790,7 +899,8 @@ def fit_force_constants( "disp_cut": disp_cut, "imaginary_tol": imaginary_tol, # "max_n_imaginary": max_n_imaginary, - "best": None, + "best_cutoff": None, + "best_rmse": np.inf } best_fit = { @@ -843,7 +953,8 @@ def fit_force_constants( # and result["n_imaginary"] < best_fit["n_imaginary"] ): best_fit.update(result) - fitting_data["best"] = result["cutoffs"] + fitting_data["best_cutoff"] = result["cutoffs"] + fitting_data["best_rmse"] = result["rmse_test"] logger.info("Finished fitting force constants.") @@ -1230,10 +1341,10 @@ def gruneisen( dLfrac = cte*temperature else: dLfrac = thermal_expansion(temperature, cte) - logger.info('Gruneisen: \n {}'.format(grun_tot)) - logger.info('Coefficient of Thermal Expansion: \n {}'.format(cte)) - logger.info('Linear Expansion Fraction: \n {}'.format(dLfrac)) - + logger.info(f"Gruneisen: \n {grun_tot}") + logger.info(f"Coefficient of Thermal Expansion: \n {cte}") + logger.info(f"Linear Expansion Fraction: \n {dLfrac}") + return grun_tot, cte, dLfrac @@ -1304,7 +1415,7 @@ def run_thermal_cond_solver( structure_data = loadfn("structure_data.json") structure = structure_data["structure"] supercell_matrix = structure_data["supercell_matrix"] - + structure = SpacegroupAnalyzer(structure).find_primitive() #TODO refactor this later logger.info(f"Temperature = {temperature}") @@ -1327,7 +1438,7 @@ def run_thermal_cond_solver( else: raise ValueError("Unsupported temperature type, must be float or dict") - logger.info(f"Creating control dict") + logger.info("Creating control dict") control_dict = { "scalebroad": 0.5, From 60e9a7b7440eb0ac8fa4f3630008009e80804bf8 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Tue, 30 Apr 2024 09:11:51 -0700 Subject: [PATCH 13/14] adding fw configs to individual jobs --- src/atomate2/common/flows/elastic.py | 1 + src/atomate2/common/flows/hiphive.py | 492 +++++++++++++++------------ src/atomate2/common/jobs/elastic.py | 1 + src/atomate2/common/jobs/hiphive.py | 361 +++++++++++++++----- src/atomate2/common/jobs/phonons.py | 2 + src/atomate2/vasp/flows/hiphive.py | 123 +++++-- 6 files changed, 647 insertions(+), 333 deletions(-) diff --git a/src/atomate2/common/flows/elastic.py b/src/atomate2/common/flows/elastic.py index 795a092659..bbda7b1b2c 100644 --- a/src/atomate2/common/flows/elastic.py +++ b/src/atomate2/common/flows/elastic.py @@ -109,6 +109,7 @@ def make( if self.prev_calc_dir_argname is not None: bulk_kwargs[self.prev_calc_dir_argname] = prev_dir bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) + bulk.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) #TODO remove this line jobs.append(bulk) structure = bulk.output.structure prev_dir = bulk.output.dir_name diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index 5ed1c6705f..5431bda9d0 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -127,12 +127,12 @@ class BaseHiphiveMaker(Maker, ABC): i * 100 for i in range(21) ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) T_RENORM: ClassVar[list[int]] = [ - 10 # 300, 500, 600, 700, 800, 900, 1000, 1500, 2500, 2700, 3000 + 100, 500, 1000, 1500 # 300, 500, 600, 700, 800, 900, 1000, 1500, 2500, 2700, 3000 ] # [i*100 for i in range(0,16)] # Temp. at which renorm. is to be performed # If renormalization is performed, # T_RENORM overrides T_KLAT for lattice thermal conductivity T_KLAT: ClassVar[dict] = {"min":100,"max":1000,"step":100} #[i*100 for i in range(0,11)] - FIT_METHOD = "rfe" + FIT_METHOD = "omp" #least-squares #omp #rfe #elasticnet RENORM_METHOD = "least_squares" # pseudoinverse RENORM_NCONFIG = 5 # Changed from 50 RENORM_CONV_THRESH = 0.1 # meV/atom @@ -228,7 +228,7 @@ def make( if self.prev_calc_dir_argname is not None: bulk_kwargs[self.prev_calc_dir_argname] = prev_dir bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) - bulk.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + bulk.update_config({"manager_config": {"_fworker": "gpu_reg_fworker"}}) jobs.append(bulk) outputs.append(bulk.output) structure = bulk.output.structure @@ -249,11 +249,18 @@ def make( } ) - # 2. if supercell_matrix is None, supercell size will be determined after relax - # maker to ensure that cell lengths are really larger than threshold. - # then, perturbations will be generated based on the supercell size. - # STATIC calculations will then be run on the perturbed structures, and the - # forces and perturbed structures will be aggregated. + # # 2. if supercell_matrix is None, supercell size will be determined after relax + # # maker to ensure that cell lengths are really larger than threshold. + # # then, perturbations will be generated based on the supercell size. + # # STATIC calculations will then be run on the perturbed structures, and the + # # forces and perturbed structures will be aggregated. + # from pymatgen.core.structure import IStructure + # file_path = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-03-30-01-25-24-547097/launcher_2024-03-31-20-40-35-997684/launcher_2024-03-31-20-43-19-730610/CONTCAR" + # prev_dir = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/bulk_relax_2_4622_96_VASP/launcher_2024-03-31-20-43-19-730610/CONTCAR" + # structure = IStructure.from_file(file_path) + # prev_dir = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-03-30-01-25-24-547097/launcher_2024-03-31-20-40-35-997684/launcher_2024-03-31-20-43-19-730610" + # prev_dir = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/tight_relax2_4622_VASP_6_phonopy_exact/launcher_2024-04-14-02-09-48-767506" + # prev_dir = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-04-10-16-08-15-920932/launcher_2024-04-14-01-38-32-712261/launcher_2024-04-14-02-09-48-767506" static_calcs = hiphive_static_calcs( structure=structure, supercell_matrix=supercell_matrix, @@ -268,52 +275,78 @@ def make( ) jobs.append(static_calcs) - - # 3. Hiphive Fitting of FCPs upto 4th order - if n_structures >= 10: - fit_force_constant = run_hiphive_individually( - mpid = mpid, - cutoffs = cutoffs, - fit_method = fit_method, - disp_cut = disp_cut, - bulk_modulus = bulk_modulus, - temperature_qha = temperature_qha, - imaginary_tol = imaginary_tol, - prev_dir_json_saver = static_calcs.output["current_dir"], - # prev_dir_json_saver = prev_dir_json_saver, - loop = loops, - ) - else: - fit_force_constant = run_hiphive( - fit_method=fit_method, - disp_cut=disp_cut, - bulk_modulus=bulk_modulus, - temperature_qha=temperature_qha, - mesh_density=mesh_density, - imaginary_tol=imaginary_tol, - prev_dir_json_saver=static_calcs.output["current_dir"], - # prev_dir_json_saver=prev_dir_json_saver, - loop=loops, - cutoffs=cutoffs - ) - fit_force_constant.name += f" {loops}" - fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) - jobs.append(fit_force_constant) - outputs.append(fit_force_constant.output) - fit_force_constant.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"fit_force_constant_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) + # # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_6350" + # # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_disp_sampling_updated/launcher_2024-04-07-12-23-44-934527" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_1/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_2/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_4/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_5/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_50_disp_sampling_updated_only_harmonic/launcher_2024-04-08-10-11-41-526328" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_6/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_7/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/pscratch/sd/h/hrushi99/atomate2/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/Rohith_Mo2Zr_6715_VASP/launcher_2024-04-10-03-20-44-201391/launcher_2024-04-10-03-24-30-451833" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_5_less_than_0,01/launcher_2024-04-11-11-54-49-707881" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_8/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_9/launcher_2024-04-01-08-09-15-480067" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_6_phonopy_INCAR/launcher_2024-04-13-21-57-43-344701" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_10/launcher_2024-04-01-08-09-15-480067" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_6_phonopy_exact_INCAR/launcher_2024-04-14-19-15-11-794304" + # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_12_phonopy_exact_INCAR/launcher_2024-04-15-17-38-04-444124" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_41_VASP/launcher_2024-04-05-15-03-16-098315" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1565_VASP" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_131_VASP/launcher_2024-04-06-01-48-46-676261" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1265_VASP/launcher_2024-03-23-20-19-18-215601" + # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_2605_VASP" + # # 3. Hiphive Fitting of FCPs upto 4th order + # if n_structures >= 10: + # fit_force_constant = run_hiphive_individually( + # mpid = mpid, + # cutoffs = cutoffs, + # fit_method = fit_method, + # disp_cut = disp_cut, + # bulk_modulus = bulk_modulus, + # temperature_qha = temperature_qha, + # imaginary_tol = imaginary_tol, + # # prev_dir_json_saver = static_calcs.output["current_dir"], + # prev_dir_json_saver = prev_dir_json_saver, + # loop = loops, + # ) + # else: + # fit_force_constant = run_hiphive( + # fit_method=fit_method, + # disp_cut=disp_cut, + # bulk_modulus=bulk_modulus, + # temperature_qha=temperature_qha, + # # mesh_density=mesh_density, + # imaginary_tol=imaginary_tol, + # # prev_dir_json_saver=static_calcs.output["current_dir"], + # prev_dir_json_saver=prev_dir_json_saver, + # loop=loops, + # cutoffs=cutoffs + # ) + # fit_force_constant.name += f" {loops}" + # fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) + # jobs.append(fit_force_constant) + # outputs.append(fit_force_constant.output) + # fit_force_constant.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"bulk_modulus={bulk_modulus}", + # f"cutoffs={cutoffs}", + # f"fit_force_constant_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) #TODO: Implement Quality Control Job # # # # # # 7. Quality Control Job to check if the desired Test RMSE is achieved, @@ -356,172 +389,195 @@ def make( # # # # # } # # # # # ) - # 4. Perform phonon renormalization to obtain temperature-dependent - # force constants using hiPhive - outputs_renorm = [] - if renormalize: - for temperature in renormalize_temperature: - nconfig = renormalize_nconfig * (1 + temperature // 100) - renormalization = run_hiphive_renormalization( - temperature=temperature, - renorm_method=renormalize_method, - nconfig=nconfig, - renorm_TE_iter=renormalize_thermal_expansion_iter, - bulk_modulus=bulk_modulus, - prev_dir_hiphive=fit_force_constant.output["current_dir"], - # prev_dir_hiphive=prev_dir_hiphive, - loop=loops, - ) - renormalization.name += f" {temperature} {loops}" - renormalization.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) - jobs.append(renormalization) - outputs_renorm.append(renormalization.output) - outputs.append(renormalization.output) - renormalization.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"run_renormalization_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) - - - # 5. Extract Phonon Band structure & DOS from FC - # for 0K - fc_pdos_pb_to_db = run_fc_to_pdos( - renormalized=renormalize, - mesh_density=mesh_density, - prev_dir_json_saver=fit_force_constant.output["current_dir"], - # prev_dir_json_saver=prev_dir_hiphive, - loop=loops, - ) - fc_pdos_pb_to_db.name += f" {loops} 0K" - jobs.append(fc_pdos_pb_to_db) - outputs.append(fc_pdos_pb_to_db.output) - fc_pdos_pb_to_db.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"temperature=0K" - f"fc_pdos_pb_to_db_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) - # for finite temperatures - if renormalize: - for i, temperature in enumerate(renormalize_temperature): - fc_pdos_pb_to_db = run_fc_to_pdos( - renormalized=renormalize, - mesh_density=mesh_density, - prev_dir_json_saver=outputs_renorm[i][0], - loop=loops, - ) - fc_pdos_pb_to_db.name += f" {loops} {temperature}K" - jobs.append(fc_pdos_pb_to_db) - outputs.append(fc_pdos_pb_to_db.output) - fc_pdos_pb_to_db.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"temperature={temperature}K" - f"fc_pdos_pb_to_db_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) - - - # 6. Lattice thermal conductivity calculation using Sheng BTE - if calculate_lattice_thermal_conductivity: - if renormalize: - temperatures = renormalize_temperature - else: - temperatures = thermal_conductivity_temperature - # Because of the way ShengBTE works, a temperature array that is not - # evenly spaced out (T_step) requires submission for each temperature - if not renormalize: - if isinstance(temperatures, dict): - pass - elif type(temperatures) in [list, np.ndarray]: - assert all(np.diff(temperatures) == np.diff(temperatures)[0]) - lattice_thermal_conductivity = run_lattice_thermal_conductivity( - renormalized=renormalize, - temperature=temperatures, - loop=loops, - prev_dir_hiphive=fit_force_constant.output["current_dir"], - therm_cond_solver= self.THERM_COND_SOLVER - ) - lattice_thermal_conductivity.name += f" {temperatures} {loops}" - lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) - jobs.append(lattice_thermal_conductivity) - outputs.append(lattice_thermal_conductivity.output) - lattice_thermal_conductivity.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"run_lattice_thermal_conductivity_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) - else: - for t, T in enumerate(temperatures): - if T == 0: - continue - lattice_thermal_conductivity = run_lattice_thermal_conductivity( - renormalized=renormalize, - temperature=T, - loop=loops, - prev_dir_hiphive=outputs_renorm[t][0], - therm_cond_solver= self.THERM_COND_SOLVER - ) - - lattice_thermal_conductivity.name += f" {T} {loops}" - lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) - jobs.append(lattice_thermal_conductivity) - outputs.append(lattice_thermal_conductivity.output) - lattice_thermal_conductivity.metadata.update( - { - "tag": [ - f"mp_id={mpid}", - f"bulk_modulus={bulk_modulus}", - f"run_lattice_thermal_conductivity_{loops}", - f"nConfigsPerStd={n_structures}", - f"fixedDispls={fixed_displs}", - f"dispCut={disp_cut}", - f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", - f"supercell_matrix={supercell_matrix}", - f"loop={loops}", - ] - } - ) + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1479_displ_sampling" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-59-31-454856-98976" # fix 3 + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-05-05-27-911205-37950" # fix 5 + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-10-19-33-51-707117-99373" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-10-22-08-52-739207-65556" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-03-39-27-546384-72231" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-05-08-16-334659-40420" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-20-30-11-696345-68541" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-20-19-30-722449-33890" + # # 4. Perform phonon renormalization to obtain temperature-dependent + # # force constants using hiPhive + # outputs_renorm = [] + # if renormalize: + # for temperature in renormalize_temperature: + # nconfig = renormalize_nconfig * (1 + temperature // 100) + # renormalization = run_hiphive_renormalization( + # temperature=temperature, + # renorm_method=renormalize_method, + # nconfig=nconfig, + # renorm_TE_iter=renormalize_thermal_expansion_iter, + # bulk_modulus=bulk_modulus, + # # prev_dir_hiphive=fit_force_constant.output["current_dir"], + # prev_dir_hiphive=prev_dir_hiphive, + # loop=loops, + # ) + # renormalization.name += f" {temperature} {loops}" + # renormalization.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) + # jobs.append(renormalization) + # outputs_renorm.append(renormalization.output) + # outputs.append(renormalization.output) + # renormalization.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"bulk_modulus={bulk_modulus}", + # f"run_renormalization_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) + + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-06-04-29-02-534935-41708" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-16-14-51-625183-73398" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-18-38-06-547834-98127" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-19-14-17-420858-98262" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-24-38-567687-98068" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-59-31-454856-98976" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-22-03-35-538120-20851" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-05-05-27-911205-37950" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-14-52-07-269546-30380" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-16-51-44-504507-24744" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-18-24-57-922525-28171" + # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-15-01-03-47-788448-49252" + # # 5. Extract Phonon Band structure & DOS from FC + # # for 0K + # fc_pdos_pb_to_db = run_fc_to_pdos( + # renormalized=renormalize, + # mesh_density=mesh_density, + # prev_dir_json_saver=fit_force_constant.output["current_dir"], + # # prev_dir_json_saver=prev_dir_hiphive, + # loop=loops, + # ) + # fc_pdos_pb_to_db.name += f" {loops} 0K" + # jobs.append(fc_pdos_pb_to_db) + # outputs.append(fc_pdos_pb_to_db.output) + # fc_pdos_pb_to_db.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"cutoffs={cutoffs}", + # f"bulk_modulus={bulk_modulus}", + # f"temperature=0K" + # f"fc_pdos_pb_to_db_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) + # # for finite temperatures + # if renormalize: + # for i, temperature in enumerate(renormalize_temperature): + # fc_pdos_pb_to_db = run_fc_to_pdos( + # renormalized=renormalize, + # mesh_density=mesh_density, + # prev_dir_json_saver=outputs_renorm[i][0], + # loop=loops, + # ) + # fc_pdos_pb_to_db.name += f" {loops} {temperature}K" + # jobs.append(fc_pdos_pb_to_db) + # outputs.append(fc_pdos_pb_to_db.output) + # fc_pdos_pb_to_db.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"bulk_modulus={bulk_modulus}", + # f"temperature={temperature}K" + # f"fc_pdos_pb_to_db_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) + + + # # 6. Lattice thermal conductivity calculation using Sheng BTE + # if calculate_lattice_thermal_conductivity: + # if renormalize: + # temperatures = renormalize_temperature + # else: + # temperatures = thermal_conductivity_temperature + # # Because of the way ShengBTE works, a temperature array that is not + # # evenly spaced out (T_step) requires submission for each temperature + # if not renormalize: + # if isinstance(temperatures, dict): + # pass + # elif type(temperatures) in [list, np.ndarray]: + # assert all(np.diff(temperatures) == np.diff(temperatures)[0]) + # lattice_thermal_conductivity = run_lattice_thermal_conductivity( + # renormalized=renormalize, + # temperature=temperatures, + # loop=loops, + # prev_dir_hiphive=fit_force_constant.output["current_dir"], + # therm_cond_solver= self.THERM_COND_SOLVER + # ) + # lattice_thermal_conductivity.name += f" {temperatures} {loops}" + # lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + # jobs.append(lattice_thermal_conductivity) + # outputs.append(lattice_thermal_conductivity.output) + # lattice_thermal_conductivity.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"bulk_modulus={bulk_modulus}", + # f"run_lattice_thermal_conductivity_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) + # else: + # for t, T in enumerate(temperatures): + # if T == 0: + # continue + # lattice_thermal_conductivity = run_lattice_thermal_conductivity( + # renormalized=renormalize, + # temperature=T, + # loop=loops, + # prev_dir_hiphive=outputs_renorm[t][0], + # therm_cond_solver= self.THERM_COND_SOLVER + # ) + + # lattice_thermal_conductivity.name += f" {T} {loops}" + # lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + # jobs.append(lattice_thermal_conductivity) + # outputs.append(lattice_thermal_conductivity.output) + # lattice_thermal_conductivity.metadata.update( + # { + # "tag": [ + # f"mp_id={mpid}", + # f"bulk_modulus={bulk_modulus}", + # f"run_lattice_thermal_conductivity_{loops}", + # f"nConfigsPerStd={n_structures}", + # f"fixedDispls={fixed_displs}", + # f"dispCut={disp_cut}", + # f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + # f"supercell_matrix={supercell_matrix}", + # f"loop={loops}", + # ] + # } + # ) return Flow(jobs=jobs, output=outputs, name=f"{mpid}_{self.THERM_COND_SOLVER}_" - f"{fixed_displs}_" + f"{disp_cut}_" + f"{cutoffs}_" + f"seperate_fit_24_harm_0.03_12_anharm_0.08" f"{self.name}") @property diff --git a/src/atomate2/common/jobs/elastic.py b/src/atomate2/common/jobs/elastic.py index f4606e78e7..88d5898ca1 100644 --- a/src/atomate2/common/jobs/elastic.py +++ b/src/atomate2/common/jobs/elastic.py @@ -146,6 +146,7 @@ def run_elastic_deformations( elastic_job_kwargs[prev_dir_argname] = prev_dir # create the job relax_job = elastic_relax_maker.make(deformed_structure, **elastic_job_kwargs) + relax_job.update_config({"manager_config": {"_fworker": "gpu_reg_fworker"}}) #TODO remove this line relax_job.append_name(f" {idx + 1}/{len(deformations)}") relaxations.append(relax_job) diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index 864b5cf70d..0f6aeb807f 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -6,6 +6,7 @@ from __future__ import annotations import json +import math import os import shlex import subprocess @@ -175,6 +176,19 @@ def hiphive_static_calcs( return Response(replace=jobs, output=outputs) + +def find_smallest_x(n_structures: int) -> int: + # Starting from x = 3 because x - 2 in the denominator implies x > 2 + x = 3 + while True: + # Calculate the value of the expression for the current x + expression_value = (n_structures - 4) / (x - 2) + + # Check if the expression value satisfies the inequality + if 5 <= expression_value <= 6: + return math.floor(x) # Explicitly return the floor of x + x += 0.1 + @job def generate_hiphive_displacements( structure: Structure, @@ -185,32 +199,91 @@ def generate_hiphive_displacements( ) -> list[Structure]: """Job generates the perturbed structures for hiPhive fitting.""" if fixed_displs is None: + # # fixed_displs = [0.01, 0.03, 0.08, 0.1] #TODO, this list is used in + # # the paper + # smallest_disp = 0.01 - # fixed_displs = [0.01, 0.03, 0.08, 0.1] #TODO, this list is used in - # the paper - smallest_disp = 0.01 + # # dictionary of largest displacements for each period + # largest_displacement = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} - # dictionary of largest displacements for each period - largest_displacement = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} + # row = int( + # np.around(np.array([s.row for s in structure.species]).mean(), 0) + # ) - row = int( - np.around(np.array([s.row for s in structure.species]).mean(), 0) - ) + # largest_disp = largest_displacement[row] + + # ratio = int(largest_disp/smallest_disp) + # if 60 >= n_structures >= 4: + # logger.info(f"n_structures inside if statement >= 4 is {n_structures}") + # factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) + # elif n_structures < 4 : + # factors = np.sqrt(np.linspace(1,ratio**2,4)) + # else: + # factors = np.sqrt(np.linspace(1,ratio**2,60)) + # amplitudes = (smallest_disp*factors).round(3) + # logger.info(f"amplitudes = {amplitudes}") + # fixed_displs = amplitudes.tolist() + # logger.info(f"list_amplitudes = {fixed_displs}") + + if n_structures < 4: + n_structures = 4 + else: + pass + logger.info(f"n_structures = {n_structures}") + # # Ask the user for n_structures + # n_structures = int(input("Enter the value of n_structures: ")) - largest_disp = largest_displacement[row] + # # Find the smallest integer x that satisfies the inequality + # smallest_x = find_smallest_x(n_structures) + + smallest_x = 2 + + # Output the result + logger.info(f"The smallest integer x that satisfies the inequality is: {smallest_x}") + + # Create structures at equal interval between disp = 0.01 & 0.055 + # equal_interval_structures = np.linspace(0.01, 0.055, smallest_x) + equal_interval_structures = np.linspace(0.001, 0.005, smallest_x) + # create a new np array with 3 elements equaling 0.03 and 0.08 respectively + equal_interval_structures = np.array([0.03, 0.08]) + + # # Remaining number of structures + # remaining_structures = n_structures - smallest_x + + + # # structure_data = loadfn(f"/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96/launcher_2024-04-01-08-09-15-480067/structure_data_1.json") + # # structure = structure_data["structure"] + + # largest_period_default = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} + + # row = int( + # np.around(np.array([s.row for s in structure.species]).mean(), 0) + # ) + # logger.info(f"row = {row}") + + # largest = largest_period_default[row] + # logger.info(f"largest = {largest}") + # smallest = 0.056 + + # ratio = largest/smallest + # logger.info(f"ratio = {ratio}") + # logger.info(f"smallest = {smallest}") + # logger.info(f"largest = {largest}") + # logger.info(f"remaining_structures = {remaining_structures}") + # factors = np.sqrt(np.linspace(1,ratio**2,remaining_structures)) + # logger.info(f"factors = {factors}") + # amplitudes = (smallest*factors).round(3) + # logger.info(f"amplitudes = {amplitudes}") + + # # Combine the two sets of structures + # all_amplitudes = np.concatenate((equal_interval_structures, amplitudes)) + + # fixed_displs = all_amplitudes.tolist() + # logger.info(f"all_amplitudes = {all_amplitudes}") + + fixed_displs = equal_interval_structures.tolist() + logger.info(f"all_amplitudes = {fixed_displs}") - ratio = int(largest_disp/smallest_disp) - if 60 >= n_structures >= 4: - logger.info(f"n_structures inside if statement >= 4 is {n_structures}") - factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) - elif n_structures < 4 : - factors = np.sqrt(np.linspace(1,ratio**2,4)) - else: - factors = np.sqrt(np.linspace(1,ratio**2,60)) - amplitudes = (smallest_disp*factors).round(3) - logger.info(f"amplitudes = {amplitudes}") - fixed_displs = amplitudes.tolist() - logger.info(f"list_amplitudes = {fixed_displs}") logger.info(f"supercell_matrix = {supercell_matrix}") supercell_structure = SupercellTransformation( @@ -225,7 +298,7 @@ def generate_hiphive_displacements( for i in range(len(fixed_displs)): supercell_ase = AseAtomsAdaptor.get_atoms(supercell_structure) structures_ase = generate_displaced_structures( - atoms=supercell_ase, n_structures=1, distance=fixed_displs[i], loop=loop + atoms=supercell_ase, n_structures=n_structures, distance=fixed_displs[i], loop=loop ) structures_ase_all.extend(structures_ase) @@ -460,7 +533,13 @@ def run_hiphive_individually( supercell_structure = structure_data["supercell_structure"] if cutoffs is None: - cutoffs = get_cutoffs(supercell_structure) + # cutoffs = get_cutoffs(supercell_structure) + # cutoffs = [[8.0]] + cutoffs = [[4, 3, 2.66]] + cutoffs = [[5, 3.5]] + cutoffs = [[3, 2.5]] + cutoffs = [[5, 2.5]] + # cutoffs = [[5, 2.5]] logger.info(f"cutoffs is {cutoffs}") else: pass @@ -579,10 +658,15 @@ def collect_hiphive_outputs( fitting_data["best_cutoff"] = result["fitting_data"]["cutoffs"][0] fitting_data["best_rmse"] = result["fitting_data"]["rmse_test"][0] best_fit["rmse_test"] = result["fitting_data"]["rmse_test"][0] + # following line is commented out only for testing purposes fitting_data["n_imaginary"] = result["fitting_data"]["n_imaginary"] + # following line is commented out only for testing purposes ends copy_hiphive_outputs(best_fit["directory"]) + # following line is commented out only for testing purposes thermal_data = loadfn("thermal_data.json") + # thermal_data = {} + # following line is commented out only for testing purposes ends dumpfn(fitting_data, "fitting_data.json") param = np.loadtxt("parameters.txt") @@ -639,7 +723,39 @@ def run_hiphive( if cutoffs is None: cutoffs = get_cutoffs(supercell_structure) - cutoffs = [[6, 3.5, 3]] + logger.info(f"cutoffs is {cutoffs}") + # cutoffs = [[6, 3.5, 3]] + # cutoffs = [[9.0, 5.5, 4.25]] + # cutoffs = [[8.0]] + # cutoffs = [[9.0, 6.5, 5]] + # cutoffs = [[9.0, 5.5, 3]] + # cutoffs = [[9, 9.25, 5]] + # cutoffs = [[9, 8.3125, 5]] + # cutoffs = [[4, 4, 4]] + # cutoffs = [[3, 2.5]] + # cutoffs = [[2, 2.7]] + # cutoffs = [[6, 4, 2]] + # cutoffs = [[1.8635, 4]] + # cutoffs = [[2.559, 4]] # results in 2nd, 3rd order -- 124, 7379 DOFs -- imaginary modes = 8 + # cutoffs = [[2.54, 4]] # results in 2nd, 3rd order -- 106, 7379 DOFs -- imaginary modes = 2 + # cutoffs = [[2.5, 4]] # results in 2nd, 3rd order -- 88, 7379 DOFs -- imaginary modes = 0 + # cutoffs = [[2.52, 4]] # results in 2nd, 3rd order -- 97, 7379 DOFs -- imaginary modes = 2 + # cutoffs = [[2.52, 3]] # results in 2nd, 3rd order -- 97, 857 DOFs -- imaginary modes = 0 + # cutoffs = [[2.52, 4.6]] # results in 2nd, 3rd order -- 97, 20804 DOFs -- imaginary modes = 2 + # cutoffs = [[2.5, 4.6]] # results in 2nd, 3rd order -- 88, 20804 DOFs -- imaginary modes = 0 + # cutoffs = [[2.5, 5.6]] # results in 2nd, 3rd order -- 88, 67511 DOFs -- imaginary modes = ? + cutoffs = [[9.0, 5.5, 3.5]] + cutoffs = [[9, 8.3125, 5]] + cutoffs = [[6, 5, 4]] # results in {2: 36, 3: 263, 4: 327} DOFs -- imaginary modes = 18, rmse = 132meV/atom + cutoffs = [[8, 5, 4]] # results in {2: 79, 3: 263, 4: 327} DOFs -- imaginary modes = 174, rmse = 68.8meV/atom + cutoffs = [[10, 5, 4]] # results in {2: 151, 3: 263, 4: 327} DOFs -- imaginary modes = 150, rmse = 40.5meV/atom + cutoffs = [[10, 7, 4]] # results in {2: 151, 3: 2149, 4: 327} DOFs -- imaginary modes = 150, rmse = 40.1meV/atom + cutoffs = [[10, 7, 5]] # results in {2: 151, 3: 2149, 4: 2876} DOFs -- imaginary modes = 150, rmse = 38.03meV/atom + cutoffs = [[9, 6, 4]] # {2: 102, 3: 775, 4: 428} -- imaginary modes = 114, rmse = 28.9meV/atom -- MgO 1265 + cutoffs = [[9, 7, 4.3]] # {2: 102, 3: 2593, 4: 806} -- imaginary modes = 114, rmse = 29.2meV/atom -- MgO 1265 + cutoffs = [[9, 6, 4.5]] # {2: 74, 3: 597, 4: 428} -- imaginary modes = 0, rmse = 13.3meV/atom -- CaO 2605 + cutoffs = [[10, 6, 4.5]] # {2: 88, 3: 597, 4: 428} -- imaginary modes = 0, rmse = 9.6meV/atom -- CaO 2605 + cutoffs = [[10, 7, 5]] # {2: 88, 3: 775, 4: 806} -- imaginary modes = 0, rmse = 9.6meV/atom -- CaO 2605 logger.info(f"cutoffs is {cutoffs}") else: pass @@ -730,26 +846,31 @@ def run_hiphive( else: logger.info("cs is not an instance of ClusterSpace") - if fitting_data["n_imaginary"] == 0: - logger.info("No imaginary modes! Writing ShengBTE files") - atoms = AseAtomsAdaptor.get_atoms(parent_structure) - # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) - fcs.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") - ForceConstants.write_to_phonopy(fcs, "fc2.hdf5", "hdf5") - ForceConstants.write_to_phono3py(fcs, "fc3.hdf5", "hdf5") - ### detour from hdf5 - supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) - FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) - FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) - else: - logger.info("ShengBTE files not written due to imaginary modes.") - logger.info("You may want to perform phonon renormalization.") + # # following code is commented only for testing puropose + # if fitting_data["n_imaginary"] == 0: + # # following code is commented only for testing puropose ends + # # if True: # change this back to above if statement + # logger.info("No imaginary modes! Writing ShengBTE files") + # atoms = AseAtomsAdaptor.get_atoms(parent_structure) + # ### fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) + # fcs.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") + # # following code is commented only for testing puropose + # ForceConstants.write_to_phonopy(fcs, "fc2.hdf5", "hdf5") + # ForceConstants.write_to_phono3py(fcs, "fc3.hdf5", "hdf5") + # ### detour from hdf5 + # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + # FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + # FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) + # # following code is commented only for testing puropose ends + # else: + # logger.info("ShengBTE files not written due to imaginary modes.") + # logger.info("You may want to perform phonon renormalization.") current_dir = os.getcwd() outputs: dict[str, list] = { - "thermal_data": thermal_data, - "anharmonic_data": anharmonic_data, + "thermal_data": thermal_data, # replace with "thermal_data" + "anharmonic_data": thermal_data, # replace with "anharmonic_data" "fitting_data": fitting_data, "param": param, "current_dir": current_dir @@ -798,25 +919,65 @@ def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: ------- A list of trial cutoffs. """ - # indexed as min_cutoffs[order][period] - # DO NOT CHANGE unless you know what you are doing + # # indexed as min_cutoffs[order][period] + # # DO NOT CHANGE unless you know what you are doing + # min_cutoffs = { + # 2: {1: 5.0, 2: 6.0, 3: 7.0, 4: 8.0, 5: 9.0, 6: 10.0, 7: 11.0}, + # 3: {1: 3.0, 2: 3.5, 3: 4.5, 4: 5.5, 5: 6.0, 6: 6.5, 7: 7.0}, + # 4: {1: 2.5, 2: 3.0, 3: 3.5, 4: 4.0, 5: 4.5, 6: 5.0, 7: 5.5}, + # } + # inc = {2: 2, 3: 1.5, 4: 0.6} + # steps = {2: 1, 3: 0.75, 4: 0.6} + + # row = int( + # np.around(np.array([s.row for s in supercell_structure.species]).mean(), 0) + # ) + # factor = row / 4 + # mins = {2: min_cutoffs[2][row], 3: min_cutoffs[3][row], 4: min_cutoffs[4][row]} + + # range_two = np.arange( + # mins[2], mins[2] + factor * (inc[2] + steps[2]), factor * steps[2] + # ) + # range_three = np.arange( + # mins[3], mins[3] + factor * (inc[3] + steps[3]), factor * steps[3] + # ) + # range_four = np.arange( + # mins[4], mins[4] + factor * (inc[4] + steps[4]), factor * steps[4] + # ) + + # cutoffs = np.array(list(map(list, product(range_two, range_three, range_four)))) + # max_cutoff = estimate_maximum_cutoff(AseAtomsAdaptor.get_atoms(supercell_structure)) + # cutoffs[cutoffs > max_cutoff] = max_cutoff + # logger.info(f"CUTOFFS \n {cutoffs}") + # logger.info(f"MAX_CUTOFF \n {max_cutoff}") + # good_cutoffs = np.all(cutoffs < max_cutoff - 0.1, axis=1) + # logger.info(f"GOOD CUTOFFS \n{good_cutoffs}") + + + min_cutoffs = { + # 2: {1: 7.0, 2: 8.0, 3: 9.0, 4: 10.0, 5: 11.0, 6: 12.0, 7: 13.0}, 2: {1: 5.0, 2: 6.0, 3: 7.0, 4: 8.0, 5: 9.0, 6: 10.0, 7: 11.0}, - 3: {1: 3.0, 2: 3.5, 3: 4.5, 4: 5.5, 5: 6.0, 6: 6.5, 7: 7.0}, - 4: {1: 2.5, 2: 3.0, 3: 3.5, 4: 4.0, 5: 4.5, 6: 5.0, 7: 5.5}, + # 2: {1: 6.0, 2: 7.0, 3: 8.0, 4: 9.0, 5: 10.0, 6: 11.0, 7: 12.0}, + 3: {1: 2.5, 2: 3.0, 3: 4.0, 4: 5.0, 5: 5.5, 6: 6.0, 7: 6.5}, + # 3: {1: 3.0, 2: 3.5, 3: 4.5, 4: 5.5, 5: 6.0, 6: 6.5, 7: 7.0}, + 4: {1: 1.5, 2: 2.0, 3: 2.5, 4: 3.0, 5: 3.5, 6: 4.0, 7: 4.5}, + # 4: {1: 2.5, 2: 3.0, 3: 3.5, 4: 4.0, 5: 4.5, 6: 5.0, 7: 5.5}, } - inc = {2: 2, 3: 1.5, 4: 0.6} - steps = {2: 1, 3: 0.75, 4: 0.6} + inc = {2: 0, 3: 2.0, 4: 1.2} + steps = {2: 0, 3: 0.75, 4: 0.6} row = int( np.around(np.array([s.row for s in supercell_structure.species]).mean(), 0) ) + logger.info(f"row = {row}") factor = row / 4 + logger.info(f"factor = {factor}") + mins = {2: min_cutoffs[2][row], 3: min_cutoffs[3][row], 4: min_cutoffs[4][row]} - range_two = np.arange( - mins[2], mins[2] + factor * (inc[2] + steps[2]), factor * steps[2] - ) + # create an NDArray of 2nd order cutofss with only one entry -> mins[2], and 3rd and 4th order cutoffs with a range + range_two = np.array([mins[2]]) range_three = np.arange( mins[3], mins[3] + factor * (inc[3] + steps[3]), factor * steps[3] ) @@ -825,12 +986,18 @@ def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: ) cutoffs = np.array(list(map(list, product(range_two, range_three, range_four)))) + logger.info(f"cutoffs = {cutoffs}") + max_cutoff = estimate_maximum_cutoff(AseAtomsAdaptor.get_atoms(supercell_structure)) cutoffs[cutoffs > max_cutoff] = max_cutoff logger.info(f"CUTOFFS \n {cutoffs}") logger.info(f"MAX_CUTOFF \n {max_cutoff}") good_cutoffs = np.all(cutoffs < max_cutoff - 0.1, axis=1) logger.info(f"GOOD CUTOFFS \n{good_cutoffs}") + + cutoffs_used = cutoffs[good_cutoffs].tolist() + logger.info(f"cutoffs_used = {cutoffs_used}") + return cutoffs[good_cutoffs].tolist() @@ -980,12 +1147,13 @@ def _run_cutoffs( cs = ClusterSpace(supercell_atoms, cutoffs, symprec=1e-3, acoustic_sum_rules=True) logger.debug(cs.__repr__()) - n2nd = cs.get_n_dofs_by_order(2) + n2nd = cs.get_n_dofs_by_order(2) # change it back to cs.get_n_dofs_by_order(2) nall = cs.n_dofs logger.info("Fitting harmonic force constants separately") separate_fit = True logger.info(f"disp_cut = {disp_cut}") + # commenting it out only for testing purposes sc = get_structure_container( cs, structures, separate_fit, disp_cut, ncut=n2nd, param2=None ) @@ -1014,8 +1182,10 @@ def _run_cutoffs( logger.info(f"omega = {omega}") imaginary = np.any(omega < -1e-3) logger.info(f"imaginary = {imaginary}") + # commenting it out only for testing purposes if imaginary: + # if True: # only for testing purposes logger.info( "Imaginary modes found! Fitting anharmonic force constants separately" ) @@ -1043,6 +1213,18 @@ def _run_cutoffs( parameters = opt.parameters logger.info(f"Training complete for cutoff: {i}, {cutoffs}") + # # only for test purpose. Remove this later + # logger.info("No imaginary modes! Fitting all force constants in one shot") + # separate_fit = False + # sc = get_structure_container( + # cs, structures, separate_fit, disp_cut=None, ncut=None, param2=None + # ) + # opt = Optimizer(sc.get_fit_data(), fit_method, [0, nall], **fit_kwargs) + # opt.train() + # parameters = opt.parameters + # logger.info(f"Training complete for cutoff: {i}, {cutoffs}") + # # only for test purpose. Remove this later + logger.info(f"parameters before enforcing sum rules {parameters}") logger.info(f"Memory use: {psutil.virtual_memory().percent} %") parameters = enforce_rotational_sum_rules(cs, parameters, ["Huang", "Born-Huang"]) @@ -1089,6 +1271,7 @@ def get_structure_container( A hiPhive StructureContainer. """ sc = StructureContainer(cs) + logger.info(f"initial shape of fit matrix = {sc.data_shape}") saved_structures = [] for _, structure in enumerate(structures): displacements = structure.get_array("displacements") @@ -1109,6 +1292,7 @@ def get_structure_container( sc.add_structure(structure) saved_structures.append(structure) + logger.info(f"final shape of fit matrix (total # of atoms in all added supercells, n_dofs) = (rows, columns) = {sc.data_shape}") logger.info("We have completed adding structures") logger.info(sc.get_fit_data()) @@ -1192,7 +1376,7 @@ def harmonic_properties( # entropy = entropy[0] # heat_capacity = heat_capacity[0] logger.info(f"Heat_capacity_harmonic_property[0]: {heat_capacity}") - temperature = temperature[0] + # temperature = temperature[0] return { "temperature": temperature, "free_energy": free_energy, @@ -1536,6 +1720,7 @@ def run_fc_to_pdos( logger.info("Finished inserting force constants and phonon data") else: + # following line is commented out only for testing purposes renorm_thermal_data = loadfn("thermal_data.json") # renorm_thermal_data.json fcs = ForceConstants.read("force_constants.fcs") T = renorm_thermal_data["temperature"] @@ -1551,6 +1736,8 @@ def run_fc_to_pdos( logger.info( f"Finished inserting renormalized force constants and phonon data at {T} K" ) + # following line is commented out only for testing purposes ends + # pass return uniform_bs, lm_bs, dos, prev_dir_json_saver @@ -1635,6 +1822,7 @@ def run_hiphive_renormalization( TD_structure_data["structure"] = parent_structure TD_structure_data["supercell_structure"] = supercell_structure + logger.info("Renormalization is now completeed") # Additional renormalization with thermal expansion - # optional - just single "iteration" for now if renorm_TE_iter: @@ -1663,27 +1851,32 @@ def run_hiphive_renormalization( TD_structure_data["structure"] = parent_structure_TD TD_structure_data["supercell_structure"] = supercell_structure_TD - # Thermodynamic integration for anharmonic free energy - TD_data = thermodynamic_integration_ifc( - TD_data, # everything TD - fcs, # original - param, # original - imag_modes_bool, # if False, only uses lambda=0 - ) + # # Thermodynamic integration for anharmonic free energy + # TD_data = thermodynamic_integration_ifc( + # TD_data, # everything TD + # fcs, # original + # param, # original + # imag_modes_bool, # if False, only uses lambda=0 + # ) # write results logger.info("Writing renormalized results") fcs_TD = TD_data['fcs'] fcs_TD.write("force_constants.fcs") # if "n_imaginary" in TD_data: if TD_data["n_imaginary"] != 0: + # thermal_keys = ["temperature","free_energy","entropy","heat_capacity", + # "free_energy_correction_S","free_energy_correction_SC", + # "free_energy_correction_TI"] thermal_keys = ["temperature","free_energy","entropy","heat_capacity", - "free_energy_correction_S","free_energy_correction_SC", - "free_energy_correction_TI"] + "free_energy_correction_S","free_energy_correction_SC"] else: + # thermal_keys = ["temperature","free_energy","entropy","heat_capacity", + # "gruneisen","thermal_expansion","expansion_fraction", + # "free_energy_correction_S","free_energy_correction_SC", + # "free_energy_correction_TI"] thermal_keys = ["temperature","free_energy","entropy","heat_capacity", "gruneisen","thermal_expansion","expansion_fraction", - "free_energy_correction_S","free_energy_correction_SC", - "free_energy_correction_TI"] + "free_energy_correction_S","free_energy_correction_SC"] TD_thermal_data = {key: [] for key in thermal_keys} for key in thermal_keys: TD_thermal_data[key].append(TD_data[key]) @@ -1695,28 +1888,40 @@ def run_hiphive_renormalization( else: logger.info("No imaginary modes! Writing ShengBTE files") - parent_atoms_TD = copy(parent_atoms) - logger.info(f"TD_data exp frac: {TD_data['expansion_fraction']}") - logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,0]}") - logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,1]}") - logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,2]}") - new_cell = Cell(np.transpose([parent_atoms_TD.get_cell()[:,i]*(1+TD_data["expansion_fraction"][0,i]) for i in range(3)])) - parent_atoms_TD.set_cell(new_cell,scale_atoms=True) + # parent_atoms_TD = copy(parent_atoms) + # logger.info(f"TD_data exp frac: {TD_data['expansion_fraction']}") + # logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,0]}") + # logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,1]}") + # logger.info(f"TD_data exp frac 0: {TD_data['expansion_fraction'][0,2]}") + # new_cell = Cell(np.transpose([parent_atoms_TD.get_cell()[:,i]*(1+TD_data["expansion_fraction"][0,i]) for i in range(3)])) + # parent_atoms_TD.set_cell(new_cell,scale_atoms=True) + + # prim_TD_phonopy = PhonopyAtoms(symbols=parent_atoms_TD.get_chemical_symbols(), + # scaled_positions=parent_atoms_TD.get_scaled_positions(), + # cell=parent_atoms_TD.cell) + # phonopy_TD = Phonopy(prim_TD_phonopy, supercell_matrix=supercell_matrix, primitive_matrix=None) + + # atoms = AseAtomsAdaptor.get_atoms(parent_structure_TD) + # # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) + # fcs_TD.write_to_phonopy(f"FORCE_CONSTANTS_2ND_{temperature}", format="text") + # ForceConstants.write_to_phonopy(fcs_TD, "fc2.hdf5", "hdf5") + # ForceConstants.write_to_phono3py(fcs_TD, "fc3.hdf5", "hdf5") + # ### detour from hdf5 + # supercell_atoms = phonopy_atoms_to_ase(phonopy_TD.supercell) + # FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + # FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD_{temperature}", atoms, order=3, fc_tol=1e-4) - prim_TD_phonopy = PhonopyAtoms(symbols=parent_atoms_TD.get_chemical_symbols(), - scaled_positions=parent_atoms_TD.get_scaled_positions(), - cell=parent_atoms_TD.cell) - phonopy_TD = Phonopy(prim_TD_phonopy, supercell_matrix=supercell_matrix, primitive_matrix=None) - atoms = AseAtomsAdaptor.get_atoms(parent_structure_TD) + atoms = AseAtomsAdaptor.get_atoms(parent_structure) + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) - fcs_TD.write_to_phonopy(f"FORCE_CONSTANTS_2ND_{temperature}", format="text") + fcs_TD.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") ForceConstants.write_to_phonopy(fcs_TD, "fc2.hdf5", "hdf5") ForceConstants.write_to_phono3py(fcs_TD, "fc3.hdf5", "hdf5") ### detour from hdf5 - supercell_atoms = phonopy_atoms_to_ase(phonopy_TD.supercell) + # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) - FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD_{temperature}", atoms, order=3, fc_tol=1e-4) + FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) dumpfn(TD_structure_data, "structure_data.json") dumpfn(TD_thermal_data, "thermal_data.json") @@ -1821,7 +2026,7 @@ def thermodynamic_integration_ifc( cs = TD_data['cs'] fcs_TD = TD_data["fcs"] param_TD = TD_data["param"] - T = TD_data['temperature'] + T = TD_data['temperature'][0] logger.info(f"Temperature = {T}") supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) renorm = Renormalization(cs, supercell_atoms, param, fcs, T, 'least_squares', 'rfe', param_TD, fcs_TD) diff --git a/src/atomate2/common/jobs/phonons.py b/src/atomate2/common/jobs/phonons.py index 6134601745..f3cadbf534 100644 --- a/src/atomate2/common/jobs/phonons.py +++ b/src/atomate2/common/jobs/phonons.py @@ -316,6 +316,7 @@ def run_phonon_displacements( phonon_job.update_maker_kwargs( {"_set": {"write_additional_data->phonon_info:json": info}}, dict_mod=True ) + phonon_job.update_config({"manager_config": {"_fworker": "gpu_reg_fworker"}}) # change to gpu_reg_fworker phonon_jobs.append(phonon_job) outputs["displacement_number"] = list(range(len(displacements))) outputs["uuids"] = [phonon_job.output.uuid] * len(displacements) @@ -328,6 +329,7 @@ def run_phonon_displacements( phonon_job = phonon_maker.make(displacement, prev_dir=prev_dir) else: phonon_job = phonon_maker.make(displacement) + phonon_job.update_config({"manager_config": {"_fworker": "gpu_reg_fworker"}}) # change to gpu_fworker phonon_job.append_name(f" {idx + 1}/{len(displacements)}") # we will add some meta data diff --git a/src/atomate2/vasp/flows/hiphive.py b/src/atomate2/vasp/flows/hiphive.py index 56e5c884bd..93dd1af11d 100644 --- a/src/atomate2/vasp/flows/hiphive.py +++ b/src/atomate2/vasp/flows/hiphive.py @@ -17,7 +17,7 @@ # Atomate2 packages from atomate2.vasp.jobs.core import StaticMaker, TightRelaxMaker from atomate2.vasp.jobs.phonons import PhononDisplacementMaker -from atomate2.vasp.sets.core import StaticSetGenerator +from atomate2.vasp.sets.core import StaticSetGenerator, TightRelaxSetGenerator if TYPE_CHECKING: from atomate2.vasp.jobs.base import BaseVaspMaker @@ -71,52 +71,101 @@ class HiphiveMaker(BaseHiphiveMaker): """ name: str = "Lattice-Dynamics-VASP" - static_energy_maker: BaseVaspMaker | None = field( - default_factory=lambda: StaticMaker( - input_set_generator=StaticSetGenerator(auto_ispin=True) - ) - ) + # static_energy_maker: BaseVaspMaker | None = field( + # default_factory=lambda: StaticMaker( + # input_set_generator=StaticSetGenerator(auto_ispin=True) + # ) + # ) bulk_relax_maker: DoubleRelaxMaker = field( - default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker()) + # default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker( + # input_set_generator=TightRelaxSetGenerator( + # user_incar_settings={ + # "ADDGRID": True, + # "ALGO": "Normal", + # "EDIFF": 1e-07, + # "SIGMA": 0.01, # changed from 0.1 + # # "EDIFFG": -1.000000e-08 + # } + # ) + # )) + default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker( + input_set_generator=TightRelaxSetGenerator( + user_incar_settings={ + "PREC": "Accurate", + "GGA": "PS", + "IBRION": 2, + "NSW": 20, + "NELMIN": 5, + "ISIF": 3, + "ENCUT": 648.744200, + "EDIFF": 1.000000e-06, + # "EDIFFG": -1.000000e-08, + "ISMEAR": 0, + "SIGMA": 1.000000e-02, + "IALGO": 38, + "LREAL": ".FALSE.", + "ADDGRID": ".TRUE.", + "LWAVE": ".FALSE.", + "LCHARG": ".FALSE.", + "NPAR": 4 + } + ) + )) ) phonon_displacement_maker: BaseVaspMaker | None = field( + # default_factory=lambda:PhononDisplacementMaker( + # input_set_generator = StaticSetGenerator( + # user_kpoints_settings={"reciprocal_density": 500}, + # user_incar_settings={ + # "ADDGRID": True, + # "ALGO": "Normal", + # "EDIFF": 1e-07, # changed from 1e-6 + # # "EDIFFG": -1.000000e-08, + # # "ENCUT": 600, + # "GGA": "PS", + # "IBRION": -1, + # "ISIF": 3, + # "ISMEAR": 0, + # "ISPIN": 2, + # "LAECHG": False, + # "LASPH": True, + # "LCHARG": False, + # "LORBIT": 11, + # "LREAL": "Auto", + # "LVHAR": False, + # "LVTOT": False, + # "LWAVE": False, + # # "MAGMOM": 250*0.6, + # "NCORE": 6, + # "NELM": 100, + # "NSW": 0, + # "PREC": "Accurate", + # "SIGMA": 0.01, # changed from 0.1 + # }, + # # auto_ispin=True, + # ) + # ) default_factory=lambda:PhononDisplacementMaker( input_set_generator = StaticSetGenerator( user_kpoints_settings={"reciprocal_density": 500}, user_incar_settings={ - # "IBRION": 2, - # "ISIF": 3, - # "ENCUT": 600, # Changed this from 600 - # "EDIFF": 1e-7, - # "LAECHG": False, - # "ALGO": "Normal", - # "NSW": 0, - # "LCHARG": False, - # "LREAL": True - - "ADDGRID": True, - "ALGO": "Normal", - "EDIFF": 1e-06, - "ENCUT": 600, + "PREC": "Accurate", + "SIGMA": 0.01, # changed from 0.1 + "PREC": "Accurate", "GGA": "PS", "IBRION": -1, - "ISIF": 3, + "NELMIN": 5, + "ENCUT": 648.744200, + "EDIFF": 1.000000e-06, "ISMEAR": 0, - "ISPIN": 2, - "LAECHG": False, - "LASPH": True, - "LCHARG": False, - "LORBIT": 11, - "LREAL": "Auto", - "LVHAR": False, - "LVTOT": False, - "LWAVE": False, - # "MAGMOM": 250*0.6, - "NCORE": 6, - "NELM": 100, - "NSW": 0, - "PREC": "Accurate", - "SIGMA": 0.1, + "SIGMA": 1.000000e-02, + "IALGO": 38, + "LREAL": ".FALSE.", + "ADDGRID": ".TRUE.", + "LWAVE": ".FALSE.", + "LCHARG": ".FALSE.", + "NPAR": 4, + }, # auto_ispin=True, ) From a18c0b71950835a87421367a9c41816aa4e7a1a6 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe Date: Mon, 27 May 2024 15:56:00 -0700 Subject: [PATCH 14/14] reorganizing code --- src/atomate2/common/flows/hiphive.py | 532 ++++++++++++--------------- src/atomate2/common/jobs/hiphive.py | 478 ++++++++++-------------- 2 files changed, 422 insertions(+), 588 deletions(-) diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index 5431bda9d0..fe7e48da22 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -127,7 +127,7 @@ class BaseHiphiveMaker(Maker, ABC): i * 100 for i in range(21) ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) T_RENORM: ClassVar[list[int]] = [ - 100, 500, 1000, 1500 # 300, 500, 600, 700, 800, 900, 1000, 1500, 2500, 2700, 3000 + 100 # 300, 500, 600, 700, 800, 900, 1000, 1500, 2500, 2700, 3000 ] # [i*100 for i in range(0,16)] # Temp. at which renorm. is to be performed # If renormalization is performed, # T_RENORM overrides T_KLAT for lattice thermal conductivity @@ -249,18 +249,11 @@ def make( } ) - # # 2. if supercell_matrix is None, supercell size will be determined after relax - # # maker to ensure that cell lengths are really larger than threshold. - # # then, perturbations will be generated based on the supercell size. - # # STATIC calculations will then be run on the perturbed structures, and the - # # forces and perturbed structures will be aggregated. - # from pymatgen.core.structure import IStructure - # file_path = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-03-30-01-25-24-547097/launcher_2024-03-31-20-40-35-997684/launcher_2024-03-31-20-43-19-730610/CONTCAR" - # prev_dir = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/bulk_relax_2_4622_96_VASP/launcher_2024-03-31-20-43-19-730610/CONTCAR" - # structure = IStructure.from_file(file_path) - # prev_dir = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-03-30-01-25-24-547097/launcher_2024-03-31-20-40-35-997684/launcher_2024-03-31-20-43-19-730610" - # prev_dir = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/tight_relax2_4622_VASP_6_phonopy_exact/launcher_2024-04-14-02-09-48-767506" - # prev_dir = "/pscratch/sd/h/hrushi99/atomate2/MgO_Zhuoying_LAC_NRC/block_2024-04-10-16-08-15-920932/launcher_2024-04-14-01-38-32-712261/launcher_2024-04-14-02-09-48-767506" + # 2. if supercell_matrix is None, supercell size will be determined after relax + # maker to ensure that cell lengths are really larger than threshold. + # then, perturbations will be generated based on the supercell size. + # STATIC calculations will then be run on the perturbed structures, and the + # forces and perturbed structures will be aggregated. static_calcs = hiphive_static_calcs( structure=structure, supercell_matrix=supercell_matrix, @@ -275,121 +268,56 @@ def make( ) jobs.append(static_calcs) - # # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_6350" - # # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_disp_sampling_updated/launcher_2024-04-07-12-23-44-934527" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_1/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_2/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_4/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_5/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_50_disp_sampling_updated_only_harmonic/launcher_2024-04-08-10-11-41-526328" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_6/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_7/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/pscratch/sd/h/hrushi99/atomate2/hiphive_4622_VASP_96_fix_3/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/Rohith_Mo2Zr_6715_VASP/launcher_2024-04-10-03-20-44-201391/launcher_2024-04-10-03-24-30-451833" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_5_less_than_0,01/launcher_2024-04-11-11-54-49-707881" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_8/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_9/launcher_2024-04-01-08-09-15-480067" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_6_phonopy_INCAR/launcher_2024-04-13-21-57-43-344701" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96_fix_10/launcher_2024-04-01-08-09-15-480067" - # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_6_phonopy_exact_INCAR/launcher_2024-04-14-19-15-11-794304" - # # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_12_phonopy_exact_INCAR/launcher_2024-04-15-17-38-04-444124" - # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_41_VASP/launcher_2024-04-05-15-03-16-098315" - # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1565_VASP" - # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_131_VASP/launcher_2024-04-06-01-48-46-676261" - # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1265_VASP/launcher_2024-03-23-20-19-18-215601" # prev_dir_json_saver = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_2605_VASP" - # # 3. Hiphive Fitting of FCPs upto 4th order - # if n_structures >= 10: - # fit_force_constant = run_hiphive_individually( - # mpid = mpid, - # cutoffs = cutoffs, - # fit_method = fit_method, - # disp_cut = disp_cut, - # bulk_modulus = bulk_modulus, - # temperature_qha = temperature_qha, - # imaginary_tol = imaginary_tol, - # # prev_dir_json_saver = static_calcs.output["current_dir"], - # prev_dir_json_saver = prev_dir_json_saver, - # loop = loops, - # ) - # else: - # fit_force_constant = run_hiphive( - # fit_method=fit_method, - # disp_cut=disp_cut, - # bulk_modulus=bulk_modulus, - # temperature_qha=temperature_qha, - # # mesh_density=mesh_density, - # imaginary_tol=imaginary_tol, - # # prev_dir_json_saver=static_calcs.output["current_dir"], - # prev_dir_json_saver=prev_dir_json_saver, - # loop=loops, - # cutoffs=cutoffs - # ) - # fit_force_constant.name += f" {loops}" - # fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) - # jobs.append(fit_force_constant) - # outputs.append(fit_force_constant.output) - # fit_force_constant.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"bulk_modulus={bulk_modulus}", - # f"cutoffs={cutoffs}", - # f"fit_force_constant_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) - - #TODO: Implement Quality Control Job - # # # # # # 7. Quality Control Job to check if the desired Test RMSE is achieved, - # # # # # # if not, then increase the number of structures -- - # # # # # # Using "addintion" feature of jobflow - # # # # # loops += 1 - # # # # # n_structures += 1 - # # # # # logger.info(f"Number of structures increased to {n_structures}") - # # # # # logger.info(f"loop = {loops}") - # # # # # error_check_job = quality_control( - # # # # # rmse_test=fit_force_constant.output[5], - # # # # # n_structures=n_structures, - # # # # # fixedDispls=fixed_displs, - # # # # # loop=loops, - # # # # # fit_method=fit_method, - # # # # # disp_cut=disp_cut, - # # # # # bulk_modulus=bulk_modulus, - # # # # # temperature_qha=temperature_qha, - # # # # # mesh_density=mesh_density, - # # # # # imaginary_tol=imaginary_tol, - # # # # # prev_dir_json_saver=static_calcs.output["current_dir"], - # # # # # prev_dir=prev_dir, - # # # # # supercell_matrix=supercell_matrix, - # # # # # # supercell_matrix_kwargs=supercell_matrix_kwargs, - # # # # # ) - # # # # # error_check_job.name += f" {loops}" - # # # # # jobs.append(error_check_job) - # # # # # outputs.append(error_check_job.output) - # # # # # error_check_job.metadata.update( - # # # # # { - # # # # # "tag": [ - # # # # # f"error_check_job_{loops}", - # # # # # f"nConfigsPerStd={n_structures}", - # # # # # f"fixedDispls={fixed_displs}", - # # # # # f"dispCut={disp_cut}", - # # # # # # f"supercell_matrix_kwargs={supercell_matrix_kwargs}", - # # # # # f"supercell_matrix={supercell_matrix}", - # # # # # f"loop={loops}", - # # # # # ] - # # # # # } - # # # # # ) - - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1479_displ_sampling" + # 3. Hiphive Fitting of FCPs upto 4th order + if n_structures >= 10: + fit_force_constant = run_hiphive_individually( + mpid = mpid, + cutoffs = cutoffs, + fit_method = fit_method, + disp_cut = disp_cut, + bulk_modulus = bulk_modulus, + temperature_qha = temperature_qha, + imaginary_tol = imaginary_tol, + prev_dir_json_saver = static_calcs.output["current_dir"], + # prev_dir_json_saver = prev_dir_json_saver, + loop = loops, + ) + else: + fit_force_constant = run_hiphive( + fit_method=fit_method, + disp_cut=disp_cut, + bulk_modulus=bulk_modulus, + temperature_qha=temperature_qha, + # mesh_density=mesh_density, + imaginary_tol=imaginary_tol, + prev_dir_json_saver=static_calcs.output["current_dir"], + # prev_dir_json_saver=prev_dir_json_saver, + loop=loops, + cutoffs=cutoffs + ) + fit_force_constant.name += f" {loops}" + fit_force_constant.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) + jobs.append(fit_force_constant) + outputs.append(fit_force_constant.output) + fit_force_constant.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"cutoffs={cutoffs}", + f"fit_force_constant_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + + + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_1479_displ_sampling" # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-59-31-454856-98976" # fix 3 # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-05-05-27-911205-37950" # fix 5 # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-10-19-33-51-707117-99373" @@ -398,181 +326,181 @@ def make( # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-05-08-16-334659-40420" # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-20-30-11-696345-68541" # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-16-20-19-30-722449-33890" - # # 4. Perform phonon renormalization to obtain temperature-dependent - # # force constants using hiPhive - # outputs_renorm = [] - # if renormalize: - # for temperature in renormalize_temperature: - # nconfig = renormalize_nconfig * (1 + temperature // 100) - # renormalization = run_hiphive_renormalization( - # temperature=temperature, - # renorm_method=renormalize_method, - # nconfig=nconfig, - # renorm_TE_iter=renormalize_thermal_expansion_iter, - # bulk_modulus=bulk_modulus, - # # prev_dir_hiphive=fit_force_constant.output["current_dir"], - # prev_dir_hiphive=prev_dir_hiphive, - # loop=loops, - # ) - # renormalization.name += f" {temperature} {loops}" - # renormalization.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) - # jobs.append(renormalization) - # outputs_renorm.append(renormalization.output) - # outputs.append(renormalization.output) - # renormalization.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"bulk_modulus={bulk_modulus}", - # f"run_renormalization_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) - - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-06-04-29-02-534935-41708" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-16-14-51-625183-73398" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-18-38-06-547834-98127" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-19-14-17-420858-98262" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-24-38-567687-98068" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-59-31-454856-98976" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-22-03-35-538120-20851" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-05-05-27-911205-37950" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-14-52-07-269546-30380" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-16-51-44-504507-24744" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-18-24-57-922525-28171" - # # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-15-01-03-47-788448-49252" - # # 5. Extract Phonon Band structure & DOS from FC - # # for 0K - # fc_pdos_pb_to_db = run_fc_to_pdos( - # renormalized=renormalize, - # mesh_density=mesh_density, - # prev_dir_json_saver=fit_force_constant.output["current_dir"], - # # prev_dir_json_saver=prev_dir_hiphive, - # loop=loops, - # ) - # fc_pdos_pb_to_db.name += f" {loops} 0K" - # jobs.append(fc_pdos_pb_to_db) - # outputs.append(fc_pdos_pb_to_db.output) - # fc_pdos_pb_to_db.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"cutoffs={cutoffs}", - # f"bulk_modulus={bulk_modulus}", - # f"temperature=0K" - # f"fc_pdos_pb_to_db_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) - # # for finite temperatures - # if renormalize: - # for i, temperature in enumerate(renormalize_temperature): - # fc_pdos_pb_to_db = run_fc_to_pdos( - # renormalized=renormalize, - # mesh_density=mesh_density, - # prev_dir_json_saver=outputs_renorm[i][0], - # loop=loops, - # ) - # fc_pdos_pb_to_db.name += f" {loops} {temperature}K" - # jobs.append(fc_pdos_pb_to_db) - # outputs.append(fc_pdos_pb_to_db.output) - # fc_pdos_pb_to_db.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"bulk_modulus={bulk_modulus}", - # f"temperature={temperature}K" - # f"fc_pdos_pb_to_db_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) - - - # # 6. Lattice thermal conductivity calculation using Sheng BTE - # if calculate_lattice_thermal_conductivity: - # if renormalize: - # temperatures = renormalize_temperature - # else: - # temperatures = thermal_conductivity_temperature - # # Because of the way ShengBTE works, a temperature array that is not - # # evenly spaced out (T_step) requires submission for each temperature - # if not renormalize: - # if isinstance(temperatures, dict): - # pass - # elif type(temperatures) in [list, np.ndarray]: - # assert all(np.diff(temperatures) == np.diff(temperatures)[0]) - # lattice_thermal_conductivity = run_lattice_thermal_conductivity( - # renormalized=renormalize, - # temperature=temperatures, - # loop=loops, - # prev_dir_hiphive=fit_force_constant.output["current_dir"], - # therm_cond_solver= self.THERM_COND_SOLVER - # ) - # lattice_thermal_conductivity.name += f" {temperatures} {loops}" - # lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) - # jobs.append(lattice_thermal_conductivity) - # outputs.append(lattice_thermal_conductivity.output) - # lattice_thermal_conductivity.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"bulk_modulus={bulk_modulus}", - # f"run_lattice_thermal_conductivity_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) - # else: - # for t, T in enumerate(temperatures): - # if T == 0: - # continue - # lattice_thermal_conductivity = run_lattice_thermal_conductivity( - # renormalized=renormalize, - # temperature=T, - # loop=loops, - # prev_dir_hiphive=outputs_renorm[t][0], - # therm_cond_solver= self.THERM_COND_SOLVER - # ) - - # lattice_thermal_conductivity.name += f" {T} {loops}" - # lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) - # jobs.append(lattice_thermal_conductivity) - # outputs.append(lattice_thermal_conductivity.output) - # lattice_thermal_conductivity.metadata.update( - # { - # "tag": [ - # f"mp_id={mpid}", - # f"bulk_modulus={bulk_modulus}", - # f"run_lattice_thermal_conductivity_{loops}", - # f"nConfigsPerStd={n_structures}", - # f"fixedDispls={fixed_displs}", - # f"dispCut={disp_cut}", - # f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", - # f"supercell_matrix={supercell_matrix}", - # f"loop={loops}", - # ] - # } - # ) + # 4. Perform phonon renormalization to obtain temperature-dependent + # force constants using hiPhive + outputs_renorm = [] + if renormalize: + for temperature in renormalize_temperature: + nconfig = renormalize_nconfig * (1 + temperature // 100) + renormalization = run_hiphive_renormalization( + temperature=temperature, + renorm_method=renormalize_method, + nconfig=nconfig, + renorm_TE_iter=renormalize_thermal_expansion_iter, + bulk_modulus=bulk_modulus, + prev_dir_hiphive=fit_force_constant.output["current_dir"], + # prev_dir_hiphive=prev_dir_hiphive, + loop=loops, + ) + renormalization.name += f" {temperature} {loops}" + renormalization.update_config({"manager_config": {"_fworker": "cpu_reg_fworker"}}) + jobs.append(renormalization) + outputs_renorm.append(renormalization.output) + outputs.append(renormalization.output) + renormalization.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"run_renormalization_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-06-04-29-02-534935-41708" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-16-14-51-625183-73398" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-18-38-06-547834-98127" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-19-14-17-420858-98262" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-24-38-567687-98068" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-20-59-31-454856-98976" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-07-22-03-35-538120-20851" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-05-05-27-911205-37950" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-14-52-07-269546-30380" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-16-51-44-504507-24744" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-08-18-24-57-922525-28171" + # prev_dir_hiphive = "/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/HPS_hiphive/job_2024-04-15-01-03-47-788448-49252" + # 5. Extract Phonon Band structure & DOS from FC + # for 0K + fc_pdos_pb_to_db = run_fc_to_pdos( + renormalized=renormalize, + mesh_density=mesh_density, + prev_dir_json_saver=fit_force_constant.output["current_dir"], + # prev_dir_json_saver=prev_dir_hiphive, + loop=loops, + ) + fc_pdos_pb_to_db.name += f" {loops} 0K" + jobs.append(fc_pdos_pb_to_db) + outputs.append(fc_pdos_pb_to_db.output) + fc_pdos_pb_to_db.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"cutoffs={cutoffs}", + f"bulk_modulus={bulk_modulus}", + f"temperature=0K" + f"fc_pdos_pb_to_db_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + # for finite temperatures + if renormalize: + for i, temperature in enumerate(renormalize_temperature): + fc_pdos_pb_to_db = run_fc_to_pdos( + renormalized=renormalize, + mesh_density=mesh_density, + prev_dir_json_saver=outputs_renorm[i][0], + loop=loops, + ) + fc_pdos_pb_to_db.name += f" {loops} {temperature}K" + jobs.append(fc_pdos_pb_to_db) + outputs.append(fc_pdos_pb_to_db.output) + fc_pdos_pb_to_db.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"temperature={temperature}K" + f"fc_pdos_pb_to_db_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + + + # 6. Lattice thermal conductivity calculation using Sheng BTE + if calculate_lattice_thermal_conductivity: + if renormalize: + temperatures = renormalize_temperature + else: + temperatures = thermal_conductivity_temperature + # Because of the way ShengBTE works, a temperature array that is not + # evenly spaced out (T_step) requires submission for each temperature + if not renormalize: + if isinstance(temperatures, dict): + pass + elif type(temperatures) in [list, np.ndarray]: + assert all(np.diff(temperatures) == np.diff(temperatures)[0]) + lattice_thermal_conductivity = run_lattice_thermal_conductivity( + renormalized=renormalize, + temperature=temperatures, + loop=loops, + prev_dir_hiphive=fit_force_constant.output["current_dir"], + therm_cond_solver= self.THERM_COND_SOLVER + ) + lattice_thermal_conductivity.name += f" {temperatures} {loops}" + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + jobs.append(lattice_thermal_conductivity) + outputs.append(lattice_thermal_conductivity.output) + lattice_thermal_conductivity.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"run_lattice_thermal_conductivity_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) + else: + for t, T in enumerate(temperatures): + if T == 0: + continue + lattice_thermal_conductivity = run_lattice_thermal_conductivity( + renormalized=renormalize, + temperature=T, + loop=loops, + prev_dir_hiphive=outputs_renorm[t][0], + therm_cond_solver= self.THERM_COND_SOLVER + ) + + lattice_thermal_conductivity.name += f" {T} {loops}" + lattice_thermal_conductivity.update_config({"manager_config": {"_fworker": "gpu_fworker"}}) + jobs.append(lattice_thermal_conductivity) + outputs.append(lattice_thermal_conductivity.output) + lattice_thermal_conductivity.metadata.update( + { + "tag": [ + f"mp_id={mpid}", + f"bulk_modulus={bulk_modulus}", + f"run_lattice_thermal_conductivity_{loops}", + f"nConfigsPerStd={n_structures}", + f"fixedDispls={fixed_displs}", + f"dispCut={disp_cut}", + f"supercell_matrix_kwargs={self.supercell_matrix_kwargs}", + f"supercell_matrix={supercell_matrix}", + f"loop={loops}", + ] + } + ) return Flow(jobs=jobs, output=outputs, name=f"{mpid}_{self.THERM_COND_SOLVER}_" f"{disp_cut}_" diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index 0f6aeb807f..5f6ff6b9de 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -36,6 +36,10 @@ from hiphive.structure_generation.random_displacement import ( generate_displaced_structures, ) +from hiphive.structure_generation.rattle import ( + generate_mc_rattled_structures, + generate_rattled_structures, +) from hiphive.utilities import get_displacements # Jobflow packages @@ -65,7 +69,10 @@ from atomate2.settings import Atomate2Settings from atomate2.utils.log import initialize_logger from atomate2.vasp.files import copy_hiphive_outputs - +from atomate2.forcefields.jobs import ( + CHGNetStaticMaker, + ForceFieldStaticMaker, +) if TYPE_CHECKING: from ase.atoms import Atoms from pymatgen.core.structure import Structure @@ -83,7 +90,7 @@ # T_KLAT = [300] # [i*100 for i in range(0,11)] T_KLAT = {"min":100,"max":1000,"step":100} #[i*100 for i in range(0,11)] T_THERMAL_CONDUCTIVITY = [0, 100, 200, 300] # [i*100 for i in range(0,16)] -IMAGINARY_TOL = 0.025 # in THz +IMAGINARY_TOL = 0.1 # in THz # changed from 0.025 FIT_METHOD = "rfe" eV2J = sp.constants.elementary_charge @@ -112,7 +119,7 @@ def hiphive_static_calcs( min_length: float | None = None, prefer_90_degrees: bool = True, n_structures: int | None = None, - # fixed_displs: list[float] = [0.01, 0.03, 0.08, 0.1], + fixed_displs: list[float] = [0.01, 0.03, 0.08, 0.1], loops: int | None = None, prev_dir: str | None = None, phonon_displacement_maker: BaseVaspMaker | None = None, @@ -145,7 +152,8 @@ def hiphive_static_calcs( structure=structure, supercell_matrix=supercell_matrix, n_structures=n_structures, - # fixed_displs=fixed_displs, + fixed_displs=fixed_displs, + prev_dir=prev_dir, loop=loops, ) jobs.append(displacements) @@ -176,129 +184,56 @@ def hiphive_static_calcs( return Response(replace=jobs, output=outputs) - -def find_smallest_x(n_structures: int) -> int: - # Starting from x = 3 because x - 2 in the denominator implies x > 2 - x = 3 - while True: - # Calculate the value of the expression for the current x - expression_value = (n_structures - 4) / (x - 2) - - # Check if the expression value satisfies the inequality - if 5 <= expression_value <= 6: - return math.floor(x) # Explicitly return the floor of x - x += 0.1 - @job def generate_hiphive_displacements( structure: Structure, + displacements: list[Structure] | None = None, supercell_matrix: list[list[int]] | None = None, n_structures: int | None = None, fixed_displs: list[float] | None = None, + prev_dir: str | None = None, loop: int | None = None, + high_force_atoms: list | None = None, ) -> list[Structure]: """Job generates the perturbed structures for hiPhive fitting.""" if fixed_displs is None: - # # fixed_displs = [0.01, 0.03, 0.08, 0.1] #TODO, this list is used in - # # the paper - # smallest_disp = 0.01 - - # # dictionary of largest displacements for each period - # largest_displacement = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} - - # row = int( - # np.around(np.array([s.row for s in structure.species]).mean(), 0) - # ) - - # largest_disp = largest_displacement[row] - - # ratio = int(largest_disp/smallest_disp) - # if 60 >= n_structures >= 4: - # logger.info(f"n_structures inside if statement >= 4 is {n_structures}") - # factors = np.sqrt(np.linspace(1,ratio**2,n_structures)) - # elif n_structures < 4 : - # factors = np.sqrt(np.linspace(1,ratio**2,4)) - # else: - # factors = np.sqrt(np.linspace(1,ratio**2,60)) - # amplitudes = (smallest_disp*factors).round(3) - # logger.info(f"amplitudes = {amplitudes}") - # fixed_displs = amplitudes.tolist() - # logger.info(f"list_amplitudes = {fixed_displs}") - - if n_structures < 4: - n_structures = 4 - else: - pass - logger.info(f"n_structures = {n_structures}") - # # Ask the user for n_structures - # n_structures = int(input("Enter the value of n_structures: ")) - - # # Find the smallest integer x that satisfies the inequality - # smallest_x = find_smallest_x(n_structures) - - smallest_x = 2 - - # Output the result - logger.info(f"The smallest integer x that satisfies the inequality is: {smallest_x}") - - # Create structures at equal interval between disp = 0.01 & 0.055 - # equal_interval_structures = np.linspace(0.01, 0.055, smallest_x) - equal_interval_structures = np.linspace(0.001, 0.005, smallest_x) - # create a new np array with 3 elements equaling 0.03 and 0.08 respectively - equal_interval_structures = np.array([0.03, 0.08]) - - # # Remaining number of structures - # remaining_structures = n_structures - smallest_x - - - # # structure_data = loadfn(f"/Users/HPSahasrabuddhe/Desktop/Acads/3rd_sem/MSE 299/Hiphive_Atomate2_integration/hiphive_4622_VASP_96/launcher_2024-04-01-08-09-15-480067/structure_data_1.json") - # # structure = structure_data["structure"] - - # largest_period_default = {1: 0.1, 2: 0.2, 3: 0.2, 4: 0.3, 5: 0.4, 6: 0.5, 7: 0.5} - - # row = int( - # np.around(np.array([s.row for s in structure.species]).mean(), 0) - # ) - # logger.info(f"row = {row}") - - # largest = largest_period_default[row] - # logger.info(f"largest = {largest}") - # smallest = 0.056 - - # ratio = largest/smallest - # logger.info(f"ratio = {ratio}") - # logger.info(f"smallest = {smallest}") - # logger.info(f"largest = {largest}") - # logger.info(f"remaining_structures = {remaining_structures}") - # factors = np.sqrt(np.linspace(1,ratio**2,remaining_structures)) - # logger.info(f"factors = {factors}") - # amplitudes = (smallest*factors).round(3) - # logger.info(f"amplitudes = {amplitudes}") - - # # Combine the two sets of structures - # all_amplitudes = np.concatenate((equal_interval_structures, amplitudes)) - - # fixed_displs = all_amplitudes.tolist() - # logger.info(f"all_amplitudes = {all_amplitudes}") - - fixed_displs = equal_interval_structures.tolist() - logger.info(f"all_amplitudes = {fixed_displs}") - + fixed_displs = [0.01, 0.03, 0.08, 0.1] #TODO, this list is used in the paper + fixed_displs = [0.01, 0.03, 0.4, 0.6] #TODO, this list is used in the paper + fixed_displs = [0.2, 0.4, 0.6, 0.8] #TODO, this list is used in the paper + logger.info(f"default fixed_displs = {fixed_displs}") + logger.info(f"default fixed_displs = {fixed_displs}") logger.info(f"supercell_matrix = {supercell_matrix}") + # supercell_matrix = [[3, 0, 0], [0, 3, 0], [0, 0, 3]] supercell_structure = SupercellTransformation( scaling_matrix=supercell_matrix ).apply_transformation(structure) logger.info(f"supercell_structure = {supercell_structure}") + structure_data = { + "structure": structure, + "supercell_structure": supercell_structure, + "supercell_matrix": supercell_matrix, + } + + dumpfn(structure_data, f"structure_data_{loop}.json") # Generate the rattled structures #### structures_ase_all = [] logger.info(f"n_structures = {n_structures}") + structures_ase_all_prev = [] # Convert to ASE atoms + if high_force_atoms is not None: + logger.info(f"len_high_force_atoms = {len(high_force_atoms)}") + logger.info(f"high_force_atoms = {high_force_atoms}") + for pymatgen_structure in displacements: + atom = AseAtomsAdaptor.get_atoms(pymatgen_structure) + structures_ase_all_prev.append(atom) for i in range(len(fixed_displs)): supercell_ase = AseAtomsAdaptor.get_atoms(supercell_structure) structures_ase = generate_displaced_structures( - atoms=supercell_ase, n_structures=n_structures, distance=fixed_displs[i], loop=loop + atoms=supercell_ase, n_structures=n_structures, distance=fixed_displs[i], + loop=loop, high_force_atoms=high_force_atoms, displ_number=i, + structures_ase_all_prev=structures_ase_all_prev ) structures_ase_all.extend(structures_ase) @@ -316,6 +251,7 @@ def generate_hiphive_displacements( for i in range(len(structures_pymatgen)): structures_pymatgen[i].to(f"POSCAR_{i}", "poscar") + dumpfn(structures_pymatgen, f"perturbed_structures_{loop}.json") return structures_pymatgen @job @@ -342,6 +278,7 @@ def collect_perturbed_structures( """ logger.info(f"scaling_matrix = {supercell_matrix}") logger.info(f"structure = {structure}") + # supercell_matrix = [[3, 0, 0], [0, 3, 0], [0, 0, 3]] # TODO: remove this line later supercell = SupercellTransformation(scaling_matrix=supercell_matrix).apply_transformation(structure) structure_data = { "structure": structure, @@ -405,114 +342,6 @@ def collect_perturbed_structures( return [all_structures, all_forces, structure_data, current_dir] - -# @job -# def quality_control( -# rmse_test: float, -# n_structures: int, -# fixed_displs: list[float], -# loop: int, -# fit_method: str, -# disp_cut: float, -# bulk_modulus: float, -# temperature_qha: float, -# mesh_density: float, -# imaginary_tol: float, -# prev_dir_json_saver: str, -# prev_dir: str, -# supercell_matrix_kwargs: list[list[int]], -# ): -# """ -# Check if the desired Test RMSE is achieved. - -# If not, then increase the number of structures -# """ -# if rmse_test > 0.010: -# return Response( -# addition=quality_control_job( -# rmse_test, -# n_structures, -# fixed_displs, -# loop, -# fit_method, -# disp_cut, -# bulk_modulus, -# temperature_qha, -# mesh_density, -# imaginary_tol, -# prev_dir_json_saver, -# prev_vasp_dir, -# supercell_matrix_kwargs, -# ) -# ) -# return None - - -# @job -# def quality_control_job( -# rmse_test, -# n_structures: int, -# fixed_displs: List[float], -# loop: int, -# fit_method: str, -# disp_cut: float, -# bulk_modulus: float, -# temperature_qha: float, -# mesh_density: float, -# imaginary_tol: float, -# prev_dir_json_saver: str, -# prev_vasp_dir: str, -# supercell_matrix_kwargs: List[List[int]], -# ): -# """Increases the number of structures if the desired Test RMSE is not achieved.""" -# jobs = [] -# outputs = [] - -# # 4. Quality Control Job to check if the desired Test RMSE is achieved, -# # if not, then increase the number of structures -- -# # Using "addition" feature of jobflow -# loop += 1 -# n_structures += 1 -# # prev_dir_json_saver="/pscratch/sd/h/hrushi99/atomate2/InAs/ -# # block_2023-06-16-04-09-51-792824/launcher_2023-06-23-23-58-57-102993/ -# # launcher_2023-06-23-23-59-34-157381" -# error_check_job = quality_control( -# rmse_test=fw_fit_force_constant.output[5], -# n_structures=n_structures, -# fixedDispls=fixed_displs, -# loop=loop, -# fit_method=fit_method, -# disp_cut=disp_cut, -# bulk_modulus=bulk_modulus, -# temperature_qha=temperature_qha, -# mesh_density=mesh_density, -# imaginary_tol=imaginary_tol, -# prev_dir_json_saver=json_saver.output[3], -# # prev_dir_json_saver = prev_dir_json_saver, -# prev_vasp_dir=prev_vasp_dir, -# supercell_matrix_kwargs=supercell_matrix_kwargs, -# ) -# error_check_job.name += f" {loop}" -# jobs.append(error_check_job) -# outputs.append(error_check_job.output) -# error_check_job.metadata.update( -# { -# "tag": [ -# f"error_check_job_{loop}", -# f"nConfigsPerStd={n_structures}", -# f"fixedDispls={fixed_displs}", -# f"dispCut={disp_cut}", -# f"supercell_matrix_kwargs={supercell_matrix_kwargs}", -# f"loop={loop}", -# ] -# } -# ) - -# flow = Flow(jobs=jobs, output=outputs) - -# quality_control_job.name = f"quality_control_job {loop}" - -# return Response(addition=flow) @job def run_hiphive_individually( mpid: str = None, @@ -712,50 +541,22 @@ def run_hiphive( copy_hiphive_outputs(prev_dir_json_saver) all_structures = loadfn(f"perturbed_structures_{loop}.json") + # all_structures = loadfn(f"perturbed_structures.json") all_forces = loadfn(f"perturbed_forces_{loop}_new.json") + # all_forces = loadfn(f"perturbed_forces.json") structure_data = loadfn(f"structure_data_{loop}.json") + # structure_data = loadfn(f"structure_data.json") parent_structure = structure_data["structure"] supercell_structure = structure_data["supercell_structure"] supercell_matrix = np.array(structure_data["supercell_matrix"]) - parent_structure = SpacegroupAnalyzer(parent_structure).find_primitive() #TODO refactor this later + # parent_structure = SpacegroupAnalyzer(parent_structure).find_primitive() #TODO refactor this later if cutoffs is None: cutoffs = get_cutoffs(supercell_structure) logger.info(f"cutoffs is {cutoffs}") - # cutoffs = [[6, 3.5, 3]] - # cutoffs = [[9.0, 5.5, 4.25]] - # cutoffs = [[8.0]] - # cutoffs = [[9.0, 6.5, 5]] - # cutoffs = [[9.0, 5.5, 3]] - # cutoffs = [[9, 9.25, 5]] - # cutoffs = [[9, 8.3125, 5]] - # cutoffs = [[4, 4, 4]] - # cutoffs = [[3, 2.5]] - # cutoffs = [[2, 2.7]] - # cutoffs = [[6, 4, 2]] - # cutoffs = [[1.8635, 4]] - # cutoffs = [[2.559, 4]] # results in 2nd, 3rd order -- 124, 7379 DOFs -- imaginary modes = 8 - # cutoffs = [[2.54, 4]] # results in 2nd, 3rd order -- 106, 7379 DOFs -- imaginary modes = 2 - # cutoffs = [[2.5, 4]] # results in 2nd, 3rd order -- 88, 7379 DOFs -- imaginary modes = 0 - # cutoffs = [[2.52, 4]] # results in 2nd, 3rd order -- 97, 7379 DOFs -- imaginary modes = 2 - # cutoffs = [[2.52, 3]] # results in 2nd, 3rd order -- 97, 857 DOFs -- imaginary modes = 0 - # cutoffs = [[2.52, 4.6]] # results in 2nd, 3rd order -- 97, 20804 DOFs -- imaginary modes = 2 - # cutoffs = [[2.5, 4.6]] # results in 2nd, 3rd order -- 88, 20804 DOFs -- imaginary modes = 0 - # cutoffs = [[2.5, 5.6]] # results in 2nd, 3rd order -- 88, 67511 DOFs -- imaginary modes = ? - cutoffs = [[9.0, 5.5, 3.5]] - cutoffs = [[9, 8.3125, 5]] - cutoffs = [[6, 5, 4]] # results in {2: 36, 3: 263, 4: 327} DOFs -- imaginary modes = 18, rmse = 132meV/atom - cutoffs = [[8, 5, 4]] # results in {2: 79, 3: 263, 4: 327} DOFs -- imaginary modes = 174, rmse = 68.8meV/atom - cutoffs = [[10, 5, 4]] # results in {2: 151, 3: 263, 4: 327} DOFs -- imaginary modes = 150, rmse = 40.5meV/atom - cutoffs = [[10, 7, 4]] # results in {2: 151, 3: 2149, 4: 327} DOFs -- imaginary modes = 150, rmse = 40.1meV/atom - cutoffs = [[10, 7, 5]] # results in {2: 151, 3: 2149, 4: 2876} DOFs -- imaginary modes = 150, rmse = 38.03meV/atom - cutoffs = [[9, 6, 4]] # {2: 102, 3: 775, 4: 428} -- imaginary modes = 114, rmse = 28.9meV/atom -- MgO 1265 - cutoffs = [[9, 7, 4.3]] # {2: 102, 3: 2593, 4: 806} -- imaginary modes = 114, rmse = 29.2meV/atom -- MgO 1265 - cutoffs = [[9, 6, 4.5]] # {2: 74, 3: 597, 4: 428} -- imaginary modes = 0, rmse = 13.3meV/atom -- CaO 2605 - cutoffs = [[10, 6, 4.5]] # {2: 88, 3: 597, 4: 428} -- imaginary modes = 0, rmse = 9.6meV/atom -- CaO 2605 - cutoffs = [[10, 7, 5]] # {2: 88, 3: 775, 4: 806} -- imaginary modes = 0, rmse = 9.6meV/atom -- CaO 2605 + cutoffs = [[9, 4, 3]] # mp-23339 {2: 1513, 3: 787, 4: 210} logger.info(f"cutoffs is {cutoffs}") else: pass @@ -785,16 +586,22 @@ def run_hiphive( atoms.positions = supercell_atoms.get_positions() structures.append(atoms) + # Calculate mean displacements mean_displacements = np.linalg.norm(displacements, axis=1).mean() logger.info(f"Mean displacements while reading individual displacements: " f"{mean_displacements}") + # Calculate standard deviation of displacements + std_displacements = np.linalg.norm(displacements, axis=1).std() + logger.info(f"Standard deviation of displacements while reading individual displacements: " + f"{std_displacements}") all_cutoffs = cutoffs logger.info(f"all_cutoffs is {all_cutoffs}") - fcs, param, cs, fitting_data = fit_force_constants( + fcs, param, cs, fitting_data, fcp = fit_force_constants( parent_structure=parent_structure, supercell_matrix=supercell_matrix, + supercell_structure=supercell_structure, structures=structures, all_cutoffs=all_cutoffs, disp_cut=disp_cut, @@ -802,9 +609,6 @@ def run_hiphive( fit_method=fit_method, ) - if fcs is None: - raise RuntimeError("Could not find a force constant solution") - logger.info("Saving Harmonic props") thermal_data, phonopy = harmonic_properties( parent_structure, supercell_matrix, fcs, t_qha, imaginary_tol @@ -818,17 +622,9 @@ def run_hiphive( thermal_data["n_imaginary"], bulk_modulus, ) - logger.info("Saving phonopy_params") - phonopy.save("phonopy_params.yaml") - fitting_data["n_imaginary"] = thermal_data.pop("n_imaginary") - thermal_data.update(anharmonic_data) - logger.info("Saving fitting_data") - dumpfn(fitting_data, "fitting_data.json") - logger.info("Saving thermal_data") - dumpfn(thermal_data, "thermal_data.json") - logger.info("Writing cluster space and force_constants") - logger.info(f"{type(fcs)}") + if fcs is None: + raise RuntimeError("Could not find a force constant solution") if isinstance(fcs, ForceConstants): logger.info("Writing force_constants") @@ -836,6 +632,10 @@ def run_hiphive( else: logger.info("fcs is not an instance of ForceConstants") + if isinstance(fcp, ForceConstantPotential): + logger.info("Writing force_constants_potential") + fcp.write("force_constants_potential.fcp") + logger.info("Saving parameters") np.savetxt("parameters.txt", param) @@ -846,6 +646,18 @@ def run_hiphive( else: logger.info("cs is not an instance of ClusterSpace") + logger.info("Saving phonopy_params") + phonopy.save("phonopy_params.yaml") + fitting_data["n_imaginary"] = thermal_data.pop("n_imaginary") + thermal_data.update(anharmonic_data) + logger.info("Saving fitting_data") + dumpfn(fitting_data, "fitting_data.json") + logger.info("Saving thermal_data") + dumpfn(thermal_data, "thermal_data.json") + + logger.info("Writing cluster space and force_constants") + logger.info(f"{type(fcs)}") + # # following code is commented only for testing puropose # if fitting_data["n_imaginary"] == 0: # # following code is commented only for testing puropose ends @@ -956,16 +768,16 @@ def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: min_cutoffs = { - # 2: {1: 7.0, 2: 8.0, 3: 9.0, 4: 10.0, 5: 11.0, 6: 12.0, 7: 13.0}, - 2: {1: 5.0, 2: 6.0, 3: 7.0, 4: 8.0, 5: 9.0, 6: 10.0, 7: 11.0}, + 2: {1: 7.0, 2: 8.0, 3: 9.0, 4: 10.0, 5: 11.0, 6: 12.0, 7: 13.0}, + # 2: {1: 5.0, 2: 6.0, 3: 7.0, 4: 8.0, 5: 9.0, 6: 10.0, 7: 11.0}, # 2: {1: 6.0, 2: 7.0, 3: 8.0, 4: 9.0, 5: 10.0, 6: 11.0, 7: 12.0}, 3: {1: 2.5, 2: 3.0, 3: 4.0, 4: 5.0, 5: 5.5, 6: 6.0, 7: 6.5}, # 3: {1: 3.0, 2: 3.5, 3: 4.5, 4: 5.5, 5: 6.0, 6: 6.5, 7: 7.0}, 4: {1: 1.5, 2: 2.0, 3: 2.5, 4: 3.0, 5: 3.5, 6: 4.0, 7: 4.5}, # 4: {1: 2.5, 2: 3.0, 3: 3.5, 4: 4.0, 5: 4.5, 6: 5.0, 7: 5.5}, } - inc = {2: 0, 3: 2.0, 4: 1.2} - steps = {2: 0, 3: 0.75, 4: 0.6} + inc = {2: 2, 3: 2.0, 4: 1.2} + steps = {2: 1, 3: 0.75, 4: 0.6} row = int( np.around(np.array([s.row for s in supercell_structure.species]).mean(), 0) @@ -977,7 +789,10 @@ def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: mins = {2: min_cutoffs[2][row], 3: min_cutoffs[3][row], 4: min_cutoffs[4][row]} # create an NDArray of 2nd order cutofss with only one entry -> mins[2], and 3rd and 4th order cutoffs with a range - range_two = np.array([mins[2]]) + # range_two = np.array([mins[2]]) + range_two = np.arange( + mins[2], mins[2] + factor * (inc[2] + steps[2]), factor * steps[2] + ) range_three = np.arange( mins[3], mins[3] + factor * (inc[3] + steps[3]), factor * steps[3] ) @@ -1004,6 +819,7 @@ def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: def fit_force_constants( parent_structure: Structure, supercell_matrix: np.ndarray, + supercell_structure: Structure, structures: list[Atoms], all_cutoffs: list[list[float]], # separate_fit: bool, @@ -1100,13 +916,14 @@ def fit_force_constants( cutoff_results = Parallel(n_jobs=min(os.cpu_count(),len(all_cutoffs)), backend="multiprocessing")(delayed(_run_cutoffs)( - i, cutoffs, n_cutoffs, parent_structure, structures, supercell_matrix, - fit_method, disp_cut, fit_kwargs) for i, cutoffs in enumerate(all_cutoffs)) + i, cutoffs, n_cutoffs, parent_structure, supercell_structure, structures, + supercell_matrix, fit_method, disp_cut, fit_kwargs) for i, cutoffs in enumerate(all_cutoffs)) for result in cutoff_results: if result is None: + print("result is None") continue - + print(f"result = {result}") fitting_data["cutoffs"].append(result["cutoffs"]) fitting_data["rmse_test"].append(result["rmse_test"]) # fitting_data["n_imaginary"].append(result["n_imaginary"]) @@ -1125,13 +942,14 @@ def fit_force_constants( logger.info("Finished fitting force constants.") - return best_fit["force_constants"], best_fit["parameters"], best_fit["cluster_space"], fitting_data + return best_fit["force_constants"], best_fit["parameters"], best_fit["cluster_space"], fitting_data, best_fit["force_constants_potential"] def _run_cutoffs( i: int, cutoffs: list[float], n_cutoffs: int, parent_structure: Structure, + supercell_structure: Structure, structures: list[Atoms], supercell_matrix: np.ndarray, # shape=(3, 3), dtype='intc', order='C'., fit_method: str, @@ -1139,13 +957,16 @@ def _run_cutoffs( fit_kwargs: dict[str, Any], ) -> dict[str, Any]: logger.info(f"Testing cutoffs {i+1} out of {n_cutoffs}: {cutoffs}") - supercell_atoms = structures[0] + supercell_atoms = structures[0] #TODO: only for testing purposes + # supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + logger.info(f"supercell_atoms = {supercell_atoms}") if not is_cutoff_allowed(supercell_atoms, max(cutoffs)): logger.info("Skipping cutoff due as it is not commensurate with supercell size") return {} cs = ClusterSpace(supercell_atoms, cutoffs, symprec=1e-3, acoustic_sum_rules=True) + # cs = ClusterSpace(supercell_atoms, cutoffs, symprec=1e-1, acoustic_sum_rules=True) logger.debug(cs.__repr__()) n2nd = cs.get_n_dofs_by_order(2) # change it back to cs.get_n_dofs_by_order(2) nall = cs.n_dofs @@ -1179,12 +1000,25 @@ def _run_cutoffs( phonopy.run_mesh(mesh, with_eigenvectors=False, is_mesh_symmetry=False) omega = phonopy.mesh.frequencies # THz omega = np.sort(omega.flatten()) - logger.info(f"omega = {omega}") + logger.info(f"omega_one_shot_fit = {omega}") imaginary = np.any(omega < -1e-3) - logger.info(f"imaginary = {imaginary}") + logger.info(f"imaginary_one_shot_fit = {imaginary}") # commenting it out only for testing purposes + # Phonopy's way of calculating phonon frequencies + structure_phonopy = get_phonopy_structure(parent_structure) + phonon = Phonopy(structure_phonopy, supercell_matrix=supercell_matrix) + phonon.set_force_constants(fcs.get_fc_array(2)) + phonon.run_mesh(mesh=100.0, is_mesh_symmetry=False, is_gamma_center=True) + mesh = phonon.get_mesh_dict() + omega = mesh["frequencies"] + omega = np.sort(omega.flatten()) + logger.info(f"omega_phonopy_one_shot_fitting = {omega}") + imaginary = np.any(omega < -1e-3) + logger.info(f"imaginary_phonopy_one_shot_fitting = {imaginary}") + if imaginary: + # if False: # if True: # only for testing purposes logger.info( "Imaginary modes found! Fitting anharmonic force constants separately" @@ -1202,7 +1036,40 @@ def _run_cutoffs( assert(nall==len(parameters)) logger.info(f"Training complete for cutoff: {i}, {cutoffs}") + fcp = ForceConstantPotential(cs, parameters) + logger.info(f"supercell atoms = {supercell_atoms}") + fcs = fcp.get_force_constants(supercell_atoms) + logger.info("Did you get the large Condition number error?") + + parent_phonopy = get_phonopy_structure(parent_structure) + phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) + phonopy.primitive.get_number_of_atoms() + mesh = supercell_matrix.diagonal() * 2 + phonopy.set_force_constants(fcs.get_fc_array(2)) + phonopy.set_mesh( + mesh, is_eigenvectors=False, is_mesh_symmetry=False + ) # run_mesh(is_gamma_center=True) + phonopy.run_mesh(mesh=100.0, with_eigenvectors=False, is_mesh_symmetry=False) + omega = phonopy.mesh.frequencies # THz + omega = np.sort(omega.flatten()) + logger.info(f"omega_seperate_fit = {omega}") + imaginary = np.any(omega < -1e-3) + logger.info(f"imaginary_seperate_fit = {imaginary}") + + # Phonopy's way of calculating phonon frequencies + structure_phonopy = get_phonopy_structure(parent_structure) + phonon = Phonopy(structure_phonopy, supercell_matrix=supercell_matrix) + phonon.set_force_constants(fcs.get_fc_array(2)) + phonon.run_mesh(mesh, is_mesh_symmetry=False, is_gamma_center=True) + mesh = phonon.get_mesh_dict() + omega = mesh["frequencies"] + omega = np.sort(omega.flatten()) + logger.info(f"omega_phonopy_seperate_fit = {omega}") + imaginary = np.any(omega < -1e-3) + logger.info(f"imaginary_phonopy_seperate_fit = {imaginary}") + else: + # if True: logger.info("No imaginary modes! Fitting all force constants in one shot") separate_fit = False sc = get_structure_container( @@ -1271,12 +1138,18 @@ def get_structure_container( A hiPhive StructureContainer. """ sc = StructureContainer(cs) + logger.info(f"sc = {sc}") logger.info(f"initial shape of fit matrix = {sc.data_shape}") saved_structures = [] for _, structure in enumerate(structures): displacements = structure.get_array("displacements") + # Calculate mean displacements mean_displacements = np.linalg.norm(displacements, axis=1).mean() logger.info(f"Mean displacements: {mean_displacements}") + # Calculate standard deviation of displacements + std_displacements = np.linalg.norm(displacements, axis=1).std() + logger.info(f"Standard deviation of displacements: " + f"{std_displacements}") if not separate_fit: # fit all sc.add_structure(structure) else: # fit separately @@ -1337,6 +1210,8 @@ def harmonic_properties( frequency at Gamma, and the free energy, entropy, and heat capacity """ + # fcs = ForceConstants.read("force_constants.fcs") + logger.info('Evaluating harmonic properties...') fcs2 = fcs.get_fc_array(2) fcs3 = fcs.get_fc_array(3) @@ -1353,18 +1228,41 @@ def harmonic_properties( logger.info('Thermal properties successfully run!') _, free_energy, entropy, heat_capacity = phonopy.get_thermal_properties() - free_energy *= 1000/sp.constants.Avogadro/eV2J/natom # kJ/mol to eV/atom - entropy *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom - heat_capacity *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom + # free_energy *= 1000/sp.constants.Avogadro/eV2J/natom # kJ/mol to eV/atom + # entropy *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom + # heat_capacity *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom logger.info(f"Heat_capacity_harmonic_property: {heat_capacity}") freq = phonopy.mesh.frequencies # in THz logger.info(f'Frequencies: {freq}') + logger.info(f"freq_flatten = {np.sort(freq.flatten())}") # find imaginary modes at gamma # phonopy.run_qpoints([0, 0, 0]) # gamma_eigs = phonopy.get_qpoints_dict()["frequencies"] n_imaginary = int(np.sum(freq < -np.abs(imaginary_tol))) min_freq = np.min(freq) + + # # Phonopy's way of calculating phonon frequencies + # # delete this block of code later. Only for testing purposes + # logger.info(f"Mesh for calculating phonon frequencies using Phonopy's way = {mesh}") + # mesh = 100 + # logger.info(f"Mesh for calculating phonon frequencies using Phonopy's way = {mesh}") + # # fcs = ForceConstants.read("force_constants.fcs") + # fcs2 = fcs.get_fc_array(2) + # structure_phonopy = get_phonopy_structure(structure) + # phonon = Phonopy(structure_phonopy, supercell_matrix=supercell_matrix) + # phonon.set_force_constants(fcs2) + # phonon.run_mesh(mesh, is_mesh_symmetry=False, is_gamma_center=True) + # mesh = phonon.get_mesh_dict() + + # omega = mesh["frequencies"] + # omega = np.sort(omega.flatten()) + # logger.info(f"omega_phonopy_while_harmonic_prop = {omega}") + # imaginary = np.any(omega < -1e-3) + # logger.info(f"imaginary_phonopy_while_harmonic_prop = {imaginary}") + # n_imaginary = int(np.sum(omega < -np.abs(imaginary_tol))) + # logger.info(f"n_imaginary_phonopy_while_harmonic_prop = {n_imaginary}") + if n_imaginary == 0: logger.info('No imaginary modes!') else: # do not calculate these if imaginary modes exist @@ -1383,7 +1281,7 @@ def harmonic_properties( "entropy": entropy, "heat_capacity": heat_capacity, "n_imaginary": n_imaginary - }, phonopy + }, phonopy #phonon def anharmonic_properties( @@ -1485,7 +1383,7 @@ def gruneisen( fcs2: np.ndarray, fcs3: np.ndarray, temperature: List, - heat_capacity: np.ndarray, # in eV/K/atom + heat_capacity: np.ndarray, # in J/K/mol bulk_modulus: float = None # in GPa ) -> Tuple[List,List]: @@ -1506,7 +1404,8 @@ def gruneisen( cte = None dLfrac = None else: - heat_capacity *= eV2J*phonopy.primitive.get_number_of_atoms() # eV/K/atom to J/K + # heat_capacity *= eV2J*phonopy.primitive.get_number_of_atoms() # eV/K/atom to J/K + heat_capacity *= 1/sp.constants.Avogadro # J/K/mol to J/K # to convert from J/K/atom multiply by phonopy.primitive.get_number_of_atoms() # Convert heat_capacity to an array if it's a scalar # heat_capacity = np.array([heat_capacity]) logger.info(f"heat capacity = {heat_capacity}") @@ -1702,6 +1601,7 @@ def run_fc_to_pdos( mesh_density = mesh_density if mesh_density else 100.0 structure_data = loadfn(f"structure_data_{loop}.json") + # structure_data = loadfn(f"structure_data.json") structure = structure_data["structure"] structure_data["supercell_structure"] supercell_matrix = structure_data["supercell_matrix"] @@ -1739,6 +1639,12 @@ def run_fc_to_pdos( # following line is commented out only for testing purposes ends # pass + # convert uniform_bs to dict + uniform_bs_dict = uniform_bs.as_dict() + + # dump uniform_bs_dict to file + dumpfn(uniform_bs_dict, "uniform_bs.json") + return uniform_bs, lm_bs, dos, prev_dir_json_saver @@ -1912,16 +1818,16 @@ def run_hiphive_renormalization( # FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD_{temperature}", atoms, order=3, fc_tol=1e-4) - atoms = AseAtomsAdaptor.get_atoms(parent_structure) - supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) - # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) - fcs_TD.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") - ForceConstants.write_to_phonopy(fcs_TD, "fc2.hdf5", "hdf5") - ForceConstants.write_to_phono3py(fcs_TD, "fc3.hdf5", "hdf5") - ### detour from hdf5 - # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) - FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) - FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) + # atoms = AseAtomsAdaptor.get_atoms(parent_structure) + # supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + # # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3) + # fcs_TD.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") + # ForceConstants.write_to_phonopy(fcs_TD, "fc2.hdf5", "hdf5") + # ForceConstants.write_to_phono3py(fcs_TD, "fc3.hdf5", "hdf5") + # ### detour from hdf5 + # # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + # FCS = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + # FCS.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) dumpfn(TD_structure_data, "structure_data.json") dumpfn(TD_thermal_data, "thermal_data.json")