Timing results for the CMM molecular dynamics calculations are shown in Figure , in terms of cpu seconds on one processor of an SGI 4D/380 workstation. The total charge, dipole, and quadrupole of each cell,
collectively termed the ``farfield,'' can be recalculated every step or can be considered constant for a number of steps. These two cases are labeled ``Update1'' and ``Update50,'' the latter referring to calculations in which the farfield was updated only every 50 steps. Also shown is the difference between calculations using only the nonbonds of the three-protein asymmetric unit, labeled ``NB3,'' and those including interactions with the entire 180-protein capsid, labeled ``NB180.'' It is interesting to note that the Update50 calculations are actually faster for NB180 than NB3 (37.2 versus 41.0 s). This is because the calculation is dominated by the nearfield interactions. As noted above, the average cell in the NB180 case has 5.4 atoms versus 6.4 for the NB3 case. This means fewer pairwise interactions need to be calculated. The farfield interaction takes longer to calculate in the NB180 case, but the effect is not significant because of the hierarchical approach described above. Only when the farfield itself must be updated every step, i.e., the Update1 calculations, does the size of the system make a significant difference. In such calculations, including the entire capsid increases the time of the nonbond calculations from 49.8 s to 92.6 s.
Figure also shows the total time required for NEIMO calculations, including the time required to calculate accelerations and velocities, and to update the system coordinates. Unlike other internal-coordinate methods, NEIMO dynamics are computationally proportional to , the number of internal degrees of freedom. Therefore, NEIMO can be used to calculate internal coordinate dynamics for a system this large, an almost impossible task for other internal-coordinate methods, which typically have computational costs proportional to . As is clear in Figure , the linear scaling with respect to makes NEIMO entirely feasible for simulations of TBSV. NEIMO calculations for the three protein chains in the asymmetric unit ( = 4335) require 5.3 s per dynamics step. The primary advantage of internal-coordinate methods is that large timesteps can be used; whereas Cartesian dynamics simulations typically require timesteps of 1 fs for accuracy, timesteps as large as 15 fs can be used reliably for small peptides. Although the current NEIMO implementation does not allow timesteps of this magnitude for large proteins, timesteps larger than 1 fs are practical. As is indicated in Figure , the overhead for NEIMO calculations is very small compared to the time required for the nonbond calculations. Therefore, an increase in step size just to 2 fs will nearly double the speed of the simulation with respect to a 1 fs Cartesian simulation.
Figure shows the average scaled energy fluctuation,
*, versus timestep for calculations on the PH7 model of TBSV. * is the average value of during a simulation, divided by the number of degrees of freedom. is determined over 0.100 ps intervals using the equation:
where is the total energy (potential plus kinetic), is the Boltzmann constant, and is the temperature of the simulation. The figure indicates that no choice for nonbonds gives superior results at all timesteps. Fluctuations using the entire capsid, with the farfield updated every step, i.e., NB180/Update1, were not measured but are unlikely to provide a substantial improvement. In every case, there was typically a 4- to 5-fold increase in the fluctuation for each 1 fs increase in the timestep. The average scaled fluctuation, *, shown in Figure , is approximately the same for 1 fs NEIMO and 1 fs Cartesian dynamics. However, the unscaled fluctuation, i.e., the average value of not scaled by the number of degrees of freedom, is of the same order of magnitude for 2 fs NEIMO and 1 fs Cartesian dynamics.
The two model systems, PH7 and NOCA, were initially optimized using Cartesian space conjugate gradients minimization. The nonbond calculations included the entire capsid, with the farfield updated every 50 steps. The radius of gyration of the entire 180 protein system was calculated every 50 steps, when the farfield was updated. The radius of gyration is defined as:
where the coordinates and mass of each particle are and , respectively, and the coordinates of the center of mass are ). This is shown in Figure
for the first 1000 steps of minimization. PH7 and NOCA actually took 1300 and 1176 steps to minimize, respectively. Both structures contracted by about 0.08%during the minimization with their radii following almost identical curves. The contraction rate remained constant after the first 100 steps, indicating it would have continued had the minimization gone longer. However, curves of the potential energy during minimization, shown in Figure , clearly indicate
that the energy had converged. The PH7 and NOCA energies also had similar curves. However, the PH7 structure is almost 700 kcal/mol lower in energy than NOCA, despite containing identical numbers of atoms. This energy difference is almost entirely due to the electrostatic energy term and indicates the large stabilizing energy of the Ca ions.
Molecular dynamics calculations were carried out on the different minimized structures, again using nonbonds from the full capsid. The farfield was updated every 0.100 ps. Two different Cartesian dynamics methods were tried: microcanonical ensemble (NVE) with temperature scaling and canonical ensemble (NVT). In addition, NEIMO dynamics, which are NVE, were used with a timestep of 2 fs. Figure shows
the radius of gyration of PH7 during the first 2.0 ps of dynamics. All three methods give different curves but there are some similarities. Both Cartesian simulations have an initial expansion phase followed by a longer contraction. The NEIMO simulation has no expansion phase, but its contraction phase closely resembles that of the Cartesian NVE simulation, with roughly the same slope, contracting until approximately 1.8 ps, when it levels out. The canonical dynamics simulation, in contrast, shows no similar leveling out through the first 2.0 ps.
Longer simulations were run using Cartesian canonical dynamics and NEIMO dynamics on both the PH7 and NOCA models. The NEIMO dynamics simulations were twice as fast, since 2 fs timesteps were used. Therefore, longer simulations could be run. Figure shows the radius of
gyration of the PH7 and NOCA models during the first 4.0 ps of NEIMO and Cartesian NVT simulations. In the NEIMO simulations, the NOCA model undergoes a rapid expansion during the first 2.0 ps, then an even sharper contraction. The PH7 model has no such expansion phase but does have a gradually increasing contraction rate. Both of these simulations show far more variation in their radius of gyration than the corresponding Cartesian simulations. However, in both simulations, the NOCA model initially has a larger radius of gyration than for PH7, but eventually becomes smaller. Although the curves of the radii cross after about 3.9 ps in the Cartesian simulation, the energy curves do not cross, as shown in Figure .
The energy plot indicates that the NOCA model is less stable, as it is undergoes larger energy fluctuations after 3.0 ps. These fluctuations continue until the end of a 5.0 ps simulation (data not shown). The PH7 model is relatively stable. For the NEIMO simulations, the contraction rate is comparatively exaggerated, but the energies do not show such large fluctuations. The NOCA model has a potential energy around -3850 kcal/mol while the PH7 model's potential remains near -4550 kcal/mol. Note that these energies are substantially lower than in the Cartesian simulation because the numerous bond and angle degrees of freedom are held at their minimum potential energy values. Therefore, the approximately 700 kcal/mol differential between PH7 and NOCA is relatively constant, even though the radii change significantly.