X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1d6e0b7dec9e3032df43f33baa692ba078bb5f7a;hp=16c811ed5e0001d905fb7760239fa0b71a055305;hb=HEAD;hpb=25c22fe95e80365056c6a7fadc548119360ca8ce diff --git a/moldyn.c b/moldyn.c index 16c811e..1d6e0b7 100644 --- a/moldyn.c +++ b/moldyn.c @@ -13,27 +13,95 @@ #include #include #include +#include +#include #include +#include + +#ifdef PARALLEL +#include +#endif + +#if defined PTHREADS || defined VISUAL_THREAD +#include +#endif + #include "moldyn.h" #include "report/report.h" +/* potential includes */ +#include "potentials/harmonic_oscillator.h" +#include "potentials/lennard_jones.h" +#include "potentials/albe.h" +#include "potentials/albe_orig.h" +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else +#include "potentials/tersoff.h" +#endif + +/* pse */ +#define PSE_MASS +#define PSE_NAME +#define PSE_COL +#include "pse.h" +#undef PSE_MASS +#undef PSE_NAME +#undef PSE_COL + +#ifdef PTHREADS +/* global mutexes */ +pthread_mutex_t *amutex; +pthread_mutex_t emutex; +#endif + +/* fully constrained relaxation technique - global pointers */ +u8 crtt; +u8 *constraints; +double *trafo_angle; + +/* + * the moldyn functions + */ + 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; + moldyn->args=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)); @@ -63,6 +131,7 @@ 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); @@ -87,6 +156,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); @@ -122,15 +223,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->vis.dim.z=z; } - moldyn->dv=0.000001*moldyn->volume; - printf("[moldyn] dimensions in A and A^3 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); printf(" y: %f\n",moldyn->dim.y); printf(" z: %f\n",moldyn->dim.z); printf(" volume: %f\n",moldyn->volume); printf(" visualize simulation box: %s\n",visualize?"yes":"no"); - printf(" delta volume (pressure calc): %f\n",moldyn->dv); return 0; } @@ -162,66 +260,55 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func) { - - moldyn->func1b=func; - - return 0; -} - -int set_potential2b(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func2b=func; - - return 0; -} - -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j1=func; - - return 0; -} - -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j2=func; - - return 0; -} - -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j3=func; - - return 0; -} - -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k1=func; - - return 0; -} - -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k2=func; - - return 0; -} - -int set_potential_params(t_moldyn *moldyn,void *params) { +int set_potential(t_moldyn *moldyn,u8 type) { - moldyn->pot_params=params; + switch(type) { + case MOLDYN_POTENTIAL_TM: + //moldyn->func1b=tersoff_mult_1bp; + moldyn->func3b_j1=tersoff_mult_3bp_j1; + moldyn->func3b_k1=tersoff_mult_3bp_k1; + moldyn->func3b_j2=tersoff_mult_3bp_j2; + moldyn->func3b_k2=tersoff_mult_3bp_k2; + moldyn->check_2b_bond=tersoff_mult_check_2b_bond; + break; + case MOLDYN_POTENTIAL_AO: + moldyn->func_j1=albe_orig_mult_3bp_j1; + moldyn->func_j1_k0=albe_orig_mult_3bp_k1; + moldyn->func_j1c=albe_orig_mult_3bp_j2; + moldyn->func_j1_k1=albe_orig_mult_3bp_k2; + moldyn->check_2b_bond=albe_orig_mult_check_2b_bond; + break; + case MOLDYN_POTENTIAL_AM: + moldyn->func_i0=albe_mult_i0; + moldyn->func_j0=albe_mult_i0_j0; + moldyn->func_j0_k0=albe_mult_i0_j0_k0; + moldyn->func_j0e=albe_mult_i0_j1; + moldyn->func_j1=albe_mult_i0_j2; + moldyn->func_j1_k0=albe_mult_i0_j2_k0; + moldyn->func_j1c=albe_mult_i0_j3; + moldyn->check_2b_bond=albe_mult_check_2b_bond; + break; + case MOLDYN_POTENTIAL_HO: + moldyn->func_j0=harmonic_oscillator; + moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; + break; + case MOLDYN_POTENTIAL_LJ: + moldyn->func_j0=lennard_jones; + moldyn->check_2b_bond=lennard_jones_check_2b_bond; + break; + default: + printf("[moldyn] set potential: unknown type %02x\n", + type); + return -1; + } return 0; } -int set_mean_skip(t_moldyn *moldyn,int skip) { +int set_avg_skip(t_moldyn *moldyn,int skip) { printf("[moldyn] skip %d steps before starting average calc\n",skip); - moldyn->mean_skip=skip; + moldyn->avg_skip=skip; return 0; } @@ -301,13 +388,26 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { dprintf(moldyn->tfd,"# temperature log file\n"); printf("temperature (%d)\n",timer); break; + case LOG_VOLUME: + moldyn->vwrite=timer; + snprintf(filename,127,"%s/volume",moldyn->vlsdir); + moldyn->vfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->vfd<0) { + perror("[moldyn] volume log file\n"); + return moldyn->vfd; + } + dprintf(moldyn->vfd,"# volume log file\n"); + printf("volume (%d)\n",timer); + break; case SAVE_STEP: moldyn->swrite=timer; printf("save file (%d)\n",timer); break; case VISUAL_STEP: - moldyn->vwrite=timer; - ret=visual_init(&(moldyn->vis),moldyn->vlsdir); + moldyn->awrite=timer; + ret=visual_init(moldyn,moldyn->vlsdir); if(ret<0) { printf("[moldyn] visual init failure\n"); return ret; @@ -422,7 +522,6 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { moldyn->vlsdir); system(sc); } - if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; } @@ -431,23 +530,48 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { * creating lattice functions */ -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) { +int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin, + t_part_params *p_params,t_defect_params *d_params, + t_offset_params *o_params) { 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; + /* defects */ + if(d_params->type) { + switch(d_params->stype) { + case DEFECT_STYPE_DB_X: + case DEFECT_STYPE_DB_Y: + case DEFECT_STYPE_DB_Z: + case DEFECT_STYPE_DB_R: + new*=2; + break; + default: + printf("[moldyn] WARNING: cl unknown defect\n"); + break; + } + } + /* allocate space for atoms */ ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom)); if(!ptr) { @@ -457,6 +581,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; @@ -471,20 +605,29 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case CUBIC: + if(o_params->use) + v3_add(&orig,&orig,&(o_params->o)); set_nn_dist(moldyn,lc); - ret=cubic_init(a,b,c,lc,atom,&orig); + ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params); + strcpy(name,"cubic"); break; case FCC: if(!origin) v3_scale(&orig,&orig,0.5); + if(o_params->use) + v3_add(&orig,&orig,&(o_params->o)); 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_params,d_params); + strcpy(name,"fcc"); break; case DIAMOND: if(!origin) v3_scale(&orig,&orig,0.25); + if(o_params->use) + v3_add(&orig,&orig,&(o_params->o)); 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_params,d_params); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -493,24 +636,62 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, /* debug */ if(ret!=new) { - printf("[moldyn] creating lattice failed\n"); + printf("[moldyn] creating %s lattice (lc=%f) incomplete\n", + name,lc); + printf(" (ignore for partial lattice creation)\n"); printf(" amount of atoms\n"); printf(" - expected: %d\n",new); printf(" - created: %d\n",ret); + } + + moldyn->count+=ret; + if(ret==new) + printf("[moldyn] created %s lattice with %d atoms\n",name,ret); + + for(new=0;newtype) { + new+=1; + atom[new].element=d_params->element; + atom[new].mass=pse_mass[d_params->element]; + atom[new].attr=d_params->attr; + atom[new].brand=d_params->brand; + atom[new].tag=count+new; + check_per_bound(moldyn,&(atom[new].r)); + atom[new].r_0=atom[new].r; +#ifdef PTHREADS + pthread_mutex_init(&(mutex[new]),NULL); +#endif + } + } + + /* fix allocation */ + ptr=realloc(moldyn->atom,moldyn->count*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (create lattice - alloc fix)"); return -1; } + moldyn->atom=ptr; - moldyn->count+=new; - printf("[moldyn] created lattice with %d atoms\n",new); +// WHAT ABOUT AMUTEX !!!! - for(ret=0;retlc.subcell->list,moldyn->count*sizeof(int)); + if(!ptr) { + perror("[moldyn] list realloc (create lattice)"); + return -1; } + moldyn->lc.subcell->list=ptr; +#endif /* update total system mass */ total_mass_calc(moldyn); @@ -518,13 +699,166 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } +int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, + t_3dvec *r,t_3dvec *v) { + + t_atom *atom; + void *ptr; + int count; + + atom=moldyn->atom; + count=(moldyn->count)++; // asshole style! + + ptr=realloc(atom,(count+1)*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (add atom)"); + return -1; + } + moldyn->atom=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 */ + memset(&(atom[count]),0,sizeof(t_atom)); + atom[count].r=*r; + atom[count].v=*v; + atom[count].element=element; + atom[count].mass=pse_mass[element]; + 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; +#if defined LOWMEM_LISTS || defined PTHREADS + void *ptr; +#endif + + 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); + +#ifdef LOWMEM_LISTS + ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int)); + if(!ptr) { + perror("[moldyn] list realloc (del atom)"); + return -1; + } + moldyn->lc.subcell->list=ptr; + // update lists + link_cell_update(moldyn); +#endif + +#ifdef PTHREADS + ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + pthread_mutex_destroy(&(amutex[moldyn->count+1])); +#endif + + + return 0; +} + +#define set_atom_positions(pos) \ + if(d_params->type) {\ + d_o.x=0; d_o.y=0; d_o.z=0;\ + d_d.x=0; d_d.y=0; d_d.z=0;\ + switch(d_params->stype) {\ + case DEFECT_STYPE_DB_X:\ + d_o.x=d_params->od;\ + d_d.x=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_Y:\ + d_o.y=d_params->od;\ + d_d.y=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_Z:\ + d_o.z=d_params->od;\ + d_d.z=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_R:\ + break;\ + default:\ + printf("[moldyn] WARNING: unknown defect\n");\ + break;\ + }\ + v3_add(&dr,&pos,&d_o);\ + v3_copy(&(atom[count].r),&dr);\ + count+=1;\ + v3_add(&dr,&pos,&d_d);\ + v3_copy(&(atom[count].r),&dr);\ + count+=1;\ + }\ + else {\ + v3_copy(&(atom[count].r),&pos);\ + count+=1;\ + } + /* 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, + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec r; int i,j,k; t_3dvec o; + t_3dvec dist; + t_3dvec p; + t_3dvec d_o; + t_3dvec d_d; + t_3dvec dr; + + p.x=0; p.y=0; p.z=0; count=0; if(origin) @@ -532,14 +866,54 @@ 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_params->type) { + p.x=p_params->p.x+(a*lc)/2.0; + p.y=p_params->p.y+(b*lc)/2.0; + p.z=p_params->p.z+(c*lc)/2.0; + } + r.x=o.x; for(i=0;itype) { + case PART_INSIDE_R: + v3_sub(&dist,&r,&p); + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(r); + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&r,&p); + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(r); + } + break; + case PART_INSIDE_D: + v3_sub(&dist,&r,&p); + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(r); + } + break; + case PART_OUTSIDE_D: + v3_sub(&dist,&r,&p); + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(r); + } + break; + default: + set_atom_positions(r); + break; + } r.z+=lc; } r.y+=lc; @@ -557,12 +931,18 @@ 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, + t_part_params *p_params,t_defect_params *d_params) { int count; int i,j,k,l; t_3dvec o,r,n; t_3dvec basis[3]; + t_3dvec dist; + t_3dvec p; + t_3dvec d_d,d_o,dr; + + p.x=0; p.y=0; p.z=0; count=0; if(origin) @@ -579,6 +959,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_params->type) { + p.x=p_params->p.x+(a*lc)/2.0; + p.y=p_params->p.y+(b*lc)/2.0; + p.z=p_params->p.z+(c*lc)/2.0; + } + /* fill up the room */ r.x=o.x; for(i=0;itype) { + case PART_INSIDE_R: + v3_sub(&dist,&r,&p); + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(r); + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&r,&p); + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(r); + } + break; + case PART_INSIDE_D: + v3_sub(&dist,&r,&p); + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(r); + } + break; + case PART_OUTSIDE_D: + v3_sub(&dist,&r,&p); + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(r); + } + break; + default: + set_atom_positions(r); + break; + } /* the three face centered atoms */ for(l=0;l<3;l++) { v3_add(&n,&r,&basis[l]); - v3_copy(&(atom[count].r),&n); - count+=1; + switch(p_params->type) { + case PART_INSIDE_R: + v3_sub(&dist,&n,&p); + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(n); + } + break; + case PART_OUTSIDE_R: + v3_sub(&dist,&n,&p); + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(n); + } + break; + case PART_INSIDE_D: + v3_sub(&dist,&n,&p); + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(n); + } + break; + case PART_OUTSIDE_D: + v3_sub(&dist,&n,&p); + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(n); + } + break; + default: + set_atom_positions(n); + break; + } } + r.z+=lc; } r.y+=lc; } @@ -612,12 +1065,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, + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec o; - count=fcc_init(a,b,c,lc,atom,origin); + count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params); o.x=0.25*lc; o.y=0.25*lc; @@ -625,43 +1079,11 @@ 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_params,d_params); 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; - - /* update total system mass */ - total_mass_calc(moldyn); - - return 0; -} - int destroy_atoms(t_moldyn *moldyn) { if(moldyn->atom) free(moldyn->atom); @@ -735,13 +1157,9 @@ 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); - - if(moldyn->total_stepsmean_skip) - return 0; - - moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/(moldyn->total_steps+1-moldyn->mean_skip); + if(moldyn->count) + moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); + else moldyn->t=0.0; return moldyn->t; } @@ -795,7 +1213,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 */ @@ -816,77 +1234,139 @@ 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) */ + moldyn->virial=0.0; + moldyn->vir.xx=0.0; + moldyn->vir.yy=0.0; + moldyn->vir.zz=0.0; + moldyn->vir.xy=0.0; + moldyn->vir.xz=0.0; + moldyn->vir.yz=0.0; + for(i=0;icount;i++) { + virial=&(moldyn->atom[i].virial); + moldyn->virial+=(virial->xx+virial->yy+virial->zz); + moldyn->vir.xx+=virial->xx; + moldyn->vir.yy+=virial->yy; + moldyn->vir.zz+=virial->zz; + moldyn->vir.xy+=virial->xy; + moldyn->vir.xz+=virial->xz; + moldyn->vir.yz+=virial->yz; + } + + /* 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); - } - - /* virial sum and mean virial */ - moldyn->virial_sum+=v; - if(moldyn->total_steps>=moldyn->mean_skip) - moldyn->mean_v=moldyn->virial_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); + /* assume up to date virial & up to date kinetic energy */ - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; + /* pressure (atom virials) */ + moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); - if(moldyn->total_steps>=moldyn->mean_skip) { - moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - } - /* pressure from 'absolute coordinates' virial */ - virial=&(moldyn->virial); - v=virial->xx+virial->yy+virial->zz; - moldyn->gp=2.0*moldyn->ekin+v; - moldyn->gp/=(3.0*moldyn->volume); - if(moldyn->total_steps>=moldyn->mean_skip) { - moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - } + //moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx; + //moldyn->px/=moldyn->volume; + //moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy; + //moldyn->py/=moldyn->volume; + //moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz; + //moldyn->pz/=moldyn->volume; + + /* pressure (absolute coordinates) */ + //moldyn->gp=2.0*moldyn->ekin+moldyn->gv; + //moldyn->gp/=(3.0*moldyn->volume); return moldyn->p; } -int energy_fluctuation_calc(t_moldyn *moldyn) { +int average_reset(t_moldyn *moldyn) { + + printf("[moldyn] average reset\n"); + + /* update skip value */ + moldyn->avg_skip=moldyn->total_steps; + + /* kinetic energy */ + moldyn->k_sum=0.0; + moldyn->k2_sum=0.0; + + /* potential energy */ + moldyn->v_sum=0.0; + moldyn->v2_sum=0.0; + + /* temperature */ + moldyn->t_sum=0.0; + + /* virial */ + moldyn->virial_sum=0.0; + //moldyn->gv_sum=0.0; + + /* pressure */ + moldyn->p_sum=0.0; + //moldyn->gp_sum=0.0; + moldyn->tp_sum=0.0; + + return 0; +} + +int average_and_fluctuation_calc(t_moldyn *moldyn) { - if(moldyn->total_stepsmean_skip) + int denom; + + if(moldyn->total_stepsavg_skip) return 0; - /* assume up to date energies */ + 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_mean=moldyn->k_sum/(moldyn->total_steps+1-moldyn->mean_skip); - moldyn->k2_mean=moldyn->k2_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean); + 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_mean=moldyn->v_sum/(moldyn->total_steps+1-moldyn->mean_skip); - moldyn->v2_mean=moldyn->v2_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean); + 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; + moldyn->tp_sum+=moldyn->tp; + moldyn->tp_avg=moldyn->tp_sum/denom; return 0; } @@ -896,13 +1376,13 @@ int get_heat_capacity(t_moldyn *moldyn) { double temp2,ighc; /* averages needed for heat capacity calc */ - if(moldyn->total_stepsmean_skip) + if(moldyn->total_stepsavg_skip) return 0; /* (temperature average)^2 */ - temp2=moldyn->mean_t*moldyn->mean_t; + temp2=moldyn->t_avg*moldyn->t_avg; printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", - moldyn->mean_t); + moldyn->t_avg); /* ideal gas contribution */ ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; @@ -910,25 +1390,27 @@ int get_heat_capacity(t_moldyn *moldyn) { ighc/moldyn->mass*KILOGRAM/JOULE); /* specific heat for nvt ensemble */ - moldyn->c_v_nvt=moldyn->dv2_mean/(K_BOLTZMANN*temp2)+ighc; + 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_mean/(ighc*K_BOLTZMANN*temp2))); + 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_mean,1.5*moldyn->count*K_B2*moldyn->mean_t*moldyn->mean_t*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); +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_up,u_down,dv; - double scale,p; + t_3dvec dim; + //t_3dvec *tp; + double h,dv; + double y0,y1; + double su,sd; t_atom *store; /* @@ -938,54 +1420,56 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * */ - scale=0.00001; - dv=8*scale*scale*scale*moldyn->volume; - + /* store atomic configuration + dimension */ store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); return -1; } - - /* save unscaled potential energy + atom/dim configuration */ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); dim=moldyn->dim; + /* x1, y1 */ + sd=0.00001; + h=(1.0-sd)*(1.0-sd)*(1.0-sd); + su=pow(2.0-h,ONE_THIRD)-1.0; + dv=(1.0-h)*moldyn->volume; + /* scale up dimension and atom positions */ - scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_up=moldyn->energy; + y1=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* scale down dimension and atom positions */ - scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_down=moldyn->energy; + y0=moldyn->energy; /* calculate pressure */ - p=-(u_up-u_down)/dv; -printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR); + moldyn->tp=-(y1-y0)/(2.0*dv); - /* restore atomic configuration + dim */ + /* restore atomic configuration */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - - /* restore energy */ - potential_force_calc(moldyn); - link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); + //potential_force_calc(moldyn); - return p; + /* free store buffer */ + if(store) + free(store); + + return moldyn->tp; } double get_pressure(t_moldyn *moldyn) { @@ -1034,11 +1518,40 @@ int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) { return 0; } +int scale_atoms_ind(t_moldyn *moldyn,double x,double y,double z) { + + int i; + t_3dvec *r; + + for(i=0;icount;i++) { + r=&(moldyn->atom[i].r); + r->x*=x; + r->y*=y; + r->z*=z; + } + + return 0; +} + +int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) { + + t_3dvec *dim; + + dim=&(moldyn->dim); + + dim->x*=x; + dim->y*=y; + dim->z*=z; + + return 0; +} + int scale_volume(t_moldyn *moldyn) { t_3dvec *dim,*vdim; double scale; t_linkcell *lc; + //double sx,sy,sz; vdim=&(moldyn->vis.dim); dim=&(moldyn->dim); @@ -1046,17 +1559,28 @@ 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 { scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); } -moldyn->debug=scale; + + + /* + sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau; + sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau; + sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau; + sx=pow(sx,ONE_THIRD); + sy=pow(sy,ONE_THIRD); + sz=pow(sz,ONE_THIRD); + */ /* scale the atoms and dimensions */ scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + //scale_atoms_ind(moldyn,sx,sy,sz); + //scale_dim_ind(moldyn,sx,sy,sz); /* visualize dimensions */ if(vdim->x!=0) { @@ -1078,6 +1602,9 @@ moldyn->debug=scale; lc->x*=scale; lc->y*=scale; lc->z*=scale; + //lc->x*=sx; + //lc->y*=sx; + //lc->z*=sy; } return 0; @@ -1091,9 +1618,17 @@ double e_kin_calc(t_moldyn *moldyn) { atom=moldyn->atom; moldyn->ekin=0.0; + //moldyn->ekinx=0.0; + //moldyn->ekiny=0.0; + //moldyn->ekinz=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; + //moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x; + //moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y; + //moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z; + } return moldyn->ekin; } @@ -1140,7 +1675,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); @@ -1151,61 +1688,147 @@ 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*)); +#elif LOWMEM_LISTS + lc->subcell=malloc(sizeof(t_lowmem_list)); +#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] FATAL: less then 27 subcells! (%d)\n", + lc->cells); if(vol) { - printf("[moldyn] initializing linked cells (%d)\n",lc->cells); +#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); +#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); + */ + } +#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])); +#endif + + /* update the list */ + link_cell_update(moldyn); - link_cell_update(moldyn); - return 0; } int link_cell_update(t_moldyn *moldyn) { int count,i,j,k; - int nx,ny; + int nx,nxy; t_atom *atom; t_linkcell *lc; - double x,y,z; +#ifdef STATIC_LISTS + int p; +#elif LOWMEM_LISTS + int p; +#endif atom=moldyn->atom; lc=&(moldyn->lc); nx=lc->nx; - ny=lc->ny; - - x=moldyn->dim.x/2; - y=moldyn->dim.y/2; - z=moldyn->dim.z/2; + nxy=nx*lc->ny; for(i=0;icells;i++) +#ifdef STATIC_LISTS + 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 + 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(&(lc->subcell[i+j*nx+k*nx*ny]), + +#ifdef STATIC_LISTS + p=0; + while(lc->subcell[i+j*nx+k*nxy][p]!=-1) + p++; + + if(p>=MAX_ATOMS_PER_LIST) { + printf("[moldyn] FATAL: amount of atoms too high!\n"); + return -1; + } + + 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*nxy]), &(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 +#elif LOWMEM_LISTS + int *cell +#else + t_list *cell +#endif + ) { t_linkcell *lc; int a; @@ -1223,7 +1846,15 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { count2=27; a=nx*ny; + if(i>=nx||j>=ny||k>=nz) + printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n", + i,nx,j,ny,k,nz); + +#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; @@ -1247,10 +1878,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } 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 } } } @@ -1263,13 +1903,28 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { int link_cell_shutdown(t_moldyn *moldyn) { +#ifndef LOWMEM_LISTS int i; +#endif t_linkcell *lc; lc=&(moldyn->lc); - for(i=0;inx*lc->ny*lc->nz;i++) - list_destroy_f(&(moldyn->lc.subcell[i])); +#if LOWMEM_LISTS + free(lc->subcell->head); + free(lc->subcell->list); + +#else + + 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 + } +#endif free(lc->subcell); @@ -1327,7 +1982,7 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { int i; - unsigned int e,m,s,v,p,t; + unsigned int e,m,s,v,p,t,a; t_3dvec momentum; t_moldyn_schedule *sched; t_atom *atom; @@ -1335,8 +1990,20 @@ int moldyn_integrate(t_moldyn *moldyn) { char dir[128]; double ds; double energy_scale; + 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; @@ -1348,51 +2015,73 @@ int moldyn_integrate(t_moldyn *moldyn) { m=moldyn->mwrite; s=moldyn->swrite; v=moldyn->vwrite; + a=moldyn->awrite; p=moldyn->pwrite; t=moldyn->twrite; /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; - /* 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; +//return 0; #endif /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) - printf("[moldyn] warning: cutoff > 0.5 x dim.x\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n"); if(moldyn->cutoff>0.5*moldyn->dim.y) - printf("[moldyn] warning: cutoff > 0.5 x dim.y\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n"); if(moldyn->cutoff>0.5*moldyn->dim.z) - printf("[moldyn] warning: cutoff > 0.5 x dim.z\n"); - ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass; - if(ds>0.05*moldyn->nnd) - printf("[moldyn] warning: forces too high / tau too small!\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n"); + 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"); /* 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++) { @@ -1403,10 +2092,16 @@ return 0; /* calculate kinetic energy, temperature and pressure */ e_kin_calc(moldyn); temperature_calc(moldyn); + virial_sum(moldyn); pressure_calc(moldyn); - energy_fluctuation_calc(moldyn); - //tp=thermodynamic_pressure_calc(moldyn); -//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR); +#ifdef PDEBUG + thermodynamic_pressure_calc(moldyn); + printf("\n\nDEBUG: numeric pressure calc: %f\n\n", + moldyn->tp/BAR); +#endif + + /* calculate fluctuations + averages */ + average_and_fluctuation_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1416,15 +2111,16 @@ return 0; /* check for log & visualization */ if(e) { - if(!(i%e)) + if(!(moldyn->total_steps%e)) dprintf(moldyn->efd, - "%f %f %f %f\n", + "%f %f %f %f %f %f\n", moldyn->time,moldyn->ekin/energy_scale, moldyn->energy/energy_scale, - get_total_energy(moldyn)/energy_scale); + get_total_energy(moldyn)/energy_scale, + moldyn->ekin/EV,moldyn->energy/EV); } 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, @@ -1433,25 +2129,37 @@ return 0; } } 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); + "%f %f %f %f %f %f %f\n",moldyn->time, + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->tp/BAR,moldyn->tp_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(v) { + if(!(moldyn->total_steps%v)) { + dprintf(moldyn->vfd, + "%f %f %f %f %f\n",moldyn->time, + moldyn->volume, + moldyn->dim.x, + moldyn->dim.y, + moldyn->dim.z); } } if(s) { - if(!(i%s)) { - snprintf(dir,128,"%s/s-%07.f.save", + if(!(moldyn->total_steps%s)) { + snprintf(dir,128,"%s/s-%08.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)); @@ -1461,25 +2169,66 @@ return 0; close(fd); } } - if(v) { - if(!(i%v)) { - visual_atoms(&(moldyn->vis),moldyn->time, - moldyn->atom,moldyn->count); + 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 */ +#ifndef PDEBUG 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); +#endif + /* get current time */ + gettimeofday(&t2,NULL); + +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, +#ifndef PDEBUG + moldyn->p/BAR,moldyn->p_avg/BAR, +#else + moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR, +#endif + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); + fflush(stdout); -printf("\n"); -get_heat_capacity(moldyn); + + /* copy over time */ + t1=t2; +#ifndef PDEBUG } +#endif /* increase absolute time */ moldyn->time+=moldyn->tau; @@ -1488,13 +2237,69 @@ get_heat_capacity(moldyn); } /* check for hooks */ - if(sched->count+1total_sched) - if(sched->hook) - sched->hook(moldyn,sched->hook_params); + 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; + + } + + /* writing a final save file! */ + if(s) { + snprintf(dir,128,"%s/s-final.save",moldyn->vlsdir); + 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)); + write(fd,moldyn->atom, + moldyn->count*sizeof(t_atom)); + } + close(fd); + } + /* writing a final visual file! */ + if(a) + visual_atoms(moldyn); + + return 0; +} + +/* basis trafo */ + +#define FORWARD 0 +#define BACKWARD 1 + +int basis_trafo(t_3dvec *r,u8 dir,double z,double x) { + t_3dvec tmp; + + if(dir==FORWARD) { + if(z!=0.0) { + v3_copy(&tmp,r); + r->x=cos(z)*tmp.x-sin(z)*tmp.y; + r->y=sin(z)*tmp.x+cos(z)*tmp.y; + } + if(x!=0.0) { + v3_copy(&tmp,r); + r->y=cos(x)*tmp.y-sin(x)*tmp.z; + r->z=sin(x)*tmp.y+cos(x)*tmp.z; + } + } + else { + if(x!=0.0) { + v3_copy(&tmp,r); + r->y=cos(-x)*tmp.y-sin(-x)*tmp.z; + r->z=sin(-x)*tmp.y+cos(-x)*tmp.z; + } + if(z!=0.0) { + v3_copy(&tmp,r); + r->x=cos(-z)*tmp.x-sin(-z)*tmp.y; + r->y=sin(-z)*tmp.x+cos(-z)*tmp.y; + } } return 0; @@ -1515,12 +2320,47 @@ int velocity_verlet(t_moldyn *moldyn) { tau_square=moldyn->tau_square; for(i=0;ifunc3b_j1==albe_mult_3bp_j1) + // albe_potential_force_calc(moldyn); + //else + potential_force_calc(moldyn); +#else + albe_potential_force_calc(moldyn); +#endif for(i=0;icount; 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 */ +#ifdef PARALLEL + i=omp_get_thread_num(); + #pragma omp parallel for private(virial) +#endif for(i=0;ifunc1b) - moldyn->func1b(moldyn,&(itom[i])); + if(moldyn->func_i0) + moldyn->func_i0(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; @@ -1617,57 +2514,163 @@ int potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; +#ifndef STATIC_LISTS + /* check for later 3 body interaction */ + if(itom[i].attr&ATOM_ATTR_3BP) + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + /* first loop over atoms j */ - if(moldyn->func2b) { + if(moldyn->func_j0) { for(j=0;j<27;j++) { + bc_ij=(jsubcell->list[p]; +#else this=&(neighbour_i[j]); list_reset_f(this); if(this->start==NULL) continue; - bc_ij=(jcurrent->data; +#endif if(jtom==&(itom[i])) continue; + /* reset 3bp run */ + moldyn->run3bp=1; + if((jtom->attr&ATOM_ATTR_2BP)& (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + moldyn->func_j0(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + + /* 3 body potential/force */ + + /* in j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) + continue; + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; + + if(moldyn->func_j0_k0==NULL) + continue; + + /* first loop over atoms k in first j loop */ + 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->func_j0_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS } +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } } - /* 3 body potential/force */ + /* continued 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 */ +#elif LOWMEM_LISTS + /* no copy needed for lowmem 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=(jsubcell->list[p]; +#else this=&(neighbour_i[j]); list_reset_f(this); if(this->start==NULL) continue; - bc_ij=(jcurrent->data; +#endif if(jtom==&(itom[i])) continue; @@ -1678,38 +2681,52 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset 3bp run */ moldyn->run3bp=1; - if(moldyn->func3b_j1) - moldyn->func3b_j1(moldyn, - &(itom[i]), - jtom, - bc_ij); + if(moldyn->func_j1) + moldyn->func_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); - /* in first j loop, 3bp run can be skipped */ + /* in j loop, 3bp run can be skipped */ if(!(moldyn->run3bp)) continue; - /* first loop over atoms k */ - if(moldyn->func3b_k1) { + /* first loop over atoms k in second j loop */ + if(moldyn->func_j1_k0) { for(k=0;k<27;k++) { + bc_ik=(ksubcell->list[q]; +#else that=&(neighbour_i2[k]); list_reset_f(that); if(that->start==NULL) continue; - bc_ik=(kcurrent->data; +#endif if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; @@ -1720,69 +2737,101 @@ int potential_force_calc(t_moldyn *moldyn) { ktom, bc_ik|bc_ij); +#ifdef STATIC_LISTS + } +#elif LOWMEM_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); + /* continued j loop after first k loop */ + if(moldyn->func_j1c) + moldyn->func_j1c(moldyn, + &(itom[i]), + jtom, + bc_ij); /* second loop over atoms k */ - if(moldyn->func3b_k2) { + if(moldyn->func_j1_k1) { for(k=0;k<27;k++) { + bc_ik=(ksubcell->list[q]; +#else that=&(neighbour_i2[k]); list_reset_f(that); if(that->start==NULL) continue; - bc_ik=(kcurrent->data; +#endif if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k2(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); +#endif } } - /* 2bp post function */ - if(moldyn->func3b_j3) { - moldyn->func3b_j3(moldyn, - &(itom[i]), - jtom,bc_ij); + /* finish of j loop after second k loop */ + if(moldyn->func_j1e) { + moldyn->func_j1e(moldyn, + &(itom[i]), + jtom,bc_ij); } - +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif } - + #ifdef DEBUG //printf("\n\n"); #endif @@ -1793,17 +2842,33 @@ 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[DATOM].f.x); + printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y); + printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z); + } #endif - /* calculate global virial */ + /* some postprocessing */ +#ifdef PARALLEL + #pragma omp parallel for +#endif 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; + /* calculate global virial */ + moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x; + moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y; + moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z; + moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x; + moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x; + moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y; + + /* check forces regarding the given timestep */ + if(v3_norm(&(itom[i].f))>\ + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + printf("[moldyn] WARNING: pfc (high force: atom %d)\n", + i); } return 0; @@ -1858,6 +2923,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { return 0; } +int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) { + + double x,y,z; + t_3dvec *dim; + + dim=&(moldyn->dim); + + x=dim->x/2; + y=dim->y/2; + z=dim->z/2; + + if(moldyn->status&MOLDYN_STAT_PBX) { + if(a->r.x>=x) { + a->pbc[0]+=1; + a->r.x-=dim->x; + } + else if(-a->r.x>x) { + a->pbc[0]-=1; + a->r.x+=dim->x; + } + } + if(moldyn->status&MOLDYN_STAT_PBY) { + if(a->r.y>=y) { + a->pbc[1]+=1; + a->r.y-=dim->y; + } + else if(-a->r.y>y) { + a->pbc[1]-=1; + a->r.y+=dim->y; + } + } + if(moldyn->status&MOLDYN_STAT_PBZ) { + if(a->r.z>=z) { + a->pbc[2]+=1; + a->r.z-=dim->z; + } + else if(-a->r.z>z) { + a->pbc[2]-=1; + a->r.z+=dim->z; + } + } + + return 0; +} + /* * debugging / critical check functions */ @@ -1910,6 +3020,255 @@ 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; + int fsize; + int corr; + + fd=open(file,O_RDONLY); + if(fd<0) { + perror("[moldyn] load save file open"); + return fd; + } + + fsize=lseek(fd,0,SEEK_END); + lseek(fd,0,SEEK_SET); + + 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); + + /* correcting possible atom data offset */ + corr=0; + if(fsize!=sizeof(t_moldyn)+size) { + corr=fsize-sizeof(t_moldyn)-size; + printf("[moldyn] WARNING: lsf (illegal file size)\n"); + printf(" modifying offset:\n"); + printf(" - current pos: %d\n",sizeof(t_moldyn)); + printf(" - atom size: %d\n",size); + printf(" - file size: %d\n",fsize); + printf(" => correction: %d\n",corr); + lseek(fd,corr,SEEK_CUR); + } + + moldyn->atom=(t_atom *)malloc(size); + if(moldyn->atom==NULL) { + perror("[moldyn] load save file malloc (atoms)"); + return -1; + } + + while(size) { + cnt=read(fd,moldyn->atom,size); + if(cnt<0) { + perror("[moldyn] load save file read (atoms)"); + return cnt; + } + size-=cnt; + } + +#ifdef PTHREADS + amutex=malloc(moldyn->count*sizeof(pthread_mutex_t)); + if(amutex==NULL) { + perror("[moldyn] load save file (mutexes)"); + return -1; + } + for(cnt=0;cntcount;cnt++) + pthread_mutex_init(&(amutex[cnt]),NULL); +#endif + + // hooks etc ... + + return 0; +} + +int moldyn_free_save_file(t_moldyn *moldyn) { + + free(moldyn->atom); + + return 0; +} + +int moldyn_load(t_moldyn *moldyn) { + + // later ... + + return 0; +} + +/* + * function to find/callback all combinations of 2 body bonds + */ + +int process_2b_bonds(t_moldyn *moldyn,void *data, + int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc)) { + + 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 *itom,*jtom; + int i,j; + + lc=&(moldyn->lc); + 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); + + for(j=0;j<27;j++) { + + bc=(jdnlc)?0:1; + +#ifdef STATIC_LISTS + 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); + + if(this->start==NULL) + continue; + + do { + + jtom=this->current->data; +#endif + + /* process bond */ + process(moldyn,&(itom[i]),jtom,data,bc); + +#ifdef STATIC_LISTS + } +#elif LOWMEM_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + } + } + + return 0; + +} + +/* + * 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 */ @@ -1925,10 +3284,508 @@ 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; } } +int pair_correlation_init(t_moldyn *moldyn,double dr) { + + + return 0; +} + +int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) { + + int i; + t_atom *atom; + t_3dvec dist; + t_3dvec final_r; + double d2; + int a_cnt; + int b_cnt; + + atom=moldyn->atom; + dc[0]=0; + dc[1]=0; + dc[2]=0; + a_cnt=0; + b_cnt=0; + + for(i=0;icount;i++) { + + /* care for pb crossing */ + final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x; + final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y; + final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z; + /* calculate distance */ + v3_sub(&dist,&final_r,&(atom[i].r_0)); + d2=v3_absolute_square(&dist); + + if(atom[i].brand) { + b_cnt+=1; + dc[1]+=d2; + } + else { + a_cnt+=1; + dc[0]+=d2; + } + + dc[2]+=d2; + } + + dc[0]*=(1.0/(6.0*moldyn->time*a_cnt)); + dc[1]*=(1.0/(6.0*moldyn->time*b_cnt)); + dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count)); + + return 0; +} + +int calculate_msd(t_moldyn *moldyn,double *msd) { + + int i; + t_atom *atom; + t_3dvec dist; + t_3dvec final_r; + double d2; + int a_cnt; + int b_cnt; + + atom=moldyn->atom; + msd[0]=0; + msd[1]=0; + msd[2]=0; + a_cnt=0; + b_cnt=0; + + for(i=0;icount;i++) { + + /* care for pb crossing */ + if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) { + printf("[moldyn] msd pb crossings for atom %d\n",i); + printf(" x: %d y: %d z: %d\n", + atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]); + } + final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x; + final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y; + final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z; + /* calculate distance */ + v3_sub(&dist,&final_r,&(atom[i].r_0)); + d2=v3_absolute_square(&dist); + + if(atom[i].brand) { + b_cnt+=1; + msd[1]+=d2; + } + else { + a_cnt+=1; + msd[0]+=d2; + } + + msd[2]+=d2; + } + + msd[0]/=a_cnt; + msd[1]/=b_cnt; + msd[2]/=moldyn->count; + + return 0; +} + +int bonding_analyze(t_moldyn *moldyn,double *cnt) { + + return 0; +} + +int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, + t_atom *jtom,void *data,u8 bc) { + + t_3dvec dist; + double d; + int s; + t_pcc *pcc; + + /* only count pairs once, + * skip same atoms */ + if(itom->tag>=jtom->tag) + return 0; + + /* + * pair correlation calc + */ + + /* get pcc data */ + pcc=data; + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater cutoff */ + if(d>moldyn->cutoff_square) + return 0; + + /* fill the slots */ + d=sqrt(d); + s=(int)(d/pcc->dr); + + /* should never happen but it does 8) - + * related to -ffloat-store problem! */ + if(s>=pcc->o1) { + printf("[moldyn] WARNING: pcc (%d/%d)", + s,pcc->o1); + printf("\n"); + s=pcc->o1-1; + } + + if(itom->brand!=jtom->brand) { + /* mixed */ + pcc->stat[s]+=1; + } + else { + /* type a - type a bonds */ + if(itom->brand==0) + pcc->stat[s+pcc->o1]+=1; + else + /* type b - type b bonds */ + pcc->stat[s+pcc->o2]+=1; + } + + return 0; +} + +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { + + t_pcc pcc; + double norm; + int i; + + pcc.dr=dr; + pcc.o1=moldyn->cutoff/dr; + pcc.o2=2*pcc.o1; + + if(pcc.o1*dr<=moldyn->cutoff) + printf("[moldyn] WARNING: pcc (low #slots)\n"); + + printf("[moldyn] pair correlation calc info:\n"); + printf(" time: %f\n",moldyn->time); + printf(" count: %d\n",moldyn->count); + printf(" cutoff: %f\n",moldyn->cutoff); + printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg); + + if(ptr!=NULL) { + pcc.stat=(double *)ptr; + } + else { + pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double)); + if(pcc.stat==NULL) { + perror("[moldyn] pair correlation malloc"); + return -1; + } + } + + memset(pcc.stat,0,3*pcc.o1*sizeof(double)); + + /* process */ + process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process); + + /* normalization */ + for(i=1;i 2 pi r r dr + // ... and actually it's a constant times r^2 + norm=i*i*dr*dr; + pcc.stat[i]/=norm; + pcc.stat[pcc.o1+i]/=norm; + pcc.stat[pcc.o2+i]/=norm; + } + /* */ + + if(ptr==NULL) { + /* todo: store/print pair correlation function */ + free(pcc.stat); + } + + return 0; +} + +int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { + + t_ba *ba; + t_3dvec dist; + double d; + + if(itom->tag>=jtom->tag) + return 0; + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater or equal cutoff */ + if(d>moldyn->cutoff_square) + return 0; + + /* check for potential bond */ + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + /* now count this bonding ... */ + ba=data; + + /* increase total bond counter + */ + ba->tcnt+=1; + + if(itom->brand==0) + ba->acnt[jtom->tag]+=1; + else + ba->bcnt[jtom->tag]+=1; + + if(jtom->brand==0) + ba->acnt[itom->tag]+=1; + else + ba->bcnt[itom->tag]+=1; + + return 0; +} + +int bond_analyze(t_moldyn *moldyn,double *quality) { + + int qcnt4; + int qcnt3; + int ccnt4; + int ccnt3; + int bcnt; + t_ba ba; + int i; + t_atom *atom; + + ba.acnt=malloc(moldyn->count*sizeof(int)); + if(ba.acnt==NULL) { + perror("[moldyn] bond analyze malloc (a)"); + return -1; + } + memset(ba.acnt,0,moldyn->count*sizeof(int)); + + ba.bcnt=malloc(moldyn->count*sizeof(int)); + if(ba.bcnt==NULL) { + perror("[moldyn] bond analyze malloc (b)"); + return -1; + } + memset(ba.bcnt,0,moldyn->count*sizeof(int)); + + ba.tcnt=0; + qcnt4=0; qcnt3=0; + ccnt4=0; ccnt3=0; + bcnt=0; + + atom=moldyn->atom; + + process_2b_bonds(moldyn,&ba,bond_analyze_process); + + for(i=0;icount;i++) { + if(atom[i].brand==0) { + if((ba.acnt[i]==0)&(ba.bcnt[i]==4)) + qcnt4+=4; + if((ba.acnt[i]==0)&(ba.bcnt[i]==3)) + qcnt3+=3; + } + else { + if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { + qcnt4+=4; + ccnt4+=1; + } + if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) { + qcnt3+=4; + ccnt3+=1; + } + bcnt+=1; + } + } + + if(quality) { + quality[0]=1.0*ccnt4/bcnt; + quality[1]=1.0*ccnt3/bcnt; + } + else { + printf("[moldyn] bond analyze: %f %f\n", + 1.0*ccnt4/bcnt,1.0*ccnt3/bcnt); + } + + return 0; +} + +/* + * visualization code + */ + +int visual_init(t_moldyn *moldyn,char *filebase) { + + strncpy(moldyn->vis.fb,filebase,128); + + 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 VISUAL_THREAD +void *visual_atoms(void *ptr) { +#else +int visual_atoms(t_moldyn *moldyn) { +#endif + + int i; + char file[128+64]; + t_3dvec dim; + double help; + t_visual *v; + t_atom *atom; + t_vb vb; + t_3dvec strain; +#ifdef VISUAL_THREAD + t_moldyn *moldyn; + + moldyn=ptr; +#endif + + v=&(moldyn->vis); + dim.x=v->dim.x; + dim.y=v->dim.y; + dim.z=v->dim.z; + atom=moldyn->atom; + + help=(dim.x+dim.y); + + 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 %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++) { + v3_sub(&strain,&(atom[i].r),&(atom[i].r_0)); + check_per_bound(moldyn,&strain); + // atom type, positions, color and kinetic energy + 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); + sqrt(v3_absolute_square(&strain))); + } + + // bonds between atoms +#ifndef VISUAL_THREAD + process_2b_bonds(moldyn,&vb,visual_bonds_process); +#endif + + // boundaries + if(dim.x) { + 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(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(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(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(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(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(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(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(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(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(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(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(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; +} +