2.調和振動子のMD計算プログラム

調和振動子の運動をVerlet法で数値的に解くプログラムを以下に示す（VC++）。

#include "stdafx.h"
#include "stdlib.h"

int _tmain(int argc, _TCHAR* argv[])
{
  double x;       // current position
  double xb;      // position at the previous step
  double f;       // force
  double v;       // velocity components
  double m;       // mass
  double C;       // spring constant of harmonic oscillator
  double tmp;     // temporal memory
  double t, dt;   // duration, time step
  double KE;      // kinetic energy
  double PE;      // potential energy
  FILE *fp;

  fopen_s(&fp,"harmonic.csv", "w");
  if(fp == NULL){
    printf("ERROR: File open failed.\n");  exit(1);
  }
  fprintf(fp,"time,kinetic energy,potential energy,Hamiltonian\n");

  // Initialize
  dt = 0.001; m  = 0.5;  C  = 1.0; x  = 1.0; xb = 1.0;

  for(t=0.0; t<5.0; t+=dt){
    // force and potential energy calculation
    f = -C * x;
    PE = 0.5 * C * x * x;
    // the Verlet algorithm
    tmp = x;
    x = 2.0 * x - xb + dt * dt * f / m;
    xb = tmp;
    // kinetic energy calculation
    v = (x - xb) / dt;
    KE = 0.5 * m * v * v;
    // print time,kinetic energy, potentialenergy, and Hamiltonian
    fprintf(fp, "%f,%f,%f,%f\n",t, KE, PE, KE+PE);
  }
  fclose(fp);
  return 0;
}