X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=edda007d9290ff6cfe1001d4105c3f8a9ab12863;hb=e1080fc0dd66b0cf5b7715c5e99e7a34ac04a8cf;hp=4d736179550d68c7402e89b4412dec3dfc9d40bb;hpb=0b96eb313c9bfec6272b1f8de0d99c4ce26d1686;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 4d73617..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,6 +126,18 @@ int set_cutoff(t_moldyn *moldyn,double cutoff) { return 0; } +int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) { + + moldyn->bondlen[0]=b0*b0; + moldyn->bondlen[1]=b1*b1; + if(bm<0) + moldyn->bondlen[2]=b0*b1; + else + moldyn->bondlen[2]=bm*bm; + + return 0; +} + int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; @@ -218,6 +287,14 @@ int set_potential_params(t_moldyn *moldyn,void *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; +} + int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { strncpy(moldyn->vlsdir,dir,127); @@ -299,7 +376,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { break; case VISUAL_STEP: moldyn->vwrite=timer; - ret=visual_init(&(moldyn->vis),moldyn->vlsdir); + ret=visual_init(moldyn,moldyn->vlsdir); if(ret<0) { printf("[moldyn] visual init failure\n"); return ret; @@ -404,14 +481,16 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { if(moldyn->rfd) { dprintf(moldyn->rfd,report_end); close(moldyn->rfd); - 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); + snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1", + 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; } @@ -421,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; @@ -447,24 +526,33 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom=&(moldyn->atom[count]); /* no atoms on the boundaries (only reason: it looks better!) */ - origin.x=0.5*lc; - origin.y=0.5*lc; - origin.z=0.5*lc; + 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: set_nn_dist(moldyn,lc); - ret=cubic_init(a,b,c,lc,atom,&origin); + ret=cubic_init(a,b,c,lc,atom,&orig); break; case FCC: - v3_scale(&origin,&origin,0.5); + 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,&origin); + ret=fcc_init(a,b,c,lc,atom,&orig); break; case DIAMOND: - v3_scale(&origin,&origin,0.25); + 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,&origin); + ret=diamond_init(a,b,c,lc,atom,&orig); break; default: printf("unknown lattice type (%02x)\n",type); @@ -490,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) { @@ -607,35 +762,6 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { return count; } -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; - - return 0; -} - int destroy_atoms(t_moldyn *moldyn) { if(moldyn->atom) free(moldyn->atom); @@ -693,13 +819,23 @@ 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 */ moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); - moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/moldyn->total_steps; return moldyn->t; } @@ -774,48 +910,133 @@ 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; + /* 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 + - * W = 1/3 sum_i f_i r_i + * with W = 1/3 sum_i f_i r_i (- skipped!) * virial = sum_i f_i r_i * * => P = (2 Ekin + virial) / (3V) */ - v=0.0; - for(i=0;icount;i++) { - virial=&(moldyn->atom[i].virial); - v+=(virial->xx+virial->yy+virial->zz); - } + /* assume up to date virial & up to date kinetic energy */ - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+v; + /* pressure (atom virials) */ + moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); - moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; - /* pressure from 'absolute coordinates' virial */ - virial=&(moldyn->virial); - v=virial->xx+virial->yy+virial->zz; - moldyn->gp=2.0*moldyn->ekin+v; + /* pressure (absolute coordinates) */ + moldyn->gp=2.0*moldyn->ekin+moldyn->gv; moldyn->gp/=(3.0*moldyn->volume); - moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps; 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,dv; + t_3dvec dim; + //t_3dvec *tp; + double u_up,u_down,dv; + double scale,p; t_atom *store; /* @@ -823,13 +1044,11 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * * => p = - dU/dV * - * dV: dx,y,z = 0.001 x,y,z */ - scale=1.00000000000001; -printf("\n\nP-DEBUG:\n"); + scale=0.00001; + dv=8*scale*scale*scale*moldyn->volume; - tp=&(moldyn->tp); store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); @@ -837,59 +1056,44 @@ printf("\n\nP-DEBUG:\n"); } /* 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_dim(moldyn,scale,TRUE,0,0); - scale_atoms(moldyn,scale,TRUE,0,0); - dv=0.00000000000001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z; + /* 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,QUIET); potential_force_calc(moldyn); - tp->x=(moldyn->energy-u)/dv; - p=tp->x*tp->x; + u_up=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - /* derivative with respect to y direction */ - scale_dim(moldyn,scale,0,TRUE,0); - scale_atoms(moldyn,scale,0,TRUE,0); - dv=0.00000000000001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z; + /* 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,QUIET); potential_force_calc(moldyn); - tp->y=(moldyn->energy-u)/dv; - p+=tp->y*tp->y; - - /* restore atomic configuration + dim */ - memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); - moldyn->dim=dim; - - /* derivative with respect to z direction */ - scale_dim(moldyn,scale,0,0,TRUE); - scale_atoms(moldyn,scale,0,0,TRUE); - dv=0.00000000000001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y; - link_cell_shutdown(moldyn); - link_cell_init(moldyn,QUIET); - potential_force_calc(moldyn); - tp->z=(moldyn->energy-u)/dv; - p+=tp->z*tp->z; + 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,QUIET); - return sqrt(p); + return p; } double get_pressure(t_moldyn *moldyn) { @@ -898,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; @@ -911,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; @@ -947,8 +1163,8 @@ int scale_volume(t_moldyn *moldyn) { moldyn->debug=scale; /* scale the atoms and dimensions */ - scale_atoms(moldyn,scale,TRUE,TRUE,TRUE); - scale_dim(moldyn,scale,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); /* visualize dimensions */ if(vdim->x!=0) { @@ -984,8 +1200,10 @@ double e_kin_calc(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; } @@ -1043,25 +1261,57 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { 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"); if(vol) { - printf("[moldyn] initializing linked cells (%d)\n",lc->cells); +#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; } @@ -1071,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); @@ -1079,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); + +#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; @@ -1160,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); @@ -1227,6 +1510,7 @@ int moldyn_integrate(t_moldyn *moldyn) { char dir[128]; double ds; double energy_scale; + struct timeval t1,t2; //double tp; sched=&(moldyn->schedule); @@ -1247,11 +1531,14 @@ int moldyn_integrate(t_moldyn *moldyn) { 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); /* calculate initial forces */ potential_force_calc(moldyn); +#ifdef DEBUG +//return 0; +#endif /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) @@ -1275,13 +1562,17 @@ int moldyn_integrate(t_moldyn *moldyn) { 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++) { @@ -1292,9 +1583,9 @@ int moldyn_integrate(t_moldyn *moldyn) { /* calculate kinetic energy, temperature and pressure */ e_kin_calc(moldyn); temperature_calc(moldyn); + virial_sum(moldyn); pressure_calc(moldyn); - //tp=thermodynamic_pressure_calc(moldyn); -//printf("thermodynamic p: %f %f %f - %f\n",moldyn->tp.x/BAR,moldyn->tp.y/BAR,moldyn->tp.z/BAR,tp/BAR); + average_and_fluctuation_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1304,7 +1595,7 @@ int moldyn_integrate(t_moldyn *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, @@ -1312,7 +1603,7 @@ int moldyn_integrate(t_moldyn *moldyn) { get_total_energy(moldyn)/energy_scale); } if(m) { - if(!(i%m)) { + if(!(moldyn->total_steps%m)) { momentum=get_total_p(moldyn); dprintf(moldyn->mfd, "%f %f %f %f %f\n",moldyn->time, @@ -1321,25 +1612,26 @@ int moldyn_integrate(t_moldyn *moldyn) { } } if(p) { - if(!(i%p)) { + if(!(moldyn->total_steps%p)) { dprintf(moldyn->pfd, "%f %f %f %f %f\n",moldyn->time, - moldyn->p/BAR,moldyn->mean_p/BAR, - moldyn->gp/BAR,moldyn->mean_gp/BAR); + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->gp/BAR,moldyn->gp_avg/BAR); } } if(t) { - if(!(i%t)) { + if(!(moldyn->total_steps%t)) { dprintf(moldyn->tfd, "%f %f %f\n", - moldyn->time,moldyn->t,moldyn->mean_t); + 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)); @@ -1350,22 +1642,28 @@ int moldyn_integrate(t_moldyn *moldyn) { } } if(v) { - if(!(i%v)) { - visual_atoms(&(moldyn->vis),moldyn->time, - moldyn->atom,moldyn->count); + if(!(moldyn->total_steps%v)) { + visual_atoms(moldyn); } } /* display progress */ - if(!(i%10)) { - printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f", - sched->count,i, - moldyn->mean_t, - moldyn->mean_p/BAR, - moldyn->mean_gp/BAR, - moldyn->volume); + //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; @@ -1374,11 +1672,15 @@ int moldyn_integrate(t_moldyn *moldyn) { } /* 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; } @@ -1443,21 +1745,30 @@ 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->virial),0,sizeof(t_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; @@ -1505,14 +1817,33 @@ int potential_force_calc(t_moldyn *moldyn) { 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; - bc_ij=(jcurrent->data; @@ -1527,6 +1858,7 @@ int potential_force_calc(t_moldyn *moldyn) { bc_ij); } } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } } @@ -1537,21 +1869,34 @@ int potential_force_calc(t_moldyn *moldyn) { 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; @@ -1577,17 +1922,24 @@ int potential_force_calc(t_moldyn *moldyn) { for(k=0;k<27;k++) { + bc_ik=(kstart==NULL) continue; - bc_ik=(kcurrent->data; +#endif if(!(ktom->attr&ATOM_ATTR_3BP)) continue; @@ -1603,9 +1955,12 @@ int potential_force_calc(t_moldyn *moldyn) { jtom, ktom, bc_ik|bc_ij); - +#ifdef STATIC_LISTS + } +#else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); +#endif } @@ -1622,17 +1977,24 @@ int potential_force_calc(t_moldyn *moldyn) { for(k=0;k<27;k++) { + bc_ik=(kstart==NULL) continue; - bc_ik=(kcurrent->data; +#endif if(!(ktom->attr&ATOM_ATTR_3BP)) continue; @@ -1649,8 +2011,12 @@ int potential_force_calc(t_moldyn *moldyn) { ktom, bc_ik|bc_ij); +#ifdef STATIC_LISTS + } +#else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); +#endif } @@ -1662,8 +2028,11 @@ int potential_force_calc(t_moldyn *moldyn) { &(itom[i]), jtom,bc_ij); } - +#ifdef STATIC_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } @@ -1677,17 +2046,23 @@ int potential_force_calc(t_moldyn *moldyn) { } #ifdef DEBUG - printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z); + //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;ivirial.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; - moldyn->virial.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; - moldyn->virial.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; - moldyn->virial.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; - moldyn->virial.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; - moldyn->virial.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; + moldyn->gvir.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; @@ -1711,7 +2086,7 @@ 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) { @@ -1793,3 +2168,413 @@ int moldyn_bc_check(t_moldyn *moldyn) { return 0; } + +/* + * restore function + */ + +int moldyn_read_save_file(t_moldyn *moldyn,char *file) { + + int fd; + int cnt,size; + + fd=open(file,O_RDONLY); + if(fd<0) { + perror("[moldyn] load save file open"); + return fd; + } + + size=sizeof(t_moldyn); + + 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); + + moldyn->atom=(t_atom *)malloc(size); + if(moldyn->atom==NULL) { + perror("[moldyn] load save file malloc (atoms)"); + return -1; + } + + while(size) { + cnt=read(fd,moldyn->atom,size); + if(cnt<0) { + perror("[moldyn] load save file read (atoms)"); + return cnt; + } + size-=cnt; + } + + // hooks etc ... + + return 0; +} + +int moldyn_load(t_moldyn *moldyn) { + + // later ... + + 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 pair_correlation_init(t_moldyn *moldyn,double dr) { + + + return 0; +} + +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { + + 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; + + lc=&(moldyn->lc); + + slots=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; + + 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); + + /* brand of atom i */ + ibrand=itom[i].brand; + + for(j=0;j<27;j++) { + + bc=(jdnlc)?0:1; + +#ifdef STATIC_LISTS + p=0; + + while(neighbour[j][p]!=0) { + + jtom=&(moldyn->atom[neighbour[j][p]]); + p++; +#else + this=&(neighbour[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + do { + + jtom=this->current->data; +#endif + /* only count pairs once, + * skip same atoms */ + if(itom[i].tag>=jtom->tag) + continue; + + /* + * pair correlation calc + */ + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom[i].r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater or equal cutoff */ + if(d>=moldyn->cutoff_square) + continue; + + /* 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 + } + } + + /* 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; + } + */ + + if(ptr==NULL) { + /* todo: store/print pair correlation function */ + free(stat); + } + + free(moldyn->atom); + + link_cell_shutdown(moldyn); + + return 0; +} + +int analyze_bonds(t_moldyn *moldyn) { + + + + + return 0; +} + +/* + * visualization code + */ + +int visual_init(t_moldyn *moldyn,char *filebase) { + + strncpy(moldyn->vis.fb,filebase,128); + + return 0; +} + +int visual_atoms(t_moldyn *moldyn) { + + int i,j,fd; + char file[128+64]; + t_3dvec dim; + double help; + t_visual *v; + t_atom *atom; + t_atom *btom; + t_linkcell *lc; +#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; + lc=&(moldyn->lc); + + help=(dim.x+dim.y); + + sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time); + fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(fd<0) { + perror("open visual save file fd"); + return -1; + } + + /* write the actual data file */ + + // povray header + dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help); + + // atomic configuration + for(i=0;icount;i++) { + // atom type, positions, color and kinetic energy + dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element], + atom[i].r.x, + atom[i].r.y, + atom[i].r.z, + pse_col[atom[i].element], + atom[i].ekin); + + /* + * bond detection should usually be done by potential + * functions. brrrrr! EVIL! + * + * todo: potentials need to export a 'find_bonds' function! + */ + + // bonds between atoms + if(!(atom[i].attr&ATOM_ATTR_VB)) + continue; + link_cell_neighbour_index(moldyn, + (atom[i].r.x+moldyn->dim.x/2)/lc->x, + (atom[i].r.y+moldyn->dim.y/2)/lc->y, + (atom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour); + for(j=0;j<27;j++) { + 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 + } + } + + // 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; +} +