X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=d715c52a11f91aa9f43244b014c281f3c2a00cc4;hb=fe84cdc3fc897581636aff89b6f4fd12db1d2718;hp=7698c24767b90de44a165602077a46714d86d476;hpb=a16579f010cb5d53b7bf44e854089019ab9b1df0;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 7698c24..d715c52 100644 --- a/moldyn.c +++ b/moldyn.c @@ -55,6 +55,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 */ @@ -258,7 +263,7 @@ 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; @@ -516,7 +521,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { 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_part_params *p_params,t_defect_params *d_params, + t_offset_params *o_params) { int new,count; int ret; @@ -588,6 +594,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, 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,p_params,d_params); strcpy(name,"cubic"); @@ -595,6 +603,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, 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,p_params,d_params); strcpy(name,"fcc"); @@ -602,6 +612,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, 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,p_params,d_params); strcpy(name,"diamond"); @@ -2071,11 +2083,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); @@ -2182,15 +2194,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)); @@ -2198,7 +2215,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; @@ -2219,6 +2238,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; } @@ -2237,15 +2309,35 @@ 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 @@ -2270,6 +2366,21 @@ int velocity_verlet(t_moldyn *moldyn) { /* check whether fixed atom */ if(atom[i].attr&ATOM_ATTR_FP) continue; + + /* constraint relaxation */ + if(crtt) { + 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]); + } + /* 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); @@ -2713,6 +2824,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 */ @@ -2803,7 +2959,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); @@ -3048,6 +3204,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; @@ -3061,8 +3218,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) { @@ -3084,6 +3245,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;