]> hackdaworld.org Git - physik/posic.git/commitdiff
more updates, now get the code running .... :/
authorhackbard <hackbard>
Tue, 28 Nov 2006 10:35:42 +0000 (10:35 +0000)
committerhackbard <hackbard>
Tue, 28 Nov 2006 10:35:42 +0000 (10:35 +0000)
Makefile
moldyn.c
moldyn.h

index 3c8fcd1157eeaf5184c809ac70b4fd906e2c1c95..b0374046069a109e4440430976f0237538df15e2 100644 (file)
--- 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
index fe54330f002f414d345295ff73dd4d1e6b6b557d..c0b99eb70f72e6e103fecaa5157ce5679573e4b4 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
 #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;
index a7fa60a2311a4487712ef3f646adbb950430c180..dafd86cb67d415636833bafc147d860dd595850f 100644 (file)
--- 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