Andersen Thermostat

#include <stdio.h>
#include <stdlib.h>


double Ep,Ek,E;
double x=1.0,v=1.0,f=0.0;
double dt=0.005,k=1.0,m=1.0;

   /*calculate the force and potential energy of the system.*/
void force()
{
	Ep=0.5*k*x*x;
	f=-k*x;
}
//calculate the actual motion of the oscillator.
void md()
{
	v+=(f/m)*dt*0.5;
	x+=v*dt;
	force();
	v+=(f/m)*dt*0.5;
	Ek=0.5*m*v*v;
	E=Ek+Ep;
}
int main(void)
{
	FILE* fp;
	fp=fopen("erkka2.txt", "w"); //write output to "erkka2.txt"
	int step;
	for(step=0;step<10000;step++)
	{
		md();
		if(step %100 ==0)
			fprintf(fp,"%.8lf %.8lf\n",x,v);
	}
	fclose(fp);
}