X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1d6e0b7dec9e3032df43f33baa692ba078bb5f7a;hp=bcc78c011a41a646623a6191ca091f18899ed302;hb=HEAD;hpb=36923a44342c3b07e84baec93a30b22c792f74fc diff --git a/moldyn.c b/moldyn.c index bcc78c0..1d6e0b7 100644 --- a/moldyn.c +++ b/moldyn.c @@ -34,6 +34,7 @@ #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 @@ -41,9 +42,11 @@ #endif /* pse */ +#define PSE_MASS #define PSE_NAME #define PSE_COL #include "pse.h" +#undef PSE_MASS #undef PSE_NAME #undef PSE_COL @@ -53,6 +56,11 @@ 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 */ @@ -256,26 +264,36 @@ int set_potential(t_moldyn *moldyn,u8 type) { switch(type) { case MOLDYN_POTENTIAL_TM: - moldyn->func1b=tersoff_mult_1bp; + //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->func3b_j1=albe_mult_3bp_j1; - moldyn->func3b_k1=albe_mult_3bp_k1; - moldyn->func3b_j2=albe_mult_3bp_j2; - moldyn->func3b_k2=albe_mult_3bp_k2; + 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->func2b=harmonic_oscillator; + moldyn->func_j0=harmonic_oscillator; moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; break; case MOLDYN_POTENTIAL_LJ: - moldyn->func2b=lennard_jones; + moldyn->func_j0=lennard_jones; moldyn->check_2b_bond=lennard_jones_check_2b_bond; break; default: @@ -512,8 +530,10 @@ 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; @@ -537,6 +557,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) { @@ -570,22 +605,28 @@ 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: @@ -595,28 +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); - return -1; } - moldyn->count+=new; - printf("[moldyn] created %s lattice with %d atoms\n",name,new); + moldyn->count+=ret; + if(ret==new) + printf("[moldyn] created %s lattice with %d atoms\n",name,ret); - for(ret=0;rettype) { + 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[ret]),NULL); + 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; + +// WHAT ABOUT AMUTEX !!!! + +#ifdef LOWMEM_LISTS + ptr=realloc(moldyn->lc.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); @@ -624,7 +699,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; @@ -667,7 +742,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; @@ -684,6 +759,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; @@ -706,16 +784,81 @@ 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) { +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) @@ -723,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; @@ -748,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) @@ -770,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; } @@ -803,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; @@ -816,7 +1079,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { if(origin) v3_add(&o,&o,origin); - count+=fcc_init(a,b,c,lc,&atom[count],&o); + count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params); return count; } @@ -894,7 +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->count) + moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); + else moldyn->t=0.0; return moldyn->t; } @@ -994,8 +1259,8 @@ double virial_sum(t_moldyn *moldyn) { } /* global virial (absolute coordinates) */ - virial=&(moldyn->gvir); - moldyn->gv=virial->xx+virial->yy+virial->zz; + //virial=&(moldyn->gvir); + //moldyn->gv=virial->xx+virial->yy+virial->zz; return moldyn->virial; } @@ -1016,9 +1281,16 @@ 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; + /* pressure (absolute coordinates) */ - moldyn->gp=2.0*moldyn->ekin+moldyn->gv; - moldyn->gp/=(3.0*moldyn->volume); + //moldyn->gp=2.0*moldyn->ekin+moldyn->gv; + //moldyn->gp/=(3.0*moldyn->volume); return moldyn->p; } @@ -1043,11 +1315,11 @@ int average_reset(t_moldyn *moldyn) { /* virial */ moldyn->virial_sum=0.0; - moldyn->gv_sum=0.0; + //moldyn->gv_sum=0.0; /* pressure */ moldyn->p_sum=0.0; - moldyn->gp_sum=0.0; + //moldyn->gp_sum=0.0; moldyn->tp_sum=0.0; return 0; @@ -1085,14 +1357,14 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) { /* 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; + //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->gp_sum+=moldyn->gp; + //moldyn->gp_avg=moldyn->gp_sum/denom; moldyn->tp_sum+=moldyn->tp; moldyn->tp_avg=moldyn->tp_sum/denom; @@ -1246,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); @@ -1265,9 +1566,21 @@ int scale_volume(t_moldyn *moldyn) { 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); /* visualize dimensions */ if(vdim->x!=0) { @@ -1289,6 +1602,9 @@ int scale_volume(t_moldyn *moldyn) { lc->x*=scale; lc->y*=scale; lc->z*=scale; + //lc->x*=sx; + //lc->y*=sx; + //lc->z*=sy; } return 0; @@ -1302,10 +1618,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; 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; @@ -1740,6 +2062,11 @@ int moldyn_integrate(t_moldyn *moldyn) { } #endif +#ifdef PTHREADS + printf("##################\n"); + printf("# USING PTHREADS #\n"); + printf("##################\n"); +#endif /* tell the world */ printf("[moldyn] integration start, go get a coffee ...\n"); @@ -1767,11 +2094,11 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); - /* +#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); @@ -1786,10 +2113,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)) { @@ -1805,7 +2133,7 @@ int moldyn_integrate(t_moldyn *moldyn) { dprintf(moldyn->pfd, "%f %f %f %f %f %f %f\n",moldyn->time, moldyn->p/BAR,moldyn->p_avg/BAR, - moldyn->gp/BAR,moldyn->gp_avg/BAR, + moldyn->p/BAR,moldyn->p_avg/BAR, moldyn->tp/BAR,moldyn->tp_avg/BAR); } } @@ -1819,7 +2147,11 @@ int moldyn_integrate(t_moldyn *moldyn) { if(v) { if(!(moldyn->total_steps%v)) { dprintf(moldyn->vfd, - "%f %f\n",moldyn->time,moldyn->volume); + "%f %f %f %f %f\n",moldyn->time, + moldyn->volume, + moldyn->dim.x, + moldyn->dim.y, + moldyn->dim.z); } } if(s) { @@ -1873,15 +2205,20 @@ int moldyn_integrate(t_moldyn *moldyn) { } /* display progress */ +#ifndef PDEBUG if(!(i%10)) { +#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, - //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/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)); @@ -1889,7 +2226,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n", /* copy over time */ t1=t2; +#ifndef PDEBUG } +#endif /* increase absolute time */ moldyn->time+=moldyn->tau; @@ -1910,6 +2249,59 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n", } + /* 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; } @@ -1928,15 +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 @@ -1961,6 +2389,33 @@ int velocity_verlet(t_moldyn *moldyn) { /* check whether fixed atom */ if(atom[i].attr&ATOM_ATTR_FP) continue; + + /* constraint relaxation */ + if(crtt) { + // forces + basis_trafo(&(atom[i].f),FORWARD, + trafo_angle[2*i],trafo_angle[2*i+1]); + if(constraints[3*i]) + atom[i].f.x=0; + if(constraints[3*i+1]) + atom[i].f.y=0; + if(constraints[3*i+2]) + atom[i].f.z=0; + basis_trafo(&(atom[i].f),BACKWARD, + trafo_angle[2*i],trafo_angle[2*i+1]); + // velocities + basis_trafo(&(atom[i].v),FORWARD, + trafo_angle[2*i],trafo_angle[2*i+1]); + if(constraints[3*i]) + atom[i].v.x=0; + if(constraints[3*i+1]) + atom[i].v.y=0; + if(constraints[3*i+2]) + atom[i].v.z=0; + basis_trafo(&(atom[i].v),BACKWARD, + trafo_angle[2*i],trafo_angle[2*i+1]); + } + /* again velocities [actually v(t+tau)] */ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass); v3_add(&(atom[i].v),&(atom[i].v),&delta); @@ -2043,8 +2498,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) - if(moldyn->func1b) - 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; @@ -2059,8 +2514,14 @@ 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=(jrun3bp=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 } @@ -2110,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) { } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; @@ -2163,18 +2681,18 @@ 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++) { @@ -2207,8 +2725,8 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; @@ -2232,14 +2750,15 @@ int potential_force_calc(t_moldyn *moldyn) { } - 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++) { @@ -2272,17 +2791,17 @@ int potential_force_calc(t_moldyn *moldyn) { 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 } @@ -2297,11 +2816,11 @@ int potential_force_calc(t_moldyn *moldyn) { } - /* 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 } @@ -2312,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) { #endif } - + #ifdef DEBUG //printf("\n\n"); #endif @@ -2404,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 */ @@ -2494,7 +3058,7 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { if(fsize!=sizeof(t_moldyn)+size) { corr=fsize-sizeof(t_moldyn)-size; printf("[moldyn] WARNING: lsf (illegal file size)\n"); - printf(" moifying offset:\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); @@ -2517,6 +3081,16 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { 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; @@ -2729,6 +3303,7 @@ 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; @@ -2742,8 +3317,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) { for(i=0;icount;i++) { - v3_sub(&dist,&(atom[i].r),&(atom[i].r_0)); - check_per_bound(moldyn,&dist); + /* 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) { @@ -2765,6 +3344,57 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) { 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; @@ -2910,9 +3540,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; @@ -2929,10 +3558,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; @@ -2952,9 +3582,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; @@ -2963,27 +3593,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; } - cset+=1; + if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) { + qcnt3+=4; + ccnt3+=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; @@ -3034,6 +3667,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; @@ -3064,14 +3698,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