Pure Chromium
Contents

Multiscale Study of Pure Chromium
Abstract
In this work, we apply multiscale modeling techniques associated with Integrated Computational Materials Engineering (ICME) in order to capture the plastic behavior of pure chromium (Cr). We used Density Functional Theory (DFT) in the open source software Quantum Espresso (QE) to generate energy versus volume (EV), energy versus atomic separation (EA), and generalized stacking fault energy (GSFE) curves along with the lattice parameter, bulk modulus, and cohesive energy at equilibrium. The Modified Embedded Atom Method (MEAM) Parameter Calibration tool along with the DFT data was utilized to create a MEAM potential for Cr. Both sets of simulations were verified by comparison to values from literature, and a sensitivity analysis was conducted for the MEAM potential. These lower scale simulations are upscaled to Dislocation Dynamics (DD) computations. Information is bridged from atomistics to the microscale by tracking dislocation mobility and quantifying the stressstrain behavior. A convergence study is conducted on different sized crystal structures and dislocation mobility calculations with the Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) program. The mobility is calculated from seven simulations at varying applied shear stresses from 250 MPa to 3000 MPa, or 2500 to 30000 bar, respectively. These calculations were run with three sets of MEAM parameters to find an upper and lower bound for our data. The postprocessing program OVITO was used to visualize the centrosymmetry parameter, thermal equilibrium, and track dislocation mobility. Then, we generated position versus time plots for each simulation to calculate the dislocation velocity, dislocation mobility, and drag coefficient. Three sets of stressstrain curves were found from a single FrankRead source in Multiscale Dislocation Dynamics Plasticity (MDDP) computations.
Keywords: Chromium (Cr), Modified Embedded Atom Method (MEAM), Density Functional Theory (DFT), Integrated Computational Materials Engineering (ICME), Generalized Stacking Fault Energy (GSFE), Molecular Dynamics (MD), Dislocation Dynamics (DD), Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS), Multiscale Dislocation Dynamics Plasticity (MDDP)
Author(s): Cagle, M. S., Fonville, T. R., Kazandjian, S. L., and Sprow, B. T.
Introduction
Horstemeyer (Horstemeyer, 2018) proposed the ICME paradigm as a means to capture the multiscale phenomenon of a material by bridging information between length scales from the electron scale to the structural (continuum) scale. Multiscale modeling with the ICME paradigm offers the accuracy required to predict failure, reduce design costs, and reduce timetomarket with the simulationbased design paradigm. The fundamental steps required to achieve multiscale modeling with the ICME paradigm are as follows:
 Step 1: Downscaling, where we define the desired higherlevel “effects,” or information needed.
 Step 2: Modeling and simulation at lower length scales (including calibration and verification)
 Step 3: Upscaling, where we pass information up to meet the requirements set at the higher length scale in step 1. Iteration is required to ensure a proper connection.
 Step 4: Validation at the higher length scale through experiment or simulation.
In this work, we select pure Cr as our material of interest where we wish to model the plastic behavior at the structural scale. First, we downscaled (Step 1) the endgoal of modeling plastic behavior to the electronic scale, where we identified EV and GSFE curves as the quantities necessary to bridge information up (upscaling) to the nanoscale. We calculated the EV and GSFE curves with ab initio DFT simulations that we ran in QE. At the nanoscale, we calibrated the MEAM potential, derived from the Embedded Atom Method (EAM), using the EV and GSFE curves. The MEAM potential accounts for the angular forces between atoms and provides electronic structure and dislocation mobility information that we will need to conduct dislocation dynamics simulations at the microscale. MEAM potentials are applicable with many kinds of crystal lattices including the primary crystal structure of chromium that is body centered cubic (BCC). We used the MEAM Parameter Calibration (MPC) program with MATLAB interface created by Horstemeyer et al. (Horstemeyer et al, 2014) in combination with the molecular dynamics (MD) code known as LAMMPS to calibrate our MEAM potential. For both the DFT and MEAM simulations, we verify the results with values from the literature. In the next section, we give an overview of the theoretical models used in this work.
According to the ICME methodology (Horstemeyer, 2012), capturing plastic deformation of Cr at the macroscale requires the linking of atomistic and DD simulations at lower length scales. A convergence study is performed to determine the appropriate atomic size for our DD simulations using the calibrated MEAM potential and an applied shear stress of 1500 MPa (15000 bar). Then, we vary the MEAM potential to create an upper and lower bound. We used the three MEAM potentials to run Molecular Dynamics (MD) simulations within the LAMMPS program. Seven simulations were performed per MEAM potential with an applied shear stress ranging from 250 to 3000 MPa (2500 to 30000 bar). We visualized the simulation results with the open source visualization tool “Ovito.” In Ovito, we can study the centrosymmetry parameter, crystal structures, thermal equilibrium, and dislocation behavior. The dislocation behavior was captured and used to create position vs. time plots that we use to calculate dislocation velocity and dislocation mobility. The dislocation mobility was then passed on to run MDDP calculations. We can use the MDDP calculations with a single FrankRead source to describe the stressstrain behavior. Three stressstrain curves are generated from the MDDP calculations utilizing the three MEAM parameters used to calculate an upper and lower bound for the dislocation mobility.
Materials and Methods
Density Functional Theory (DFT)
Energy vs. Volume (EV) Curve
We plotted initial energy versus lattice parameter (EA) and energy versus volume (EV) curves from points calculated with Quantum Espresso (QE). The range of values was 80 – 120% of the assumed lattice constant of chromium (2.880). We set the initial Kpoint grid to 3 and assumed a constant energy cutoff value of 30 Rydberg. We show the initial EA and EV curves below in Figures 1 and 2, respectively.
KPOINT Convergence
Next, we performed a convergence study on a single atom of BCC chromium to determine which Kpoint grid offered a good tradeoff between lattice parameter, bulk modulus, minimum cohesive energy, and CPU time. We varied the Kpoints from 1 to 10 and plotted the results against derived parameters. We held energy cutoff constant at 30 Rydberg through each Kpoint convergence study. We show the Kpoint versus lattice parameter, bulk modulus, minimum cohesive energy, and CPU time plots below in Figures 36. We determine a Kpoint grid of 16 constituted convergence of lattice parameter (± 0.00086 A), bulk modulus (± 1.963958 Kbar), and minimum cohesive energy (± 0.000753 eV) without requiring more CPU time. Using our converged Kpoint grid, we conducted a second convergence study on a single atom of BCC chromium to determine which energy cutoff value offered a good tradeoff between lattice parameter, bulk modulus, and CPU time. We varied the energy curoff values from 20 to 100 Rydberg and then plotted the results against the derived parameters. We show the energy cutoff versus lattice parameter, bulk modulus, minimum cohesive energy, and CPU time plots below in Figures 710.
Equation of State (EOS)
We determined an energy cutoff of 60 Rydberg constituted convergence of lattice parameter (± 0.000144 A), bulk modulus (± 1.623263 Kbar), and minimum cohesive energy (± 0.000896 eV) without requiring more CPU time. Next, we compared two equations of state (EOS), namely, Birch and Murnaghan. We use an EOS to derive the lattice parameter, bulk modulus, and minimum cohesive energy. We show the results for lattice parameter, bulk modulus, and minimum cohesive energy using the Birch EOS below in Figures 1113. We do not find good convergence for our lattice parameter, bulk modulus, or minimum cohesive energy using the Birch EOS. We show the results for lattice parameter, bulk modulus, and minimum cohesive energy using the Murnaghan EOS below in Figures 1416. Again, we do not find good convergence for our lattice parameter, bulk modulus, or minimum cohesive energy using the Murnaghan EOS. However, we select the Birch EOS for because the results for the three parameters were closer to experimental values found in the literature.
Final Converged EA and EV Plots
Finally, using our converged Kpoint grid (16 points), and converged energy cutoff (60 Rydberg) values, we created a set of energy versus lattice parameter (EA) and energy versus volume (EV) curves for the primary chromium crystal structure (BCC), as well as the secondary (FCC) and tertiary (HCP) crystal structures. We show the final EA and EV curves for BCC, FCC, and HCP chromium below in Figures 1722.
Generalized Stacking Fault Energy (GSFE) Curves
We ran DFT simulations in QE to generate the GSFE curves for pure chromium using a Python script available on the ICME website. This script can generate GSFE curves along the three BCC slip systems presented in Table 1.
Name  Slip Plane  Slip Direction 

Full  (110)  [111] 
Partial  (110)  [001] 
Longpartial  (110)  [110] 
Table 2 shows the inputs to the code including the converged lattice parameter, energy cutoff, and number of kpoints from the previous energyvolume curve study. We did not conduct a convergence study for the GSFE curves assuming that the converged values from previous DFT simulations produced reliable results. These calculations also employed the same pseudopotential as the previous DFT calculations.
Parameter  Value 

Lattice Parameter  2.8497783 Å 
Energy Cutoff  60 Ry 
Number of Kpoints  16 
Element Weight  51.9961 atomic units^{[1]} 
Smear  0.06 
Vacuum Size  20.0 
Number of Stacking Layers  10 
GSFE Curves
We generated GSFE curves for pure chromium along three BCC slip systems: (110)[111] (full), (110)[001] (partial), and (110) [110] (longpartial). We computed sixteen GSFE data points for all three slip systems. Data for the longpartial slip system only spanned half of the slip distance; therefore, we assumed this curve was symmetric based on the symmetry of the BCC structure. We show the GSFE curves (change in energy versus upper atom displacement) in Figure 23 below. The Yaxis represents change in energy (mJ/m2), and the Xaxis represents displacement normalized by lattice parameter (2.8497783 Å).
Based on the GSFE curves, the full slip system presents the lowest energy barrier to slip. Therefore, slip is most likely to occur along the full slip system instead of either the partial or longpartial systems.
In Figure 24, we compare data from our full slip system GSFE curve with DFT from Yang et. al.^{[2]} The axes in Figure 24 are the same as in Figure 23 except both curves are scaled such that the displacement spans zero to one. Our DFT calculations show close agreement with Yang et. al.
Modified Embedded Atom Method Potential Calibration
MEAM potential calibration has been a manual process since the creation of the MEAM model. Horstemeyer et al. (Horstemeyer et al, 2014) presented a simple MEAM potential calibration methodology with the MPC software tool. In this work, we follow a similar approach where we calibrated BCC, FCC, and HCP energy versus lattice spacing using energy versus atomic volume (EV) data produced in QE. After loading the EV curve results from DFT, we calibrate the MEAM model parameter values of equilibrium lattice constant (alat), energy per unit volume (esub), and the constant related to the bulk modulus (alpha) for LAMMPS to match the minimum BCC energy as shown in Figure 25 below.
Sensitivity Analysis
We calculate values far away from equilibrium using the attract and repuls parameters proposed by Hughes et al. (Hughes et al., 2014) that represent the cubic attraction and repulsion terms in the UEOS. We use single crystal elastic moduli experimental values from Yang et al. (Yang et al., 2018) (Table 3) to calibrate the MEAM potential elastic constants by varying the screening parameter, Cmin, scaling factor for embedding energy, asub, and exponential decay factor for atomic density, b0. These variables were shown to have the largest effect on the EV curve calibration based on parametric sensitivity studies.
Elastic Constant  Experimental  MEAM 

C11  391.0  391.0 
C12  89.6  89.6 
C44  103.2  103.2 
Furthermore, calibration to the GSFE curve from DFT is achieved by adjusting the t3 parameter (weighting factor for atomic density). In Figure 26 below, we compare the GSFE curves obtained from DFT and the MPC tool and in Table 4, we summarize the MEAM potential constants after calibration to the DFT results.
Name  Calibrated Value  Name  Calibrated Value 

alat  2.8864  t0  1.0000 
esub  4.0750  t1  0.2075 
asub  1.0045  t2  12.2600 
alpha  5.1748  t3  0.3500 
b0  3.0545  Cmin  1.8996 
b1  1.0000  attract  0.0000 
b2  1.0000  repuls  0.0000 
b3  1.0000 
Molecular Dynamics
Atomic Position Generation
We generated the atom positions files for dislocation mobility simulations via a Fortran routine available on the ICME website Atomistic Dislocation Generation. To setup the code, we selected the body centered cubic (BCC) crystal structure to match the crystal structure of Cr. Then, we selected to study an edge dislocation for a Periodic Array of Dislocations (PAD). The atom positions application does not allow for custom lattice parameters with BCC crystal structures. For this reason, we selected the preprogrammed iron, which has the closest lattice parameter (2.85 Å) to chromium. However, the first step of each LAMMPS simulation equilibrates the system using the userprovided MEAM potential, which moves the atoms to the lowest energy state. This allows for a slight tolerance in the initial lattice parameter; simulations can be visually inspected after equilibrium to check that the initial lattice parameter does not adversely affect the resulting structure as was the case in our simulations.
Dislocation Mobility
Dislocation Velocity vs. Applied Shear Stress
To study the effect of number of atoms on dislocation mobility, we conducted six LAMMPS simulations with a constant applied shear stress of 1500 MPa (15000 bar). We varied the number of unit cells in the X and Y directions to produce the different volumes. We held the number of unit cells in Z constant at 3 unit cells because the edge dislocation is assumed to be infinitely long in this direction. Therefore, Z dimension can be accurately modeled with only a few unit cells. The number of unit cells in X and Y were selected and varied consistently to prevent edge and size effects. The X direction was the largest dimension because the dislocation translates in the X direction. In Table 5 below, we show the atom position sizes and their resulting dislocation velocity.
Simulation  (X, Y, Z) units  # of atoms  Velocity (nm/ps) 

1  11 x 10 x 3  330  0.416 
2  31 x 14 x 3  1302  0.243 
3  41 x 20 x 3  2460  0.133 
4  55 x 36 x 3  5940  0.102 
5  73 x 34 x 3  7446  0.104 
6  97 x 46 x 3  13386  0.104 
We can see the change in dislocation velocity is below 5% after simulation number 4, or 5940 atoms, therefore we consider the result to be converged. For all future calculations, we select simulation 5 as our baseline atom position file having the dimensions of 73 x 34 x 3 and a total of 7446 atoms.
MEAM Parameter Bounds
Previously, we calibrated a MEAM parameter with respect to energy vs. lattice parameter, energy vs. volume, elastic constants, and the Generalized Stacking Fault Energy (GSFE) curve. In this section, we use the calibrated MEAM parameter as a mean value and generate a set of upper and lower parameters by varying the GSFE curve. We found the t3 parameter had the largest influence on the GSFE curve; therefore, we conducted a parametric study on this parameter to find a GSFE curve that has a peak value approximately 10% higher and lower than the mean. We show the parametric study and results from the MEAM Parameter Calibration (MPC) tool in Figure 27 below.
We can see a good variation in the GSFE curve by ranging the t3 parameter from 1 to 1, where the original t3 parameter was approximately 0.435. From this parametric study, we could further refine the parameters and we selected t3=1.3 for our new lower GSFE curve, and t3=0.7 for our new upper GSFE curve. The three MEAM parameters were then combined with the atom positions file to run the DD simulations in LAMMPS.
Dislocation Velocity vs. Applied Shear Stress
In this section we study the effect of applied shear stress on dislocation velocity for the three sets of MEAM parameters to determine the dislocation mobility. We conducted seven simulations for each MEAM parameter set to study the effect of applied shear stress on dislocation velocity. The applied shear stresses ranged from 250 MPa to 3000 MPa, or 2500 to 30000 bar, respectively. From the position vs. time results, we can calculate dislocation velocity as the slope of each line. Then we convert angstroms to nanometer by multiplying each result by 0.1. In Figure 28 below, we show the results for all 21 simulations velocity versus applied shear stress.
We can calculate the dislocation mobility from the linear region of the velocity vs. applied shear stress curve where the slope is equivalent to the burgers vector divided by the drag coefficient. We know that dislocation mobility is equivalent to the inverse of the drag coefficient ^{[3]}, therefore, we must only divide the slope of the linear region by the burgers vector to obtain the dislocation mobility. From these results shown in Figure 28, we take the linear region to lie between 0 and 1000 MPa. We can then find the slope of this line and then divide the slope by the burgers vector (2.47e10) to calculate the mobility. We can also find where this line intersects the XAxis, which is the Peierls stress. In Table 6 below, we show the dislocation mobility, Peireles Stress, and drag coefficient (inverse of mobility) for each of our MEAM parameter sets.
MEAM Parameter  Mobility (1/Pa*s)  Peierls Stress (MPa)  Drag Coefficient (Pa*s) 

Lower  726  204.34  0.00138 
Mean  277  167.41  0.0036 
Upper  70.7  174.82  0.0141 
Dislocation Dynamics
MDDP with Single FrankRead Source
We generated a single Frank Read Source (FRS) pair within BCC chromium using Multiscale Dislocation Dynamics Plasticity (MDDP). This software illustrates dislocation movement throughout a crystal. The input files included values shown in Table 7. Dislocation mobility was calculated from previous LAMMPS simulations for a lower, mean, and upper MEAM parameter set.
Parameter  Value 

Crystal Size (Burger's vectors)  10000 X 10000 X 10000 
Xdirection  [1 0 0] 
Ydirection  [0 1 0] 
Zdirection  [0 0 1] 
Strain Rate  1e5 
Loading Direction  Xaxis 
Density (kg/m3)  7190 
Shear Modulus (GPa)  110 
Poisson’s Ratio  0.31 
Burger’s Vector Length (m)  2.47E10 
Temperature (K)  300 
Stacking Fault Energy (J/m2)  1.338 
Lower Bound Dislocation Mobility (1/Pa*s)  726 
Mean Dislocation Mobility (1/Pa*s)  277 
Upper Bound Dislocation Mobility (1/Pa*s)  70.7 
We illustrated dislocation loops from the Frank Read Source pairs using TecPlot 360. Several intervals are depicted in Figures 29 – 40. We ran MDDP until the FRS had completed one loop.
MDDP with multiple FrankRead Sources
We generated a multiple Frank Read Source (FRS) pairs within BCC chromium using Multiscale Dislocation Dynamics Plasticity (MDDP). This software illustrates dislocation movement throughout a crystal. The input files included values shown in Table 8. Dislocation mobility was calculated from previous LAMMPS simulations for a lower, mean, and upper MEAM parameter set.
Parameter  Value 

Crystal Size (Burger's vectors)  10000 X 10000 X 10000 
Xdirection  [1 0 0] 
Ydirection  [0 1 0] 
Zdirection  [0 0 1] 
Strain Rate  1e5 
Loading Direction  Xaxis 
# Frank Read Source Pairs  10 
Density (kg/m3)  7190 
Shear Modulus (GPa)  110 
Poisson’s Ratio  0.31 
Burger’s Vector Length (m)  2.47E10 
Temperature (K)  300 
Stacking Fault Energy (J/m2)  1.338 
Lower Bound Dislocation Mobility (1/Pa*s)  726 
Mean Dislocation Mobility (1/Pa*s)  277 
Upper Bound Dislocation Mobility (1/Pa*s)  70.7 
Determination of Voce Hardening Law
We calculated the Voce hardening law parameters by using Equation 1 [2].
Where K is hardening, α is a constant relating to strength (assumed to be 0.35 [1]), μ is the shear modulus (110 GPa for chromium), b is the Burger’s vector (2.47E10 m for chromium), and ρ is the dislocation density. We determined dislocation densities from MDDP and plugged those values into Equation 1 to calculate kappa (K). We plotted K values against time, and determined initial strength (K0), saturation strength (Ks), and initial hardening rate (h0) from the decaying exponential function (Equation 2 [2]).
K0 is the starting value of kappa versus time, and Ks is the asymptotic value in which the function converges. The hardening rate (h0) is the fitted slope in Excel of the decaying exponential graph. Figures 41, 42, and 43 shows the fitted curves to the lower bound, mean, and upper bound data from MDDP. Table 9 lists the resulting values.
Parameter  Value 

Lower Bound K0 (MPa)  52 
Mean K0 (MPa)  45 
Upper Bound K0 (MPa)  37 
Lower Bound KS (MPa)  150 
Mean KS (MPa)  600 
Upper Bound KS (MPa)  700 
Lower Bound h0 (GPa)  27.5 
Mean h0 (GPa)  25 
Upper Bound h0 (GPa)  10 
To investigate the effect of different numbers of Frank Read Source (FRS) pairs on the Voce hardening relationship, we ran several simulations with MDDP. Figures 4446 depict the change in K0, Ks, and h0 parameters with respect to varying FRS pairs. K0 seems to increase steadily while Ks remains the same. The h0 value seems to converge after 8 FRS pairs. We attempted a simulation with zero FRS; however, MDDP was not able to run. With no dislocation source initially, a single crystal has theoretical infinite strength and produces no stressstrain response.
Crystal Plasticity
Single Crystal Stress Strain Behavior
The Voce hardening law material constants are upscaled to the mesoscale via MDDP computations using lower scale mobility information. We simulate hardening of Cr within a single Crystal Plasticity Finite Element Method (CPFEM) ABAQUS model. By understanding dislocation interaction at the microscale, we can determine hardening and plastic deformation at higher length scales. The mechanical response for each MEAM parameter set was nearly identical. Figure 47 depicts the stressstrain data provided from MDDP. Each MDDP calculation ran for 30 time steps. The plotted data shows mostly linear elastic response, but dislocations are present within the crystal.
Figure 48 shows the stressstrain response from ten Frank Read Source pairs within BCC chromium from three different mobility inputs. The upper bound resulted in the largest yield point, which is expected from the lower dislocation mobility. The mean and lower bound simulations yield much sooner due to higher dislocation mobility. The yield point is much higher than chromium’s yield point at macroscale due to size effects at the lower length scales.
A usermaterial subroutine is implemented in ABAQUS to perform a single C3D8R brick element simulation outputting plots of the normal stress (S33) and logarithmic strain (LE33) components along the zaxis for uniaxial tension and compression. For the case of simple shear, the Von Mises stress is plotted against the logarithmic strain (LE13) component acting in the z coordinate direction on the plane whose outward normal is parallel to the xaxis. The dislocation drag coefficient and hardening parameters for Cr are specified within ABAQUS input files along with other quantities such as Euler angles dictating the crystal orientation. Computations are performed in uniaxial tension, uniaxial compression, and simple shear for lower bound, upper bound, and mean DD simulation results. In the results section we show the stressstrain behaviors for singlecrystal chromium evaluated at various levels of plastic strain.
Polycrystal Stress Strain Behavior
In this section we used Abaqus to generate true stresstrue strain curves from three different hardening laws for tension, compression, and torsion (simple shear) at a strain rate of roughly 1/s. Theses curves were generated for both a single crystal, single element simulations and a polycrystalline (500 grains), single element simulations. We used an elastic modulus of 210000 GPa and poisson’s ratio of 0.30 for all the polycrystalline simulations and the singlecrystal torsion simulations. Figure 49, Figure 50, and Figure 51 present the tension, compression, and torsion stressstrain curves respectively. Tensile and compressive stressstrain curves are presented as axial stresses and strains, but torsional stressstrain curves are presented as Von Mises equivalent stresses and strains. Furthermore, the shear strain in the torsion (simple shear) simulations was assumed to be roughly equivalent to the Von Mises equivalent strain.
For every stress state, the lower set of Voce hardening parameters produced the highest yield stress and vice versa. Both the tension and compression simulations were carried out to relatively high strains; however, the physical behavior of Chromium is not necessarily expected to be accurately captured at these high strains. We compared the polycrystalline simulations to the singlegrain simulations in Figure 52, Figure 53, and Figure 54.
The single crystal stressstrain curves yielded at lower stresses in tension and compression and higher stresses in torsion. This behavior may be due to the orientation of the singlegrain; more singlegrain simulations should be run with various Euler angles to understand this factor. The Hall Petch effect may also play a role in the polycrystalline simulations increased strength. In the next section we calibrate an Internal State Variable (ISV) model to the polycrystalline results to model the macroscopic stressstrain behavior of pure Chromium.
Internal State Variable Model
Calibration
In this section we use the stressstrain behavior from the polycrystalline single element CPFEM to calibrate an ISV model with the DMGfit tool and then verify our results with another set of single element analysis in Abaqus. A tutorial for using the MSU DMGfit tool to calibrate an ISV model to experimental data can be found here File:DMGfitTutorial.zip. With the tool, we can find useful material constants for the ISV model that are used in a finite element calculation as a User Material subroutine (UMAT) in Abaqus. In this work, we fit three ISV models, one for the lower, mean, and upper sets of data we have been collecting. Ideally, the ISV model would capture temperature effects, strain rate dependencies, damage, hardening, and recovery for a range of loading conditions. However, in this work, we are limited to a single strain rate, a constant temperature, and disregard cyclic hardening/recovery and damage. Upon opening the DMGfit tool, we load in a set of tension, compression, and torsion data points and then list out our material dataset parameters. The only parameters we adjusted in the dataset are listed below in Table 10. The period was calculated automatically from the assigned strain rate.
dataset parameter  value 

Units  5 = MPa 
Strain Rate (1/s)  1 
Tinit (initial temerature)  300 K 
G (Shear Modulus (MPa))  110000 
Bulk (Bulk Modulus (MPa))  251000 
Tmelt (melting temerature (K))  2180 
Finally, we imported the stress strain data from the polycrystal single element analysis for calibration. This included data for compression, tension, and torsion up to 30% strain. As we mentioned earlier, in this work we only calibrate the ISV model for the yield stress, hardening, recovery, and loading conditions; therefore, the only parameters we optimized were the Yield Stress (C3), Transition (C5), Isotropic Hardening Modulus (C15), and Isotropic Dynamic Recovery (C13). We list the calibrated parameters for these variables in Table 11 below.
Material Constant  Lower ISV  Mean ISV  Upper ISV 

Yield Stress (C3)  158.52  148.98  115.22 
Transition (C5)  1.036  1.036  1.036 
Kinematic anisotropic hardening (C9)  0  0  127.78 
Kinematic dynamic recovery (C7)  0.00217  0.00223  0.00223 
Isotropic hardening modulus (C15)  227.95  264.91  66.85 
Isotropic dynamic recovery (C13)  0.0034  0.0034  0.5146 
Verification
To verify the ISV model calibration, we updated the input decks from the original polycrystalline CPFEM single element analysis with the new set of dependent variables and material constants. Additionally, to avoid hourglassing errors in the analysis, we changed the element type from C3D8R to a fully integrated C3D8 element. We submitted nine simulations to verify the three ISV models for tension, compression, and torsional loading. In Figures 55, 56, and 57, we show the overlaid stressstrain results for the original single element analysis, the ISV model fit, and the single element analysis with the ISV user material subroutine.
The results were similar for the lower and upper datasets and are not reported here. From Figure 57, we see we need to work on the ISV fit for the torsion including material constants Ca and Cb. However, we are satisfied with the results and believe we could move forward to validate the model by experimentally testing a piece of pure Chromium for comparison to ISV model.
Theoretical Models
MEAM
We use atomistic analysis methods, suggested by Lee et al. (Lee et al, 2010), to describe the energy of a system related to interatomic forces and the relative position of atoms. According to the EAM developed by Baskes (Baskes, 1997), the total energy of a system is composed of the embedding function F_i and pair interaction ϕ terms shown in Equation 3 below:
where E is energy, ρ_i is background electron density, ϕ is atomic pair interaction, and R_ij is distance. We describe the atomic pair interactions (ϕ) at a distance R in Equation 4 below:
where Z is the number of first neighbors, E^u (R) is the background electronic density, and ρ ̅^0 (R) is energy per atom for the reference structure. We use the universal equation of state (UEOS) or first principles computations to obtain the atomic energy (E^u (R)). We calculate the UEOS parameters of cohesive energy (E_c), nearest neighbor distance (r_e), atomic volume (Ω), and bulk modulus (B) at equilibrium in the reference structure in Equations 5, 6, and 7, respectively.
The MEAM embedding function comprises the aforementioned terms and the adjustable parameter A as shown in Equation 8 below:
The background electron density is the combination of the spherically symmetric (ρ^(0)) and angular (ρ^(1)),ρ^(2),ρ^(3)) partial electron density terms:
where the terms α,β,γ represent summations over the x, y, and z coordinate directions. Atomic electron density ρ^a(h) decreases exponentially with increasing distance r^i between an atom i and the site of interest as illustrated in Equation 10 that contains constant decay lengths, β^h.
The equations for the partial electron densities are coordinate invariant and equal to zero for cubic symmetric crystals. We represent the angular contributions from the partial electron densities in simplified form in Equation 10 (with constants t^h) and the equation for background electron density in Equation 11.
In order to obtain accurate material properties and avoid imaginary electron densities for Г<1, we use the functional form of G (shown in Equation 11) where G(0)=1/2. When G(0)=1, the background electron density is non angular dependent and is equivalent to the EAM expression for ρ (Baskes, 1997).
Voce Law
We define the Voce law for materials exhibiting isotropic workhardening or stress saturation at elevated stresses or strains in Equation 2 below. This equation represents how the flow stress evolves with plastic work.
where κ_S is the saturation strength. The initial strength and hardening rate is given by κ_0 and h_0, respectively. Moreover, the constant plastic shear strain rate, C, is computed with DD (Palaparti et al.; Horstemeyer, 2012).
Conclusion
In this work, we present DFT and MEAM simulations for pure chromium. We started by generating EA and EV curves and computing the bulk modulus, lattice parameter, and minimum cohesive energy. These values converged at 8 K point grid and energy cutoff of 60 Rydberg. A bulk modulus of 251 GPa, lattice parameter of 2.845 Å, and minimum cohesive energy of 4.10 eV was found for BCC chromium. DFT simulation data is in good agreement with existing literature. First principles computations were linked to atomistics through the calibration of MEAM model constants to DFT results and elastic moduli found in literature. Information obtained using atomistic methods is passed on to dislocation dynamics simulations at the microscale. To visualize dislocation movement within chromium, we produced a structure in LAMMPS with several atoms and a dislocation present. Once we had a converged system without too much image effects, we sheared the top and bottom rows of atoms and received dislocation position and time information. From those results, we calculated the dislocation velocity and calculated the mobility from three MEAM parameter sets. The stressstrain response and dislocation loop movement from MDDP gave us crucial information to feed a crystal plasticity model. After analyzing dislocation formation and movement within a crystal given a mechanical stimulus, we create our crystal plasticity model of Cr through hardening law development.
References
 ↑ Meija, Juris, et al. "Atomic weights of the elements 2013 (IUPAC Technical Report)." Pure and Applied Chemistry 88.3 (2016): 265291.
 ↑ Yang, Keliang, et al. "Modified analytic embedded atom method potential for chromium." Modelling and Simulation in Materials Science and Engineering 26.6 (2018): 065001.
 ↑ Horstemeyer, M.F., 2012. Integrated Computational Materials Engineering (ICME) for Metals: Using Multiscale Modeling to Invigorate Engineering Design with Science. John Wiley & Sons.
[4] Atomistic Dislocation Generation, Engineering Virtual Organization for CyberDesign, March 2017. [Online]. Available: https://icme.hpc.msstate.edu/mediawiki/index.php/Atomistic_Dislocation_Generation [Accessed: 4 April. 2019].
[5] Baskes, M. (1997). Determination of modified embedded atom method parameters for nickel. Materials Chemistry and Physics, 50(2), 152158. doi:10.1016/s02540584(97)802520.
[6] B. Huddleston, GSFE Curve Generation Script, Engineering Virtual Organization for CyberDesign, March 2015. [Online]. Available: https://icme.hpc.msstate.edu/mediawiki/index.php/GSFE_Curve_Generation_Script [Accessed: 5 March. 2019].
[7] Horstemeyer, Mark F. Integrated Computational Materials Engineering (ICME) for Metals: Concepts and Case Studies. John Wiley & Sons, Inc., 2018.
[8] Horstemeyer, M. F., Hughes, J. M., Sukhija, N., Lawrimore, W. B., Kim, S., Carino, R., & Baskes, M. I. (2014). Hierarchical Bridging Between Ab Initio and Atomistic Level Computations: Calibrating the Modified Embedded Atom Method (MEAM) Potential (Part A). Jom,67(1), 143147. doi:10.1007/s1183701412440.
[9] J. Meija et al, “Atomic weights of the elements 2013 (IUPAC Technical Report),” Pure and Applied Chemistry, vol. 88, no. 3, Feb 2016, pp. 265–91.
[10] K. Yang et al, “Modified analytic embedded atom method potential for chromium,” Modelling Simul. Mater. Sci. Eng., vol. 26, no. 6, June 2018, pp. 1–15.
[11] Hughes, J. M., Horstemeyer, M. F., Carino, R., Sukhija, N., Lawrimore, W. B., Kim, S., & Baskes, M. I. (2014). Hierarchical Bridging Between Ab Initio and Atomistic Level Computations: Sensitivity and Uncertainty Analysis for the Modified EmbeddedAtom Method (MEAM) Potential (Part B). Jom, 67(1), 148153. doi:10.1007/s1183701412057.
[12] K. Yang, L. Lang, H. Deng, F. Gao, and W. Hu, “Modified analytic embedded atom method potential for chromium,” Model. Simul. Mater. Sci. Eng., vol. 26, no. 6, p. 065001, Jun. 2018.
Lee, B., Ko, W., Kim, H., & Kim, E. (2010). The modified embeddedatom method interatomic potentials and recent progress in atomistic simulations. Calphad,34(4), 510522. doi:10.1016/j.calphad.2010.10.007.
Palaparti, D. P. Rao, et al. “Influence of Temperature on Tensile Stress–Strain and Work Hardening Behaviour of Forged Thick Section 9Cr–1Mo Ferritic Steel.” Transactions of the Indian Institute of Metals, vol. 65, no. 5, 2012, pp. 413–423., doi:10.1007/s1266601201465.