From 11bf064825efc474fa93c7270870856769a63e84 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Tue, 28 Nov 2006 10:35:42 +0000
Subject: [PATCH] more updates, now get the code running .... :/

---
 Makefile |   6 +-
 moldyn.c | 298 +++++++++++++++++++++++--------------------------------
 moldyn.h |  54 +++++++---
 3 files changed, 167 insertions(+), 191 deletions(-)

diff --git a/Makefile b/Makefile
index 3c8fcd1..b037404 100644
--- a/Makefile
+++ b/Makefile
@@ -5,10 +5,10 @@ LDFLAGS=-lm
 OBJS=init/init.o visual/visual.o math/math.o random/random.o list/list.o
 OBJS+=moldyn.o
 
-all: posic
+all: sic
 
-posic: $(OBJS)
+sic: $(OBJS)
 
 .PHONY:clean
 clean:
-	rm -f *.o posic */*.o
+	rm -f *.o sic */*.o
diff --git a/moldyn.c b/moldyn.c
index fe54330..c0b99eb 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -23,206 +23,155 @@
 #include "visual/visual.h"
 #include "list/list.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("-V <steps> <filebase> (rasmol file)\n");
-	printf("--- physics options ---\n");
-	printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
-	printf("-t <timestep tau> [s] (%.15f)\n",MOLDYN_TAU);
-	printf("-C <cutoff radius> [m] (%.15f)\n",MOLDYN_CUTOFF);
-	printf("-R <runs> (%d)\n",MOLDYN_RUNS);
-	printf(" -- integration algo --\n");
-	printf("  -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
-	printf("     0: velocity verlet\n");
-	printf(" -- potential --\n");
-	printf("  -P <number> <param1 param2 ...>\n");
-	printf("     0: harmonic oscillator\n");
-	printf("        param1: spring constant\n");
-	printf("        param2: equilibrium distance\n");
-	printf("     1: lennard jones\n");
-	printf("        param1: epsilon\n");
-	printf("        param2: sigma\n");
-	printf("\n");
+
+int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
+
+	//int ret;
+
+	//ret=moldyn_parse_argv(moldyn,argc,argv);
+	//if(ret<0) return ret;
+
+	memset(moldyn,0,sizeof(t_moldyn));
+
+	rand_init(&(moldyn->random),NULL,1);
+	moldyn->random.status|=RAND_STAT_VERBOSE;
 
 	return 0;
 }
 
-int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
+int moldyn_shutdown(t_moldyn *moldyn) {
 
-	int i;
+	link_cell_shutdown(moldyn);
+	moldyn_log_shutdown(moldyn);
+	rand_close(&(moldyn->random));
+	free(moldyn->atom);
 
-	memset(moldyn,0,sizeof(t_moldyn));
+	return 0;
+}
 
-	/* default values */
-	moldyn->t=MOLDYN_TEMP;
-	moldyn->tau=MOLDYN_TAU;
-	moldyn->time_steps=MOLDYN_RUNS;
-	moldyn->integrate=velocity_verlet;
-
-	/* 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 'S':
-					moldyn->swrite=atoi(argv[++i]);
-					strncpy(moldyn->sfb,argv[++i],64);
-					break;
-				case 'V':
-					moldyn->vwrite=atoi(argv[++i]);
-					strncpy(moldyn->vfb,argv[++i],64);
-					break;
-				case 'T':
-					moldyn->t=atof(argv[++i]);
-					break;
-				case 't':
-					moldyn->tau=atof(argv[++i]);
-					break;
-				case 'C':
-					moldyn->cutoff=atof(argv[++i]);
-					break;
-				case 'R':
-					moldyn->time_steps=atoi(argv[++i]);
-					break;
-				case 'I':
-	/* integration algorithm */
-	switch(atoi(argv[++i])) {
-		case MOLDYN_INTEGRATE_VERLET:
-			moldyn->integrate=velocity_verlet;
-			break;
-		default:
-			printf("unknown integration algo %s\n",argv[i]);
-			moldyn_usage(argv);
-			return -1;
-	}
+int set_int_alg(t_moldyn *moldyn,u8 algo) {
 
-				case 'P':
-	/* potential + params */
-	switch(atoi(argv[++i])) {
-		case MOLDYN_POTENTIAL_HO:
-			hop.spring_constant=atof(argv[++i]);
-			hop.equilibrium_distance=atof(argv[++i]);
-			moldyn->pot_params=malloc(sizeof(t_ho_params));
-			memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
-			moldyn->potential_force_function=harmonic_oscillator;
-			break;
-		case MOLDYN_POTENTIAL_LJ:
-			e=atof(argv[++i]);
-			s=atof(argv[++i]);
-			ljp.epsilon4=4*e;
-			ljp.sigma6=s*s*s*s*s*s;
-			ljp.sigma12=ljp.sigma6*ljp.sigma6;
-			moldyn->pot_params=malloc(sizeof(t_lj_params));
-			memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
-			moldyn->potential_force_function=lennard_jones;
+	switch(alg) {
+		case 'MOLDYN_INTEGRATE_VERLET':
+			moldyn->integrate=velocity_verlet;
 			break;
 		default:
-			printf("unknown potential %s\n",argv[i]);
-			moldyn_usage(argv);
-			return -1;
-	}
-
-				default:
-					printf("unknown option %s\n",argv[i]);
-					moldyn_usage(argv);
-					return -1;
-			}
-		} else {
-			moldyn_usage(argv);
+			printf("unknown integration algorithm: %02x\",alg);
 			return -1;
-		}
 	}
 
 	return 0;
 }
 
-int moldyn_log_init(t_moldyn *moldyn) {
+int set_cutoff(t_moldyn *moldyn,double cutoff) {
 
-	moldyn->lvstat=0;
-	t_visual *vis;
+	moldyn->cutoff=cutoff;
 
-	vis=&(moldyn->vis);
+	return 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;
-	}
+int set_temperature(t_moldyn *moldyn,double t) {
+	
+	moldyn->t=t;
 
-	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;
-	}
+	return 0;
+}
 
-	if(moldyn->swrite)
-		moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
+int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 
-	if((moldyn->vwrite)&&(vis)) {
-		moldyn->visual=vis;
-		visual_init(vis,moldyn->vfb);
-		moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
-	}
+	moldyn->dim.x=x;
+	moldyn->dim.y=y;
+	moldyn->dim.z=z;
 
-	moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
+	if(visualize) {
+		moldyn->vis.dim.x=x;
+		moldyn->vis.dim.y=y;
+		moldyn->vis.dim.z=z;
+	}
 
 	return 0;
 }
 
-int moldyn_log_shutdown(t_moldyn *moldyn) {
+int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 
-	if(moldyn->efd) close(moldyn->efd);
-	if(moldyn->mfd) close(moldyn->efd);
-	if(moldyn->dfd) close(moldyn->efd);
-	if(moldyn->visual) visual_tini(moldyn->visual);
+	if(x)
+		moldyn->status|=MOLDYN_STAT_PBX;
+
+	if(y)
+		moldyn->status|=MOLDYN_STAT_PBY;
+
+	if(z)
+		moldyn->status|=MOLDYN_STAT_PBZ;
 
 	return 0;
 }
 
-int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
+int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params) {
 
-	int ret;
+	switch(type) {
+		case MOLDYN_1BP:
+			moldyn->pf_func1b=func;
+			moldyn->pot1b_params=params;
+			break;
+		case MOLDYN_2BP:
+			moldyn->pf_func2b=func;
+			moldyn->pot2b_params=params;
+			break;
+		case MOLDYN_3BP:
+			moldyn->pf_func3b=func;
+			moldyn->pot3b_params=params;
+			break;
+		default:
+			printf("unknown potential type: %02x\n",type);
+			return -1;
+	}
 
-	ret=moldyn_parse_argv(moldyn,argc,argv);
-	if(ret<0) return ret;
+	return 0;
+}
 
-	ret=moldyn_log_init(moldyn);
-	if(ret<0) return ret;
+int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) {
 
-	rand_init(&(moldyn->random),NULL,1);
-	moldyn->random.status|=RAND_STAT_VERBOSE;
-
-	moldyn->status=0;
+	switch(type) {
+		case LOG_TOTAL_ENERGY:
+			moldyn->ewrite=timer;
+			moldyn->efd=open(fb,O_WRONLY|O_CREAT|O_TRUNC);
+			if(moldyn->efd<0) {
+				perror("[moldyn] efd open");
+				return moldyn->efd;
+			}
+			dprintf("# moldyn total energy log file\n");
+			break;
+		case LOG_TOTAL_MOMENTUM:
+			moldyn->mwrite=timer;
+			moldyn->mfd=open(fb,O_WRONLY|O_CREAT|O_TRUNC);
+			if(moldyn->mfd<0) {
+				perror("[moldyn] mfd open");
+				return moldyn->mfd;
+			}
+			dprintf("# moldyn total momentum log file\n");
+			break;
+		case SAVE_STEP:
+			moldyn->swrite=timer;
+			strncpy(moldyn->sfb,fb,63);
+			break;
+		case VISUAL_STEP:
+			moldyn->mwrite=timer;
+			strncpy(moldyn->vfb,fb,63);
+			visual_init(&(moldyn->vis),fb);
+			break;
+		default:
+			printf("unknown log mechanism: %02x\n",type);
+			return -1;
+	}
 
 	return 0;
 }
 
-int moldyn_shutdown(t_moldyn *moldyn) {
+int moldyn_log_shutdown(t_moldyn *moldyn) {
 
-	moldyn_log_shutdown(moldyn);
-	rand_close(&(moldyn->random));
-	free(moldyn->atom);
+	if(moldyn->efd) close(moldyn->efd);
+	if(moldyn->mfd) close(moldyn->mfd);
+	if(moldyn->visual) visual_tini(moldyn->visual);
 
 	return 0;
 }
@@ -239,6 +188,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	atom=moldyn->atom;
 
 	if(type==FCC) count*=4;
+
 	if(type==DIAMOND) count*=8;
 
 	atom=malloc(count*sizeof(t_atom));
@@ -269,6 +219,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 		return -1;
 	}
 
+	moldyn->count=count;
+
 	while(count) {
 		atom[count-1].element=element;
 		atom[count-1].mass=mass;
@@ -382,18 +334,18 @@ int scale_velocity(t_moldyn *moldyn) {
 	return 0;
 }
 
-double get_e_kin(t_atom *atom,int count) {
+double get_e_kin(t_moldyn *moldyn) {
 
 	int i;
-	double e;
+	t_atom *atom;
 
-	e=0.0;
+	atom=moldyn->atom;
+	moldyn->ekin=0.0;
 
-	for(i=0;i<count;i++) {
-		e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
-	}
+	for(i=0;i<moldyn->count;i++)
+		moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 
-	return e;
+	return moldyn->ekin;
 }
 
 double get_e_pot(t_moldyn *moldyn) {
@@ -403,18 +355,16 @@ double get_e_pot(t_moldyn *moldyn) {
 
 double get_total_energy(t_moldyn *moldyn) {
 
-	double e;
-
-	e=get_e_kin(moldyn->atom,moldyn->count);
-	e+=get_e_pot(moldyn);
-
-	return e;
+	return(get_e_kin(moldyn)+get_e_pot(moldyn));
 }
 
-t_3dvec get_total_p(t_atom *atom, int count) {
+t_3dvec get_total_p(t_moldyn *moldyn) {
 
 	t_3dvec p,p_total;
 	int i;
+	t_atom *atom;
+
+	atom=moldyn->atom;
 
 	v3_zero(&p_total);
 	for(i=0;i<count;i++) {
@@ -577,7 +527,7 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 	return 0;
 }
 
-int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau ) {
+int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
 
 	int count;
 	void *ptr;
diff --git a/moldyn.h b/moldyn.h
index a7fa60a..dafd86c 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -58,7 +58,7 @@ typedef struct s_moldyn_schedule {
 	int content_count;
 	int *runs;
 	double *tau;
-	int (*hook)(t_moldyn *,void *);
+	int (*hook)(void *moldyn,void *hook);
 	void *hook_params;
 } t_moldyn_schedule;
 
@@ -96,20 +96,20 @@ typedef struct s_moldyn {
 	double tau_square;	/* delta t squared */
 	double elapsed;		/* total elapsed time */
 
-	double energy;		/* energy */
+	double energy;		/* potential energy */
+	double ekin;		/* kinetic energy */
 
 	t_visual vis;		/* visualization/log/save interface structure */
 	u8 lvstat;		/* log & vis properties */
 	unsigned int ewrite;	/* how often to log energy */
-	char efb[64];		/* energy log filename */
 	int efd;		/* fd for energy log */
 	unsigned int mwrite;	/* how often to log momentum */
-	char mfb[64];		/* momentum log filename */
 	int mfd;		/* fd for momentum log */
 	unsigned int vwrite;	/* how often to visualize atom information */
 	char vfb[64];		/* visualization file name base */
 	void *visual;		/* pointer (hack!) */
 	unsigned int swrite;	/* how often to create a save file */
+	char sfb[64];		/* visualization file name base */
 
 	u8 status;		/* general moldyn properties */
 
@@ -126,6 +126,10 @@ typedef struct s_moldyn {
 #define MOLDYN_STAT_PBY			0x10	/* y */
 #define MOLDYN_STAT_PBZ			0x20	/* and z direction */
 
+#define MOLDYN_1BP			0x00
+#define MOLDYN_2BP			0x01
+#define MOLDYN_3BP			0x02
+
 
 /*
  *
@@ -236,6 +240,9 @@ typedef struct s_tersoff_mult_params {
 #define MOLDYN_SET_
 #define MOLDYN_SET_
 
+#define TRUE				1
+#define FALSE				0
+
 /*
  *
  * phsical values / constants
@@ -264,21 +271,32 @@ typedef struct s_tersoff_mult_params {
  *
  */
 
-int moldyn_usage(char **argv);
-int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv);
-int moldyn_log_init(t_moldyn *moldyn);
 int moldyn_init(t_moldyn *moldyn,int argc,char **argv);
 int moldyn_shutdown(t_moldyn *moldyn);
 
-int create_lattice(u8 type,int element,double mass,double lc,
-                   int a,int b,int c,t_atom **atom);
-int destroy_lattice(t_atom *atom);
+int set_int_alg(t_moldyn *moldyn,u8 algo);
+int set_cutoff(t_moldyn *moldyn,double cutoff);
+int set_temperature(t_moldyn *moldyn,double t);
+int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
+int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z);
+int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params);
+
+int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer);
+nt moldyn_log_shutdown(t_moldyn *moldyn);
+
+int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
+                   u8 attr,u8 bnum,int a,int b,int c);
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+             t_3dvec r,t_3dvec v);
+int destroy_atoms(t_moldyn *moldyn);
+
 int thermal_init(t_moldyn *moldyn);
 int scale_velocity(t_moldyn *moldyn);
-double get_e_kin(t_atom *atom,int count);
+
+double get_e_kin(t_moldyn *moldyn);
 double get_e_pot(t_moldyn *moldyn);
 double get_total_energy(t_moldyn *moldyn);
-t_3dvec get_total_p(t_atom *atom,int count);
+t_3dvec get_total_p(t_moldyn *moldyn);
 
 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t);
 
@@ -287,10 +305,18 @@ int link_cell_update(t_moldyn *moldyn);
 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell);
 int link_cell_shutdown(t_moldyn *moldyn);
 
+int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau);
+int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params);
+
 int moldyn_integrate(t_moldyn *moldyn);
 int velocity_verlet(t_moldyn *moldyn);
 
-int harmonic_oscillator(t_moldyn *moldyn);
-int lennard_jones(t_moldyn *moldyn);
+int potential_force_calc(t_moldyn *moldyn);
+int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
+int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_Atom *aj,u8 bc);
+int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
+int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
 
 #endif
-- 
2.20.1