|
|
| A two-dimensional array of boxes. As molecule 1 moves from the central box into box B it is replaced by it's image which moves from box F into the central box. This movement is replicated across all the boxes. |
(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright Ross Walker 2004
5) Running Minimisation and MD (in explicit solvent)
Up to this point we have been running simulations based on our in vacuo 10-mer polyA-polyT DNA structure. With this structure we have run both in vacuo simulations and also generalised Born implicit solvent simulations. In the next part of this tutorial we are going to run simulations in explicit solvent with periodic boundaries. This significantly increases the complexity of the simulations increasing both the time required to run the simulations and also the way in which we set up the simulations and visualise the results.
5.1) Equilibrating and running MD with solvated poly(A)-poly(T)
The problems with bad van der Waal (non bond) and electrostatic interactions in the initial structure, that led to the need to minimise our in vacuo structure, are considerably magnified when we start to run simulations in explicit solvent. Our initial structure was a gas phase DNA molecule with no protons. xleap added these for us at pre-determined positions. We then added 18 sodium ions at positions of high negative electric potential around our DNA molecule in order to neutralise it. To this system we then added a truncated octahedral box of pre-equilibrated TIP3P water. This water has not felt the influence of the solute or charges and moreover there may be gaps between the solvent and solute and solvent and box edges. If we are not careful such holes can lead to "vacuum" bubbles forming and subsequently an instability in our molecular dynamics simulation. Thus we need to run careful minimisation before slowly heating our system to 300 K. It is also a good idea to allow the water box to relax during a MD equilibration stage prior to running production MD. In this phase it is a good idea, since we will be using periodic boundaries, to keep the pressure constant and so allow the volume of the box to change. This approach will allow the water to equilibrate around the solute and come to an equilibrium density. It is essential that we monitor this equilibrium phase in order to be certain our solvated system has reached equilibrium before we start obtaining results (production data) from our MD simulation.
In this section will be presented one possible path towards equilibrating the solvated DNA structure such that it is ready for use in production MD. For this we will use periodic boundary simulations based on the particle mesh Ewald (PME) method. After equilibration has been demonstrated and deemed to have been successful, we will set up a "production" MD run. Please note, there are many differing opinions on how one should equilibrate a simulation, herein is presented only one, based on my experience with performing molecular dynamics simulations on biological systems.
Note: some of the calculations, particularly the equilibration phase, may take a very long time (several days) to run unless you have access to a computer cluster for parallel simulation. For the purposes of this simulation I will provide copies of all of the input files along with example output files run on four 1.3GHz CPUs of a SGI Altix. I will make it clear where the simulations are likely to take a long time such that you can just use the provided output files.
r + ix + jy + kz (i,j,k -> -inf, +inf)
where x, y and z are the vectors corresponding to the box edges. During the simulation only one of the particles is represented, but the effects are reproduced over all the image particles with each particle not only interacting with the other particles but also with their images in neighbouring boxes. Particles that leave one side of the box re-enter from the opposite side as their image. In this way the total number of particles in the central box remains constant.Upon initial inspection such a method would appear to be very computationally intensive requiring the evaluation of an infinite number of interacting pairs. However, by choosing the non-bonded cutoff distance such that each atom ‘sees’ only one image of all of the other atoms the complexity is greatly reduced. Hence as long as the box size is more than twice the cut-off distance it is impossible for a particle to interact with any two images of the same particle simultaneously. This cutoff requirement is important for us to remember when we choose what cutoff to use later on. It must not be larger than half the smallest box dimension. If you look back to when we initially created our solvated structure (here), you will see that the box size reported by xleap was 60.281 x 60.281 x 60.281 angstroms. Thus our maximum cutoff size is approximately 30.1 angstroms.
|
|
| A two-dimensional array of boxes. As molecule 1 moves from the central box into box B it is replaced by it's image which moves from box F into the central box. This movement is replicated across all the boxes. |
The choice of cutoff is therefore an important factor in carrying out a MD simulation. The figure below shows the variation in total system energy after 100 steps of steepest descent minimisation of our solvated DNA 10-mer as a function of cutoff distance. It can be seen that the total energy converges quickly and that by 10 angstroms the gain in going to large cutoffs is small and no longer justifies the extra computational cost. Hence for the solvated simulations that will be carried out in the remainder of this tutorial a cutoff of 10 angstroms will be utilised.
 |
|
Energy of 10-mer DNA after 100 steps of steepest descent minimisation (A), and processor time (1 x 3 GHz P4) in seconds (B), as a function of cutoff distance. |
There are two options for keeping the DNA atoms fixed during the initial minimisation, one method is to use the "BELLY" option in which all the selected atoms have the forces on them zeroed at each step. This results in what is essentially a partial minimisation. This method used to be the favoured method but it can, especially if used with constant pressure periodic boundary simulations, lead to instabilities and strange artefacts. As a result of this it is no longer the suggested method. Instead we shall use "positional restraints" on each of the DNA atoms to keep them essentially fixed in the same position. Such restraints work by specifying a reference structure, in this case our starting structure, and then restraining the selected atoms to conform to this structure via the use of a force constant. This can be visualised as basically attaching a spring to each of the solute atoms that connects it to its initial position. Moving the atom from this position results in a force which acts to restore it to the initial position. By varying the magnitude of the force constant so the effect can be increased or decreased. This can be especially useful with structures based on homology models which may be a long way from equilibrium and so require a number of stages of minimisation / MD with the force constant being reduced each time.
Here is the input file we shall use for our initial minimisation of solvent and ions:
|
polyAT_wat_min1.in |
polyA-polyT 10-mer: initial minimisation solvent + ions &cntrl imin = 1, maxcyc = 1000, ncyc = 500, ntb = 1, ntr = 1, cut = 10 / Hold the DNA fixed 500.0 RES 1 20 END END |
The meaning of each of the terms are as follows:
Hold the DNA fixed
500.0
RES 1 20
END
ENDNote that whenever you run using the GROUP option in the input you should carefully check the top of the output file to make sure you've selected as many atoms as you thought you did. (Note also that 500 kcal/mol-A**2 is a very large force constant, much larger than is needed. You can use this for minimization, but for dynamics, stick to much smaller values, say around 10.)
----- READING GROUP 1; TITLE: Hold the DNA fixed GROUP 1 HAS HARMONIC CONSTRAINTS 500.00000 GRP 1 RES 1 TO 20 Number of atoms in this group = 638 ----- END OF GROUP READ -----
Ok, we are ready so let's run the minimisation. Note that we have an extra option on the command line (-ref). This specifies the structure to which we want to restrain the atoms, in this case our initial structure in the inpcrd file. We will be using the solvated prmtop and inpcrd files we created at the beginning of this tutorial.
$AMBERHOME/exe/sander -O -i polyAT_wat_min1.in -o polyAT_wat_min1.out -p polyAT_wat.prmtop -c polyAT_wat.inpcrd -r polyAT_wat_min1.rst -ref polyAT_wat.inpcrd
Input files:
polyAT_wat_min1.in,
polyAT_wat.prmtop,
polyAT_wat.inpcrd
Output files:
polyAT_wat_min1.out,
polyAT_wat_min1.rst
This takes about 350 s (5 mins, 50 s) on a 3 GHz Pentium 4, so you may want to go get a coffee.
Take a look at the output file from this job. You should notice that there is initially a rather high van der Waals (VDWAALS and 1-4 VDW) energy which represents bad contacts in both the water and DNA. You should also see that the electrostatic energy (EEL) drops very rapidly. This is due to optimising of the favourable water-water and water-DNA electrostatic interactions as the waters re-orientate themselves into a lower energy geometry. The RESTRAINT energy represents the effect of the harmonic restraints we have imposed. It rises rapidly but then levels off.
Now we have minimised the water and ions the next stage of our minimisation is to minimise the entire system. In this case we will run 2,500 steps of minimisation without the restraints this time. Here is the input file:
|
polyAT_wat_min2.in |
polyA-polyT 10-mer: initial minimisation whole system &cntrl imin = 1, maxcyc = 2500, ncyc = 1000, ntb = 1, ntr = 0, cut = 10 / |
Note, choosing the number of minimisation steps to run is a bit of a black art. Running too few can lead to instabilities when you start running MD. Running too many will not do any real harm since we will just get ever closer to the nearest local minima. It will, however, waste cpu time. Since minimisation is generally very quick with respect to the molecular dynamics it is often a good idea to run more minimisation steps than really necessary. Here, 2,500 should be more than enough.
We run this using the following command. Note, we use the restrt file from the previous run, which contains the last structure from the first stage of minimisation, as the input structure for stage 2 of our minimisation. Note, we also no longer need the -ref switch:
$AMBERHOME/exe/sander -O -i polyAT_wat_min2.in -o polyAT_wat_min2.out -p polyAT_wat.prmtop -c polyAT_wat_min1.rst -r polyAT_wat_min2.rst
Input files:
polyAT_wat_min2.in,
polyAT_wat.prmtop,
polyAT_wat_min1.rst
Output files:
polyAT_wat_min2.out,
polyAT_wat_min2.rst
This takes approximately 854.0 s (14 min 14 s) on a 3 GHz Pentium 4.
Lets take a quick look at the structure, in order to ensure it still sensible, before proceeding:
$AMBERHOME/exe/ambpdb -p polyAT_wat.prmtop < polyAT_wat.inpcrd > polyAT_wat.pdb
$AMBERHOME/exe/ambpdb -p polyAT_wat.prmtop < polyAT_wat_min2.rst > polyAT_wat_min2.pdb
We shall use VMD to visualise the 2 structures (polyAT_wat.pdb, polyAT_wat_min2.pdb) and compare the structure of the DNA at the beginning and after our minimisation:
vmd
File -> New Molecule (load the first pdb polyAT_wat.pdb)
Lets make the DNA stand out from the water better:
Graphics -> Representations
Hit "Create Rep" This will create a second representation of our molecule
In the "Selected Atoms" box type: all not water
Then under "Drawing Method" select "Ribbons"
Next select the representation that still has the style "Lines"
And change the drawing method for this to "Points"
And there you have it, a much clearer representation of our initial DNA structure in the truncated octahedral water box:
Now, lets compare the starting structure and the minimised structure. This time we will instruct VMD not to draw the solvent. Quickest way to do this is to close VMD (File->Quit), and then re-load it. We will specify the first pdb file on the command line so we only need to load the second one manually:
vmd polyAT_wat.pdb
File -> New Molecule Load the polyAT_wat_min2.pdb file
Turn off the water to make it easier to compare the DNA's:
Select Graphics->Representations
In the Selected Atoms box enter: all not water
Click "Apply"
Now under the Selected Molecule drop down box at the top of the Graphical Representations window select molecule 0: polyAT_wat.pdb and repeat the process above.
You can see that the structure has changed slightly but there are no major problems so it would appear that the minimisation was successful. Note: The non-bonded blue dots are the sodium ions.
Now that we have successfully minimised our system the next stage in our equilibration protocol is allow our system to heat up from 0 K to 300 K. In order to ensure this happens without any wild fluctuations in our solute we will use a weak restraint, as we did in stage 1 of the minimisation, on the solute (DNA) atoms. Prior to AMBER 8 the recommended method for maintaining temperature was to use the Berendsen thermostat (NTT=1). This method is not very efficient at ensuring the temperature remains even across the system and so one would typically have to use NMR restraints in order to ensure that the heating occurred very gradually over a timescale of about 20 ps. This was essential in order to avoid problems with hot solvent cold solute etc. AMBER 8 now supports the new Langevin temperature equilibration scheme (NTT=3) which is significantly better at maintaining and equalising the system temperature. As a result the use of NMR restraints is no longer necessary. In fact it may even be possible to simply start our system at 300 K but we shall remain cautious and start from 0 K.
Our final aim is to run production dynamics at constant temperature and pressure since this more closely resembles laboratory conditions. However, at low temperatures, as our system will be for the first few ps the calculation of pressure is very inaccurate and so using constant pressure periodic boundaries in this situation can lead to problems. Using constant pressure with restraints can also cause problems so initially we will run 20 ps of MD at constant volume. Once our system has equilibrated over approximately 20 ps we will then switch off the restraints and change to constant pressure before running a further 100 ps of equilibration at 300 K.
Since these simulations are significantly more computationally expensive than the in vacuo and implicit solvent simulations it is essential that we try to reduce the computational complexity as much as possible. One way to do this is to use triangulated water, that is water in which the angle between the hydrogens is kept fixed. One such model is the TIP3P water model with which we solvated our system. When using such a water model it is essential that the hydrogen atom motion of the water also be fixed since failure to do this can lead to very large inaccuracies in the calculation of the densities etc. Since the hydrogen atom motion in our DNA is unlikely to effect its large scale dynamics we can fix these hydrogens as well. One method of doing this is to use the SHAKE algorithm in which all bonds involving hydrogen are constrained (NTC=2). This method of removing hydrogen motion has the fortunate effect of removing the highest frequency oscillation in the system, that of the hydrogen vibrations. Since it is the highest frequency oscillation that determines the maximum size of our time step (1 fs for the in vacuo runs) by removing the hydrogen motion we can safely increase our time step to 2 fs without introducing any instability into our MD simulation. This has the effect of allowing us to cover a given amount of phase space in half as much time since we need only 50,000 steps to cover 100 ps in time as opposed to the 100,000 we required with a 1 fs time step.
Note, if your system is very inhomogeneous you may need to run short MD runs with smaller time steps, say 0.5 fs, before conducting long MD runs with a 2 fs time step. This approach allows the system to rapidly move out of an unfavourable geometry without wild oscillations in the energy.
So, stage 1, 20 ps of MD with weak positional restraints on the DNA. Here is the input file:
| polyAT_wat_md1.in |
polyA-polyT 10-mer: 20ps MD with res on DNA &cntrl imin = 0, irest = 0, ntx = 1, ntb = 1, cut = 10, ntr = 1, ntc = 2, ntf = 2, tempi = 0.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 10000, dt = 0.002 ntpr = 100, ntwx = 100, ntwr = 1000 / Keep DNA fixed with weak restraints 10.0 RES 1 20 END END |
The meaning of each of the terms are as follows:
We run this using the following command. Note, we use the restrt file from the second stage of our minimisation since this contains the final minimised structure. We also use this as the reference structure with which to restrain the DNA to:
$AMBERHOME/exe/sander -O -i polyAT_wat_md1.in -o polyAT_wat_md1.out -p polyAT_wat.prmtop -c polyAT_wat_min2.rst -r polyAT_wat_md1.rst -x polyAT_wat_md1.mdcrd -ref polyAT_wat_min2.rst
Input files:
polyAT_wat_md1.in,
polyAT_wat.prmtop,
polyAT_wat_min2.rst
Output files: polyAT_wat_md1.out,
polyAT_wat_md1.rst,(polyAT_wat_md1.mdcrd.gz
[8.2 MB])
Note: This simulation may take a long time. In parallel on four 1.3GHz processors of a SGI ALTIX it takes approximately 682 s (11 mins 22 s). You can run this yourself if you want but expect it to take at least an hour on a single processor Pentium 4 machine. Alternatively you can just use the output files provided in the links above. For information on running sander in parallel please refer to the users manual.
Now that we have successfully heated our system at constant volume with weak restraints on the DNA the next stage is to switch to using constant pressure, so that the density of our water can relax. At the same time, now we are at 300 K we can safely remove the restraints on our DNA.
We will run this equilibration for 100 ps to give our system plenty of time to relax.
| polyAT_wat_md2.in |
polyA-polyT 10-mer: 100ps MD &cntrl imin = 0, irest = 1, ntx = 7, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 10, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 50000, dt = 0.002, ntpr = 100, ntwx = 100, ntwr = 1000 / |
The meaning of each of the terms are as follows:
We run this using the following command. Note, we use the restrt file from the first stage of our MD since this contains the final MD frame, velocities and box information.
$AMBERHOME/exe/sander -O -i polyAT_wat_md2.in -o polyAT_wat_md2.out -p polyAT_wat.prmtop -c polyAT_wat_md1.rst -r polyAT_wat_md2.rst -x polyAT_wat_md2.mdcrd
Input files:
polyAT_wat_md2.in,
polyAT_wat.prmtop,
polyAT_wat_md1.rst
Output files: polyAT_wat_md2.out,
polyAT_wat_md2.rst,(polyAT_wat_md2.mdcrd.gz
[42.8 MB])
Note: This simulation may take considerable amount of time to run (about 5 times more than the 20 ps run above), in parallel on four 1.3 GHz processors of a SGI ALTIX it takes approximately 3670 s (1 hr 1 min 10 s). You can run this yourself if you want but expect it to take several hours on a single processor Pentium 4 machine, alternatively you can just use the output files provided in the links above.
Now that we have equilibrated our system it is essential that we check that equilibrium has successfully been achieved before moving on to running any production MD simulations in which we would hope to learn something new about our molecule.
There are a number of system properties that we should monitor to check the quality of our equilibrium. These include:
The following steps will illustrate how to extract this data and plot it and also discuss what we expect to see for a successful equilibration.
In this section we will look at the system properties that can be extracted from the data written to the mdout files during our 120ps of MD simulation (polyAT_wat_md1.out, polyAT_wat_md2.out). This includes the potential, kinetic and total energies, the temperature, pressure, volume and density.
Lets start by extracting the various properties from our two output files. For this we will use the process_mdout perl script that was introduced earlier. This script can be given any number of files to process and it will append the results to single summary files:
mkdir analysis
cd analysis
process_mdout.perl ../polyAT_wat_md1.out ../polyAT_wat_md2.out
This should give us a series of summary files in our analysis directory (polyAT_wat_md1+2.summary.tar.gz). Lets plot some of these to see if we have successfully reached an equilibrium. First of all the energies, potential, kinetic and total (= potential + kinetic):
xmgr summary.EPTOT summary.EKTOT summary.ETOT
Here the red line, which is positive, is our kinetic energy, the black line is the potential energy (negative) and the green line is our total energy. The important points to note with this plot is that all of the energies increased during the first few ps, corresponding to our heating from 0 K to 300 K. The kinetic energy then remained constant for the remainder of the simulation implying that our temperature thermostat, which acts on the kinetic energy, was working correctly. The potential energy, and consequently the total energy since total energy = potential energy + kinetic energy, initially increased and then plateaued during the constant volume stage (0 to 20 ps) before decreasing as our system relaxed when we switched off the DNA restraints and moved to constant pressure (20 to 40 ps). The potential energy then levelled off and remained constant for the remainder of our simulation (40 to 120 ps) indicating that the relaxation was completed and that we reached an equilibrium.
Next we plot the temperature:
xmgr summary.TEMP
Here we see the expected behaviour, our system temperature started at 0 K and then increased to 300 K over a period of about 5 ps. The temperature then remained more or less constant for the remainder of the simulation indicating the use of Langevin dynamics for temperature regulation was successful.
Next the pressure:
xmgr summary.PRES
The pressure plot is slightly different than the previous plots. For the first 20 ps the pressure is zero. This is to be expected since we were running a constant volume simulation in which the pressure wasn't evaluated. At 20 ps we switched to constant pressure, allowing the volume of our box to change, at which point the pressure drops sharply becoming negative. The negative pressures correspond to a "force" acting to decrease the box size and the positive pressures to a "force" acting to increase the box size. The important point here is that while the pressure graph seems to show that the pressure fluctuated wildly during the simulation the mean pressure stabilised around 1 atm after about 50 ps of simulation. This is sufficient to indicate successful equilibration.
Next we will consider the volume. Unfortunately we cannot simply plot the summary.VOLUME file since the first 20 ps do not contain any volume information since during constant volume simulations the volume of the box is not written to the output file. So, we need to remove the first 101 lines of the file. You can use your favourite text editor to do this - I like vi:
vi summary.VOLUME
d100 (removes the current (first) line plus the next 100 lines)
<esc>:wq summary.VOLUME.modified
Here's the modified file (summary.VOLUME.modified).
xmgr summary.VOLUME.modified
Notice how the volume of our system (in angstroms3) initially decreases as our water box relaxes and reaches and equilibrated density, and thus volume. The smooth transitions in this plot followed by the oscillations about a mean value suggest our equilibration has been successful.
Finally, lets take a look at the density, we expect this to mirror the volume. We have the same problem with this summary file as we did with the volume since the density is not written to the output during constant volume simulations. Here is the modified file with the top 101 lines removed (summary.DENSITY.modified).
This is as expected. Our system has equilibrated at a density of approximately 1.04 g cm-3. This seems reasonable. The density of pure liquid water at 300 K is approximately 1.00 g cm-3 so adding a 10-mer of DNA and the associated charges has increased the density by around 4 %.
Now we have analysed our output files and seen that the main system properties suggest that our equilibration has been successful the next step is to consider the structures. Have they remained reasonable? One useful measure to consider is the root mean square deviation (RMSd) from the starting structure. We can use ptraj to calculate this for us as a function of time. In this situation we are interested in the backbone of our DNA so we shall consider just the main backbone atoms, P, O3*, O5* C3*, C4* and C5* (P, O3', O5'. C3', C4' and C5' in our prmtop). We will conduct a standard (not mass weighted) RMS fit to the first structure of our MD (the final structure from the minimisation). Note, the time is set to 0.2 this time around since although we wrote to the trajectory file every 100 steps the time step was 2 fs - hence 100 steps = 0.2 ps. Here is the input file for ptraj:
| polyAT_wat_calc_backbone_rms.in |
| trajin polyAT_wat_md1.mdcrd trajin polyAT_wat_md2.mdcrd rms first out polyAT_wat_backbone.rms @P,O3',O5',C3',C4',C5' time 0.2 |
Note: in early versions of amber it was often necessary to re-image molecules back into the primary box before calculating RMSd fits. This is no longer necessary since by default AMBER tracks the original molecule, even if it moves out of the central box during the course of the simulation. In early versions of AMBER you would see the parts of the molecule that moved out of the central box re-appear at the opposite edge. Calculating and RMSd of such a structure would result in huge (>20 angstrom) jumps in the results. This is no longer a problem. However, it you want to calculate certain parameters, such as densities you will need to re-image everything back into the central box. This is also necessary if you want movies of the trajectory to look "good" since otherwise it will look like your water molecules are flying off into space. In reality the periodic boundary approach results in them being re-introduced on the opposite side of the box, this is just not shown in the trajectory file. How to re-image a trajectory file is covered in the next stage of the tutorial.
Lets run ptraj to calculate the RMSd of the DNA backbone atoms (P, O3*, O5* C3*, C4* and C5*):
$AMBERHOME/exe/ptraj polyAT_wat.prmtop < polyAT_wat_calc_backbone_rms.in
Here is the output file (polyAT_wat_backbone.rms).
xmgr polyAT_wat_backbone.rms
Here we see that the RMSd of the DNA backbone atoms remained low for the first 20 ps, this was due to the restraints. Upon removing the restraints the RMSd shot up as the DNA relaxed within the solvent. After that our RMSd is fairly stable with no wild oscillations.
Lets take a look at the trajectory using VMD. First of all we will do this without re-imaging back to the primary box. We will then use ptraj to re-image our trajectory and then look at it again. So, first of all lets load both trajectory files into vmd:
Step 1 - load the prmtop file (polyAT_wat.prmtop)
remembering to select "parm7" Now - rewind to the first frame and press play:
There you have it, a water droplet flying through space. Notice what happens to the water as we move through the simulation (see below).
|
It appears like the water box is breaking down and the waters are boiling off into space. This would indeed be true if we were not using periodic boundaries. However, as mentioned above the atom coordinates written to the trajectory files are those of the original atom, even if it migrates outside of the original box and is replaced by one of its images. We can "fix" this by re-imaging our trajectory as we will now do.
Lets use ptraj to re-image our trajectory files so that only the atoms in the primary box are shown. Here is the ptraj input file:
| polyAT_wat_reimage.ptraj |
trajin polyAT_wat_md1.mdcrd trajin polyAT_wat_md2.mdcrd trajout polyAT_wat_md_reimaged.mdcrd center :1-20 image familiar go |
$AMBERHOME/exe/ptraj polyAT_wat.prmtop < polyAT_wat_reimage.ptraj
This will append the md2 trajectory file to the md1 file, change the trajectory such that the centre of geometry of residues 1 to 20 is placed at the origin and then image all of the coordinates outside of the primary box back to their images inside the box. The familiar keyword makes sure ptraj writes out the box as a familiar truncated octahedron as opposed to the default, triclinic imaging. Finally the converted trajectory will be written to (polyAT_wat_md_reimaged.mdcrd.gz [50.1 MB]).
You can now open this trajectory file in VMD where you will find our flying water droplet is much more intact. Have a play with the representations in vmd, remove the water and look at just the DNA etc. Observe how the dynamics of the DNA molecule change after the first 20ps when the restraints are switched off. Also, observe how much more stable our DNA chain is now than it was in the in vacuo simulations.
You have now covered the techniques necessary to create an initial input structure for a small DNA 10-mer, modify the pdb file to make it compatible with LEaP and then load this into xleap. We have covered how to charge neutralise this structure, solvate it in water and then create the necessary input files for running MD using sander. What to do when you have residues that are not in the default AMBER templates will be covered in later tutorials.
You should now know how to use these input files to run simulations in vacuo, in implict solvent and in explicit solvent. The concept of minimisation and equilibration has been covered and finally you now know how to perform simple analysis of the results and view trajectory files.
The final stage of this tutorial, which is optional, is to use these newly learned techniques in a practical example. One might ask the question: What happens if I start the simulations using the structure for A-DNA? In the last stage of this tutorial we shall try just that.
(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright Ross Walker 2004
