From 792f14f95b47989f7f12df0ea70b54619be016ee Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Tue, 4 Apr 2006 10:09:32 +0000
Subject: [PATCH] checkin before cafeteria visit ;)

---
 moldyn.c | 200 +++++++++++++++++++++++++++++++++++++++++++++++++++----
 moldyn.h |  45 +++++++++++--
 posic.c  | 159 ++++---------------------------------------
 3 files changed, 240 insertions(+), 164 deletions(-)

diff --git a/moldyn.c b/moldyn.c
index aec87f6..f2ac637 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -5,17 +5,136 @@
  *
  */
 
-#include "moldyn.h"
-
+#define _GNU_SOURCE
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
 #include <math.h>
 
+#include "moldyn.h"
+
 #include "math/math.h"
 #include "init/init.h"
 #include "random/random.h"
 #include "visual/visual.h"
 
+int moldyn_usage(char **argv) {
+
+	printf("\n%s usage:\n\n",argv[0]);
+	printf("--- general options ---\n");
+	printf("-E <steps> <file> (log total energy)\n");
+	printf("-M <steps> <file> (log total momentum)\n");
+	printf("-D <steps> <file> (dump total information)\n");
+	printf("-S <steps> <filebase> (single save file)\n");
+	printf("--- physics options ---\n");
+	printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
+	printf("-t <timestep tau> [s] (%f)\n",MOLDYN_TAU);
+	printf("-R <runs> (%d)\n",MOLDYN_RUNS);
+	printf("\n");
+
+	return 0;
+}
+
+int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
+
+	int i;
+
+	memset(moldyn,0,sizeof(t_moldyn));
+
+	/* default values */
+	moldyn->t=MOLDYN_TEMP;
+	moldyn->tau=MOLDYN_TAU;
+	moldyn->time_steps=MOLDYN_RUNS;
+
+	/* parse argv */
+	for(i=1;i<argc;i++) {
+		if(argv[i][0]=='-') {
+			switch(argv[i][1]){
+				case 'E':
+					moldyn->ewrite=atoi(argv[++i]);
+					strncpy(moldyn->efb,argv[++i],64);
+					break;
+				case 'M':
+					moldyn->mwrite=atoi(argv[++i]);
+					strncpy(moldyn->mfb,argv[++i],64);
+					break;
+				case 'D':
+					moldyn->dwrite=atoi(argv[++i]);
+					strncpy(moldyn->dfb,argv[++i],64);
+					break;
+				case 'S':
+					moldyn->swrite=atoi(argv[++i]);
+					strncpy(moldyn->sfb,argv[++i],64);
+					break;
+				case 'T':
+					break;
+					moldyn->t=atof(argv[++i]);
+				case 't':
+					moldyn->tau=atof(argv[++i]);
+					break;
+				case 'R':
+					moldyn->time_steps=atoi(argv[++i]);
+					break;
+				default:
+					printf("unknown option %s\n",argv[i]);
+					moldyn_usage(argv);
+					return -1;
+			}
+		} else {
+			moldyn_usage(argv);
+			return -1;
+		}
+	}
+
+	return 0;
+}
+
+int moldyn_log_init(t_moldyn *moldyn) {
+
+	moldyn->lvstat=0;
+
+	if(moldyn->ewrite) {
+		moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
+		if(moldyn->efd<0) {
+			perror("[moldyn] efd open");
+			return moldyn->efd;
+		}
+		dprintf(moldyn->efd,"# moldyn total energy logfile\n");
+		moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
+	}
+
+	if(moldyn->mwrite) {
+		moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
+		if(moldyn->mfd<0) {
+			perror("[moldyn] mfd open");
+			return moldyn->mfd;
+		}
+		dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
+		moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
+	}
+
+	if(moldyn->swrite)
+		moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
+
+	if(moldyn->dwrite) {
+		moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
+                if(moldyn->dfd<0) {
+			perror("[moldyn] dfd open");
+			return moldyn->dfd;
+		}
+		write(moldyn->dfd,moldyn,sizeof(t_moldyn));
+		moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
+	}
+
+	if(moldyn->dwrite)
+		moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
+
+	return 0;
+}
 
 int create_lattice(unsigned char type,int element,double mass,double lc,
                    int a,int b,int c,t_atom **atom) {
@@ -73,7 +192,7 @@ int destroy_lattice(t_atom *atom) {
 	return 0;
 }
 
-int thermal_init(t_atom *atom,t_random *random,int count,double t) {
+int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
 
 	/*
 	 * - gaussian distribution of velocities
@@ -84,11 +203,14 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) {
 	int i;
 	double v,sigma;
 	t_3dvec p_total,delta;
+	t_atom *atom;
+
+	atom=moldyn->atom;
 
 	/* gaussian distribution of velocities */
 	v3_zero(&p_total);
 	for(i=0;i<count;i++) {
-		sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
+		sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
 		/* x direction */
 		v=sigma*rand_get_gauss(random);
 		atom[i].v.x=v;
@@ -111,15 +233,18 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) {
 	}
 
 	/* velocity scaling */
-	scale_velocity(atom,count,t);
+	scale_velocity(moldyn,count);
 
 	return 0;
 }
 
-int scale_velocity(t_atom *atom,int count,double t) {
+int scale_velocity(t_moldyn *moldyn,int count) {
 
 	int i;
 	double e,c;
+	t_atom *atom;
+
+	atom=moldyn->atom;
 
 	/*
 	 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
@@ -127,7 +252,7 @@ int scale_velocity(t_atom *atom,int count,double t) {
 	e=0.0;
 	for(i=0;i<count;i++)
 		e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
-	c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
+	c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*moldyn->t));
 	for(i=0;i<count;i++)
 		v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 
@@ -202,9 +327,23 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
 int moldyn_integrate(t_moldyn *moldyn) {
 
 	int i;
-	int write;
+	unsigned int e,m,s,d,v;
+	unsigned char lvstat;
+	t_3dvec p;
+
+	int fd;
+	char fb[128];
 
-	write=moldyn->write;
+	e=moldyn->ewrite;
+	m=moldyn->mwrite;
+	s=moldyn->swrite;
+	d=moldyn->dwrite;
+	v=moldyn->vwrite;
+
+	if(!(lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
+		printf("[moldyn] warning, lv system not initialized\n");
+		return -1;
+	}
 
 	/* calculate initial forces */
 	moldyn->force(moldyn);
@@ -213,11 +352,44 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		/* integration step */
 		moldyn->integrate(moldyn);
 
-		/* check for visualiziation */
-		if(!(i%write)) {
-			visual_atoms(moldyn->visual,i*moldyn->tau,
-			             moldyn->atom,moldyn->count);
-			printf("finished %d / %d\n",i,moldyn->time_steps);
+		/* check for log & visualiziation */
+		if(e) {
+			if(!(i%e))
+				dprintf(moldyn->efd,
+				        "%.15f %.45f\n",i*moldyn->tau,
+				        get_total_energy(moldyn));
+		}
+		if(m) {
+			if(!(i%m)) {
+				p=get_total_p(moldyn->atom,moldyn->count);
+				dprintf(moldyn->mfd,
+				        "%.15f %.45f\n",i*moldyn->tau,
+				        v3_norm(&p));
+			}
+		}
+		if(s) {
+			if(!(i%s)) {
+				snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
+				         moldyn->t,i*moldyn->tau);
+				fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
+				if(fd<0) perror("[moldyn] save fd open");
+				else {
+					write(fd,moldyn,sizeof(t_moldyn));
+					write(fd,moldyn->atom,
+					      moldyn->count*sizeof(t_atom));
+				}
+			}	
+		}
+		if(d) {
+			if(!(i%d))
+				write(moldyn->dfd,moldyn->atom,
+				      moldyn->count*sizeof(t_atom));
+
+		}
+		if(v) {
+			if(!(i%v))
+				visual_atoms(moldyn->visual,i*moldyn->tau,
+				             moldyn->atom,moldyn->count);
 		}
 	}
 
diff --git a/moldyn.h b/moldyn.h
index a2633ba..cbfc173 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -24,19 +24,39 @@ typedef struct s_atom {
 } t_atom;
 
 typedef struct s_moldyn {
+	/* atoms, amount, dimensions */
 	int count;
 	t_atom *atom;
+	t_3dvec dim;
+	/* potential, force & parameters */
 	double (*potential)(struct s_moldyn *moldyn);
-	void *pot_params;
 	int (*force)(struct s_moldyn *moldyn);
+	void *pot_params;
 	double cutoff;
 	double cutoff_square;
-	t_3dvec dim;
+	/* temperature */
+	double t;
+	/* integration of newtons equations */
 	int (*integrate)(struct s_moldyn *moldyn);
 	int time_steps;
 	double tau;
+	/* logging & visualization */
+	unsigned char lvstat;
+	unsigned int ewrite;
+	char efb[64];
+	int efd;
+	unsigned int mwrite;
+	char mfb[64];
+	int mfd;
+	unsigned int swrite;
+	char sfb[64];
+	int sfd;
+	unsigned int dwrite;
+	char dfb[64];
+	int dfd;
+	unsigned int vwrite;
 	void *visual;
-	int write;
+	/* moldyn general status */
 	unsigned char status;
 } t_moldyn;
 
@@ -57,9 +77,20 @@ typedef struct s_lj_params {
 
 /* general defines */
 
+#define MOLDYN_LVSTAT_TOTAL_E		0x01
+#define MOLDYN_LVSTAT_TOTAL_M		0x02
+#define MOLDYN_LVSTAT_SAVE		0x04
+#define MOLDYN_LVSTAT_DUMP		0x08
+#define MOLDYN_LVSTAT_VISUAL		0x10
+#define MOLDYN_LVSTAT_INITIALIZED	0x80
+
 #define MOLDYN_STAT_POTENTIAL		0x01
 #define MOLDYN_STAT_FORCE		0x02
 
+#define MOLDYN_TEMP			273.0
+#define MOLDYN_TAU			1.0e-15
+#define MOLDYN_RUNS			1000000
+
 /* phsical values */
 
 #define K_BOLTZMANN		1.3807e-27			/* Nm/K */
@@ -79,11 +110,15 @@ typedef struct s_lj_params {
 
 /* function prototypes */
 
+int moldyn_usage(char **argv);
+int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv);
+int moldyn_log_init(t_moldyn *moldyn);
+
 int create_lattice(unsigned char type,int element,double mass,double lc,
                    int a,int b,int c,t_atom **atom);
 int destroy_lattice(t_atom *atom);
-int thermal_init(t_atom *atom,t_random *random,int count,double t);
-int scale_velocity(t_atom *atom,int count,double t);
+int thermal_init(t_moldyn *moldyn,t_random *random,int count);
+int scale_velocity(t_moldyn *moldyn,int count);
 double get_e_kin(t_atom *atom,int count);
 double get_e_pot(t_moldyn *moldyn);
 double get_total_energy(t_moldyn *moldyn);
diff --git a/posic.c b/posic.c
index 5ff3bc1..7c1c418 100644
--- a/posic.c
+++ b/posic.c
@@ -31,7 +31,6 @@ int main(int argc,char **argv) {
 	int count;
 
 	t_lj_params lj;
-	t_ho_params ho;
 
 	char fb[32]="saves/lj_test";
 
@@ -51,15 +50,18 @@ int main(int argc,char **argv) {
 	b=LEN_Y;
 	c=LEN_Z;
 
-	vis.dim.x=a*LC_SI;
-	vis.dim.y=b*LC_SI;
-	vis.dim.z=c*LC_SI;
+	/* set for 'bounding atoms' */
+	//vis.dim.x=a*LC_SI;
+	//vis.dim.y=b*LC_SI;
+	//vis.dim.z=c*LC_SI;
 
 	t=TEMPERATURE;
 
 	printf("placing silicon atoms ... ");
-	//count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
-	//printf("(%d) ok!\n",count);
+	count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
+	printf("(%d) ok!\n",count);
+
+	/* testing purpose
 	count=2;
 	si=malloc(2*sizeof(t_atom));
 	si[0].r.x=0.16e-9;
@@ -72,11 +74,12 @@ int main(int argc,char **argv) {
 	si[1].r.z=0;
 	si[1].element=SI;
 	si[1].mass=M_SI;
+	*/
 
 	printf("setting thermal fluctuations\n");
-	//thermal_init(si,&random,count,t);
-	v3_zero(&(si[0].v));
-	v3_zero(&(si[1].v));
+	thermal_init(si,&random,count,t);
+	//v3_zero(&(si[0].v));
+	//v3_zero(&(si[1].v));
 
 	/* check kinetic energy */
 
@@ -92,14 +95,10 @@ int main(int argc,char **argv) {
 	md.count=count;
 	md.atom=si;
 	md.potential=potential_lennard_jones;
-	//md.potential=potential_harmonic_oscillator;
 	md.force=force_lennard_jones;
-	//md.force=force_harmonic_oscillator;
-	//md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
-	md.cutoff=(0.4e-9);
-	md.cutoff_square=(0.6e-9*0.6e-9);
+	md.cutoff=R_CUTOFF;
+	md.cutoff_square=(R_CUTOFF*R_CUTOFF);
 	md.pot_params=&lj;
-	//md.pot_params=&ho;
 	md.integrate=velocity_verlet;
 	md.time_steps=RUNS;
 	md.tau=TAU;
@@ -113,9 +112,6 @@ int main(int argc,char **argv) {
 	lj.sigma12=lj.sigma6*lj.sigma6;
 	lj.epsilon=LJ_EPSILON_SI;
 
-	ho.equilibrium_distance=0.3e-9;
-	ho.spring_constant=1.0e-9;
-
 	u=get_e_pot(&md);
 
 	printf("potential energy: %.40f [J]\n",u);
@@ -147,133 +143,6 @@ int main(int argc,char **argv) {
 
 	rand_close(&random);
 	
-
-	//printf("starting velocity verlet: ");
-	//fflush(stdout);
-
-	//for(runs=0;runs<RUNS;runs++) {
-
-	/* 
-	 * velocity verlet
-	 *
-	 * r(t+h) = r(t) + h * dr/dt|t + h^2/2m * F(t)
-	 * dr/dt|(t+h) = dr/dt|t + h/2m * (F(t) + F(t+h))
-	 *
-	 */
-	//for(i=0;i<amount_si;i++) {
-//		/* calculation of new positions r(t+h) */
-//		si[i].x+=si[i].vx*tau;
-//		si[i].y+=si[i].vy*tau;
-//		si[i].z+=si[i].vz*tau;
-//		si[i].x+=(tau2*si[i].fx/m2);
-//		if(si[i].x>LX) si[i].x-=LEN_X;
-//		else if(si[i].x<-LX) si[i].x+=LEN_X;
-//		si[i].y+=(tau2*si[i].fy/m2);
-//		if(si[i].y>LY) si[i].y-=LEN_Y;
-//		else if(si[i].y<-LY) si[i].y+=LEN_Y;
-//		si[i].z+=(tau2*si[i].fz/m2);
-//		if(si[i].z>LZ) si[i].z-=LEN_Z;
-//		else if(si[i].z<-LZ) si[i].z+=LEN_Z;
-//		/* calculation of velocities v(t+h/2) */
-//		si[i].vx+=(tau*si[i].fx/m2);
-//		si[i].vy+=(tau*si[i].fy/m2);
-//		si[i].vz+=(tau*si[i].fz/m2);
-//		/* reset of forces */
-//		si[i].fx=.0;
-//		si[i].fy=.0;
-//		si[i].fz=.0;
-//	}
-//	for(i=0;i<amount_si;i++) {
-//		/* calculation of forces at new positions r(t+h) */
-//		for(j=0;j<i;j++) {
-//			deltax=si[i].x-si[j].x;
-//			if(deltax>LX) deltax-=LEN_X;
-//			else if(-deltax>LX) deltax+=LEN_X;
-//			deltax2=deltax*deltax;
-//			deltay=si[i].y-si[j].y;
-//			if(deltay>LY) deltay-=LEN_Y;
-//			else if(-deltay>LY) deltay+=LEN_Y;
-//			deltay2=deltay*deltay;
-//			deltaz=si[i].z-si[j].z;
-//			if(deltaz>LZ) deltaz-=LEN_Z;
-//			else if(-deltaz>LZ) deltaz+=LEN_Z;
-//			deltaz2=deltaz*deltaz;
-//			distance=deltax2+deltay2+deltaz2;
-//			if(distance<=R2_CUTOFF) {
-//				tmp=1.0/distance; // 1/r^2
-//				lj1=tmp; // 1/r^2
-//				tmp*=tmp; // 1/r^4
-//				lj1*=tmp; // 1/r^6
-//				tmp*=tmp; // 1/r^8
-//				lj2=tmp; // 1/r^8
-//				lj1*=tmp; // 1/r^14
-//				lj1*=LJ_SIGMA_12;
-//				lj2*=LJ_SIGMA_06;
-//				lj=-2*lj1+lj2;
-//				si[i].fx-=lj*deltax;
-//				si[i].fy-=lj*deltay;
-//				si[i].fz-=lj*deltaz;
-//				si[j].fx+=lj*deltax;
-//				si[j].fy+=lj*deltay;
-//				si[j].fz+=lj*deltaz;
-//			}
-//		}
-//	}
-//	for(i=0;i<amount_si;i++) {
-//		/* calculation of new velocities v(t+h) */
-//		si[i].vx+=(tau*si[i].fx/m2);
-//		si[i].vy+=(tau*si[i].fy/m2);
-//		si[i].vz+=(tau*si[i].fz/m2);
-//	}
-//
-//	if(!(runs%150)) {
-//
-//	/* rasmol script & xyz file */
-//	sprintf(xyz,"./saves/si-%.15f.xyz",time);
-//	sprintf(ppm,"./video/si-%.15f.ppm",time);
-//	fd1=open(xyz,O_WRONLY|O_CREAT|O_TRUNC);
-//	if(fd1<0) {
-//		perror("rasmol xyz file open");
-//		return -1;
-//	}
-//	dprintf(fd2,"load xyz %s\n",xyz);
-//	dprintf(fd2,"spacefill 200\n");
-//	dprintf(fd2,"rotate x 11\n");
-//	dprintf(fd2,"rotate y 13\n");
-//	dprintf(fd2,"set ambient 20\n");
-//	dprintf(fd2,"set specular on\n");
-//	dprintf(fd2,"zoom 400\n");
-//	dprintf(fd2,"write ppm %s\n",ppm);
-//	dprintf(fd2,"zap\n");
-//	dprintf(fd1,"%d\nsilicon\n",amount_si+9);
-//	for(i=0;i<amount_si;i++)
-//		dprintf(fd1,"Si %f %f %f %f\n",
-//		        si[i].x,si[i].y,si[i].z,time);
-//	dprintf(fd1,"H 0.0 0.0 0.0 %f\n",time);
-//	dprintf(fd1,"He %f %f %f %f\n",LX,LY,LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",LX,LY,-LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,-LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,-LZ,time);
-//	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,-LZ,time);
-//	close(fd1);
-//
-//	}
-//
-//	/* increase time */
-//	time+=tau;
-//	printf(".");
-//	fflush(stdout);
-//
-//	}
-//
-//	printf(" done\n");
-//	close(fd2);
-//	free(si);
-//
-
 	return 0;
 }
 
-- 
2.39.2