PFC3D Microscopic Simulations

• PFC3D is a computational framework that models atomistic dynamics and defect nucleation in FCC and BCC crystals using a rescaled density field and phase-field-crystal formulation.
• It employs both diffusive and inertia-enhanced (MPFC) dynamical equations with spectral solvers to capture glide and climb dislocation mechanisms effectively.
• The approach bridges MD and continuum models, enabling microsecond-scale simulations that efficiently study complex three-dimensional defect morphologies.

Three-dimensional phase-field-crystal (PFC3D) microscopic simulations constitute a computational framework for modeling the atomistic dynamics of crystalline defects in face-centered cubic (FCC) and body-centered cubic (BCC) crystals over diffusive time scales. PFC3D bridges the gap between molecular dynamics (MD) simulations, which directly resolve atomic motions but are limited to nanoseconds, and continuum models, which cannot capture true atomistic phenomena. The PFC approach efficiently captures both conservative (glide) and non-conservative (climb) dislocation mechanisms, nucleation dynamics at interfaces and inclusions, and the evolution of complex three-dimensional defect configurations within a coarse-grained atomistic description (Berry et al., 2014).

1. PFC3D Free-Energy Formulation and Correlation Kernel Structure

The PFC model is built upon a rescaled atomic density field

$$n(\mathbf r) = \frac{\rho(\mathbf r)}{\rho_\ell^0} - 1$$

where $\rho(\mathbf r)$ is the local number density and $\rho_\ell^0$ is a reference liquid density. In three dimensions, the standard free energy functional for FCC symmetry in dimensionless form is

$$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$

with phenomenological parameters $w$, $u$, and a two-point correlation kernel $C_2(r)$. Its Fourier transform for FCC is

$$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$

with $\tilde k = k/(2\pi)$, $r \propto T$, $\rho(\mathbf r)$0 controlling elastic moduli, and the Gaussian term stabilizing stacking faults (parameters: $\rho(\mathbf r)$1, $\rho(\mathbf r)$2, $\rho(\mathbf r)$3). For BCC crystals, a structural/XPFC kernel is used:

$\rho(\mathbf r)$4

with each lattice plane family $\rho(\mathbf r)$5 assigned peak position $\rho(\mathbf r)$6, width $\rho(\mathbf r)$7, planar density $\rho(\mathbf r)$8, multiplicity $\rho(\mathbf r)$9, and temperature parameter $\rho_\ell^0$0 (Berry et al., 2014).

2. Dynamical Equations and Numerical Strategy

The dynamical evolution of the density field in PFC3D uses two principal schemes:

• Model B (purely diffusive):

$\rho_\ell^0$1

• MPFC (inertia-enhanced):

$\rho_\ell^0$2

where $\rho_\ell^0$3 is the damping rate, and $\rho_\ell^0$4 sets an elastic/sound-speed scale.

Spatial discretization is performed with a spectral solver. Linear terms are treated implicitly, and nonlinearities explicitly. Parameter choices for FCC and BCC are tuned to target lattice constants, elastic moduli, and stacking-fault energies. For FCC: $\rho_\ell^0$5, $\rho_\ell^0$6, $\rho_\ell^0$7, $\rho_\ell^0$8, $\rho_\ell^0$9, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$0, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$1, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$2, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$3, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$4, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$5. For BCC: $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$6, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$7, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$8, $$\tilde F[n] = \int d\mathbf r \left\{\frac{1}{2} n^2 - \frac{w}{6} n^3 + \frac{u}{12} n^4\right\} - \frac{1}{2} \int d\mathbf r\,d\mathbf r'\, n(\mathbf r) C_2(|\mathbf r-\mathbf r'|) n(\mathbf r')$$9, $w$0, $w$1, $w$2, $w$3, MPFC or Model B dynamics.

Typical simulation parameters are:

• FCC: $w$4, $w$5
• BCC: $w$6, $w$7 with system sizes up to $w$8 density peaks under periodic boundary conditions.

Strain is imposed by a penalty-region ("quasi-liquid") approach, through an added term:

$w$9

where $u$0 is large in boundary slabs and $u$1 is a rigid translation of the crystal, allowing controlled movement to impose shear or uniaxial strain (Berry et al., 2014).

3. Modeling Defect Sources and Nucleation Phenomena

PFC3D simulates both conservative and non-conservative dislocation sources:

• Frank–Read (glide) sources: A rectangular $u$2 prismatic loop is constructed, with sessile and glissile segments. Under applied shear ($u$3), glissile partials bow out; pinch-off and repeated nucleation are observed in three dimensions. Additionally, in polycrystallites generated via Voronoi tessellation, grain boundaries act as dislocation sources under strain, yielding stress–strain responses consistent with MD simulations and experimental findings, including reverse Hall–Petch effects for small grains.
• Bardeen–Herring (climb) sources: Spherical inclusions of higher modulus are modeled via penalty fields. Uniaxial strain triggers climb-nucleated $u$4 loops—detachment occurs through vacancy diffusion. The resultant loop geometries depend on strain orientation, inclusion size $u$5, and lattice anisotropy, producing serpentine, linked, and dual-loop structures.

Spherical inclusions or precipitates are introduced by applying a uniform penalty proportional to $u$6 within a defined spherical region—no explicit "void" reconstruction is necessary (Berry et al., 2014).

4. Critical Strain, Complex Defect Transformations, and 3D Structures

The critical strain for loop nucleation from an inclusion in continuum elasticity is given by

$u$7

where $u$8 is the Burgers vector magnitude, $u$9 the Poisson ratio, $C_2(r)$0 the bulk modulus ratio, and $C_2(r)$1 the sphere radius. PFC3D results match this analytic threshold when a finite-size strain offset $C_2(r)$2 is included, with best-fit $C_2(r)$3 and $C_2(r)$4 (Berry et al., 2014).

A newly observed dipole-to-quadrupole transformation occurs for rectangular sources with short glissile and long sessile segments: under shear, sessile segments stretch, misalign from their glide plane, nucleate jogs, and dissociate onto adjacent $C_2(r)$5 planes. This process results in a quadrupole (two dipoles) emission and subsequently returns to the original loop geometry, depending on activation stresses and segment mobilities.

Stacking-fault tetrahedra (SFTs), formed via Frank triangular loop nucleation and relaxation, are captured in PFC3D using the Silcox–Hirsch mechanism. Interactions between SFTs and gliding $C_2(r)$6 dislocations reproduce molecular dynamics results, including Lomer–Cottrell stair-rod pinning, Orowan loop formation, cut-through events, and ledge formation on SFT faces (Berry et al., 2014).

5. Timescale, Computational Efficiency, and Trade-Offs

PFC3D dynamics operate naturally on diffusive timescales. By calibrating the PFC vacancy diffusivity $C_2(r)$7 to experimental reference values (e.g., copper at $C_2(r)$8), the model achieves effective time steps in the microsecond regime ($C_2(r)$9–$$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$0s), compared to nanoseconds for classical MD. Strain rates are $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$1–$$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$2 s$$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$3, more than $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$4-fold slower than in MD.

A system of $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$5 atoms over 3 ms simulation time requires approximately 57 hours on 48 CPUs. Molecular dynamics for an equivalent system and time would require over $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$6 CPU hours. The PFC approach does not explicitly resolve phonons or deliver quantitative elastic constants but allows the simulation of long-term dynamics and arbitrary defect morphologies including climb (Berry et al., 2014).

6. Simulation Best Practices and Interpretive Guidelines

Best practices for PFC3D simulations include:

• Tuning correlation kernel parameters to match the target lattice constant, elastic properties, and stacking-fault energies.
• Employing a spectral spatial solver with grid spacing $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$7 to resolve density peaks; time step $$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$8 must ensure stability ($$\hat C_2(k) = -r + 1 - B^x (1-\tilde k^2)^2 - H_0 \exp\left[-\frac{(k-k_0)^2}{2\alpha_0^2}\right]$$9–$\tilde k = k/(2\pi)$0).
• Loading via penalty-region pinning, ensuring clean stress/strain fields under periodic boundary conditions.
• Seeding defects (dislocation loops, grain boundaries, inclusions) by direct density modifications or local penalty fields, before external driving.
• Monitoring the PFC vacancy diffusion coefficient to calibrate effective, physically meaningful time units.
• Interpreting results in the context of PFC’s coarse-grained atomistic nature: core geometries of close-packed dislocations are faithfully represented, but systems with complex chemistries or highly directional bonding may require extensions of the kernel.
• For simulations exceeding millisecond timescales, using MPI-parallelized FFTs can efficiently exploit the spectral method’s strengths (Berry et al., 2014).

PFC3D enables unified, atomistically-grounded studies of dislocation dynamics (glide and climb), defect nucleation at interfaces and inclusions, grain boundary emission, stacking fault–dislocation interactions, and the evolution of complex 3D microstructures on diffusive, experimentally relevant timescales. Its computational efficiency and flexibility make it especially suitable for investigating crystalline plasticity scenarios inaccessible to standard MD (Berry et al., 2014).