From: Nick Schafer (npschafer_at_wisc.edu)
Date: Tue Jul 10 2007 - 11:01:36 CDT
Dear NAMDers -
My name is Nick Schafer and I am an undergraduate at the University of
Wisconsin, Madison. I'm interested in numerical methods for the phase
space propagation of molecular dynamics. In order to get a feel for how
the Verlet integrator stacks up against other simple methods, I've
implemented both Forward Euler and Heun's method in NAMD-Lite. As far
as I can tell, my Heun's method implementation is working. However, I'm
not very sure that I'm handling the force evaluations correctly because
the Verlet integrator that I used as a model uses only one force
evaluation, which is at the current positions. Heun's method requires a
force evaluation at a "predicted" state of the system, which makes it
just slightly more complicated. Anyway, here is some pseudo code that I
was trying to follow during my implementation, followed by my code. If
someone who knows exactly how the force evaluations in NAMD-lite work
could take a look and tell me whether they think I am doing it right, it
would be greatly appreciated.
Pseudo code:
% first, we calculate the current acceleration
current_acceleration = function(current_velocity, etc.)
% we use this to do a poor job of estimating future velocity
future_velocity = current_velocity + new_acceleration * timestep;
% use (poorly predicted) future velocity to calc future acceleration
future_acceleration = function(future_velocity, etc.)
% find the average acceleration
acceleration = 0.5 * (current_acceleration + future_acceleration)
% use the average acceleration
new_velocity = current_velocity + acceleration * timestep;
Inside runge.c, equivalent to verlet.c :
int step_compute_runge(Step *s, StepSystem *sys, int32 numsteps)
{
const double dt = s->param->timestep;
const double half_dt = 0.5 * dt;
double konst;
double pe;
void *f_obj = s->param->force_object;
int32 (*f_eval)(void *, double *, MD_Dvec *, MD_Dvec *)
= s->param->force_compute;
MD_Dvec *f = sys->force;
MD_Dvec *vel = sys->vel;
MD_Dvec *pos = sys->pos;
const double *scal_inv_mass = s->scal_inv_mass;
const int32 natoms = s->param->natoms;
int32 i;
MD_Dvec *pred_pos;
pred_pos = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec));
MD_Dvec *pred_vel;
pred_vel = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec));
MD_Dvec *pred_force;
pred_force = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec));
double temppe = 0;
/* perform two step Runge Kutta integration of system */
/* compute force at current positions */
if (f_eval(f_obj, &pe, f, pos)) return STEP_FAIL;
/* predict velocities, store in temporary array */
for (i = 0; i < natoms; i++) {
konst = dt * scal_inv_mass[i];
pred_vel[i].x = vel[i].x + f[i].x * konst;
pred_vel[i].y = vel[i].y + f[i].y * konst;
pred_vel[i].z = vel[i].z + f[i].z * konst;
}
/* predict new positions, store in temporary array */
for (i = 0; i < natoms; i++) {
konst = dt * scal_inv_mass[i];
pred_pos[i].x = pos[i].x + pred_vel[i].x * dt;
pred_pos[i].y = pos[i].y + pred_vel[i].y * dt;
pred_pos[i].z = pos[i].z + pred_vel[i].z * dt;
}
/* finished computing rough approximation of position at the new
time step */
/* copy old force array into pred_force array */
for (i = 0; i < natoms; i++) {
pred_force[i].x = f[i].x;
pred_force[i].y = f[i].y;
pred_force[i].z = f[i].z;
}
/* compute NEW force array at updated positions, saving current pe */
temppe = pe;
if (f_eval(f_obj, &pe, f, pred_pos)) return STEP_FAIL;
pe = temppe;
/* update velocities using average of two force arrays */
for (i = 0; i < natoms; i++) {
konst = half_dt * scal_inv_mass[i];
vel[i].x += konst * (f[i].x+pred_force[i].x);
vel[i].y += konst * (f[i].y+pred_force[i].y);
vel[i].z += konst * (f[i].z+pred_force[i].z);
}
/* update positions using new velocities */
for (i = 0; i < natoms; i++) {
pos[i].x += dt * vel[i].x;
pos[i].y += dt * vel[i].y;
pos[i].z += dt * vel[i].z;
}
/* free pred_pos, pred_vel and pred_force */
free(pred_pos);
free(pred_vel);
free(pred_force);
return 0;
}
Thanks much!
Nick Schafer
http://sbel.wisc.edu/People/schafer/mainframeset.html
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:57 CST