X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=8451c4436f43db2213bec83677c359fb4062f133;hb=787848dc9bf3daa2ccb76296de904688cc504a45;hp=832d171be44b345c61a6f05306bfe1b12e9c8c3d;hpb=3592f3b0cf641729566674a7cdfd8eff3f61f59a;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 832d171..8451c44 100644 --- a/moldyn.c +++ b/moldyn.c @@ -17,6 +17,16 @@ #include #include +#include + +#ifdef PARALLEL +#include +#endif + +#if defined PTHREADS || defined VISUAL_THREAD +#include +#endif + #include "moldyn.h" #include "report/report.h" @@ -30,100 +40,18 @@ #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", -}; - -/* -static double pse_mass[]={ - 0, - 0, - 0, - 0, - 0, - 0, - M_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - M_SI, - 0, - 0, - 0, - 0, -}; - -static double pse_lc[]={ - 0, - 0, - 0, - 0, - 0, - 0, - LC_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - LC_SI, - 0, - 0, - 0, - 0, -}; -*/ +/* pse */ +#define PSE_NAME +#define PSE_COL +#include "pse.h" +#undef PSE_NAME +#undef PSE_COL + +#ifdef PTHREADS +/* global mutexes */ +pthread_mutex_t *amutex; +pthread_mutex_t emutex; +#endif /* * the moldyn functions @@ -133,6 +61,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; @@ -141,13 +72,28 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { rand_init(&(moldyn->random),NULL,1); moldyn->random.status|=RAND_STAT_VERBOSE; +#ifdef PTHREADS + pthread_mutex_init(&emutex,NULL); +#endif + return 0; } int moldyn_shutdown(t_moldyn *moldyn) { +#ifdef PTHREADS + int i; +#endif + printf("[moldyn] shutdown\n"); +#ifdef PTHREADS + for(i=0;icount;i++) + pthread_mutex_destroy(&(amutex[i])); + free(amutex); + pthread_mutex_destroy(&emutex); +#endif + moldyn_log_shutdown(moldyn); link_cell_shutdown(moldyn); rand_close(&(moldyn->random)); @@ -208,7 +154,7 @@ int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { moldyn->pt_scale|=ptype; moldyn->p_tc=ptc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] p scaling:\n"); printf(" p: %s",ptype?"yes":"no "); if(ptype) @@ -224,7 +170,7 @@ int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { moldyn->pt_scale|=ttype; moldyn->t_tc=ttc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] t scaling:\n"); printf(" t: %s",ttype?"yes":"no "); if(ttype) @@ -315,7 +261,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; @@ -567,18 +513,27 @@ 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,t_3dvec *origin) { + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin, + u8 p_type,t_part_vals *p_vals) { int new,count; int ret; t_3dvec orig; void *ptr; t_atom *atom; + char name[16]; +#ifdef PTHREADS + pthread_mutex_t *mutex; +#endif 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; @@ -592,6 +547,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, moldyn->atom=ptr; atom=&(moldyn->atom[count]); +#ifdef PTHREADS + ptr=realloc(amutex,(count+new)*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + mutex=&(amutex[count]); +#endif + /* no atoms on the boundaries (only reason: it looks better!) */ if(!origin) { orig.x=0.5*lc; @@ -607,19 +572,22 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case CUBIC: set_nn_dist(moldyn,lc); - ret=cubic_init(a,b,c,lc,atom,&orig); + ret=cubic_init(a,b,c,lc,atom,&orig,p_type,p_vals); + 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); + ret=fcc_init(a,b,c,lc,atom,&orig,p_type,p_vals); + 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); + ret=diamond_init(a,b,c,lc,atom,&orig,p_type,p_vals); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -636,7 +604,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 + +#ifdef PTHREADS + ptr=realloc(amutex,(count+1)*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + pthread_mutex_init(&(amutex[count]),NULL); +#endif + atom=moldyn->atom; /* initialize new atom */ @@ -721,12 +711,14 @@ int del_atom(t_moldyn *moldyn,int tag) { } /* cubic init */ -int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { +int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, + u8 p_type,t_part_vals *p_vals) { int count; t_3dvec r; int i,j,k; t_3dvec o; + t_3dvec dist; count=0; if(origin) @@ -734,14 +726,39 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { else v3_zero(&o); + /* shift partition values */ + if(p_type) { + p_vals->p.x+=(a*lc)/2.0; + p_vals->p.y+=(b*lc)/2.0; + p_vals->p.z+=(c*lc)/2.0; + } + r.x=o.x; for(i=0;ip)); + if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { + v3_copy(&(atom[count].r),&r); + count+=1; + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&r,&(p_vals->p)); + if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { v3_copy(&(atom[count].r),&r); count+=1; + } + break; + default: + v3_copy(&(atom[count].r),&r); + count+=1; + break; + } r.z+=lc; } r.y+=lc; @@ -759,12 +776,14 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { } /* fcc lattice init */ -int fcc_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, + u8 p_type,t_part_vals *p_vals) { int count; int i,j,k,l; t_3dvec o,r,n; t_3dvec basis[3]; + t_3dvec dist; count=0; if(origin) @@ -781,6 +800,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { basis[2].y=0.5*lc; basis[2].z=0.5*lc; + /* shift partition values */ + if(p_type) { + p_vals->p.x+=(a*lc)/2.0; + p_vals->p.y+=(b*lc)/2.0; + p_vals->p.z+=(c*lc)/2.0; + } + /* fill up the room */ r.x=o.x; for(i=0;ip)); + if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { v3_copy(&(atom[count].r),&r); count+=1; - r.z+=lc; + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&r,&(p_vals->p)); + if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { + v3_copy(&(atom[count].r),&r); + count+=1; + } + break; + default: + v3_copy(&(atom[count].r),&r); + count+=1; + break; + } /* the three face centered atoms */ for(l=0;l<3;l++) { v3_add(&n,&r,&basis[l]); + switch(p_type) { + case PART_INSIDE_R: + v3_sub(&dist,&n,&(p_vals->p)); + if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { + v3_copy(&(atom[count].r),&r); + count+=1; + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&n,&(p_vals->p)); + if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { + v3_copy(&(atom[count].r),&r); + count+=1; + } + break; + default: v3_copy(&(atom[count].r),&n); count+=1; + break; + } } + r.z+=lc; } r.y+=lc; } @@ -814,12 +876,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { return count; } -int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { +int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, + u8 p_type,t_part_vals *p_vals) { int count; t_3dvec o; - count=fcc_init(a,b,c,lc,atom,origin); + count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals); o.x=0.25*lc; o.y=0.25*lc; @@ -827,7 +890,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { if(origin) v3_add(&o,&o,origin); - count+=fcc_init(a,b,c,lc,&atom[count],&o); + count+=fcc_init(a,b,c,lc,&atom[count],&o,p_type,p_vals); return count; } @@ -959,7 +1022,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 */ @@ -1269,7 +1332,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 { @@ -1364,7 +1427,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); @@ -1379,6 +1444,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 @@ -1389,12 +1456,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); @@ -1418,6 +1489,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])); @@ -1432,22 +1514,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 @@ -1459,7 +1545,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) { @@ -1467,9 +1553,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) @@ -1485,6 +1575,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 @@ -1510,7 +1602,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; @@ -1534,10 +1630,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 } } } @@ -1550,11 +1655,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]); @@ -1563,6 +1676,7 @@ int link_cell_shutdown(t_moldyn *moldyn) { list_destroy_f(&(lc->subcell[i])); #endif } +#endif free(lc->subcell); @@ -1631,6 +1745,17 @@ int moldyn_integrate(t_moldyn *moldyn) { struct timeval t1,t2; //double tp; +#ifdef VISUAL_THREAD + 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; @@ -1665,17 +1790,35 @@ 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 VISUAL_THREAD + 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 + +#ifdef PTHREADS + printf("##################\n"); + printf("# USING PTHREADS #\n"); + printf("##################\n"); +#endif /* tell the world */ printf("[moldyn] integration start, go get a coffee ...\n"); @@ -1703,7 +1846,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); @@ -1756,7 +1903,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); @@ -1771,19 +1918,49 @@ int moldyn_integrate(t_moldyn *moldyn) { } if(a) { if(!(moldyn->total_steps%a)) { +#ifdef VISUAL_THREAD + /* 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); -printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", +printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n", 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)); @@ -1791,7 +1968,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; @@ -1853,7 +2030,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); @@ -1983,6 +2166,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1994,6 +2178,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 @@ -2008,6 +2197,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 @@ -2019,10 +2210,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); @@ -2063,10 +2261,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); @@ -2092,8 +2297,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); @@ -2118,10 +2326,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); @@ -2150,6 +2365,8 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@ -2167,6 +2384,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 @@ -2193,6 +2412,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; @@ -2204,7 +2426,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); } @@ -2405,18 +2627,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++) { @@ -2434,10 +2656,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); @@ -2455,6 +2684,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 @@ -2465,6 +2696,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 */ @@ -2676,8 +2985,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; @@ -2793,7 +3100,11 @@ int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, return 0; } +#ifdef VISUAL_THREAD +void *visual_atoms(void *ptr) { +#else int visual_atoms(t_moldyn *moldyn) { +#endif int i; char file[128+64]; @@ -2802,6 +3113,11 @@ int visual_atoms(t_moldyn *moldyn) { t_visual *v; t_atom *atom; t_vb vb; +#ifdef VISUAL_THREAD + t_moldyn *moldyn; + + moldyn=ptr; +#endif v=&(moldyn->vis); dim.x=v->dim.x; @@ -2811,17 +3127,19 @@ int visual_atoms(t_moldyn *moldyn) { help=(dim.x+dim.y); - sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time); + 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 VISUAL_THREAD return -1; +#endif } /* write the actual data file */ // povray header - dprintf(vb.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 @@ -2835,7 +3153,9 @@ int visual_atoms(t_moldyn *moldyn) { atom[i].ekin); // bonds between atoms +#ifndef VISUAL_THREAD process_2b_bonds(moldyn,&vb,visual_bonds_process); +#endif // boundaries if(dim.x) { @@ -2881,6 +3201,32 @@ int visual_atoms(t_moldyn *moldyn) { close(vb.fd); +#ifdef VISUAL_THREAD + 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; }