<jupyter_start><jupyter_text>Create an array of elongation we will investigate.<jupyter_code>ds_list = np.linspace(-1, 6, 51)<jupyter_output><empty_output><jupyter_text>An "interactive job" keeps the same Lammps instance runnning, while we feed it new structures on the fly.  We use it here to supply the strained structures on after the others.<jupyter_code>job = pr.create_job(pr.job_type.Lammps, "strain",True)
job.structure = model()
job.potential = potential
job.interactive_open()
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(ds=ds)
    job.run()
    
job.interactive_close()<jupyter_output>The job strain was saved and received the ID: 10051545<jupyter_text>`job.output` collects all the output from lammps for each of the structures that we gave during the interactive run.<jupyter_code>plt.plot(ds_list,job.output.energy_tot,'b', linewidth = 3,label='total energy')
plt.title('elongation with no seperation')
plt.xlabel('Elongation [$\AA$]')
plt.ylabel('$E$ [eV]')
plt.legend(loc = 'best', fontsize = 14)<jupyter_output><empty_output><jupyter_text>Now repeat the same calculations, but introduce a "vaccum" separation between one crystal plane, instead of elastically stretching it.  We just perform a static energy calculation here, so this corresponds to the "rigid" work of separation as in Rebecca's talk.<jupyter_code>job = pr.create_job(pr.job_type.Lammps, "vac", True)
job.structure = model()
job.potential = potential
job.interactive_open()
job.calc_static()
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(ds=ds, vac=True)
    job.run()
    
job.interactive_close()<jupyter_output>The job vac was saved and received the ID: 10069276<jupyter_text>In the decohesion zone model, we assume that crack formation is governed by a simple energy balance between the energy stored in the strain field versus the energy of the two surface energy.  Once two surfaces are more favorable  crack should form.  Below we plot both energies to find that point.<jupyter_code>strain = pr.load('strain')
vac = pr.load('vac')
plt.plot(ds_list, strain['output/generic/energy_tot'], label=strain.name, linewidth = 3)
plt.plot(ds_list, vac['output/generic/energy_tot'], label="crack", linewidth = 3)
    
plt.legend()
plt.title('comparison of two mechanisms')
plt.xlabel('Elongation [$\AA$]')
plt.ylabel('$E$ [eV]')
plt.legend(loc = 'best', fontsize = 14)<jupyter_output><empty_output><jupyter_text>Size ConvergenceThe argumentation above has one problem, that the strain energy goes with the volume of the sample, but the crack formation energy only on the surface area.  Hence as the total volume increases the energy barrier per atom will decrease.  We will repeat both caluclations, but with half the system size in the strain direction.<jupyter_code>n_layer = 3
job = pr.create_job(pr.job_type.Lammps, f"strain_l{n_layer}",True)
job.structure = model()
job.potential = potential
job.interactive_open()
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(ds=ds, repeat = [2, 2, n_layer])
    job.run()
    
job.interactive_close()
job = pr.create_job(pr.job_type.Lammps, f"vac_l{n_layer}", True)
job.structure = model()
job.potential = potential
job.interactive_open()
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(ds=ds, vac=True, repeat=[2, 2, n_layer])
    job.run()
    
job.interactive_close()<jupyter_output>The job vac_l3 was saved and received the ID: 10069609<jupyter_text>The system sizes are stored for each caluclation step and we can use to normalize the separation to strains.  Similarly the atomic structure is saved and we can use its length to normalize to energy per atom.<jupyter_code>job['output/generic/cells'][:, 2,2]
label_mat=['strain_l6','vac_l6','strain_l3','vac_l3'];
i=0
jobs = [
    pr.load('strain'),
    pr.load('vac'),
    pr.load('strain_l3'),
    pr.load('vac_l3'),
]
for job in jobs:
    L =job['output/generic/cells'][:, 2,2]
    N = job['output/generic/positions'].shape[1]
    plt.plot(ds_list / L, job['output/generic/energy_tot'] / N, '+',
             label="{} N={}".format(job.name.split('_')[0],  N), linewidth = 2)

plt.legend()
plt.title('Size dependence of the results')
plt.xlabel('Strain [1]')
plt.ylabel('$E$ [eV/atom]')<jupyter_output><empty_output><jupyter_text>We can also look at the stress-strain curves  to see different maximum stresses the two mechanism support.  Unsurprisingly the crack mode supports a much lower stress.<jupyter_code>for job in jobs:
    L =job['output/generic/cells'][:, 2,2]
    N = job['output/generic/positions'].shape[1]
    plt.plot(ds_list / L, -job['output/generic/pressures'][:, 2, 2], '+',
             label="{} N={}".format(job.name.split('_')[0],  N), linewidth = 2)

plt.legend()
plt.title('size dependence of the cohesive stress')
plt.xlabel('strain [1]')
plt.ylabel('$\sigma$ [GPa]')<jupyter_output><empty_output><jupyter_text>We want to see now the effect of adding hydrogen to our structure.  We modify our `model()` to support this by copying the original structure, shifting the positions by (fractional coordinates) 1/4,1/4,1/4 and then adding hydrogen to a random subset of this new lattice.  The shift ensures that we place hydrogen on the tetrahydral interstitials.<jupyter_code>def model(el='Al', strain=None, ds=0, vac=False, n_H=0, repeat=(2,2,6), solute='H'):
    structure = pr.create_ase_bulk(el, cubic=True).repeat(repeat)
    if n_H > 0:
        h_structure = structure.copy()
        h_structure[:] = solute
        pos = h_structure.get_scaled_positions()
        h_structure.set_scaled_positions(pos + 1/4/np.array(repeat))
        structure += h_structure[np.random.choice(len(structure), replace=False, size=n_H)]
    
    if not vac:
        a3_0 = structure.cell[2, 2]
        if strain is None:
            strain = ds/a3_0

        return structure.apply_strain([0, 0, strain], return_box=True)
    else:
        structure.cell[2, 2] += ds
        return structure
model(n_H=10).plot3d()<jupyter_output><empty_output><jupyter_text>We now repeat the calculations in Nickel, just strain, strain with additional hydrogen and the vacuum with no hydrogen.<jupyter_code>n_H = 5
job = pr.create_job(pr.job_type.Lammps, f"ni_strain_l{n_layer}", True)
job.structure = model(el='Ni', n_H=0, repeat=[4,4,n_layer])
job.potential = potential
job.interactive_open()
job.calc_minimize(n_print=1)
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(el='Ni', n_H=0, ds=ds, repeat=[4, 4, n_layer])
    job.run()
    
job.interactive_close()
job.animate_structure()
plt.plot(job.output.energy_tot,linewidth=3)
plt.title('Ni-H, $N_H=5$')
plt.xlabel('separation [A]')
plt.ylabel('E [eV/atom]')
n_H = 5
job = pr.create_job(pr.job_type.Lammps, f"ni_strain_h{n_H}_l{n_layer}", True)
job.structure = model(el='Ni', n_H=n_H, repeat=[4,4,n_layer])
job.potential = potential
job.interactive_open()
job.calc_minimize(n_print=1)
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(el='Ni', n_H=n_H, ds=ds, repeat=[4, 4, n_layer])
    job.run()
    
job.interactive_close()
pr['ni_strain_h5_l3'].animate_structure()
plt.plot(job.output.energy_tot,linewidth=3)
plt.title('Ni-H, $N_H=5$')
plt.xlabel('separation [A]')
plt.ylabel('E [eV/atom]')
job = pr.create_job(pr.job_type.Lammps, f"ni_vac_l{n_layer}", True)
job.structure = model(el='Ni', repeat=[4, 4, n_layer])
job.potential = potential
job.interactive_open()
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(el='Ni', ds=ds, vac=True, repeat=[4, 4, n_layer])
    job.run()
    
job.interactive_close()
plt.plot(ds_list, pr.load('ni_strain_h5_l3').output.energy_tot, label='H')
plt.plot(ds_list, pr.load('ni_vac_l3').output.energy_tot, label='no H, vacuum')
plt.plot(ds_list, pr.load('ni_strain_l3').output.energy_tot, label='no H, elastic')
plt.legend()
plt.legend()
plt.title('Effect of H on the decohesion in nickel')
plt.xlabel('Separation [A]')
plt.ylabel('E [eV]')
plt.plot(ds_list, -pr.load('ni_strain_h5_l3').output.pressures[:, 2,2], label='H',linewidth=3)
plt.plot(ds_list, -pr.load('ni_vac_l3').output.pressures[:, 2,2], label = 'no H, vacuum',linewidth=3)
plt.plot(ds_list, -pr.load('ni_strain_l3').output.pressures[:, 2,2], label = 'no H, elastic',linewidth=3)
plt.title('Effect of H on the cohesive stress in nickel')
plt.xlabel('ds (A)')
plt.ylabel('$\sigma$ (GPa)')
plt.legend()
plt.ylim(0, None)<jupyter_output><empty_output><jupyter_text>Here we just add another set of calculations with a aluminum solute, to see if the effect is hydrogen specific.<jupyter_code>n_H = 5
job = pr.create_job(pr.job_type.Lammps, f"ni_vac_al{n_H}_l{n_layer}", True)
job.structure = model(el='Ni', n_H=n_H, repeat=[4,4,n_layer], solute='Al')
job.potential = potential
job.interactive_open()
job.calc_minimize(n_print=1)
job.interactive_enforce_structure_reset=True
for ds in ds_list:
    job.structure = model(el='Ni', n_H=n_H, ds=ds, vac=True, repeat=[4, 4, n_layer], solute='Al')
    job.run()
    
job.interactive_close()
plt.plot(ds_list, pr.load('ni_vac_h5_l3').output.energy_tot, label='Ni-H',linewidth=3)
plt.plot(ds_list, pr.load('ni_vac_al5_l3').output.energy_tot, label='Ni-Al',linewidth=3)
plt.plot(ds_list, pr.load('ni_strain_l3').output.energy_tot, label='no H, elastic')
plt.plot(ds_list, pr.load('ni_vac_l3').output.energy_tot, label='pure Ni',linewidth=3)
plt.title('Effect of H and Al soluteson the decohesion in nickel')
plt.xlabel('ds (A)')
plt.ylabel('$E (eV)')
plt.legend()
plt.plot(ds_list, -pr.load('ni_vac_h5_l3').output.pressures[:, 2,2], label='Ni-H',linewidth=3)
plt.plot(ds_list, -pr.load('ni_vac_al5_l3').output.pressures[:, 2,2], label='Ni-Al')
plt.plot(ds_list, -pr.load('ni_vac_l3').output.pressures[:, 2,2], label = 'pure Ni')
plt.plot(ds_list, -pr.load('ni_strain_l3').output.pressures[:, 2,2], label = 'pure Ni, elastic',linewidth=3)
plt.legend()
plt.title('Effect of H and Al solutes on the cohesive stress in nickel')
plt.xlabel('ds (A)')
plt.ylabel('$\sigma$ (GPa)')
plt.legend()
plt.ylim(0, None)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14566156262269336, 'fraction_numerical': 0.016195524146054182, 'mean_word_length': 5.565077319587629, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '76985', 'license': 'MIT', 'n_tokens_mistral': 3897, 'n_tokens_neox': 3632, 'n_words': 875, 'path': 'notebooks/HandsOnDecohesion.ipynb', 'repo_name': 'srmnitc/structural_analysis_ADIS_2020', 'size': '10188'}