X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=edda007d9290ff6cfe1001d4105c3f8a9ab12863;hb=e1080fc0dd66b0cf5b7715c5e99e7a34ac04a8cf;hp=881f0b698565fe7495ee151ecca63239d32e74cf;hpb=c0ddf2bdd8067456f39f6b63fe2261624ebde6b7;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 881f0b6..edda007 100644 --- a/moldyn.c +++ b/moldyn.c @@ -13,17 +13,74 @@ #include #include #include +#include +#include #include #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"); memset(moldyn,0,sizeof(t_moldyn)); + moldyn->argc=argc; + moldyn->args=argv; + rand_init(&(moldyn->random),NULL,1); moldyn->random.status|=RAND_STAT_VERBOSE; @@ -69,11 +126,23 @@ 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; - printf("[moldyn] temperature: %f\n",moldyn->t_ref); + printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref); return 0; } @@ -82,7 +151,7 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { moldyn->p_ref=p_ref; - printf("[moldyn] pressure: %f\n",moldyn->p_ref); + printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR); return 0; } @@ -162,34 +231,66 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) { +int set_potential1b(t_moldyn *moldyn,pf_func1b func) { moldyn->func1b=func; - moldyn->pot1b_params=params; return 0; } -int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { +int set_potential2b(t_moldyn *moldyn,pf_func2b func) { moldyn->func2b=func; - moldyn->pot2b_params=params; return 0; } -int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) { +int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { + + 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->func2b_post=func; - moldyn->pot2b_params=params; + moldyn->func3b_k2=func; return 0; } -int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { +int set_potential_params(t_moldyn *moldyn,void *params) { - moldyn->func3b=func; - moldyn->pot3b_params=params; + moldyn->pot_params=params; + + return 0; +} + +int set_avg_skip(t_moldyn *moldyn,int skip) { + + printf("[moldyn] skip %d steps before starting average calc\n",skip); + moldyn->avg_skip=skip; return 0; } @@ -243,13 +344,39 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { dprintf(moldyn->efd,"# total momentum log file\n"); printf("total momentum (%d)\n",timer); break; + case LOG_PRESSURE: + moldyn->pwrite=timer; + snprintf(filename,127,"%s/pressure",moldyn->vlsdir); + moldyn->pfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->pfd<0) { + perror("[moldyn] pressure log file\n"); + return moldyn->pfd; + } + dprintf(moldyn->pfd,"# pressure log file\n"); + printf("pressure (%d)\n",timer); + break; + case LOG_TEMPERATURE: + moldyn->twrite=timer; + snprintf(filename,127,"%s/temperature",moldyn->vlsdir); + moldyn->tfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->tfd<0) { + perror("[moldyn] temperature log file\n"); + return moldyn->tfd; + } + dprintf(moldyn->tfd,"# temperature log file\n"); + printf("temperature (%d)\n",timer); + break; case SAVE_STEP: moldyn->swrite=timer; printf("save file (%d)\n",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; @@ -265,18 +392,52 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { perror("[moldyn] report fd open"); return moldyn->rfd; } - snprintf(filename,127,"%s/plot.scr",moldyn->vlsdir); - moldyn->pfd=open(filename, - O_WRONLY|O_CREAT|O_EXCL, - S_IRUSR|S_IWUSR); - if(moldyn->pfd<0) { - perror("[moldyn] plot fd open"); - return moldyn->pfd; + printf("report -> "); + if(moldyn->efd) { + snprintf(filename,127,"%s/e_plot.scr", + moldyn->vlsdir); + moldyn->epfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->epfd<0) { + perror("[moldyn] energy plot fd open"); + return moldyn->epfd; + } + dprintf(moldyn->epfd,e_plot_script); + close(moldyn->epfd); + printf("energy "); + } + if(moldyn->pfd) { + snprintf(filename,127,"%s/pressure_plot.scr", + moldyn->vlsdir); + moldyn->ppfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->ppfd<0) { + perror("[moldyn] p plot fd open"); + return moldyn->ppfd; + } + dprintf(moldyn->ppfd,pressure_plot_script); + close(moldyn->ppfd); + printf("pressure "); + } + if(moldyn->tfd) { + snprintf(filename,127,"%s/temperature_plot.scr", + moldyn->vlsdir); + moldyn->tpfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->tpfd<0) { + perror("[moldyn] t plot fd open"); + return moldyn->tpfd; + } + dprintf(moldyn->tpfd,temperature_plot_script); + close(moldyn->tpfd); + printf("temperature "); } dprintf(moldyn->rfd,report_start, moldyn->rauthor,moldyn->rtitle); - dprintf(moldyn->pfd,plot_script); - close(moldyn->pfd); + printf("\n"); break; default: printf("unknown log type: %02x\n",type); @@ -291,21 +452,45 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { char sc[256]; printf("[moldyn] log shutdown\n"); - if(moldyn->efd) close(moldyn->efd); + if(moldyn->efd) { + close(moldyn->efd); + if(moldyn->rfd) { + dprintf(moldyn->rfd,report_energy); + snprintf(sc,255,"cd %s && gnuplot e_plot.scr", + moldyn->vlsdir); + system(sc); + } + } if(moldyn->mfd) close(moldyn->mfd); + if(moldyn->pfd) { + close(moldyn->pfd); + if(moldyn->rfd) + dprintf(moldyn->rfd,report_pressure); + snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr", + moldyn->vlsdir); + system(sc); + } + if(moldyn->tfd) { + close(moldyn->tfd); + if(moldyn->rfd) + dprintf(moldyn->rfd,report_temperature); + snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr", + moldyn->vlsdir); + system(sc); + } if(moldyn->rfd) { dprintf(moldyn->rfd,report_end); close(moldyn->rfd); - snprintf(sc,255,"cd %s && gnuplot plot.scr",moldyn->vlsdir); + snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); - snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); + snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); - snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); - system(sc); - snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir); + snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); } - if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; } @@ -315,11 +500,11 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { */ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, - u8 attr,u8 brand,int a,int b,int c) { + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) { int new,count; int ret; - t_3dvec origin; + t_3dvec orig; void *ptr; t_atom *atom; @@ -339,21 +524,35 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->atom=ptr; atom=&(moldyn->atom[count]); - - v3_zero(&origin); + + /* no atoms on the boundaries (only reason: it looks better!) */ + if(!origin) { + orig.x=0.5*lc; + orig.y=0.5*lc; + orig.z=0.5*lc; + } + else { + orig.x=origin->x; + orig.y=origin->y; + orig.z=origin->z; + } switch(type) { case CUBIC: - origin.x=0.5*lc; - origin.y=0.5*lc; - origin.z=0.5*lc; - ret=cubic_init(a,b,c,lc,atom,&origin); + set_nn_dist(moldyn,lc); + ret=cubic_init(a,b,c,lc,atom,&orig); break; case FCC: - ret=fcc_init(a,b,c,lc,atom,NULL); + if(!origin) + v3_scale(&orig,&orig,0.5); + set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); + ret=fcc_init(a,b,c,lc,atom,&orig); break; case DIAMOND: - ret=diamond_init(a,b,c,lc,atom,&origin); + if(!origin) + v3_scale(&orig,&orig,0.25); + set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); + ret=diamond_init(a,b,c,lc,atom,&orig); break; default: printf("unknown lattice type (%02x)\n",type); @@ -379,11 +578,78 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom[ret].brand=brand; atom[ret].tag=count+ret; check_per_bound(moldyn,&(atom[ret].r)); + atom[ret].r_0=atom[ret].r; } + /* update total system mass */ + total_mass_calc(moldyn); + return ret; } +int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + t_3dvec *r,t_3dvec *v) { + + t_atom *atom; + void *ptr; + int count; + + atom=moldyn->atom; + count=(moldyn->count)++; + + ptr=realloc(atom,(count+1)*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (add atom)"); + return -1; + } + moldyn->atom=ptr; + + atom=moldyn->atom; + atom[count].r=*r; + atom[count].v=*v; + atom[count].element=element; + atom[count].mass=mass; + atom[count].brand=brand; + atom[count].tag=count; + atom[count].attr=attr; + check_per_bound(moldyn,&(atom[count].r)); + atom[count].r_0=atom[count].r; + + /* update total system mass */ + total_mass_calc(moldyn); + + 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) { @@ -426,65 +692,55 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { int count; - int i,j; + int i,j,k,l; t_3dvec o,r,n; t_3dvec basis[3]; - double help[3]; - double x,y,z; - - x=a*lc; - y=b*lc; - z=c*lc; - if(origin) v3_copy(&o,origin); - else v3_zero(&o); + count=0; + if(origin) + v3_copy(&o,origin); + else + v3_zero(&o); /* construct the basis */ - for(i=0;i<3;i++) { - for(j=0;j<3;j++) { - if(i!=j) help[j]=0.5*lc; - else help[j]=.0; - } - v3_set(&basis[i],help); - } + memset(basis,0,3*sizeof(t_3dvec)); + basis[0].x=0.5*lc; + basis[0].y=0.5*lc; + basis[1].x=0.5*lc; + basis[1].z=0.5*lc; + basis[2].y=0.5*lc; + basis[2].z=0.5*lc; - v3_zero(&r); - count=0; - /* fill up the room */ r.x=o.x; - while(r.xatom; - count=(moldyn->count)++; - - ptr=realloc(atom,(count+1)*sizeof(t_atom)); - if(!ptr) { - perror("[moldyn] realloc (add atom)"); - return -1; - } - moldyn->atom=ptr; - - atom=moldyn->atom; - atom[count].r=*r; - atom[count].v=*v; - atom[count].element=element; - atom[count].mass=mass; - atom[count].brand=brand; - atom[count].tag=count; - atom[count].attr=attr; - - return 0; -} - int destroy_atoms(t_moldyn *moldyn) { if(moldyn->atom) free(moldyn->atom); @@ -592,6 +819,18 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) { return 0; } +double total_mass_calc(t_moldyn *moldyn) { + + int i; + + moldyn->mass=0.0; + + for(i=0;icount;i++) + moldyn->mass+=moldyn->atom[i].mass; + + return moldyn->mass; +} + double temperature_calc(t_moldyn *moldyn) { /* assume up to date kinetic energy, which is 3/2 N k_B T */ @@ -671,39 +910,145 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) { return p; } -double pressure_calc(t_moldyn *moldyn) { +double virial_sum(t_moldyn *moldyn) { int i; double v; t_virial *virial; - /* - * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i ) - * - * virial = f_i r_i - */ - + /* virial (sum over atom virials) */ v=0.0; for(i=0;icount;i++) { virial=&(moldyn->atom[i].virial); v+=(virial->xx+virial->yy+virial->zz); } + moldyn->virial=v; + + /* global virial (absolute coordinates) */ + virial=&(moldyn->gvir); + moldyn->gv=virial->xx+virial->yy+virial->zz; + + return moldyn->virial; +} + +double pressure_calc(t_moldyn *moldyn) { + + /* + * PV = NkT + + * with W = 1/3 sum_i f_i r_i (- skipped!) + * virial = sum_i f_i r_i + * + * => P = (2 Ekin + virial) / (3V) + */ - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+v; + /* assume up to date virial & up to date kinetic energy */ + + /* pressure (atom virials) */ + moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); + /* pressure (absolute coordinates) */ + moldyn->gp=2.0*moldyn->ekin+moldyn->gv; + moldyn->gp/=(3.0*moldyn->volume); + return moldyn->p; -} +} + +int average_and_fluctuation_calc(t_moldyn *moldyn) { + + int denom; + + if(moldyn->total_stepsavg_skip) + return 0; + + denom=moldyn->total_steps+1-moldyn->avg_skip; + + /* assume up to date energies, temperature, pressure etc */ + + /* kinetic energy */ + moldyn->k_sum+=moldyn->ekin; + moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin); + moldyn->k_avg=moldyn->k_sum/denom; + moldyn->k2_avg=moldyn->k2_sum/denom; + moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg); + + /* potential energy */ + moldyn->v_sum+=moldyn->energy; + moldyn->v2_sum+=(moldyn->energy*moldyn->energy); + moldyn->v_avg=moldyn->v_sum/denom; + moldyn->v2_avg=moldyn->v2_sum/denom; + moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg); + + /* temperature */ + moldyn->t_sum+=moldyn->t; + moldyn->t_avg=moldyn->t_sum/denom; + + /* virial */ + moldyn->virial_sum+=moldyn->virial; + moldyn->virial_avg=moldyn->virial_sum/denom; + moldyn->gv_sum+=moldyn->gv; + moldyn->gv_avg=moldyn->gv_sum/denom; + + /* pressure */ + moldyn->p_sum+=moldyn->p; + moldyn->p_avg=moldyn->p_sum/denom; + moldyn->gp_sum+=moldyn->gp; + moldyn->gp_avg=moldyn->gp_sum/denom; + + return 0; +} + +int get_heat_capacity(t_moldyn *moldyn) { + + double temp2,ighc; + + /* averages needed for heat capacity calc */ + if(moldyn->total_stepsavg_skip) + return 0; + + /* (temperature average)^2 */ + temp2=moldyn->t_avg*moldyn->t_avg; + printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", + moldyn->t_avg); + + /* ideal gas contribution */ + ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; + printf(" ideal gas contribution: %f\n", + ighc/moldyn->mass*KILOGRAM/JOULE); + + /* specific heat for nvt ensemble */ + moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc; + moldyn->c_v_nvt/=moldyn->mass; + + /* specific heat for nve ensemble */ + moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2))); + moldyn->c_v_nve/=moldyn->mass; + + printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE); + printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE); +printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); + + return 0; +} double thermodynamic_pressure_calc(t_moldyn *moldyn) { - t_3dvec dim,*tp; - double u,p; - double scale; + t_3dvec dim; + //t_3dvec *tp; + double u_up,u_down,dv; + double scale,p; t_atom *store; - tp=&(moldyn->tp); + /* + * dU = - p dV + * + * => p = - dU/dV + * + */ + + scale=0.00001; + dv=8*scale*scale*scale*moldyn->volume; + store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); @@ -711,77 +1056,44 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { } /* save unscaled potential energy + atom/dim configuration */ - u=moldyn->energy; memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); dim=moldyn->dim; - /* derivative with respect to x direction */ - scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z); - scale_dim(moldyn,scale,TRUE,0,0); - scale_atoms(moldyn,scale,TRUE,0,0); - link_cell_shutdown(moldyn); - link_cell_init(moldyn); - potential_force_calc(moldyn); - tp->x=(moldyn->energy-u)/moldyn->dv; - p=tp->x*tp->x; - - /* restore atomic configuration + dim */ - memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); - moldyn->dim=dim; - - /* derivative with respect to y direction */ - scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z); - scale_dim(moldyn,scale,0,TRUE,0); - scale_atoms(moldyn,scale,0,TRUE,0); + /* scale up dimension and atom positions */ + scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->y=(moldyn->energy-u)/moldyn->dv; - p+=tp->y*tp->y; + u_up=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - /* derivative with respect to z direction */ - scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y); - scale_dim(moldyn,scale,0,0,TRUE); - scale_atoms(moldyn,scale,0,0,TRUE); - link_cell_shutdown(moldyn); - link_cell_init(moldyn); - potential_force_calc(moldyn); - tp->z=(moldyn->energy-u)/moldyn->dv; - p+=tp->z*tp->z; - - /* restore atomic configuration + dim */ - memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); - moldyn->dim=dim; - - printf("dU/dV komp addiert = %f %f %f\n",tp->x,tp->y,tp->z); - - scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD); - -printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); - scale_dim(moldyn,scale,1,1,1); - scale_atoms(moldyn,scale,1,1,1); + /* scale down dimension and atom positions */ + scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); -printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); - - printf("dU/dV einfach = %f\n",((moldyn->energy-u)/moldyn->dv)/ATM); + u_down=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); /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* restore energy */ - moldyn->energy=u; + potential_force_calc(moldyn); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); - return sqrt(p); + return p; } double get_pressure(t_moldyn *moldyn) { @@ -790,12 +1102,18 @@ double get_pressure(t_moldyn *moldyn) { } -int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { +int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) { t_3dvec *dim; dim=&(moldyn->dim); + if(dir==SCALE_UP) + scale=1.0+scale; + + if(dir==SCALE_DOWN) + scale=1.0-scale; + if(x) dim->x*=scale; if(y) dim->y*=scale; if(z) dim->z*=scale; @@ -803,11 +1121,17 @@ int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { return 0; } -int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { +int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) { int i; t_3dvec *r; + if(dir==SCALE_UP) + scale=1.0+scale; + + if(dir==SCALE_DOWN) + scale=1.0-scale; + for(i=0;icount;i++) { r=&(moldyn->atom[i].r); if(x) r->x*=scale; @@ -820,68 +1144,55 @@ int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { int scale_volume(t_moldyn *moldyn) { - t_atom *atom; t_3dvec *dim,*vdim; - double scale,v; - t_virial virial; + double scale; t_linkcell *lc; - int i; - atom=moldyn->atom; - dim=&(moldyn->dim); vdim=&(moldyn->vis.dim); + dim=&(moldyn->dim); lc=&(moldyn->lc); - memset(&virial,0,sizeof(t_virial)); - - for(i=0;icount;i++) { - virial.xx+=atom[i].virial.xx; - virial.yy+=atom[i].virial.yy; - virial.zz+=atom[i].virial.zz; - virial.xy+=atom[i].virial.xy; - virial.xz+=atom[i].virial.xz; - virial.yz+=atom[i].virial.yz; - } - - /* just a guess so far ... */ - v=virial.xx+virial.yy+virial.zz; - -printf("%f\n",v); - /* get pressure from virial */ - moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v; - moldyn->p/=moldyn->volume; -printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM); - - /* scale factor */ - if(moldyn->pt_scale&P_SCALE_BERENDSEN) - scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc); - else - /* should actually never be used */ - scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0); - -printf("scale = %f\n",scale); - /* actual scaling */ - dim->x*=scale; - dim->y*=scale; - dim->z*=scale; - if(vdim->x) vdim->x=dim->x; - if(vdim->y) vdim->y=dim->y; - if(vdim->z) vdim->z=dim->z; - moldyn->volume*=(scale*scale*scale); - - /* check whether we need a new linkcell init */ - if((dim->x/moldyn->cutoff!=lc->nx)|| - (dim->y/moldyn->cutoff!=lc->ny)|| - (dim->z/moldyn->cutoff!=lc->nx)) { + /* scaling factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) { + 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); + scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + + /* visualize dimensions */ + if(vdim->x!=0) { + vdim->x=dim->x; + vdim->y=dim->y; + vdim->z=dim->z; + } + + /* recalculate scaled volume */ + moldyn->volume=dim->x*dim->y*dim->z; + + /* adjust/reinit linkcell */ + if(((int)(dim->x/moldyn->cutoff)!=lc->nx)|| + ((int)(dim->y/moldyn->cutoff)!=lc->ny)|| + ((int)(dim->z/moldyn->cutoff)!=lc->nx)) { link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); + } else { + lc->x*=scale; + lc->y*=scale; + lc->z*=scale; } return 0; } -double get_e_kin(t_moldyn *moldyn) { +double e_kin_calc(t_moldyn *moldyn) { int i; t_atom *atom; @@ -889,17 +1200,14 @@ double get_e_kin(t_moldyn *moldyn) { atom=moldyn->atom; moldyn->ekin=0.0; - for(i=0;icount;i++) - moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + for(i=0;icount;i++) { + atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + moldyn->ekin+=atom[i].ekin; + } return moldyn->ekin; } -double update_e_kin(t_moldyn *moldyn) { - - return(get_e_kin(moldyn)); -} - double get_total_energy(t_moldyn *moldyn) { return(moldyn->ekin+moldyn->energy); @@ -939,7 +1247,7 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { /* linked list / cell method */ -int link_cell_init(t_moldyn *moldyn) { +int link_cell_init(t_moldyn *moldyn,u8 vol) { t_linkcell *lc; int i; @@ -953,20 +1261,57 @@ int link_cell_init(t_moldyn *moldyn) { lc->y=moldyn->dim.y/lc->ny; lc->nz=moldyn->dim.z/moldyn->cutoff; lc->z=moldyn->dim.z/lc->nz; - lc->cells=lc->nx*lc->ny*lc->nz; + +#ifdef STATIC_LISTS + lc->subcell=malloc(lc->cells*sizeof(int*)); +#else lc->subcell=malloc(lc->cells*sizeof(t_list)); +#endif + + if(lc->subcell==NULL) { + perror("[moldyn] cell init (malloc)"); + return -1; + } if(lc->cells<27) printf("[moldyn] FATAL: less then 27 subcells!\n"); - printf("[moldyn] initializing linked cells (%d)\n",lc->cells); + if(vol) { +#ifdef STATIC_LISTS + printf("[moldyn] initializing 'static' linked cells (%d)\n", + lc->cells); +#else + printf("[moldyn] initializing 'dynamic' linked cells (%d)\n", + lc->cells); +#endif + printf(" x: %d x %f A\n",lc->nx,lc->x); + printf(" y: %d x %f A\n",lc->ny,lc->y); + printf(" z: %d x %f A\n",lc->nz,lc->z); + } +#ifdef STATIC_LISTS + /* list init */ + for(i=0;icells;i++) { + lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int)); + if(lc->subcell[i]==NULL) { + perror("[moldyn] list init (malloc)"); + return -1; + } + /* + if(i==0) + printf(" ---> %d malloc %p (%p)\n", + i,lc->subcell[0],lc->subcell); + */ + } +#else for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); +#endif + + /* update the list */ + link_cell_update(moldyn); - link_cell_update(moldyn); - return 0; } @@ -976,7 +1321,9 @@ int link_cell_update(t_moldyn *moldyn) { int nx,ny; t_atom *atom; t_linkcell *lc; - double x,y,z; +#ifdef STATIC_LISTS + int p; +#endif atom=moldyn->atom; lc=&(moldyn->lc); @@ -984,25 +1331,50 @@ int link_cell_update(t_moldyn *moldyn) { nx=lc->nx; ny=lc->ny; - x=moldyn->dim.x/2; - y=moldyn->dim.y/2; - z=moldyn->dim.z/2; - for(i=0;icells;i++) +#ifdef STATIC_LISTS + memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int)); +#else list_destroy_f(&(lc->subcell[i])); - +#endif + for(count=0;countcount;count++) { 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); - list_add_immediate_f(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]), + +#ifdef STATIC_LISTS + p=0; + while(lc->subcell[i+j*nx+k*nx*ny][p]!=0) + p++; + + if(p>=MAX_ATOMS_PER_LIST) { + printf("[moldyn] FATAL: amount of atoms too high!\n"); + return -1; + } + + lc->subcell[i+j*nx+k*nx*ny][p]=count; +#else + list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]), &(atom[count])); + /* + if(j==0&&k==0) + printf(" ---> %d %d malloc %p (%p)\n", + i,count,lc->subcell[i].current,lc->subcell); + */ +#endif } return 0; } -int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { +int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, +#ifdef STATIC_LISTS + int **cell +#else + t_list *cell +#endif + ) { t_linkcell *lc; int a; @@ -1065,8 +1437,14 @@ int link_cell_shutdown(t_moldyn *moldyn) { lc=&(moldyn->lc); - for(i=0;inx*lc->ny*lc->nz;i++) - list_destroy_f(&(moldyn->lc.subcell[i])); + for(i=0;icells;i++) { +#ifdef STATIC_LISTS + free(lc->subcell[i]); +#else + //printf(" ---> %d free %p\n",i,lc->subcell[i].start); + list_destroy_f(&(lc->subcell[i])); +#endif + } free(lc->subcell); @@ -1124,38 +1502,43 @@ 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; - t_3dvec p; + unsigned int e,m,s,v,p,t; + t_3dvec momentum; t_moldyn_schedule *sched; t_atom *atom; int fd; char dir[128]; double ds; double energy_scale; + struct timeval t1,t2; + //double tp; sched=&(moldyn->schedule); atom=moldyn->atom; /* initialize linked cell method */ - link_cell_init(moldyn); + link_cell_init(moldyn,VERBOSE); /* logging & visualization */ e=moldyn->ewrite; m=moldyn->mwrite; s=moldyn->swrite; v=moldyn->vwrite; + p=moldyn->pwrite; + t=moldyn->twrite; /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; - /* energy scaling factor */ - energy_scale=moldyn->count*EV; + /* get current time */ + gettimeofday(&t1,NULL); -printf("debug: %f\n",moldyn->atom[0].f.x); /* calculate initial forces */ potential_force_calc(moldyn); -printf("debug: %f\n",moldyn->atom[0].f.x); +#ifdef DEBUG +//return 0; +#endif /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) @@ -1170,6 +1553,7 @@ printf("debug: %f\n",moldyn->atom[0].f.x); /* zero absolute time */ moldyn->time=0.0; + moldyn->total_steps=0; /* debugging, ignore */ moldyn->debug=0; @@ -1178,13 +1562,17 @@ printf("debug: %f\n",moldyn->atom[0].f.x); printf("[moldyn] integration start, go get a coffee ...\n"); /* executing the schedule */ - for(sched->count=0;sched->counttotal_sched;sched->count++) { + sched->count=0; + while(sched->counttotal_sched) { /* setting amount of runs and finite time step size */ moldyn->tau=sched->tau[sched->count]; moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->time_steps=sched->runs[sched->count]; + /* energy scaling factor (might change!) */ + energy_scale=moldyn->count*EV; + /* integration according to schedule */ for(i=0;itime_steps;i++) { @@ -1192,20 +1580,22 @@ printf("debug: %f\n",moldyn->atom[0].f.x); /* integration step */ moldyn->integrate(moldyn); + /* calculate kinetic energy, temperature and pressure */ + e_kin_calc(moldyn); + temperature_calc(moldyn); + virial_sum(moldyn); + pressure_calc(moldyn); + average_and_fluctuation_calc(moldyn); + /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) scale_volume(moldyn); - update_e_kin(moldyn); - temperature_calc(moldyn); - pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); - /* check for log & visualization */ if(e) { - if(!(i%e)) + if(!(moldyn->total_steps%e)) dprintf(moldyn->efd, "%f %f %f %f\n", moldyn->time,moldyn->ekin/energy_scale, @@ -1213,17 +1603,35 @@ printf("debug: %f\n",moldyn->atom[0].f.x); get_total_energy(moldyn)/energy_scale); } if(m) { - if(!(i%m)) { - p=get_total_p(moldyn); + if(!(moldyn->total_steps%m)) { + momentum=get_total_p(moldyn); dprintf(moldyn->mfd, - "%f %f\n",moldyn->time,v3_norm(&p)); + "%f %f %f %f %f\n",moldyn->time, + momentum.x,momentum.y,momentum.z, + v3_norm(&momentum)); + } + } + if(p) { + if(!(moldyn->total_steps%p)) { + dprintf(moldyn->pfd, + "%f %f %f %f %f\n",moldyn->time, + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->gp/BAR,moldyn->gp_avg/BAR); + } + } + if(t) { + if(!(moldyn->total_steps%t)) { + dprintf(moldyn->tfd, + "%f %f %f\n", + moldyn->time,moldyn->t,moldyn->t_avg); } } if(s) { - if(!(i%s)) { + if(!(moldyn->total_steps%s)) { snprintf(dir,128,"%s/s-%07.f.save", moldyn->vlsdir,moldyn->time); - fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT); + fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT, + S_IRUSR|S_IWUSR); if(fd<0) perror("[moldyn] save fd open"); else { write(fd,moldyn,sizeof(t_moldyn)); @@ -1234,26 +1642,45 @@ printf("debug: %f\n",moldyn->atom[0].f.x); } } if(v) { - if(!(i%v)) { - visual_atoms(&(moldyn->vis),moldyn->time, - moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, debug: %f", - sched->count,i,moldyn->p/ATM); - fflush(stdout); + if(!(moldyn->total_steps%v)) { + visual_atoms(moldyn); } } + /* display progress */ + //if(!(moldyn->total_steps%10)) { + /* get current time */ + gettimeofday(&t2,NULL); + +printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", + sched->count,i,moldyn->total_steps, + moldyn->t,moldyn->t_avg, + moldyn->p_avg/BAR,moldyn->gp_avg/BAR, + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); + + fflush(stdout); + + /* copy over time */ + t1=t2; + //} + /* increase absolute time */ moldyn->time+=moldyn->tau; + moldyn->total_steps+=1; } /* check for hooks */ - if(sched->hook) + if(sched->hook) { + printf("\n ## schedule hook %d start ##\n", + sched->count); sched->hook(moldyn,sched->hook_params); + printf(" ## schedule hook end ##\n"); + } - /* get a new info line */ - printf("\n"); + /* increase the schedule counter */ + sched->count+=1; } @@ -1318,19 +1745,31 @@ int potential_force_calc(t_moldyn *moldyn) { t_atom *itom,*jtom,*ktom; t_virial *virial; t_linkcell *lc; +#ifdef STATIC_LISTS + int *neighbour_i[27]; + int p,q; + t_atom *atom; +#else t_list neighbour_i[27]; t_list neighbour_i2[27]; t_list *this,*that; +#endif u8 bc_ij,bc_ik; int dnlc; count=moldyn->count; itom=moldyn->atom; lc=&(moldyn->lc); +#ifdef STATIC_LISTS + atom=moldyn->atom; +#endif /* reset energy */ moldyn->energy=0.0; + /* reset global virial */ + memset(&(moldyn->gvir),0,sizeof(t_virial)); + /* reset force, site energy and virial of every atom */ for(i=0;ifunc1b(moldyn,&(itom[i])); + if(moldyn->func1b) + moldyn->func1b(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; @@ -1371,54 +1813,188 @@ int potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; + /* first loop over atoms j */ + if(moldyn->func2b) { + for(j=0;j<27;j++) { + + bc_ij=(jattr&ATOM_ATTR_2BP)& + (itom[i].attr&ATOM_ATTR_2BP)) { + moldyn->func2b(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + } +#else + this=&(neighbour_i[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + do { + jtom=this->current->data; + + if(jtom==&(itom[i])) + continue; + + if((jtom->attr&ATOM_ATTR_2BP)& + (itom[i].attr&ATOM_ATTR_2BP)) { + moldyn->func2b(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + + } + } + + /* 3 body potential/force */ + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + /* copy the neighbour lists */ +#ifdef STATIC_LISTS + /* no copy needed for static lists */ +#else + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + + /* second loop over atoms j */ for(j=0;j<27;j++) { + bc_ij=(jstart==NULL) continue; - bc_ij=(jcurrent->data; +#endif if(jtom==&(itom[i])) continue; - if((jtom->attr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; - /* 3 body potential/force */ + /* reset 3bp run */ + moldyn->run3bp=1; - if(!(itom[i].attr&ATOM_ATTR_3BP)|| - !(jtom->attr&ATOM_ATTR_3BP)) - continue; + if(moldyn->func3b_j1) + moldyn->func3b_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); - /* copy the neighbour lists */ - memcpy(neighbour_i2,neighbour_i, - 27*sizeof(t_list)); + /* in first j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) + continue; + + /* first loop over atoms k */ + if(moldyn->func3b_k1) { - /* get neighbours of i */ for(k=0;k<27;k++) { + bc_ik=(kstart==NULL) continue; + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + if(ktom==jtom) + continue; + + if(ktom==&(itom[i])) + continue; + + moldyn->func3b_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); +#endif + + } + + } + + if(moldyn->func3b_j2) + moldyn->func3b_j2(moldyn, + &(itom[i]), + jtom, + bc_ij); + + /* second loop over atoms k */ + if(moldyn->func3b_k2) { + + for(k=0;k<27;k++) { + bc_ik=(kstart==NULL) + continue; + do { ktom=that->current->data; +#endif if(!(ktom->attr&ATOM_ATTR_3BP)) continue; @@ -1429,37 +2005,66 @@ int potential_force_calc(t_moldyn *moldyn) { if(ktom==&(itom[i])) continue; - moldyn->func3b(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func3b_k2(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS + } +#else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); +#endif + } + } /* 2bp post function */ - if(moldyn->func2b_post) { - moldyn->func2b_post(moldyn, - &(itom[i]), - jtom,bc_ij); + if(moldyn->func3b_j3) { + moldyn->func3b_j3(moldyn, + &(itom[i]), + jtom,bc_ij); } - +#ifdef STATIC_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } + +#ifdef DEBUG + //printf("\n\n"); +#endif +#ifdef VDEBUG + printf("\n\n"); +#endif } #ifdef DEBUG -printf("\n\n"); -#endif -#ifdef VDEBUG -printf("\n\n"); + //printf("\nATOM 0: %f %f %f\n\n",itom->f.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); + } #endif + /* calculate global virial */ + 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; + } + return 0; } @@ -1467,7 +2072,8 @@ printf("\n\n"); * virial calculation */ -inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { +//inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { +int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { a->virial.xx+=f->x*d->x; a->virial.yy+=f->y*d->y; @@ -1480,10 +2086,11 @@ inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { } /* - * periodic boundayr checking + * periodic boundary checking */ -inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { +//inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { +int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { double x,y,z; t_3dvec *dim; @@ -1510,826 +2117,464 @@ inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { return 0; } - /* - * example potentials + * debugging / critical check functions */ -/* harmonic oscillator potential and force */ - -int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - - t_ho_params *params; - t_3dvec force,distance; - double d,f; - double sc,equi_dist; - - params=moldyn->pot2b_params; - sc=params->spring_constant; - equi_dist=params->equilibrium_distance; - - if(air),&(ai->r)); - - if(bc) check_per_bound(moldyn,&distance); - d=v3_norm(&distance); - if(d<=moldyn->cutoff) { - moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); - /* f = -grad E; grad r_ij = -1 1/r_ij distance */ - f=sc*(1.0-equi_dist/d); - v3_scale(&force,&distance,f); - v3_add(&(ai->f),&(ai->f),&force); - virial_calc(ai,&force,&distance); - virial_calc(aj,&force,&distance); /* f and d signe switched */ - v3_scale(&force,&distance,-f); - v3_add(&(aj->f),&(aj->f),&force); - } + t_atom *atom; + t_3dvec *dim; + int i; + double x; + u8 byte; + int j,k; - return 0; -} + atom=moldyn->atom; + dim=&(moldyn->dim); + x=dim->x/2; -/* lennard jones potential & force for one sort of atoms */ - -int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - - t_lj_params *params; - t_3dvec force,distance; - double d,h1,h2; - double eps,sig6,sig12; - - params=moldyn->pot2b_params; - eps=params->epsilon4; - sig6=params->sigma6; - sig12=params->sigma12; - - if(air),&(ai->r)); - if(bc) check_per_bound(moldyn,&distance); - d=v3_absolute_square(&distance); /* 1/r^2 */ - if(d<=moldyn->cutoff_square) { - d=1.0/d; /* 1/r^2 */ - h2=d*d; /* 1/r^4 */ - h2*=d; /* 1/r^6 */ - h1=h2*h2; /* 1/r^12 */ - moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); - h2*=d; /* 1/r^8 */ - h1*=d; /* 1/r^14 */ - h2*=6*sig6; - h1*=12*sig12; - d=+h1-h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(aj->f),&(aj->f),&force); - v3_scale(&force,&force,-1.0); /* f = - grad E */ - v3_add(&(ai->f),&(ai->f),&force); - virial_calc(ai,&force,&distance); - virial_calc(aj,&force,&distance); /* f and d signe switched */ + for(i=0;icount;i++) { + if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) { + printf("FATAL: atom %d: x: %.20f (%.20f)\n", + i,atom[i].r.x,dim->x/2); + printf("diagnostic:\n"); + printf("-----------\natom.r.x:\n"); + for(j=0;j<8;j++) { + memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1); + for(k=0;k<8;k++) + printf("%d%c", + ((byte)&(1<=dim->y/2||-atom[i].r.y>dim->y/2) + printf("FATAL: atom %d: y: %.20f (%.20f)\n", + i,atom[i].r.y,dim->y/2); + if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2) + printf("FATAL: atom %d: z: %.20f (%.20f)\n", + i,atom[i].r.z,dim->z/2); } return 0; } /* - * tersoff potential & force for 2 sorts of atoms + * restore function */ -/* create mixed terms from parameters and set them */ -int tersoff_mult_complete_params(t_tersoff_mult_params *p) { - - printf("[moldyn] tersoff parameter completion\n"); - p->S2[0]=p->S[0]*p->S[0]; - p->S2[1]=p->S[1]*p->S[1]; - p->Smixed=sqrt(p->S[0]*p->S[1]); - p->S2mixed=p->Smixed*p->Smixed; - p->Rmixed=sqrt(p->R[0]*p->R[1]); - p->Amixed=sqrt(p->A[0]*p->A[1]); - p->Bmixed=sqrt(p->B[0]*p->B[1]); - p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); - p->mu_m=0.5*(p->mu[0]+p->mu[1]); - - printf("[moldyn] tersoff mult parameter info:\n"); - printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed); - printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed); - printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV); - printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); - printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], - p->lambda_m); - printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); - printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]); - printf(" n | %f | %f\n",p->n[0],p->n[1]); - printf(" c | %f | %f\n",p->c[0],p->c[1]); - printf(" d | %f | %f\n",p->d[0],p->d[1]); - printf(" h | %f | %f\n",p->h[0],p->h[1]); - printf(" chi | %f \n",p->chi); - - return 0; -} +int moldyn_read_save_file(t_moldyn *moldyn,char *file) { -/* tersoff 1 body part */ -int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { - - int brand; - t_tersoff_mult_params *params; - t_tersoff_exchange *exchange; - - brand=ai->brand; - params=moldyn->pot1b_params; - exchange=&(params->exchange); + int fd; + int cnt,size; - /* - * simple: point constant parameters only depending on atom i to - * their right values - */ + fd=open(file,O_RDONLY); + if(fd<0) { + perror("[moldyn] load save file open"); + return fd; + } - exchange->beta_i=&(params->beta[brand]); - exchange->n_i=&(params->n[brand]); - exchange->c_i=&(params->c[brand]); - exchange->d_i=&(params->d[brand]); - exchange->h_i=&(params->h[brand]); + size=sizeof(t_moldyn); - exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i)); - exchange->ci2=params->c[brand]*params->c[brand]; - exchange->di2=params->d[brand]*params->d[brand]; - exchange->ci2di2=exchange->ci2/exchange->di2; + while(size) { + cnt=read(fd,moldyn,size); + if(cnt<0) { + perror("[moldyn] load save file read (moldyn)"); + return cnt; + } + size-=cnt; + } - return 0; -} - -/* tersoff 2 body part */ -int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - - t_tersoff_mult_params *params; - t_tersoff_exchange *exchange; - t_3dvec dist_ij,force; - double d_ij,d_ij2; - double A,B,R,S,S2,lambda,mu; - double f_r,df_r; - double f_c,df_c; - int brand; - double s_r; - double arg; - - params=moldyn->pot2b_params; - brand=aj->brand; - exchange=&(params->exchange); - - /* clear 3bp and 2bp post run */ - exchange->run3bp=0; - exchange->run2bp_post=0; - - /* reset S > r > R mark */ - exchange->d_ij_between_rs=0; - - /* - * calc of 2bp contribution of V_ij and dV_ij/ji - * - * for Vij and dV_ij we need: - * - f_c_ij, df_c_ij - * - f_r_ij, df_r_ij - * - * for dV_ji we need: - * - f_c_ji = f_c_ij, df_c_ji = df_c_ij - * - f_r_ji = f_r_ij; df_r_ji = df_r_ij - * - */ + size=moldyn->count*sizeof(t_atom); - /* constants */ - if(brand==ai->brand) { - S=params->S[brand]; - S2=params->S2[brand]; - R=params->R[brand]; - A=params->A[brand]; - B=params->B[brand]; - lambda=params->lambda[brand]; - mu=params->mu[brand]; - exchange->chi=1.0; + moldyn->atom=(t_atom *)malloc(size); + if(moldyn->atom==NULL) { + perror("[moldyn] load save file malloc (atoms)"); + return -1; } - else { - S=params->Smixed; - S2=params->S2mixed; - R=params->Rmixed; - A=params->Amixed; - B=params->Bmixed; - lambda=params->lambda_m; - mu=params->mu_m; - params->exchange.chi=params->chi; - } - - /* dist_ij, d_ij */ - v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij2=v3_absolute_square(&dist_ij); - - /* if d_ij2 > S2 => no force & potential energy contribution */ - if(d_ij2>S2) - return 0; - /* now we will need the distance */ - //d_ij=v3_norm(&dist_ij); - d_ij=sqrt(d_ij2); - - /* save for use in 3bp */ - exchange->d_ij=d_ij; - exchange->d_ij2=d_ij2; - exchange->dist_ij=dist_ij; - - /* more constants */ - exchange->beta_j=&(params->beta[brand]); - exchange->n_j=&(params->n[brand]); - exchange->c_j=&(params->c[brand]); - exchange->d_j=&(params->d[brand]); - exchange->h_j=&(params->h[brand]); - if(brand==ai->brand) { - exchange->betajnj=exchange->betaini; - exchange->cj2=exchange->ci2; - exchange->dj2=exchange->di2; - exchange->cj2dj2=exchange->ci2di2; - } - else { - exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j)); - exchange->cj2=params->c[brand]*params->c[brand]; - exchange->dj2=params->d[brand]*params->d[brand]; - exchange->cj2dj2=exchange->cj2/exchange->dj2; + while(size) { + cnt=read(fd,moldyn->atom,size); + if(cnt<0) { + perror("[moldyn] load save file read (atoms)"); + return cnt; + } + size-=cnt; } - /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */ - f_r=A*exp(-lambda*d_ij); - df_r=lambda*f_r/d_ij; - - /* f_a, df_a calc (again, same for ij and ji) | save for later use! */ - exchange->f_a=-B*exp(-mu*d_ij); - exchange->df_a=mu*exchange->f_a/d_ij; + // hooks etc ... - /* f_c, df_c calc (again, same for ij and ji) */ - if(d_ij r > R */ - exchange->d_ij_between_rs=1; - } + return 0; +} - /* add forces of 2bp (ij, ji) contribution - * dVij = dVji and we sum up both: no 1/2) */ - v3_add(&(ai->f),&(ai->f),&force); +int moldyn_load(t_moldyn *moldyn) { - /* virial */ - ai->virial.xx-=force.x*dist_ij.x; - ai->virial.yy-=force.y*dist_ij.y; - ai->virial.zz-=force.z*dist_ij.z; - ai->virial.xy-=force.x*dist_ij.y; - ai->virial.xz-=force.x*dist_ij.z; - ai->virial.yz-=force.y*dist_ij.z; + // later ... -#ifdef DEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVij, dVji (2bp) contrib:\n"); - printf("%f | %f\n",force.x,ai->f.x); - printf("%f | %f\n",force.y,ai->f.y); - printf("%f | %f\n",force.z,ai->f.z); -} -#endif -#ifdef VDEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVij, dVji (2bp) contrib:\n"); - printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx); - printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy); - printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz); + return 0; } -#endif - /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */ - moldyn->energy+=(0.5*f_r*f_c); +/* + * post processing functions + */ - /* save for use in 3bp */ - exchange->f_c=f_c; - exchange->df_c=df_c; +int get_line(int fd,char *line,int max) { - /* enable the run of 3bp function and 2bp post processing */ - exchange->run3bp=1; - exchange->run2bp_post=1; + int count,ret; - /* reset 3bp sums */ - exchange->zeta_ij=0.0; - exchange->zeta_ji=0.0; - v3_zero(&(exchange->dzeta_ij)); - v3_zero(&(exchange->dzeta_ji)); + count=0; - return 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; + } } -/* tersoff 2 body post part */ +int pair_correlation_init(t_moldyn *moldyn,double dr) { -int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + return 0; +} - /* - * here we have to allow for the 3bp sums - * - * that is: - * - zeta_ij, dzeta_ij - * - zeta_ji, dzeta_ji - * - * to compute the 3bp contribution to: - * - Vij, dVij - * - dVji - * - */ +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { - t_tersoff_mult_params *params; - t_tersoff_exchange *exchange; + int slots; + double *stat; + int i,j; + t_linkcell *lc; +#ifdef STATIC_LISTS + int *neighbour[27]; + int p; +#else + t_list neighbour[27]; +#endif + t_atom *itom,*jtom; + t_list *this; + unsigned char bc; + t_3dvec dist; + double d; + //double norm; + int o,s; + unsigned char ibrand; - t_3dvec force,temp; - t_3dvec *dist_ij; - double b,db,tmp; - double f_c,df_c,f_a,df_a; - double chi,ni,betaini,nj,betajnj; - double zeta; + lc=&(moldyn->lc); - params=moldyn->pot2b_params; - exchange=&(params->exchange); + slots=moldyn->cutoff/dr; + o=2*slots; - /* we do not run if f_c_ij was detected to be 0! */ - if(!(exchange->run2bp_post)) - return 0; + 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); - f_c=exchange->f_c; - df_c=exchange->df_c; - f_a=exchange->f_a; - df_a=exchange->df_a; - betaini=exchange->betaini; - betajnj=exchange->betajnj; - ni=*(exchange->n_i); - nj=*(exchange->n_j); - chi=exchange->chi; - dist_ij=&(exchange->dist_ij); - - /* Vij and dVij */ - zeta=exchange->zeta_ij; - if(zeta==0.0) { - moldyn->debug++; /* just for debugging ... */ - b=chi; - v3_scale(&force,dist_ij,df_a*b*f_c); + if(ptr!=NULL) { + stat=(double *)ptr; } else { - tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ - b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ - db=chi*pow(b,-1.0/(2*ni)-1); /* x(...)^(-1/2n - 1) */ - b=db*b; /* b_ij */ - db*=-0.5*tmp; /* db_ij */ - v3_scale(&force,&(exchange->dzeta_ij),f_a*db); - v3_scale(&temp,dist_ij,df_a*b); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,f_c); - } - v3_scale(&temp,dist_ij,df_c*b*f_a); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,-0.5); - - /* add force */ - v3_add(&(ai->f),&(ai->f),&force); + stat=(double *)malloc(3*slots*sizeof(double)); + if(stat==NULL) { + perror("[moldyn] pair correlation malloc"); + return -1; + } + } - /* virial */ - ai->virial.xx-=force.x*dist_ij->x; - ai->virial.yy-=force.y*dist_ij->y; - ai->virial.zz-=force.z*dist_ij->z; - ai->virial.xy-=force.x*dist_ij->y; - ai->virial.xz-=force.x*dist_ij->z; - ai->virial.yz-=force.y*dist_ij->z; + memset(stat,0,3*slots*sizeof(double)); -#ifdef DEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVij (3bp) contrib:\n"); - printf("%f | %f\n",force.x,ai->f.x); - printf("%f | %f\n",force.y,ai->f.y); - printf("%f | %f\n",force.z,ai->f.z); -} -#endif -#ifdef VDEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVij (3bp) contrib:\n"); - printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); - printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); - printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); -} -#endif + link_cell_init(moldyn,VERBOSE); - /* add energy of 3bp sum */ - moldyn->energy+=(0.5*f_c*b*f_a); + itom=moldyn->atom; + + for(i=0;icount;i++) { + /* neighbour indexing */ + link_cell_neighbour_index(moldyn, + (itom[i].r.x+moldyn->dim.x/2)/lc->x, + (itom[i].r.y+moldyn->dim.y/2)/lc->x, + (itom[i].r.z+moldyn->dim.z/2)/lc->x, + neighbour); - /* dVji */ - zeta=exchange->zeta_ji; - if(zeta==0.0) { - moldyn->debug++; - b=chi; - v3_scale(&force,dist_ij,df_a*b*f_c); - } - else { - tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ - b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ - db=chi*pow(b,-1.0/(2*nj)-1); /* x(...)^(-1/2n - 1) */ - b=db*b; /* b_ij */ - db*=-0.5*tmp; /* db_ij */ - v3_scale(&force,&(exchange->dzeta_ji),f_a*db); - v3_scale(&temp,dist_ij,df_a*b); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,f_c); - } - v3_scale(&temp,dist_ij,df_c*b*f_a); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,-0.5); - - /* add force */ - v3_add(&(ai->f),&(ai->f),&force); - - /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ -// TEST ... with a minus instead - ai->virial.xx-=force.x*dist_ij->x; - ai->virial.yy-=force.y*dist_ij->y; - ai->virial.zz-=force.z*dist_ij->z; - ai->virial.xy-=force.x*dist_ij->y; - ai->virial.xz-=force.x*dist_ij->z; - ai->virial.yz-=force.y*dist_ij->z; + /* brand of atom i */ + ibrand=itom[i].brand; + + for(j=0;j<27;j++) { -#ifdef DEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVji (3bp) contrib:\n"); - printf("%f | %f\n",force.x,ai->f.x); - printf("%f | %f\n",force.y,ai->f.y); - printf("%f | %f\n",force.z,ai->f.z); -} -#endif -#ifdef VDEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVji (3bp) contrib:\n"); - printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); - printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); - printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); -} -#endif + bc=(jdnlc)?0:1; - return 0; -} +#ifdef STATIC_LISTS + p=0; -/* tersoff 3 body part */ - -int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { - - t_tersoff_mult_params *params; - t_tersoff_exchange *exchange; - t_3dvec dist_ij,dist_ik,dist_jk; - t_3dvec temp1,temp2; - t_3dvec *dzeta; - double R,S,S2,s_r; - double B,mu; - double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2; - double rr,dd; - double f_c,df_c; - double f_c_ik,df_c_ik,arg; - double f_c_jk; - double n,c,d,h; - double c2,d2,c2d2; - double cos_theta,d_costheta1,d_costheta2; - double h_cos,d2_h_cos2; - double frac,g,zeta,chi; - double tmp; - int brand; - - params=moldyn->pot3b_params; - exchange=&(params->exchange); - - if(!(exchange->run3bp)) - return 0; + while(neighbour[j][p]!=0) { - /* - * calc of 3bp contribution of V_ij and dV_ij/ji/jk & - * 2bp contribution of dV_jk - * - * for Vij and dV_ij we still need: - * - b_ij, db_ij (zeta_ij) - * - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk - * - * for dV_ji we still need: - * - b_ji, db_ji (zeta_ji) - * - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik - * - * for dV_jk we need: - * - f_c_jk - * - f_a_jk - * - db_jk (zeta_jk) - * - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki - * - */ + jtom=&(moldyn->atom[neighbour[j][p]]); + p++; +#else + this=&(neighbour[j]); + list_reset_f(this); - /* - * get exchange data - */ + if(this->start==NULL) + continue; - /* dist_ij, d_ij - this is < S_ij ! */ - dist_ij=exchange->dist_ij; - d_ij=exchange->d_ij; - d_ij2=exchange->d_ij2; + do { - /* f_c_ij, df_c_ij (same for ji) */ - f_c=exchange->f_c; - df_c=exchange->df_c; + jtom=this->current->data; +#endif + /* only count pairs once, + * skip same atoms */ + if(itom[i].tag>=jtom->tag) + continue; - /* - * calculate unknown values now ... - */ + /* + * pair correlation calc + */ - /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */ + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom[i].r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); - /* dist_ik, d_ik */ - v3_sub(&dist_ik,&(ak->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ik); - d_ik2=v3_absolute_square(&dist_ik); + /* ignore if greater or equal cutoff */ + if(d>=moldyn->cutoff_square) + continue; - /* ik constants */ - brand=ai->brand; - if(brand==ak->brand) { - R=params->R[brand]; - S=params->S[brand]; - S2=params->S2[brand]; - } - else { - R=params->Rmixed; - S=params->Smixed; - S2=params->S2mixed; - } - - /* zeta_ij/dzeta_ij contribution only for d_ik < S */ - if(d_ik2n_i); - c=*(exchange->c_i); - d=*(exchange->d_i); - h=*(exchange->h_i); - c2=exchange->ci2; - d2=exchange->di2; - c2d2=exchange->ci2di2; - - /* cosine of theta_ijk by scalaproduct */ - rr=v3_scalar_product(&dist_ij,&dist_ik); - dd=d_ij*d_ik; - cos_theta=rr/dd; - - /* d_costheta */ - tmp=1.0/dd; - d_costheta1=cos_theta/d_ij2-tmp; - d_costheta2=cos_theta/d_ik2-tmp; - - /* some usefull values */ - h_cos=(h-cos_theta); - d2_h_cos2=d2+(h_cos*h_cos); - frac=c2/(d2_h_cos2); - - /* g(cos_theta) */ - g=1.0+c2d2-frac; - - /* d_costheta_ij and dg(cos_theta) - needed in any case! */ - v3_scale(&temp1,&dist_ij,d_costheta1); - v3_scale(&temp2,&dist_ik,d_costheta2); - v3_add(&temp1,&temp1,&temp2); - v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ - - /* f_c_ik & df_c_ik + {d,}zeta contribution */ - dzeta=&(exchange->dzeta_ij); - if(d_ik f_c_ik=1.0; - // => df_c_ik=0.0; of course we do not set this! - - /* zeta_ij */ - exchange->zeta_ij+=g; - - /* dzeta_ij */ - v3_add(dzeta,dzeta,&temp1); - } - else { - /* {d,}f_c_ik */ - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); - - /* zeta_ij */ - exchange->zeta_ij+=f_c_ik*g; - - /* dzeta_ij */ - v3_scale(&temp1,&temp1,f_c_ik); - v3_scale(&temp2,&dist_ik,g*df_c_ik); - v3_add(&temp1,&temp1,&temp2); - v3_add(dzeta,dzeta,&temp1); + /* fill the slots */ + d=sqrt(d); + s=(int)(d/dr); + + /* should never happen but it does 8) - + * related to -ffloat-store problem! */ + if(s>=slots) s=slots-1; + + 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; + } +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } } - /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */ - - /* dist_jk, d_jk */ - v3_sub(&dist_jk,&(ak->r),&(aj->r)); - if(bc) check_per_bound(moldyn,&dist_jk); - d_jk2=v3_absolute_square(&dist_jk); + /* 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; + } + */ - /* jk constants */ - brand=aj->brand; - if(brand==ak->brand) { - R=params->R[brand]; - S=params->S[brand]; - S2=params->S2[brand]; - B=params->B[brand]; - mu=params->mu[brand]; - chi=1.0; + if(ptr==NULL) { + /* todo: store/print pair correlation function */ + free(stat); } - else { - R=params->Rmixed; - S=params->Smixed; - S2=params->S2mixed; - B=params->Bmixed; - mu=params->mu_m; - chi=params->chi; - } - - /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */ - if(d_jk2n_j); - c=*(exchange->c_j); - d=*(exchange->d_j); - h=*(exchange->h_j); - c2=exchange->cj2; - d2=exchange->dj2; - c2d2=exchange->cj2dj2; - - /* cosine of theta_jik by scalaproduct */ - rr=-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji */ - dd=d_ij*d_jk; - cos_theta=rr/dd; - - /* d_costheta */ - d_costheta1=1.0/dd; - d_costheta2=cos_theta/d_ij2; - - /* some usefull values */ - h_cos=(h-cos_theta); - d2_h_cos2=d2+(h_cos*h_cos); - frac=c2/(d2_h_cos2); - - /* g(cos_theta) */ - g=1.0+c2d2-frac; - - /* d_costheta_jik and dg(cos_theta) - needed in any case! */ - v3_scale(&temp1,&dist_jk,d_costheta1); - v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */ - //v3_add(&temp1,&temp1,&temp2); - v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */ - v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ - - /* store dg in temp2 and use it for dVjk later */ - v3_copy(&temp2,&temp1); - - /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ - dzeta=&(exchange->dzeta_ji); - if(d_jkzeta_ji+=g; - - /* dzeta_ji */ - v3_add(dzeta,dzeta,&temp1); - } - else { - /* f_c_jk */ - s_r=S-R; - arg=M_PI*(d_jk-R)/s_r; - f_c_jk=0.5+0.5*cos(arg); - - /* zeta_ji */ - exchange->zeta_ji+=f_c_jk*g; - - /* dzeta_ji */ - v3_scale(&temp1,&temp1,f_c_jk); - v3_add(dzeta,dzeta,&temp1); - } - /* dV_jk stuff | add force contribution on atom i immediately */ - if(exchange->d_ij_between_rs) { - zeta=f_c*g; - v3_scale(&temp1,&temp2,f_c); - v3_scale(&temp2,&dist_ij,df_c*g); - v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */ - } - else { - zeta=g; - // dzeta_jk is simply dg, which is stored in temp2 - } - /* betajnj * zeta_jk ^ nj-1 */ - tmp=exchange->betajnj*pow(zeta,(n-1.0)); - tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp; - v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); - v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ - /* scaled with 0.5 ^ */ - - /* virial */ - ai->virial.xx-=temp2.x*dist_jk.x; - ai->virial.yy-=temp2.y*dist_jk.y; - ai->virial.zz-=temp2.z*dist_jk.z; - ai->virial.xy-=temp2.x*dist_jk.y; - ai->virial.xz-=temp2.x*dist_jk.z; - ai->virial.yz-=temp2.y*dist_jk.z; + free(moldyn->atom); -#ifdef DEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVjk (3bp) contrib:\n"); - printf("%f | %f\n",temp2.x,ai->f.x); - printf("%f | %f\n",temp2.y,ai->f.y); - printf("%f | %f\n",temp2.z,ai->f.z); -} -#endif -#ifdef VDEBUG -if(ai==&(moldyn->atom[0])) { - printf("dVjk (3bp) contrib:\n"); - printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx); - printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy); - printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz); + link_cell_shutdown(moldyn); + + return 0; } -#endif - } +int analyze_bonds(t_moldyn *moldyn) { + + + return 0; } - /* - * debugging / critical check functions + * visualization code */ -int moldyn_bc_check(t_moldyn *moldyn) { +int visual_init(t_moldyn *moldyn,char *filebase) { - t_atom *atom; - t_3dvec *dim; - int i; - double x; - u8 byte; - int j,k; + 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; +#ifdef STATIC_LISTS + int *neighbour[27]; + int p; +#else + t_list neighbour[27]; +#endif + 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; - dim=&(moldyn->dim); - x=dim->x/2; + 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++) { - if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) { - printf("FATAL: atom %d: x: %.20f (%.20f)\n", - i,atom[i].r.x,dim->x/2); - printf("diagnostic:\n"); - printf("-----------\natom.r.x:\n"); - for(j=0;j<8;j++) { - memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1); - for(k=0;k<8;k++) - printf("%d%c", - ((byte)&(1<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 } - if(atom[i].r.x==x) printf("the same!\n"); - else printf("different!\n"); +#else + } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT); +#endif } - if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2) - printf("FATAL: atom %d: y: %.20f (%.20f)\n", - i,atom[i].r.y,dim->y/2); - if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2) - printf("FATAL: atom %d: z: %.20f (%.20f)\n", - i,atom[i].r.z,dim->z/2); } + // 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; } +