X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=7698c24767b90de44a165602077a46714d86d476;hb=d01a01f5bebe3c2b2d5f5700e5b8b4d42ae41207;hp=a9341941909c452a613ff62877e6b016bbf10571;hpb=d843bfcc38ac4ad1ceffbb82a1aba8096f791b41;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index a934194..7698c24 100644 --- a/moldyn.c +++ b/moldyn.c @@ -41,9 +41,11 @@ #endif /* pse */ +#define PSE_MASS #define PSE_NAME #define PSE_COL #include "pse.h" +#undef PSE_MASS #undef PSE_NAME #undef PSE_COL @@ -512,9 +514,9 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { * creating lattice functions */ -int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, +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, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int new,count; int ret; @@ -538,6 +540,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, 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) { @@ -572,21 +589,21 @@ 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,p_type,p_vals); + 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); set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); - ret=fcc_init(a,b,c,lc,atom,&orig,p_type,p_vals); + 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); set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); - ret=diamond_init(a,b,c,lc,atom,&orig,p_type,p_vals); + ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params); strcpy(name,"diamond"); break; default: @@ -598,7 +615,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, if(ret!=new) { printf("[moldyn] creating %s lattice (lc=%f) incomplete\n", name,lc); - printf(" (ignore in case of partial lattice creation)\n"); + printf(" (ignore for partial lattice creation)\n"); printf(" amount of atoms\n"); printf(" - expected: %d\n",new); printf(" - created: %d\n",ret); @@ -610,7 +627,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, 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 */ @@ -629,6 +659,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->atom=ptr; +// WHAT ABOUT AMUTEX !!!! + #ifdef LOWMEM_LISTS ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int)); if(!ptr) { @@ -644,7 +676,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, +int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, t_3dvec *r,t_3dvec *v) { t_atom *atom; @@ -687,7 +719,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, atom[count].r=*r; atom[count].v=*v; atom[count].element=element; - atom[count].mass=mass; + atom[count].mass=pse_mass[element]; atom[count].brand=brand; atom[count].tag=count; atom[count].attr=attr; @@ -704,6 +736,9 @@ 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; @@ -726,12 +761,69 @@ int del_atom(t_moldyn *moldyn,int tag) { 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, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec r; @@ -739,6 +831,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, 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; @@ -749,10 +844,10 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, v3_zero(&o); /* shift partition values */ - if(p_type) { - p.x=p_vals->p.x+(a*lc)/2.0; - p.y=p_vals->p.y+(b*lc)/2.0; - p.z=p_vals->p.z+(c*lc)/2.0; + 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; @@ -761,42 +856,39 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, for(j=0;jtype) { case PART_INSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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: - v3_copy(&(atom[count].r),&r); - count+=1; + set_atom_positions(r); break; } r.z+=lc; @@ -817,7 +909,7 @@ 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, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int count; int i,j,k,l; @@ -825,6 +917,7 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, 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; @@ -844,10 +937,10 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, basis[2].z=0.5*lc; /* shift partition values */ - if(p_type) { - p.x=p_vals->p.x+(a*lc)/2.0; - p.y=p_vals->p.y+(b*lc)/2.0; - p.z=p_vals->p.z+(c*lc)/2.0; + 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 */ @@ -858,83 +951,77 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, r.z=o.z; for(k=0;ktype) { case PART_INSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + 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: - v3_copy(&(atom[count].r),&r); - count+=1; + set_atom_positions(r); break; } /* the three face centered atoms */ for(l=0;l<3;l++) { v3_add(&n,&r,&basis[l]); - switch(p_type) { + switch(p_params->type) { case PART_INSIDE_R: v3_sub(&dist,&n,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&n); - count+=1; + 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_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&n); - count+=1; + 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)) { - v3_copy(&(atom[count].r),&n); - count+=1; + 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_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&n); - count+=1; + 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: - v3_copy(&(atom[count].r),&n); - count+=1; + set_atom_positions(n); break; } } @@ -956,12 +1043,12 @@ int fcc_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) { + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec o; - count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals); + count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params); o.x=0.25*lc; o.y=0.25*lc; @@ -969,7 +1056,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,p_type,p_vals); + count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params); return count; } @@ -1171,12 +1258,12 @@ double pressure_calc(t_moldyn *moldyn) { moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); - 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; + //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; @@ -1439,16 +1526,15 @@ int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) { int scale_volume(t_moldyn *moldyn) { t_3dvec *dim,*vdim; - //double scale; + double scale; t_linkcell *lc; - double sx,sy,sz; + //double sx,sy,sz; vdim=&(moldyn->vis.dim); dim=&(moldyn->dim); lc=&(moldyn->lc); /* scaling factor */ - /* if(moldyn->pt_scale&P_SCALE_BERENDSEN) { scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau; scale=pow(scale,ONE_THIRD); @@ -1456,20 +1542,22 @@ int scale_volume(t_moldyn *moldyn) { else { scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); } - */ + + /* 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); + 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) { @@ -1488,12 +1576,12 @@ int scale_volume(t_moldyn *moldyn) { link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); } else { - //lc->x*=scale; - //lc->y*=scale; - //lc->z*=scale; - lc->x*=sx; - lc->y*=sx; - lc->z*=sy; + lc->x*=scale; + lc->y*=scale; + lc->z*=scale; + //lc->x*=sx; + //lc->y*=sx; + //lc->z*=sy; } return 0; @@ -1507,16 +1595,16 @@ 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; + //moldyn->ekinx=0.0; + //moldyn->ekiny=0.0; + //moldyn->ekinz=0.0; 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; + //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; @@ -2002,10 +2090,11 @@ int moldyn_integrate(t_moldyn *moldyn) { if(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(!(moldyn->total_steps%m)) { @@ -3140,9 +3229,8 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, ba=data; /* increase total bond counter - * ... double counting! */ - ba->tcnt+=2; + ba->tcnt+=1; if(itom->brand==0) ba->acnt[jtom->tag]+=1; @@ -3159,10 +3247,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, int bond_analyze(t_moldyn *moldyn,double *quality) { - // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total - - int qcnt; - int ccnt,cset; + int qcnt4; + int qcnt3; + int ccnt4; + int ccnt3; + int bcnt; t_ba ba; int i; t_atom *atom; @@ -3182,9 +3271,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { memset(ba.bcnt,0,moldyn->count*sizeof(int)); ba.tcnt=0; - qcnt=0; - ccnt=0; - cset=0; + qcnt4=0; qcnt3=0; + ccnt4=0; ccnt3=0; + bcnt=0; atom=moldyn->atom; @@ -3193,27 +3282,30 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { for(i=0;icount;i++) { if(atom[i].brand==0) { if((ba.acnt[i]==0)&(ba.bcnt[i]==4)) - qcnt+=4; + qcnt4+=4; + if((ba.acnt[i]==0)&(ba.bcnt[i]==3)) + qcnt3+=3; } else { if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { - qcnt+=4; - ccnt+=1; + qcnt4+=4; + ccnt4+=1; + } + if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) { + qcnt3+=4; + ccnt3+=1; } - cset+=1; + bcnt+=1; } } - printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset); - printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt); - if(quality) { - quality[0]=1.0*ccnt/cset; - quality[1]=1.0*qcnt/ba.tcnt; + quality[0]=1.0*ccnt4/bcnt; + quality[1]=1.0*ccnt3/bcnt; } else { - printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt); - printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt); + printf("[moldyn] bond analyze: %f %f\n", + 1.0*ccnt4/bcnt,1.0*ccnt3/bcnt); } return 0; @@ -3264,6 +3356,7 @@ int visual_atoms(t_moldyn *moldyn) { t_visual *v; t_atom *atom; t_vb vb; + t_3dvec strain; #ifdef VISUAL_THREAD t_moldyn *moldyn; @@ -3294,14 +3387,18 @@ int visual_atoms(t_moldyn *moldyn) { moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help); // atomic configuration - for(i=0;icount;i++) + for(i=0;icount;i++) { + 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); + //atom[i].ekin); + sqrt(v3_absolute_square(&strain))); + } // bonds between atoms #ifndef VISUAL_THREAD