""" Utilities for multi scale modelling """ import numpy as np from matplotlib import pyplot as plt from aiida.orm import StructureData, Dict, Code, load_node, Group, load_group, ArrayData from aiida.engine import submit from aiida_kkr.workflows import kkr_scf_wc, kkr_dos_wc from aiida_kkr.tools import plot_kkr, kkrparams, find_parent_structure from masci_tools.util.constants import BOHR_A from aiida.engine import calcfunction from tqdm import tqdm def load_or_submit_jij_calc(scf, uuid=None, options=None, code=None): """Submission helper function to submit Jij calculations""" if uuid is not None: return load_node(uuid) from aiida_kkr.calculations import KkrCalculation kkr_calc_converged = scf.outputs.last_RemoteData.get_incoming(node_class=KkrCalculation).first().node Jij_params = kkrparams(**kkr_calc_converged.inputs.parameters.get_dict()) # reuse old settings # add JIJRAD (remember: in alat units) if len(find_parent_structure(kkr_calc_converged).sites)==1: Jij_params.set_multiple_values(JIJRAD=5.0, NATOMIMPD=800) nsym = 48 else: Jij_params.set_multiple_values(JIJRAD=2.5, NATOMIMPD=1600) nsym = 8 # add 'XCPL' runopt to list of runopts runopts = Jij_params.get_value('RUNOPT') runopts.append('XCPL ') Jij_params.set_value('RUNOPT', runopts) # increase k-point sampling Jij_params.set_multiple_values(BZDIVIDE=[100, 100, 100], KPOIBZ=int(100**3/nsym*1.5)) # now use updated parameters and continue from converged calculation builder = kkr_calc_converged.get_builder_restart() builder.metadata.label = 'Jij_calculation_noSOC_GGA' if options is not None: for k,v in options.items(): builder.metadata.options[k] = v if code is not None: builder.code = code builder.parent_folder = kkr_calc_converged.outputs.remote_folder builder.parameters = Dict(dict=Jij_params) jij_calc = submit(builder) print('submitted', jij_calc) return jij_calc def make_jij_builder(scf, options=None, code=None): """prepare process bulder for Jij calculation""" from aiida_kkr.calculations import KkrCalculation kkr_calc_converged = scf.outputs.last_RemoteData.get_incoming(node_class=KkrCalculation).first().node Jij_params = kkrparams(**kkr_calc_converged.inputs.parameters.get_dict()) # reuse old settings # add JIJRAD (remember: in alat units) Jij_params.set_multiple_values(JIJRAD=5.0, NATOMIMPD=800) # add 'XCPL' runopt to list of runopts runopts = Jij_params.get_value('RUNOPT') runopts.append('XCPL ') Jij_params.set_value('RUNOPT', runopts) # increase k-point sampling Jij_params.set_multiple_values(BZDIVIDE=[100, 100, 100], KPOIBZ=int(100**3/48*1.5)) # now use updated parameters and continue from converged calculation builder = kkr_calc_converged.get_builder_restart() if options is not None: for k,v in options.items(): builder.metadata.options[k] = v if code is not None: builder.code = code builder.parent_folder = kkr_calc_converged.outputs.remote_folder builder.parameters = Dict(dict=Jij_params) jij_calc = submit(builder) return builder def make_positions(spirit_calc, flatten_array=True): """construct positions of the spins in the spirit calculation cell""" # get structure information from spirit input structure struc = spirit_calc.inputs.structure cell = np.array(struc.cell) pos_cell = np.array([i.position for i in struc.sites]) n_basis_cells = spirit_calc.inputs.parameters.get_dict().get('n_basis_cells', [5, 5, 5]) # set positions and direciton of the spins positions = [] for iz in range(n_basis_cells[2]): for iy in range(n_basis_cells[1]): for ix in range(n_basis_cells[0]): for p in pos_cell: # set position positions.append(p + ix * cell[0] + iy * cell[1] + iz * cell[2]) # make flattened array positions = np.array(positions).reshape(-1, 3) if not flatten_array: positions = positions.reshape( n_basis_cells[0], n_basis_cells[1], n_basis_cells[2], len(pos_cell), 3 ) return positions def make_AFM_start(spirit_calc): p = make_positions(spirit_calc, flatten_array=False) m = np.zeros_like(p) m[:,:,:,0] = [0,0,1] m[:,:,:,1] = [0,0,-1] m[:,:,:,2] = [0,0,1] m[:,:,:,3] = [0,0,-1] init_AFM = ArrayData() init_AFM.set_array('initial_state', m.reshape(-1,3)) return init_AFM def prepare_spirit_calc_4(scf, jij_data, spirit_struct, label='', uuid=None, FM=False, **kwargs): """.""" if uuid is not None: return load_node(uuid) # parameters that go into the spirit configuration file input_para = { ### Hamiltonian # impurity spin moments in mu_Bohr # list of spin moments, needed to apply an external field 'mu_s': [scf.outputs.last_calc_out['magnetism_group']['spin_moment_per_atom'][0] for i in range(4)], # Spirit supercell 'boundary_conditions': [False, False, False], 'n_basis_cells': [40, 40, 40], # external magnetic field in T (point in +z direction) 'external_field_magnitude': 0.0, 'external_field_normal': [0.0, 0.0, 1.0], # single-ion anisotropy in meV 'anisotropy_magnitude': 0., 'anisotropy_normal': [0.0, 0.0, 1.0], ### LLG settings 'llg_n_iterations': 100000, 'llg_n_iterations_log': 1000, 'llg_temperature': 0, # temperature noise (K) 'llg_dt': 1.0E-3, # LLG time step (ps) 'llg_force_convergence': 1e-7, } for k,v in kwargs.items(): if k not in ['solver']: input_para[k] = v parameters = Dict(dict=input_para) # run modes of spirit that go into the run_spirit.py script ropt = { 'simulation_method': 'llg', 'solver': 'vp', } if 'solver' in kwargs: ropt['solver'] = kwargs['solver'] if FM: ropt['configuration'] = {'plus_z': True} run_opts = Dict(dict=ropt) #set up a spirit calculations from aiida.engine import submit from aiida.plugins import CalculationFactory builder = CalculationFactory('spirit').get_builder() builder.structure = spirit_struct builder.jij_data = jij_data builder.parameters = parameters builder.run_options = run_opts # code and computer options (queue name etc.) builder.code = Code.get_from_string('spirit@iffslurm') builder.metadata.options = { 'withmpi': False, # Spirit does not have MPI but by default used OpenMP if that is available 'resources': {'num_machines': 1, 'tot_num_mpiprocs': 1}, 'queue_name': 'oscar', # use oscar partition 'max_wallclock_seconds': 10*3600 # 10 h max runtime } builder.metadata.label = label return builder @calcfunction def expand_Jijs_big_cell(jij_data, struc): """ . """ j1 = jij_data.get_array('Jij_expanded') p1 = jij_data.get_array('positions_expanded') cell = np.array(struc.cell) a = cell[0,0]*2 cell4=np.array([[a,0,0], [0,a,0], [0,0,a]]) pos4 = np.array([ [0,0,0], cell[0], cell[1], cell[2] ]) struc4 = StructureData(cell = cell4) struc4.append_atom(position=pos4[0], symbols='Fe') struc4.append_atom(position=pos4[1], symbols='Fe') struc4.append_atom(position=pos4[2], symbols='Fe') struc4.append_atom(position=pos4[3], symbols='Fe') jij_big = [] for da in range(-5,5): for db in range(-5,5): for dc in range(-5,5): for i in range(4): for j in range(4): p = pos4[i] - pos4[j] + da*cell4[0] + db*cell4[1] + dc*cell4[2] jij_big.append([i, j, da, db, dc, 0, p[0], p[1], p[2]]) jij_big = np.array(jij_big) rmax = 10 jij_big = jij_big[np.linalg.norm(jij_big[:,-3:], axis=1)<=rmax] # limit max distance jij_big = jij_big[np.linalg.norm(jij_big[:,-3:], axis=1)>1e-5] # throw away zero distances jij_big = jij_big[np.linalg.norm(jij_big[:,-3:], axis=1).argsort()] # sort by distance # map by finding the distance for ii, p in enumerate(jij_big): d = np.linalg.norm(p[-3:]-p1, axis=1) j = np.where(d==d.min())[0] if len(j)==1: jij_big[ii,5] = j1[j[0],5] # # test plot # plt.plot(np.linalg.norm(jij_big[:,-3:], axis=1), jij_big[:,5], 'o--') jij_data_big = ArrayData() jij_data_big.set_array('Jij_expanded', jij_big[:,:-3]) jij_data_big.set_array('positions_expanded', jij_big[:,-3:]) return {'jij_data': jij_data_big, 'structure': struc4} def fft_central_spins(spirit_calc): """FFT of s_z along x=0, y=0 and along x=0, z=0 lines""" positions = make_positions(spirit_calc) m = spirit_calc.outputs.magnetization.get_array('final') alat = spirit_calc.inputs.structure.cell[0][0] # find central position diff = np.linalg.norm(positions - positions.mean(axis=0), axis=1) p0 = positions[diff==diff.min()][0] # pick x=0 mask = np.where(abs(positions[:,0]-p0[0])==abs(positions[:,0]-p0[0]).min()) p1, m1 = positions[mask], m[mask] # add second sub-lattice mask = np.where(abs(positions[:,0]-p0[0]-alat/2)==abs(positions[:,0]-p0[0]-alat/2).min()) p1 = np.array(list(p1) + list(positions[mask])) m1 = np.array(list(m1) + list(m[mask])) # along y=0 mask = np.where(abs(p1[:,1]-p0[1])==abs(p1[:,1]-p0[1]).min()) p2, m2 = p1[mask], m1[mask] isort = p2[:,2].argsort() p2, m2 = p2[isort], m2[isort] # along z=0 mask = np.where(abs(p1[:,2]-p0[2])==abs(p1[:,2]-p0[2]).min()) p3, m3 = p1[mask], m1[mask] isort = p3[:,1].argsort() p3, m3 = p3[isort], m3[isort] # fft only of z-component ffz = abs(np.fft.fftshift(np.fft.fft(m2[:,2]))) # positions along z=0, x=0, i.e. in z-direction ffy = abs(np.fft.fftshift(np.fft.fft(m3[:,2]))) # positions along y=0, x=0, i.e. in y-direction # pick y=0 mask = np.where(abs(positions[:,1]-p0[1])==abs(positions[:,1]-p0[1]).min()) p4, m4 = positions[mask], m[mask] # add second sub-lattice mask = np.where(abs(positions[:,1]-p0[1]-alat/2)==abs(positions[:,1]-p0[1]-alat/2).min()) p4 = np.array(list(p4) + list(positions[mask])) m4 = np.array(list(m4) + list(m[mask])) # along z=0 mask = np.where(abs(p4[:,2]-p0[2])==abs(p4[:,2]-p0[2]).min()) p5, m5 = p4[mask], m4[mask] isort = p5[:,0].argsort() p5, m5 = p5[isort], m5[isort] # fft only of z-component ffx = abs(np.fft.fftshift(np.fft.fft(m5[:,2]))) # positions along y=0, z=0, i.e. in x-direction # fft frequency x = np.fft.fftshift(np.fft.fftfreq(len(ffy))) * 2 return x, ffx, ffy, ffz, p5[:,0], m5, p3[:,1], m3, p2[:,2], m2