CLAMPS Methods
This page describes the algorithms used in the various simulation methods with some references, and a description of the output of each method. Note that the different methods keep separate averages of such things as the potential energy. For example, the potential energy for molecular dynamics appears in the averages under MD_PE and for Monte Carlo as MC_PE. The equality of these two quantities provides a powerful check for the program.
MD Molecular dynamics
The routine MOLDY(ISTEPS,TBATH) executes ISTEPS integration- steps starting from the position XOLD=X(O) and velocity VELO=V(0). Upon each entry to MOLDY the next position XNOW=X(1) is obtained from the routine DSTART by using a second order Taylor series expansion about the point XOLD. Then MOLDY iterates the following sequence ISTEPS times.
X(n+l)=2*X(n)-X(n-l)+TAU**2*F(n)/M
V(n)=(X(n+l)-X(n-1))/2*TAU
On exit XOLD and VELO contain positions and velocities advanced a time ISTEPS*TAU
MD computes the following scalars:
MD_EN Total energy/particle
MD_POT Potential energy/particle
mom**2 Total momentum squared (should be zero)
VBONDD Total polymer potential energy
Note that the ENERGY should be constant and that MOM**2 should remain very small since initially the momentum is set to zero.
If TBATH is nonzero there are two additional forces felt by each particle 1) a viscous drag equal to -M*V/(TBATH*XINV) where XINV is the inverse friction coefficient for this particle 2) a Gaussian random force, mean zero and variance equal to 2*M/(BETA*TAU*TBATH*XINV). This force is such as to maintain energy balance at a temperature equal to 1/BETA. The integration is performed with the Verlet algorithm. This is equivalent to simulating the Langevin equation, but is also useful in bringing a molecular dynamics trajectory into a desired temperature regime. The kinetic energy will relax with a time constant equal to TBATH*XINV. Hence this should be chosen to be somewhat less than the duration of the dynamical trajectory.
The algorithm used is due to L. Verlet , Phys. Rev. 159,98(1967).
MC : Metropolis Monte Carlo
A call to MONTE(ISTEPS) causes each particle to be moved sequentially ISTEPS times in the following manner. Say particle i is being moved from its present position XOLD(i). A trial position XNEW is sampled uniformly inside a cube (square in 2D) of side TAU and center XOLD (i). The new position is accepted with probability p where:
P = MIN (1,EXP (-BETA*(V(XNEW)-i(XOLD))))
where V(x) is the total potential energy felt by the ith particle at position x.
MC computes the following scalars:
MC_ACC acceptance ratio/move
MC_PE total potential energy/particle
VBONDM total polymer potential energy (IFPOLY.NE.0)
Note that TAU should be adjusted so that MC_ACC is between 0.1 and 0.5.
The algorithm used is : M. Metropolis, A.W. Rosenbluth, M. N. Rosenbluth, A. N. Teller and E.Teller, J. Chem. Phys. 21,1087 (1953).
PR : Polymer reptation
The algorithm used is the following (). A call to SNAKE(ISTEPS) causes NCHAINS*ISTEPS elementary trial moves. In such a move the chain number and end of the chain to consider as head are picked at random. Suppose the polymer chain has coordinates xl,x2 ... xn and xn is the head. A new trial head position xl=xn + dx is chosen, where dx is a vector, randomly chosen (by the routine RNBOND using a rejection technique) from the probability distribution
EXP(-BETA*VBOND(dx))*normalization
VBOND(dx) is the polymer potential energy of a bond stretched a distance dx (see INPUT 8 ). The trial chain then has coordinates x2,x3,...,xn,x' This configuration is accepted with probability p where:
p=MIN(l,EXP(-BETA*(VP(x')-VP(xl))))
where VPi(x) non-bonding potential only (e.g. Lennard-Jones of particle i at position x. if the move is accepted then xl=x' , particle I becomes the new head and particle 2 the new tail. Note that the head and tail are adjacent mod(n) but are not otherwise restricted. In this way we avoid shuffling the particle coordinates after every move.
Polymer reptation is a much quicker way of moving through the configuration space of a polymer system since a single particle move cause the chain to move significantly.
PR computes the following scalars:
PR_ACC acceptance ratio/step
ENDTOEND end-to-end distance for a polymer of type i one entry/POLYMER.distance computed with periodic boundary conditions
BD : Brownian Dynamics
This algorithm is also known as 'smart Monte Carlo' since with a large time step TAU one may sample the Boltzmann distribution quicker than with the standard Metropolis Monte Carlo. However, if TAU is chosen to be small then the algorithm simulates the diffusion through configuration space of a system where the viscous forces of the surrounding solvent bath dominate the inertial terms. This process is described by the Smoluchowski equation:
Df/Dt=gradi(DCONSTi[gradi(f)-BETA*f*Fi])
where f(R,t) is the time dependent probability distribution in phase space, Fi is the force on the ith particle, and a sum over particles i and dimensions is implied on the right hand side. Here DCONST = XINV/(BETA*MASS) is the diffusion constant, XINV is the reciprocal friction coefficient.
A call to BROWN(ISTEPS) results in moving every particle sequentially to a trial position ISTEPS times. Say particle i has the current position XOLD. Then a trial position XNEW=XOLD+DX is chosen where DX is a normally distributed random vector with mean BETA*DCONSTi*Fi(XOLD)*TAU and variance 2*NDIM*DCONSTi*TAU. This move is accepted with probability p where:
p=MIN(l,EXP(-BETA*(delta))
delta=Vi(XNEW)-vi(XOLD)
+25*(Fi(XOLD)+Fi(XNEW))*(DCONSTi*P,ETA*TAU**(Fi(XOLD)-Fi(XNEW))-2*DX))))
This is the Metropolis acceptance probability generalized to a non-uniform transition probability. As a consequence, the probability distribution will converge to the Boltzmann distribution no matter how large TAU is. BD computes the following scalars:
APTIME actual diffusion/theoretical diffusion
BD_VBND total polymer potential energy (IFPLY.NE.0)
BD_ACC acceptance ratio per step
BD_PE potential energy per particle
Because some trial moves are rejected, the diffusion constant of the center of mass, which one can show from the Smoluchoskwi equation to be equal to 2*NDIM*DCONST*time/NATOMS will decrease. APTIME gives the relative decrease of this diffusion. It is reasonable to assume all other motions will be slowed down by an equal amount. Hence the difference in time between two successive configurations is closer to APTIME*TAU then to TAU and all times should be scaled by the factor APTIME. Tests on polymeric systems show that if one chooses TAU so that ACCEPT is greater than 0.95 the truncation errors in most (time rescaled) correlation function are smaller than statistical errors.
Note that for equivalent moves, the time step of Metropolis Monte Carlo will be related to that of Brownian Dynamics by:
TAU_BD=DCONST*TAU_MC**2/24
The algorithm used here is the single particle version of that discussed in P.J. Rossky, J.D.Doll and H.L.Friedman, J.Chem.Phys. 69 , 4628 (1978).
Keywords | CLAMPS homepage | Examples
D. Ceperley 8/29/98