# -*- coding: utf-8 -*- ######################################################### #Instructions #Set up a template .inp file for the parent in the correct format with correct flags (see example .inp) #Change directories in lines below #Change to the right directory #Press run ######################################################### # JWL - updated 2016 # JSOE - updated 18/9/2020 for new isodistort format with space group numbers and origin line import os import re import math import glob import subprocess import xlsxwriter from shutil import copyfile ######################################################### # Failsafe Function ######################################################### def dead(reason): # indicates script closure with reason print reason raise SystemExit ######################################################### # User required info - Remember to cd to working directory before running # Delete any files in best_inp_files before running ######################################################### working_directory = r'C:\data\covid19_work\symmetry_descent_james' # Enter directory where trial subgroups .strs are stored topas_directory = r'C:\software\topasa_v7' # Enter directory topas where tc.exe is located thermal_prms = {"W":1.5, "O":3} # Do not have a dictionary entry if beq not to be changed seed_file = working_directory + r'/WO3_850c_start.inp' # pre setup .inp file to be automatically edited for each subgroup result_file = working_directory + r'/results_summary.txt' # New file created any old one with the same name is overwritten basefile = r"topas" # prefix to subgroups .str files generated by ISODISTORT, could be in another directory #subgroup_numbers_to_run = [1,3,58,60,72] # List of subgroup numbers (suffix of ISODISTORT output) to run subgroup_numbers_to_run = [x for x in range(3,73) if x!=31] # List of subgroup numbers (suffix of ISODISTORT output) to run number_of_steps = 2 # number of different instructions sets in the topas.inp to run for each subgroup i.e. #ifdef run_1, #ifdef run_2, etc ######################################################### # Delete any old start.inp and best.inp files in subdirectories ######################################################### best_list = glob.glob(os.path.join(working_directory,'best_inp_files','*.inp')) start_list = glob.glob(os.path.join(working_directory,'start_inp_files','*.inp')) delete_list = best_list + start_list for filePath in delete_list: try: os.remove(filePath) except: print("Error while deleting file : ", filePath) ######################################################### # Search compilers (to allow editing of ISODISTORT .str files) # Assumes fixed formats for the .str files as below ######################################################### lpa_line = re.compile(r"\s+a (.+?)+-?") # anything then string starting with a then 4 spaces then any number (positive or negative) lpb_line = re.compile(r"\s+b (.+?)+-?") lpc_line = re.compile(r"\s+c (.+?)+-?") be_line = re.compile(r"\s+be (.+?)+-?\d+?") al_line = re.compile(r"\s+al (.+?)+-?\d+?") ga_line = re.compile(r"\s+ga (.+?)+-?\d+?") Space_group_line = re.compile(r"\s+space_group\s+?") mode_site_line = re.compile("\t\tsite") # looks for two tabs then site displacive_line = re.compile(r"\s+prm\s+!a\d+\s+-?\d+\.?\d+'?") origin_site_line =re.compile(" origin ") #looks for origin line and deletes ######################################################### # Main Inp File Creation Function ######################################################### def inp_maker(basefile, file_number, thermal_prms): # this function converts all str files into input files in the seed format new_inp_file = (basefile + "%04d" % (file_number) + '.inp') structural_file = (basefile + "%04d" % (file_number) + '.str') new_inp_file_stream = open(new_inp_file,'w') structural_file_stream = open(structural_file, 'r') structure_lines = structural_file_stream.readlines() structural_file_stream.close() #jsoe to get rid of the origin line here displacive_modes_in_this_str = 0 for line in structure_lines: if displacive_line.search(line): displacive_modes_in_this_str += 1 space_group, space_group_name, crystal_class = get_space_group(structure_lines,subgroup) section_flag = 0 for i, line in enumerate(open(seed_file).readlines()): if re.search("""'!End_Lattice_Prm_Section!""", line) or re.search("""'!End_Insert_Sites_Section!""", line): section_flag = 0 elif section_flag == 1: continue elif re.search("""'!Lattice_Prm_Section!""", line): refine_lps(new_inp_file_stream, structure_lines, crystal_class,space_group,space_group_name) section_flag = 1 elif re.search("""'!Insert_Sites_Section!""", line): add_structural_info(new_inp_file_stream, structure_lines, thermal_prms) section_flag = 1 elif re.search("""'!Insert_Penalties_Section!""", line): for i in range(0, displacive_modes_in_this_str): new_inp_file_stream.write("aP(a%d)\n" % (i+1)) else: new_inp_file_stream.write(line) #print(line) new_inp_file_stream.close() ########################################################### # Supplementary Editing Functions ########################################################### def add_structural_info(new_inp_file_stream, structure_lines, thermal_prms): # Add mode definitions and structural sites for i, line in enumerate(structure_lines[13:]): site = None if displacive_line.search(line): new_inp_file_stream.write(line.replace("prm !a", "mode_onoff(a").replace("0.00000", ", 0.00000) drand '")) elif mode_site_line.search(line) and not origin_site_line.search(line): #elif mode_site_line.search(line): for atom in thermal_prms: if ("site %s" % (atom)) in line: site = atom break new_inp_file_stream.write(line.replace("beq 0.0", "beq %f" % (thermal_prms[site]))) else: new_inp_file_stream.write(line) #jsoe to just work with space group number from str file here so can ignore the search of sgcom5 #will just be able to say k=int(space_group) #then get the class ranges from isodistort_tree.py def get_space_group(structure_lines,subgroup): # determine space group and crystal class (in order to determine which lps to refine) for i,line in enumerate(structure_lines): if Space_group_line.search(line): space_group = line.split()[1] space_group_name=re.sub("'","",structure_lines[i-1].split()[0]) k=int(space_group) print "Subgroup number: "+str(subgroup)+" Space group number: "+str(k)+" Space group name: ",space_group_name if k==1 or k==2: print 'Triclinic space group detected!' crystal_class = "triclinic" elif k>=3 and k<=15: print 'Monoclinic space group detected!' crystal_class = "monoclinic" elif k>=16 and k<=74: print 'Orthorhombic space group detected!' crystal_class = "orthorhombic" elif k>=75 and k<=142: print 'Tetragonal space group detected!' crystal_class = "tetragonal" elif k>=143 and k<=166: print 'Trigonal space group detected! Warning: assuming hexagonal cell setting' crystal_class = "trigonal" elif k>=167 and k<=194: print 'Hexagonal space group detected!' crystal_class = "hexagonal" elif k>=195 and k<=230: print 'Cubic space group detected!' crystal_class = "cubic" else: dead('space group not assigned') return space_group, space_group_name, crystal_class def refine_lps(new_inp_file_stream, structure_lines, crystal_class,space_group,space_group_name): # Add lattice prms with appropriate refinement flags if crystal_class == "triclinic": refinement_flags = ["lpa", "lpb", "lpc", "alph", "bet", "gam"] elif crystal_class == "monoclinic": refinement_flags = ["lpa", "lpb", "lpc", "!alph", "bet", "!gam"] elif crystal_class == "orthorhombic": refinement_flags = ["lpa", "lpb", "lpc", "!alph", "!bet", "!gam"] elif crystal_class == "tetragonal": refinement_flags = ["lpa", "lpa", "lpc", "!alph", "!bet", "!gam"] elif crystal_class == "trigonal": refinement_flags = ["lpa", "lpa", "lpc", "!alph", "!bet", "!gam"] # In ISODISTORT it appears all trigonal lattices are in a hexagonal setting elif crystal_class == "trigonal_hexagonal_setting": refinement_flags = ["lpa", "lpa", "lpc", "!alph", "!bet", "!gam"] elif crystal_class == "hexagonal": refinement_flags = ["lpa", "lpa", "lpc", "!alph", "!bet", "!gam"] elif crystal_class == "cubic": refinement_flags = ["lpa", "lpa", "lpa", "!alph", "!bet", "!gam"] #new_inp_file_stream.write((structure_lines[4])) new_inp_file_stream.write("\t\tspace_group "+space_group_name+" '"+space_group+"\n") new_inp_file_stream.write((structure_lines[6].replace("\n", "").replace("a", "prm %s_dummy" % (refinement_flags[0]))) + " del = 0.0001 Val; val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*a_start; min = 0.95*a_start; max = 1.05*a_start;\n") if refinement_flags[1] != "lpa": new_inp_file_stream.write((structure_lines[7].replace("\n", "").replace("b", "prm %s_dummy" % (refinement_flags[1]))) + " del = 0.0001 Val; val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*b_start; min = 0.95*b_start; max = 1.05*b_start;\n") if refinement_flags[2] != "lpa": new_inp_file_stream.write((structure_lines[8].replace("\n", "").replace("c", "prm %s_dummy" % (refinement_flags[2]))) + " del = 0.0001 Val; val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*c_start; min = 0.95*c_start; max = 1.05*c_start;\n") new_inp_file_stream.write((structure_lines[9].replace("\n", "").replace("al", "prm %s_dummy" % (refinement_flags[3]))) + " del 0.01 val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*al_start; min = 0.95*al_start; max = 1.05*al_start;\n") new_inp_file_stream.write((structure_lines[10].replace("\n", "").replace("be", "prm %s_dummy" % (refinement_flags[4]))) + " del 0.01 val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*be_start; min = 0.95*be_start; max = 1.05*be_start;\n") new_inp_file_stream.write((structure_lines[11].replace("\n", "").replace("ga", "prm %s_dummy" % (refinement_flags[5]))) + " del 0.01 val_on_continue = Rand((100-cell_perc)/100,(100+cell_perc)/100)*ga_start; min = 0.95*ga_start; max = 1.05*ga_start;\n") new_inp_file_stream.write("a %s = %s_dummy;" % (refinement_flags[0].replace("!",""), refinement_flags[0].replace("!","")) + " min = 0.95*a_start; max = 1.05*a_start;\n") new_inp_file_stream.write("b %s = %s_dummy;" % (refinement_flags[1].replace("!",""), refinement_flags[1].replace("!","")) + " min = 0.95*b_start; max = 1.05*b_start;\n") new_inp_file_stream.write("c %s = %s_dummy;" % (refinement_flags[2].replace("!",""), refinement_flags[2].replace("!","")) + " min = 0.95*c_start; max = 1.05*c_start;\n") new_inp_file_stream.write("al %s = %s_dummy;" % (refinement_flags[3].replace("!",""), refinement_flags[3].replace("!","")) + " min = 0.95*al_start; max = 1.05*al_start;\n") new_inp_file_stream.write("be %s = %s_dummy;" % (refinement_flags[4].replace("!",""), refinement_flags[4].replace("!","")) + " min = 0.95*be_start; max = 1.05*be_start;\n") new_inp_file_stream.write("ga %s = %s_dummy;" % (refinement_flags[5].replace("!",""), refinement_flags[5].replace("!","")) + " min = 0.95*ga_start; max = 1.05*ga_start;\n") new_inp_file_stream.write("volume vol 0\n") def ExcelResults(result_file): result_stream = open(result_file, 'r') result_lines = result_stream.readlines() result_stream.close() workbook = xlsxwriter.Workbook(r'paste_me.xlsx') worksheet = workbook.add_worksheet() for i, line in enumerate(result_lines): if i != 0: worksheet.write(i, 2, float(line.split()[1])) worksheet.write(i, 4, float(line.split()[2])) worksheet.write(i, 5, float(line.split()[3])) worksheet.write(i, 6, float(line.split()[4])) worksheet.write(i, 7, float(line.split()[5])) worksheet.write(i, 8, float(line.split()[6])) worksheet.write(i, 9, float(line.split()[7])) worksheet.write(i, 10, float(line.split()[8])) workbook.close() ##################################################################### # Execution Section ##################################################################### print "Initialising refinements of all subgroups in subgroup_numbers_to_run" ##################################################################### # Create Result Files and Folders (in Working Directory) ##################################################################### result_file_stream = open(result_file, 'w') result_file_stream.write("subgroup rwp a b c al be ga volume\n") result_file_stream.close() Folders_to_create = ["start_inp_files\\", "best_inp_files\\"] for folder in Folders_to_create: if not os.path.exists(os.path.dirname(working_directory + r"\\%s" % (folder))): os.makedirs(os.path.dirname(working_directory + r"\\%s" % (folder))) ##################################################################### # Subgroup Tests ##################################################################### for subgroup in subgroup_numbers_to_run: ##################################################################### # Make inp file ##################################################################### inp_maker(basefile, subgroup, thermal_prms) copyfile(working_directory + r'/' + basefile + "%04d" % (subgroup) + r""".inp""", working_directory + r'\start_inp_files/' + basefile + "%04d_start" % (subgroup) + r""".inp""") # save original (unrefined) input file (copy) for run in range(1,number_of_steps+1): # differences between runs are found in ifdef run_x statements in seed file ##################################################################### # Send to topas ##################################################################### print "Running subgroup number %d - Run %d..." % (subgroup, run) command = topas_directory + r"""\tc.exe """ + basefile + "%04d" % (subgroup) + r""".inp""" + r""" "#define GUI_LINES #define run_%d macro subgroup_number {prm !subgroup %d}""" % (run, subgroup) a = subprocess.call(command, creationflags=0x08000000) ##################################################################### # Save ouputs ##################################################################### copyfile(working_directory + r'/' + basefile + "%04d" % (subgroup) + r""".out""", working_directory + r'/' + basefile + "%04d" % (subgroup) + r""".inp""") # turn this step .out into next step .inp (copy and overwrite) os.rename(working_directory + r'/' + basefile + "%04d" % (subgroup) + r""".out""", working_directory + r'\best_inp_files/' + basefile + "%04d_best" % (subgroup) + r""".inp""") # save best (refined) input file after all steps (cut and paste) os.remove(working_directory + r'/' + basefile + "%04d" % (subgroup) + r""".inp""") ExcelResults(result_file) print "All Runs Complete!"