``` ##────────────────────────────────────────────────────────────────────── ## ## ## ## [ VASP - SCF File Generation] ## ## 1. Call structure from Materials Project API ## ## 2. Convert to Primitive structure ## ## 3. Create POSCAR, INCAR for each structure ## ## 4. Change some INCAR options to user settings ## ## 5. Keep POTCAR writing order same as POSCAR ## ## ## ##=======================================================================## ## ## ## [ Changelog] ## ## 2025.04.23: Modified some INCAR options to user settings ## ## 2025.04.23: Improved KPOINTS generation density from 100 -> 400 ## ## 2025.05.05: Modified POTCAR reading order to normal, _pv, _sv ## ## 2025.05.08: Changed POTCAR to VASP recommended (removed) ## ## 2025.05.08: Changed POTCAR to MPStaticSet recommended ## ## 2025.07.01: Modified MAGMOM change mechanism ## ## 2025.07.01: Overall code structure and detail modifications ## ## 2025.07.05: Re-modified MAGMOM setting method ## ## 1) If CONFIG not used, use site property ## ## 2) Site property loaded from primitive and attached ## ## 3) Confirmed MAGMOM matches site property, ## ## but note there may be errors ## ## 2025.08.29: Sorted atom order in POTCAR writing to match POSCAR order## ## 2025.09.01: Added printing of total material count and progress ## ## 2025.09.15: Added writing Magnetization data to structure.txt file ## ## 2025.09.24: Added exception handling for cases when KPOINTS not generated ## ## 2025.09.25: Added writing calculation origin (GGA, GGA+U, r2SCAN) of Magnetism data ## ## 2025.09.25: Added writing calculation origin (GGA, GGA+U, r2SCAN) read from ## ## Origin for Structure, Energy, and Magnetism data ## ## 2025.09.29: Reverted from r2SCAN to GGA calculation, commented out r2SCAN options ## ## ## ##=======================================================================## ## ## ## Written by G. Maeng. ## ## Last Update : 2025.09.25 ## ## ## ##───────────────────────────────────────────────────────────────────────## import os # Basic module for file path and directory management import subprocess # Module for executing external commands (used here for merging POTCAR via 'cat' command) import traceback from mp_api.client import MPRester # Client for accessing Materials Project API from pymatgen.io.vasp.sets import MPRelaxSet, MPStaticSet, MPScanStaticSet # MP's standard VASP input set classes from pymatgen.analysis.bond_valence import BVAnalyzer # Oxidation state analyzer (import added) # ------------------------------------------------------------------ # 0. Basic settings # ------------------------------------------------------------------ API_KEY = "-" # Next-Gen MP API, inu e-mail key (not the MDT Lab e-mail key) POTCAR_DIR = "-" # final_structure = True && config_magmom = True : Set to value defined by MAGMOM + set to 0 if MAGMOM not specified # final_structure = True && config_magmom = False : Preferentially use Site property # final_structure = False && config_magmom = True : Output only values set by MAGMOM to INCAR # final Structure = False && config_magmom = False : Unify MAGMOM to 0.6 use_final_structure = True # True: final structure, False: initial structure USE_CONFIG_MAGMOM = False # True: force CONFIG_MAGMOM, False: automatic # ------------------------------------------------------------------ # 1. CONFIG_MAGMOM table # ------------------------------------------------------------------ CONFIG_MAGMOM = { 'Ce':5,'Ce3+':1,'Co':0.6,'Co3+':0.6,'Co4+':1,'Cr':5,'Dy3+':5,'Er3+':3,'Eu':10,'Eu2+':7, 'Eu3+':6,'Fe':5,'Gd3+':7,'Ho3+':4,'La3+':0.6,'Lu3+':0.6,'Mn':5,'Mn3+':4,'Mn4+':3,'Mo':5, 'Nd3+':3,'Ni':5,'Pm3+':4,'Pr3+':2,'Sm3+':5,'Tb3+':6,'Tm3+':2,'V':5,'W':5,'Yb3+':1 } # ------------------------------------------------------------------ # 2. Recommended POTCAR list (MPStaticSet recommended POTCAR) # ------------------------------------------------------------------ RECOMMENDED_POTCARS = { 'Ac':'Ac','Ag':'Ag','Al':'Al','Ar':'Ar','As':'As','Au':'Au','B':'B','Ba':'Ba_sv','Be':'Be_sv', 'Bi':'Bi','Br':'Br','C':'C','Ca':'Ca_sv','Cd':'Cd','Ce':'Ce','Cl':'Cl','Co':'Co','Cr':'Cr_pv', 'Cs':'Cs_sv','Cu':'Cu_pv','Dy':'Dy_3','Er':'Er_3','Eu':'Eu','F':'F','Fe':'Fe_pv','Ga':'Ga_d', 'Gd':'Gd','Ge':'Ge_d','H':'H','He':'He','Hf':'Hf_pv','Hg':'Hg','Ho':'Ho_3','I':'I','In':'In_d', 'Ir':'Ir','K':'K_sv','Kr':'Kr','La':'La','Li':'Li_sv','Lu':'Lu_3','Mg':'Mg_pv','Mn':'Mn_pv', 'Mo':'Mo_pv','N':'N','Na':'Na_pv','Nb':'Nb_pv','Nd':'Nd_3','Ne':'Ne','Ni':'Ni_pv','Np':'Np', 'O':'O','Os':'Os_pv','P':'P','Pa':'Pa','Pb':'Pb_d','Pd':'Pd','Pm':'Pm_3','Pr':'Pr_3','Pt':'Pt', 'Pu':'Pu','Rb':'Rb_sv','Re':'Re_pv','Rh':'Rh_pv','Ru':'Ru_pv','S':'S','Sb':'Sb','Sc':'Sc_sv', 'Se':'Se','Si':'Si','Sm':'Sm_3','Sn':'Sn_d','Sr':'Sr_sv','Ta':'Ta_pv','Tb':'Tb_3','Tc':'Tc_pv', 'Te':'Te','Th':'Th','Ti':'Ti_pv','Tl':'Tl_d','Tm':'Tm_3','U':'U','V':'V_pv','W':'W_pv','Xe':'Xe', 'Y':'Y_sv','Yb':'Yb_2','Zn':'Zn','Zr':'Zr_sv' } # ------------------------------------------------------------------ # 3. Utility: Oxidation state / POTCAR # ------------------------------------------------------------------ def add_oxi(structure): """Add oxidation state to site_properties using BVAnalyzer.""" try: bva = BVAnalyzer() oxi_structure = bva.get_oxi_state_decorated_structure(structure) print("Successfully added oxidation state information") return oxi_structure except Exception as error: print(f"Failed to add oxidation state via BVAnalyzer: {error}") return structure # ------------------------------------------------------------------ # 4. Merge and generate POTCAR files in order of elements contained in the structure. Preferentially use recommended POTCAR # ------------------------------------------------------------------ # ★ Modified: function now directly receives vasp_set object def generate_potcar(structure, potdir): """ Merges and generates POTCAR file in the order of elements in the POSCAR contained in vasp_set. Searches for recommended POTCAR first, and if not found, searches for all possible alternative POTCARs. """ # ★ Modified: use vasp_set.poscar.site_symbols instead of structure.composition.elements # This ensures the actual element order of the POSCAR file elems = vasp_set.poscar.site_symbols paths, used_potcars = [], [] for el in elems: found = False # 1. Search for recommended POTCAR if el in RECOMMENDED_POTCARS: tag = RECOMMENDED_POTCARS[el] # Check recommended POTCAR directory with various case formats for name_case in (tag, tag.upper(), tag.lower(), tag.capitalize()): p = os.path.join(potdir, name_case, "POTCAR") if os.path.exists(p): paths.append(p) used_potcars.append(f"{el}:{tag}") print(f"✓ {el}: using recommended POTCAR '{tag}'") found = True break # If recommended POTCAR is specified but the actual file does not exist, add to list if not found and (el, tag) not in not_found_recommended_elements: not_found_recommended_elements.append((el, tag)) # 2. If recommended POTCAR was not found or not in the list, search for alternative POTCAR if not found: print(f"✗ {el}: recommended POTCAR not found, searching for alternative POTCAR...") # Expand search scope to include all cases like capitalize, upper, lower as before search_suffixes = ["", "_pv", "_sv", "_d", "_h", "_s"] candidate_paths = [] for suffix in search_suffixes: for case_fn in (lambda x: x, str.capitalize, str.upper, str.lower): dir_name = case_fn(el) + suffix candidate_paths.append(os.path.join(potdir, dir_name, "POTCAR")) for p in candidate_paths: if os.path.exists(p): dir_name = os.path.basename(os.path.dirname(p)) paths.append(p) used_potcars.append(f"{el}:{dir_name}") print(f"✓ {el}: using alternative POTCAR '{dir_name}'") found = True break if not found: print(f"✗ Error: could not find POTCAR for {el}!") return False # 3. Merge POTCAR files try: subprocess.run(f"cat {' '.join(paths)} > POTCAR", shell=True, check=True) print(f"POTCAR generation complete: {', '.join(used_potcars)}") return True except subprocess.CalledProcessError as e: print(f"POTCAR merge error: {e}") return False # ------------------------------------------------------------------ # 5. User default INCAR / KPOINTS options # ------------------------------------------------------------------ user_incar_settings = { # "KPAR": 4, # same as nk 4 option in QE # "EDIFF": 1e-6, # electronic convergence criterion (use if needed) # "ISMEAR": 0, # smearing type # "SIGMA": 0.05, # smearing width # "IBRION": -1, # no ion update # "NSW": 0, # number of ion relaxation steps (0 means no movement) # "LREAL": Auto, # projection operator setting (False: reciprocal, True: real) # "LWAVE": False, # do not save WAVECAR # "LCHARG": False, # do not save CHGCAR # "ALGO": "Fast", # ALGO setting # Additional options can be extended! # Add values to use for r2SCAN below # "LVHAR": False, # whether to record local potential in LOCPOT; if LVHAR=True, LVTOT is changed to False # "LVTOT": False, # whether to record Total local potential in LOCPOT } user_kpoints_settings = { # "reciprocal_density": 100, # "kpts": [[4, 4, 4]], # directly specify KPOINTS grid (use when uncommented) # "gamma_centered": True, # whether gamma-centered # "shift": [0.0, 0.0, 0.0], # k-point shift } # ------------------------------------------------------------------ # 6. MP search / input generation loop # ------------------------------------------------------------------ with MPRester(API_KEY) as mpr, open("structure.txt","w") as file: # --- Materials Project condition-based search settings --- structures = mpr.materials.summary.search( # material_ids = ["mp-13", "mp-23", "mp-30", "mp-33", "mp-54", "mp-1042545", "mp-20722"]) # can directly specify multiple structures by MPID # elements=["Co"], # required included atom exclude_elements=["Co"], # exclude alloys containing Co chemsys="*-*", # specify chemical formula format # formula_anonymous = "AB", # search by chemical formula composition ratio family # band_gap = (0.5, 1.0), # Band Gap min, max values energy_above_hull=(0, 0.02), # energy stability (eV) num_sites=(0, 20), # number of atoms per cell is_metal=True, # metallic materials only ) file.write(f"Number of documents: {len(structures)}\n\n") if not structures: file.write("No results found.\n") raise SystemExit("No search results found") # Determine total number of materials | added 2025.09.01 total_materials = len(structures) print(f"Found {total_materials} materials. Starting individual data download.") # Save script execution location (use os.chdir(dirname), os.chdir(home) to move to mp-id folder and back to parent folder) home = os.getcwd() for idx, structure_summary in enumerate(structures,1): try: # --- 1. Collect data --- material_id = structure_summary.material_id # MPID formula = structure_summary.formula_pretty # chemical formula energy_above_hull = structure_summary.energy_above_hull # Hull energy # is_metal = structure_summary.is_metal # whether metal # is_stable = structure_summary.is_stable # whether stable phase # crystal = structure_summary.structure # crystal structure symmetry = structure_summary.symmetry # sym structure mag_order = structure_summary.ordering # magnetism such as Ferro, Anti-Ferro etc. mag_total = structure_summary.total_magnetization # total magnetization mag_norm_fu = structure_summary.total_magnetization_normalized_formula_units # Fetch structure, energy, magnetism information origin from origin (task_id search) structure_task_id = None energy_task_id = None magnetism_task_id = None if structure_summary.origins: for origin in structure_summary.origins: if origin.name == "structure": structure_task_id = origin.task_id elif origin.name == "energy": energy_task_id = origin.task_id elif origin.name == "magnetism": magnetism_task_id = origin.task_id # Efficiently look up run_type of all extracted Task IDs with a single API call structure_run_type = "N/A" energy_run_type = "N/A" magnetism_run_type = "N/A" # Make a deduplicated list of IDs to query (excluding None) ids_to_query = list(set(filter(None, [structure_task_id, energy_task_id, magnetism_task_id]))) if ids_to_query: try: # Fetch information of all tasks with a single API call task_docs = mpr.tasks.search(task_ids=ids_to_query, fields=["task_id", "run_type"]) # Convert result into a dictionary (map) for fast lookup run_type_map = {doc.task_id: doc.run_type for doc in task_docs} # Find and assign each run_type from the map if structure_task_id: structure_run_type = run_type_map.get(structure_task_id, "Not Found") if energy_task_id: energy_run_type = run_type_map.get(energy_task_id, "Not Found") if magnetism_task_id: magnetism_run_type = run_type_map.get(magnetism_task_id, "Not Found") except IndexError: print(f" - Warning: Could not retrieve run_type for some Task IDs: {ids_to_query}") # Print progress status when processing each material (print total material count and the number of currently processed material) | added 2025.09.01 print(f"[{idx} / {total_materials}] Processing: {formula} ({material_id})") # --- 2. Get structure and preprocess --- structure = mpr.get_structure_by_material_id(material_id, final=use_final_structure) if isinstance(structure, list): structure = structure[0] primitive_cell = structure.get_primitive_structure() # Add oxidation state primitive_cell = add_oxi(primitive_cell) # Get Full Formula full_formula = primitive_cell.composition.formula # Full formula full_formula_no_spaces = "".join(full_formula.split()) # remove spaces # Get Cohesive energy cohesive_energy = mpr.get_cohesive_energy( material_ids=[material_id], normalization="atom" # normalize by atom count )[material_id] # result is returned as a dictionary, so access via material_id # Generate VASP input (for each directory) directory_name =f"{''.join(primitive_cell.composition.formula.split())}_{material_id}" os.makedirs(directory_name, exist_ok = True) os.chdir(directory_name) # --- 3. INCAR / MAGMOM settings--- incar = user_incar_settings.copy() if USE_CONFIG_MAGMOM: # 1) CONFIG_MAGMOM forced mode: use only CONFIG_MAGMOM values unconditionally # if magmom exists in site_properties, delete it # (remove pymatgen auto magmom) if "magmom" in primitive_cell.site_properties: del primitive_cell.site_properties["magmom"] # if "magmom" in conventional_cell.site_properties: # << for Conventional structure # del conventional_cell.site_properties["magmom"] # << for Conventional structure # based on species_string, use CONFIG_MAGMOM.get(..., 0.0) mag_dict = {} for site in primitive_cell: key = site.species_string # e.g., "Fe3+" mag_dict[key] = CONFIG_MAGMOM.get(key, 0.0) incar["MAGMOM"] = mag_dict mag_info_for_log = f"config: {mag_dict}" print(f"[{material_id}] Applied MAGMOM(config) dict → {mag_dict}") else: # 2) site_property mode: use only site_properties["magmom"] if "magmom" not in primitive_cell.site_properties: # if site_property does not exist, warn the user msg = f"[{material_id}] ERROR: site_properties['magmom'] not found!" print(msg) mag_info_for_log = "ERROR: no site_properties['magmom']" else: sp_mag = primitive_cell.site_properties["magmom"] incar["MAGMOM"] = sp_mag mag_info_for_log = f"site_properties: {sp_mag}" print(f"[{material_id}] Applied MAGMOM(site_properties) → {sp_mag}") # --- 4. Generate VASP input files --- vasp_set = MPStaticSet( primitive_cell, user_incar_settings=user_incar_settings, user_kpoints_settings=user_kpoints_settings ) # Reflect incar value determined above into MPStaticSet internally vasp_set.incar.update(incar) vasp_set.poscar.write_file("POSCAR") vasp_set.incar.write_file("INCAR") # Added 2025.09.24: issue where MPScanStaticSet uses KSPACING and does not generate KPOINTS file # Added exception handling in case KPOINTS is not generated if vasp_set.kpoints: # same as if vasp_set.kopints is not None vasp_set.kpoints.write_file("KPOINTS") # write KPOINTS file if it exists else: kpoint_error_msg = f"[{material_id}] Warning: KPOINTS not automatically generated" print(kpoint_error_msg) # print error message if KPOINTS generation fails file.write(kpoint_error_msg + "\n") # keep log for tracking if KPOINTS generation fails # Generate POTCAR # ★ Modified: pass vasp_set object instead of primitive_cell if not generate_potcar(vasp_set, POTCAR_DIR): print(f"[{material_id}] POTCAR generation failed") # --- 5. Save detailed log --- file.write(f"--- Document {idx} ---\n") file.write(f"Material ID : {material_id}\n") file.write(f"Full Formula : {full_formula}\n") file.write(f"Directory : {directory_name}\n") file.write(f"Structure from : {structure_run_type}\n") file.write(f"Energy from : {energy_run_type}\n") file.write(f"Magnetism from : {magnetism_run_type}\n") file.write(f"Energy Above Hull : {energy_above_hull}\n") # file.write(f"Chemical Formula: {formula}\n") file.write(f"Cohesive Energy : {cohesive_energy} eV/atom\n") # file.write(f"Is Metal: {is_metal}\n") # file.write(f"Stable: {is_stable}\n") # file.write(f"Crystal: {structure}\n") file.write(f"Symmetry : {symmetry}\n") file.write(f"MAGMOM mode : {mag_info_for_log}\n") file.write(f"Magnetic Ordering : {mag_order}\n") file.write(f"Total Magnetization: {mag_total}\n") file.write(f"Total Magnetization Normalized Formula Unit: {mag_norm_fu}\n") file.write(f"Primitive Structure:\n{primitive_cell}\n") # file.write(f"Primitive Structure:\n{conventional_cell}\n") file.write('-' * 40 + "\n\n") print(f"[{material_id}] Complete → {directory_name}\n") except Exception as e: # Log error and continue processing when an error occurs print(f"[{idx}, {material_id if 'material_id' in locals() else 'N/A'}] Error occurred during processing: {e}") traceback.print_exc() file.write(f"--- Document {idx} ({material_id if 'material_id' in locals() else 'N/A'}) ---\nError occurred: {e}\n\n") finally: # Always return to the original directory regardless of success/failure of the try block os.chdir(home) # After all processing is finished, print the list of elements for which recommended POTCAR could not be found if not_found_recommended_elements: print("\n===== List of elements for which recommended POTCAR was not found =====") for element,tag in not_found_recommended_elements: print(f"{element}: recommended POTCAR '{tag}' does not exist – using alternative POTCAR") print("="*50) ```