X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=eb3c4b03dfbcdb6b427c66d03aace481a6fe91ff;hb=67df6b3a722a44e36fd56f3f040c3c9726b3fc3f;hp=84cf7003c13f70cd6e7180a7a23b53ddda6722da;hpb=2346d076f3d0955fb3a7bfb47f2b17b9eae6dd5b;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 84cf700..eb3c4b0 100644 --- a/moldyn.c +++ b/moldyn.c @@ -20,53 +20,15 @@ #include "moldyn.h" #include "report/report.h" -/* - * 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", -}; +/* potential includes */ +#include "potentials/harmonic_oscillator.h" +#include "potentials/lennard_jones.h" +#include "potentials/albe.h" +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else +#include "potentials/tersoff.h" +#endif /* * the moldyn functions @@ -120,24 +82,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) { int set_cutoff(t_moldyn *moldyn,double cutoff) { moldyn->cutoff=cutoff; + moldyn->cutoff_square=cutoff*cutoff; printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff); return 0; } -int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) { - - moldyn->bondlen[0]=b0*b0; - moldyn->bondlen[1]=b1*b1; - if(bm<0) - moldyn->bondlen[2]=b0*b1; - else - moldyn->bondlen[2]=bm*bm; - - return 0; -} - int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; @@ -156,6 +107,38 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { return 0; } +int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { + + moldyn->pt_scale&=(~(P_SCALE_MASK)); + moldyn->pt_scale|=ptype; + moldyn->p_tc=ptc; + + printf("[moldyn] p scaling:\n"); + + printf(" p: %s",ptype?"yes":"no "); + if(ptype) + printf(" | type: %02x | factor: %f",ptype,ptc); + printf("\n"); + + return 0; +} + +int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { + + moldyn->pt_scale&=(~(T_SCALE_MASK)); + moldyn->pt_scale|=ttype; + moldyn->t_tc=ttc; + + printf("[moldyn] t scaling:\n"); + + printf(" t: %s",ttype?"yes":"no "); + if(ttype) + printf(" | type: %02x | factor: %f",ttype,ttc); + printf("\n"); + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); @@ -228,58 +211,37 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func) { - - moldyn->func1b=func; - - return 0; -} - -int set_potential2b(t_moldyn *moldyn,pf_func2b func) { +int set_potential(t_moldyn *moldyn,u8 type) { - moldyn->func2b=func; - - return 0; -} - -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j1=func; - - return 0; -} - -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j2=func; - - return 0; -} - -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j3=func; - - return 0; -} - -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k1=func; - - return 0; -} - -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k2=func; - - return 0; -} - -int set_potential_params(t_moldyn *moldyn,void *params) { - - moldyn->pot_params=params; + switch(type) { + case MOLDYN_POTENTIAL_TM: + 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_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->check_2b_bond=albe_mult_check_2b_bond; + break; + case MOLDYN_POTENTIAL_HO: + moldyn->func2b=harmonic_oscillator; + moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; + break; + case MOLDYN_POTENTIAL_LJ: + moldyn->func2b=lennard_jones; + moldyn->check_2b_bond=lennard_jones_check_2b_bond; + break; + default: + printf("[moldyn] set potential: unknown type %02x\n", + type); + return -1; + } return 0; } @@ -517,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; @@ -551,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); @@ -579,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;rettau_square=moldyn->tau*moldyn->tau; - moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; /* get current time */ gettimeofday(&t1,NULL); @@ -1609,13 +1578,16 @@ int moldyn_integrate(t_moldyn *moldyn) { printf("[moldyn] WARNING: cutoff > 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; @@ -1647,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); @@ -1728,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)); @@ -2358,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++) { @@ -2424,7 +2398,8 @@ int get_line(int fd,char *line,int max) { ret=read(fd,line+count,1); if(ret<=0) return ret; if(line[count]=='\n') { - line[count]='\0'; + memset(line+count,0,max-count-1); + //line[count]='\0'; return count+1; } count+=1; @@ -2508,8 +2483,8 @@ int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, if(bc) check_per_bound(moldyn,&dist); d=v3_absolute_square(&dist); - /* ignore if greater or equal 2 times cutoff */ - if(d>=4.0*moldyn->cutoff_square) + /* ignore if greater cutoff */ + if(d>moldyn->cutoff_square) return 0; /* fill the slots */ @@ -2548,7 +2523,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { int i; pcc.dr=dr; - pcc.o1=2.0*moldyn->cutoff/dr; + pcc.o1=moldyn->cutoff/dr; pcc.o2=2*pcc.o1; if(pcc.o1*dr<=moldyn->cutoff) @@ -2612,7 +2587,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, d=v3_absolute_square(&dist); /* ignore if greater or equal cutoff */ - if(d>=moldyn->cutoff_square) + if(d>moldyn->cutoff_square) + return 0; + + /* check for potential bond */ + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) return 0; /* now count this bonding ... */ @@ -2641,6 +2620,7 @@ 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; t_ba ba; int i; t_atom *atom; @@ -2660,8 +2640,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; atom=moldyn->atom; @@ -2673,17 +2654,25 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { qcnt+=4; } else { - if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) + if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { qcnt+=4; + ccnt+=1; + } + cset+=1; } } -printf("%d %d\n",qcnt,ba.tcnt); - if(quality) - *quality=1.0*qcnt/ba.tcnt; - else - printf("[moldyn] bond analyze: quality = %f\n", - 1.0*qcnt/ba.tcnt); + 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; + } + 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); + } return 0; } @@ -2699,39 +2688,48 @@ int visual_init(t_moldyn *moldyn,char *filebase) { return 0; } +int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { + + t_vb *vb; + + vb=data; + + if(itom->tag>=jtom->tag) + return 0; + + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB)) + dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n", + itom->r.x,itom->r.y,itom->r.z, + jtom->r.x,jtom->r.y,jtom->r.z); + + return 0; +} + int visual_atoms(t_moldyn *moldyn) { - int i,j,fd; + int i; char file[128+64]; t_3dvec dim; double help; t_visual *v; t_atom *atom; - t_atom *btom; - t_linkcell *lc; -#ifdef STATIC_LISTS - int *neighbour[27]; - int p; -#else - t_list neighbour[27]; -#endif - u8 bc; - t_3dvec dist; - double d2; - u8 brand; + t_vb vb; v=&(moldyn->vis); dim.x=v->dim.x; dim.y=v->dim.y; dim.z=v->dim.z; atom=moldyn->atom; - lc=&(moldyn->lc); help=(dim.x+dim.y); sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time); - fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); - if(fd<0) { + vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(vb.fd<0) { perror("open visual save file fd"); return -1; } @@ -2739,118 +2737,65 @@ int visual_atoms(t_moldyn *moldyn) { /* write the actual data file */ // povray header - dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n", 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++) // atom type, positions, color and kinetic energy - dprintf(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); - - /* - * bond detection should usually be done by potential - * functions. brrrrr! EVIL! - * - * todo: potentials need to export a 'find_bonds' function! - */ - - // bonds between atoms - if(!(atom[i].attr&ATOM_ATTR_VB)) - continue; - link_cell_neighbour_index(moldyn, - (atom[i].r.x+moldyn->dim.x/2)/lc->x, - (atom[i].r.y+moldyn->dim.y/2)/lc->y, - (atom[i].r.z+moldyn->dim.z/2)/lc->z, - neighbour); - for(j=0;j<27;j++) { - bc=jdnlc?0:1; -#ifdef STATIC_LISTS - p=0; - while(neighbour[j][p]!=0) { - btom=&(atom[neighbour[j][p]]); - p++; -#else - list_reset_f(&neighbour[j]); - if(neighbour[j].start==NULL) - continue; - do { - btom=neighbour[j].current->data; -#endif - if(btom==&atom[i]) // skip identical atoms - continue; - //if(btom<&atom[i]) // skip half of them - // continue; - v3_sub(&dist,&(atom[i].r),&(btom->r)); - if(bc) check_per_bound(moldyn,&dist); - d2=v3_absolute_square(&dist); - brand=atom[i].brand; - if(brand==btom->brand) { - if(d2>moldyn->bondlen[brand]) - continue; - } - else { - if(d2>moldyn->bondlen[2]) - continue; - } - dprintf(fd,"# [B] %f %f %f %f %f %f\n", - atom[i].r.x,atom[i].r.y,atom[i].r.z, - btom->r.x,btom->r.y,btom->r.z); -#ifdef STATIC_LISTS - } -#else - } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT); -#endif - } - } - + 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); + + // bonds between atoms + process_2b_bonds(moldyn,&vb,visual_bonds_process); + // boundaries if(dim.x) { - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,-dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,-dim.z/2, -dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,-dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,-dim.z/2, dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, dim.x/2,-dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, -dim.x/2,dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,dim.z/2, dim.x/2,-dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,dim.z/2, dim.x/2,dim.y/2,dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,-dim.y/2,dim.z/2, -dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", -dim.x/2,dim.y/2,dim.z/2, -dim.x/2,dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,-dim.y/2,dim.z/2, dim.x/2,-dim.y/2,-dim.z/2); - dprintf(fd,"# [D] %f %f %f %f %f %f\n", + dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n", dim.x/2,dim.y/2,dim.z/2, dim.x/2,dim.y/2,-dim.z/2); } - close(fd); + close(vb.fd); return 0; }