Coarse-grained molecular dynamics: Nonlinear finite elements and finite temperature

Coarse-grained molecular dynamics (CGMD) is a technique developed as a concurrent multiscale model that couples conventional molecular dynamics (MD) to a more coarse-grained description of the periphery. The coarse-grained regions are modeled on a mesh in a formulation that generalizes conventional finite element modeling (FEM) of continuum elasticity. CGMD is derived solely from the MD model, however, and has no continuum parameters. As a result, it provides a coupling that is smooth and provides control of errors that arise at the coupling between the atomistic and coarse-grained regions. In this paper, we elaborate on the formulation of CGMD, describing in detail how CGMD is applied to anharmonic solids and finite-temperature simulations. As tests of CGMD, we present in detail the calculation of the phonon spectra for solid argon and tantalum in three dimensions, demonstrating how CGMD provides a better description of the elastic waves than that provided by FEM. We also present elastic wave scattering calculations that show the elastic wave scattering is more benign in CGMD than FEM. We also discuss the dependence of scattering on the properties of the mesh. We introduce a rigid approximation to CGMD that eliminates internal relaxation, similar to the quasicontinuum technique, and compare it to the full CGMD.

Figure 1

(Color online) Schematic of a coarse-grained simulation of a NEMS silicon microresonator (Refs. 4, 5, 6). The coarse-grained region comprises most of the volume, but the molecular dynamics region contains most of the simulated degrees of freedom. Each sphere shown in the MD region represents an atom. Note that the CG mesh is refined to the atomic scale where it joins with the MD lattice.

Figure 2

Feynman diagrams of connected graphs contributing to the anharmonic CGMD internal energy. The thick external legs represent factors of $H$, which in the usual terminology correspond to dressed external fields. Physically, they are contributions from the CG fields with some zero-temperature relaxation of the internal degrees of freedom. The dashed internal lines represent factors of $G$ that arise from the fluctuations of the internal degrees of freedom. These fluctuations are temperature dependent, and each factor of $G$ comes with a factor of ${\beta }^{-1}(=kT)$.

Figure 3

The leading anharmonic contributions to the partition function affecting the thermal expansion in CGMD.

Figure 4

The anharmonic contributions to the partition function affecting the temperature dependence of the adiabatic CGMD stiffness. The Green’s function $G_1$ is represented by a dashed line with a large filled circle. It is defined in the upper box.

Figure 5

A comparison of the acoustic wave spectra for CGMD, R-CGMD, and MD with frequency vs wave number plotted. The units are $k_0=\pi ∕\left(N_{\text{per}}a\right)$ and ${\omega }_0=k_{0}a\sqrt{K∕m}$. The MD spectrum, corresponding to a simple 1D ball and spring model, is the ideal case. The CGMD and R-CGMD spectra are computed on a regular mesh with $N_{\text{per}}=32$ atoms per cell. Both full CGMD and its rigid approximation, R-CGMD, are in good agreement with the MD spectrum for $k\sim 0$, i.e., at long wavelengths. At short wavelengths near the zone boundary, the CGMD spectrum is more accurate than the R-CGMD spectrum, a property that we attribute to the relaxation effects accounted for by CGMD but eliminated in the rigid approximation.

Figure 6

A comparison of the error in the acoustic wave frequency at the zone boundary $[k=\pi ∕\left(N_{\text{per}}a\right)]$ for CGMD and R-CGMD as a function of the level of coarse graining, as expressed by $N_{\text{per}}$. At $N_{\text{per}}=1$, there is no error in either frequency. At all other values of $N_{\text{per}}$, the error in the CGMD frequency is less than that of R-CGMD, asymptotically going to 0.66% and 10.3%, respectively.

Figure 7

The phonon spectra for solid argon in the atomic limit are shown as plots of wave frequencies vs wave number along various high symmetry lines through the Brillouin zone (Ref. 3). The high symmetry points are labeled according to the standard convention for an fcc lattice (Refs. 62, 63); for example, $\Gamma$ is at the center of the zone, $\mathbf{k}=0$ and $X=(1,0,0)\pi ∕a$ where $a$ is the fcc lattice constant. The nine curves represent the three branches of the spectra for CGMD and two versions of FEM, one using the conventional distributed mass matrix and one using a diagonal lumped mass matrix. In the atomic limit, CGMD reproduces the Lennard-Jones MD spectra exactly, whereas the FEM spectra show significant error, especially with the distributed mass matrix.

Figure 8

The phonon spectra for solid argon on a mesh with slight coarse graining are shown as plots of wave frequencies vs wave number along various high symmetry lines through the Brillouin zone (cf. the caption to Fig. 7). The 12 curves represent the three branches of the spectra for MD, CGMD, and two versions of FEM, one using the conventional distributed mass matrix and one using a diagonal lumped mass matrix. Each cell of the CG mesh contains eight atoms. For this first level of coarse graining the CGMD spectra is in better agreement with the MD spectra than either of the FEM spectra. Of the two FEM spectra, the lumped mass spectrum is somewhat better.

Figure 9

The phonon spectra for solid argon on a mesh approaching the continuum limit (32 768 atoms per cell) are shown as plots of wave frequencies vs wave number along various high symmetry lines through the Brillouin zone (Ref. 3; cf. the caption to Fig. 7). The 12 curves represent the three branches of the spectra for MD, CGMD, and two versions of FEM, one using the conventional distributed mass matrix and one using a diagonal lumped mass matrix. With this significant level of coarse graining, the CGMD spectra is again in better agreement with the MD spectra than either of the FEM spectra. Of the two FEM spectra, now the distributed mass spectrum is better.

Figure 10

The room temperature phonon spectra for solid tantalum on a mesh with no coarse graining are shown as plots of wave frequencies vs wave number along various high symmetry lines through the Brillouin zone (for comparison to the solid argon case shown in Fig. 7). As in the case of argon, the CGMD phonon spectrum agrees exactly with the MD spectrum when the mesh is refined to the atomic limit. It is common practice to leave the gap between the second $N$ and $P$ points since that part of the spectrum is already represented to the left.

Figure 11

The room temperature phonon spectra for solid tantalum on a mesh with eight atoms per mesh cell are shown as plots of wave frequencies vs wave number along various high symmetry lines through the Brillouin zone (cf. Fig. 8). As in the case of argon, the CGMD phonon spectrum agrees very well with the MD spectrum in the limit of long wavelengths (near the $\Gamma$ point), and it agrees reasonably well throughout the coarse-grained Brillouin zone.

Figure 12

A comparison of the reflection of elastic waves from a CG region in three cases: CGMD and two varieties of FEM. Note that the reflection coefficient is plotted on a logarithmic scale. A similar graph plotted on a linear scale is shown in Ref. 17. The dashed line marks the natural cutoff $[k_0=\pi ∕\left(N_{\max }a\right)]$, where $N_{\max }$ is the number of atoms in the largest cells. The bumps in the curves are scattering resonances. Note that at long wavelengths, CGMD offers significantly suppressed scattering.

Figure 13

A plot of the mesh used for the scattering calculations. The atomic scale mesh, partially visible at the ends, extends infinitely to the right and left. The mesh is commensurate with the underlying atomic lattice, and the largest cell size is $N_{\max }=20$ atoms.

Figure 14

A comparison of the transmission of elastic waves through a CG region in three cases: CGMD and two varieties of FEM. The dashed line marks the natural cutoff $[k_0=\pi ∕\left(N_{\max }a\right)]$. Note that the bumps evident in the plot of the reflection coefficient are absent from the transmission coefficient.

Figure 15

A comparison of the reflection of elastic waves from a CG region for a mesh with smoothly varying cell size and one with an abrupt change in cell size, both computed in CGMD. The dashed line marks the natural cutoff $[k_0=\pi ∕\left(N_{\max }a\right)]$. Note that on the linear scale the resonances are not visible for the smooth mesh, but are quite pronounced for the abruptly changing mesh.