Multiscale Molecular Dynamics Simulations of Alkanes

We performed MD simulations of butane, pentane, and hexane in vacuum and TMC278 with AGBNP implicit solvation [30 (link)]. The numbers of dihedral angles excluding methyl groups rotations are one for butane, two for pentane, and three for hexane. The chemical structure of TMC278, which contains five rotatable bonds is shown in Fig. 1. The OPLS 2005 force field was used.[31 (link)] We employed the Nosé-Hoover thermostat,[32 –34 (link)] whose equations of motion are $${\dot{\mathit{r}}}_i=\frac{{\mathit{p}}_i}{m_i},$$
$${\dot{\mathit{p}}}_i={\mathit{F}}_i-\zeta {\mathit{p}}_i,$$
$$\dot{\zeta }=\frac{g{k}_B}{Q}(𝒯(t)-T_0),$$ where r_i, p_i, m_i, and F_i are coordinates, momentum, mass, and force of atom i, respectively. The variable ζ is a “viscosity” parameter for temperature control. The constant g is related to the number of atoms N as g = 3N − 6. The constant Q is the artificial “mass” for the thermostat. 𝒯 (t) and T₀ are the instantaneous temperature and the set temperature, respectively. We integrated the equations of motion in Eqs. (2)–(4) by the time-reversible algorithm by Martyna et al. [35 ]. We used the relation Q = gk_BT₀τ² for the Nosé-Hoover thermostat[36 ] with a relaxation time of τ = 10 fs. The MD time step was set to Δt = 0.5 fs.
Serial molecular dynamics simulations were conducted at T₀ = 300 K. We performed equilibration runs for 20 ns and then sampled conformations for 80 ns for the alkanes. For TMC 278, we performed an equilibration run for 50 ns and sampled conformations for 200 ns. Parallel replica exchange MD simulations were performed with four replicas at 300 K, 350 K, 400 K, and 450 K for the alkanes and with eight replicas at 300 K, 350 K, 400 K, 450 K, 500 K, 550 K, 600 K, and 650 K for TMC278. We performed REMD for 20 ns after 5 ns equilibration for the alkanes and for 25 ns after 6.25 ns equilibration runs for TMC278. The total sampling times are 80 ns for each alkane and 200 ns for TMC278. These sampling times are the same as for the corresponding serial MD simulations. Temperature exchanges were attempted 1 ps between adjacent temperatures and were accepted with probability $$w=\text{min}\{1,\text{exp}[-({\beta }_i-{\beta }_j)(E_i-E_j)]\},$$ where β_i = 1/k_BT_i and β_j = 1/k_BT_i are the inverse temperature before the exchange and E_i and E_j are the potential energies of replica i and j, respectively.
Dihedral angle distributions at a particular temperature can be calculated by simply binning the data from only the corresponding replica. We also calculated the dihedral angle distributions using the temperature WHAM method [29 (link), 37 ]: $$P_{\beta }(\varphi )={\displaystyle \sum _E\frac{{\displaystyle \sum _{i=1}^M{N}_i(E,\varphi )e^{-\beta E}}}{{\displaystyle \sum _{i=1}^M{n}_i{e}^{f_i-{\beta }_{i}E}}},}$$
$$e^{-f_i}={\displaystyle \sum _{\varphi }{P}_{{\beta }_i}(\varphi ),}$$ where P_β(ϕ) is the unnormalized angle distribution at inverse temperature β, M is the number of replicas, N_i(E, ϕ) is the histogram of potential energy, E, and dihedral angles ϕ = (ϕ₁, ϕ₂, …) at temperature T_i, and n_i is the total number of samples at temperature T_i. Eqs. (6) and (7) are solved iteratively [29 (link)].
Dihedral angle distributions for the alkanes were also obtained by numerical grid integration of the corresponding partition functions. Let us suppose a model of an alkane with fixed bond lengths. This model has only internal degrees of freedom of dihedral angles ϕ₀, ϕ₁, … ϕ_n+1 and bond angles θ₁, θ₂,… θ_n+1, where the number of dihedral angles n is n = 1 for butane, n = 2 for pentane, and n = 3 for hexane. ϕ₀ and ϕ_n+1 are the dihedral angles of the ends of the alkane chain, which include hydrogen atoms. Under these assumptions, the partition function Z of the molecule is given by $$Z={\displaystyle \int d{\varphi }_0\cdots d{\varphi }_{n+1}{\displaystyle \int d{\theta }_1\cdots d{\theta }_{n+1}{\displaystyle \prod _{i=1}^{n+1}\text{sin}{\theta }_i\text{exp}\{-\beta E(\varphi ,\theta )\},}}}$$ where the bold letters ϕ and θ denote the set of dihedral angles and bond angles, respectively. ${\displaystyle {\prod }_{i=1}^{n+1}\text{sin}{\theta }_i}$ is the Jacobian for the ϕ and θ internal coordinates. The dihedral angle distribution P(ϕ₁, … , ϕ_n) is calculated from the partition function Z by $$P({\varphi }_1,\cdots ,{\varphi }_n)=\frac{{\displaystyle \int d{\varphi }_{0}d{\varphi }_{n+1}{\displaystyle \int d{\theta }_1\cdots d{\theta }_{n+1}{\displaystyle \prod _{i=1}^{n+1}\text{sin}{\theta }_i\text{exp}\{-\beta E(\varphi ,\theta )\}}}}}{{\displaystyle \int d{\varphi }_0\cdots d{\varphi }_{n+1}{\displaystyle \int d{\theta }_1\cdots d{\theta }_{n+1}{\displaystyle \prod _{i=1}^{n+1}\text{sin}{\theta }_i\text{exp}\{-\beta E(\varphi ,\theta )\}}}}}.$$ The grid spacing of dihedral angles was set to 10°. Two grid points θ_eq and θ_eq + 4°, where θ_eq is the minimum energy bond angle, were employed for the bond angles integration (the potential energy of the cis conformation of butane is smallest for θ = θ_eq + 4°).