X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=d7d89bc2328506715a1d5d82e439267ce41281fa;hb=b4628acf6967b179451bd36f74c2ba312590ecbf;hp=a1bfcaea95c82580dc636f5655649a61d0e8cf8f;hpb=155e1cfea83209d09c2a06ae4fb7f5e1652fc00a;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index a1bfcae..d7d89bc 100644 --- a/moldyn.c +++ b/moldyn.c @@ -17,6 +17,16 @@ #include #include +#include + +#ifdef PARALLEL +#include +#endif + +#ifdef PTHREADS +#include +#endif + #include "moldyn.h" #include "report/report.h" @@ -30,54 +40,12 @@ #include "potentials/tersoff.h" #endif - -/* - * 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", -}; +/* pse */ +#define PSE_NAME +#define PSE_COL +#include "pse.h" +#undef PSE_NAME +#undef PSE_COL /* * the moldyn functions @@ -87,6 +55,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { printf("[moldyn] init\n"); + /* only needed if compiled without -msse2 (float-store prob!) */ + //fpu_set_rtd(); + memset(moldyn,0,sizeof(t_moldyn)); moldyn->argc=argc; @@ -131,24 +102,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) { int set_cutoff(t_moldyn *moldyn,double cutoff) { moldyn->cutoff=cutoff; + moldyn->cutoff_square=cutoff*cutoff; printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff); return 0; } -int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) { - - moldyn->bondlen[0]=b0*b0; - moldyn->bondlen[1]=b1*b1; - if(bm<0) - moldyn->bondlen[2]=b0*b1; - else - moldyn->bondlen[2]=bm*bm; - - return 0; -} - int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; @@ -167,6 +127,38 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { return 0; } +int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { + + moldyn->pt_scale&=(~(P_SCALE_MASK)); + moldyn->pt_scale|=ptype; + moldyn->p_tc=ptc; + + printf("[moldyn] p scaling:\n"); + + printf(" p: %s",ptype?"yes":"no "); + if(ptype) + printf(" | type: %02x | factor: %f",ptype,ptc); + printf("\n"); + + return 0; +} + +int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { + + moldyn->pt_scale&=(~(T_SCALE_MASK)); + moldyn->pt_scale|=ttype; + moldyn->t_tc=ttc; + + printf("[moldyn] t scaling:\n"); + + printf(" t: %s",ttype?"yes":"no "); + if(ttype) + printf(" | type: %02x | factor: %f",ttype,ttc); + printf("\n"); + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); @@ -248,7 +240,7 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func3b_k1=tersoff_mult_3bp_k1; moldyn->func3b_j2=tersoff_mult_3bp_j2; moldyn->func3b_k2=tersoff_mult_3bp_k2; - // missing: check 2b bond func + moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_AM: moldyn->func3b_j1=albe_mult_3bp_j1; @@ -507,11 +499,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, t_3dvec orig; void *ptr; t_atom *atom; + char name[16]; new=a*b*c; count=moldyn->count; /* how many atoms do we expect */ + if(type==NONE) { + new*=1; + printf("[moldyn] WARNING: create 'none' lattice called"); + } if(type==CUBIC) new*=1; if(type==FCC) new*=4; if(type==DIAMOND) new*=8; @@ -541,18 +538,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, case CUBIC: set_nn_dist(moldyn,lc); ret=cubic_init(a,b,c,lc,atom,&orig); + strcpy(name,"cubic"); break; case FCC: 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); + strcpy(name,"fcc"); break; case DIAMOND: 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); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -569,7 +569,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->count+=new; - printf("[moldyn] created lattice with %d atoms\n",new); + printf("[moldyn] created %s lattice with %d atoms\n",name,new); for(ret=0;retatom=ptr; +#ifdef LOWMEM_LISTS + ptr=realloc(moldyn->lc.subcell->list,(count+1)*sizeof(int)); + if(!ptr) { + perror("[moldyn] list realloc (add atom)"); + return -1; + } + moldyn->lc.subcell->list=ptr; +#endif + atom=moldyn->atom; /* initialize new atom */ @@ -892,7 +901,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { scale*=2.0; else if(moldyn->pt_scale&T_SCALE_BERENDSEN) - scale=1.0+(scale-1.0)/moldyn->t_tc; + scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc; scale=sqrt(scale); /* velocity scaling */ @@ -1202,7 +1211,7 @@ int scale_volume(t_moldyn *moldyn) { /* scaling factor */ if(moldyn->pt_scale&P_SCALE_BERENDSEN) { - scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc; + scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau; scale=pow(scale,ONE_THIRD); } else { @@ -1297,7 +1306,9 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { int link_cell_init(t_moldyn *moldyn,u8 vol) { t_linkcell *lc; +#ifndef LOWMEM_LISTS int i; +#endif lc=&(moldyn->lc); @@ -1312,6 +1323,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { #ifdef STATIC_LISTS lc->subcell=malloc(lc->cells*sizeof(int*)); +#elif LOWMEM_LISTS + lc->subcell=malloc(sizeof(t_lowmem_list)); #else lc->subcell=malloc(lc->cells*sizeof(t_list)); #endif @@ -1322,12 +1335,16 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { } if(lc->cells<27) - printf("[moldyn] FATAL: less then 27 subcells!\n"); + printf("[moldyn] FATAL: less then 27 subcells! (%d)\n", + lc->cells); if(vol) { #ifdef STATIC_LISTS printf("[moldyn] initializing 'static' linked cells (%d)\n", lc->cells); +#elif LOWMEM_LISTS + printf("[moldyn] initializing 'lowmem' linked cells (%d)\n", + lc->cells); #else printf("[moldyn] initializing 'dynamic' linked cells (%d)\n", lc->cells); @@ -1351,6 +1368,17 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { i,lc->subcell[0],lc->subcell); */ } +#elif LOWMEM_LISTS + lc->subcell->head=malloc(lc->cells*sizeof(int)); + if(lc->subcell->head==NULL) { + perror("[moldyn] head init (malloc)"); + return -1; + } + lc->subcell->list=malloc(moldyn->count*sizeof(int)); + if(lc->subcell->list==NULL) { + perror("[moldyn] list init (malloc)"); + return -1; + } #else for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); @@ -1365,22 +1393,26 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { int link_cell_update(t_moldyn *moldyn) { int count,i,j,k; - int nx,ny; + int nx,nxy; t_atom *atom; t_linkcell *lc; #ifdef STATIC_LISTS int p; +#elif LOWMEM_LISTS + int p; #endif atom=moldyn->atom; lc=&(moldyn->lc); nx=lc->nx; - ny=lc->ny; + nxy=nx*lc->ny; for(i=0;icells;i++) #ifdef STATIC_LISTS - memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int)); + memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int)); +#elif LOWMEM_LISTS + lc->subcell->head[i]=-1; #else list_destroy_f(&(lc->subcell[i])); #endif @@ -1392,7 +1424,7 @@ int link_cell_update(t_moldyn *moldyn) { #ifdef STATIC_LISTS p=0; - while(lc->subcell[i+j*nx+k*nx*ny][p]!=0) + while(lc->subcell[i+j*nx+k*nxy][p]!=-1) p++; if(p>=MAX_ATOMS_PER_LIST) { @@ -1400,9 +1432,13 @@ int link_cell_update(t_moldyn *moldyn) { return -1; } - lc->subcell[i+j*nx+k*nx*ny][p]=count; + lc->subcell[i+j*nx+k*nxy][p]=count; +#elif LOWMEM_LISTS + p=i+j*nx+k*nxy; + lc->subcell->list[count]=lc->subcell->head[p]; + lc->subcell->head[p]=count; #else - list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]), + list_add_immediate_f(&(lc->subcell[i+j*nx+k*nxy]), &(atom[count])); /* if(j==0&&k==0) @@ -1418,6 +1454,8 @@ int link_cell_update(t_moldyn *moldyn) { int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, #ifdef STATIC_LISTS int **cell +#elif LOWMEM_LISTS + int *cell #else t_list *cell #endif @@ -1443,7 +1481,11 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n", i,nx,j,ny,k,nz); +#ifndef LOWMEM_LISTS cell[0]=lc->subcell[i+j*nx+k*a]; +#else + cell[0]=lc->subcell->head[i+j*nx+k*a]; +#endif for(ci=-1;ci<=1;ci++) { bx=0; x=i+ci; @@ -1467,10 +1509,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, } if(!(ci|cj|ck)) continue; if(bx|by|bz) { +#ifndef LOWMEM_LISTS cell[--count2]=lc->subcell[x+y*nx+z*a]; +#else + cell[--count2]=lc->subcell->head[x+y*nx+z*a]; +#endif + } else { +#ifndef LOWMEM_LISTS cell[count1++]=lc->subcell[x+y*nx+z*a]; +#else + cell[count1++]=lc->subcell->head[x+y*nx+z*a]; +#endif } } } @@ -1483,11 +1534,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, int link_cell_shutdown(t_moldyn *moldyn) { +#ifndef LOWMEM_LISTS int i; +#endif t_linkcell *lc; lc=&(moldyn->lc); +#if LOWMEM_LISTS + free(lc->subcell->head); + free(lc->subcell->list); + +#else + for(i=0;icells;i++) { #ifdef STATIC_LISTS free(lc->subcell[i]); @@ -1496,6 +1555,7 @@ int link_cell_shutdown(t_moldyn *moldyn) { list_destroy_f(&(lc->subcell[i])); #endif } +#endif free(lc->subcell); @@ -1564,6 +1624,17 @@ int moldyn_integrate(t_moldyn *moldyn) { struct timeval t1,t2; //double tp; +#ifdef PTHREADS + u8 first,change; + pthread_t io_thread; + int ret; + t_moldyn md_copy; + t_atom *atom_copy; + + first=1; + change=0; +#endif + sched=&(moldyn->schedule); atom=moldyn->atom; @@ -1581,7 +1652,6 @@ int moldyn_integrate(t_moldyn *moldyn) { /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; /* get current time */ gettimeofday(&t1,NULL); @@ -1599,17 +1669,30 @@ int moldyn_integrate(t_moldyn *moldyn) { printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n"); if(moldyn->cutoff>0.5*moldyn->dim.z) printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n"); - ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass; - if(ds>0.05*moldyn->nnd) + if(moldyn->count) { + ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass; + if(ds>0.05*moldyn->nnd) printf("[moldyn] WARNING: forces too high / tau too small!\n"); + } /* zero absolute time */ - moldyn->time=0.0; - moldyn->total_steps=0; + // should have right values! + //moldyn->time=0.0; + //moldyn->total_steps=0; /* debugging, ignore */ moldyn->debug=0; + /* zero & init moldyn copy */ +#ifdef PTHREADS + memset(&md_copy,0,sizeof(t_moldyn)); + atom_copy=malloc(moldyn->count*sizeof(t_atom)); + if(atom_copy==NULL) { + perror("[moldyn] malloc atom copy (init)"); + return -1; + } +#endif + /* tell the world */ printf("[moldyn] integration start, go get a coffee ...\n"); @@ -1637,7 +1720,11 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); + /* + thermodynamic_pressure_calc(moldyn); + printf("\n\nDEBUG: numeric pressure calc: %f\n\n", + moldyn->tp/BAR); + */ /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); @@ -1690,7 +1777,7 @@ int moldyn_integrate(t_moldyn *moldyn) { } if(s) { if(!(moldyn->total_steps%s)) { - snprintf(dir,128,"%s/s-%07.f.save", + snprintf(dir,128,"%s/s-%08.f.save", moldyn->vlsdir,moldyn->time); fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT, S_IRUSR|S_IWUSR); @@ -1705,12 +1792,41 @@ int moldyn_integrate(t_moldyn *moldyn) { } if(a) { if(!(moldyn->total_steps%a)) { +#ifdef PTHREADS + /* check whether thread has not terminated yet */ + if(!first) { + ret=pthread_join(io_thread,NULL); + } + first=0; + /* prepare and start thread */ + if(moldyn->count!=md_copy.count) { + free(atom_copy); + change=1; + } + memcpy(&md_copy,moldyn,sizeof(t_moldyn)); + if(change) { + atom_copy=malloc(moldyn->count*sizeof(t_atom)); + if(atom_copy==NULL) { + perror("[moldyn] malloc atom copy (change)"); + return -1; + } + } + md_copy.atom=atom_copy; + memcpy(atom_copy,moldyn->atom,moldyn->count*sizeof(t_atom)); + change=0; + ret=pthread_create(&io_thread,NULL,visual_atoms,&md_copy); + if(ret) { + perror("[moldyn] create visual atoms thread\n"); + return -1; + } +#else visual_atoms(moldyn); +#endif } } /* display progress */ - //if(!(moldyn->total_steps%10)) { + if(!(i%10)) { /* get current time */ gettimeofday(&t2,NULL); @@ -1718,6 +1834,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", sched->count,i,moldyn->total_steps, moldyn->t,moldyn->t_avg, moldyn->p/BAR,moldyn->p_avg/BAR, + //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR, moldyn->volume, (int)(t2.tv_sec-t1.tv_sec)); @@ -1725,7 +1842,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", /* copy over time */ t1=t2; - //} + } /* increase absolute time */ moldyn->time+=moldyn->tau; @@ -1787,7 +1904,11 @@ int velocity_verlet(t_moldyn *moldyn) { link_cell_update(moldyn); /* forces depending on chosen potential */ +#ifndef ALBE_FAST potential_force_calc(moldyn); +#else + albe_potential_force_calc(moldyn); +#endif for(i=0;igvir),0,sizeof(t_virial)); /* reset force, site energy and virial of every atom */ +#ifdef PARALLEL + i=omp_get_thread_num(); + #pragma omp parallel for private(virial) +#endif for(i=0;iattr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } - } + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour_i[j]); list_reset_f(this); @@ -1917,6 +2040,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1928,6 +2052,11 @@ int potential_force_calc(t_moldyn *moldyn) { jtom, bc_ij); } +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@ -1942,6 +2071,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* copy the neighbour lists */ #ifdef STATIC_LISTS /* no copy needed for static lists */ +#elif LOWMEM_LISTS + /* no copy needed for lowmem lists */ #else memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); #endif @@ -1953,10 +2084,17 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS p=0; - while(neighbour_i[j][p]!=0) { + while(neighbour_i[j][p]!=-1) { jtom=&(atom[neighbour_i[j][p]]); p++; +#elif LOWMEM_LISTS + p=neighbour_i[j]; + + while(p!=-1) { + + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour_i[j]); list_reset_f(this); @@ -1997,10 +2135,17 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; +#elif LOWMEM_LISTS + q=neighbour_i[k]; + + while(q!=-1) { + + ktom=&(itom[q]); + q=lc->subcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@ -2026,8 +2171,11 @@ int potential_force_calc(t_moldyn *moldyn) { jtom, ktom, bc_ik|bc_ij); + #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@ -2052,10 +2200,17 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; +#elif LOWMEM_LISTS + q=neighbour_i[k]; + + while(q!=-1) { + + ktom=&(itom[q]); + q=lc->subcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@ -2084,6 +2239,8 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@ -2101,6 +2258,8 @@ int potential_force_calc(t_moldyn *moldyn) { } #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@ -2127,6 +2286,9 @@ int potential_force_calc(t_moldyn *moldyn) { #endif /* some postprocessing */ +#ifdef PARALLEL + #pragma omp parallel for +#endif for(i=0;igvir.xx+=itom[i].r.x*itom[i].f.x; @@ -2138,7 +2300,7 @@ int potential_force_calc(t_moldyn *moldyn) { /* check forces regarding the given timestep */ if(v3_norm(&(itom[i].f))>\ - 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) printf("[moldyn] WARNING: pfc (high force: atom %d)\n", i); } @@ -2339,18 +2501,18 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, #ifdef STATIC_LISTS int *neighbour[27]; int p; +#elif LOWMEM_LISTS + int neighbour[27]; + int p; #else t_list neighbour[27]; + t_list *this; #endif u8 bc; t_atom *itom,*jtom; int i,j; - t_list *this; lc=&(moldyn->lc); - - link_cell_init(moldyn,VERBOSE); - itom=moldyn->atom; for(i=0;icount;i++) { @@ -2368,10 +2530,17 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, #ifdef STATIC_LISTS p=0; - while(neighbour[j][p]!=0) { + while(neighbour[j][p]!=-1) { jtom=&(moldyn->atom[neighbour[j][p]]); p++; +#elif LOWMEM_LISTS + p=neighbour[j]; + + while(p!=-1) { + + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour[j]); list_reset_f(this); @@ -2389,6 +2558,8 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@ -2399,6 +2570,84 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, } +/* + * function to find neighboured atoms + */ + +int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom, + int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom, + void *data,u8 bc)) { + + t_linkcell *lc; +#ifdef STATIC_LISTS + int *neighbour[27]; + int p; +#elif LOWMEM_LISTS + int neighbour[27]; + int p; +#else + t_list neighbour[27]; + t_list *this; +#endif + u8 bc; + t_atom *natom; + int j; + + lc=&(moldyn->lc); + + /* neighbour indexing */ + link_cell_neighbour_index(moldyn, + (atom->r.x+moldyn->dim.x/2)/lc->x, + (atom->r.y+moldyn->dim.y/2)/lc->x, + (atom->r.z+moldyn->dim.z/2)/lc->x, + neighbour); + + for(j=0;j<27;j++) { + + bc=(jdnlc)?0:1; + +#ifdef STATIC_LISTS + p=0; + + while(neighbour[j][p]!=-1) { + + natom=&(moldyn->atom[neighbour[j][p]]); + p++; +#elif LOWMEM_LISTS + p=neighbour[j]; + + while(p!=-1) { + + natom=&(moldyn->atom[p]); + p=lc->subcell->list[p]; +#else + this=&(neighbour[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + do { + + natom=this->current->data; +#endif + + /* process bond */ + process(moldyn,atom,natom,data,bc); + +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + } + + return 0; + +} + /* * post processing functions */ @@ -2414,7 +2663,8 @@ int get_line(int fd,char *line,int max) { ret=read(fd,line+count,1); if(ret<=0) return ret; if(line[count]=='\n') { - line[count]='\0'; + memset(line+count,0,max-count-1); + //line[count]='\0'; return count+1; } count+=1; @@ -2609,8 +2859,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) return 0; - d=sqrt(d); - /* now count this bonding ... */ ba=data; @@ -2705,158 +2953,151 @@ int visual_init(t_moldyn *moldyn,char *filebase) { return 0; } +int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { + + t_vb *vb; + + vb=data; + + if(itom->tag>=jtom->tag) + return 0; + + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB)) + dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n", + itom->r.x,itom->r.y,itom->r.z, + jtom->r.x,jtom->r.y,jtom->r.z); + + return 0; +} + +#ifdef PTHREADS +void *visual_atoms(void *ptr) { +#else int visual_atoms(t_moldyn *moldyn) { +#endif - int i,j,fd; + int i; char file[128+64]; t_3dvec dim; double help; t_visual *v; t_atom *atom; - t_atom *btom; - t_linkcell *lc; -#ifdef STATIC_LISTS - int *neighbour[27]; - int p; -#else - t_list neighbour[27]; + t_vb vb; +#ifdef PTHREADS + t_moldyn *moldyn; + + moldyn=ptr; #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) { + sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time); + vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(vb.fd<0) { perror("open visual save file fd"); +#ifndef PTHREADS return -1; +#endif } /* write the actual data file */ // povray header - dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + dprintf(vb.fd,"# [P] %d %08.f <%f,%f,%f>\n", moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help); // atomic configuration - for(i=0;icount;i++) { + for(i=0;icount;i++) // atom type, positions, color and kinetic energy - dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element], - atom[i].r.x, - atom[i].r.y, - atom[i].r.z, - pse_col[atom[i].element], - atom[i].ekin); - - /* - * bond detection should usually be done by potential - * functions. brrrrr! EVIL! - * - * todo: potentials need to export a 'find_bonds' function! - */ - - // bonds between atoms - if(!(atom[i].attr&ATOM_ATTR_VB)) - continue; - link_cell_neighbour_index(moldyn, - (atom[i].r.x+moldyn->dim.x/2)/lc->x, - (atom[i].r.y+moldyn->dim.y/2)/lc->y, - (atom[i].r.z+moldyn->dim.z/2)/lc->z, - neighbour); - for(j=0;j<27;j++) { - bc=jdnlc?0:1; -#ifdef STATIC_LISTS - p=0; - while(neighbour[j][p]!=0) { - btom=&(atom[neighbour[j][p]]); - p++; -#else - list_reset_f(&neighbour[j]); - if(neighbour[j].start==NULL) - continue; - do { - btom=neighbour[j].current->data; -#endif - if(btom==&atom[i]) // skip identical atoms - continue; - //if(btom<&atom[i]) // skip half of them - // continue; - v3_sub(&dist,&(atom[i].r),&(btom->r)); - if(bc) check_per_bound(moldyn,&dist); - d2=v3_absolute_square(&dist); - brand=atom[i].brand; - if(brand==btom->brand) { - if(d2>moldyn->bondlen[brand]) - continue; - } - else { - if(d2>moldyn->bondlen[2]) - continue; - } - dprintf(fd,"# [B] %f %f %f %f %f %f\n", - atom[i].r.x,atom[i].r.y,atom[i].r.z, - btom->r.x,btom->r.y,btom->r.z); -#ifdef STATIC_LISTS - } -#else - } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT); -#endif - } - } - + dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element], + atom[i].r.x, + atom[i].r.y, + atom[i].r.z, + pse_col[atom[i].element], + atom[i].ekin); + + // bonds between atoms + process_2b_bonds(moldyn,&vb,visual_bonds_process); + // boundaries if(dim.x) { - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,-dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,-dim.z/2, -dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,-dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,-dim.z/2, dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, dim.x/2,-dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, -dim.x/2,dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,dim.z/2, dim.x/2,-dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,dim.z/2, dim.x/2,dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, -dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,dim.z/2, -dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,-dim.y/2,dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,dim.z/2, dim.x/2,dim.y/2,-dim.z/2); } - close(fd); + close(vb.fd); + +#ifdef PTHREADS + pthread_exit(NULL); + +} +#else + + return 0; +} +#endif + +/* + * fpu cntrol functions + */ + +// set rounding to double (eliminates -ffloat-store!) +int fpu_set_rtd(void) { + + fpu_control_t ctrl; + + _FPU_GETCW(ctrl); + + ctrl&=~_FPU_EXTENDED; + ctrl|=_FPU_DOUBLE; + + _FPU_SETCW(ctrl); return 0; }