X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=258985eb412a83ec6d3dd466c1c76c3c3a624c2c;hb=95cfeec6fbbfa975d5ac5b99ec3f7386ca3d6071;hp=415581c475fef1e8bb85b4fb8a1595135de3a88f;hpb=10645866d8c186ee58a40ff17a60ddb444fb4b82;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 415581c..258985e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -20,6 +20,58 @@ #include "moldyn.h" #include "report/report.h" +/* + * global variables, pse and atom colors (only needed here) + */ + +static char *pse_name[]={ + "*", + "H", + "He", + "Li", + "Be", + "B", + "C", + "N", + "O", + "F", + "Ne", + "Na", + "Mg", + "Al", + "Si", + "P", + "S", + "Cl", + "Ar", +}; + +static char *pse_col[]={ + "*", + "White", + "He", + "Li", + "Be", + "B", + "Gray", + "N", + "Blue", + "F", + "Ne", + "Na", + "Mg", + "Al", + "Yellow", + "P", + "S", + "Cl", + "Ar", +}; + +/* + * the moldyn functions + */ + int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { printf("[moldyn] init\n"); @@ -74,6 +126,18 @@ int set_cutoff(t_moldyn *moldyn,double cutoff) { return 0; } +int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) { + + moldyn->bondlen[0]=b0*b0; + moldyn->bondlen[1]=b1*b1; + if(bm<0) + moldyn->bondlen[2]=b0*b1; + else + moldyn->bondlen[2]=bm*bm; + + return 0; +} + int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; @@ -312,7 +376,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { break; case VISUAL_STEP: moldyn->vwrite=timer; - ret=visual_init(&(moldyn->vis),moldyn->vlsdir); + ret=visual_init(moldyn,moldyn->vlsdir); if(ret<0) { printf("[moldyn] visual init failure\n"); return ret; @@ -427,7 +491,6 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { moldyn->vlsdir); system(sc); } - if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; } @@ -558,6 +621,35 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, return 0; } +int del_atom(t_moldyn *moldyn,int tag) { + + t_atom *new,*old; + int cnt; + + old=moldyn->atom; + + new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom)); + if(!new) { + perror("[moldyn]malloc (del atom)"); + return -1; + } + + for(cnt=0;cntcount;cnt++) { + new[cnt-1]=old[cnt]; + new[cnt-1].tag=cnt-1; + } + + moldyn->count-=1; + moldyn->atom=new; + + free(old); + + return 0; +} + /* cubic init */ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { @@ -1483,8 +1575,7 @@ return 0; } if(v) { if(!(moldyn->total_steps%v)) { - visual_atoms(&(moldyn->vis),moldyn->time, - moldyn->atom,moldyn->count); + visual_atoms(moldyn); } } @@ -1975,6 +2066,157 @@ int analyze_bonds(t_moldyn *moldyn) { + return 0; +} + +/* + * visualization code + */ + +int visual_init(t_moldyn *moldyn,char *filebase) { + + strncpy(moldyn->vis.fb,filebase,128); + + return 0; +} + +int visual_atoms(t_moldyn *moldyn) { + + int i,j,fd; + char file[128+64]; + t_3dvec dim; + double help; + t_visual *v; + t_atom *atom; + t_atom *btom; + t_linkcell *lc; + t_list neighbour[27]; + u8 bc; + t_3dvec dist; + double d2; + u8 brand; + + v=&(moldyn->vis); + dim.x=v->dim.x; + dim.y=v->dim.y; + dim.z=v->dim.z; + atom=moldyn->atom; + lc=&(moldyn->lc); + + help=(dim.x+dim.y); + + sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time); + fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(fd<0) { + perror("open visual save file fd"); + return -1; + } + + /* write the actual data file */ + + // povray header + dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help); + + // atomic configuration + for(i=0;icount;i++) { + // atom type, positions, color and kinetic energy + dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element], + atom[i].r.x, + atom[i].r.y, + atom[i].r.z, + pse_col[atom[i].element], + atom[i].ekin); + + /* + * bond detection should usually be done by potential + * functions. brrrrr! EVIL! + * + * todo: potentials need to export a 'find_bonds' function! + */ + + // bonds between atoms + if(!(atom[i].attr&ATOM_ATTR_VB)) + continue; + link_cell_neighbour_index(moldyn, + (atom[i].r.x+moldyn->dim.x/2)/lc->x, + (atom[i].r.y+moldyn->dim.y/2)/lc->y, + (atom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour); + for(j=0;j<27;j++) { + list_reset_f(&neighbour[j]); + if(neighbour[j].start==NULL) + continue; + bc=jdnlc?0:1; + do { + btom=neighbour[j].current->data; + if(btom==&atom[i]) // skip identical atoms + continue; + //if(btom<&atom[i]) // skip half of them + // continue; + v3_sub(&dist,&(atom[i].r),&(btom->r)); + if(bc) check_per_bound(moldyn,&dist); + d2=v3_absolute_square(&dist); + brand=atom[i].brand; + if(brand==btom->brand) { + if(d2>moldyn->bondlen[brand]) + continue; + } + else { + if(d2>moldyn->bondlen[2]) + continue; + } + dprintf(fd,"# [B] %f %f %f %f %f %f\n", + atom[i].r.x,atom[i].r.y,atom[i].r.z, + btom->r.x,btom->r.y,btom->r.z); + } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT); + } + } + + // boundaries + if(dim.x) { + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,-dim.y/2,-dim.z/2, + dim.x/2,-dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,-dim.y/2,-dim.z/2, + -dim.x/2,dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dim.x/2,dim.y/2,-dim.z/2, + dim.x/2,-dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,dim.y/2,-dim.z/2, + dim.x/2,dim.y/2,-dim.z/2); + + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,-dim.y/2,dim.z/2, + dim.x/2,-dim.y/2,dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,-dim.y/2,dim.z/2, + -dim.x/2,dim.y/2,dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dim.x/2,dim.y/2,dim.z/2, + dim.x/2,-dim.y/2,dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,dim.y/2,dim.z/2, + dim.x/2,dim.y/2,dim.z/2); + + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,-dim.y/2,dim.z/2, + -dim.x/2,-dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + -dim.x/2,dim.y/2,dim.z/2, + -dim.x/2,dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dim.x/2,-dim.y/2,dim.z/2, + dim.x/2,-dim.y/2,-dim.z/2); + dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dim.x/2,dim.y/2,dim.z/2, + dim.x/2,dim.y/2,-dim.z/2); + } + + close(fd); + return 0; }