#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);
+ printf("unknown integration algorithm: %02x\",alg);
return -1;
}
- 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) {
+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 ret;
+int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params) {
- ret=moldyn_parse_argv(moldyn,argc,argv);
- if(ret<0) return 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_log_init(moldyn);
- if(ret<0) return ret;
+ return 0;
+}
- rand_init(&(moldyn->random),NULL,1);
- moldyn->random.status|=RAND_STAT_VERBOSE;
+int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) {
- 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;
}
-int create_lattice(u8 type,int element,double mass,double lc,
- int a,int b,int c,t_atom **atom) {
+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 count;
int ret;
t_3dvec origin;
+ t_atom *atom;
count=a*b*c;
+ atom=moldyn->atom;
if(type==FCC) count*=4;
+
if(type==DIAMOND) count*=8;
- *atom=malloc(count*sizeof(t_atom));
- if(*atom==NULL) {
+ atom=malloc(count*sizeof(t_atom));
+ if(atom==NULL) {
perror("malloc (atoms)");
return -1;
}
switch(type) {
case FCC:
- ret=fcc_init(a,b,c,lc,*atom,&origin);
+ ret=fcc_init(a,b,c,lc,atom,&origin);
break;
case DIAMOND:
- ret=diamond_init(a,b,c,lc,*atom,&origin);
+ ret=diamond_init(a,b,c,lc,atom,&origin);
break;
default:
printf("unknown lattice type (%02x)\n",type);
return -1;
}
+ moldyn->count=count;
+
while(count) {
- (*atom)[count-1].element=element;
- (*atom)[count-1].mass=mass;
+ atom[count-1].element=element;
+ atom[count-1].mass=mass;
+ atom[count-1].attr=attr;
+ atom[count-1].bnum=bnum;
count-=1;
}
return ret;
}
-int destroy_lattice(t_atom *atom) {
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+ t_3dvec r,t_3dvec v) {
- if(atom) free(atom);
+ t_atom *atom;
+ void *ptr;
+ int count;
+
+ atom=moldyn->atom;
+ count=++(moldyn->count);
+
+ ptr=realloc(atom,count*sizeof(t_atom));
+ if(!ptr) {
+ perror("[moldyn] realloc (add atom)");
+ return -1;
+ }
+
+ atom=ptr;
+ atom->r=r;
+ atom->v=v;
+ atom->element=element;
+ atom->bnum=bnum;
+ atom->attr=attr;
+
+ return 0;
+}
+
+int destroy_atoms(t_moldyn *moldyn) {
+
+ if(moldyn->atom) free(moldyn->atom);
return 0;
}
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) {
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++) {
return 0;
}
+int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
+
+ int count;
+ void *ptr;
+ t_moldyn_schedule *schedule;
+
+ schedule=moldyn->schedule;
+ count=++(schedule->content_count);
+
+ ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+ if(!ptr) {
+ perror("[moldyn] realloc (runs)");
+ return -1;
+ }
+ moldyn->schedule.runs[count-1]=runs;
+
+ ptr=realloc(schedule->tau,count*sizeof(double));
+ if(!ptr) {
+ perror("[moldyn] realloc (tau)");
+ return -1;
+ }
+ moldyn->schedule.tau[count-1]=tau;
+
+ return 0;
+}
+
+int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
+
+ moldyn->schedule.hook=hook;
+ moldyn->schedule.hook_params=hook_params;
+
+ return 0;
+}
+
/*
*
* 'integration of newtons equation' - algorithms
moldyn->potential_force_function(moldyn);
for(sched=0;sched<moldyn->schedule.content_count;sched++) {
- moldyn->tau=;
- moldyn->tau_square=;
- // hier weiter ...
+ /* setting amont of runs and finite time step size */
+ moldyn->tau=schedule->tau[sched];
+ moldyn->tau_square=moldyn->tau*moldyn->tau;
+ moldyn->timesteps=schedule->runs[sched];
+
+ /* integration according to schedule */
for(i=0;i<moldyn->time_steps;i++) {
}
}
+ /* check for hooks */
+ if(schedule->hook)
+ schedule->hook(moldyn,schedule->hook_params);
+
return 0;
}
projects / physik / posic.git / blobdiff