#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;
}
atom=moldyn->atom;
if(type==FCC) count*=4;
+
if(type==DIAMOND) count*=8;
atom=malloc(count*sizeof(t_atom));
return -1;
}
+ moldyn->count=count;
+
while(count) {
atom[count-1].element=element;
atom[count-1].mass=mass;
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 moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
int count;
void *ptr;
int content_count;
int *runs;
double *tau;
- int (*hook)(t_moldyn *,void *);
+ int (*hook)(void *moldyn,void *hook);
void *hook_params;
} t_moldyn_schedule;
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 */
#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
+
/*
*
#define MOLDYN_SET_
#define MOLDYN_SET_
+#define TRUE 1
+#define FALSE 0
+
/*
*
* phsical values / constants
*
*/
-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);
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
projects / physik / posic.git / commitdiff