diff --git a/sgrcsp/main.py b/sgrcsp/main.py index 2c95204..f6523d5 100644 --- a/sgrcsp/main.py +++ b/sgrcsp/main.py @@ -8,16 +8,16 @@ def print_logo(): print( """ - Package name is not decided yet |\__/,| (` _.|o o |_ ) ) -(((---(((-------- + Space Group Restricted Crystal Structure Prediction """ ) - print("-------------------(version", __version__, ")--------------------\n") + print("-------------------(version", __version__,")--------------------\n") print("A Python package for symmerty restricted crystal structure prediction") - print("The source code is available at ..") - print("Developed by ColdSnaaap \n\n") + print("The source code is available at https://github.com/ColdSnaap/sgrcsp.git") + print("Developed by ColdSnaap \n\n") sys.stdout.flush() diff --git a/sgrcsp/optimization.py b/sgrcsp/optimization.py index 009d5ca..fcb76fa 100644 --- a/sgrcsp/optimization.py +++ b/sgrcsp/optimization.py @@ -192,6 +192,7 @@ def simulated_annealing( print("-------------------------------------------------") print(f"Current state: {current_energy} eV, {current_energy_atom} eV/atom") print(f"New state: {next_energy} eV, {next_energy_atom} eV/atom") + sys.stdout.flush() # Determine if we should accept the new state accept = False diff --git a/sgrcsp/routine.py b/sgrcsp/routine.py index 311460b..83794e1 100644 --- a/sgrcsp/routine.py +++ b/sgrcsp/routine.py @@ -99,11 +99,10 @@ def read_continue_config(): config["count"] = int(config["count"]) config["best_energy"] = float(config["best_energy"]) config["accept_list"] = ast.literal_eval(config["accept_list"]) - current_state = Poscar.from_file(continue_path+"/POSCAR_current") + # current_state = Poscar.from_file(continue_path+"/POSCAR_current") best_state = Poscar.from_file(continue_path+"/POSCAR_best") - return config["step_current"], config["accept_list"], config["count"], config["best_energy"], \ - current_state, best_state + return config["step_current"], config["accept_list"], config["count"], config["best_energy"], best_state class Routine: @@ -142,34 +141,52 @@ def __init__(self) -> None: print(f"MaxStep: {self.max_step}\n") sys.stdout.flush() # generate initial structure - if self.initial_relax: - while True: - try: - self.structure = structure_generation( + file_path = os.getcwd() + "/Input/initial.cif" + if not self.job_contine: + if os.path.isfile(file_path): + self.structure = structure_generation( self.sg, self.mol_list, self.mol_number, self.sites, factor=self.factor ) - print("Relaxing the initial structure ...") - sys.stdout.flush() - cal = CalculateEnergy(self.structure) - cal_struc = cal.energy_calculate("VASP", True, self.sub_command)[0] - self.structure = structure_type_converter(cal_struc, "pyxtal") - print("Initial structure generated") - sys.stdout.flush() - break - except ReadSeedError: - print("Molecule not intact, trying to generate a new structure to relax") - else: - self.structure = structure_generation( - self.sg, - self.mol_list, - self.mol_number, - self.sites, - factor=self.factor - ) + else: + while True: + try: + self.structure = structure_generation( + self.sg, + self.mol_list, + self.mol_number, + self.sites, + factor=self.factor + ) + print("Relaxing the initial structure ...") + sys.stdout.flush() + cal = CalculateEnergy(self.structure) + cal_struc = cal.energy_calculate("VASP", True, self.sub_command)[0] + self.structure = structure_type_converter(cal_struc, "pyxtal") + print("Initial structure generated") + sys.stdout.flush() + break + except ReadSeedError: + print("Molecule not intact, trying to generate a new structure to relax") + + elif self.job_contine: + print("Continue ...") + file_path = os.getcwd() + "/StrucBackup/POSCAR" + self.structure = structure_type_converter(file_path, "pyxtal") + print("Backup structure read") + sys.stdout.flush() + + # write backup + backup_dir = os.getcwd() + "/StrucBackup" + if not os.path.isdir(backup_dir): + os.mkdir(backup_dir) + struc_py = structure_type_converter(self.structure, "pymatgen") + poscar = Poscar(struc_py) + poscar.write_file(f"{backup_dir}/POSCAR") + def log( self, @@ -188,6 +205,7 @@ def log( result_dir = os.getcwd() + "/Result" struc_log = os.getcwd() + "/Result/BestStrucsList" struc_file = os.getcwd() + "/Result/BestStrucs" + backup_dir = os.getcwd() + "/StrucBackup" sga = SpacegroupAnalyzer(structure, symprec=0.1) refined_structure = sga.get_refined_structure() @@ -199,8 +217,14 @@ def log( # raise NameError("Result directory not exist") if not os.path.isdir(struc_file): os.makedirs(struc_file) + if not os.path.isdir(backup_dir): + os.mkdir(backup_dir) name_en = abs(energy) + # backup structure + poscar = Poscar(structure) + poscar.write_file(f"{backup_dir}/POSCAR") + if not os.path.isfile(f"{struc_file}/{name_en}.cif"): # refined_structure.to(f'{struc_file}/{name_en}.cif', fmt="cif") cif_writer = CifWriter(refined_structure, symprec=0.1) @@ -308,32 +332,44 @@ def routine3(self): step_current = 0 accept_list = None - current_state = self.structure count = None best_energy = None best_state = None rate_scale = None else: - step_current, accept_list, count, best_energy, current_state, best_state = \ - read_continue_config() + continue_path = os.getcwd() + "/Result/Continue" + if os.path.isdir(continue_path): + step_current, accept_list, count, best_energy, best_state = read_continue_config() + else: + step_current = 0 + accept_list = None + count = None + best_energy = None + best_state = None + rate_scale = None + current_state = self.structure + for step in range(self.loop_number): - while True: - try: - current_state_ori = structure_type_converter(current_state, "pymatgen") - current_state = copy.deepcopy(current_state_ori) - cal1 = CalculateEnergy(current_state) - cal = cal1.energy_calculate("VASP", True, self.sub_command) - current_state, relax_energy, relax_energy_atom = cal - current_state = structure_type_converter(current_state, "pyxtal") - break - except ReadSeedError: - print("Molecule not intact, trying to generate a new structure to relax") - sys.stdout.flush() - for i in range(10): - current_state_ori = apply_perturbation(current_state_ori) - + try: + current_state_ori = structure_type_converter(current_state, "pymatgen") + current_state = copy.deepcopy(current_state_ori) + cal1 = CalculateEnergy(current_state) + cal = cal1.energy_calculate("VASP", True, self.sub_command) + current_state, relax_energy, relax_energy_atom = cal + current_state = structure_type_converter(current_state, "pyxtal") + except ReadSeedError: + print("Molecule not intact, use backup structure") + sys.stdout.flush() + backup_file = os.getcwd()+"/StrucBackup/POSCAR" + current_state_ori = Structure.from_file(backup_file) + current_state = copy.deepcopy(current_state_ori) + cal1 = CalculateEnergy(current_state) + cal = cal1.energy_calculate("VASP", True, self.sub_command) + current_state, relax_energy, relax_energy_atom = cal + current_state = structure_type_converter(current_state, "pyxtal") + self.log(step_current, current_state, relax_energy) # small perturbation diff --git a/sgrcsp/structure.py b/sgrcsp/structure.py index b2f42b7..bace755 100644 --- a/sgrcsp/structure.py +++ b/sgrcsp/structure.py @@ -2,6 +2,8 @@ import copy import sys import random +import numpy as np +import shutil from itertools import combinations from pathlib import Path from readconfig import ReadConfig @@ -9,6 +11,7 @@ from pymatgen.core import Structure from pymatgen.io.vasp import Poscar from pymatgen.io.cif import CifParser +from pyxtal.msg import ReadSeedError class Perturbation: @@ -21,20 +24,33 @@ def __init__(self, structure): sys.stdout.flush() + def intact_mol_check(self, structure): + struc_py = structure_type_converter(structure, "pymatgen") + try: + struc_pyxtal = structure_type_converter(struc_py, "pyxtal") + return True + except ReadSeedError: + return False + + def constant(self, d_coor, d_rot, d_lat, trial_max=1000): print("Constant perturbation") - print(f"lat, coor, rot: {d_lat}, {d_coor}, {d_rot}") + # print(f"lat, coor, rot: {d_lat}, {d_coor}, {d_rot}") sys.stdout.flush() min_dis_struc = 0.001 count = 0 + mol_intact = False - while min_dis_struc < self.min_dis_allow: + while min_dis_struc < self.min_dis_allow or mol_intact == False: count += 1 struc_trial = copy.deepcopy(self.struc_pyxtal) struc_trial.apply_perturbation(d_lat, d_coor, d_rot) + # check minimum distance between atoms struc_py = structure_type_converter(struc_trial, "pymatgen") - min_dis_struc = distance_check(struc_py) + min_dis_struc = distance_min(struc_py) + # check if molecules are intact + mol_intact = self.intact_mol_check(struc_trial) if count == trial_max: print(f"Could not find good structure") @@ -56,32 +72,48 @@ def constant(self, d_coor, d_rot, d_lat, trial_max=1000): def uniform(self, d_coor1, d_rot1, d_lat1, trial_max=1000): print("Uniform perturbation") - print(f"lat, coor, rot: 0-{d_lat1}, 0-{d_coor1}, 0-{d_rot1}") + # print(f"lat, coor, rot: 0-{d_lat1}, 0-{d_coor1}, 0-{d_rot1}") sys.stdout.flush() min_dis_struc = 0.001 + # min_dis_mol = 0.001 count = 0 - + # count2 = 0 + # mol_intact = False while min_dis_struc < self.min_dis_allow: + # while min_dis_struc < self.min_dis_allow or \ + # min_dis_mol < 2.0*self.min_dis_allow: + # mol_intact == False: count += 1 struc_trial = copy.deepcopy(self.struc_pyxtal) d_lat = round(random.uniform(0, d_lat1), 3) d_coor = round(random.uniform(0, d_coor1), 3) d_rot = round(random.uniform(0, d_rot1), 3) - print(f"lat, coor, rot: {d_lat}, {d_coor}, {d_rot}") - sys.stdout.flush() + # print(f"lat, coor, rot: {d_lat}, {d_coor}, {d_rot}") + # sys.stdout.flush() struc_trial.apply_perturbation(d_lat, d_coor, d_rot) + # check minimum distance between atoms struc_py = structure_type_converter(struc_trial, "pymatgen") - min_dis_struc = distance_check(struc_py) - + min_dis_struc = distance_min(struc_py) + # try: + # min_dis_mol = distance_min_mol(struc_trial) + # except ReadSeedError: + # x_write = Poscar(struc_py) + # x_write.write_file("POSCAR_error") + # check if molecules are intact + # mol_intact = self.intact_mol_check(struc_trial) + # print(f"min_dis_struc:{min_dis_struc}") + # print(f"min_dis_mol:{min_dis_mol}") if count == trial_max: + # count2 += 1 print(f"Could not find good structure after {trial_max} trail") print(f"d_lat, d_coor, d_rot: {d_lat}, {d_coor}, {d_rot}") sys.stdout.flush() - count = 0 + # if count2 == 5: + # raise ValueError("No good structure found after {count2} loop attemps") print(f"Trail: {count}") sys.stdout.flush() @@ -143,7 +175,7 @@ def structure_type_converter(structure, target_type, molecule=True): if file_type == "cif": parser = CifParser(structure) structure = parser.get_structures()[0] - elif file_type == "POSCAR": + elif file_type == "POSCAR" or "CONTCAR": poscar = Poscar.from_file(structure) structure = poscar.structure @@ -178,12 +210,12 @@ def structure_type_converter(structure, target_type, molecule=True): raise TypeError(f"molecule={molecule} not support yet") elif isinstance(structure, pyxtal): return structure - + else: raise NameError(f"Does not support target_type={target_type} yet") -def distance_check(structure) -> float: +def distance_min(structure) -> float: """ Return the minimum distance between atoms in a structure. """ @@ -191,7 +223,6 @@ def distance_check(structure) -> float: structure = structure_type_converter(structure, "pymatgen") min_distance = float('inf') # Initialize with infinity - # Generate all unique pairs of sites for site1, site2 in combinations(structure.sites, 2): # Calculate distance between the pair of sites directly @@ -203,6 +234,36 @@ def distance_check(structure) -> float: return min_distance +def distance_min_mol(structure) -> bool: + + struc_pymatgen = structure_type_converter(structure, "pymatgen") + struc_pyxtal = structure_type_converter(structure, "pyxtal") + + # get the elements in molecular -> list + rigid_elem_list = [] + for i, site in enumerate(struc_pyxtal.mol_sites): + if len(site.get_coords_and_species()[0]) != 1: + rigid_elem_list += site.get_coords_and_species()[1] + centroid = np.mean(site.get_coords_and_species()[0], axis =0) + struc_pymatgen.append("H", centroid) + seen = set() + unique_elements = [x for x in rigid_elem_list if not (x in seen or seen.add(x))] + + # structure with only single atoms and "H" as molecules + struc_pymatgen.remove_species(unique_elements) + # get the min distance between atom and molecules + # make sure atom is not "inside" molecules + min_distance = float('inf') + for site1, site2 in combinations(struc_pymatgen.sites, 2): + if "H" in [str(site1.specie), str(site2.specie)]: + # print(site1, site2) + distance = site1.distance(site2) + if distance < min_distance: + min_distance = distance + + return min_distance + + def structure_generation( sg, mol_list, @@ -210,19 +271,26 @@ def structure_generation( sites, factor=1.0 ) -> pyxtal: - print("Generating structure ...") - sys.stdout.flush() - while True: - try: - structure = pyxtal(molecular=True) - structure.from_random(3, sg, mol_list, mol_number, sites=sites, factor=factor) - break - except RuntimeError as e: - print(f"RuntimeError: {e}.") - finally: - factor += 0.05 - print(f"Increasing volume, volume factor: {factor}") - sys.stdout.flush() - - print("Structure generated\n") - return structure + file_path = os.getcwd() + "/Input/initial.cif" + if os.path.isfile(file_path): + print("Initial structure found") + structure = structure_type_converter(file_path, "pyxtal") + return structure + else: + print("Generating structure ...") + sys.stdout.flush() + while True: + try: + structure = pyxtal(molecular=True) + structure.from_random(3, sg, mol_list, mol_number, sites=sites, factor=factor) + break + except RuntimeError as e: + print(f"RuntimeError: {e}.") + finally: + factor += 0.05 + print(f"Increasing volume, volume factor: {factor}") + sys.stdout.flush() + + print("Structure generated\n") + return structure + diff --git a/sgrcsp/test.py b/sgrcsp/test.py index 4b40378..30386d2 100644 --- a/sgrcsp/test.py +++ b/sgrcsp/test.py @@ -22,11 +22,12 @@ from caltool import CalculateEnergy import subprocess from structure import structure_type_converter -from structure import distance_check +from structure import distance_min_mol from combination import WyckoffCombinations from analysis import vasp_ml_compare from analysis import check_trajectory from structure import Perturbation +from pymatgen.io.cif import CifParser config = ReadConfig() @@ -37,7 +38,7 @@ x = WyckoffCombinations(mol_list) for i in sg: - yy = x.mol_ratio_comb_list_sg([1, 3], i, 3) + yy = x.mol_ratio_comb_list_sg([1, 1, 1], i, 4) # per_dict = { # "d_lat": 0.0, @@ -80,8 +81,37 @@ # py = structure_type_converter(x, "pymatgen") -# poscar = Poscar.from_file("/Users/qizhang/Desktop/paper_new/PSLi/BESTgatheredPOSCARS_order") +# poscar = Poscar.from_file(os.getcwd()+"/99.77045.cif") # structure = poscar.structure + +# strcu_pyxtal = structure_type_converter(os.getcwd()+"/POSCAR", "pyxtal") +# strcu_pyxtal = structure_type_converter(os.getcwd()+"/99.77045.cif", "pyxtal") +# print(strcu_pyxtal) +# x = distance_min_mol(strcu_pyxtal) + + +# pertub = Perturbation(strcu_pyxtal) +# strcu_pyxtal = pertub.uniform(d_coor1=0.2, d_rot1=2.0, d_lat1=0.0) +# print(strcu_pyxtal.mol_sites[0].get_coords_and_species()) +# print("\n") +# coor1 = strcu_pyxtal.mol_sites[2].get_coords_and_species()[0] +# coor2 = strcu_pyxtal.mol_sites[3].get_coords_and_species()[0] +# print(strcu_pyxtal.mol_sites[0].get_distances(coor1, coor2)) +# print(str(strcu_pyxtal.mol_sites[2])) + +# for i in range(200): +# print(i) +# pertub = Perturbation(strcu_pyxtal) +# strcu_pyxtal = pertub.uniform(d_coor1=0.2, d_rot1=2.0, d_lat1=0.0) +# print("\n") + +# x_py = structure_type_converter(strcu_pyxtal, "pymatgen") +# x_write = Poscar(x_py) +# x_write.write_file("POSCAR_result") + + + +# print(x) # cif = CifWriter(py, 0.1) # cif.write_file("test.cif") diff --git a/sgrcsp/test2.py b/sgrcsp/test2.py new file mode 100644 index 0000000..74d4e49 --- /dev/null +++ b/sgrcsp/test2.py @@ -0,0 +1,19 @@ +import os +from readconfig import ReadConfig +from pymatgen.io.vasp import Poscar +from structure import apply_perturbation +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer +from pymatgen.io.cif import CifWriter + +file = "/Users/qizhang/Desktop/phonon_check/structure_check/POSCAR3" +poscar = Poscar.from_file(file) +structure = poscar.structure + +sga = SpacegroupAnalyzer(structure, symprec=0.1) +print(sga.get_space_group_number()) + +refined_structure = sga.get_refined_structure() + + +cif_writer = CifWriter(refined_structure, symprec=0.1) +cif_writer.write_file("/Users/qizhang/Desktop/phonon_check/structure_check//poscar3.cif") \ No newline at end of file