Introduction to atomman: Semidiscrete variational Peierls-Nabarro model

2.1. Disregistry and dislocation density

The disregistry of a dislocation, \(\delta_i\), is defined as the difference in the displacement just above and below the slip plane taken as a function of the coordinate \(\xi\) that is along the slip plane and perpendicular to the dislocation line.

A dislocation density, \(\rho_i\), is also defined that can be thought of as the infinitesimal Burgers vector of the dislocation at any specific \(\xi\) point, and is related to \(\delta_i\) as

For a slip plane containing a single dislocation in an otherwise infinite medium, \(\delta_i\) and \(\rho_i\) can be related to the dislocation’s Burgers vector, \(b_i\)

2.2. Semidiscrete variational Peierls-Nabarro model

The semidiscrete variational Peierls-Nabarro (SDVPN) model expresses the energy of a dislocation as a summation of different energy terms, each of which is a function of \(\delta_i\). Here, the general dislocation energy equation is

Each of the energy terms is expressed with a discrete formulation that assumes the disregistry is a step function evaluated at certain points and constant in between. The disregistry is evaluated at regular intervals, with \(\Delta \xi\) typically corresponding to a spacing of atomic columns along the \(\xi\) direction giving the calculation a discreteness corresponding to an actual atomic discreteness. The SDVPN model then solves for the disregistry at the evaluated \(\xi\) points that minimizes the dislocation’s energy.

2.3. Elastic energy

The elastic energy is due to the elastic energy of the dislocation. It contains two components: a configuration-dependent part associated with the interaction of the dislocation at one position interacting with the dislocation at another point, and a configuration-independent part associated with the dislocation’s long-range strain field. The elastic energy is expressed as

where \(i\) and \(j\) are integer indexes of the \(\xi\) evaluation points, i.e. \(\xi[i+1] = \xi[i] + \Delta \xi\), \(K_{lm}\) is the dislocation’s energy coefficient tensor, \(L\) is a cutoff radius for the long-range configuration-independent elastic energy, and \(\chi\) is a multiplicative factor related to the distance between points \(\xi[i]\) and \(\xi[j]\) given by

The dislocation density, \(\rho_l\), is obtained from the discrete \(\delta_l\) as

NOTE: slightly different forms for \(U_{elastic}\) are found in the references due to how the different authors define the energy coefficient tensor, \(K_{lm}\). Here, \(K_{lm}\) is defined such that in the isotropic case \(K_{33} = K_{screw} = \mu\), and \(K_{11} = K_{edge} = \frac{\mu}{1-\nu}\) where \(\mu\) is the shear modulus and \(\nu\) is the Poisson’s ratio.

2.4. Misfit energy

The misfit energy arises from the disregistry within the dislocation spreading coinciding with crystallographically unfavorable configurations. This energy is calculated by mapping the dislocation’s disregistry onto a generalized stacking fault energy map, i.e \(\gamma\) surface

2.5. Stress energy

The stress energy computes the energy due to an externally applied stress, \(\tau_{ij}\), interacting with the dislocation line. Two forms have been used in SDVPN models: the original full stress form

and a truncated form introduced by [8]

The two should be equivalent except for a configuration-independent constant.

2.6. Surface energy

The surface energy is a correction term introduced by [5] to account for the surface effect along the slip plane. The discrete gradient formulation by [6] is used by atomman

where \(\beta_{lj}\) is a coefficient tensor. Terms in the coefficient tensor can be

2.7. Nonlocal energy

The nonlocal energy is an alternate correction term introduced by [4] to account for the nonlocal interactions directly. It is of the form

where \(\alpha_m\) are correction coefficients corresponding to increasing separations between the evaluated \(\xi\) points. Typically, only the first coefficient \(\alpha_1\) is used and it is fitted by minimizing the discrepancy between the SDVPN disregistry and a full atomistic disregistry.

4.1. Load \(\gamma\) surface results

The atomman.defect.GammaSurface class loads and represents generalized stacking fault data.

See 4.5. Gamma surface plotting Notebook for more details on this class.

[2]:

# Load pre-calculated stacking fault energy data (same data as used in Gamma Surface Plotting Notebook)
model = Path('files', '2008--Mendelev-M-I--Al--A1--Cu--fcc--111sf.json')
gamma = am.defect.GammaSurface(model=model)

box = gamma.box
print('a1vect, a2vect used:')
print(gamma.a1vect)
print(gamma.a2vect)

a1vect, a2vect used:
[ 0.  -0.5  0.5]
[ 0.5 -0.5  0. ]

Plot smallest cell

[3]:

fig = gamma.E_gsf_surface_plot()

Plot common cell representation

[4]:

plotting_a1vect = np.array([1, 1, -2]) / 2
plotting_a2vect = np.array([1,-1, 0]) / 2
fig = gamma.E_gsf_surface_plot(a1vect = plotting_a1vect,
                               a2vect = plotting_a2vect)

4.2. Get a Volterra dislocation solution

The model uses Volterra solutions to compute the interaction stresses. A Volterra solution can easily be found using the solve_volterra_dislocation.

See 4.4. Dislocation solution and generator Notebook for more details on the solve_volterra_dislocation function.

Note: The slip plane used here must match the slip plane of the gamma surface, otherwise SDVPN will not know how to relate the two!

[5]:

# Elastic constants for Al using the 2008--Mendelev-M-I--Al potential
C11 = uc.set_in_units(105.09, 'GPa')
C12 = uc.set_in_units( 59.46, 'GPa')
C44 = uc.set_in_units( 30.66, 'GPa')
C = am.ElasticConstants(C11=C11, C12=C12, C44=C44)
print('Cij in GPa:')
print(uc.get_in_units(C.Cij, 'GPa'))

Cij in GPa:
[[105.09  59.46  59.46   0.     0.     0.  ]
 [ 59.46 105.09  59.46   0.     0.     0.  ]
 [ 59.46  59.46 105.09   0.     0.     0.  ]
 [  0.     0.     0.    30.66   0.     0.  ]
 [  0.     0.     0.     0.    30.66   0.  ]
 [  0.     0.     0.     0.     0.    30.66]]

[6]:

# Miller Burgers vector
burgers_uvw = 0.5 * np.array([1, 0, -1])

# Miller slip plane
slip_hkl = np.array([1, 1, 1])

# The dislocation's line direction (ξ_uvw) determines character
ξ_uvw = [ 1,-2, 1] # 90 degree edge

volterra = am.defect.solve_volterra_dislocation(C, burgers_uvw, ξ_uvw=ξ_uvw, slip_hkl=slip_hkl, box=box)

4.3. Initializing an SDVPN object

The Volterra dislocation solution and gamma surface provide all of the materials properties associated with the dislocation problem. During initialization, verification checks will be performed to make certain that the dislocation solution and gamma surface are compatible.

There are many parameters that can be set during initialization. The first three are the most important, and the remainder are associated with settings listed in the theory section above. All the settings parameters (tau and below) can also be changed later.

Parameters

• volterra (atomman.defect.VolterraDislocation, optional) The elastic solution for a Volterra dislocation to use as the basis of the model. Either volterra or model are required, and both cannot be given at the same time.

• gamma (atomman.defect.GammaSurface, optional) The gamma surface to use for the solution. Required unless model is given and the model content contains gamma surface data.

• model (str or DataModelDict, optional) Saved data from previous SDVPN runs to load. Either volterra or model are required, and both cannot be given at the same time.

• tau (numpy.ndarray, optional) A (3,3) array giving the stress tensor to apply to the system using the stress energy term. Only the xy, yy, and yz components are used. Default value is all zeros.

• alpha (list of float, optional) The alpha coefficient(s) used by the nonlocal energy term. Default value is [0.0].

• beta (numpy.ndarray, optional) The (3,3) array of beta coefficient(s) used by the gradient energy term. Default value is all zeros.

• cutofflongrange (float, optional) The cutoff distance to use for computing the long-range energy. Default value is 1000 Angstroms.

• fullstress (bool, optional) Flag indicating which stress energy algorithm to use. Default value is True.

• cdiffelastic (bool, optional) Flag indicating if the dislocation density for the elastic energy component is computed with central difference (True) or simply neighboring values (False). Default value is False.

• cdiffsurface (bool, optional) Flag indicating if the dislocation density for the surface energy component is computed with central difference (True) or simply neighboring values (False). Default value is True.

• cdiffstress (bool, optional) Flag indicating if the dislocation density for the stress energy component is computed with central difference (True) or simply neighboring values (False). Only matters if fullstress is True. Default value is False.

• min_method (str, optional) The scipy.optimize.minimize method to use. Default value is ‘Powell’, which seems to do well for this task.

• min_options (dict, optional) Any options to pass on to scipy.optimize.minimize. Default value is {}.

Execution tips: Tests seem to indicate that the ‘Powell’ minimization option works well for this problem. For best results, use min_options to set xtol and ftol to smaller values and run the solver multiple times in a row (minimizer may get stuck).

[7]:

min_options = {}
min_options['disp'] = True     # will display convergence info
min_options['xtol'] = 1e-6     # smaller convergence tolerance than default
min_options['ftol'] = 1e-6     # smaller convergence tolerance than default
#min_options['maxiter'] = 2    # for testing purposes

# Initialize SDVPN object from volterra and gamma
pnsolution = am.defect.SDVPN(volterra=volterra, gamma=gamma, min_options=min_options)

5.1 Initial disregistry

pn_arctan_disregistry()

The atomman.defect.pn_arctan_disregistry() function computes the classic Peierls-Nabarro arctan disregistry for an array of points x

where \(\xi\) is the dislocation’s halfwidth and \(b_i\) is the Burgers vector.

Parameters

• xmax (float or None, optional) Maximum value of x to use. Minimum value is taken as -xmax. At least 2 of xmax, xstep, and xnum must be not None. Default value is None.

• xstep (float or None, optional) Step size to use between each x value. At least 2 of xmax, xstep, and xnum must be not None. Default value is None.

• xnum (int or None, optional) Number of x values to use. At least 2 of xmax, xstep, and xnum must be not None. Default value is None.

• burgers (numpy.ndarray, optional) The Burgers vector for the dislocation. Default value is [1, 0, 0].

• halfwidth (float, optional) The dislocation halfwidth to use. Default value is 1.

• normalize (bool, optional) If True (default), the disregistry values will be scaled such that the two endpoints differ by exactly one Burgers vector.

• shift (bool, optional) If True (default), the disregistry will range [0, 0, 0] to burgers. If False, the disregistry will range from -burgers to burgers.

Returns

• x (numpy.ndarray) The x-coordinates for the disregistry values.

• disregistry (numpy.ndarray) The disregistry vector at each x-coordinate.

NOTE: the parameters normalize and shift should both be True for use with the SDVPN model.

Check burgers vector as transformed for the SDVPN solution:

• x = edge component.

• y = normal (out-of-plane) component. Should be 0.0

• z = screw component.

[8]:

b = pnsolution.burgers
print('Transformed Burgers vector =', b)

Transformed Burgers vector = [2.86043063 0.         0.        ]

Generate arctan disregistry using burgers.

[9]:

xstep = box.a * 2**0.5 / 4  # atomic column separation along [10-1] direction
xnum = 100
halfwidth = 2

x, disreg = am.defect.pn_arctan_disregistry(xstep=xstep, xnum=xnum, burgers=b, halfwidth=halfwidth)

Set x, disregistry to the SDVPN object

[10]:

pnsolution.x = x
pnsolution.disregistry = disreg

5.2 Check initial disregistry

The SDVPN class has some built-in methods for checking the disregistry profile and the corresponding energy values.

• total_energy(), misfit_energy(), elastic_energy(), longrange_energy(), stress_energy(), nonlocal_energy() and surface_energy() return the values for the different energy terms based on the current disregistry.

• check_energies() prints out the current energies for all of the different energy terms.

• disregistry_plot() creates a simple plot of the disregistry profiles.

• E_gsf_surface_plot() extends the corresponding GammaSurface method by including a line for the disregistry profile.

• E_gsf_vs_x_plot() plots the stacking fault energy associated with the disregistry values along x versus x.

[11]:

pnsolution.check_energies()

Dislocation energy terms in eV/Å:
Misfit energy =     0.15241997612248678
Elastic energy =    -0.21141715213385232
Long-range energy = 2.3123103911189946
Stress energy =     -0.0
Surface energy =    0.0
Nonlocal energy =   0.0
Total energy =      2.253313215107629

[12]:

pnsolution.disregistry_plot()
plt.show()

[13]:

pnsolution.E_gsf_surface_plot(a1vect=plotting_a1vect, a2vect=plotting_a2vect, energyperarea_unit='mJ/m^2')
plt.show()

[14]:

pnsolution.E_gsf_vs_x_plot(energyperarea_unit='mJ/m^2')
plt.show()

6.1. Call solve the first time

Since x and disregistry have been set using the identified initial values, the minimization will start from there. Also note all options are as specified above.

[15]:

pnsolution.solve()

Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 11
         Function evaluations: 31672

6.2. Check results

[16]:

pnsolution.check_energies()
pnsolution.disregistry_plot()
plt.show()

Dislocation energy terms in eV/Å:
Misfit energy =     0.15773037600126952
Elastic energy =    -0.25615441125278854
Long-range energy = 2.3123103911189946
Stress energy =     -0.0
Surface energy =    0.0
Nonlocal energy =   0.0
Total energy =      2.2138863558674755

6.3. Repeat minimization

The disregistry found after relaxing will replace the initial values in the class. Subsequent calls to solve() will then use the previous run’s final results as the initial starting place. Doing repeated solves is recommended as it can help further refine the results as the minimization algorithms (especially Powell) benefit from restarts.

Also, note that you can change the various options and settings for subsequent runs. This can save time if you wish to see relaxed disregistries for multiple different settings, like different stress states, as the two relaxed disregistry states are likely closer to each other than the arbitrary initial state.

[17]:

for i in range(10):
    pnsolution.solve()

Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 3790
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 3785
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 3710
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 3837
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 3951
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 4376
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 4493
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 4554
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 4503
Optimization terminated successfully.
         Current function value: 2.213886
         Iterations: 1
         Function evaluations: 4571

6.4. Check results

[18]:

pnsolution.check_energies()
pnsolution.disregistry_plot()
plt.show()

Dislocation energy terms in eV/Å:
Misfit energy =     0.15770921002385463
Elastic energy =    -0.2561337894922829
Long-range energy = 2.3123103911189946
Stress energy =     -0.0
Surface energy =    0.0
Nonlocal energy =   0.0
Total energy =      2.213885811650566

[19]:

pnsolution.E_gsf_surface_plot(a1vect=plotting_a1vect, a2vect=plotting_a2vect, energyperarea_unit='mJ/m^2')
plt.show()

[20]:

pnsolution.E_gsf_surface_plot(a1vect=[ 0.5, 0., -0.5], a2vect=[0.0,0.5,-0.5], energyperarea_unit='mJ/m^2')
plt.show()

[21]:

pnsolution.E_gsf_vs_x_plot(energyperarea_unit='mJ/m^2')
plt.show()