Polymer Visualization: Simulated Annealing and Molecular Dynamics

Roger Edberg and Drew Whitehouse
Supercomputer Facility Australian National University

Visualization techniques are used to follow the progress of two types of polymer calculations that are performed on the ANU's VP2200 supercomputer. First, we look at a simulated annealing calculation that unravels a random walk on a tetrahedral lattice and produces a realistic polymer conformation with excluded volume. Second, we take the resulting polymer and simulate its motion by molecular dynamics.

Random walk conformations of any length are quite easy to generate on the computer. On the other hand, conformations with excluded volume can be difficult to generate for even a small molecule of say 20 atoms (except for a few trivial conformations such as all-trans.)

Here, we start with a random walk conformation on a tetrahedral lattice, and proceed to optimize it, i.e. to remove all of the excluded volume interactions, by a simulated annealing [1] procedure. The function to be optimized is the 'cost function':

C = n0 + n1 + n2

where n0 is the number of self-intersections, n1 is the number of nearest- neighbor contacts, and n2 is the number of next-nearest-neighbor contacts. The optimization algorithm is:

Steps 2-5 are iterated until a conformation with C=0 is obtained. Note that the quantity T plays the role that temperature does in the Metropolis Monte Carlo algorithm; here it controls the acceptance of moves that result in increased energy. For our definition of the cost function, the algorithm should converge for T<1.0.

The annealing calculation shown on film was performed for a 200 atom molecule at T=0.8. Colors indicate the localization of the cost function within the molecule, i.e. excluded volume contacts. The unravelling process is not a physically realistic one because portions of the molecule may pass through one another as bond vectors are changed. Some additional work is necessary to improve the convergence of the algorithm; the cost function decreases quickly in the initial stages of the calculation, while in latter stages it takes much time to unravel the last few 'knots' in the molecule.

Once we have the unravelled conformation we can use it as a starting configuration for a molecular dynamics simulation. We choose a model with harmonic atom-atom bonds and Lennard-Jones (12-6) non-bonded interactions for atoms more than four apart on the chain. The velocity version of the Verlet algorithm is used to integrate the equations of motion [2]. Atomic velocities were randomized with an initial kinetic energy per particle of 3.0 (reduced units; initial temperature=2.0), and the system was allowed to evolve isoenergetically.

The film shows the molecule relaxing once the lattice geometry restriction is removed. The Lennard-Jones non-bonded potential causes the molecule to relax into a globular conformation. The mean squared radius of gyration decreases from its initial value of 26 to about 5.

In future work we plan to apply similar visualization techniques to Monte Carlo and molecular dynamics simulations of polymeric liquids and mixtures.

Appendix:

Lattice Parameters:

atomic diameter (D) = 3.923 Angstroms bond angle = 109.47 degrees bond length = 0.39D

MD parameters (reduced units):

timestep = 0.002 equilibrium bond length = 0.39 bond force constant = 100.0

References:

[1] S. Kirkpatrick, C.D. Gelatt, and M.P. Vecchi, "Optimization By Simulated Annealing", Science 220, 671 (1983)

[2] M.P. Allen and D. Tildesley, "Computer Simulation of Liquids", Clarendon Press, Oxford (1987)