#include <time.h>
#include <math.h>
+#include <fpu_control.h>
+
+#ifdef PARALLEL
+#include <omp.h>
+#endif
+
+#if defined PTHREADS || defined VISUAL_THREAD
+#include <pthread.h>
+#endif
+
#include "moldyn.h"
#include "report/report.h"
#include "potentials/tersoff.h"
#endif
+/* 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
*/
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;
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;i<moldyn->count;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));
*/
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;
void *ptr;
t_atom *atom;
char name[16];
+#ifdef PTHREADS
+ pthread_mutex_t *mutex;
+#endif
new=a*b*c;
count=moldyn->count;
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;
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:
atom[ret].tag=count+ret;
check_per_bound(moldyn,&(atom[ret].r));
atom[ret].r_0=atom[ret].r;
+#ifdef PTHREADS
+ pthread_mutex_init(&(mutex[ret]),NULL);
+#endif
}
/* update total system mass */
}
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 */
}
/* 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)
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;i<a;i++) {
r.y=o.y;
for(j=0;j<b;j++) {
r.z=o.z;
for(k=0;k<c;k++) {
+ switch(p_type) {
+ case PART_INSIDE_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;
+ 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;
}
/* 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)
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;i<a;i++) {
r.z=o.z;
for(k=0;k<c;k++) {
/* first atom */
+ switch(p_type) {
+ case PART_INSIDE_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;
- 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;
}
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;
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;
}
scale*=2.0;
else
if(moldyn->pt_scale&T_SCALE_BERENDSEN)
- scale=1.0+(scale-1.0)*moldyn->timestep/moldyn->t_tc;
+ scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc;
scale=sqrt(scale);
/* velocity scaling */
/* scaling factor */
if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
- scale=(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->timestep;
- scale=pow(1.0-scale,ONE_THIRD);
+ 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);
int link_cell_init(t_moldyn *moldyn,u8 vol) {
t_linkcell *lc;
+#ifndef LOWMEM_LISTS
int i;
+#endif
lc=&(moldyn->lc);
#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->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);
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;i<lc->cells;i++)
list_init_f(&(lc->subcell[i]));
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;i<lc->cells;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
#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) {
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)
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
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;
}
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
}
}
}
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;i<lc->cells;i++) {
#ifdef STATIC_LISTS
free(lc->subcell[i]);
list_destroy_f(&(lc->subcell[i]));
#endif
}
+#endif
free(lc->subcell);
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;
/* 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");
}
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);
}
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,
/* copy over time */
t1=t2;
- //}
+ }
/* increase absolute time */
moldyn->time+=moldyn->tau;
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;i<count;i++) {
/* check whether fixed atom */
int *neighbour_i[27];
int p,q;
t_atom *atom;
+#elif LOWMEM_LISTS
+ int neighbour_i[27];
+ int p,q;
#else
t_list neighbour_i[27];
t_list neighbour_i2[27];
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;i<count;i++) {
/* reset force */
#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];
- if(jtom==&(itom[i]))
- continue;
+ while(p!=-1) {
- if((jtom->attr&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);
do {
jtom=this->current->data;
+#endif
if(jtom==&(itom[i]))
continue;
jtom,
bc_ij);
}
+#ifdef STATIC_LISTS
+ }
+#elif LOWMEM_LISTS
+ }
+#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
/* 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
#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);
#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);
jtom,
ktom,
bc_ik|bc_ij);
+
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(that)!=\
L_NO_NEXT_ELEMENT);
#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);
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(that)!=\
L_NO_NEXT_ELEMENT);
}
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
#endif
/* some postprocessing */
+#ifdef PARALLEL
+ #pragma omp parallel for
+#endif
for(i=0;i<count;i++) {
/* calculate global virial */
moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
#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);
itom=moldyn->atom;
#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);
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
}
+/*
+ * 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=(j<lc->dnlc)?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
*/
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_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;
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
atom[i].ekin);
// bonds between atoms
+#ifndef VISUAL_THREAD
process_2b_bonds(moldyn,&vb,visual_bonds_process);
+#endif
// boundaries
if(dim.x) {
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;
}
projects / physik / posic.git / blobdiff