Computational Chemistry/Molecular dynamics

Previous chapter - Molecular mechanics

Introduction to Molecular dynamics

You will remember from your spectroscopy that vibrations, and that includes internal rotations, are quantized. In MD we use a classical force field and ignore this quantization putting in classical energy according to Newton's Laws. This might at first sight appear to be totally unsound in view of all we have been taught about quantum theory but in fact what it means is we are replacing a summation with a Monte Carlo integration. Only in very special cases is there a dramatic difference. Situations with strong force constants and light particles i.e. protons are likely to show this. If, and this is a big if, the trajectories of the particles are chaotic all the phase space of the molecular vibrations and in particular internal rotational conformers, will be covered. Only if this has occurred do we have a temperature. (A temperature is a single number which allows one to reproduce the occupations of the states of the system using Boltzmann statistics. If there is no Boltzmann distribution, there is no temperature.)

One focus of research here has been the development of methods which allow a reliable temperature to be created for gas phase molecular systems. (Condensed phases with plentiful collisions are self randomizing and the Boltzmann distribution is rapidly established.) Simply apply Newtons laws of motion. In one dimensional form they are:

$f=ma$

$v=u+at$

$s=ut+{\frac {1}{2}}at^{2}$

s is a new position. Force is $-{\frac {dV_{s}}{dx}}$ . $t+\delta t$ then run for n timesteps.

Verlet algorithm

$s(t+\delta t)=s(t)+v(t)\delta t+{\frac {1}{2}}a(t)\delta t^{2}+....$

$s(t-\delta t)=s(t)-v(t)\delta t+{\frac {1}{2}}a(t)\delta t^{2}-....$

add these together then :

$s(t+\delta t)=2s(t)-s(t-\delta t)+a(t)\delta t^{2}-....$

These are applied to the 3N-6 dimensional hypersurface of a single molecule, or to all free coordinates of an infinitely repeating system.

Single molecule dynamics

Applications in liquid crystals and biological chemistry often involve dynamics on a single molecule in an infinite vacuum. There are considerable technical problems here both in establishing a temperature and in ensuring the 3N-6 dimensional phase space is covered with populations according to the Boltzmann distribution. In practice this can only be assured in special circumstances and for proteins there is the continual phenomenon of sticking in a given neighbourhood of conformations. In this case the timescales of the physical processes of changing between energetically equivalent conformations may be longer than the simulation. Putting molecules in a box with inert gas bumping into it is one way of trying for conformational equilibrium.

Table - Timescales of molecular processes


Secs	name	Process

10^-1		Surface reorganizations
10^-3	milliseconds	Protein Folding
10^-8	10 nanoseconds	NMR
10^-9	to 10⁺³	Brownian dynamics
10^-9	to 10⁰	Lattice Monte-Carlo
10^-12	picoseconds	Fast vibrations
10^-15	femtoseconds	Electronic excitation
		(molecular dynamics)

A good simulation regime for a teaching exercise is 10 femtoseconds with C-H bond vibrations frozen out. Research requires considerably more care and much longer runs! It would not be a good idea to investigate protein folding using a conventional MD program and a 0.5 femtosecond time step.

The finite box size

A finite box is also a problem, but a minimum image convention allows for infinitely repeating boxes in the 3 Cartesian directions. Software for infinite boxes and the AMBER and GROMOS force field exists within CCP5 (The DLPOLY package). There are problems with electrostatic interactions. The Coulomb Interaction goes as inverse $r$ to the power one. This dies away very slowly over several images of the box before dying away to insignificance. This must be handled by a complicated procedure known as an Ewald Sum.

Liquid simulations often calculate G-of-r, not to be confused with the quantum theory matrix G-of-R. This is the radial probability density, which is experimentally accessible via neutron diffraction. It can therefore act as simulation quality check on whether the right liquid solvation shell structure is being generated.

A good description of modern simulation techniques is in Allen and Tildesley.

Problems

i) A Gas phase liquid crystal molecule

Make a liquid crystal prototype by adding a (CH2)n tail onto phenol's acidic proton. Run a dynamic simulation at various temperatures to see the side chain move about.

2) Histogramming

Make hexamine-diamine. Use ANALYZE to monitor the distance between the N-atoms. Use DYNAMIC and MonSt to create a histogram on this N to N distance. Repeat this with a flexible dihedral angle in a molecule of interest to you.

Bibliography

M. P. Allen, and D. J. Tildesley, Computer Simulation of Liquids,(Oxford University Press, Oxford,1987).

Next Chapter - Molecular quantum mechanics