X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=eb3c4b03dfbcdb6b427c66d03aace481a6fe91ff;hb=67df6b3a722a44e36fd56f3f040c3c9726b3fc3f;hp=988b4671780d6fb68a8826c92390c2cedd2ff663;hpb=d7f67c88195ab155f2737e57cc5e81973d3feb0c;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 988b467..eb3c4b0 100644 --- a/moldyn.c +++ b/moldyn.c @@ -30,99 +30,6 @@ #include "potentials/tersoff.h" #endif - -/* - * global variables, pse and atom colors (only needed here) - */ - -static char *pse_name[]={ - "*", - "H", - "He", - "Li", - "Be", - "B", - "C", - "N", - "O", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Si", - "P", - "S", - "Cl", - "Ar", -}; - -static char *pse_col[]={ - "*", - "White", - "He", - "Li", - "Be", - "B", - "Gray", - "N", - "Blue", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Yellow", - "P", - "S", - "Cl", - "Ar", -}; - -static double *pse_mass[]={ - 0, - 0, - 0, - 0, - 0, - 0, - M_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - M_SI, - 0, - 0, - 0, - 0, -}; - -static double *pse_lc[]={ - 0, - 0, - 0, - 0, - 0, - 0, - LC_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - LC_SI, - 0, - 0, - 0, - 0, -}; - /* * the moldyn functions */ @@ -202,37 +109,27 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { - moldyn->pt_scalei&=(^(P_SCALE_MASK)); + moldyn->pt_scale&=(~(P_SCALE_MASK)); moldyn->pt_scale|=ptype; moldyn->p_tc=ptc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] p scaling:\n"); printf(" p: %s",ptype?"yes":"no "); if(ptype) printf(" | type: %02x | factor: %f",ptype,ptc); printf("\n"); - printf(" t: %s",ttype?"yes":"no "); - if(ttype) - printf(" | type: %02x | factor: %f",ttype,ttc); - printf("\n"); - return 0; } int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { - moldyn->pt_scalei&=(^(T_SCALE_MASK)); + moldyn->pt_scale&=(~(T_SCALE_MASK)); moldyn->pt_scale|=ttype; moldyn->t_tc=ttc; - printf("[moldyn] p/t scaling:\n"); - - printf(" p: %s",ptype?"yes":"no "); - if(ptype) - printf(" | type: %02x | factor: %f",ptype,ptc); - printf("\n"); + printf("[moldyn] t scaling:\n"); printf(" t: %s",ttype?"yes":"no "); if(ttype) @@ -323,7 +220,7 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func3b_k1=tersoff_mult_3bp_k1; moldyn->func3b_j2=tersoff_mult_3bp_j2; moldyn->func3b_k2=tersoff_mult_3bp_k2; - // missing: check 2b bond func + moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_AM: moldyn->func3b_j1=albe_mult_3bp_j1; @@ -582,11 +479,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, t_3dvec orig; void *ptr; t_atom *atom; + char name[16]; 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; @@ -616,18 +518,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, case CUBIC: set_nn_dist(moldyn,lc); ret=cubic_init(a,b,c,lc,atom,&orig); + 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); + 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); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -644,7 +549,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->count+=new; - printf("[moldyn] created lattice with %d atoms\n",new); + printf("[moldyn] created %s lattice with %d atoms\n",name,new); for(ret=0;ret 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) + 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; @@ -1711,7 +1619,11 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); + /* + thermodynamic_pressure_calc(moldyn); + printf("\n\nDEBUG: numeric pressure calc: %f\n\n", + moldyn->tp/BAR); + */ /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); @@ -1792,6 +1704,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", sched->count,i,moldyn->total_steps, moldyn->t,moldyn->t_avg, moldyn->p/BAR,moldyn->p_avg/BAR, + //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR, moldyn->volume, (int)(t2.tv_sec-t1.tv_sec)); @@ -2422,9 +2335,6 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, t_list *this; lc=&(moldyn->lc); - - link_cell_init(moldyn,VERBOSE); - itom=moldyn->atom; for(i=0;icount;i++) { @@ -2684,8 +2594,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) return 0; - d=sqrt(d); - /* now count this bonding ... */ ba=data;