#include <math.h>  
#include <stdio.h> 
#include <stdlib.h>
#define n 100
#define k 8.61738569226e-5 // eV/K

double M=1;
double Omega=1;
double Inter=0;
double A=0.025;
double Temp=294.15;
double Delta=0.03;


double tau(double x1[],double x2[],double x3[])
{
	int i,j;
	double rij=0,interaction=0,uij=0;
	for (i=1;i<=n;i++)
	{
		for (j=1;j<=n;j++)
		{
			rij=(pow((x1[i]-x1[j]),2)+pow((x2[i]-x2[j]),2)+pow((x3[i]-x3[j]),2));
			interaction+=Inter*0.5*rij*exp(-sqrt(rij)/A);
		}
		uij+=M*Omega*Omega*0.5*(pow(x1[i],2)+pow(x2[i],2)+pow(x3[i],2))-interaction;
	}
	return uij;
}

inline double embt(double t)
{
	static double beta=1/(k*Temp);
	return (exp(-beta*t));
}

void main()
{
	FILE *stream1;
	FILE *stream;
	char datei1[30];
	printf("\nEingabedatei: ");
	scanf("%s",datei1);
	if ((stream=fopen(datei1,"r"))==NULL)
	{
		printf("Kann Datei nicht öffnen!!!");
	}
	fscanf(stream,"%lf",&M);printf("m: %1.2lf ",M);
	fscanf(stream,"%lf",&Omega);printf("omega: %1.2lf ",Omega);
	fscanf(stream,"%lf",&Inter);printf("X: %1.2lf ",Inter);
	fscanf(stream,"%lf",&A);printf("a: %1.4lf ",A);
	fscanf(stream,"%lf",&Temp);printf("T: %4.2lf ",Temp);
	fscanf(stream,"%lf",&Delta);printf("Delta: %le\n",Delta);
	int nz=1,j,count=0,wo=0;
	double t=0,tt=0;
	double s=0;
	double x1[n+1]={0};
	double xt1[n+1]={0};
	double x2[n+1]={0};
	double xt2[n+1]={0};
	double x3[n+1]={0};
	double xt3[n+1]={0};
	srand(10000);
	for (j=1;j<=n;j++) {
		x1[j]=0;
		x2[j]=0;
		x3[j]=0;}
	t=tau(x1,x2,x3);
	double emb=embt(t);
	double energie=0;
	char datei[50];
	printf("Wie oft: "); scanf("%d",&wo);
	printf("Ausgabedateiname: "); scanf("%s",datei);
	if ((stream=fopen(datei,"w"))==NULL)                 
	{                                            
		printf("Kann Datei nicht öffnen!!!");
		exit(1);                             
	}                                            
	static double beta=1/(k*Temp);
	nz=1;
	do {
		for (j=1;j<=n;j++) {
			xt1[j]=x1[j]+Delta*(double(rand())/(RAND_MAX/2)-1);
			xt2[j]=x2[j]+Delta*(double(rand())/(RAND_MAX/2)-1);
			xt3[j]=x3[j]+Delta*(double(rand())/(RAND_MAX/2)-1);}
		fprintf(stream,"%lf %lf\n",xt1[1],xt2[1]);
		tt=tau(xt1,xt2,xt3);
		//printf("%le\n",tt);
		if ((double(rand())/(RAND_MAX)) < (embt(tt)/emb))
			{
				for (j=1;j<=n;j++) {x1[j]=xt1[j];x2[j]=xt2[j];x3[j]=xt3[j];emb=embt(tt);t=tt;}
			}
		else count++;
		s+=t;
		if (nz>1) {
			energie=+1.5*n/beta+s/nz;
			}
		nz++;
	} while (nz <= wo);
	printf("Ergebnis: %1.5le bei n=%d\n",energie,(nz-1));
	printf("%d von %d Werten gleich.\n",count,(nz-1));
}