C Program for ODHO

#include 
#include 
#define NSTEP 500

/* An example demonstrating the Verlet algorthm for Molecular Dynamics
   in one dimension in a harmonic well

             |
     *       |       *      PE = 1/2 kx^2
     *       |       *
      *      |      *       KE = 1/2 mv^2
      *      |      *
       *     |     *        a  = F/m = -kx/m
        *    |   **
         **  | **
   ---------***--------
        
   Author:  G.C. Rutledge     Department of Chemical Engineering
	                      MIT
	                      March 1993
			      
*/

void odho(float dt, int flag)
{
  int i,t=0;

  float pe,ke,te,te0,dtsq,a,x_new,x_old;
  float mv,force;
  float k = 1.0, m = 1.0;
  float x = 0.;
  float v = 1.;
  float xmax,mvmax,fmax;

/*  initial state */
  pe = 0.5 * k * x * x;
  ke = 0.5 * m * v * v;
  te = pe + ke;
  te0 = te;
  mv = m*v;
  force = -k*x;
  a  = 0.;
  xmax = sqrt(2.*te0/k);
  mvmax = sqrt(2.*te0*m);
  fmax = k*xmax;

/* update position to get Verlet started */
  x_old = x;
  x     = x + dt * v;

  dtsq = dt*dt;

  for(t=1;t<=NSTEP;t++) {
    a     = -k * x / m;
    x_new = 2.0 * x - x_old + dtsq * a;
    v     = (x_new - x_old)/2.0/dt;
    
    pe    = 0.5 * k * x * x;
    ke    = 0.5 * m * v * v;
    te    = pe + ke;
    mv = m * v;
    force = -k * x;

    x_old = x;
    x     = x_new;
  }
}