diff --git a/nafuma/xrd/refinement.py b/nafuma/xrd/refinement.py index 4fe4f78..1b33bc3 100644 --- a/nafuma/xrd/refinement.py +++ b/nafuma/xrd/refinement.py @@ -16,19 +16,25 @@ import nafuma.auxillary as aux def make_initial_inp(data: dict, options={}): - required_options = ['filename', 'overwrite', 'include', 'save_results', 'save_dir', 'topas_options', 'background', 'capillary', 'th2_offset', 'fit_peak_width', 'start', 'finish', 'exclude'] + required_options = ['filename', 'overwrite', 'include', 'save_results', 'save_dir', 'instrument', 'topas_options', 'background', 'capillary', 'th2_offset', 'fit_peak_width', 'simple_axial_model', 'TCHZ_Peak_Type', 'start', 'finish', 'exclude', 'magnetic_atoms', 'radiation', 'magnetic_space_group', 'interval'] default_options = { 'filename': 'start.inp', 'overwrite': False, 'save_results': False, 'save_dir': 'results/', + 'instrument': None, + 'radiation': 'synchrotron', + 'magnetic_space_group': None, + 'magnetic_atoms': [], 'include': [], # Any files that should be included with #include 'topas_options': {}, 'background': 7, 'capillary': False, 'th2_offset': False, 'fit_peak_width': False, + 'simple_axial_model': False, + 'TCHZ_Peak_Type': False, 'interval': [None, None], # Start and finish values that TOPAS should refine on. Overrides 'start' and 'finish' 'start': None, # Start value only. Overridden by 'interval' if this is set. 'finish': None, # Finish value only. Overridden by 'interval' if this is set. @@ -72,6 +78,16 @@ def make_initial_inp(data: dict, options={}): write_str(fout=fout, data=data, options=options, index=i) + if 'peak' in data.keys(): + + fout.write('\n\'------------------------------------------------------------------------------------------------------------------------------------------\n') + fout.write('\'PEAK PHASES') + fout.write('\n\'------------------------------------------------------------------------------------------------------------------------------------------\n') + + write_peaks(fout=fout, data=data, options=options, index=i) + + + with open(options['filename'], 'r') as fin: lines = [] @@ -124,8 +140,21 @@ def write_xdd(fout, data, options): # Write wavelength and LP-factor fout.write('\n\n') + + if options['instrument'] == 'RECX2': + options['radiation'] = 'MoKa' + for line in snippets['RECX2']: + fout.write('\t'+line+'\n') + + if options['radiation'] == 'synchrotron': + fout.write('\t'+snippets['synchrotron'].format(data['wavelength'])+'\n') + if options['radiation'] == 'neutron': + fout.write('\t'+snippets['neutron'][0].format(data['wavelength']) +'\n') + fout.write('\t'+snippets['neutron'][1] + '\n') + if options['radiation'] == 'MoKa': + for line in snippets['MoKa']: + fout.write('\t'+line+'\n') - fout.write('\t'+snippets['wavelength'].format(data['wavelength'])+'\n') fout.write('\t'+snippets['lp_factor'].format(topas_options['lp_factor'])+'\n') if options['th2_offset']: @@ -190,21 +219,21 @@ def write_params(fout, data, options, index=0): # If start_values is defined: fout.write(f'#ifdef start_values_{label}\n') data['define_statements'].append(f'start_values_{label}') - lpa = f'local !lpa_{label} {a} ;: {a}' - lpb = f'local !lpb_{label} {b} ;: {b}' - lpc = f'local !lpc_{label} {c} ;: {c}' + lpa = f'local !lpa_{label} = {a} ;: {a}' + lpb = f'local !lpb_{label} = {b} ;: {b}' + lpc = f'local !lpc_{label} = {c} ;: {c}' fout.write('{: <55} {: <55} {: <55}\n'.format(lpa, lpb, lpc)) - lpal = f'local !lpal_{label} {alpha} ;: {alpha}' - lpbe = f'local !lpbe_{label} {beta} ;: {beta}' - lpga = f'local !lpga_{label} {gamma} ;: {gamma}' + lpal = f'local !lpal_{label} = {alpha} ;: {alpha}' + lpbe = f'local !lpbe_{label} = {beta} ;: {beta}' + lpga = f'local !lpga_{label} = {gamma} ;: {gamma}' fout.write('{: <55} {: <55} {: <55}\n'.format(lpal, lpbe, lpga)) # Otherwise fout.write('\n') fout.write('#else\n') - lpa = f'local !lpa_{label} {a}' - lpb = f'local !lpb_{label} {b}' - lpc = f'local !lpc_{label} {c}' + lpa = f'local !lpa_{label} = {a} ;: {a}' + lpb = f'local !lpb_{label} = {b} ;: {b}' + lpc = f'local !lpc_{label} = {c} ;: {c}' fout.write('{: <55} {: <55} {: <55}\n'.format(lpa, lpb, lpc)) lpal = f'local !lpal_{label} {alpha}' lpbe = f'local !lpbe_{label} {beta}' @@ -224,6 +253,7 @@ def write_params(fout, data, options, index=0): # WRITE SITE PARAMETERS + mag = [] for site in sites: params = [] @@ -232,13 +262,30 @@ def write_params(fout, data, options, index=0): value = atoms["atoms"][site][attr].split("(")[0] value = value if value != '.' else 0. + + params.append('{: <20} {: <20}'.format(f'local !{attrs[attr]}_{site}_{label}', f' = {value} ;: {value}')) - #fout.write(f'local {attrs[attr]}_{site}_{label}\t\t =\t {value} \t ;= \t\t\t\t\t') + + if options['magnetic_space_group']: + for mag_atom in options['magnetic_atoms']: + if mag_atom in site: + mag.append(site) fout.write('{: <55} {: <55} {: <55} {: <55} {: <55}\n'.format(*params)) fout.write('\n') + if options['magnetic_space_group']: + mag = list(dict.fromkeys(mag)) # remove duplicates + for m in mag: + fout.write('{: <55} {: <55} {: <55}\n'.format( + f'local !ml_x_{m}_{label}_XXXX = 0 ;: 0', + f'local !ml_y_{m}_{label}_XXXX = 0 ;: 0', + f'local !ml_z_{m}_{label}_XXXX = 0 ;: 0' + )) + + fout.write('\n') + fout.write('{: <55} {: <55} {: <55} {: <55}\n'.format( f'local !csgc_{label}_XXXX = 200 ;: 200', f'local !cslc_{label}_XXXX = 200 ;: 200', @@ -277,7 +324,10 @@ def write_str(fout, data, options, index=0): fout.write('\tstr\n') fout.write(f'\t\tphase_name "{atoms["_chemical_name_common"]} ({atoms["_space_group_name_H-M_alt"]})"\n'.replace('\'', '')) - fout.write(f'\t\tspace_group {atoms["_space_group_IT_number"]}\n') + if options['magnetic_space_group']: + fout.write(f'\t\tmag_space_group {options["magnetic_space_group"]}\n') + else: + fout.write(f'\t\tspace_group {atoms["_space_group_IT_number"]}\n') fout.write(f'\t\ta = lpa_{label} ;\n') fout.write(f'\t\tb = lpb_{label} ;\n') fout.write(f'\t\tc = lpc_{label} ;\n') @@ -287,8 +337,11 @@ def write_str(fout, data, options, index=0): fout.write('\n') fout.write(f'\t\tcell_volume\t vol_{label}_XXXX {atoms["_cell_volume"]}\n') fout.write(f'\t\tcell_mass\t mass_{label}_XXXX 1\n') - fout.write(f'\t\tweight_percent\t wp_{label}_XXXX 100\n\n') - fout.write('\n') + fout.write(f'\t\tweight_percent\t wp_{label}_XXXX 100\n\n\n') + if options['simple_axial_model']: + fout.write(f'\t\t{snippets["Simple_Axial_Model"]}\n') + if options['TCHZ_Peak_Type']: + fout.write(f'\t\t{snippets["TCHZ_Peak_Type"]}\n\n\n') fout.write('\t\tscale @ 1.0\n') fout.write('\t\tr_bragg 1.0\n\n') if options['fit_peak_width']: @@ -342,6 +395,11 @@ def write_str(fout, data, options, index=0): fout.write('\t\t{: <9} {: <30}; {: <30}; {: <30}; {: <30}; {: <30};'.format(site, x, y, z, occ, beq)) fout.write('\n') + if options['magnetic_space_group']: + for mag_atom in options['magnetic_atoms']: + if mag_atom in atom: + fout.write('\t\t'+snippets['magnetic_moment_str'].format(atom, label, atom, label, atom, label)+'\n') + if options['save_results']: @@ -364,11 +422,13 @@ def write_output(fout, data, options, index=0): fout.write('#ifdef output\n') - fout.write(f'\t\tOut_Riet({options["save_dir"]}/{filename}_{label}_XXXX_riet.xy)\n') - fout.write(f'\t\tOut_CIF_STR({options["save_dir"]}/{filename}_{label}_XXXX.cif)\n') - fout.write(f'\t\tOut_CIF_ADPs({options["save_dir"]}/{filename}_{label}_XXXX.cif)\n') - fout.write(f'\t\tOut_CIF_Bonds_Angles({options["save_dir"]}/{filename}_{label}_XXXX.cif)\n') - fout.write(f'\t\tCreate_hklm_d_Th2_Ip_file({options["save_dir"]}/{filename}_{label}_XXXX_hkl.dat)\n') + fout.write(f'\t\tOut_Riet({options["save_dir"]}/{filename}_{label}_riet.xy)\n') + fout.write(f'\t\tOut_CIF_STR({options["save_dir"]}/{filename}_{label}.cif)\n') + fout.write(f'\t\tOut_CIF_ADPs({options["save_dir"]}/{filename}_{label}.cif)\n') + fout.write(f'\t\tOut_CIF_Bonds_Angles({options["save_dir"]}/{filename}_{label}.cif)\n') + fout.write(f'\t\tCreate_hklm_d_Th2_Ip_file({options["save_dir"]}/{filename}_{label}_hkl.dat)\n') + if options['magnetic_space_group']: + fout.write(f'\t\tOut_CIF_mag({options["save_dir"]}/{filename}_{label}_magnetic.cif)\n') fout.write('\n') @@ -446,7 +506,26 @@ def write_output(fout, data, options, index=0): fout.write('\n\n') +def write_peaks(fout, data, options, index=0): + import nafuma.xrd as xrd + basedir = os.path.dirname(xrd.__file__) + + with open(os.path.join(basedir, 'snippets.json'), 'r') as snip: + snippets = json.load(snip) + + + if not isinstance(data['peak'], list): + data['peak'] = [data['peak']] + + for peak in data['peak']: + fout.write('\t'+snippets['peak'][0]+'\n') + fout.write('\t\t'+snippets['peak'][1].format(peak)+'\n') + fout.write('\t\t'+snippets['peak'][2]+'\n') + fout.write('\t\t'+snippets['peak'][3]+'\n') + fout.write('\t\t'+snippets['peak'][4]+'\n') + + fout.write('\n\n') def read_cif(path): @@ -600,7 +679,7 @@ def make_big_inp(data: dict, options={}): aux.backup_file(filename=options['output'], backup_dir=options['backup_dir']) - runlist = os.path.join(os.path.dirname(options['output']), 'runlist.txt') + runlist = os.path.abspath(os.path.join(os.path.dirname(options['output']), f'runlist_{os.path.basename(options["output"]).split(".")[0]}.txt')) options['include'].append(runlist) with open(options['template'], 'r') as template: @@ -615,6 +694,8 @@ def make_big_inp(data: dict, options={}): for i, path in enumerate(data['path']): + path = os.path.abspath(path) + s = make_inp_entry(template=template, xdd=path, num=i, options=options) output.write(s) output.write('\n\n') @@ -745,6 +826,7 @@ def make_inp_entry(template: str, xdd: str, num: int, options: dict) -> str: # Replace diffractogram-path s = template.replace(temp_xdd, xdd).replace('XXXX', num_str) + s = s.replace(os.path.basename(temp_xdd).split('.')[0], os.path.basename(xdd).split('.')[0]) return s @@ -867,3 +949,73 @@ def read_results(path): return results + + +def fix_magnetic_cif(path): + with open(path, 'r') as f: + lines = f.readlines() + + for i, line in enumerate(lines): + if '_magnetic_space_group_symop_operation_timereversal' in line: # last line before symmetry operations + j = 1 + line = lines[i + j] + while len(line) > 1: + lines[i+j] = f'{j} ' + line + j += 1 + line = lines[i+j] + + + + with open(path, 'w') as f: + for line in lines: + f.write(line) + + + +def get_magnetic_moment(path): + with open(path, 'r') as f: + lines = f.readlines() + + for i, line in enumerate(lines): + if '_magnetic_atom_site_moment_crystalaxis_mz' in line: # last line before magnetic moments + j = 1 + line = lines[i+j] + + magnetic_moments = {} + while line: + magnetic_moments[line.split()[0]] = (float(line.split()[1]), float(line.split()[2]), float(line.split()[3])) + + j += 1 + + if i+j < len(lines): + line = lines[i+j] + else: + break + + + + return magnetic_moments + + +def get_magnetic_moment_magnitude(magmom): + magnitudes = {} + for site in magmom.keys(): + magnitudes[site] = np.sqrt(magmom[site][0]**2+magmom[site][1]**2+magmom[site][2]**2) + + return magnitudes + + +def strip_labels_from_inp(path, save_path=None): + + if not save_path: + save_path = path + + with open(path, 'r') as fin: + lines = fin.readlines() + + for i, line in enumerate(lines): + lines[i] = line.replace('_XXXX', '') + + with open(save_path, 'w') as fout: + for line in lines: + fout.write(line) diff --git a/nafuma/xrd/snippets.json b/nafuma/xrd/snippets.json index 00266bd..60f169d 100644 --- a/nafuma/xrd/snippets.json +++ b/nafuma/xrd/snippets.json @@ -9,13 +9,36 @@ ], "gauss_fwhm": "gauss_fwhm = Sqrt({} Cos(2 * Th)^4 + {} Cos(2 * Th)^2 + {});", "lp_factor": "LP_Factor({}) 'change the LP correction or lh value if required", - "wavelength": "lam ymin_on_ymax 0.0001 la 1.0 lo {} lh 0.1", - "zero_error": "Zero_Error(zero, 0)", + "synchrotron": "lam ymin_on_ymax 0.0001 la 1.0 lo {} lh 0.1", + "neutron": [ + "lam ymin_on_ymax 0.0001 la 1.0 lo {} lh 0.5", + "neutron_data" + ], + "MoKa":[ + "lam ymin_on_ymax 0.0001", + "la 0.6533 lo 0.7093 lh 0.2695", + "la 0.3467 lo 0.713574 lh 0.2795" + ], + "RECX2": [ + "Rp 280", + "Rs 280" + ], + "zero_error": "Zero_Error(!zero, 0)", "th2_offset": [ - "prm zero \t\t\t0 \t\t\t\tmin = Max(Val - 20 Yobs_dx_at(X1), -100 Yobs_dx_at(X1)); max = Min(Val + 20 Yobs_dx_at(X2), 100 Yobs_dx_at(X2)); del = .01 Yobs_dx_at(X1); val_on_continue 0", - "prm cos_shift \t\t0 \t\t\t\tmin = Val-.8; max = Val+.8; del 0.001", - "prm sin_shift \t\t0 \t\t\t\tmin = Val-.8; max = Val+.8; del 0.001", + "prm !zero\t\t\t= 0 ;: 0 \t\t\t\tmin = Max(Val - 20 Yobs_dx_at(X1), -100 Yobs_dx_at(X1)); max = Min(Val + 20 Yobs_dx_at(X2), 100 Yobs_dx_at(X2)); del = .01 Yobs_dx_at(X1); val_on_continue 0", + "prm !cos_shift\t\t= 0 ;: 0 \t\t\t\tmin = Val-.8; max = Val+.8; del 0.001", + "prm !sin_shift\t\t= 0 ;: 0 \t\t\t\tmin = Val-.8; max = Val+.8; del 0.001", "th2_offset = (zero) + (cos_shift) Cos(Th) + (sin_shift) Sin(Th) ;" ], - "fit_peak_width": "DC1( ad, 0, bd, 0, cd, 0)" + "fit_peak_width": "DC1( ad, 0, bd, 0, cd, 0)", + "TCHZ_Peak_Type": "TCHZ_Peak_Type(pku_1, 0, pkv_1, 0,pkw_1, 0, !pkx_1, 0.0000,pky_1, 0,!pkz_1, 0.0000)", + "Simple_Axial_Model": "Simple_Axial_Model( axial_1, 0)", + "magnetic_moment_str": "mlx = ml_x_{}_{}_XXXX ; \t mly = ml_y_{}_{}_XXXX ; \t mlz = ml_z_{}_{}_XXXX ; \t MM_CrystalAxis_Display( 0, 0, 0)", + "peak": [ + "xo_Is", + "xo @ {}", + "peak_type fp", + "LVol_FWHM_CS_G_L( 1, 0, 0.89, 0,,,@, 2)", + "I @ 35.35632`" + ] } \ No newline at end of file