Example for AiiDA-Spirit paper to demonstrate the multi scale modelling, i.e. interoperability with AiiDA-KKR
Pseudo code for workflow:
# 1. DFT calculations
for alat in alat_range:
# first run DFT self-consistency
scf_all = run_KKR_scf(alat)
# then compute DOS and Jij's
dos_all = run_DOS(scf_all)
Jij_all = run_Jij(scf_all)
# 2. then Spirit calculations (LLG)
for Jij in Jij_all:
# starting from random spin configuration
LLG_random_all = run_LLG(Jij)
# starting from ferromagnetic
LLG_FM_all = run_LLG_FM(Jij)
from aiida import load_profile
load_profile()
from util1 import submit_batch_iffslurm_all_free
from util3 import *
struc = StructureData(cell=np.array([[0.5,0.5,0], [0.5,0,0.5], [0,0.5,0.5]])*6.82*BOHR_A)
struc.append_atom(position=[0,0,0], symbols='Fe')
struc.label = 'fcc_Fe_alat=Cu'
struc.description = 'Fe in fcc lattice with lattice constant of Cu (https://www.flapw.de/MaX-5.1/future/M2/)'
# plot_kkr(struc, silent=True)
builder = kkr_scf_wc.get_builder()
# scf settings
scf_settings = kkr_scf_wc.get_wf_defaults(silent=True)[0]
scf_settings['mag_init'] = True
scf_settings['check_dos'] = False
scf_settings['nsteps'] = 200
scf_settings.convergence_setting_coarse = scf_settings.convergence_setting_fine # increases energy mesh to use ranks more efficiently
# code that uses AMD nodes and has the BdG functionality
kkr_code = Code.get_from_string('kkrhost_3.5_AMD@iffslurm')
voronoi_code = Code.get_from_string('voronoi_3.5_AMD@iffslurm')
builder.kkr = kkr_code
builder.voronoi = voronoi_code
options_AMD2 = Dict(dict={
'withmpi': True,
'resources': {'num_machines': 1, 'tot_num_mpiprocs': 24},
'queue_name': 'th1-2020-64', # select AMD nodes that match the AMD codes (see above)
'custom_scheduler_command': 'export OMP_NUM_THREADS=1\n',
'max_wallclock_seconds': 3600*10 # 10 hours max runtime
})
builder.options = options_AMD2
builder.structure = struc
builder.metadata.label = 'scf:fcc_Fe_alat=Cu'
builder.calc_parameters = Dict(dict=kkrparams(NSPIN=2, LMAX=3, RMAX=10, GMAX=100, KEXCOR=4))
# scf0_GGA = submit(builder)
scf0_GGA = load_node('b79e0d1f-d2a2-4363-8bf4-613077453039') # fcc
# plot_kkr(scf0_GGA, silent=False)
print('GGA:', scf0_GGA.outputs.last_calc_out['magnetism_group']['total_spin_moment'], 'mu_B')
scale_factors = np.sort([0.8, 0.85, 0.9] + list(np.linspace(0.91, 1.09, 19)) + [1.1, 1.15, 1.2])
# run scf calculations for different lattice constants (+/-20%)
scf_group_scaled = load_group(label='fcc_Fe_GGA_scaled:scf')
builders = []
# prepare process builders for missing jobs
for scalefac in scale_factors:
label = f'fcc_Fe_GGA_scaled={scalefac:.2f}'
group_labels = [i.label for i in scf_group_scaled.nodes]
if label not in group_labels:
# scale the structure
cell = np.array(scf0_GGA.inputs.structure.cell)
asestruc = scf0_GGA.inputs.structure.get_ase()
asestruc.set_cell(cell*scalefac, scale_atoms=True)
struc_scaled = StructureData(ase=asestruc)
# submit calculation and add to group
builder = scf0_GGA.get_builder_restart()
builder.metadata.label = label
builder.structure = struc_scaled
builders.append(builder)
# now submit is any jobs are missing
if len(builders)>0:
submit_batch_iffslurm_all_free(builders, allowed_queues=['th1-2020-32', 'th1-2020-64'], dry_run=False, calculations_group=scf_group_scaled)
print('submitted', label)
from aiida_kkr.workflows import kkr_dos_wc
dos_group_scaled = load_group(label='fcc_Fe_GGA_scaled:dos')
builders = []
# prepare process builders for missing jobs
for scalefac in scale_factors:
label = f'fcc_Fe_GGA_scaled={scalefac:.2f}'
group_labels = [i.label for i in dos_group_scaled.nodes]
if label not in group_labels:
scf = [i for i in scf_group_scaled.nodes if i.label==label][0]
if scf.is_finished_ok:
builder = kkr_dos_wc.get_builder()
builder.kkr = scf.inputs.kkr
builder.options = scf.inputs.options
builder.remote_data = scf.outputs.last_RemoteData
builder.wf_parameters = Dict(dict={'nepts': 96, 'tempr': 200.0, 'emin': -10.0, 'emax': 5.0, 'kmesh': [50, 50, 50]})
builder.metadata.label = label
builders.append(builder)
# now submit is any jobs are missing
if len(builders)>0:
submit_batch_iffslurm_all_free(builders, allowed_queues=['th1-2020-32', 'th1-2020-64'], dry_run=False, calculations_group=dos_group_scaled)
# run Jij calculations
jij_calcs_group_scaled = load_group(label='fcc_Fe_GGA_scaled:Jij')
builders = []
# prepare process builders for missing jobs
for scalefac in scale_factors:
label = f'fcc_Fe_GGA_scaled={scalefac:.2f}'
group_labels = [i.label for i in jij_calcs_group_scaled.nodes]
if label not in group_labels:
scf = [i for i in scf_group_scaled.nodes if i.label==label][0]
if scf.is_finished_ok:
builder = make_jij_builder(scf)
builder.metadata.label = label
builders.append(builder)
# now submit is any jobs are missing
if len(builders)>0:
submit_batch_iffslurm_all_free(builders, allowed_queues=['th1-2020-32', 'th1-2020-64'], dry_run=False, calculations_group=jij_calcs_group_scaled)
print('submitted', label)
# extract jij nodes such that they are understood by AiiDA-Spirit
from cpa_jijs import parse_jij_calc
jij_data_group_scaled = load_group(label='fcc_Fe_GGA_scaled:jijdata')
for scalefac in scale_factors:
label = f'fcc_Fe_GGA_scaled={scalefac:.2f}'
group_labels = [i.label for i in jij_data_group_scaled.nodes]
if label+'_jijdata' not in group_labels:
# find calculation and parse the outcome
jij_calc = [i for i in jij_calcs_group_scaled.nodes if i.label==label][0]
out_dict = parse_jij_calc(jij_calc.outputs.retrieved)
# expand to rectangular structure
out_big = expand_Jijs_big_cell(jij_data=out_dict['jij_data'], struc=find_parent_structure(jij_calc))
jij_data_big, struc_big = out_big['jij_data'], out_big['structure']
# add label to find the nodes in the group
jij_data_big.label = label+'_jijdata'
struc_big.label = label+'_jij_structure'
# add to nodes group
jij_data_group_scaled.add_nodes([jij_data_big, struc_big])
from matplotlib import cm
from masci_tools.util.constants import BOHR_A
from scipy.optimize import curve_fit
def fitfunc(x, a, b, c):
return a + b*x + c*x**2
# get spin moments and scale factors
mu_s = np.array([(i.label.split('=')[1], i.outputs.output_parameters['magnetism_group']['total_spin_moment']) for i in jij_calcs_group_scaled.nodes], dtype=float)
mu_s = mu_s[mu_s[:,0].argsort()]
mu_s = mu_s[mu_s[:,0]>=0.90]
mu_s = mu_s[mu_s[:,0]<=1.1]
# find structure and from there the lattice constant
struc = find_parent_structure([i for i in jij_calcs_group_scaled.nodes if '=1.00' in i.label][0])
a_fac = struc.cell[0][0]*2
# get total energies and fit high spin and low spin parabolas
etot_scf = np.array([(i.label.split('=')[1], i.outputs.last_calc_out['energy']) for i in scf_group_scaled.nodes], dtype=float)
etot_scf = etot_scf[etot_scf[:,0].argsort()]
etot_scf[:,1] -= etot_scf[etot_scf[:,0]==1,1]
xall = etot_scf[2:-2,0]
yall = etot_scf[2:-2,1]
popt1, pcov1 = curve_fit(fitfunc, xall[xall<1]*a_fac, yall[xall<1])
popt2, pcov2 = curve_fit(fitfunc, xall[xall>=1]*a_fac, yall[xall>=1])
symbols = ['s', 'o', 'v', '^', '>', '<']
fontsize = 'large'
fontsize2 = 'large'
plt.figure(figsize=(12,6))
# total energies and fits for high spin / low spin
plt.subplot(2,2,1)
plt.plot(a_fac*etot_scf[2:-2,0], etot_scf[2:-2,1], 's', color='C1', ms=8, alpha=0.66, label='total energy', zorder=100)
x = np.linspace(0.89, 1.02, 1000)*a_fac
plt.plot(x, fitfunc(x, popt1[0], popt1[1], popt1[2]), color='C2', lw=3, label='low-spin')
x = np.linspace(0.97, 1.12, 1000)*a_fac
plt.plot(x, fitfunc(x, popt2[0], popt2[1], popt2[2]), color='C4', lw=3, label='high-spin')
plt.legend(fontsize=fontsize2)
plt.ylabel('total energy (eV)', fontsize=fontsize)
plt.xlabel('lattice constant ($\AA$)', fontsize=fontsize)
plt.xlim(3.2,4.0)
# DOS for some hi spin and low spin calculations
plt.subplot(2,2,2)
plotopt = dict(silent=True, l_channels=False, noshow=True, nofig=True)
for ii, i in enumerate(dos_group_scaled.nodes[6:17]):
if i.is_terminated:
plot_kkr(i, color=cm.viridis_r(ii/10), label=np.round(float(i.label.split('=')[1])*a_fac, 2), **plotopt)
plt.legend(ncol=3)
plt.title('')
plt.xlim(-10,4)
plt.axhline(0, color='grey', lw=0.5)
plt.ylabel('DOS (1/eV)', fontsize=fontsize)
plt.xlabel('$E-E_F$ (eV)', fontsize=fontsize)
plt.annotate("HS", fontsize='large',
xy=(-0.4, 0.5), xycoords='data',
xytext=(-2.5, 0.2), textcoords='data',
arrowprops=dict(arrowstyle="<-",
connectionstyle="arc3,rad=0.2"
),
)
plt.annotate("LS", fontsize='large',
xy=(0.8, 0.5), xycoords='data',
xytext=(2.5, 1), textcoords='data',
arrowprops=dict(arrowstyle="<-",
connectionstyle="arc3,rad=0.3",
color='grey'
),
)
#
plt.subplot(2,2,3)
plt.plot(struc.cell[0][0]*2*mu_s[:,0], mu_s[:,1], 'o--', lw=3, ms=10)
plt.xlabel('lattice constant ($\AA$)', fontsize=fontsize)
plt.ylabel('spin moment ($\mu_B$)', fontsize=fontsize)
plt.ylim(0)
plt.xlim(3.2,4.0)
plt.axvspan(3.378, 3.59, color='grey', ls='-', lw=1, alpha=0.3, zorder=-100)
plt.annotate('LS', (3.47, 2.0), zorder=100, bbox=dict(boxstyle="round,pad=0.3", fc="w", ec="w", lw=1))
plt.axvspan(3.59, 4.1, color='C3', ls='-', lw=1, alpha=0.3, zorder=-100, hatch='//')
plt.annotate('HS', (3.78, 2.0), zorder=100, bbox=dict(boxstyle="round,pad=0.3", fc="w", ec="w", lw=1))
plt.subplot(2,2,4)
isort = np.argsort([i.label for i in jij_data_group_scaled.nodes if 'jijdata' in i.label])
jijs_all = np.array([i for i in jij_data_group_scaled.nodes if 'jijdata' in i.label])[isort]
for ii, jij_big in enumerate(jijs_all):
scale = float(jij_big.label.split('=')[1].split('_')[0])
if scale>=0.94 and scale<1.05:
pos_big = jij_big.get_array('positions_expanded')
jij_big = jij_big.get_array('Jij_expanded')
plt.plot(np.linalg.norm(pos_big, axis=1), jij_big[:,5], '--', label=np.round(scale*a_fac, 2), lw=1, ms=8, color=cm.viridis_r((ii-6)/(len(mu_s)-8)), marker=symbols[ii%6])
plt.axhline(0, color='grey', lw=0.5)
plt.legend(ncol=3)
plt.ylabel('$J_{ij}/$ (meV/$\mu_B$)', fontsize=fontsize)
plt.xlabel('distance ($\AA$)', fontsize=fontsize)
plt.tight_layout()
plt.annotate('(a)', (0.005, 0.95), xycoords='figure fraction', fontsize='x-large', weight="bold")
plt.annotate('(b)', (0.505, 0.95), xycoords='figure fraction', fontsize='x-large', weight="bold")
plt.annotate('(c)', (0.005, 0.46), xycoords='figure fraction', fontsize='x-large', weight="bold")
plt.annotate('(d)', (0.505, 0.46), xycoords='figure fraction', fontsize='x-large', weight="bold")
plt.show()
starting from random spins, FM, or spin-spiral state
spirit_calcs_group_scaled = load_group(label='fcc_Fe_GGA_scaled:spirit_40x40x40_open')
from aiida.orm import ArrayData
# spirit_calc0 = [i for i in spirit_calcs_group_scaled.nodes if i.label=='fcc_Fe_GGA_scaled=1.00'][0]
# init = ArrayData()
# init.set_array('initial_state', spirit_calc0.outputs.magnetization.get_array('final'))
init = load_node('7b03e464-a6e3-4e1b-a1b5-5cd3747b32ef')
builders = []
for scalefac in tqdm(scale_factors):
label = f'fcc_Fe_GGA_scaled={scalefac:.2f}'
group_labels = [i.label for i in spirit_calcs_group_scaled.nodes]
if label not in group_labels:
# find calculation and jijdata and structure
scf = [i for i in scf_group_scaled.nodes if i.label==label][0]
jij_data_big = [i for i in jij_data_group_scaled.nodes if i.label==label+'_jijdata'][0]
struc_big = [i for i in jij_data_group_scaled.nodes if i.label==label+'_jij_structure'][0]
# prepare spirit calculation builders
if scf.outputs.last_calc_out['magnetism_group']['spin_moment_per_atom'][0]>1e-3:
builder = prepare_spirit_calc_4(scf, jij_data_big, struc_big, label=label, FM=False)
builder_FM = prepare_spirit_calc_4(scf, jij_data_big, struc_big, label=label+'_FM', FM=True)
builders.append(builder)
builders.append(builder_FM)
if label+'_SS' not in group_labels and '1.00' not in label:
# find calculation and jijdata and structure
scf = [i for i in scf_group_scaled.nodes if i.label==label][0]
jij_data_big = [i for i in jij_data_group_scaled.nodes if i.label==label+'_jijdata'][0]
struc_big = [i for i in jij_data_group_scaled.nodes if i.label==label+'_jij_structure'][0]
# prepare spirit calculation builders
if scf.outputs.last_calc_out['magnetism_group']['spin_moment_per_atom'][0]>1e-3:
builder = prepare_spirit_calc_4(scf, jij_data_big, struc_big, label=label+'_SS', FM=False)
builder.initial_state = init
builders.append(builder)
# now submit is any jobs are missing
if len(builders)>0:
submit_batch_iffslurm_all_free(builders, dry_run=False, calculations_group=spirit_calcs_group_scaled)
else:
print('nothing to submit')
# run spirit calculations with modified first Jij parameter
spirit_calcs_jij_mod = load_group(label='fcc_Fe_GGA_scaled:spirit_changed_jij')
# create stating point for AFM structure
spirit_calc0 = [i for i in spirit_calcs_group_scaled.nodes if i.label=='fcc_Fe_GGA_scaled=1.00'][0]
init_AFM = make_AFM_start(spirit_calc0)
# load original data
scf = [i for i in scf_group_scaled.nodes if i.label=='fcc_Fe_GGA_scaled=1.00'][0]
struc_big = [i for i in jij_data_group_scaled.nodes if i.label=='fcc_Fe_GGA_scaled=1.00_jij_structure'][0]
jij_data_big = [i for i in jij_data_group_scaled.nodes if i.label=='fcc_Fe_GGA_scaled=1.00_jijdata'][0]
jij = jij_data_big.get_array('Jij_expanded')
pos = jij_data_big.get_array('positions_expanded')
# prepare builders
builders = []
for jij0 in tqdm(np.linspace(-5,15,21)):
# prepare modified Jij's
jij[abs(np.linalg.norm(pos, axis=1)-np.linalg.norm(pos, axis=1).min())<1e-3,5] = jij0
jij_data = ArrayData()
jij_data.set_array('Jij_expanded', jij)
jij_data.set_array('positions_expanded', pos)
# set up the builders
group_labels = [i.label for i in spirit_calcs_jij_mod.nodes]
label = f'Jij0={jij0:.2f}meV'
if label not in group_labels:
builder = prepare_spirit_calc_4(scf, jij_data, struc_big, label=label, FM=False)
builder_FM = prepare_spirit_calc_4(scf, jij_data, struc_big, label=label+'_FM', FM=True)
builders.append(builder)
builders.append(builder_FM)
# start from spin spiral state
if label + '_SS' not in group_labels:
builder = prepare_spirit_calc_4(scf, jij_data, struc_big, label=label+'_SS', FM=False)
builder.initial_state = init
builders.append(builder)
# start from AFM state
if label + '_AFM' not in group_labels:
builder = prepare_spirit_calc_4(scf, jij_data, struc_big, label=label+'_AFM', FM=True)
# overwrite with AFM starting structure
builder.initial_state = init_AFM
builders.append(builder)
# now submit is any jobs are missing
if len(builders)>0:
submit_batch_iffslurm_all_free(builders, dry_run=False, calculations_group=spirit_calcs_jij_mod)
else:
print('nothing to submit')
def get_calc(scale_fac):
return [i for i in spirit_calcs_group_scaled.nodes if i.label==f'fcc_Fe_GGA_scaled={scale_fac:.2f}'][0]
plt.figure(figsize=(6,6))
x, fft_x, fft_y, fft_z, px, mx, py, my, pz, mz = fft_central_spins(get_calc(1.00))
plt.subplot(2,1,1)
plt.plot(my[:,2]+2, lw=2, label='y-dir', color='C1')
plt.plot(np.arange(len(mz)//2)*2, mz[::2,2]-2, '--', lw=2, color='C0')
plt.plot(np.arange(len(mz)//2)*2+1, mz[1::2,2]-2, ':', lw=2, color='C0')
plt.plot(mz[:,2], lw=2, label='z-dir', color='C0')
plt.ylabel('$S_z$', fontsize='large')
plt.xlabel('site index', fontsize='large')
plt.xlim(0,80)
plt.subplot(2,1,2)
ff = fft_y
plt.plot(x, ff/ff.max(), lw=2, label='y-dir', color='C1')
ff = fft_z
plt.plot(x, ff/ff.max(), lw=2, label='z-dir', color='C0')
plt.xlim(-1,1)
plt.ylabel('$|\mathcal{F}(q)|$', fontsize='large')
plt.xlabel('$q$ ($2\pi/a_\mathrm{lat}$)', fontsize='large')
plt.legend(loc=0, fontsize='large')
plt.tight_layout()
plt.show()
def get_scalefac(node):
return float(node.label.replace('_FM','').split('=')[1].split('_')[0])
def get_energy(node):
return node.outputs.energies.get_array('energies')[-1,1]
# start FM
nodes = [i for i in spirit_calcs_group_scaled.nodes if '_FM' in i.label and i.is_finished_ok]
e_FM = np.array([(get_scalefac(i), get_energy(i)) for i in nodes]); e_FM = e_FM[e_FM[:,0].argsort()]
e_FM = e_FM[e_FM[:,0]<=1.1]; e_FM = e_FM[e_FM[:,0]>=0.9]
# start random
nodes = [i for i in spirit_calcs_group_scaled.nodes if '_FM' not in i.label and 'SS' not in i.label and i.is_finished_ok]
e = np.array([(get_scalefac(i), get_energy(i)) for i in nodes]); e = e[e[:,0].argsort()]
e = e[e[:,0]<=1.1]; e = e[e[:,0]>=0.9]
# spin spiral ground state
nodes = [i for i in spirit_calcs_group_scaled.nodes if 'SS' in i.label and i.is_finished_ok] + [i for i in spirit_calcs_group_scaled.nodes if '1.00'==i.label[-4:]]
e_SS = np.array([(get_scalefac(i), get_energy(i)) for i in nodes]); e_SS = e_SS[e_SS[:,0].argsort()]
e_SS = e_SS[e_SS[:,0]<=1.1]; e_SS = e_SS[e_SS[:,0]>=0.9]
plt.figure(figsize=(5.7,4))
plt.plot(e_FM[:,0]*a_fac, e_FM[:,1], 'o-', label='$E_\mathrm{FM}$')
plt.plot(e[:,0]*a_fac, e[:,1], 'o--', label='$E$')
plt.ylabel('energy (meV / spin)')
plt.xlabel('lattice constant ($\AA$)', )
plt.ylim(-175, 25)
legend = plt.legend(loc=5)
legend.get_frame().set_alpha(1)
plt.axvspan(a_fac*0.975, a_fac*0.995, color='grey', ls='-', lw=1, alpha=0.3, zorder=-100)
plt.axvspan(a_fac*1.005, 4, color='C3', ls='-', lw=1, alpha=0.3, hatch='//', zorder=-100)
plt.xlim(3.4,4)
plt.annotate('SS', (3.45, -160), zorder=100)
plt.annotate('AFM', (3.525, -160), zorder=100)
plt.annotate('FM', (3.65, -160), zorder=100, bbox=dict(boxstyle="round,pad=0.3", fc="w", ec="w", lw=1, alpha=0.8))
plt.tight_layout()
plt.show()
ys = np.array([i for i in scale_factors if i>=0.94 and i<=1.03])
fs = []
for s in tqdm(ys):
x, fft_x, fft_y, fft_z, px, mx, py, my, pz, mz = fft_central_spins(get_calc(s))
fs.append(fft_x + fft_y + fft_z)
fs = np.array(fs)
plt.figure()
for i, s in enumerate(ys[::-1]):
ff = fs[::-1][i]
ff = (ff-ff.min())/(ff.max()-ff.min())
ls = '-'
if i in [0,1,2]:
ls=':'
elif i in [4,5]:
ls='--'
name = str(np.round(a_fac*s,2))
offset = 1*(100-s*100)
dat = ff
plt.plot(x, offset + dat, label=name, ls=ls, lw=2)
plt.annotate(name, (0.45, offset+0.25), fontsize='large')
if i==3:
xpos = -0.179; ypos = 0.20
plt.annotate('', (xpos, offset+ypos), (xpos, offset+1+ypos), arrowprops={'arrowstyle': '->', })
if i==6:
xpos = -0.128; ypos = 0.6
plt.annotate('', (xpos, offset+ypos), (xpos, offset+1+ypos), arrowprops={'arrowstyle': '->', })
if i==7:
xpos = -0.153; ypos = 0.5
plt.annotate('', (xpos, offset+ypos), (xpos, offset+1+ypos), arrowprops={'arrowstyle': '->', })
if i==8:
xpos = -0.179; ypos = 0.5
plt.annotate('', (xpos, offset+ypos), (xpos, offset+1+ypos), arrowprops={'arrowstyle': '->', })
if i==9:
xpos = -0.195; ypos = 0.8
plt.annotate('', (xpos, offset+ypos), (xpos, offset+1+ypos), arrowprops={'arrowstyle': '->' })
plt.annotate('$a_{\mathrm{lat}}$', (0.45, 7.3), fontsize='x-large')
plt.xlim(-1.05,1.05)
plt.ylim(-4,8.5)
plt.xlabel('$q$ $(2\pi/a_\mathrm{lat})$', fontsize='large')
plt.ylabel('normalized FFT amplitude', fontsize='large')
plt.show()
nodes_FM = [i for i in spirit_calcs_jij_mod.nodes if '_FM' in i.label and i.is_finished_ok]
e_FM = np.array([(float(i.label.split('=')[1].split('meV')[0]), get_energy(i)) for i in nodes_FM])
e_FM = e_FM[e_FM[:,0].argsort()]
nodes_AFM = [i for i in spirit_calcs_jij_mod.nodes if 'AFM' in i.label and i.is_finished_ok]
e_AFM = np.array([(float(i.label.split('=')[1].split('meV')[0]), get_energy(i)) for i in nodes_AFM])
e_AFM = e_AFM[e_AFM[:,0].argsort()]
nodes = [i for i in spirit_calcs_jij_mod.nodes if 'SS' in i.label and i.is_finished_ok]
e_SS = np.array([(float(i.label.split('=')[1].split('meV')[0]), get_energy(i)) for i in nodes])
e_SS = e_SS[e_SS[:,0].argsort()]
nodes = [i for i in spirit_calcs_jij_mod.nodes if 'FM' not in i.label and 'SS' not in i.label and i.is_finished_ok]
e = np.array([(float(i.label.split('=')[1].split('meV')[0]), get_energy(i)) for i in nodes])
e = e[e[:,0].argsort()]
lw=2; ms=5
def make_plot1(e, e_FM, e_AFM, e_SS):
plt.plot(e_FM[:,0], e_FM[:,1], 'o--', ms=ms, lw=lw, color='C0', label='$E_\mathrm{FM}$')
plt.plot(e_AFM[:,0], e_AFM[:,1], 'o--', ms=ms, lw=lw, color='C1', label='$E_\mathrm{AFM}$')
plt.plot(e_SS[:,0], e_SS[:,1], 'o--', ms=ms, lw=lw, color='C2', label='$E_\mathrm{SS}$')
efill = []
for ie in e:
iess = e_SS[e_SS[:,0]==ie[0]]
iefm = e_FM[e_FM[:,0]==ie[0]]
ieafm = e_AFM[e_AFM[:,0]==ie[0]]
if len(iess)>0 and len(iefm)>0 and len(ieafm)>0:
emin = np.min([ie[1], iess[0][1], iefm[0][1], ieafm[0][1]])
emin2 = np.min([iefm[0][1], ieafm[0][1]])
efill.append([ie[0], emin, emin2])
efill.append([7.28, -83.3, -75.8])
efill = np.array(efill)
efill = efill[efill[:,0].argsort()]
plt.fill_between(efill[:,0], y1=efill[:,1], y2=efill[:,2], color='k', alpha=0.3)
plt.ylabel('energy (meV / spin)', fontsize='large')
plt.xlabel('$J_1$ (meV)', fontsize='large')
plt.figure(figsize=(6,6))
plt.subplot(2,1,1)
make_plot1(e, e_FM, e_AFM, e_SS)
plt.legend(fontsize='large')
plt.plot((2.5, 2.5, 10, 10, 2.5), (-100, -60, -60, -100, -100), color='k')
plt.subplot(2,1,2)
make_plot1(e, e_FM, e_AFM, e_SS)
plt.xlim(2.5, 10)
plt.ylim(-100, -60)
plt.tight_layout()
plt.show()
J0s = np.array([float(i.label.split('=')[1].split('meV')[0]) for i in spirit_calcs_jij_mod.nodes if 'SS' in i.label])
fs = []
for calc in tqdm([i for i in spirit_calcs_jij_mod.nodes if 'SS' in i.label]):
x, fft_x, fft_y, fft_z, px, mx, py, my, pz, mz = fft_central_spins(calc)
fs.append(fft_x + fft_y + fft_z)
fs = np.array(fs)
isort = J0s.argsort()
J0s = J0s[isort]
fs = fs[isort]
fig = plt.figure(figsize=(6,6))
for i, f in enumerate(fs[2:]):
# normalize
f = (f-f.min())/(f.max()-f.min()) * 1.5
# now plot with shifted curves
ls = '-'
plt.plot(x, f+i*0.5, ls=ls, lw=3, zorder=-i, alpha=0.7)
plt.annotate(f'{J0s[2:][i]}', (0.45, i*0.5+0.15), fontsize='large')
plt.xlim(-1.01, 1.01)
plt.ylim(0,10.6)
plt.annotate(f'$J_1$ (meV)', (0.4, 9.8), fontsize='x-large')
plt.xlabel('q ($2\pi/a_\mathrm{lat}$)', fontsize='large')
plt.ylabel('normalized FFT amplitude', fontsize='large')
plt.tight_layout()
plt.show()
builder = spirit_calc_GGA.get_builder_restart()
params = builder.parameters.get_dict()
params['boundary_conditions'] = [False, False, False]
params['n_basis_cells'] = [20,20,20]
builder.parameters = Dict(dict=params)
builder.metadata.label = 'test_20x20x20_open'
builder.metadata.options['max_wallclock_seconds'] = 3600*2 # increase wall clock time (full run doesn't fit)
# spirit_calc_20_open = submit(builder)
spirit_calc_20_open = load_node('bd80f3d8-7c48-460f-80d3-7cdff39e0c01')
spirit_calc_20_open
builder = spirit_calc_GGA.get_builder_restart()
params = builder.parameters.get_dict()
params['boundary_conditions'] = [False, False, False]
params['n_basis_cells'] = [50,50,50]
builder.parameters = Dict(dict=params)
builder.run_options = Dict(
dict={'simulation_method': 'llg',
'solver': 'vp', #'depondt',
# 'configuration': {'plus_z': True}
})
builder.metadata.label = 'test_50x50x50_open'
builder.metadata.options['max_wallclock_seconds'] = 3600*10 # increase wall clock time (full run doesn't fit)
# spirit_calc_50_open = submit(builder)
# spirit_calc_50_open = load_node('7cb1577f-fa50-4c8e-9f6e-e2d48059b089') # Depondt solver
# spirit_calc_50_open = load_node('d2f5f6a0-24fa-423d-9e31-68fd5a0393d0') # VP solver
# spirit_calc_50_open = submit_batch_iffslurm_all_free([builder])[0]
# spirit_calc_50_open = load_node('3d424517-12ed-4f7d-95f1-ca0faa2db533') # VP solver, submitted to next free node with submit_batch_iffslurm_all_free
spirit_calc_50_open = load_node('d685e76c-30e8-4901-850a-9f7d2349b2bb') # VP solver, submitted to next free node with submit_batch_iffslurm_all_free
spirit_calc_50_open
builder = spirit_calc_GGA.get_builder_restart()
params = builder.parameters.get_dict()
params['boundary_conditions'] = [False, False, False]
params['n_basis_cells'] = [100 ,100 ,100]
builder.parameters = Dict(dict=params)
builder.run_options = Dict(
dict={'simulation_method': 'llg',
'solver': 'vp', #'depondt',
# 'configuration': {'plus_z': True}
})
builder.metadata.label = 'test_100x100x100_open'
builder.metadata.options['max_wallclock_seconds'] = 3600*10 # increase wall clock time (full run doesn't fit)
# spirit_calc_100_open = submit_batch_iffslurm_all_free([builder])[0]
spirit_calc_100_open = load_node('1cd1b9cc-8b43-4f72-9d8f-c3735ae503c8') # VP solver, submitted to next free node with submit_batch_iffslurm_all_free
spirit_calc_100_open
builder = spirit_calc_GGA4.get_builder_restart()
params = builder.parameters.get_dict()
params['boundary_conditions'] = [False, False, False]
params['n_basis_cells'] = [20 ,20 ,20]
builder.parameters = Dict(dict=params)
builder.run_options = Dict(
dict={'simulation_method': 'llg',
'solver': 'vp', #'depondt',
# 'configuration': {'plus_z': True}
})
builder.metadata.label = 'test4_20x20x20_open'
builder.metadata.options['max_wallclock_seconds'] = 3600*10 # increase wall clock time (full run doesn't fit)
# spirit_calc4_20_open = submit_batch_iffslurm_all_free([builder])[0]
# spirit_calc4_20_open = load_node('9c70ca37-fcfb-4993-8a2f-17ce23d6cb70') # VP solver, submitted to next free node with submit_batch_iffslurm_all_free, wrong jij's
# spirit_calc4_20_open = load_node('31d9bd91-fd85-4b4b-b412-a6de1e54d068') # VP solver, new jij's, still wrong
spirit_calc4_20_open = load_node('9a0ad4a7-5a31-4f63-bf3b-fe5a4effefba') # VP solver, new jij's, maybe now
spirit_calc4_20_open
# p[:,0] *= mask
# p[:,1] *= mask
# p[:,2] *= mask
np.fft.fftshift(np.fft.fftfreq(30, d=pp[1]/30)), x
p = make_positions(spirit_calc_50_open)
m = spirit_calc_50_open.outputs.magnetization.get_array('final')
plt.figure(figsize=(12,12))
p1 = p[abs(p[:,0])==np.median(abs(p[:,0]))]
m1 = m[abs(p[:,0])==np.median(abs(p[:,0]))]
# print( len(set(abs(p1[:,2]))), [len(m1[abs(p1[:,2])==ip]) for ip in np.sort(list(set(abs(p1[:,2]))))] )
# plt.hist([len(m1[abs(p1[:,2])==ip]) for ip in np.sort(list(set(abs(p1[:,2]))))])
ii = 0
for ip in np.sort(list(set(abs(p1[:,2])))):
if len(m1[abs(p1[:,2])==ip])>=20:
ii += 1
# m1 = m1[abs(p1[:,2])==np.median(abs(p1[:,2]))]#[10:-10]
# p1 = p1[abs(p1[:,2])==np.median(abs(p1[:,2]))]#[10:-10]
m2 = m1[abs(p1[:,2])==ip]
p2 = p1[abs(p1[:,2])==ip]
# print(m2.shape)
pp = p2[:,1] - p2[:,1].min()
# plt.plot(p2[:,1], p2[:,2]); plt.axis('equal')
# plt.quiver(p2[:,1], p2[:,2], m2[:,1], m2[:,2], m2[:,0], pivot='mid', scale=50 )
# plt.plot( np.linspace(-np.pi/pp[-1], np.pi/pp[-1], len(pp)), np.fft.fftshift( abs(np.fft.fft(m1[:,2])) ) )
x = np.linspace(-np.pi/pp.max() * len(pp) / np.pi*pp[1], np.pi/pp.max() * len(pp) / np.pi*pp[1], len(pp), endpoint=False)
f = np.fft.fftshift( abs(np.fft.fft(m2[:,2])) )
f1 = np.fft.fftshift( abs(np.fft.fft(m2[:,1])) )
f2 = np.fft.fftshift( abs(np.fft.fft(m2[:,0])) )
# plt.plot(m1[:,2])
# plt.plot(m1[:,1])
# plt.plot(m1[:,0])
# plt.twiny()
# plt.plot(x, f, 'C1')
# plt.plot(x, f1, 'C2')
# plt.plot(x, f2, 'C3')
# plt.plot(x, f+f1+f2, 'C4')
plt.plot(x, f+f1+f2, )
print(ii)
# np.fft.fftfreq(len(pp), d=1/pp[1])
# xm = x[f==f.max()]
# plt.axvline(xm[0]); plt.axvline(xm[1])
# pp[1], xm, (xm/np.pi*pp[1]), 1/xm, np.pi*2 / (30/5)
for ii in range(10):
p1 = p[abs(p[:,0])==np.sort(list(set(list(abs(p[:,0])))))[30+ii]]
m1 = m[abs(p[:,0])==np.sort(list(set(list(abs(p[:,0])))))[30+ii]]
m1 = m1[abs(p1[:,2])==np.sort(list(set(list(abs(p1[:,2])))))[30+ii]]
p1 = p1[abs(p1[:,2])==np.sort(list(set(list(abs(p1[:,2])))))[30+ii]]
# print(ii, p1)
if len(p1)>10:
pp = p1[:,1] - p1[:,1].min()
# plt.plot(p1[:,1], p1[:,2]); plt.axis('equal')
# plt.quiver(p1[:,1], p1[:,2], m1[:,1], m1[:,2], m1[:,0], pivot='mid', scale=50 )
# plt.plot( np.linspace(-np.pi/pp[-1], np.pi/pp[-1], len(pp)), np.fft.fftshift( abs(np.fft.fft(m1[:,2])) ) )
x = np.linspace(-np.pi/pp[-1] * len(pp) / np.pi*pp[1], np.pi/pp[-1] * len(pp) / np.pi*pp[1], len(pp), endpoint=False)
f = np.fft.fftshift( abs(np.fft.fft(m1[:,2])) )
f1 = np.fft.fftshift( abs(np.fft.fft(m1[:,1])) )
f2 = np.fft.fftshift( abs(np.fft.fft(m1[:,0])) )
plt.plot(x, f+f1+f2, )
xm = x[f==f.max()][0]
# print(ii, xm, pp.shape)
def make_mask(spirit_calc):
nborder = [10, 10, 10]
p = make_positions(spirit_calc, flatten_array=False)
mask = np.zeros_like(p[:,:,:,0])
print(mask.shape, nborder[0],-nborder[0], nborder[1],-nborder[1], nborder[2],-nborder[2])
mask[nborder[0]:-nborder[0], nborder[1]:-nborder[1], nborder[2]:-nborder[2]] = 1
mask = mask.reshape(-1)
return mask
m = spirit_calc_50_open.outputs.magnetization.get_array('final')
mask = make_mask(spirit_calc_50_open)
m[:,0] *= mask
m[:,1] *= mask
m[:,2] *= mask
plt.plot( np.linalg.norm(m, axis=1) , '.')
spirit_calc_50_open.is_finished, spirit_calc_50_open.exit_message
!verdi process status 42589
!verdi process report 42589
init_spinview(vfr_frame_id='2')
# show_spins?
# show_spins(spirit_calc_50_open, vfr_frame_id='2', mask=mask)
# show_spins(spirit_calc4_20_open, vfr_frame_id='2')
# show_spins(spirit_calc_GGA4, vfr_frame_id='2')
show_spins(spirit_calc4_20_open, vfr_frame_id='2')
f = np.fft.fftn( spirit_calc4_20_open.outputs.magnetization.get_array('final').reshape(20,20,20,4,3) )
x = np.fft.fftshift( 2 * np.fft.fftfreq(20, 1) )
# ii = 0
# plt.figure(figsize=(16,12))
# for i in range(4):
# for j in range(3):
# ii += 1
# plt.subplot(4,3,ii)
# plt.pcolormesh( np.fft.fftshift(abs(f[0,:,:,i,j])) )
for i in range(3):
plt.figure()
plt.pcolormesh(x, x, np.fft.fftshift(np.sum(np.sum(np.sum(abs(f[:,:,:,:,:]), axis=i), axis=-1), axis=-1)), shading='gouraud', cmap='binary' ); plt.colorbar(); plt.axis('equal')
spirit_calc_group = load_group(label='spirit_calcs_fccFe_GGA')
5**3 * 10 / 3600.
125