Optimization of partial multicanonical molecular dynamics simulations applied to an alanine dipeptide in explicit water solvent

Abstract

The partial multicanonical algorithm for molecular dynamics and Monte Carlo simulations samples a wide range of an important part of the potential energy. Although it is a strong technique for structure prediction of biomolecules, the choice of the partial potential energy has not been optimized. In order to find the best choice, partial multicanonical molecular dynamics simulations of an alanine dipeptide in explicit water solvent were performed with 15 trial choices for the partial potential energy. The best choice was found to be the sum of the electrostatic, Lennard-Jones, and torsion-angle potential energies between solute atoms. In this case, the partial multicanonical simulation sampled all of the local-minimum free-energy states of the P_II, C₅, α_R, α_P, α_L, and C^ax₇ states and visited these states most frequently. Furthermore, backbone dihedral angles ϕ and ψ rotated very well. It is also found that the most important term among these three terms is the electrostatic potential energy and that the Lennard-Jones term also helps the simulation to overcome the steric restrictions. On the other hand, multicanonical simulation sampled all of the six states, but visited these states fewer times. Conventional canonical simulation sampled only four of the six states: The P_II, C₅, α_R, and α_P states.

---

I. Introduction

Biomolecules like proteins and peptides have complicated free-energy surfaces with many local minima. An important problem is to predict the structures of these biomolecules theoretically or computationally. Conventional molecular dynamics (MD) and Monte Carlo (MC) simulations in physical ensembles, such as the canonical ensemble,^1–4 tend to get trapped in these local-minimum states and give inaccurate conformational results. Generalized-ensemble algorithms^5–9 are powerful techniques to avoid this difficulty. One of the most well-known generalized-ensemble algorithms is the multicanonical ensemble algorithm.^10–13 A non-Boltzmann weight factor is employed in this ensemble so that a free one-dimensional random walk can be realized in the potential-energy space. Thus, a simulation with this algorithm can escape from the free-energy-minimum states and sample a wide range of the conformational space. Because of this advantage, the multicanonical algorithm has frequently been applied to a biomolecule, which has a free-energy surface with many local minima. Modification of the multicanonical algorithm has also been proposed.^14–33

The non-Boltzmann weight factor in the multicanonical algorithm is not a priori known and has to be determined by a preliminary simulation.^34–36 As the system size increases, the distribution of the total potential energy P(E) gets narrower in proportion to , where N is the number of atoms. Thus, the determination of the multicanonical weight factor to give a flat distribution becomes difficult in a large system. New techniques are desired to reduce the effort for the weight-factor determination with keeping the high sampling efficiency. The multicanonical algorithm covers a wide range of the total potential energy. The total potential energy of a biomolecular system includes bond-length, bond-angle, torsion-angle, improper torsion-angle, Lennard-Jones, and electrostatic energy terms for most of the force fields.^37–41 However, it is redundant to sample a wide range of all these terms.

The author proposed recently the partial multicanonical ensemble algorithm for MD and MC simulations,²⁵ in which only an important part of potential-energy terms, instead of the total potential energy, are sampled widely. The partial potential energy terms can, in principle, be chosen arbitrarily. The purpose of this article is to optimize the choice of the partial potential energy terms. For this purpose, the partial multicanonical MD simulations of an alanine dipeptide are performed in explicit water and many choices of the partial potential energy are tested. Conventional canonical and multicanonical MD simulations are performed, too. Their sampling efficiencies are compared with the partial multicanonical MD simulations. The author also discusses which potential energy term contributes most to sample many conformations. The alanine dipeptide has been studied theoretically extensively as a model biomolecule because of its simplicity (see, for example, ref. 23–25 and 42–49).

In section II the partial multicanonical algorithm is explained. The computational details of the alanine dipeptide simulations are described in section III. The results and discussion are presented in section IV. Section V is devoted to conclusions.

II. Methods

In the canonical ensemble, the distribution P_NVT(E) of the total potential energy E is given by

P_NVT(E) = n(E)e^−β₀E, (1)

where n(E) is the density of states, β₀ is the inverse of the product of the Boltzmann constant k_B and the absolute temperature T₀, at which simulations are performed. This ensemble has a bell-shaped distribution in the potential-energy space.

In the multicanonical ensemble,^10–13 every state is sampled with a weight factor W_muca(E) ≡ exp{−β₀E_muca(E)} so that a uniform distribution of E may be obtained:

P_muca(E) = n(E) W_muca(E) = n(E) exp{−β₀E_muca(E)} = constant, (2)

where E_muca(E) and W_muca(E) are referred to as the multicanonical energy and the multicanonical weight factor, respectively. The difference between E_muca and E is written as δE(E): E_muca(E) ≡ E + δE(E). The difference δE(E) is additional potential energy which is introduced to obtain a flat distribution of E in eqn (2). The case of δE(E) = 0 gives the regular canonical ensemble. The multicanonical simulation covers a wide range of E and escapes from local-minimum free-energy states. The expectation value of a physical quantity A at temperature T is obtained by reweighting techniques⁵⁰ fromwhere 〈⋯〉_muca is the multicanonical ensemble average, which is defined as 〈⋯〉_muca ≡ ∫dr⋯e^{−β₀{E(r)+δE(E(r))}}/∫dre^{−β₀{E(r)+δE(E(r))}}. The vector r is the set of coordinates r ≡ {r₁,…,r_N}. Because of the random walks in the potential-energy space, we can calculate physical quantities in a wide range of T.

Although the multicanonical algorithm is a powerful technique, the multicanonical weight factor W_muca(E) or the multicanonical energy E_muca(E) is not a priori known and has to be determined in advance of the simulation for data collection, for example, by the iterations of short simulations.^34–36 The determination of the multicanonical weight factor is time consuming for a larger system. In order to alleviate the effort to determine the weight factor, the total potential energy is divided into two terms in the partial multicanonical algorithm:

E = E₁ + E₂, (4)

where E₁ is the part of the total potential energy which needs to be sampled for a wide range of conformational space. The partial potential energy E₁ can be composed of any of the components of E, in principle. The second term E₂ is the rest of the potential energy. The partial multicanonical algorithm ensures that the wide range of E₁ is efficiently sampled. The partial multicanonical simulation is performed with a weight factor W_pmuca(E₁, E₂) ≡ exp{−β₀E_pmuca(E₁, E₂)} so that a uniform distribution of E₁ may be obtained:where E_pmuca(E₁, E₂) and W_pmuca(E₁, E₂) are referred to as the partial multicanonical energy and the partial multicanonical weight factor, respectively. This equation means that P_pmuca(E₁, E₂) depends on E₂ but is independent of E₁. The difference between E_pmuca and E depends only on E₁ and is written as δE₁(E₁):

E_pmuca(E₁, E₂) ≡ E₁ + E₂ + δE₁(E₁) = E + δE₁(E₁). (6)

The difference δE₁(E₁) is characteristic of the partial multicanonical simulation. If we define here a distribution P_pmuca(E₁) of E₁ and the density of states n(E₁) of E₁ as

P_pmuca(E₁) ≡ ∫dE₂P_pmuca(E₁, E₂), (7)

n(E₁) ≡ ∫dE₂n(E₁, E₂)exp(−β₀E₂), (8)

Eqn (5) is rewritten asThe expectation value of a physical quantity A is obtained by the reweighting techniques⁵⁰ fromwhere 〈⋯〉_pmuca is the partial multicanonical ensemble average, which is defined as 〈⋯〉_pmuca ≡ ∫dr⋯e^{−β₀{E(r)+δE₁(E₁(r))}}/∫dre^{−β₀{E(r)+δE₁(E₁(r))}}. Note that 〈A〉_NVT can be calculated correctly at a temperature T only in the vicinity of T₀, although 〈A〉_NVT can be calculated in a wide range of T in the multicanonical algorithm. It is because a limited range of E₂ is sampled with the normal Boltzmann weight factor exp(−β₀E₂) in the partial multicanonical algorithm.

The scheme of the partial multicanonical algorithm includes the canonical and the multicanonical algorithms. In the case of E₁ = 0 and E₂ = E, this algorithm becomes the regular canonical algorithm as a limiting case. This is because eqn (5) becomes the canonical distribution as P_pmuca(E) = n(E) exp(−β₀E). In the case of E₁ = E and E₂ = 0, the multicanonical algorithm is obtained because eqn (5) becomes the multicanonical distribution as P_pmuca(E) = n(E) exp{−β₀(E + δE)} = constant. Furthermore, if we take E₁ as the sum of the potential energy between solute atoms E^slt–slt and that of the solute–solvent interaction E^slt–svn (i.e., E₁ = E^slt–slt + E^slt–svn) and E₂ as the potential energy between solvent atoms E^svn–svn (i.e., E₂ = E^svn–svn), the partial multicanonical algorithm gives the selectively enhanced multicanonical algorithm.¹⁷ The partial multicanonical algorithm can be applied not only to an explicit solvent model but also to an implicit solvent model, in a similar way to the combination of the multicanonical algorithm and the RISM theory.⁵¹

A constant-temperature MD method is employed to perform a partial multicanonical MD simulation. One of the representative constant-temperature methods is the Nosé thermostat.^1,2 The partial multicanonical Hamiltonian based on the Nosé thermostat is given by

where is the conjugate momentum for r_i. The real momentum p_i is related to the virtual momentum by . The variable s is the additional degree of freedom for the Nosé thermostat. The variable P_s is the conjugate momentum for s. The constant m_i is the mass of atom i and Q is the artificial “mass” related to s. The real time interval Δt is related to the virtual time interval Δt′ by Δt = Δt′/s. In the case of a system consisting of N atoms, g equals 3N if the time development is performed in the real time t, or g equals 3N + 1 if the time development is performed in virtual time t′.

The equations of motion in the real time development are derived from the Hamiltonian in eqn (11). In the Nosé-Hoover formalism,^1–3 a variable ζ ≡ ṡ/s is introduced and the equations of motion are given by

where the dot (˙) stands for the real time derivative d/dt and forces are defined byIn the case of a rigid-body molecule system, the equations of motion in the partial multicanonical ensemble are given bywhere subscript α is an index of the rigid-body molecules, r_α is the coordinate of the center of mass of molecule α, p_α is the momentum of molecule α, m_α is the mass of molecule α, and q_α is a quaternion of molecule α, which indicates the orientation of the rigid-body molecule. Here, the quaternion q = (q₀, q₁, q₂, q₃)^T is related to the Euler angle (ϕ, θ, ψ) as follows:The elements of the matrix are given byω and ω⁽⁴⁾ are three-dimensional and four-dimensional angular velocity, respectively, which are given by

ω = (ω₁, ω₂, ω₃)^T, (27)

ω⁽⁴⁾ = (0, ω₁, ω₂, ω₃)^T, (28)

where ω₁, ω₂, and ω₃ are the angular velocities along each of the corresponding principal axes. The vector F_α is the force acting on molecule α, which is related to the force F_i acting on atom i in molecule α byThe torque N_α is calculated bywhere r_i and F_i are the coordinate and force of atom i, respectively, in a rigid-body-fixed coordinate system for molecule α. The matrix I_α is a 3 × 3 diagonal matrix consisting of the principle moments of inertia of molecule α.

The equations of motion (13), (18), and (20) are obtained from the regular canonical-ensemble equations of motion by modifying the force from

F_i = F_1i + F_2i (31)

toThe partial multicanonical MD formalism has been described here based on the Nosé-Hoover thermostat. The partial multicanonical MD simulation can also be performed in a similar way based on the Gaussian thermostat by modifying the force from eqn (31) to (32).

Instead of MD simulations, it is also possible to perform the partial multicanonical simulations using the MC method. In the partial multicanonical MC simulations, all atoms are updated from r_i to as trial moves by uniform random numbers. Let us suppose that the partial multicanonical energy is consequently changed from E_pmuca(r) to by these trial moves. The Metropolis criterion⁴ is employed here with the weight factor of exp(−β₀E_pmuca) instead of the Boltzmann weight factor exp(−β₀E) in the canonical ensemble. These trial moves are accepted with the probability

The partial multicanonical probability distribution P_pmuca(E₁) can be obtained either by the MD simulation in eqn (12)–(21) or by the MC simulation in eqn (33).

III. Computational details

Molecular dynamics simulations were performed for a system consisting of one alanine dipeptide molecule ((S)-2-(acetylamino)-N-methylpropanamide) and 67 water molecules. The N-terminus and the C-terminus of the alanine residue were blocked by the acetyl group and the N-metyl group, respectively. Enough water molecules are used so that the alanine dipeptide molecule was always held perfectly within the simulation box. A cubic unit cell was employed with periodic boundary conditions. The side length was set to be 13.39 Å. Thus, the volume V of the unit cell was 2400 ų. The initial condition of the backbone dihedral angles ϕ and ψ were ϕ = ψ = 180°, as shown in Fig. 1. The AMBER parm96 force field³⁸ was used for the alanine dipeptide molecule and the TIP3P rigid-body model⁵³ was employed for the water molecules. The electrostatic potential energy was calculated using the Ewald method. The Lennard-Jones interaction was calculated for all pairs of the atoms within the minimum image convention instead of introducing the spherical potential cutoff.

Fig. 1 Initial conformation of alanine dipeptide for the MD simulations. The dihedral angles are ϕ = ψ = 180°. The figure was created with RasMol.⁵²

The combination⁵⁴ of the Nosé-Hoover thermostat^1–3 and the symplectic quaternion scheme for the rigid-body water molecules⁵⁵ was used. The time step was taken to be Δt = 0.5 fs. The temperature was set to T₀ = 300 K.

The total potential energy is the sum of the potential energy between solute atoms E^slt–slt, that of the solute–solvent interaction E^slt–svn, and that between solvent atoms E^svn–svn as follows:

E = E^slt–slt + E^slt–svn + E^svn–svn. (34)

Each term for the AMBER force field^37,38 is given by

E^slt–svn = E^slt–svn_elec + E^slt–svn_LJ, (36)

E^svn–svn = E^svn–svn_elec + E^svn–svn_LJ, (37)

where the subscripts bond, angle, torsion, imp, elec, and LJ stand for the bond-length, bond-angle, torsion-angle, improper torsion-angle, electrostatic, and Lennard-Jones potential energy, respectively. Although these terms may vary for different force fields, the partial multicanonical algorithm can be applied in essentially the same way.

The purpose of this article is to find the best choice of the partial potential energy E₁. For this purpose, the following 15 cases of the E₁ choices were tested.

case 1: E₁ = 0, (38)

case 2: E₁ = E^slt–slt + E^slt–svn + E^svn–svn = E, (39)

case 3: E₁ = E^slt–slt + E^slt–svn, (40)

case 4: E₁ = E^slt–slt, (41)

case 5: E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, (42)

case 6: E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp, (43)

case 7: E₁ = E^slt–slt_tor + E^slt–slt_elec, (44)

case 8: E₁ = E^slt–slt_tor + E^slt–slt_LJ, (45)

case 9: E₁ = E^slt–slt_elec + E^slt–slt_LJ, (46)

case 10: E₁ = E^slt–slt_tor, (47)

case 11: E₁ = E^slt–slt_elec, (48)

case 12: E₁ = E^slt–slt_LJ, (49)

case 13: E₁ = E^slt–slt_bond, (50)

case 14: E₁ = E^slt–slt_angle, (51)

case 15: E₁ = E^slt–slt_imp. (52)

Cases 1, 2, and 3 correspond to the conventional canonical MD simulation, the multicanonical MD simulation, and the selectively enhanced multicanonical MD simulation, respectively. In case 4, E₁ includes all of the solute–solute interactions. Case 5 is the same choice as in the preceding article (see ref. 25). At that time the author expected this choice may work well because of the following reasons: even if we sample wide ranges of bond-length, bond-angle, and improper torsion-angle energy in E^slt–slt , the conformation of the solute molecule is not expected to be changed very much. Sampling a wide range of solvent energy E^svn–svn would not change the solute structure very much, either. Sampling a wide range of the solute–solvent potential energy E^slt–svn may sample a wide range of the interface area between the solute molecule and solvent molecules. The solute molecule of the alanine dipeptide, however, is small and this effect may not necessarily have to be taken into account in this system. In case 6, E₁ includes the solute–solute energy terms which were not expected to be important in the preceding article.²⁵ Cases 7–9 include two of the three energy terms in case 5. Cases 10–15 include only one of the solute–solute energy terms.

To obtain an optimal partial multicanonical weight factor W_pmuca(E₁, E₂), both the iterative method^34–36 and the Wang-Landau techniques⁵⁶ were employed. This is because this combination is more efficient to determine W_pmuca(E₁, E₂) than employing only the iterative method. This technique updates W_pmuca(E₁, E₂) at every MD step n as W_pmuca(E₁, E₂, n) ≡ exp[−{E₁ + E₂ + δE₁(E₁, n)}/k_BT₀] to lower energy barriers. The difference δE₁(E₁, n) at time step n in the i-th iteration was updated using a histogram of E₁ and is given by

where P⁽ⁱ⁾(E₁, n) is the histogram of E₁ accumulated until one time step before n, that is, until time step (n − 1). Here, the time step in the i-th iteration takes the values n = 1,2,…,N⁽ⁱ⁾_t. For the first iteration (i = 1), the partial multicanonical energy starts from the normal potential energy with

δE⁽⁰⁾₁(E₁) = 0. (54)

In the case of i ≥ 2, δE⁽ⁱ⁻¹⁾₁(E₁) was given by the histogram P⁽ⁱ⁻¹⁾(E₁) of the partial potential energy E₁ accumulated until the last time step in the (i − 1)-th iteration as follows:where log P⁽ⁱ⁻¹⁾(E₁) was shifted by log P⁽ⁱ⁻¹⁾(E^min₁) for the partial potential energy larger than a threshold E^min₁(E₁ ≥ E^min₁) to make the potential-energy distribution wider towards the larger value of E₁ to overcome free-energy barriers. The threshold E^min₁ was determined as the minimum value of E₁ in the canonical MD simulation (MD simulation of the case 1). This procedure was iterated 5 times (i = 1–5). In the first iteration the MD simulation was performed for 100 ps (N⁽¹⁾_t = 200 000). At later iterations, the MD simulations we carried out for longer times, e.g., for 500 ps (N⁽⁵⁾_t = 1 000 000) in the fifth iteration. An optimal weight factor, W_pmuca(E₁, E₂), was obtained after these preliminary simulations for 1.5 ns in total.

A partial multicanonical MD simulation was performed for an equilibration process for 20 ps before sampling. A long partial multicanonical MD simulation was then performed for 10 ns for data collection. This simulation length is longer than the preceding partial multicanonical article,²⁵ in which the data were collected for 6 ns.

For convenience, E₁ was discretized into bins and histograms used to calculate δE₁(E₁). The bin size was set as ΔE₁ = 5 kcal mol⁻¹ for the multicanonical (case 2) and selectively enhanced multicanonical (case 3) simulation, ΔE₁ = 2 kcal mol⁻¹ for the cases 4, 6, and 13 simulation, and ΔE₁ = 1 kcal mol⁻¹ for the other partial multicanonical simulations. In order to calculate the derivative of δE₁(E₁), that is, ∂δE₁(E₁)/∂E₁, it is necessary to interpolate the histogram. The histogram was interpolated by a third-order polynomial in E₁ to make the derivative continuous at the boundaries between two adjacent bins (see ref. 22). Deviation of the value of the Hamiltonian from its initial value was suppressed by this polynomial.

IV. Results and discussion

The time series of the total potential energy E are shown in Fig. 2. Only four results are shown here for simplicity, which were obtained by the conventional canonical (E₁ = 0), multicanonical (E₁ = E), selectively enhanced multicanonical (E₁ = E^slt–slt + E^slt–svn), and partial multicanonical (E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ) MD simulations. The potential energy E fluctuated in a narrow range in the canonical MD simulation, as shown in Fig. 2(a). The fluctuation ranges in the selectively enhanced multicanonical and partial multicanonical MD simulations were almost the same as that in the canonical simulation, as shown in Fig. 2(c) and (d). The energy fluctuation range of the selectively enhanced multicanonical and that of the partial multicanonical was 1.33 times wider and 1.10 times wider, respectively, than that of the canonical MD simulations. On the other hand, Fig. 2(b) shows that the multicanonical MD simulation covered a wide range of E. The energy range was 1.87 times wider than that of the canonical MD simulations.

Fig. 2 The time series of total potential energy E obtained by (a) the canonical (E₁ = 0), (b) the multicanonical (MUCA, E₁ = E), (c) the selectively enhanced multicanonical (SEMUCA, E₁ = E^slt–slt + E^slt–svn), and (d) the partial multicanonical (PMUCA, E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ) MD simulation.

Fig. 3 shows the time series of E^slt–slt + E^slt–svn. As expected, the selectively enhanced multicanonical MD simulation, of which E₁ was set as E₁ = E^slt–slt + E^slt–svn , performed a random walk in a wide range of E^slt–slt + E^slt–svn, as shown in Fig. 3(c). The fluctuation range was 2.06 times wider in the selectively enhanced multicanonical MD simulation than that in the canonical MD simulation. On the other hand, this energy fluctuated in a narrow range in the canonical, multicanonical, and partial multicanonical simulations. The energy ranges in the multicanonical and partial multicanonical simulations were only 1.25 times wider and 1.26 times wider, respectively, than that in the canonical MD simulation.

Fig. 3 The time series of partial potential energy E^slt–slt + E^slt–svn. See the caption of Fig. 2 for further details.

Fig. 4 shows the time series of E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ. The partial multicanonical MD simulation in the case of E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ performed a random walk in a wide range of this energy, as shown in Fig. 4(d). This fluctuation range was 2.65 times wider than that in the canonical MD simulation. On the other hand, this energy fluctuated in a narrow range in the canonical, multicanonical, and selectively enhanced multicanonical simulations. The ranges in the multicanonical and selectively enhanced multicanonical simulations were only 1.19 times wider and 1.39 times wider, respectively, than the canonical MD simulation.

Fig. 4 Time series of the partial potential energy E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ. See the caption of Fig. 2 for further details.

The probability distributions of the potential energy are shown in Fig. 5. The distribution of the total potential energy E in the selectively enhanced multicanonical and partial multicanonical simulation was slightly wider than that in the canonical simulation and shifted toward a higher value. The distribution of E in the multicanonical simulation was much wider than those in the other two simulations. Similar situations can be seen in E^slt–slt + E^slt–svn and E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ. The selectively enhanced multicanonical algorithm sampled a wide range of E^slt–slt + E^slt–svn, whereas its distributions in the multicanonical and partial multicanonical simulations were almost the same as that in the canonical simulation. The partial multicanonical simulation covered a wider range of E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ than the canonical, multicanonical and selectively enhanced multicanonical simulations. Fig. 2–5 show that these generalized-ensemble MD simulations actually sampled a wide range of the potential energy which and only which are pre-designated to be sampled widely.

Fig. 5 Logarithms of the probability distributions of (a) total potential energy E, (b) partial potential energy E^slt–slt + E^slt–svn, and (c) partial potential energy E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ. The dotted lines were obtained by the canonical (E₁ = 0), the short-dashed lines were obtained by the multicanonical (E₁ = E), the long-dashed lines were obtained by the selectively enhanced multicanonical (E₁ = E^slt–slt + E^slt–svn), and the solid lines were obtained by the partial multicanonical (E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ) MD simulation.

The time series of the dihedral angle ϕ are shown in Fig. 6. The canonical MD simulation did not make a complete rotation, as shown in Fig. 6(a). It covered only 65% of the whole range of ϕ. There is potential energy barriers around ϕ = 0° and ϕ = 120° which cannot be overcome by the canonical simulation (these energy barriers correspond to the boundaries between the local-minimum free-energy states defined in Table 2 below). It means that the canonical MD simulation was trapped in the local-minimum free-energy states.

Fig. 6 The time series of the backbone dihedral angle ϕ obtained by (a) the canonical (E₁ = 0), (b) the multicanonical (MUCA, E₁ = E), (c) the selectively enhanced multicanonical (SEMUCA, E₁ = E^slt–slt + E^slt–svn), the partial multicanonical (PMUCA) MD simulations in the cases of (d) E₁ = E^slt–slt, (e) E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, (f) E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp, (g) E₁ = E^slt–slt_tor + E^slt–slt_elec, (h) E₁ = E^slt–slt_tor + E^slt–slt_LJ, (i) E₁ = E^slt–slt_elec + E^slt–slt_LJ, (j) E₁ = E^slt–slt_tor, (k) E₁ = E^slt–slt_elec, (l) E₁ = E^slt–slt_LJ, (m) E₁ = E^slt–slt_bond, (n) E₁ = E^slt–slt_angle, and (o) E₁ = E^slt–slt_imp.

As well as the canonical MD simulation, the MD simulations did not overcome these energy barriers at ϕ = 0° and ϕ = 120°, if E₁ was selected as E₁ = E^slt–slt (Fig. 6(d)), E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp (Fig. 6(f)), or E₁ = E^slt–slt_angle (Fig. 6(n)). It means that these choices of E₁ do not contribute to enhance the sampling efficiency.

The multicanonical MD simulation did not make a complete rotation, either. It overcame the potential energy barrier around ϕ = 0°, but did not overcome that around ϕ = 120° and sampled 91% of the whole range of ϕ, as shown in Fig. 6(b). The multicanonical MD simulation surely have higher sampling efficiency than the canonical MD simulation. Note that these simulation results do not mean that the multicanonical algorithm cannot sample the whole range of ϕ. The results mean that longer preliminary simulations are required to determine a better multicanonical weight factor W_muca(E). If a better weight factor W_muca(E) were obtained with the longer preliminary simulations, the whole range of ϕ could be sampled. In this article, the multicanonical and the partial multicanonical algorithms are compared in the same condition of the effort to determine the weight factors for 1.5 ns.

The partial multicanonical MD simulations with E₁ = E^slt–slt_bond (Fig. 6(m)) and E₁ = E^slt–slt_imp (Fig. 6(o)) also sampled the range of 0° < ϕ < 120°, but only once. The sampling efficiency of these choices might be slightly better than the canonical MD simulation, but worse than the multicanonical MD simulation.

The whole range of ϕ was sampled in the cases of E₁ = E^slt–slt + E^slt–svn (Fig. 6(c)), E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ (Fig. 6(e)), E₁ = E^slt–slt_tor + E^slt–slt_elec (Fig. 6(g)), E₁ = E^slt–slt_tor + E^slt–slt_LJ (Fig. 6(h)), E₁ = E^slt–slt_elec + E^slt–slt_LJ (Fig. 6(i)), E₁ = E^slt–slt_tor (Fig. 6(j)), E₁ = E^slt–slt_elec (Fig. 6(k)), and E₁ = E^slt–slt_LJ (Fig. 6(l)). These MD simulations overcame both energy barriers and covered the whole range of ϕ. These results mean that the partial multicanonical algorithm with appropriate choice of E₁ shows higher sampling efficiency than the multicanonical algorithm. Table 1 shows how many times ϕ underwent a complete rotation (covering 360°) in these MD simulations for 10 ns. Here, one turn of rotation is defined when ϕ undergoes all values between −180 and 180° and count another turn when ϕ undergoes all values after that. The dihedral angle ϕ rotated particularly well if E₁ was set as E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, E₁ = E^slt–slt_tor + E^slt–slt_LJ, E₁ = E^slt–slt_elec + E^slt–slt_LJ, or E₁ = E^slt–slt_LJ. In these choices of E₁, ϕ made seven turns or more. The Lennard-Jones potential energy is commonly included in these choices of E₁. It indicates that sampling a higher value of the Lennard-Jones potential energy helps the MD simulation to overcome steric restrictions when the solute atoms collide with one another.

Table 1 Numbers of the complete rotations of the backbone dihedral angles ϕ and ψ during the 10 ns MD simulations. The numbers in parentheses are the ratio of the dihedral angle range sampled by the MD simulation to the whole range of the dihedral angle (360°) if the whole range was not sampled

Choice of E1	ϕ	ψ
E 1 = 0	(0.65)	27
E 1 = E	(0.91)	54
E 1 = E^slt–slt + E^slt–svn	3	68
E 1 = E^slt–slt	(0.49)	(0.92)
E 1 = E^slt–slttor + E^slt–sltelec + E^slt–sltLJ	8	46
E 1 = E^slt–sltbond + E^slt–sltangle + E^slt–sltimp	(0.56)	9
E 1 = E^slt–slttor + E^slt–sltelec	2	6
E 1 = E^slt–slttor + E^slt–sltLJ	10	89
E 1 = E^slt–sltelec + E^slt–sltLJ	7	53
E 1 = E^slt–slttor	2	37
E 1 = E^slt–sltelec	1	31
E 1 = E^slt–sltLJ	9	57
E 1 = E^slt–sltbond	(0.82)	37
E 1 = E^slt–sltangle	(0.56)	18
E 1 = E^slt–sltimp	(0.82)	31

---

Fig. 7 shows the time series of the other backbone dihedral angle ψ. All MD simulations sampled the whole range of ψ except the case of E₁ = E^slt–slt. The dihedral angle ψ rotated more easily than ϕ even in the canonical MD simulation. Table 1 also shows how many times ψ realized a complete rotation. The following choices of E₁ made ψ rotate well (more than 40 turns): the multicanonical (E₁ = E), selectively enhanced multicanonical (E₁ = E^slt–slt + E^slt–svn), partial multicanonical with E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, E₁ = E^slt–slt_tor + E^slt–slt_LJ, E₁ = E^slt–slt_elec + E^slt–slt_LJ, and E₁ = E^slt–slt_LJ. These choices of E₁ also include the Lennard-Jones potential energy. This means again that sampling of a wide range of the Lennard-Jones energy helps the simulations to overcome the free-energy barriers.

Fig. 7 Time series of the backbone dihedral angle ψ. See the caption of Fig. 6 for further details.

Fig. 8 shows the potential of the mean force or the free-energy surface A(ϕ,ψ) as a function of ϕ and ψ at T = 300 K. Fig. 8 also shows the borders between the adjacent states. These borders are listed in Table 2, too. The potential of mean force A(ϕ, ψ) was calculated from the probability distributions P_NVT(ϕ, ψ) of ϕ and ψ in the canonical ensemble by

A(ϕ, ψ) = −k_BT log P_NVT(ϕ, ψ). (56)

In the canonical MD simulation, the distribution P_NVT(ϕ,ψ) is easily calculated. On the other hand, to obtain the distribution P_NVT(ϕ, ψ) from the multicanonical simulation and that from the partial multicanonical simulation, the reweighting techniques in eqn (3) and (10) were used, respectively. For example, P_NVT(ϕ, ψ) was calculated from the partial multicanonical MD simulation as follows: fist, an unnormalized histogram 〈N(ϕ, ψ)〉_NVT of ϕ and ψ in the canonical ensemble was given by

〈N(ϕ,ψ)〉_NVT = 〈N(ϕ,ψ; r)exp[β₀{E(r) + δE₁(E₁(r))}]exp{−βE(r)}〉_pmuca, (57)

where ϕ and ψ were discretized using the bin sizes of Δϕ = Δψ = 10° to calculate the histogram N(ϕ, ψ; r) andwhere ϕ(r) and ψ(r) stand for ϕ and ψ, respectively, for the coordinates r. Eqn (57) means that 〈N(ϕ, ψ)〉_NVT is obtained by summing up the values of exp[β₀{E(r) + δE₁(E₁(r))}]exp{−βE(r)} at (ϕ(r), ψ(r)). The probability distribution P_NVT(ϕ, ψ) was then obtained by normalizing 〈N(ϕ, ψ)〉:The peak position of P(ϕ, ψ) of each state obtained by the partial multicanonical MD simulation with E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ is listed in Table 3.

Fig. 8 Contour maps of the potential of mean force A(ϕ, ψ) as a function of the backbone dihedral angles ϕ and ψ obtained by the reweighting techniques at T = 300 K. See the caption of Fig. 6 for further details.

Table 2 Dihedral angle ranges of (ϕ, ψ) for each state. The state labels are taken from ref. 46

State	ϕ	ψ
PII	(−100°, 0°)	(30°, −120°)
C5	(120°, −100°)	(30°, −120°)
α R	(−100°, 0°)	(−120°, 30°)
α P	(120°, −100°)	(−120°, 30°)
α L	(0°, 120°)	(0°, 140°)
C^ax7	(0°, 120°)	(140°, 0°)

---

Table 3 Peak positions (ϕ, ψ) of the probability distribution P_NVT(ϕ, ψ) for six states at T = 300 K. The results were obtained by the reweighting techniques from the partial multicanonical MD simulation in the case of E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp

State	(ϕ, ψ)
PII	(−65°, 145°)
C5	(−145°, 155°)
α R	(−65°, −45°)
α P	(−145°, −55°)
α L	(45°, 65°)
C^ax7	(55°, −85°)

---

The volume under the surface of P(ϕ, ψ) around each peak corresponds to the population W of each state. In order to calculate W, the whole (ϕ, ψ) plane was divided into six states as shown in Fig. 8 and Table 2. For example, the population W_PII of the P_II state is calculated by the integral of P(ϕ, ψ) in the area in which ϕ and ψ take the P_II conformation:

where the integration area of (ϕ, ψ) means the range for the corresponding state in Table 2. The population of each state at T = 300 K calculated by the reweighting techniques from the results of the partial multicanonical MD simulation is listed in Table 4. The statistical uncertainties were estimated by the jackknife method⁵⁷ in which the production run was divided into ten segments. The populations of the P_II and C₅ states are higher than those of the α_R and α_P states. In the case of longer peptides or proteins, the α_R state corresponds to an α-helix structure and the P_II and C₅ states correspond to a β-strand structure. It is known that, in general, the AMBER parm96 force field tends to form a β-strand structure.⁵⁸ The potential of mean force A(ϕ, ψ) in Fig. 8 and Table 4 is consistent with this feature and P_NVT(ϕ, ψ) calculated by other methods such as the multibaric-multithermal MD algorithm.²⁴

Table 4 Population W of each state at T = 300 K obtained by the reweighting techniques from the results of the partial multicanonical MD simulation in the case of E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp

State	W
PII	0.409(18)
C5	0.470(21)
α R	0.067(7)
α P	0.051(7)
α L	0.00062(12)
C^ax7	0.0020(14)

---

Fig. 8 shows that the P_II, C₅, α_R, and α_P states were sampled in the canonical MD simulation (Fig. 8(a)) and the partial multicanonical MD simulations in the cases of E₁ = E^slt–slt (Fig. 8(d)), E₁ = E^slt–slt_bond + E^slt–slt_angle + E^slt–slt_imp (Fig. 8(f)), and E₁ = E^slt–slt_angle (Fig. 8(n)). Neither the α_L state nor the C^ax₇ state were sampled in these simulations because the range of ϕ between 0 and 120° was not sampled in Fig. 6. These choices of E₁ are not appropriate to enhance the sampling efficiency.

All of the six states of the P_II, C₅, α_R, α_P, α_L, and C^ax₇ states were sampled in the multicanonical (E₁ = E, Fig. 8(b)), selectively enhanced multicanonical (E₁ = E^slt–slt + E^slt–svn, Fig. 8(c)), and partial multicanonical MD simulations in the cases of E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ (Fig. 8(e)), E₁ = E^slt–slt_tor + E^slt–slt_elec (Fig. 8(g)), E₁ = E^slt–slt_tor + E^slt–slt_LJ (Fig. 8(h)), E₁ = E^slt–slt_elec + E^slt–slt_LJ (Fig. 8(i)), E₁ = E^slt–slt_tor (Fig. 8(j)), E₁ = E^slt–slt_elec (Fig. 8(k)), E₁ = E^slt–slt_LJ (Fig. 8(l)), E₁ = E^slt–slt_bond (Fig. 8(m)), and E₁ = E^slt–slt_imp (Fig. 8(o)).

Although these choices of E₁ sampled all states, not all of these choices were good. In order to judge which choice of E₁ works well, the author calculated how many times each simulation visits each state, as listed in Table 5. The number of visiting times is good information to find the best choice of E₁ as well as the number of rotations of ϕ and ψ. The simulation, however, enters into the neighboring state's area moving by a few degrees of ϕ or ψ, but soon comes back to the original state's area. For example, a small part of the foot of the P_II state peak extends over the area of the α_L state in the canonical ensemble, as shown in Fig. 8(a), but this simulation essentially did not sample the α_L state. To avoid counting such a situation as a visit, when the simulation comes into the area listed in Table 2 for more than 10 degrees, it counts for one visit. For example, when the simulation enters into the area of 10° < ϕ < 110° and 10° < ψ < 130°, it is regarded as one visit to the α_L state, although the α_L state area is originally 0° < ϕ < 120° and 0° < ψ < 140°. The numbers of visiting times are listed in Table 5. The P_II, C₅, α_R, and α_P states were already sampled well in the canonical MD simulation. Important states here are the α_L and C^ax₇ states, because they were not sampled in the canonical MD simulation. These states were visited many times when E₁ was chosen as E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, E₁ = E^slt–slt_tor + E^slt–slt_elec, E₁ = E^slt–slt_elec + E^slt–slt_LJ, or E₁ = E^slt–slt_elec. In particular, when E₁ was chosen as E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ, all of the P_II, C₅, α_R, α_P, α_L and C^ax₇ states were visited more than 400 times. It means that the simulation did not get trapped in any local-minimum free-energy states and visited all states many times. The electrostatic potential energy between the solute atoms E^slt–slt_elec was commonly chosen in these four choices of E₁. We can see that covering higher energy values of E^slt–slt_elec contributes the most to sampling so many states.

Table 5 The number of visiting times of each state

Choice of E1	PII	C5	α R	α P	α L	C^ax7
E 1 = 0	1318	1417	194	291	0	0
E 1 = E	1305	1423	266	382	8	9
E 1 = E^slt–slt + E^slt–svn	1421	1497	306	366	11	21
E 1 = E^slt–slt	205	208	2	3	0	0
E 1 = E^slt–slttor + E^slt–sltelec + E^slt–sltLJ	1019	452	810	515	1141	862
E 1 = E^slt–sltbond + E^slt–sltangle + E^slt–sltimp	1150	1170	112	131	0	0
E 1 = E^slt–slttor + E^slt–sltelec	1017	115	834	53	541	402
E 1 = E^slt–slttor + E^slt–sltLJ	1745	1390	1172	779	99	121
E 1 = E^slt–sltelec + E^slt–sltLJ	1102	691	592	282	816	732
E 1 = E^slt–slttor	1667	1661	740	759	46	24
E 1 = E^slt–sltelec	495	249	459	240	995	970
E 1 = E^slt–sltLJ	1509	1057	848	389	27	38
E 1 = E^slt–sltbond	1208	1309	282	377	1	1
E 1 = E^slt–sltangle	1246	1300	151	205	0	0
E 1 = E^slt–sltimp	1040	1108	145	216	2	1

---

On the other hand, the Lennard-Jones potential energy between solute atoms E^slt–slt_LJ, instead of E^slt–slt_elec, was most important for the numbers of the complete rotations of ϕ and ψ, as listed in Table 1. It depends on whether the simulation can overcome the free energy barrier at ϕ = 120° or not. The backbone dihedral angle ϕ made nine turns of the rotation when E₁ was chosen as E₁ = E^slt–slt_LJ, as shown in Fig. 6(l). It is more than one turn of the rotation in the case of E₁ = E^slt–slt_elec. However, the simulation visited the range of 0° < ϕ < 120° more frequently if E₁ = E^slt–slt_elec than E₁ = E^slt–slt_LJ. The reason why the number of rotation in the case of E₁ = E^slt–slt_elec is less than that in the case of E^slt–slt_LJ is the simulation of E₁ = E^slt–slt_elec did not overcome the free energy barrier at ϕ = 120°.

In order to discuss the relation between the free energy barrier and the atomic conformation, the potential of the mean force and some typical atomic conformations, of which the dihedral angle ψ is larger than zero (ψ > 0°), are shown in Fig. 9. Both the potential of the mean force and the atomic conformations were obtained from the partial multicanonical MD simulation in the case of E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ. The distance r between the O atom of the acetyl group and the O atom of the alanine at the P_II state is r = 3.92 Å. Here, the Lennard-Jones pair potential energy u_LJ is given by

where σ is the distance at which the interatomic potential is zero and ε is the depth of the potential well. For these O atoms, σ = 2.96 Å and the Lennard-Jones energy of this pair is u_LJ = −0.13 kcal mol⁻¹. The distance between the H atom which has a covalent bond with the N atom of the alanine and the O atom of the alanine is r = 3.26 Å. The Lennard-Jones parameter of this pair is σ = 2.01 Å and the Lennard-Jones energy is u_LJ = −0.012 kcal mol⁻¹. These facts mean that both pairs of the atoms at the P_II state do not collide with each other. The electrostatic potential u_elec between these pairs are u_elec = 27.3 kcal mol⁻¹ for the pair of the O atom of the acetyl group and the O atom of the alanine and u_elec = 15.7 kcal mol⁻¹ for the pair of the H atom and the O atom of the alanine.

Fig. 9 Typical conformations of the alanine dipeptide and contour maps of the potential of mean force A(ϕ, ψ) at T = 300 K obtained by the partial multicanonical MD simulation in the case of E₁ = E^slt–slt_tor + E^slt–slt_elec + E^slt–slt_LJ.

At the C₅ state, the distance between the H atom and the O atom of the alanine is r = 2.26 Å. It is relatively close to that of the P_II state, but they do not collide with each other yet, because σ = 2.01 Å and it is still shorter than r. Thus, the Lennard-Jones energy is low enough as u_LJ = −0.057 kcal mol⁻¹ and the electrostatic energy is u_elec = −22.7 kcal mol⁻¹. Because the distances of any pairs of atoms at the C₅ states are not too close with each other, the C₅ state is stable and its population is high as that of the P_II state, as listed in Table 4.

The population of the α_L state is much less than the P_II and C₅ states. The distance between the O atom of the acetyl group and the O atom of the alanine is 2.77 Å. This distance is shorter than the parameter σ = 2.96 Å of this pair. It means that these O atoms overlap slightly with each other, thus the Lennard-Jones energy of this pair is relatively high as u_LJ = 0.61 kcal mol⁻¹ and the electrostatic potential is u_elec = 38.7 kcal mol⁻¹. Both u_LJ and u_elec are higher than those of the P_II and C₅ states. This is the reason why the population of the α_L state is less than the P_II and C₅ states.

The distance between the O atom of the acetyl group and the O atom of the alanine of the transition state at ϕ = 0° between the P_II state and the α_L state is 2.57 Å. They overlap each other more than those of the α_L state. Thus, both Lennard-Jones energy and the electrostatic energy are higher than the α_L state as u_LJ = 2.62 kcal mol⁻¹ and u_elec = 41.7 kcal mol⁻¹. The reason for the free energy barrier at ϕ = 0° is both the Lennard-Jones energy and the electrostatic energy.

The distance between the O atom of the acetyl group and the H atom of the alanine side chain of the transition state between the α_L state and the C₅ state at ϕ = 120° is 2.20 Å. These atoms overlap most among any non-bonding atom pairs shown in Fig. 9. The parameter σ of this pair is σ = 2.80 Å and the Lennard-Jones energy is highest as u_LJ = 3.25 kcal mol⁻¹. On the other hand, the electrostatic energy is as low as u_elec = −5.2 kcal mol⁻¹ because the charges of these atoms have opposite signs. The reason for the free energy barrier at ϕ = 120° is solely the Lennard-Jones energy, not the electrostatic energy. It means that the Lennard-Jones term E^slt–slt_LJ should be included in E₁ to overcome the free energy barrier at ϕ = 120°.

The free energy barrier at ϕ = 0° is also caused by the Lennard-Jones repulsive force but it is not as high as that at ϕ = 120°. The electrostatic interaction also causes the free energy barrier at ϕ = 0°. It is why the simulation in which E₁ includes only the electrostatic interaction (E₁ = E^slt–slt_elec) can overcome the free energy barrier at ϕ = 0°. To overcome both free energy barriers at ϕ = 0° and ϕ = 120° efficiently, both E^slt–slt_LJ and E^slt–slt_elec should be included in E₁. In fact, both simulations in the cases of E₁ = E^slt–slt_elec + E^slt–slt_LJ + E^slt–slt_tor and E₁ = E^slt–slt_elec + E^slt–slt_LJ visited all states many times. Strictly speaking, the simulation sampling the wide range of E₁ = E^slt–slt_elec + E^slt–slt_LJ + E^slt–slt_tor visited all states slightly more than that of E₁ = E^slt–slt_elec + E^slt–slt_LJ, as listed in Table 5. This fact indicates that sampling the wide range of the torsion angle potential energy E^slt–slt_tor also contributes somewhat to enhancing the sampling efficiency.

V. Conclusions

In the multicanonical ensemble, a free one-dimensional random walk is realized in the potential-energy space and a simulation does not get trapped in the local-minimum free-energy states in biomolecules. But still there is the problem of the multicanonical algorithm that the determination of the multicanonical weight factor becomes difficult for a large system. It is because the multicanonical simulation samples a wide range of total potential energy and thus a redundant effort is required to determine the weight factor. In order to alleviate this difficulty, the author recently proposed the partial multicanonical algorithm for MD and MC simulations.²⁵ The partial multicanonical simulation samples a wide range of only a part of the potential energy terms which is necessary to sample the conformational space widely, so that one can save the effort to determine the weight factor. The partial multicanonical MD simulation can not only sample the local-minimum free-energy states but also overcome the free-energy barriers. Furthermore, its sampling efficiency is higher than the multicanonical algorithm by focusing to sample the wide range of the important energy terms.

To find the best choice of the partial potential energy E₁ which should be sampled widely, 15 choices of E₁ were tried in this article. The canonical, multicanonical, selectively enhanced multicanonical, and partial multicanonical MD simulations of an alanine dipeptide were performed in explicit water. The AMBER parm96 force field³⁸ was used for the peptide molecule and the TIP3P rigid-body model⁵³ was employed for the water molecules.

The canonical MD simulation sampled only four states of the P_II, C₅, α_R, and α_P states. The multicanonical and selectively enhanced multicanonical MD simulations also sampled the α_L and C^ax₇ states. Although these simulations sampled all of the six states, the α_L and C^ax₇ states were not visited frequently and their sampling was not enough.

The best choice of E₁ in this system is the sum of the electrostatic, Lennard-Jones, and torsion-angle terms between the solute atoms. In this case, all states were sampled sufficiently many times. The most important term is the electrostatic potential energy, which makes the simulation visit all states many times. The Lennard-Jones potential is also important for overcoming the steric restriction. In the case of an alanine dipeptide, there are free energy barriers at ϕ = 0° and ϕ = 120° due to the steric restriction. Including the torsion-angle potential energy also contributes to enhancing the sampling efficiency.

The partial potential energy E₁ in the multicanonical and selectively enhanced multicanonical algorithms also includes the electrostatic, Lennard-Jones, and torsion-angle terms. But these algorithms do not concentrate on widening the sampling energy ranges of these terms. Most of the effort to widen the sampling energy-ranges is spent on other unimportant energy terms. This is why the sampling efficiency of the multicanonical and selectively enhanced multicanonical algorithms is not as good as the partial multicanonical algorithm in the cases of E₁ = E^slt–slt_elec + E^slt–slt_LJ + E^slt–slt_tor, E₁ = E^slt–slt_elec + E^slt–slt_LJ, and E₁ = E^slt–slt_elec.

The best choice of E₁ may depend on the system and force field. For example, in the case of aquaporin, a protein which includes water molecules, the solute–solvent interaction may be important and should be included in E₁. However, the detailed discussion on the choice of E₁ in this article will be useful for other systems. The partial multicanonical algorithm, of course, will be also useful for a MC simulation. This algorithm will be a great use for MD and MC simulations of a large system like a protein beyond the limit of the multicanonical algorithm. Furthermore the partial multicanonical algorithm can be employed for free-energy sampling problems⁵⁹ in general such as supercooled liquids and glasses.

This idea, in which a simulation covers a wide range of only important potential energy terms, can be utilized not only for the multicanonical algorithm but also for multidimensional extensions^{14,15,18–24,30,31} of the multicanonical algorithm such as the multibaric–multithermal^18–24 and multicanonical–multioverlap³⁰ algorithms. Combining these algorithms, we can develop new stronger simulation techniques.

References

1. S. Nosé, Mol. Phys., 1984, 52, 255.
2. S. Nosé, J. Chem. Phys., 1984, 81, 511.
3. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695.
4. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087.
5. U. H. E. Hansmann and Y. Okamoto, in Ann. Rev. Comput. Phys. VI, ed. D. Stauffer, World Scientific, Singapore, 1999, p. 129.
6. A. Mitsutake, Y. Sugita and Y. Okamoto, Biopolymers, 2001, 60, 96.
7. B. A. Berg, Comput. Phys. Commun., 2002, 147, 52.
8. Y. Okamoto, J. Mol. Graphics Modell., 2004, 22, 425.
9. S. G. Itoh, H. Okumura and Y. Okamoto, Mol. Simul., 2007, 33, 47.
10. B. A. Berg and T. Neuhaus, Phys. Lett. B, 1991, 267, 249.
11. B. A. Berg and T. Neuhaus, Phys. Rev. Lett., 1992, 68, 9.
12. U. H. E. Hansmann, Y. Okamoto and F. Eisenmenger, Chem. Phys. Lett., 1996, 259, 321.
13. N. Nakajima, H. Nakamura and A. Kidera, J. Phys. Chem. B, 1997, 101, 817.
14. C. Bartels and M. Karplus, J. Comput. Chem., 1997, 18, 1450.
15. J. Higo, N. Nakajima, H. Shirai, A. Kidera and H. Nakamura, J. Comput. Chem., 1997, 18, 2086.
16. Y. Iba, G. Chikenji and M. Kikuchi, J. Phys. Soc. Jpn., 1998, 67, 3327.
17. N. Nakajima, Chem. Phys. Lett., 1998, 288, 319.
18. H. Okumura and Y. Okamoto, Chem. Phys. Lett., 2004, 383, 391.
19. H. Okumura and Y. Okamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 026702.
20. H. Okumura and Y. Okamoto, J. Phys. Soc. Jpn., 2004, 73, 3304.
21. H. Okumura and Y. Okamoto, Chem. Phys. Lett., 2004, 391, 248.
22. H. Okumura and Y. Okamoto, J. Comput. Chem., 2006, 27, 379.
23. H. Okumura and Y. Okamoto, Bull. Chem. Soc. Jpn., 2007, 80, 1114.
24. H. Okumura and Y. Okamoto, J. Phys. Chem. B, 2008, 112, 12038.
25. H. Okumura, J. Chem. Phys., 2008, 129, 124116.
26. B. A. Berg, H. Noguchi and Y. Okamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 036126.
27. S. G. Itoh and Y. Okamoto, Chem. Phys. Lett., 2004, 400, 308.
28. S. G. Itoh and Y. Okamoto, J. Chem. Phys., 2006, 124, 104103.
29. S. G. Itoh and Y. Okamoto, Mol. Simul., 2007, 33, 83.
30. S. G. Itoh and Y. Okamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 026705.
31. Y. Yoshimoto, J. Chem. Phys., 2006, 125, 184103.
32. T. Yamaguchi, T. Kiuchi, T. Matsuoka and S. Koda, Bull. Chem. Soc. Jpn., 2005, 78, 2098.
33. T. Morishita and M. Mikami, J. Chem. Phys., 2007, 127, 034104.
34. B. A. Berg and T. Celik, Phys. Rev. Lett., 1992, 69, 2292.
35. J. Lee, Phys. Rev. Lett., 1993, 71, 211; J. Lee, Phys. Rev. Lett., 1993, 71, 2353 (E).
36. Y. Okamoto and U. H. E. Hansmann, J. Phys. Chem., 1995, 99, 11276.
37. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179.
38. P. A. Kollman, R. Dixon, W. Cornell, T. Fox, C. Chipot and A. Pohorille, in Computer Simulation of Biomolecular Systems 3, ed. A. Wilkinson, P. Weiner and W. F. van Gunsteren, Elsevier, Dordrecht, 1997, p. 83.
39. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225.
40. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2001, 105, 6474.
41. A. D. MacKerell, Jr., D. Bashford, M. Bellott, R. L. Dunbrack, Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, J. Phys. Chem. B, 1998, 102, 3586.
42. I. R. Gould and P. A. Kollman, J. Phys. Chem., 1992, 96, 9255.
43. J. Apostolakis, P. Ferrara and A. Caflisch, J. Chem. Phys., 1999, 110, 2099.
44. P. E. Smith, J. Chem. Phys., 1999, 111, 5568.
45. A. D. MacKerell, Jr., M. Feig and C. L. Brooks III, J. Comput. Chem., 2004, 25, 1400.
46. J. D. Chodera, W. C. Swope, J. W. Pitera and K. A. Dill, Multiscale Model. Simul., 2006, 5, 1214.
47. D. Branduardi, F. L. Gervasio and M. Parrinello, J. Chem. Phys., 2007, 126, 054103.
48. D. Branduardi, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 020603.
49. M. Bonomi, A. Barducci and M. Parrinello, J. Comput. Chem., 2009, 30, 1615.
50. A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1988, 61, 2635; A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1989, 63, 1658 (E).
51. A. Mitsutake, M. Kinoshita, Y. Okamoto and F. Hirata, Chem. Phys. Lett., 2000, 329, 295.
52. R. A. Sayle and E. J. Milner-White, Trends Biochem. Sci., 1995, 20, 374.
53. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926.
54. H. Okumura, S. G. Itoh and Y. Okamoto, J. Chem. Phys., 2007, 126, 084103.
55. T. F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns and G. J. Martyna, J. Chem. Phys., 2002, 116, 8649.
56. F. Wang and D. P. Landau, Phys. Rev. Lett., 2001, 86, 2050.
57. B. A. Berg, Introduction to Monte Carlo Simulations and Their Statistical Analysis, World Scientific, Singapore, 2004.
58. T. Yoda, Y. Sugita and Y. Okamoto, Chem. Phys. Lett., 2004, 386, 460.
59. J. C. Schön, J. Chem. Phys., 1996, 105, 10072.

---

Footnote

† This article was submitted as part of a Themed Issue on solid state and cluster structure prediction. Other papers on this topic can be found in issue 30 of vol. 12 (2010). This issue can be found from the PCCP homepage [http://www.rsc.org/pccp]

---