-
Notifications
You must be signed in to change notification settings - Fork 1
An example of a model script
canrong qiu edited this page May 29, 2021
·
2 revisions
Here is one example of a completed model script, which is generated from the script generating wizard mentioned above. The script is tagged using specially designed tags for code block modification as well as a better clarity to read.
import os
import models.structure_tools.sxrd_dafy as model
from models.utils import UserVars
import numpy as np
import batchfile.locate_path as batch_path
import dump_files.locate_path as output_path
import models.domain_creator as domain_creator
from models.structure_tools import tool_box
from models.structure_tools import sorbate_tool_beta as sorbate_tool
import models.setup_domain_hematite_rcut as setup_domain_hematite_rcut
from accessory_functions.data_formating.data_formating import format_hkl
from UtilityFunctions import config_file_parser_bv, update_O_NUMBER, setup_raxr_pars_new
from FitEnginePool import bond_valence_constraint
model_type = 'ctr'
#/raxs/begin#
RAXS_EL = 'As'
RAXS_FIT_MODE = 'MI'
NUMBER_SPECTRA = 9
E0 = 11867
F1F2_FILE = '/Users/canrong/apps/DaFy/examples/data_examples/As_K_edge_March28_2018.f1f2'
#/raxs/end#
if NUMBER_SPECTRA!=0:
rgh_raxs,F1F2 = setup_raxr_pars_new(NUMBER_SPECTRA, F1F2_FILE)
else:
rgh_raxs, F1F2 = None, []
#--global settings--#
#/globalsetting/begin#
#/path/begin#
batch_path_head=batch_path.module_path_locator()
output_file_path=output_path.module_path_locator()
#/path/end#
#/bv/begin#
USE_BV = True
#/bv/end#
#/wavelength/begin#
wal = 1.4922
#/wavelength/end#
#/slabnumber/begin#
num_surface_slabs = 1
num_sorbate_slabs = 2
#/slabnumber/end#
#--set unitcell--#
#/unitcell/begin#
lat_pars = [5.038, 5.434, 7.3707, 90, 90, 90]
unitcell = model.UnitCell(*lat_pars)
#/unitcell/end#
#/expconstant/begin#
re = 2.818e-5#electron radius
kvect=2*np.pi/wal#k vector
Egam = 6.626*(10**-34)*3*(10**8)/wal*10**10/1.602*10**19#energy in ev
LAM=1.5233e-22*Egam**6 - 1.2061e-17*Egam**5 + 2.5484e-13*Egam**4 + 1.6593e-10*Egam**3 + 1.9332e-06*Egam**2 + 1.1043e-02*Egam
exp_const = 4*kvect/LAM
auc=unitcell.a*unitcell.b*np.sin(unitcell.gamma)
#/expconstant/end#
#globalsetting/end#
#--set instrument--#
#/instrument/begin#
inst = model.Instrument(wavel = wal, alpha = 2.0)
#/instrument/end#
#--set bulk slab--#
#/bulk/begin#
bulk = model.Slab(T_factor = 'B')
tool_box.add_atom_in_slab(bulk,'/Users/canrong/apps/DaFy/util/batchfile/hematite_rcut/bulk.str')
#/bulk/end#
#--set surface slabs--#
#/surfaceslab/begin#
surface_1 = model.Slab(T_factor = 'B')
tool_box.add_atom_in_slab(surface_1, '/Users/canrong/apps/DaFy/util/batchfile/hematite_rcut/half_layer2.str')
#/surfaceslab/end#
#--set sorbate properties--#
#/sorbateproperties/begin#
sorbate_instance_1 = sorbate_tool.Tetrahedra.build_instance(substrate_domain = surface_1,anchored_ids={'attach_atm_ids': ['O1_1_0', 'O1_2_0'], 'offset': [None, None], 'anchor_ref': None, 'anchor_offset': None},binding_mode = {'mode': 'BD'},structure_pars_dict={'top_angle_offset': 0.0, 'angle_offset': [0, 0], 'phi': 0.0, 'edge_offset': [0, 0]},anchor_id = 'As1',ids = ['As1', 'O1', 'O2'], els = ['As', 'O', 'O'], lat_pars = [5.038, 5.434, 7.3707, 90, 90, 90], T = unitcell.lattice.RealTM, T_INV = unitcell.lattice.RealTMInv)
domain_sorbate_1 = sorbate_instance_1.domain
rgh_sorbate_1 = sorbate_instance_1.rgh
atm_gp_sorbate_1 = sorbate_instance_1.make_atm_group(instance_name = 'atm_gp_sorbate_1')
sorbate_instance_2 = sorbate_tool.Tetrahedra.build_instance(substrate_domain = surface_1,anchored_ids={'attach_atm_ids': ['O1_1_0', 'O1_2_0', 'O1_3_0'], 'offset': [None, None, None], 'anchor_ref': None, 'anchor_offset': None},binding_mode = {'mode': 'TD', 'mirror': True},structure_pars_dict={'center_offset': 0.0, 'edge_offset': 0},anchor_id = 'As1',ids = ['As1', 'O1'], els = ['As', 'O'], lat_pars = [5.038, 5.434, 7.3707, 90, 90, 90], T = unitcell.lattice.RealTM, T_INV = unitcell.lattice.RealTMInv)
domain_sorbate_2 = sorbate_instance_2.domain
rgh_sorbate_2 = sorbate_instance_2.rgh
atm_gp_sorbate_2 = sorbate_instance_2.make_atm_group(instance_name = 'atm_gp_sorbate_2')
#/sorbateproperties/end#
#/sorbatestructure/begin#
#
#/sorbatestructure/end#
#--set rgh--#
#/rgh/begin#
#/rgh/global/begin#
rgh = UserVars()
rgh.new_var('beta', 0.0)#roughness factor
rgh.new_var('mu',1)#liquid film thickness
scales=['scale_nonspecular_rods','scale_specular_rod']
for scale in scales:
rgh.new_var(scale,0.0005)
#/rgh/global/end#
#/rgh/layer_water/begin#
rgh_lw = UserVars()
pars=['u0','ubar','d_w','first_layer_height', 'density_w']
rgh_lw.new_var('u0', 2)
rgh_lw.new_var('ubar', 2)
rgh_lw.new_var('first_layer_height', 2)
rgh_lw.new_var('d_w', 2)
rgh_lw.new_var('density_w', 0.033)
#/rgh/layer_water/end#
#/rgh/domain_weight/begin#
rgh_wt = UserVars()
for i in range(num_surface_slabs):
rgh_wt.new_var('wt_domain{}'.format(i+1),1)
#/rgh/domain_weight/end#
#/rgh/end#
#/atmgroup/begin#
atm_gp_1 = model.AtomGroup(instance_name = 'atm_gp_1')
atm_gp_1.add_atom(surface_1,'O1_2_0',matrix=[1, 0, 0, 0, 1, 0, 0, 0, 1])
atm_gp_1.add_atom(surface_1,'O1_1_0',matrix=[1, 0, 0, 0, 1, 0, 0, 0, 1])
atm_gp_2 = model.AtomGroup(instance_name = 'atm_gp_2')
atm_gp_2.add_atom(surface_1,'O1_4_0',matrix=[1, 0, 0, 0, 1, 0, 0, 0, 1])
atm_gp_2.add_atom(surface_1,'O1_3_0',matrix=[1, 0, 0, 0, 1, 0, 0, 0, 1])#/atmgroup/end#
#/sorbatesym/begin#
sorbate_syms_1=[model.SymTrans([[1,0],[0,1]],t=[0,0]),model.SymTrans([[-1, 0], [0, 1]],t=[0.5, 0.5])]
sorbate_syms_2=[model.SymTrans([[1,0],[0,1]],t=[0,0]),model.SymTrans([[-1, 0], [0, 1]],t=[0.5, 0.5])]
#/sorbatesym/end#
#/sample/begin#
surface_parms = {'delta1': 0.0, 'delta2': 0.1391}
domains = {}
for i in range(num_surface_slabs):
domains['domain{}'.format(i+1)] = {}
domains['domain{}'.format(i+1)]['slab'] = globals()['surface_{}'.format(i+1)]
domains['domain{}'.format(i+1)]['sorbate'] = [globals()['domain_sorbate_{}'.format(j+1)] for j in range(num_sorbate_slabs)]
domains['domain{}'.format(i+1)]['wt'] = getattr(globals()['rgh_wt'],'wt_domain{}'.format(i+1))
domains['domain{}'.format(i+1)]['sorbate_sym'] = globals()['sorbate_syms_{}'.format(i+1)]
domains['domain{}'.format(i+1)]['layered_water'] = rgh_lw
sample = model.Sample(inst, bulk, domains, unitcell, surface_parms = surface_parms)
setattr(sample, 'rgh_raxs', rgh_raxs)
setattr(sample, 'E0', E0)
setattr(sample, 'f1f2', F1F2)
setattr(sample, 'res_el', RAXS_EL)
setattr(sample, 'mode', RAXS_FIT_MODE)
#/sample/end#
#setup bond valence attributes
locals().update(config_file_parser_bv(os.path.join(batch_path_head,'bv_data_base','config_bond_valence_db.ini')))
for i in range(num_surface_slabs):
vars()['bv_constraint_domain{}'.format(i+1)] = bond_valence_constraint.factory_function_bv(r0_container = R0_BV,\
domain_list= [domains['domain{}'.format(i+1)]['slab']] + [globals()['domain_sorbate_{}'.format(j+1)] for j in range(num_sorbate_slabs)],\
TM = unitcell.lattice.RealTM)
#a long list storing the info whether each dataset is used (flatten to the full length of each dataset)
data_use_array = np.array(sum([[each_set.use]*len(each_set.x) for each_set in data],[]))
def Sim(data,VARS=vars(),kwargs = {}):
F =[]
fom_scaler=[]
beta=rgh.beta
#/update_sorbate/begin#
sorbate_instance_2.update_structure()
sorbate_instance_1.update_structure()
#/update_sorbate/end#
#normalize the domain weight to make total = 1
wt_list = [getattr(rgh_wt, 'wt_domain{}'.format(i+1)) for i in range(num_surface_slabs)]
total_wt = sum(wt_list)
for i in range(num_surface_slabs):
sample.domain['domain{}'.format(i+1)]['wt']=wt_list[i]/total_wt
#this is used to generate dummy data from GUI dialog
if 'ctr' in kwargs:
h_ctr, k_ctr, l_ctr = kwargs['ctr'][:,0], kwargs['ctr'][:,1], kwargs['ctr'][:,2]
f_ctr = sample.calc_f_all(h_ctr, k_ctr, l_ctr)
F_ctr = abs(f_ctr*f_ctr)
Ferr_ctr = F_ctr*0.01
y_ctr = F_ctr*0
LB_ctr = [2]*len(F_ctr)
dL_ctr = [2]*len(F_ctr)
dummy_data_ctr = np.array([l_ctr,h_ctr,k_ctr,y_ctr,F_ctr,Ferr_ctr,dL_ctr,LB_ctr]).T
if 'raxs' in kwargs:
h_rs, k_rs, l_rs, E_rs = kwargs['raxs'][:,0], kwargs['raxs'][:,1], kwargs['raxs'][:,2],kwargs['raxs'][:,3]
f_rs = sample.calc_f_all_RAXS(h_rs, k_rs, l_rs, E_rs)
F_rs = abs(f_rs*f_rs)
LB_rs = [2]*len(F_rs)
dL_rs = [2]*len(F_rs)
dummy_data_rs = np.array([E_rs,h_rs,k_rs,l_rs,F_rs, F_rs*0.001, dL_rs, LB_rs]).T
kwargs['func'](np.vstack((dummy_data_ctr,dummy_data_rs)))
else:
kwargs['func'](dummy_data_ctr)
#faster solution(a factor of two faster than using loop)
#ctr datasets
condition_ctr = data.ctr_data_all[:,-1]<100
condition_used_ctr = data.binary_comparison_and(data_use_array, condition_ctr)
if True in list(condition_used_ctr):
h_, k_, x_,LB_,dL_ = data.ctr_data_all[condition_used_ctr][:,0], data.ctr_data_all[condition_used_ctr][:,1], data.ctr_data_all[condition_used_ctr][:,2],data.ctr_data_all[condition_used_ctr][:,4],data.ctr_data_all[condition_used_ctr][:,5]
rough_ = (1-beta)/((1-beta)**2 + 4*beta*np.sin(np.pi*(x_-LB_)/dL_)**2)**0.5
f_ = rough_*sample.calc_f_all(h_, k_, x_)
F_ = abs(f_*f_)#either f_*f_ for intensity or f_ only structure factor
F_all = data.ctr_data_all[condition_ctr][:,0]*0
sub_sets, data_info = data.split_used_dataset(F_, data_type = 'CTR')
F_ = data.insert_datasets(full_set = F_all, sub_sets = sub_sets, data_info = data_info, data_type = 'CTR')
#you need to edit the list of extra scaling factor accordingly
scaling_factors = [[rgh.scale_nonspecular_rods, rgh.scale_specular_rod][int(each=='specular_rod')] for each in data.scaling_tag]
F_ctr = data.split_fullset(F_,scaling_factors)
else:
F_ctr = data.split_fullset(data.ctr_data_all[condition_ctr][:,0]*0,scale_factors=1)
#raxs datasets
condition_raxs = data.ctr_data_all[:,-1]>=100
condition_used_raxs = data.binary_comparison_and(data_use_array, condition_raxs)
if True in list(condition_used_raxs):
h_, k_, E_, l_, LB_,dL_, tag= data.ctr_data_all[condition_used_raxs][:,0], data.ctr_data_all[condition_used_raxs][:,1], data.ctr_data_all[condition_used_raxs][:,2],data.ctr_data_all[condition_used_raxs][:,3], data.ctr_data_all[condition_used_raxs][:,4],data.ctr_data_all[condition_used_raxs][:,5],data.ctr_data_all[condition_used_raxs][:,-1]-100
rough_ = (1-beta)/((1-beta)**2 + 4*beta*np.sin(np.pi*(l_-LB_)/dL_)**2)**0.5
f_ = rough_*sample.calc_f_all_RAXS(h_, k_, l_, E_, tag)
F_ = abs(f_*f_)#either f_*f_ for intensity or f_ only structure factor
F_all = data.ctr_data_all[condition_raxs][:,0]*0
sub_sets, data_info = data.split_used_dataset(F_, data_type = 'RAXS')
F_ = data.insert_datasets(full_set = F_all, sub_sets = sub_sets, data_info = data_info, data_type = 'RAXS')
#you need to edit the list of extra scaling factor accordingly
scaling_factors = [[rgh.scale_nonspecular_rods, rgh.scale_specular_rod][int(each=='specular_rod')] for each in data.scaling_tag_raxs]
F_raxs = data.split_fullset(F_,scale_factors=scaling_factors, data_type = 'RAXS')
else:
F_raxs=data.split_fullset(data.ctr_data_all[condition_raxs][:,0]*0,scale_factors=1, data_type = 'RAXS')
#Now merge both datasets together
F = data.merge_datasets(ctr_datasets = F_ctr, raxs_datasets = F_raxs)
fom_scaler = [1]*len(F)
#calculate bv panelty factor
bv = 0
if USE_BV:
for i in range(num_surface_slabs):
bv += VARS['bv_constraint_domain{}'.format(i+1)].cal_distance()
return F,1+bv,fom_scaler