X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=b52d51d4b790728b9e14f1fd38f48ebd41c8a7b2;hb=61d24d027511c3e96b2f94558dc1b31c72725767;hp=56817795283d87b38f648283cfb2d4eeac47306f;hpb=e84546faff88301a68958b58f66f546b03294503;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 5681779..b52d51d 100644 --- a/moldyn.c +++ b/moldyn.c @@ -20,6 +20,17 @@ #include "moldyn.h" #include "report/report.h" +/* potential includes */ +#include "potentials/harmonic_oscillator.h" +#include "potentials/lennard_jones.h" +#include "potentials/albe.h" +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else +#include "potentials/tersoff.h" +#endif + + /* * global variables, pse and atom colors (only needed here) */ @@ -120,24 +131,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) { int set_cutoff(t_moldyn *moldyn,double cutoff) { moldyn->cutoff=cutoff; + moldyn->cutoff_square=cutoff*cutoff; printf("[moldyn] cutoff [A]: %f\n",moldyn->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; @@ -191,15 +191,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->vis.dim.z=z; } - moldyn->dv=0.000001*moldyn->volume; - printf("[moldyn] dimensions in A and A^3 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); printf(" y: %f\n",moldyn->dim.y); printf(" z: %f\n",moldyn->dim.z); printf(" volume: %f\n",moldyn->volume); printf(" visualize simulation box: %s\n",visualize?"yes":"no"); - printf(" delta volume (pressure calc): %f\n",moldyn->dv); return 0; } @@ -231,58 +228,37 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func) { - - moldyn->func1b=func; - - return 0; -} - -int set_potential2b(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func2b=func; - - return 0; -} - -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { +int set_potential(t_moldyn *moldyn,u8 type) { - moldyn->func3b_j1=func; - - return 0; -} - -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j2=func; - - return 0; -} - -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j3=func; - - return 0; -} - -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k1=func; - - return 0; -} - -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k2=func; - - return 0; -} - -int set_potential_params(t_moldyn *moldyn,void *params) { - - moldyn->pot_params=params; + switch(type) { + case MOLDYN_POTENTIAL_TM: + moldyn->func1b=tersoff_mult_1bp; + moldyn->func3b_j1=tersoff_mult_3bp_j1; + moldyn->func3b_k1=tersoff_mult_3bp_k1; + moldyn->func3b_j2=tersoff_mult_3bp_j2; + moldyn->func3b_k2=tersoff_mult_3bp_k2; + // missing: check 2b bond func + break; + case MOLDYN_POTENTIAL_AM: + moldyn->func3b_j1=albe_mult_3bp_j1; + moldyn->func3b_k1=albe_mult_3bp_k1; + moldyn->func3b_j2=albe_mult_3bp_j2; + moldyn->func3b_k2=albe_mult_3bp_k2; + moldyn->check_2b_bond=albe_mult_check_2b_bond; + break; + case MOLDYN_POTENTIAL_HO: + moldyn->func2b=harmonic_oscillator; + moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; + break; + case MOLDYN_POTENTIAL_LJ: + moldyn->func2b=lennard_jones; + moldyn->check_2b_bond=lennard_jones_check_2b_bond; + break; + default: + printf("[moldyn] set potential: unknown type %02x\n", + type); + return -1; + } return 0; } @@ -370,12 +346,25 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { dprintf(moldyn->tfd,"# temperature log file\n"); printf("temperature (%d)\n",timer); break; + case LOG_VOLUME: + moldyn->vwrite=timer; + snprintf(filename,127,"%s/volume",moldyn->vlsdir); + moldyn->vfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->vfd<0) { + perror("[moldyn] volume log file\n"); + return moldyn->vfd; + } + dprintf(moldyn->vfd,"# volume log file\n"); + printf("volume (%d)\n",timer); + break; case SAVE_STEP: moldyn->swrite=timer; printf("save file (%d)\n",timer); break; case VISUAL_STEP: - moldyn->vwrite=timer; + moldyn->awrite=timer; ret=visual_init(moldyn,moldyn->vlsdir); if(ret<0) { printf("[moldyn] visual init failure\n"); @@ -595,7 +584,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, int count; atom=moldyn->atom; - count=(moldyn->count)++; + count=(moldyn->count)++; // asshole style! ptr=realloc(atom,(count+1)*sizeof(t_atom)); if(!ptr) { @@ -605,6 +594,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, moldyn->atom=ptr; atom=moldyn->atom; + + /* initialize new atom */ + memset(&(atom[count]),0,sizeof(t_atom)); atom[count].r=*r; atom[count].v=*v; atom[count].element=element; @@ -913,16 +905,26 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) { double virial_sum(t_moldyn *moldyn) { int i; - double v; t_virial *virial; /* virial (sum over atom virials) */ - v=0.0; + moldyn->virial=0.0; + moldyn->vir.xx=0.0; + moldyn->vir.yy=0.0; + moldyn->vir.zz=0.0; + moldyn->vir.xy=0.0; + moldyn->vir.xz=0.0; + moldyn->vir.yz=0.0; for(i=0;icount;i++) { virial=&(moldyn->atom[i].virial); - v+=(virial->xx+virial->yy+virial->zz); + moldyn->virial+=(virial->xx+virial->yy+virial->zz); + moldyn->vir.xx+=virial->xx; + moldyn->vir.yy+=virial->yy; + moldyn->vir.zz+=virial->zz; + moldyn->vir.xy+=virial->xy; + moldyn->vir.xz+=virial->xz; + moldyn->vir.yz+=virial->yz; } - moldyn->virial=v; /* global virial (absolute coordinates) */ virial=&(moldyn->gvir); @@ -954,12 +956,44 @@ double pressure_calc(t_moldyn *moldyn) { return moldyn->p; } +int average_reset(t_moldyn *moldyn) { + + printf("[moldyn] average reset\n"); + + /* update skip value */ + moldyn->avg_skip=moldyn->total_steps; + + /* kinetic energy */ + moldyn->k_sum=0.0; + moldyn->k2_sum=0.0; + + /* potential energy */ + moldyn->v_sum=0.0; + moldyn->v2_sum=0.0; + + /* temperature */ + moldyn->t_sum=0.0; + + /* virial */ + moldyn->virial_sum=0.0; + moldyn->gv_sum=0.0; + + /* pressure */ + moldyn->p_sum=0.0; + moldyn->gp_sum=0.0; + moldyn->tp_sum=0.0; + + return 0; +} + int average_and_fluctuation_calc(t_moldyn *moldyn) { + int denom; + if(moldyn->total_stepsavg_skip) return 0; - int denom=moldyn->total_steps+1-moldyn->avg_skip; + denom=moldyn->total_steps+1-moldyn->avg_skip; /* assume up to date energies, temperature, pressure etc */ @@ -992,6 +1026,8 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) { moldyn->p_avg=moldyn->p_sum/denom; moldyn->gp_sum+=moldyn->gp; moldyn->gp_avg=moldyn->gp_sum/denom; + moldyn->tp_sum+=moldyn->tp; + moldyn->tp_avg=moldyn->tp_sum/denom; return 0; } @@ -1031,9 +1067,11 @@ printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->coun double thermodynamic_pressure_calc(t_moldyn *moldyn) { - t_3dvec dim,*tp; - double u_up,u_down,dv; - double scale,p; + t_3dvec dim; + //t_3dvec *tp; + double h,dv; + double y0,y1; + double su,sd; t_atom *store; /* @@ -1043,54 +1081,56 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * */ - scale=0.00001; - dv=8*scale*scale*scale*moldyn->volume; - + /* store atomic configuration + dimension */ store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); return -1; } - - /* save unscaled potential energy + atom/dim configuration */ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); dim=moldyn->dim; + /* x1, y1 */ + sd=0.00001; + h=(1.0-sd)*(1.0-sd)*(1.0-sd); + su=pow(2.0-h,ONE_THIRD)-1.0; + dv=(1.0-h)*moldyn->volume; + /* scale up dimension and atom positions */ - scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_up=moldyn->energy; + y1=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* scale down dimension and atom positions */ - scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_down=moldyn->energy; + y0=moldyn->energy; /* calculate pressure */ - p=-(u_up-u_down)/dv; -printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR); + moldyn->tp=-(y1-y0)/(2.0*dv); - /* restore atomic configuration + dim */ + /* restore atomic configuration */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - - /* restore energy */ - potential_force_calc(moldyn); - link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); + //potential_force_calc(moldyn); - return p; + /* free store buffer */ + if(store) + free(store); + + return moldyn->tp; } double get_pressure(t_moldyn *moldyn) { @@ -1151,13 +1191,12 @@ int scale_volume(t_moldyn *moldyn) { /* scaling factor */ if(moldyn->pt_scale&P_SCALE_BERENDSEN) { - scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc; + scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc; scale=pow(scale,ONE_THIRD); } else { scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); } -moldyn->debug=scale; /* scale the atoms and dimensions */ scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); @@ -1339,7 +1378,7 @@ int link_cell_update(t_moldyn *moldyn) { i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x); j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y); k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z); - + #ifdef STATIC_LISTS p=0; while(lc->subcell[i+j*nx+k*nx*ny][p]!=0) @@ -1389,6 +1428,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, count2=27; a=nx*ny; + if(i>=nx||j>=ny||k>=nz) + printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n", + i,nx,j,ny,k,nz); + cell[0]=lc->subcell[i+j*nx+k*a]; for(ci=-1;ci<=1;ci++) { bx=0; @@ -1436,7 +1479,6 @@ int link_cell_shutdown(t_moldyn *moldyn) { for(i=0;icells;i++) { #ifdef STATIC_LISTS - //printf(" ---> %d free %p\n",i,lc->subcell[i]); free(lc->subcell[i]); #else //printf(" ---> %d free %p\n",i,lc->subcell[i].start); @@ -1500,7 +1542,7 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { int i; - unsigned int e,m,s,v,p,t; + unsigned int e,m,s,v,p,t,a; t_3dvec momentum; t_moldyn_schedule *sched; t_atom *atom; @@ -1522,12 +1564,12 @@ int moldyn_integrate(t_moldyn *moldyn) { m=moldyn->mwrite; s=moldyn->swrite; v=moldyn->vwrite; + a=moldyn->awrite; p=moldyn->pwrite; t=moldyn->twrite; /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; /* get current time */ gettimeofday(&t1,NULL); @@ -1540,14 +1582,14 @@ int moldyn_integrate(t_moldyn *moldyn) { /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) - printf("[moldyn] warning: cutoff > 0.5 x dim.x\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n"); if(moldyn->cutoff>0.5*moldyn->dim.y) - printf("[moldyn] warning: cutoff > 0.5 x dim.y\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n"); if(moldyn->cutoff>0.5*moldyn->dim.z) - printf("[moldyn] warning: cutoff > 0.5 x dim.z\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n"); ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass; if(ds>0.05*moldyn->nnd) - printf("[moldyn] warning: forces too high / tau too small!\n"); + printf("[moldyn] WARNING: forces too high / tau too small!\n"); /* zero absolute time */ moldyn->time=0.0; @@ -1583,6 +1625,9 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); + //thermodynamic_pressure_calc(moldyn); + + /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); /* p/t scaling */ @@ -1612,9 +1657,10 @@ int moldyn_integrate(t_moldyn *moldyn) { if(p) { if(!(moldyn->total_steps%p)) { dprintf(moldyn->pfd, - "%f %f %f %f %f\n",moldyn->time, + "%f %f %f %f %f %f %f\n",moldyn->time, moldyn->p/BAR,moldyn->p_avg/BAR, - moldyn->gp/BAR,moldyn->gp_avg/BAR); + moldyn->gp/BAR,moldyn->gp_avg/BAR, + moldyn->tp/BAR,moldyn->tp_avg/BAR); } } if(t) { @@ -1624,6 +1670,12 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->time,moldyn->t,moldyn->t_avg); } } + if(v) { + if(!(moldyn->total_steps%v)) { + dprintf(moldyn->vfd, + "%f %f\n",moldyn->time,moldyn->volume); + } + } if(s) { if(!(moldyn->total_steps%s)) { snprintf(dir,128,"%s/s-%07.f.save", @@ -1639,8 +1691,8 @@ int moldyn_integrate(t_moldyn *moldyn) { close(fd); } } - if(v) { - if(!(moldyn->total_steps%v)) { + if(a) { + if(!(moldyn->total_steps%a)) { visual_atoms(moldyn); } } @@ -1650,13 +1702,14 @@ int moldyn_integrate(t_moldyn *moldyn) { /* get current time */ gettimeofday(&t2,NULL); - printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", - sched->count,i, - moldyn->t,moldyn->t_avg, - moldyn->p_avg/BAR,moldyn->gp_avg/BAR, - moldyn->volume, - (int)(t2.tv_sec-t1.tv_sec)); - fflush(stdout); +printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", + sched->count,i,moldyn->total_steps, + moldyn->t,moldyn->t_avg, + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); + + fflush(stdout); /* copy over time */ t1=t2; @@ -1670,8 +1723,8 @@ int moldyn_integrate(t_moldyn *moldyn) { /* check for hooks */ if(sched->hook) { - printf("\n ## schedule hook %d/%d start ##\n", - sched->count+1,sched->total_sched-1); + printf("\n ## schedule hook %d start ##\n", + sched->count); sched->hook(moldyn,sched->hook_params); printf(" ## schedule hook end ##\n"); } @@ -1699,6 +1752,9 @@ int velocity_verlet(t_moldyn *moldyn) { tau_square=moldyn->tau_square; for(i=0;if.x,itom->f.y,itom->f.z); if(moldyn->time>DSTART&&moldyn->timeatom[5832].f.x); - printf(" y: %0.40f\n",moldyn->atom[5832].f.y); - printf(" z: %0.40f\n",moldyn->atom[5832].f.z); + printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x); + printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y); + printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z); } #endif - /* calculate global virial */ + /* some postprocessing */ for(i=0;igvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; - moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; - moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; - moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; - moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; - moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; + /* calculate global virial */ + moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x; + moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y; + moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z; + moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x; + moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x; + moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y; + + /* check forces regarding the given timestep */ + if(v3_norm(&(itom[i].f))>\ + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + printf("[moldyn] WARNING: pfc (high force: atom %d)\n", + i); } return 0; @@ -2174,6 +2243,8 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { int fd; int cnt,size; + int fsize; + int corr; fd=open(file,O_RDONLY); if(fd<0) { @@ -2181,25 +2252,48 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { return fd; } + fsize=lseek(fd,0,SEEK_END); + lseek(fd,0,SEEK_SET); + size=sizeof(t_moldyn); - cnt=read(fd,moldyn,size); - if(cnt!=size) { - perror("[moldyn] load save file read (moldyn)"); - return cnt; + + while(size) { + cnt=read(fd,moldyn,size); + if(cnt<0) { + perror("[moldyn] load save file read (moldyn)"); + return cnt; + } + size-=cnt; } size=moldyn->count*sizeof(t_atom); + /* correcting possible atom data offset */ + corr=0; + if(fsize!=sizeof(t_moldyn)+size) { + corr=fsize-sizeof(t_moldyn)-size; + printf("[moldyn] WARNING: lsf (illegal file size)\n"); + printf(" moifying offset:\n"); + printf(" - current pos: %d\n",sizeof(t_moldyn)); + printf(" - atom size: %d\n",size); + printf(" - file size: %d\n",fsize); + printf(" => correction: %d\n",corr); + lseek(fd,corr,SEEK_CUR); + } + moldyn->atom=(t_atom *)malloc(size); if(moldyn->atom==NULL) { perror("[moldyn] load save file malloc (atoms)"); return -1; } - cnt=read(fd,moldyn->atom,size); - if(cnt!=size) { - perror("[moldyn] load save file read (atoms)"); - return cnt; + while(size) { + cnt=read(fd,moldyn->atom,size); + if(cnt<0) { + perror("[moldyn] load save file read (atoms)"); + return cnt; + } + size-=cnt; } // hooks etc ... @@ -2207,46 +2301,28 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { return 0; } -int moldyn_load(t_moldyn *moldyn) { +int moldyn_free_save_file(t_moldyn *moldyn) { - // later ... + free(moldyn->atom); return 0; } -/* - * post processing functions - */ - -int get_line(int fd,char *line,int max) { - - int count,ret; - - count=0; - - while(1) { - if(count==max) return count; - ret=read(fd,line+count,1); - if(ret<=0) return ret; - if(line[count]=='\n') { - line[count]='\0'; - return count+1; - } - count+=1; - } -} +int moldyn_load(t_moldyn *moldyn) { -int pair_correlation_init(t_moldyn *moldyn,double dr) { + // later ... - return 0; } -int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { +/* + * function to find/callback all combinations of 2 body bonds + */ + +int process_2b_bonds(t_moldyn *moldyn,void *data, + int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc)) { - int slots; - double *stat; - int i,j; t_linkcell *lc; #ifdef STATIC_LISTS int *neighbour[27]; @@ -2254,38 +2330,13 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { #else t_list neighbour[27]; #endif + u8 bc; t_atom *itom,*jtom; + int i,j; t_list *this; - unsigned char bc; - t_3dvec dist; - double d,norm; - int o,s; - unsigned char ibrand; lc=&(moldyn->lc); - slots=(int)(moldyn->cutoff/dr); - o=2*slots; - - printf("[moldyn] pair correlation calc info:\n"); - printf(" time: %f\n",moldyn->time); - printf(" count: %d\n",moldyn->count); - printf(" cutoff: %f\n",moldyn->cutoff); - printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg); - - if(ptr!=NULL) { - stat=(double *)ptr; - } - else { - stat=(double *)malloc(3*slots*sizeof(double)); - if(stat==NULL) { - perror("[moldyn] pair correlation malloc"); - return -1; - } - } - - memset(stat,0,3*slots*sizeof(double)); - link_cell_init(moldyn,VERBOSE); itom=moldyn->atom; @@ -2298,9 +2349,6 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { (itom[i].r.z+moldyn->dim.z/2)/lc->x, neighbour); - /* brand of atom i */ - ibrand=itom[i].brand; - for(j=0;j<27;j++) { bc=(jdnlc)?0:1; @@ -2324,78 +2372,313 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { jtom=this->current->data; #endif - if(jtom==&(itom[i])) - continue; + /* process bond */ + process(moldyn,&(itom[i]),jtom,data,bc); - /* only count pairs once */ - if(itom[i].tag>jtom->tag) - continue; +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + } + } - /* - * pair correlation calc - */ + return 0; - /* distance */ - v3_sub(&dist,&(jtom->r),&(itom[i].r)); - if(bc) check_per_bound(moldyn,&dist); - d=v3_absolute_square(&dist); +} - /* ignore if greater cutoff */ - if(d>moldyn->cutoff_square) - continue; +/* + * post processing functions + */ - /* fill the slots */ - d=sqrt(d); - s=(int)(d/dr); +int get_line(int fd,char *line,int max) { - if(ibrand!=jtom->brand) { - /* mixed */ - stat[s]+=1; - } - else { - /* type a - type a bonds */ - if(ibrand==0) - stat[s+slots]+=1; - else - /* type b - type b bonds */ - stat[s+o]+=1; - } + int count,ret; -#ifdef STATIC_LISTS - } -#else - } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); -#endif + count=0; + + while(1) { + if(count==max) return count; + ret=read(fd,line+count,1); + if(ret<=0) return ret; + if(line[count]=='\n') { + memset(line+count,0,max-count-1); + //line[count]='\0'; + return count+1; + } + count+=1; + } +} + +int pair_correlation_init(t_moldyn *moldyn,double dr) { + + + return 0; +} + +int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) { + + int i; + t_atom *atom; + t_3dvec dist; + double d2; + int a_cnt; + int b_cnt; + + atom=moldyn->atom; + dc[0]=0; + dc[1]=0; + dc[2]=0; + a_cnt=0; + b_cnt=0; + + for(i=0;icount;i++) { + + v3_sub(&dist,&(atom[i].r),&(atom[i].r_0)); + check_per_bound(moldyn,&dist); + d2=v3_absolute_square(&dist); + + if(atom[i].brand) { + b_cnt+=1; + dc[1]+=d2; + } + else { + a_cnt+=1; + dc[0]+=d2; + } + + dc[2]+=d2; + } + + dc[0]*=(1.0/(6.0*moldyn->time*a_cnt)); + dc[1]*=(1.0/(6.0*moldyn->time*b_cnt)); + dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count)); + + return 0; +} + +int bonding_analyze(t_moldyn *moldyn,double *cnt) { + + return 0; +} + +int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, + t_atom *jtom,void *data,u8 bc) { + + t_3dvec dist; + double d; + int s; + t_pcc *pcc; + + /* only count pairs once, + * skip same atoms */ + if(itom->tag>=jtom->tag) + return 0; + + /* + * pair correlation calc + */ + + /* get pcc data */ + pcc=data; + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater cutoff */ + if(d>moldyn->cutoff_square) + return 0; + + /* fill the slots */ + d=sqrt(d); + s=(int)(d/pcc->dr); + + /* should never happen but it does 8) - + * related to -ffloat-store problem! */ + if(s>=pcc->o1) { + printf("[moldyn] WARNING: pcc (%d/%d)", + s,pcc->o1); + printf("\n"); + s=pcc->o1-1; + } + + if(itom->brand!=jtom->brand) { + /* mixed */ + pcc->stat[s]+=1; + } + else { + /* type a - type a bonds */ + if(itom->brand==0) + pcc->stat[s+pcc->o1]+=1; + else + /* type b - type b bonds */ + pcc->stat[s+pcc->o2]+=1; + } + + return 0; +} + +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { + + t_pcc pcc; + double norm; + int i; + + pcc.dr=dr; + pcc.o1=moldyn->cutoff/dr; + pcc.o2=2*pcc.o1; + + if(pcc.o1*dr<=moldyn->cutoff) + printf("[moldyn] WARNING: pcc (low #slots)\n"); + + printf("[moldyn] pair correlation calc info:\n"); + printf(" time: %f\n",moldyn->time); + printf(" count: %d\n",moldyn->count); + printf(" cutoff: %f\n",moldyn->cutoff); + printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg); + + if(ptr!=NULL) { + pcc.stat=(double *)ptr; + } + else { + pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double)); + if(pcc.stat==NULL) { + perror("[moldyn] pair correlation malloc"); + return -1; } } - /* normalization - for(i=1;i 2 pi r r dr - norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr; - stat[i]/=norm; - stat[slots+i]/=norm; - stat[o+i]/=norm; + // ... and actually it's a constant times r^2 + norm=i*i*dr*dr; + pcc.stat[i]/=norm; + pcc.stat[pcc.o1+i]/=norm; + pcc.stat[pcc.o2+i]/=norm; } - */ + /* */ if(ptr==NULL) { /* todo: store/print pair correlation function */ - free(stat); + free(pcc.stat); } - free(moldyn->atom); + return 0; +} - link_cell_shutdown(moldyn); +int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { + + t_ba *ba; + t_3dvec dist; + double d; + + if(itom->tag>=jtom->tag) + return 0; + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater or equal cutoff */ + if(d>moldyn->cutoff_square) + return 0; + + /* check for potential bond */ + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + d=sqrt(d); + + /* now count this bonding ... */ + ba=data; + + /* increase total bond counter + * ... double counting! + */ + ba->tcnt+=2; + + if(itom->brand==0) + ba->acnt[jtom->tag]+=1; + else + ba->bcnt[jtom->tag]+=1; + + if(jtom->brand==0) + ba->acnt[itom->tag]+=1; + else + ba->bcnt[itom->tag]+=1; return 0; } -int analyze_bonds(t_moldyn *moldyn) { +int bond_analyze(t_moldyn *moldyn,double *quality) { - - + // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total + + int qcnt; + int ccnt,cset; + t_ba ba; + int i; + t_atom *atom; + + ba.acnt=malloc(moldyn->count*sizeof(int)); + if(ba.acnt==NULL) { + perror("[moldyn] bond analyze malloc (a)"); + return -1; + } + memset(ba.acnt,0,moldyn->count*sizeof(int)); + + ba.bcnt=malloc(moldyn->count*sizeof(int)); + if(ba.bcnt==NULL) { + perror("[moldyn] bond analyze malloc (b)"); + return -1; + } + memset(ba.bcnt,0,moldyn->count*sizeof(int)); + + ba.tcnt=0; + qcnt=0; + ccnt=0; + cset=0; + + atom=moldyn->atom; + + process_2b_bonds(moldyn,&ba,bond_analyze_process); + + for(i=0;icount;i++) { + if(atom[i].brand==0) { + if((ba.acnt[i]==0)&(ba.bcnt[i]==4)) + qcnt+=4; + } + else { + if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { + qcnt+=4; + ccnt+=1; + } + cset+=1; + } + } + + printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset); + printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt); + + if(quality) { + quality[0]=1.0*ccnt/cset; + quality[1]=1.0*qcnt/ba.tcnt; + } + else { + printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt); + printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt); + } return 0; } @@ -2411,39 +2694,48 @@ int visual_init(t_moldyn *moldyn,char *filebase) { return 0; } +int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { + + t_vb *vb; + + vb=data; + + if(itom->tag>=jtom->tag) + return 0; + + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB)) + dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n", + itom->r.x,itom->r.y,itom->r.z, + jtom->r.x,jtom->r.y,jtom->r.z); + + return 0; +} + int visual_atoms(t_moldyn *moldyn) { - int i,j,fd; + int i; char file[128+64]; t_3dvec dim; double help; t_visual *v; t_atom *atom; - t_atom *btom; - t_linkcell *lc; -#ifdef STATIC_LISTS - int *neighbour[27]; - int p; -#else - t_list neighbour[27]; -#endif - u8 bc; - t_3dvec dist; - double d2; - u8 brand; + t_vb vb; 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) { + vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(vb.fd<0) { perror("open visual save file fd"); return -1; } @@ -2451,118 +2743,65 @@ int visual_atoms(t_moldyn *moldyn) { /* write the actual data file */ // povray header - dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + dprintf(vb.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++) { + 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++) { - bc=jdnlc?0:1; -#ifdef STATIC_LISTS - p=0; - while(neighbour[j][p]!=0) { - btom=&(atom[neighbour[j][p]]); - p++; -#else - list_reset_f(&neighbour[j]); - if(neighbour[j].start==NULL) - continue; - do { - btom=neighbour[j].current->data; -#endif - 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); -#ifdef STATIC_LISTS - } -#else - } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT); -#endif - } - } - + dprintf(vb.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); + + // bonds between atoms + process_2b_bonds(moldyn,&vb,visual_bonds_process); + // boundaries if(dim.x) { - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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", + dprintf(vb.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); + close(vb.fd); return 0; }