Ubiquitin in a Water Sphere: Simulation with Non-Periodic Boundary Conditions

Next: Ubiquitin in a Water Up: Basics of NAMD Previous: Solvating the Protein

Ubiquitin in a Water Sphere: Simulation with Non-Periodic Boundary Conditions

In this section, you will examine the minimization and equilibration of ubiquitin in a water sphere placed in vacuum.

Configuration File

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...our file is correct before you submit a simulation to be run!} \end{minipage} }$

$% latex2html id marker 4320 \fbox{ \begin{minipage}{.2\textwidth} \includegra... ...hort on time, and examine the contents of this section later.} \end{minipage} }$

1: Go to your 1-2-sphere directory by typing cd ..\1-2-sphere. Here, you will find a configuration file for the minimization and equilibration of ubiquitin in a water sphere.

All output files for the minimization and equilibration of your ubiquitin in a water sphere system will be placed in this directory. The configuration file is the only input file placed here, since it is particular to this simulation. The pdb, psf, and parameter files, which may be used by many other simulations, are placed in your common directory and are called by each respective configuration file.

Open WordPad by clicking Start $\rightarrow$ Programs $\rightarrow$ Accessories $\rightarrow$ WordPad. Open the configuration file, ubq_ws_eq.conf in your 1-2-sphere directory.

The configuration file may seem complex at first, but it will be examined line by line to determine its function in your simulation.

Note that when ``#" appears at the beginning of a line, the entire line is treated as a comment and ignored by NAMD. In the middle of a line, ``;#" is used to comment out the remainder of the line.

2

The ``Job Description" section contains only comments, and its purpose is to inform those who view the configuration file about what it is meant for. Your comment should read

# Minimization and Equilibration of

# Ubiquitin in a Water Sphere

3

The ``Adjustable Parameters" section contains five commands:

structure: calls the psf file describing the system ( ..\common\ubq_ws.psf).
coordinates: calls the initial coordinate data from the file listed next to it (in this case ..\common\ubq_ws.pdb).
set temperature: creates a variable called ``temperature" in which to store a value for the initial temperature of the system. If you place the text ``$temperature" in the configuration file, NAMD simply reads it as a label for the number ``310". (Creating variables is useful since you can alter the value of that variable in a single place if it is listed many times in the configuration file.)
set outputname: creates a variable called ``outputname" in which to store a generic name for output files. If you place the text ``$outputname" in the configuration file, NAMD simply reads it as a label for ``ubq_ws_eq".
firsttimestep: simply sets a number value for the first time step of the simulation. It is typically useful when restarting a simulation. For instance, if a previous simulation ended at time step 552, the command firsttimestep 553 would be used.

4

The ``Simulation Parameters" section contains many commands, commented into different categories:

Input
- paraTypeCharmm: indicates whether or not the parameter file is in the format used by the CHARMM force field. on indicates that it is; off indicates that it is not. (If this command is not specified, NAMD assumes the file is in X-PLOR format by default.)
- parameters: calls the force field parameters from the file listed next to it (in this case ..\common\par_all27_prot_lipid.inp).
- temperature: sets the initial temperature of the system in Kelvin (K), with the value listed next to it (in this case $temperature, or 310). This is done by assigning random velocities to atoms picked from a Maxwell distribution such that their average kinetic energy accurately represents the given temperature.
Force-Field Parameters
- exclude: specifies which atomic interactions are to be excluded from consideration. See Figure 4 for general atom labels. The scaled1-4 value indicates that interactions between any such atoms 1 and 2 and 1 and 3 are neglected and interactions between atoms 1 and 4 are weakened. The van der Waals interaction for ``1-4" atoms are modified using special 1-4 parameters defined in the parameter files, and electrostatic interaction is modified as shown in the next command.
  
  Figure 4: Number labels for atoms which are a given amount of bonds away from one another
  $\begin{figure}\begin{center} \par \par \latex{ \includegraphics[scale=0.5]{pictures/tut_unit01_bonds} } \end{center} \end{figure}$
- 1-4scaling: specifies the degree to which the electrostatic interaction between 1-4 atoms is to be taken into account. It may be a decimal between 0 and 1 (here, it is 1) and indicates how much the interaction is ``turned off" or ``on", respectively.
- cutoff: indicates the distance in Å beyond which electrostatic and van der Waals interactions are cut-off. Otherwise, those interactions are considered over the entire volume of your system. This can be computationally too costly. When also employing a fast solver for electrostatics such as Particle Mesh Ewald (Section 1.5.1) or Multilevel Summation Method (Section 1.4.1), the cutoff parameter instead defines the splitting distance for the interaction potential, in which a short-ranged part is evaluated exactly between atoms within the cutoff distance, and the remaining long-ranged part is approximated by the solver.
- switching: indicates whether or not switching functions are used to smoothly take electrostatic and van der Waals interactions to zero at the cutoff distance. You may list ``on" or ``off" next to it as a yes/no answer.
- switchdist: indicates the distance in Å at which the functional form of electrostatic and van der Waals interactions is modified to allow their values to approach zero at the cutoff distance. A visual explanation of this is useful and may be found in Figure 5.
- pairlistdist: is designed to make computation faster. It specifies a distance in Å. NAMD will only search within this distance for atoms which may interact by electrostatic or van der Waals interactions. This way, NAMD does not have to search the entire system. The distance must be greater than the cutoff distance, and the list must be updated during the simulation. See Figure 5.
Figure 5: Cutoff and switching distances indicated on the left. Pair list distance indicated on the right.
$\begin{figure}\begin{center} \par \par \latex{ \includegraphics[scale=0.5]{pictures/tut_unit01_cutoff} } \end{center} \end{figure}$
Integrator Parameters
- timestep: indicates the value of the time step size used in the simulation. MD simulations solve Newton's laws in a discrete approximation to determine the trajectories of atoms. The time step tells NAMD how to discretize the particle dynamics. It is specified in femtoseconds (here, 2 fs).
- rigidBonds: specifies which bonds involving hydrogen are considered to be rigid (non-vibrating). The value all specifies all linear bonds involving hydrogen and any other atoms.
- nonbondedFreq: specifies in number of time steps how often nonbonded interactions should be calculated. It is useful for saving computational time.
- fullElectFrequency: specifies in number of time steps how often full electrostatic interactions should be calculated.
- stepspercycle: Atoms are reassigned pair list identities (as explained above) once every cycle. This command specifies how long one cycle lasts, i.e. the number of time steps in one cycle.
Constant Temperature Control
- langevin: indicates whether or not the simulation uses Langevin dynamics; uses values on and off. See the science box below for more on Langevin dynamics.
- langevinDamping: sets the value of the Langevin coupling coefficient, which quantifies the friction applied to the system, removing energy from the system, slowing atoms down, etc. It is specified in ps $^{-1}$ .
- langevinTemp: Langevin dynamics may be applied to all atoms or only non-hydrogen atoms in the system. This command specifies the temperature at which to keep those atoms, even though friction and random forces will be acting on them. (Remember $temperature is a variable for the value 310.)
- langevinHydrogen: indicates whether or not Langevin dynamics will be applied to hydrogen atoms in the simulation; uses values on and off.
Electrostatics with MSM
Multilevel Summation Method (MSM) is a useful method for dealing with electrostatic interactions in a system that does not have fully periodic boundary conditions along all three dimensions, in other words, when employing non-periodic or semi-periodic boundaries. Although MSM also supports fully periodic boundary conditions, the Particle Mesh Ewald (PME) method discussed in Section 1.5 is more efficient for smaller periodic systems, such as ubiquitin simulated in a water box. In the present simulation involving a water sphere, MSM is the best alternative for electrostatics.
MSM uses particle meshes, similar to PME, arranged as 3-D grids for evaluating long-range electrostatic interactions. Charges are assigned to nearby grid points, and the interaction of these gridded charges result in gridded potentials that are interpolated back to the atoms to determine the long-range atomic forces. Unlike PME which uses a single grid resolution, MSM uses a nesting of grids having coarser resolution at each successive ``level'' in real space for resolving the long-ranged tail of the Coulomb potential. Rather than evaluating the Ewald sum as PME does, MSM uses a continued splitting of the interaction potential and performs interpolation from these nested grid levels.
- MSM: indicates whether or not the simulation uses the Multilevel Summation Method; uses values yes and no.
- MSMGridSpacing: The grid spacing between the finest level grid points determines, in part, the accuracy and efficiency of MSM. According to an error versus cost analysis, the optimal value is for the grid spacing to be a bit larger than the average interatomic spacing; since the performance is extremely sensitive to the grid spacing, it is best to leave this parameter set to its default value of 2.5Å.
- MSMxmin, MSMxmax, MSMymin, MSMymax, MSMzmin, MSMzmax: For non-periodic simulation, MSM needs to know the expected maximum and minimum coordinates so as to establish grids that will contain all atoms throughout the entire simulation. To keep the atoms from leaving the simulated space, it might be necessary to use a restraining potential, such as the spherical boundary restraints used for the water sphere simulation.
Additional MSM parameters that further control accuracy and performance are documented in the NAMD User's Guide.
Output
- outputName: Several types of output data may be written by NAMD for any simulation. This command specifies the prefix for output filenames. NAMD always returns the following two output files for every simulation: a pdb file containing the final coordinates of all atoms in the system and a pdb file containing the final velocities of all atoms in the system; the extensions for these files are .coor and .vel, respectively. Thus, this configuration file will write two files named ubq_ws_eq.coor and ubq_ws_eq.vel
- restartfreq: During the simulation, NAMD can also create restart files, one of which is a pdb file which stores atomic coordinates, and the other of which stores atomic velocities. This command specifies the amount of time steps between writing to the restart file (here, every 500 steps, or 1000fs, or 1 ps). If this command is not set, NAMD will not create a restart file. Furthermore, NAMD will store the file from the previous cycle each time it writes a new file. The filename is appended with a .old extension; it is created in case NAMD fails in writing the new restart file.
- dcdfreq: The dcd file contains only atomic coordinates, and they are written to the file several times over the course of a simulation. Thus, it provides a trajectory of the system over the runtime. This command specifies the number of time steps between writing new coordinates to the dcd output file. If this command is not set, NAMD will not create a dcd file.
  
  In addition to the output files described, NAMD also prints a log of the simulation (which we will redirect into a .log file). This output is explained further below.
- outputEnergies: specifies the number of time steps between each output of system energies (for various force field interactions) into the .log file (here, every 100 steps, or 200 fs).

**Figure 4:** Number labels for atoms which are a given amount of bonds away from one another
$\begin{figure}\begin{center} \par \par \latex{ \includegraphics[scale=0.5]{pictures/tut_unit01_bonds} } \end{center} \end{figure}$

**Figure 5:** Cutoff and switching distances indicated on the left. Pair list distance indicated on the right.
$\begin{figure}\begin{center} \par \par \latex{ \includegraphics[scale=0.5]{pictures/tut_unit01_cutoff} } \end{center} \end{figure}$

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...ture. \textbf{A value of 1.0 is often a good starting point.}} \end{minipage} }$

$\framebox[\textwidth]{ \begin{minipage}{.2\textwidth} \includegraphics[width=2... ...se water molecules have been parametrized as rigid molecules.} \end{minipage} }$

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ... simulations using non-periodic or semi-periodic boundaries.} \end{minipage} }$

5

The ``Extra Parameters" section contains commands which are applicable to more specific simulations. Included here are commands which characterize the spherical boundary conditions on the water sphere. These conditions prevent the sphere from undergoing evaporation or diffusion.

Spherical Boundary Conditions
- sphericalBC: indicates whether or not the simulation uses spherical boundary conditions; uses values on and off.
- sphericalBCcenter: sets the , , coordinates of the center of the sphere to which spherical boundary conditions are applied. This is why you need to record the center of your water sphere after you solvate!
The next three commands are parameters which must be specified for spherical boundary conditions to work properly.
- sphericalBCr1: sets the distance in Å at which the first boundary potential begins to act.
- sphericalBCk1: The first potential applied by spherical boundary conditions to keep the sphere together (or pull it apart) is a harmonic one. This command sets the value of the force constant for that potential. Its value is specified in kcal/mol $\cdot$ Å $^{2}$ .
- sphericalBCexp1: sets the value of the exponent used to formulate the potential; must be a positive, even integer.

6

The ``Execution Script" section contains three commands, the first two of which apply to minimization and the last one of which applies to equilibration.

Minimization
- minimize: sets the number of iterations over which to vary atom positions to search for a local minimum in the potential (in this case 100).
- reinitvels: Minimization is performed on the system after all atomic velocities have been set to zero. This command resets atomic velocities such that the system starts at the temperature specified (in this case, $temperature, or 310K).
run: sets the number of time steps over which to run the MD equilibration (in this case 2500, which corresponds to 5,000 fs or 5 ps, since a 2 fs time step has been used).

7: Now, close the configuration file by clicking File $\rightarrow$ Exit.

$% latex2html id marker 4604 \framebox[\textwidth]{ \begin{minipage}{.2\textwid... ...rent/ug/}{http://www.ks.uiuc.edu/Research/namd/current/ug/}}.} \end{minipage} }$

$\framebox[\textwidth]{ \begin{minipage}{.2\textwidth} \includegraphics[width=2... ...owing: {\tt temperature, velocities,} or {\tt binvelocities}.} \end{minipage} }$

Run your Simulation

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...ow: \\ \\ {\tt namd2 ubq\_ws\_eq.conf > ubq\_ws\_eq.log \&} } \end{minipage} }$

$\framebox[\textwidth]{ \begin{minipage}{.2\textwidth} \includegraphics[width=2... ...are in, you should be able to type {\tt namd2} and run it. } \end{minipage} }$

Your simulation should take about 20 minutes to complete, and NAMD will produce the output files described in Section 1.7.

Next: Ubiquitin in a Water Up: Basics of NAMD Previous: Solvating the Protein

namd@ks.uiuc.edu

`# Minimization and Equilibration of`
`# Ubiquitin in a Water Sphere`