From: hackbard Date: Mon, 20 Feb 2012 15:55:56 +0000 (+0100) Subject: Merge branch 'leadoff' X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9;hp=-c Merge branch 'leadoff' Conflicts: Makefile moldyn.c potentials/albe.c potentials/albe.h runmd ... hopefully fixed! Just a test to merge leadoff back to master! :) --- 6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9 diff --combined mdrun.c index 907388a,fb1cbe0..26fbc87 --- a/mdrun.c +++ b/mdrun.c @@@ -11,7 -11,6 +11,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 @@@ -94,9 -93,18 +94,18 @@@ int add_stage(t_mdrun *mdrun,u8 type,vo case STAGE_DISPLACE_ATOM: psize=sizeof(t_displace_atom_params); break; + case STAGE_DEL_ATOMS: + psize=sizeof(t_del_atoms_params); + break; + case STAGE_MODIFY_ATOMS: + psize=sizeof(t_modify_atoms_params); + break; case STAGE_INSERT_ATOMS: psize=sizeof(t_insert_atoms_params); break; + case STAGE_INSERT_MIXED_ATOMS: + psize=sizeof(t_insert_mixed_atoms_params); + break; case STAGE_CONTINUE: psize=sizeof(t_continue_params); break; @@@ -109,6 -117,21 +118,21 @@@ case STAGE_CHSATTR: psize=sizeof(t_chsattr_params); break; + case STAGE_SET_TEMP: + psize=sizeof(t_set_temp_params); + break; + case STAGE_SET_TIMESTEP: + psize=sizeof(t_set_timestep_params); + break; + case STAGE_FILL: + psize=sizeof(t_fill_params); + break; + case STAGE_THERMAL_INIT: + psize=0; + break; + case STAGE_CRT: + psize=sizeof(t_crt_params); + break; default: printf("%s unknown stage type: %02x\n",ME,type); return -1; @@@ -142,16 -165,23 +166,23 @@@ int mdrun_parse_config(t_mdrun *mdrun) char error[128]; char line[128]; char *wptr; - char word[16][64]; + char word[32][64]; int wcnt; int i,o; t_displace_atom_params dap; + t_modify_atoms_params map; t_insert_atoms_params iap; + t_insert_mixed_atoms_params imp; t_continue_params cp; t_anneal_params ap; t_chaattr_params cap; t_chsattr_params csp; + t_set_temp_params stp; + t_set_timestep_params stsp; + t_fill_params fp; + t_del_atoms_params delp; + t_crt_params crtp; /* open config file */ fd=open(mdrun->cfile,O_RDONLY); @@@ -176,10 -206,17 +207,17 @@@ // reset memset(&iap,0,sizeof(t_insert_atoms_params)); + memset(&map,0,sizeof(t_modify_atoms_params)); + memset(&imp,0,sizeof(t_insert_mixed_atoms_params)); memset(&cp,0,sizeof(t_continue_params)); memset(&ap,0,sizeof(t_anneal_params)); memset(&cap,0,sizeof(t_chaattr_params)); memset(&csp,0,sizeof(t_chsattr_params)); + memset(&stp,0,sizeof(t_set_temp_params)); + memset(&stsp,0,sizeof(t_set_timestep_params)); + memset(&fp,0,sizeof(t_fill_params)); + memset(&delp,0,sizeof(t_del_atoms_params)); + memset(&crtp,0,sizeof(t_crt_params)); // get command + args wcnt=0; @@@ -197,8 -234,6 +235,8 @@@ if(!strncmp(word[0],"potential",9)) { if(!strncmp(word[1],"albe",4)) mdrun->potential=MOLDYN_POTENTIAL_AM; + if(!strncmp(word[1],"albe_o",6)) + mdrun->potential=MOLDYN_POTENTIAL_AO; if(!strncmp(word[1],"tersoff",7)) mdrun->potential=MOLDYN_POTENTIAL_TM; if(!strncmp(word[1],"ho",2)) @@@ -250,21 -285,133 +288,133 @@@ mdrun->lattice=FCC; if(!strncmp(word[1],"diamond",7)) mdrun->lattice=DIAMOND; + if(!strncmp(word[1],"none",4)) + mdrun->lattice=NONE; + if(wcnt==3) + mdrun->lc=atof(word[2]); } else if(!strncmp(word[0],"element1",8)) { mdrun->element1=atoi(word[1]); - mdrun->m1=pse_mass[mdrun->element1]; } else if(!strncmp(word[0],"element2",8)) { mdrun->element2=atoi(word[1]); - mdrun->m2=pse_mass[mdrun->element2]; } else if(!strncmp(word[0],"fill",6)) { - // only lc mode by now - mdrun->lx=atoi(word[2]); - mdrun->ly=atoi(word[3]); - mdrun->lz=atoi(word[4]); - mdrun->lc=atof(word[5]); + // default values + fp.fill_element=mdrun->element1; + fp.fill_brand=0; + fp.lattice=mdrun->lattice; + fp.p_params.type=0; + fp.d_params.type=0; + fp.d_params.stype=0; + // parse fill command + i=1; + while(ilc=fp.lc; + continue; + } + if(!strncmp(word[i],"eb",2)) { + fp.fill_element=atoi(word[++i]); + fp.fill_brand=atoi(word[++i]); + continue; + } + if(word[i][0]=='p') { + i+=1; + switch(word[i][0]) { + case 'i': + if(word[i][1]=='r') + fp.p_params.type=PART_INSIDE_R; + else + fp.p_params.type=PART_INSIDE_D; + break; + case 'o': + if(word[i][1]=='r') + fp.p_params.type=PART_OUTSIDE_R; + else + fp.p_params.type=PART_OUTSIDE_D; + break; + default: + break; + } + if((fp.p_params.type==PART_INSIDE_R)|| + (fp.p_params.type==PART_OUTSIDE_R)) { + fp.p_params.r=atof(word[++i]); + fp.p_params.p.x=atof(word[++i]); + fp.p_params.p.y=atof(word[++i]); + fp.p_params.p.z=atof(word[++i]); + } + if((fp.p_params.type==PART_INSIDE_D)|| + (fp.p_params.type==PART_OUTSIDE_D)) { + fp.p_params.p.x=atof(word[++i]); + fp.p_params.p.y=atof(word[++i]); + fp.p_params.p.z=atof(word[++i]); + fp.p_params.d.x=atof(word[++i]); + fp.p_params.d.y=atof(word[++i]); + fp.p_params.d.z=atof(word[++i]); + } + continue; + } + if(word[i][0]=='d') { + switch(word[++i][0]) { + case '0': + + fp.d_params.type=DEFECT_TYPE_0D; + if(!strncmp(word[i+1],"dbx",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_X; + } + if(!strncmp(word[i+1],"dby",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_Y; + } + if(!strncmp(word[i+1],"dbz",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_Z; + } + if(!strncmp(word[i+1],"dbr",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_R; + } + i+=1; + fp.d_params.od=atof(word[++i]); + fp.d_params.dd=atof(word[++i]); + fp.d_params.element=atoi(word[++i]); + fp.d_params.brand=atoi(word[++i]); + // parsed in future + fp.d_params.attr=ATOM_ATTR_HB|ATOM_ATTR_VA; + fp.d_params.attr|=ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP; + break; + + case '1': + fp.d_params.type=DEFECT_TYPE_1D; + break; + case '2': + fp.d_params.type=DEFECT_TYPE_2D; + break; + case '3': + fp.d_params.type=DEFECT_TYPE_3D; + break; + default: + break; + } + continue; + + } + // offset + if(word[i][0]=='o') { + fp.o_params.o.x=atof(word[++i])*fp.lc; + fp.o_params.o.y=atof(word[++i])*fp.lc; + fp.o_params.o.z=atof(word[++i])*fp.lc; + fp.o_params.use=1; + continue; + } + i+=1; + } + add_stage(mdrun,STAGE_FILL,&fp); + } + else if(!strncmp(word[0],"thermal_init",12)) { + add_stage(mdrun,STAGE_THERMAL_INIT,NULL); } else if(!strncmp(word[0],"aattr",5)) { // for aatrib line we need a special stage @@@ -280,24 -427,30 +430,30 @@@ case 'e': cap.type|=CHAATTR_ELEMENT; break; + case 'n': + cap.type|=CHAATTR_NUMBER; default: break; } } i=2; if(cap.type&CHAATTR_REGION) { - cap.x0=atof(word[1]); - cap.y0=atof(word[2]); - cap.z0=atof(word[3]); - cap.x1=atof(word[4]); - cap.y1=atof(word[5]); - cap.z1=atof(word[6]); + cap.x0=atof(word[2]); + cap.y0=atof(word[3]); + cap.z0=atof(word[4]); + cap.x1=atof(word[5]); + cap.y1=atof(word[6]); + cap.z1=atof(word[7]); i+=6; } if(cap.type&CHAATTR_ELEMENT) { cap.element=atoi(word[i]); i+=1; } + if(cap.type&CHAATTR_NUMBER) { + cap.element=atoi(word[i]); + i+=1; + } for(o=0;o0) + csp.ptau=0.01/(csp.ptau*GPA); csp.type|=CHSATTR_PCTRL; } if(!strncmp(word[i],"tctrl",5)) { @@@ -384,6 -539,29 +542,29 @@@ dap.dz=atof(word[5]); add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap); } + else if(!strncmp(word[1],"del_atoms",9)) { + delp.o.x=atof(word[2]); + delp.o.y=atof(word[3]); + delp.o.z=atof(word[4]); + delp.r=atof(word[5]); + add_stage(mdrun,STAGE_DEL_ATOMS,&delp); + } + else if(!strncmp(word[1],"mod_atoms",8)) { + i=2; + while(istage.current->data; + delp=stage->params; + + if(delp->r<0) + outer=1; + + for(i=0;icount;i++) { + v3_sub(&dist,&(delp->o),&(moldyn->atom[i].r)); + //printf("%d ----> %f %f %f = %f | %f\n",i,dist.x,dist.y,dist.z,v3_absolute_square(&dist),delp->r*delp->r); + if(v3_absolute_square(&dist)<=(delp->r*delp->r)) { + if(!outer) { + del_atom(moldyn,moldyn->atom[i].tag); + printf("%s atom deleted: %d %d %d\n",ME, + moldyn->atom[i].tag, + moldyn->atom[i].element, + moldyn->atom[i].brand); + } + } + else { + if(outer) { + del_atom(moldyn,moldyn->atom[i].tag); + printf("%s atom deleted: %d %d %d\n",ME, + moldyn->atom[i].tag, + moldyn->atom[i].element, + moldyn->atom[i].brand); + } + } + } + + return 0; + } + + int modify_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_modify_atoms_params *map; + t_stage *stage; + t_atom *atom; + t_3dvec v; + int i; + + atom=moldyn->atom; + stage=mdrun->stage.current->data; + map=stage->params; + v.x=0.0; v.y=0.0; v.z=0.0; + + for(i=0;icount;i++) { + if(atom[i].tag==map->tag) { + v.x=sqrt(2.0*fabs(map->ekin.x)/atom[i].mass); + if(map->ekin.x<0.0) + v.x=-v.x; + v.y=sqrt(2.0*fabs(map->ekin.y)/atom[i].mass); + if(map->ekin.y<0.0) + v.y=-v.y; + v.z=sqrt(2.0*fabs(map->ekin.z)/atom[i].mass); + if(map->ekin.z<0.0) + v.z=-v.z; + v3_copy(&(atom[i].v),&v); + printf("%s atom modified: v = (%f %f %f)\n", + ME,v.x,v.y,v.z); + } + } + + return 0; + } + int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { t_insert_atoms_params *iap; @@@ -573,7 -917,7 +920,7 @@@ z0=-z/2.0; cr_check=TRUE; break; - case INS_REGION: + case INS_RECT: x=iap->x1-iap->x0; x0=iap->x0; y=iap->y1-iap->y0; @@@ -582,7 -926,17 +929,17 @@@ z0=iap->z0; cr_check=TRUE; break; + case INS_SPHERE: + x=2.0*iap->x1; + x0=iap->x0-iap->x1; + y=x; + y0=iap->y0-iap->x1; + z=x; + z0=iap->z0-iap->x1; + cr_check=TRUE; + break; case INS_POS: + case INS_RELPOS: x0=iap->x0; y0=iap->y0; z0=iap->z0; @@@ -598,7 -952,7 +955,7 @@@ while(cntins_atoms) { run=1; while(run) { - if(iap->type!=INS_POS) { + if((iap->type!=INS_POS)&&(iap->type!=INS_RELPOS)) { r.x=rand_get_double(&(moldyn->random))*x; r.y=rand_get_double(&(moldyn->random))*y; r.z=rand_get_double(&(moldyn->random))*z; @@@ -608,9 -962,16 +965,16 @@@ r.y=0.0; r.z=0.0; } - r.x+=x0; - r.y+=y0; - r.z+=z0; + if(iap->type==INS_RELPOS) { + r.x+=x0*mdrun->lc; + r.y+=y0*mdrun->lc; + r.z+=z0*mdrun->lc; + } + else { + r.x+=x0; + r.y+=y0; + r.z+=z0; + } // offset if(iap->type!=INS_TOTAL) { r.x+=o; @@@ -633,19 -994,146 +997,146 @@@ dmin=d; } } + if(iap->type==INS_SPHERE) { + if((r.x-iap->x0)*(r.x-iap->x0)+ + (r.y-iap->y0)*(r.y-iap->y0)+ + (r.z-iap->z0)*(r.z-iap->z0)> + (iap->x1*iap->x1)) { + run=1; + } + } } - add_atom(moldyn,iap->element,pse_mass[iap->element], + add_atom(moldyn,iap->element, iap->brand,iap->attr,&r,&v); printf("%s atom inserted (%d/%d): %f %f %f\n", ME,(iap->cnt_steps+1)*iap->ins_atoms, iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z); - printf(" -> d2 = %f/%f\n",dmin,iap->cr*iap->cr); + printf(" attributes: "); + if(iap->attr&ATOM_ATTR_VB) + printf("b "); + if(iap->attr&ATOM_ATTR_HB) + printf("h "); + if(iap->attr&ATOM_ATTR_VA) + printf("v "); + if(iap->attr&ATOM_ATTR_FP) + printf("f "); + if(iap->attr&ATOM_ATTR_1BP) + printf("1 "); + if(iap->attr&ATOM_ATTR_2BP) + printf("2 "); + if(iap->attr&ATOM_ATTR_3BP) + printf("3 "); + printf("\n"); + printf(" d2 = %f/%f\n",dmin,iap->cr*iap->cr); cnt+=1; } return 0; } + int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_stage *stage; + t_insert_mixed_atoms_params *imp; + t_atom *atom; + double x,x0,y,y0,z,z0; + double dmin,d,cmin,cmax; + u8 retry; + t_3dvec r,v,dist; + int i; + + + stage=mdrun->stage.current->data; + imp=stage->params; + + x=moldyn->dim.x; + x0=-x/2.0; + y=moldyn->dim.y; + y0=-y/2.0; + z=moldyn->dim.z; + z0=-z/2.0; + + v.x=0.0; + v.y=0.0; + v.z=0.0; + + cmin=imp->crmin*imp->crmin; + cmax=imp->crmax*imp->crmax; + + while(imp->amount1|imp->amount2) { + if(imp->amount1) { + retry=1; + while(retry) { + retry=0; + r.x=rand_get_double(&(moldyn->random))*x; + r.y=rand_get_double(&(moldyn->random))*y; + r.z=rand_get_double(&(moldyn->random))*z; + r.x+=x0; + r.y+=y0; + r.z+=z0; + dmin=1000.0; // for sure too high! + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); + v3_sub(&dist,&(atom->r),&r); + check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + if(dcmax) + retry=1; + } + add_atom(moldyn,imp->element1, + imp->brand1,imp->attr1,&r,&v); + printf("%s (mixed) atom inserted (%d): %f %f %f\n", + ME,imp->amount1,r.x,r.y,r.z); + printf(" -> d2 = %f/%f/%f\n",dmin,cmin,cmax); + imp->amount1-=1; + } + if(imp->amount2) { + retry=1; + while(retry) { + retry=0; + r.x=rand_get_double(&(moldyn->random))*x; + r.y=rand_get_double(&(moldyn->random))*y; + r.z=rand_get_double(&(moldyn->random))*z; + r.x+=x0; + r.y+=y0; + r.z+=z0; + dmin=1000.0; // for sure too high! + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); + v3_sub(&dist,&(atom->r),&r); + check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + if(dcmax) + retry=1; + } + add_atom(moldyn,imp->element2, + imp->brand2,imp->attr2,&r,&v); + printf("%s (mixed) atom inserted (%d): %f %f %f\n", + ME,imp->amount2,r.x,r.y,r.z); + printf(" -> d2 = %f/%f/%f\n",dmin,cmin,cmax); + imp->amount2-=1; + } + } + + return 0; + } + int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) { t_stage *stage; @@@ -662,23 -1150,32 +1153,32 @@@ if(cap->element!=atom->element) continue; } + if(cap->type&CHAATTR_NUMBER) { + if(cap->element!=atom->tag) + continue; + } if(cap->type&CHAATTR_REGION) { - if(cap->x0r.x) + if(cap->x0>atom->r.x) continue; - if(cap->y0r.y) + if(cap->y0>atom->r.y) continue; - if(cap->z0r.z) + if(cap->z0>atom->r.z) continue; - if(cap->x1>atom->r.x) + if(cap->x1r.x) continue; - if(cap->y1>atom->r.y) + if(cap->y1r.y) continue; - if(cap->z1>atom->r.z) + if(cap->z1r.z) continue; } + if(!(cap->type&CHAATTR_TOTALV)) + printf(" changing attributes of atom %d (0x%x)\n", + i,cap->attr); atom->attr=cap->attr; } + printf("\n\n"); + return 0; } @@@ -694,13 -1191,13 +1194,13 @@@ int chsattr(t_moldyn *moldyn,t_mdrun *m if(csp->ptau>0) set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau); else - set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau); + set_p_scale(moldyn,P_SCALE_NONE,1.0); } if(csp->type&CHSATTR_TCTRL) { if(csp->ttau>0) set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau); else - set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau); + set_t_scale(moldyn,T_SCALE_NONE,1.0); } if(csp->type&CHSATTR_PRELAX) { if(csp->dp<0) @@@ -729,6 -1226,138 +1229,138 @@@ return 0; } + int crt(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_stage *stage; + t_crt_params *crtp; + + int fd; + char line[128]; + char *wptr; + int acount; + int ret; + void *ptr; + + t_atom *atom; + t_3dvec disp; + double frac; + int i; + + stage=mdrun->stage.current->data; + crtp=stage->params; + + acount=0; + + /* initial stuff */ + + if(crtp->count==0) { + printf(" crt init\n"); + // read final positions, constraints and do the alloc + fd=open(crtp->file,O_RDONLY); + if(fd<0) { + perror("[mdrun] FATAL reading constraints file"); + return fd; + } + while(1) { + ret=get_line(fd,line,128); + // check for end of file + if(ret<=0) { + printf(" read %d atom positions\n",acount); + if(acount!=moldyn->count) + printf(" atom count mismatch!!!\n"); + printf("\n"); + break; + } + // ignore # lines and \n + if((line[0]=='#')|(ret==1)) + continue; + // allocate new memory + ptr=realloc(crtp->r_fin,(acount+1)*sizeof(t_3dvec)); + if(ptr==NULL) { + perror("[mdrun] FATAL realloc crt positions"); + return -1; + } + crtp->r_fin=ptr; + ptr=realloc(constraints,(acount+1)*3*sizeof(u8)); + if(ptr==NULL) { + perror("[mdrun] FATAL realloc crt constraints"); + return -1; + } + constraints=ptr; + // ignore type + wptr=strtok(line," \t"); + // read x y z + wptr=strtok(NULL," \t"); + crtp->r_fin[acount].x=atof(wptr); + wptr=strtok(NULL," \t"); + crtp->r_fin[acount].y=atof(wptr); + wptr=strtok(NULL," \t"); + crtp->r_fin[acount].z=atof(wptr); + // read constraints + wptr=strtok(NULL," \t"); + constraints[3*acount]=atoi(wptr); + wptr=strtok(NULL," \t"); + constraints[3*acount+1]=atoi(wptr); + wptr=strtok(NULL," \t"); + constraints[3*acount+2]=atoi(wptr); + // done reading + acount+=1; + } + close(fd); + // allocate trafo angles + trafo_angle=malloc(acount*2*sizeof(double)); + if(trafo_angle==NULL) { + perror("[mdrun] FATAL alloc trafo angles"); + return -1; + } + // set crt mode + crtt=crtp->type; + } + + /* write a save file s-crt_xofy.save */ + snprintf(line,128,"%s/s-crt_%03dof%03d.save", + moldyn->vlsdir,crtp->count,crtp->steps); + fd=open(line,O_WRONLY|O_TRUNC|O_CREAT,S_IRUSR|S_IWUSR); + if(fd<0) perror("[mdrun] crt save fd open"); + else { + write(fd,moldyn,sizeof(t_moldyn)); + write(fd,moldyn->atom, + moldyn->count*sizeof(t_atom)); + } + close(fd); + /* visualize atoms */ + visual_atoms(moldyn); + + /* output energy */ + printf(" crt energy: %d - %f\n\n", + crtp->count,(moldyn->ekin+moldyn->energy)/EV); + + /* crt routines: calculate displacement + set individual constraints */ + + printf(" crt step %d of %d in total\n\n",crtp->count+1,crtp->steps); + + if((crtp->type==1)|(crtp->count==0)) + printf(" crt angle update\n\n"); + + for(i=0;icount;i++) { + // calc displacements + atom=moldyn->atom; + v3_sub(&disp,&(crtp->r_fin[i]),&(atom[i].r)); + // angles + if((crtp->type==1)|(crtp->count==0)) { + trafo_angle[2*i]=atan2(disp.x,disp.y); + trafo_angle[2*i+1]=-atan2(disp.z, + sqrt(disp.x*disp.x+disp.y*disp.y)); + } + // move atoms + frac=1.0/(crtp->steps-crtp->count); + v3_scale(&disp,&disp,frac); + v3_add(&(atom[i].r),&(atom[i].r),&disp); + } + + return 0; + } + #define stage_print(m) if(!(stage->executed)) \ printf("%s",m) @@@ -739,12 -1368,17 +1371,17 @@@ int mdrun_hook(void *ptr1,void *ptr2) t_stage *stage; t_list *sl; int steps; - double tau; u8 change_stage; + t_3dvec o; t_insert_atoms_params *iap; + t_insert_mixed_atoms_params *imp; t_continue_params *cp; t_anneal_params *ap; + t_set_temp_params *stp; + t_set_timestep_params *stsp; + t_fill_params *fp; + t_crt_params *crtp; moldyn=ptr1; mdrun=ptr2; @@@ -760,9 -1394,8 +1397,8 @@@ /* get stage description */ stage=sl->current->data; - /* default steps and tau values */ + /* steps in next schedule */ steps=mdrun->relax_steps; - tau=mdrun->timestep; /* check whether relaxation steps are necessary */ if((check_pressure(moldyn,mdrun)==TRUE)&\ @@@ -780,6 -1413,16 +1416,16 @@@ displace_atom(moldyn,mdrun); change_stage=TRUE; break; + case STAGE_DEL_ATOMS: + stage_print(" -> del atoms\n\n"); + del_atoms(moldyn,mdrun); + change_stage=TRUE; + break; + case STAGE_MODIFY_ATOMS: + stage_print(" -> modify atoms\n\n"); + modify_atoms(moldyn,mdrun); + change_stage=TRUE; + break; case STAGE_INSERT_ATOMS: stage_print(" -> insert atoms\n\n"); iap=stage->params; @@@ -790,6 -1433,18 +1436,18 @@@ insert_atoms(moldyn,mdrun); iap->cnt_steps+=1; break; + + + + case STAGE_INSERT_MIXED_ATOMS: + stage_print(" -> insert mixed atoms\n\n"); + imp=stage->params; + insert_mixed_atoms(moldyn,mdrun); + change_stage=TRUE; + break; + + + case STAGE_CONTINUE: stage_print(" -> continue\n\n"); if(stage->executed==TRUE) { @@@ -810,17 -1465,103 +1468,103 @@@ set_temperature(moldyn, moldyn->t_ref+ap->dt); ap->count+=1; + steps=ap->interval; break; case STAGE_CHAATTR: - stage_print(" -> chaattr\n\n"); + stage_print(" -> change atom attributes\n\n"); chaatr(moldyn,mdrun); change_stage=TRUE; break; case STAGE_CHSATTR: - stage_print(" -> chsattr\n\n"); + stage_print(" -> change sys attributes\n\n"); chsattr(moldyn,mdrun); change_stage=TRUE; break; + case STAGE_SET_TEMP: + stage_print(" -> set temperature\n\n"); + stp=stage->params; + if(stp->type&SET_TEMP_CURRENT) { + set_temperature(moldyn,moldyn->t_avg); + } + else { + set_temperature(moldyn,stp->val); + } + change_stage=TRUE; + break; + case STAGE_SET_TIMESTEP: + stage_print(" -> set timestep\n\n"); + stsp=stage->params; + mdrun->timestep=stsp->tau; + change_stage=TRUE; + break; + case STAGE_FILL: + stage_print(" -> fill lattice\n\n"); + fp=stage->params; + switch(fp->lattice) { + case ZINCBLENDE: + + o.x=0.5*0.25*fp->lc; + o.y=o.x; + o.z=o.x; + create_lattice(moldyn, + FCC,fp->lc, + mdrun->element1, + DEFAULT_ATOM_ATTR,0, + fp->lx,fp->ly,fp->lz, + &o, + &(fp->p_params), + &(fp->d_params), + &(fp->o_params)); + o.x+=0.25*fp->lc; + o.y=o.x; + o.z=o.x; + create_lattice(moldyn, + FCC,fp->lc, + mdrun->element2, + DEFAULT_ATOM_ATTR,1, + fp->lx,fp->ly,fp->lz, + &o, + &(fp->p_params), + &(fp->d_params), + &(fp->o_params)); + break; + + default: + + create_lattice(moldyn, + fp->lattice,fp->lc, + fp->fill_element, + DEFAULT_ATOM_ATTR, + fp->fill_brand, + fp->lx,fp->ly,fp->lz, + NULL, + &(fp->p_params), + &(fp->d_params), + &(fp->o_params)); + break; + } + moldyn_bc_check(moldyn); + change_stage=TRUE; + break; + case STAGE_THERMAL_INIT: + stage_print(" -> thermal init\n\n"); + thermal_init(moldyn,TRUE); + change_stage=TRUE; + break; + case STAGE_CRT: + stage_print(" -> constraint relaxation"); + stage_print(" technique\n\n"); + crtp=stage->params; + if(crtp->count==crtp->steps) { + free(constraints); + free(trafo_angle); + free(crtp->r_fin); + change_stage=TRUE; + break; + } + crt(moldyn,mdrun); + crtp->count+=1; + break; default: printf("%s unknwon stage type\n",ME); break; @@@ -837,7 -1578,6 +1581,6 @@@ return 0; } steps=0; - tau=mdrun->timestep; } } @@@ -850,7 -1590,7 +1593,7 @@@ } /* continue simulation */ - moldyn_add_schedule(moldyn,steps,tau); + moldyn_add_schedule(moldyn,steps,mdrun->timestep); return 0; } @@@ -859,12 -1599,17 +1602,17 @@@ int main(int argc,char **argv) t_mdrun mdrun; t_moldyn moldyn; - t_3dvec o; + //t_3dvec o; /* clear structs */ memset(&mdrun,0,sizeof(t_mdrun)); memset(&moldyn,0,sizeof(t_moldyn)); + /* init crt variables */ + crtt=0; + constraints=NULL; + trafo_angle=NULL; + /* parse arguments */ if(mdrun_parse_argv(&mdrun,argc,argv)<0) return -1; @@@ -908,11 -1653,6 +1656,11 @@@ mdrun.element1, mdrun.element2); break; + case MOLDYN_POTENTIAL_AO: + albe_orig_mult_set_params(&moldyn, + mdrun.element1, + mdrun.element2); + break; case MOLDYN_POTENTIAL_TM: tersoff_mult_set_params(&moldyn, mdrun.element1, @@@ -939,26 -1679,33 +1687,33 @@@ /* initial lattice and dimensions */ set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis); set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz); + /* replaced by fill stage !! switch(mdrun.lattice) { case FCC: - create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1, - mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,NULL); + create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element, + mdrun.m1,DEFAULT_ATOM_ATTR, + mdrun.fill_brand,mdrun.lx, + mdrun.ly,mdrun.lz,NULL,0,NULL); break; case DIAMOND: - create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1, - mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,NULL); + create_lattice(&moldyn,DIAMOND,mdrun.lc, + mdrun.fill_element, + mdrun.m1,DEFAULT_ATOM_ATTR, + mdrun.fill_brand,mdrun.lx, + mdrun.ly,mdrun.lz,NULL,0,NULL); break; case ZINCBLENDE: o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x; create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1, mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,&o); + mdrun.ly,mdrun.lz,&o,0,NULL); o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x; create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2, mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx, - mdrun.ly,mdrun.lz,&o); + mdrun.ly,mdrun.lz,&o,0,NULL); + break; + case NONE: + set_nn_dist(&moldyn,mdrun.nnd); break; default: printf("%s unknown lattice type: %02x\n", @@@ -966,11 -1713,14 +1721,14 @@@ return -1; } moldyn_bc_check(&moldyn); + replaced by fill stage !! */ /* temperature and pressure */ set_temperature(&moldyn,mdrun.temperature); set_pressure(&moldyn,mdrun.pressure); + /* replaced thermal_init stage thermal_init(&moldyn,TRUE); + */ addsched: /* first schedule */ diff --combined moldyn.c index 010b8ed,33224a6..1d6e0b7 --- a/moldyn.c +++ b/moldyn.c @@@ -17,6 -17,16 +17,16 @@@ #include #include + #include + + #ifdef PARALLEL + #include + #endif + + #if defined PTHREADS || defined VISUAL_THREAD + #include + #endif + #include "moldyn.h" #include "report/report.h" @@@ -24,13 -34,32 +34,33 @@@ #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 #include "potentials/tersoff.h" #endif + /* pse */ + #define PSE_MASS + #define PSE_NAME + #define PSE_COL + #include "pse.h" + #undef PSE_MASS + #undef PSE_NAME + #undef PSE_COL + + #ifdef PTHREADS + /* global mutexes */ + 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 */ @@@ -39,6 -68,9 +69,9 @@@ int moldyn_init(t_moldyn *moldyn,int ar printf("[moldyn] init\n"); + /* only needed if compiled without -msse2 (float-store prob!) */ + //fpu_set_rtd(); + memset(moldyn,0,sizeof(t_moldyn)); moldyn->argc=argc; @@@ -47,13 -79,28 +80,28 @@@ rand_init(&(moldyn->random),NULL,1); moldyn->random.status|=RAND_STAT_VERBOSE; + #ifdef PTHREADS + pthread_mutex_init(&emutex,NULL); + #endif + return 0; } int moldyn_shutdown(t_moldyn *moldyn) { + #ifdef PTHREADS + int i; + #endif + printf("[moldyn] shutdown\n"); + #ifdef PTHREADS + for(i=0;icount;i++) + pthread_mutex_destroy(&(amutex[i])); + free(amutex); + pthread_mutex_destroy(&emutex); + #endif + moldyn_log_shutdown(moldyn); link_cell_shutdown(moldyn); rand_close(&(moldyn->random)); @@@ -216,36 -263,26 +264,36 @@@ int set_potential(t_moldyn *moldyn,u8 t switch(type) { case MOLDYN_POTENTIAL_TM: - moldyn->func_i0=tersoff_mult_1bp; - moldyn->func_j1=tersoff_mult_3bp_j1; - moldyn->func_j1_k0=tersoff_mult_3bp_k1; - moldyn->func_j1c=tersoff_mult_3bp_j2; - moldyn->func_j1_k1=tersoff_mult_3bp_k2; + //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: @@@ -482,8 -519,10 +530,10 @@@ int moldyn_log_shutdown(t_moldyn *moldy * 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; @@@ -491,15 -530,37 +541,37 @@@ void *ptr; t_atom *atom; char name[16]; + #ifdef PTHREADS + pthread_mutex_t *mutex; + #endif 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; + /* 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) { @@@ -509,6 -570,16 +581,16 @@@ moldyn->atom=ptr; atom=&(moldyn->atom[count]); + #ifdef PTHREADS + ptr=realloc(amutex,(count+new)*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + mutex=&(amutex[count]); + #endif + /* no atoms on the boundaries (only reason: it looks better!) */ if(!origin) { orig.x=0.5*lc; @@@ -523,22 -594,28 +605,28 @@@ 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: @@@ -548,25 -625,62 +636,62 @@@ /* 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); + } + + moldyn->count+=ret; + if(ret==new) + printf("[moldyn] created %s lattice with %d atoms\n",name,ret); + + 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 */ + ptr=realloc(moldyn->atom,moldyn->count*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (create lattice - alloc fix)"); return -1; } + moldyn->atom=ptr; - moldyn->count+=new; - printf("[moldyn] created %s lattice with %d atoms\n",name,new); + // WHAT ABOUT AMUTEX !!!! - for(ret=0;retlc.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); @@@ -574,7 -688,7 +699,7 @@@ 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; @@@ -591,6 -705,25 +716,25 @@@ } moldyn->atom=ptr; + #ifdef LOWMEM_LISTS + ptr=realloc(moldyn->lc.subcell->list,(count+1)*sizeof(int)); + if(!ptr) { + perror("[moldyn] list realloc (add atom)"); + return -1; + } + moldyn->lc.subcell->list=ptr; + #endif + + #ifdef PTHREADS + ptr=realloc(amutex,(count+1)*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + pthread_mutex_init(&(amutex[count]),NULL); + #endif + atom=moldyn->atom; /* initialize new atom */ @@@ -598,7 -731,7 +742,7 @@@ 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; @@@ -615,6 -748,9 +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; @@@ -637,16 -773,81 +784,81 @@@ 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) @@@ -654,14 -855,54 +866,54 @@@ 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; @@@ -679,12 -920,18 +931,18 @@@ } /* 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) @@@ -701,6 -948,13 +959,13 @@@ 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; } @@@ -734,12 -1054,13 +1065,13 @@@ 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; @@@ -747,7 -1068,7 +1079,7 @@@ 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; } @@@ -825,7 -1146,9 +1157,9 @@@ double temperature_calc(t_moldyn *moldy /* 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; } @@@ -879,7 -1202,7 +1213,7 @@@ int scale_velocity(t_moldyn *moldyn,u8 scale*=2.0; else if(moldyn->pt_scale&T_SCALE_BERENDSEN) - scale=1.0+(scale-1.0)/moldyn->t_tc; + scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc; scale=sqrt(scale); /* velocity scaling */ @@@ -925,8 -1248,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; } @@@ -947,9 -1270,16 +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; } @@@ -974,11 -1304,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; @@@ -1016,14 -1346,14 +1357,14 @@@ int average_and_fluctuation_calc(t_mold /* 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; @@@ -1177,11 -1507,40 +1518,40 @@@ int scale_atoms(t_moldyn *moldyn,u8 dir 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); @@@ -1189,16 -1548,28 +1559,28 @@@ /* scaling factor */ if(moldyn->pt_scale&P_SCALE_BERENDSEN) { - scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc; + scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau; scale=pow(scale,ONE_THIRD); } 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); /* visualize dimensions */ if(vdim->x!=0) { @@@ -1220,6 -1591,9 +1602,9 @@@ lc->x*=scale; lc->y*=scale; lc->z*=scale; + //lc->x*=sx; + //lc->y*=sx; + //lc->z*=sy; } return 0; @@@ -1233,10 -1607,16 +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; @@@ -1284,7 -1664,9 +1675,9 @@@ double estimate_time_step(t_moldyn *mol int link_cell_init(t_moldyn *moldyn,u8 vol) { t_linkcell *lc; + #ifndef LOWMEM_LISTS int i; + #endif lc=&(moldyn->lc); @@@ -1299,6 -1681,8 +1692,8 @@@ #ifdef STATIC_LISTS lc->subcell=malloc(lc->cells*sizeof(int*)); + #elif LOWMEM_LISTS + lc->subcell=malloc(sizeof(t_lowmem_list)); #else lc->subcell=malloc(lc->cells*sizeof(t_list)); #endif @@@ -1309,12 -1693,16 +1704,16 @@@ } if(lc->cells<27) - printf("[moldyn] FATAL: less then 27 subcells!\n"); + printf("[moldyn] FATAL: less then 27 subcells! (%d)\n", + lc->cells); if(vol) { #ifdef STATIC_LISTS printf("[moldyn] initializing 'static' linked cells (%d)\n", lc->cells); + #elif LOWMEM_LISTS + printf("[moldyn] initializing 'lowmem' linked cells (%d)\n", + lc->cells); #else printf("[moldyn] initializing 'dynamic' linked cells (%d)\n", lc->cells); @@@ -1338,6 -1726,17 +1737,17 @@@ i,lc->subcell[0],lc->subcell); */ } + #elif LOWMEM_LISTS + lc->subcell->head=malloc(lc->cells*sizeof(int)); + if(lc->subcell->head==NULL) { + perror("[moldyn] head init (malloc)"); + return -1; + } + lc->subcell->list=malloc(moldyn->count*sizeof(int)); + if(lc->subcell->list==NULL) { + perror("[moldyn] list init (malloc)"); + return -1; + } #else for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); @@@ -1352,22 -1751,26 +1762,26 @@@ int link_cell_update(t_moldyn *moldyn) { int count,i,j,k; - int nx,ny; + int nx,nxy; t_atom *atom; t_linkcell *lc; #ifdef STATIC_LISTS int p; + #elif LOWMEM_LISTS + int p; #endif atom=moldyn->atom; lc=&(moldyn->lc); nx=lc->nx; - ny=lc->ny; + nxy=nx*lc->ny; for(i=0;icells;i++) #ifdef STATIC_LISTS - memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int)); + memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int)); + #elif LOWMEM_LISTS + lc->subcell->head[i]=-1; #else list_destroy_f(&(lc->subcell[i])); #endif @@@ -1379,7 -1782,7 +1793,7 @@@ #ifdef STATIC_LISTS p=0; - while(lc->subcell[i+j*nx+k*nx*ny][p]!=0) + while(lc->subcell[i+j*nx+k*nxy][p]!=-1) p++; if(p>=MAX_ATOMS_PER_LIST) { @@@ -1387,9 -1790,13 +1801,13 @@@ return -1; } - lc->subcell[i+j*nx+k*nx*ny][p]=count; + lc->subcell[i+j*nx+k*nxy][p]=count; + #elif LOWMEM_LISTS + p=i+j*nx+k*nxy; + lc->subcell->list[count]=lc->subcell->head[p]; + lc->subcell->head[p]=count; #else - list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]), + list_add_immediate_f(&(lc->subcell[i+j*nx+k*nxy]), &(atom[count])); /* if(j==0&&k==0) @@@ -1405,6 -1812,8 +1823,8 @@@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, #ifdef STATIC_LISTS int **cell + #elif LOWMEM_LISTS + int *cell #else t_list *cell #endif @@@ -1430,7 -1839,11 +1850,11 @@@ printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n", i,nx,j,ny,k,nz); + #ifndef LOWMEM_LISTS cell[0]=lc->subcell[i+j*nx+k*a]; + #else + cell[0]=lc->subcell->head[i+j*nx+k*a]; + #endif for(ci=-1;ci<=1;ci++) { bx=0; x=i+ci; @@@ -1454,10 -1867,19 +1878,19 @@@ } if(!(ci|cj|ck)) continue; if(bx|by|bz) { + #ifndef LOWMEM_LISTS cell[--count2]=lc->subcell[x+y*nx+z*a]; + #else + cell[--count2]=lc->subcell->head[x+y*nx+z*a]; + #endif + } else { + #ifndef LOWMEM_LISTS cell[count1++]=lc->subcell[x+y*nx+z*a]; + #else + cell[count1++]=lc->subcell->head[x+y*nx+z*a]; + #endif } } } @@@ -1470,11 -1892,19 +1903,19 @@@ int link_cell_shutdown(t_moldyn *moldyn) { + #ifndef LOWMEM_LISTS int i; + #endif t_linkcell *lc; lc=&(moldyn->lc); + #if LOWMEM_LISTS + free(lc->subcell->head); + free(lc->subcell->list); + + #else + for(i=0;icells;i++) { #ifdef STATIC_LISTS free(lc->subcell[i]); @@@ -1483,6 -1913,7 +1924,7 @@@ list_destroy_f(&(lc->subcell[i])); #endif } + #endif free(lc->subcell); @@@ -1551,6 -1982,17 +1993,17 @@@ int moldyn_integrate(t_moldyn *moldyn) struct timeval t1,t2; //double tp; + #ifdef VISUAL_THREAD + u8 first,change; + pthread_t io_thread; + int ret; + t_moldyn md_copy; + t_atom *atom_copy; + + first=1; + change=0; + #endif + sched=&(moldyn->schedule); atom=moldyn->atom; @@@ -1585,9 -2027,11 +2038,11 @@@ 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 */ // should have right values! @@@ -1597,6 -2041,21 +2052,21 @@@ /* debugging, ignore */ moldyn->debug=0; + /* zero & init moldyn copy */ + #ifdef VISUAL_THREAD + memset(&md_copy,0,sizeof(t_moldyn)); + atom_copy=malloc(moldyn->count*sizeof(t_atom)); + if(atom_copy==NULL) { + perror("[moldyn] malloc atom copy (init)"); + return -1; + } + #endif + + #ifdef PTHREADS + printf("##################\n"); + printf("# USING PTHREADS #\n"); + printf("##################\n"); + #endif /* tell the world */ printf("[moldyn] integration start, go get a coffee ...\n"); @@@ -1624,11 -2083,11 +2094,11 @@@ 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); @@@ -1643,10 -2102,11 +2113,11 @@@ 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)) { @@@ -1662,7 -2122,7 +2133,7 @@@ 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); } } @@@ -1676,12 -2136,16 +2147,16 @@@ 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) { if(!(moldyn->total_steps%s)) { - snprintf(dir,128,"%s/s-%07.f.save", + snprintf(dir,128,"%s/s-%08.f.save", moldyn->vlsdir,moldyn->time); fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT, S_IRUSR|S_IWUSR); @@@ -1696,20 -2160,54 +2171,54 @@@ } if(a) { if(!(moldyn->total_steps%a)) { + #ifdef VISUAL_THREAD + /* check whether thread has not terminated yet */ + if(!first) { + ret=pthread_join(io_thread,NULL); + } + first=0; + /* prepare and start thread */ + if(moldyn->count!=md_copy.count) { + free(atom_copy); + change=1; + } + memcpy(&md_copy,moldyn,sizeof(t_moldyn)); + if(change) { + atom_copy=malloc(moldyn->count*sizeof(t_atom)); + if(atom_copy==NULL) { + perror("[moldyn] malloc atom copy (change)"); + return -1; + } + } + md_copy.atom=atom_copy; + memcpy(atom_copy,moldyn->atom,moldyn->count*sizeof(t_atom)); + change=0; + ret=pthread_create(&io_thread,NULL,visual_atoms,&md_copy); + if(ret) { + perror("[moldyn] create visual atoms thread\n"); + return -1; + } + #else visual_atoms(moldyn); + #endif } } /* display progress */ - //if(!(moldyn->total_steps%10)) { + #ifndef PDEBUG + if(!(i%10)) { + #endif /* get current time */ gettimeofday(&t2,NULL); - printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", + 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)); @@@ -1717,7 -2215,9 +2226,9 @@@ /* copy over time */ t1=t2; - //} + #ifndef PDEBUG + } + #endif /* increase absolute time */ moldyn->time+=moldyn->tau; @@@ -1738,6 -2238,59 +2249,59 @@@ } + /* 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; } @@@ -1756,15 -2309,47 +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 for(i=0;igvir),0,sizeof(t_virial)); /* reset force, site energy and virial of every atom */ + #ifdef PARALLEL + i=omp_get_thread_num(); + #pragma omp parallel for private(virial) + #endif for(i=0;ifunc1b) - 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; @@@ -1876,24 -2503,25 +2514,31 @@@ 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=(jsubcell->list[p]; #else this=&(neighbour_i[j]); list_reset_f(this); @@@ -1908,86 -2536,17 +2553,74 @@@ if(jtom==&(itom[i])) continue; + /* reset 3bp run */ + moldyn->run3bp=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 } - #else - } while(list_next_f(that)!=\ - L_NO_NEXT_ELEMENT); - #endif - - } - - /* finish of first j loop after first k loop */ - if(moldyn->func_j0e) - moldyn->func_j0e(moldyn, - &(itom[i]), - jtom, - bc_ij); - #ifdef STATIC_LISTS } + #elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@@ -1995,11 -2554,20 +2628,20 @@@ } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; + /* copy the neighbour lists */ + #ifdef STATIC_LISTS + /* no copy needed for static lists */ + #elif LOWMEM_LISTS + /* no copy needed for lowmem lists */ + #else + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); + #endif + /* second loop over atoms j */ for(j=0;j<27;j++) { @@@ -2007,10 -2575,17 +2649,17 @@@ #ifdef STATIC_LISTS p=0; - while(neighbour_i[j][p]!=0) { + while(neighbour_i[j][p]!=-1) { jtom=&(atom[neighbour_i[j][p]]); p++; + #elif LOWMEM_LISTS + p=neighbour_i[j]; + + while(p!=-1) { + + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour_i[j]); list_reset_f(this); @@@ -2032,18 -2607,18 +2681,18 @@@ /* 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++) { @@@ -2051,10 -2626,17 +2700,17 @@@ #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; + #elif LOWMEM_LISTS + q=neighbour_i[k]; + + while(q!=-1) { + + ktom=&(itom[q]); + q=lc->subcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@@ -2069,19 -2651,22 +2725,22 @@@ if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func_j1_k0(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func3b_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); + #ifdef STATIC_LISTS } + #elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@@ -2091,15 -2676,14 +2750,15 @@@ } - 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++) { @@@ -2107,10 -2691,17 +2766,17 @@@ #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; + #elif LOWMEM_LISTS + q=neighbour_i[k]; + + while(q!=-1) { + + ktom=&(itom[q]); + q=lc->subcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@@ -2125,20 -2716,22 +2791,22 @@@ 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 } + #elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@@ -2148,20 -2741,22 +2816,22 @@@ } - /* 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 } + #elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif } - + #ifdef DEBUG //printf("\n\n"); #endif @@@ -2182,6 -2777,9 +2852,9 @@@ #endif /* some postprocessing */ + #ifdef PARALLEL + #pragma omp parallel for + #endif for(i=0;igvir.xx+=itom[i].r.x*itom[i].f.x; @@@ -2193,7 -2791,7 +2866,7 @@@ /* check forces regarding the given timestep */ if(v3_norm(&(itom[i].f))>\ - 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) printf("[moldyn] WARNING: pfc (high force: atom %d)\n", i); } @@@ -2250,6 -2848,51 +2923,51 @@@ int check_per_bound(t_moldyn *moldyn,t_ 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 */ @@@ -2340,7 -2983,7 +3058,7 @@@ int moldyn_read_save_file(t_moldyn *mol 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); @@@ -2363,6 -3006,16 +3081,16 @@@ 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; @@@ -2394,13 -3047,16 +3122,16 @@@ int process_2b_bonds(t_moldyn *moldyn,v #ifdef STATIC_LISTS int *neighbour[27]; int p; + #elif LOWMEM_LISTS + int neighbour[27]; + int p; #else t_list neighbour[27]; + t_list *this; #endif u8 bc; t_atom *itom,*jtom; int i,j; - t_list *this; lc=&(moldyn->lc); itom=moldyn->atom; @@@ -2420,10 -3076,17 +3151,17 @@@ #ifdef STATIC_LISTS p=0; - while(neighbour[j][p]!=0) { + while(neighbour[j][p]!=-1) { jtom=&(moldyn->atom[neighbour[j][p]]); p++; + #elif LOWMEM_LISTS + p=neighbour[j]; + + while(p!=-1) { + + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour[j]); list_reset_f(this); @@@ -2441,6 -3104,8 +3179,8 @@@ #ifdef STATIC_LISTS } + #elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@@ -2451,6 -3116,84 +3191,84 @@@ } + /* + * function to find neighboured atoms + */ + + int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom, + int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom, + void *data,u8 bc)) { + + t_linkcell *lc; + #ifdef STATIC_LISTS + int *neighbour[27]; + int p; + #elif LOWMEM_LISTS + int neighbour[27]; + int p; + #else + t_list neighbour[27]; + t_list *this; + #endif + u8 bc; + t_atom *natom; + int j; + + lc=&(moldyn->lc); + + /* neighbour indexing */ + link_cell_neighbour_index(moldyn, + (atom->r.x+moldyn->dim.x/2)/lc->x, + (atom->r.y+moldyn->dim.y/2)/lc->x, + (atom->r.z+moldyn->dim.z/2)/lc->x, + neighbour); + + for(j=0;j<27;j++) { + + bc=(jdnlc)?0:1; + + #ifdef STATIC_LISTS + p=0; + + while(neighbour[j][p]!=-1) { + + natom=&(moldyn->atom[neighbour[j][p]]); + p++; + #elif LOWMEM_LISTS + p=neighbour[j]; + + while(p!=-1) { + + natom=&(moldyn->atom[p]); + p=lc->subcell->list[p]; + #else + this=&(neighbour[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + do { + + natom=this->current->data; + #endif + + /* process bond */ + process(moldyn,atom,natom,data,bc); + + #ifdef STATIC_LISTS + } + #elif LOWMEM_LISTS + } + #else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); + #endif + } + + return 0; + + } + /* * post processing functions */ @@@ -2485,6 -3228,7 +3303,7 @@@ int calculate_diffusion_coefficient(t_m int i; t_atom *atom; t_3dvec dist; + t_3dvec final_r; double d2; int a_cnt; int b_cnt; @@@ -2498,8 -3242,12 +3317,12 @@@ 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) { @@@ -2521,6 -3269,57 +3344,57 @@@ 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; @@@ -2666,9 -3465,8 +3540,8 @@@ int bond_analyze_process(t_moldyn *mold ba=data; /* increase total bond counter - * ... double counting! */ - ba->tcnt+=2; + ba->tcnt+=1; if(itom->brand==0) ba->acnt[jtom->tag]+=1; @@@ -2685,10 -3483,11 +3558,11 @@@ 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; @@@ -2708,9 -3507,9 +3582,9 @@@ 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; @@@ -2719,27 -3518,30 +3593,30 @@@ 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; @@@ -2777,7 -3579,11 +3654,11 @@@ int visual_bonds_process(t_moldyn *mold return 0; } + #ifdef VISUAL_THREAD + void *visual_atoms(void *ptr) { + #else int visual_atoms(t_moldyn *moldyn) { + #endif int i; char file[128+64]; @@@ -2786,6 -3592,12 +3667,12 @@@ t_visual *v; t_atom *atom; t_vb vb; + t_3dvec strain; + #ifdef VISUAL_THREAD + t_moldyn *moldyn; + + moldyn=ptr; + #endif v=&(moldyn->vis); dim.x=v->dim.x; @@@ -2795,31 -3607,39 +3682,39 @@@ help=(dim.x+dim.y); - sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time); + sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time); vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); if(vb.fd<0) { perror("open visual save file fd"); + #ifndef VISUAL_THREAD return -1; + #endif } /* write the actual data file */ // povray header - dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n", + dprintf(vb.fd,"# [P] %d %08.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++) { + 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 process_2b_bonds(moldyn,&vb,visual_bonds_process); + #endif // boundaries if(dim.x) { @@@ -2865,6 -3685,32 +3760,32 @@@ close(vb.fd); + #ifdef VISUAL_THREAD + pthread_exit(NULL); + + } + #else + + return 0; + } + #endif + + /* + * fpu cntrol functions + */ + + // set rounding to double (eliminates -ffloat-store!) + int fpu_set_rtd(void) { + + fpu_control_t ctrl; + + _FPU_GETCW(ctrl); + + ctrl&=~_FPU_EXTENDED; + ctrl|=_FPU_DOUBLE; + + _FPU_SETCW(ctrl); + return 0; } diff --combined moldyn.h index b2d6a34,a1866a9..0670f3d --- a/moldyn.h +++ b/moldyn.h @@@ -45,6 -45,7 +45,7 @@@ typedef struct s_atom u8 brand; /* brand id */ int tag; /* atom unique id (number of atom) */ u8 attr; /* attributes */ + int pbc[3]; /* pb crossing in x, y and z direction */ } t_atom; #define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ @@@ -58,6 -59,12 +59,12 @@@ #define DEFAULT_ATOM_ATTR 0x74 // 1,2,3 body interaction + visualize + /* special list structure for low mem approach */ + typedef struct s_lowmem_list { + int *head; + int *list; + } t_lowmem_list; + /* cell lists */ typedef struct s_linkcell { int nx,ny,nz; /* amount of cells in x, y and z direction */ @@@ -66,6 -73,8 +73,8 @@@ double x,y,z; /* the actual cell lengthes */ #ifdef STATIC_LISTS int **subcell; /* pointer to the cell lists */ + #elif LOWMEM_LISTS + t_lowmem_list *subcell; /* low mem approach list */ #else t_list *subcell; /* pointer to the cell lists */ #endif @@@ -104,18 -113,15 +113,18 @@@ typedef struct s_moldyn double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ /* potential force function and parameter pointers */ - int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); - int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_k1)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); - int (*func3b_k2)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_i0)(struct s_moldyn *moldyn,t_atom *ai); + int (*func_j0)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j0_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j0e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1c)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k1)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); void *pot_params; unsigned char run3bp; @@@ -132,14 -138,22 +141,22 @@@ double t_sum; /* sum over all t */ double t_avg; /* average value of t */ + /* for sale! */ t_virial gvir; /* global virial (absolute coordinates) */ - double gv; - double gv_sum; - double gv_avg; - - double gp; /* pressure computed from global virial */ - double gp_sum; /* sum over all gp */ - double gp_avg; /* average value of gp */ + //double gv; + //double gv_sum; + //double gv_avg; + double sale1; + double sale2; + double sale3; + + // gp stuff exchanged by kinetic energies + //double gp; /* pressure computed from global virial */ + //double gp_sum; /* sum over all gp */ + //double gp_avg; /* average value of gp */ + double ekinx; + double ekiny; + double ekinz; t_virial vir; /* actual virial */ double virial; @@@ -245,6 -259,43 +262,43 @@@ typedef struct s_vb int fd; } t_vb; + typedef struct s_part_params { + u8 type; + double r; + t_3dvec p; + t_3dvec d; + } t_part_params; + + #define PART_INSIDE_R 1 + #define PART_OUTSIDE_R 2 + #define PART_INSIDE_D 3 + #define PART_OUTSIDE_D 4 + + typedef struct s_defect_params { + u8 type; + u8 stype; + double od; + double dd; + int element; + u8 brand; + u8 attr; + } t_defect_params; + + #define DEFECT_TYPE_0D 1 + #define DEFECT_TYPE_1D 2 + #define DEFECT_TYPE_2D 3 + #define DEFECT_TYPE_3D 4 + + #define DEFECT_STYPE_DB_X 1 + #define DEFECT_STYPE_DB_Y 2 + #define DEFECT_STYPE_DB_Z 3 + #define DEFECT_STYPE_DB_R 4 + + typedef struct s_offset_params { + t_3dvec o; + u8 use; + } t_offset_params; + /* * * defines @@@ -307,7 -358,6 +361,7 @@@ #define MOLDYN_POTENTIAL_LJ 0x01 #define MOLDYN_POTENTIAL_TM 0x02 #define MOLDYN_POTENTIAL_AM 0x03 +#define MOLDYN_POTENTIAL_AO 0x04 #define LOG_TOTAL_ENERGY 0x01 #define LOG_TOTAL_MOMENTUM 0x02 @@@ -356,12 -406,13 +410,13 @@@ #define FCC 0x02 #define DIAMOND 0x04 #define ZINCBLENDE 0x08 + #define NONE 0x80 /* * more includes */ - #include "pse.h" + //#include "pse.h" /* * @@@ -391,14 -442,19 +446,19 @@@ int moldyn_set_report(t_moldyn *moldyn, int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); int moldyn_log_shutdown(t_moldyn *moldyn); - 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 add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + 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 add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, t_3dvec *r,t_3dvec *v); int del_atom(t_moldyn *moldyn,int tag); - int cubic_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); - int diamond_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 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 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 destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); @@@ -427,6 -483,8 +487,8 @@@ int link_cell_init(t_moldyn *moldyn,u8 int link_cell_update(t_moldyn *moldyn); #ifdef STATIC_LISTS int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell); + #elif LOWMEM_LISTS + int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int *cell); #else int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell); #endif @@@ -445,6 -503,7 +507,7 @@@ int virial_calc(t_atom *a,t_3dvec *f,t_ //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) // __attribute__((always_inline)); int check_per_bound(t_moldyn *moldyn,t_3dvec *a); + int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a); //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) // __attribute__((always_inline)); @@@ -456,10 -515,15 +519,15 @@@ int moldyn_load(t_moldyn *moldyn) int process_2b_bonds(t_moldyn *moldyn,void *data, int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, void *data,u8 bc)); + int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom, + int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom, + void *data,u8 bc)); + int get_line(int fd,char *line,int max); int pair_correlation_init(t_moldyn *moldyn,double dr); int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc); + int calculate_msd(t_moldyn *moldyn,double *msd); int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, t_atom *jtom,void *data,u8 bc); int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr); @@@ -470,7 -534,13 +538,13 @@@ int bond_analyze(t_moldyn *moldyn,doubl int visual_init(t_moldyn *moldyn,char *filebase); int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, void *data,u8 bc); + #ifdef VISUAL_THREAD + void *visual_atoms(void *ptr); + #else int visual_atoms(t_moldyn *moldyn); + #endif + + int fpu_set_rtd(void); #endif diff --combined potentials/albe.c index 9d8336e,1ec3938..9f98282 --- a/potentials/albe.c +++ b/potentials/albe.c @@@ -53,10 -53,7 +53,10 @@@ int albe_mult_set_params(t_moldyn *mold p->mu[0]=ALBE_MU_SI; p->gamma[0]=ALBE_GAMMA_SI; p->c[0]=ALBE_C_SI; + p->c2[0]=p->c[0]*p->c[0]; p->d[0]=ALBE_D_SI; + p->d2[0]=p->d[0]*p->d[0]; + p->c2d2[0]=p->c2[0]/p->d2[0]; p->h[0]=ALBE_H_SI; switch(element2) { case C: @@@ -70,10 -67,7 +70,10 @@@ p->mu[1]=ALBE_MU_C; p->gamma[1]=ALBE_GAMMA_C; p->c[1]=ALBE_C_C; + p->c2[1]=p->c[1]*p->c[1]; p->d[1]=ALBE_D_C; + p->d2[1]=p->d[1]*p->d[1]; + p->c2d2[1]=p->c2[1]/p->d2[1]; p->h[1]=ALBE_H_C; /* mixed type: silicon carbide */ p->Smixed=ALBE_S_SIC; @@@ -85,10 -79,7 +85,10 @@@ p->mu_m=ALBE_MU_SIC; p->gamma_m=ALBE_GAMMA_SIC; p->c_mixed=ALBE_C_SIC; + p->c2_mixed=p->c_mixed*p->c_mixed; p->d_mixed=ALBE_D_SIC; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2_m=p->c2_mixed/p->d2_mixed; p->h_mixed=ALBE_H_SIC; break; default: @@@ -105,6 -96,15 +105,15 @@@ p->S2[0]=p->S[0]*p->S[0]; p->S2[1]=p->S[1]*p->S[1]; p->S2mixed=p->Smixed*p->Smixed; + p->c2[0]=p->c[0]*p->c[0]; + p->c2[1]=p->c[1]*p->c[1]; + p->c2_mixed=p->c_mixed*p->c_mixed; + p->d2[0]=p->d[0]*p->d[0]; + p->d2[1]=p->d[1]*p->d[1]; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2[0]=p->c2[0]/p->d2[0]; + p->c2d2[1]=p->c2[1]/p->d2[1]; + p->c2d2_m=p->c2_mixed/p->d2_mixed; printf("[albe] mult parameter info:\n"); printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed); @@@ -114,56 -114,34 +123,56 @@@ printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], p->lambda_m); printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); - printf(" gamma | %f | %f\n",p->gamma[0],p->gamma[1]); - printf(" c | %f | %f\n",p->c[0],p->c[1]); - printf(" d | %f | %f\n",p->d[0],p->d[1]); - printf(" h | %f | %f\n",p->h[0],p->h[1]); + printf(" gamma | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m); + printf(" c | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed); + printf(" d | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed); + printf(" c2 | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed); + printf(" d2 | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed); + printf(" c2d2 | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m); + printf(" h | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed); return 0; } -/* albe 3 body potential function (first ij loop) */ -int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { +/* first i loop */ +int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) { t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double S2; - t_3dvec dist_ij; - double d_ij2,d_ij; + + int i; params=moldyn->pot_params; exchange=&(params->exchange); - /* reset zeta sum */ - exchange->zeta_ij=0.0; + /* zero exchange values */ + memset(exchange->zeta,0,ALBE_MAXN*sizeof(double)); + for(i=0;idzeta[i],0,ALBE_MAXN*sizeof(t_3dvec)); + exchange->jcnt=0; + exchange->j2cnt=0; - /* - * set ij depending values - */ + return 0; +} +/* first j loop within first i loop */ +int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + double S2,S,R,d2,d,s_r,arg; + t_3dvec dist; + int j; + u8 brand; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* get j counter */ + j=exchange->jcnt; + + /* set ij depending values */ brand=ai->brand; if(brand==aj->brand) { S2=params->S2[brand]; @@@ -173,226 -151,166 +182,204 @@@ } /* dist_ij, d_ij2 */ - v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij2=v3_absolute_square(&dist_ij); + v3_sub(&dist,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist); + exchange->dist[j]=dist; + d2=v3_absolute_square(&dist); + exchange->d2[j]=d2; /* if d_ij2 > S2 => no force & potential energy contribution */ - if(d_ij2>S2) { + if(d2>S2) { moldyn->run3bp=0; + exchange->skip[j]=1; + exchange->jcnt+=1; return 0; } + exchange->skip[j]=0; + + /* more ij depending values */ + if(brand==aj->brand) { + R=params->R[brand]; + S=params->S[brand]; + /* albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma[brand]); + exchange->c_[j]=&(params->c[brand]); + exchange->d_[j]=&(params->d[brand]); + exchange->h_[j]=&(params->h[brand]); + exchange->c2_[j]=&(params->c2[brand]); + exchange->d2_[j]=&(params->d2[brand]); + exchange->c2d2_[j]=&(params->c2d2[brand]); + } + else { + R=params->Rmixed; + S=params->Smixed; + /* albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma_m); + exchange->c_[j]=&(params->c_mixed); + exchange->d_[j]=&(params->d_mixed); + exchange->h_[j]=&(params->h_mixed); + exchange->c2_[j]=&(params->c2_mixed); + exchange->d2_[j]=&(params->d2_mixed); + exchange->c2d2_[j]=&(params->c2d2_m); + } /* d_ij */ - d_ij=sqrt(d_ij2); + d=sqrt(d2); + exchange->d[j]=d; + + /* f_c, df_c */ + if(df_c[j]=1.0; + exchange->df_c[j]=0.0; + } + else { + s_r=S-R; + arg=M_PI*(d-R)/s_r; + exchange->f_c[j]=0.5+0.5*cos(arg); + exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d)); + } - /* store values */ - exchange->dist_ij=dist_ij; - exchange->d_ij2=d_ij2; - exchange->d_ij=d_ij; + /* reset k counter */ + exchange->kcnt=0; - /* reset k counter for first k loop */ - exchange->kcount=0; - return 0; } -/* albe 3 body potential function (first k loop) */ -int albe_mult_3bp_k1(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { +/* first k loop within first j loop within first i loop */ +int albe_mult_i0_j0_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double R,S,S2; - t_3dvec dist_ij,dist_ik; - double d_ik2,d_ik,d_ij; - double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg; - double f_c_ik,df_c_ik; - int kcount; + + int j,k; + t_3dvec distj,distk; + double dj,dk,djdk_inv,cos_theta; + double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j; + double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k; + t_3dvec dcosdrj,dcosdrk,tmp; + t_3dvec *dzjj,*dzkk,*dzjk,*dzkj; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - if(kcount>ALBE_MAXN) { - printf("FATAL: neighbours = %d\n",kcount); - printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag); + if(aj==ak) { + exchange->kcnt+=1; + return 0; } - /* ik constants */ - brand=ai->brand; - if(brand==ak->brand) { - R=params->R[brand]; - S=params->S[brand]; - S2=params->S2[brand]; - /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma[brand]); - exchange->c_i=&(params->c[brand]); - exchange->d_i=&(params->d[brand]); - exchange->h_i=&(params->h[brand]); - } - else { - R=params->Rmixed; - S=params->Smixed; - S2=params->S2mixed; - /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma_m); - exchange->c_i=&(params->c_mixed); - exchange->d_i=&(params->d_mixed); - exchange->h_i=&(params->h_mixed); + /* kjcnt; + k=exchange->kcnt; + if(k>=ALBE_MAXN) { + printf("FATAL: too many neighbours! (%d)\n",k); + printf(" atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag); } - exchange->ci2=*(exchange->c_i)**(exchange->c_i); - exchange->di2=*(exchange->d_i)**(exchange->d_i); - exchange->ci2di2=exchange->ci2/exchange->di2; - - /* dist_ik, d_ik2 */ - v3_sub(&dist_ik,&(ak->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ik); - d_ik2=v3_absolute_square(&dist_ik); - - /* store data for second k loop */ - exchange->dist_ik[kcount]=dist_ik; - exchange->d_ik2[kcount]=d_ik2; - - /* return if not within cutoff */ - if(d_ik2>S2) { - exchange->kcount++; + if((k>=j)|(exchange->skip[k])) { + exchange->kcnt+=1; return 0; } - /* d_ik */ - d_ik=sqrt(d_ik2); - - /* dist_ij, d_ij */ - dist_ij=exchange->dist_ij; - d_ij=exchange->d_ij; + /* distances */ + distj=exchange->dist[j]; + distk=exchange->dist[k]; + dj=exchange->d[j]; + dk=exchange->d[k]; + djdk_inv=1.0/(dj*dk); /* cos theta */ - cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - - /* g_ijk */ - h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism - d2_h_cos2=exchange->di2+(h_cos*h_cos); - frac=exchange->ci2/d2_h_cos2; - g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac); - dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. - - /* zeta sum += f_c_ik * g_ijk */ - if(d_ik<=R) { - exchange->zeta_ij+=g; - f_c_ik=1.0; - df_c_ik=0.0; + cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv; + + /* g(cos(theta)) ij and ik values */ + h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism + d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j); + frac_j=*exchange->c2_[j]/d2_h_cos2_j; + gj=1.0+*exchange->c2d2_[j]-frac_j; + gj*=*(exchange->gamma_[j]); + dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe + if(ak->brand==aj->brand) { + gk=gj; + dgk=dgj; } else { - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); - exchange->zeta_ij+=f_c_ik*g; + h_cos_k=*(exchange->h_[k])+cos_theta; + d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k); + frac_k=*exchange->c2_[k]/d2_h_cos2_k; + gk=1.0+*exchange->c2d2_[k]-frac_k; + gk*=*(exchange->gamma_[k]); + dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k; } - /* zeta - for albe: ik depending g function */ - //if(ai->tag==0) { - // printf("------> %.15f %.15f\n",dj,dk); - // printf("------> %.15f %.15f\n",dj,dk); - //} - - exchange->zeta[j]+=(exchange->f_c[k]*gk); - exchange->zeta[k]+=(exchange->f_c[j]*gj); - - /* cos theta derivatives */ - v3_scale(&dcosdrj,&distk,djdk_inv); // j - v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]); - v3_add(&dcosdrj,&dcosdrj,&tmp); - v3_scale(&dcosdrk,&distj,djdk_inv); // k - v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]); - v3_add(&dcosdrk,&dcosdrk,&tmp); + #ifdef DEBUG + if(ai==&(moldyn->atom[DATOM])) + printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik); + #endif - /* zeta derivatives */ - dzjj=&(exchange->dzeta[j][j]); - dzkk=&(exchange->dzeta[k][k]); - dzjk=&(exchange->dzeta[j][k]); - dzkj=&(exchange->dzeta[k][j]); - v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk); - v3_add(dzjj,dzjj,&tmp); // j j - v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj); - v3_add(dzkk,dzkk,&tmp); // k k - v3_scale(&tmp,&distk,-exchange->df_c[k]*gk); // dk rik = - di rik - v3_add(dzjk,dzjk,&tmp); - v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk); - v3_add(dzjk,dzjk,&tmp); // j k - v3_scale(&tmp,&distj,-exchange->df_c[j]*gj); // dj rij = - di rij - v3_add(dzkj,dzkj,&tmp); - v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj); - v3_add(dzkj,dzkj,&tmp); // k j + /* store even more data for second k loop */ + exchange->g[kcount]=g; + exchange->dg[kcount]=dg; + exchange->d_ik[kcount]=d_ik; + exchange->cos_theta[kcount]=cos_theta; + exchange->f_c_ik[kcount]=f_c_ik; + exchange->df_c_ik[kcount]=df_c_ik; /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; + + return 0; +} + +/* first j loop within first i loop */ +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->jcnt+=1; return 0; } -int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { +/* second j loop within first i loop */ +int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - t_3dvec force; - double f_a,df_a,b,db,f_c,df_c; - double f_r,df_r; - double scale; - double mu,B; - double lambda,A; - double d_ij,r0; - unsigned char brand; - double S,R,s_r,arg; + + int j; + double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db; + double A,B,mu,lambda,r0; double energy; + t_3dvec *dist,force; + double scale; + u8 brand; params=moldyn->pot_params; exchange=&(params->exchange); + /* get j counter */ + j=exchange->j2cnt; + + /* skip if j not within cutoff */ + if(exchange->skip[j]) { + moldyn->run3bp=0; + exchange->j2cnt+=1; + return 0; + } + + /* distance */ + d=exchange->d[j]; + dist=&(exchange->dist[j]); + f_c=exchange->f_c[j]; + df_c=exchange->df_c[j]; + + /* determine parameters to calculate fa, dfa, fr, dfr */ brand=aj->brand; if(brand==ai->brand) { - S=params->S[brand]; - R=params->R[brand]; B=params->B[brand]; A=params->A[brand]; r0=params->r0[brand]; @@@ -400,6 -318,8 +387,6 @@@ lambda=params->lambda[brand]; } else { - S=params->Smixed; - R=params->Rmixed; B=params->Bmixed; A=params->Amixed; r0=params->r0_mixed; @@@ -407,151 -327,206 +394,211 @@@ lambda=params->lambda_m; } - d_ij=exchange->d_ij; - - /* f_c, df_c */ - if(d_ijzeta_ij==0.0) { - b=1.0; - db=0.0; - } - else { - b=1.0/sqrt(1.0+exchange->zeta_ij); - db=-0.5*b/(1.0+exchange->zeta_ij); - } + b=1.0/sqrt(1.0+exchange->zeta[j]); + db=-0.5*b/(1.0+exchange->zeta[j]); + + /* energy contribution */ + energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism + moldyn->energy+=energy; + ai->e+=energy; - /* force contribution for atom i */ + /* force contribution for atom i due to ij bond */ scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism - v3_scale(&force,&(exchange->dist_ij),scale); + v3_scale(&force,dist,scale); v3_add(&(ai->f),&(ai->f),&force); - /* force contribution for atom j */ +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif + + /* force contribution for atom j due to ij bond */ v3_scale(&force,&force,-1.0); // dri rij = - drj rij v3_add(&(aj->f),&(aj->f),&force); - #ifdef NDEBUG - if(aj->tag==0) { - printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]); - printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); - } + /* virial */ + virial_calc(ai,&force,&(exchange->dist_ij)); + + #ifdef DEBUG + if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) { + printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + if(ai==&(moldyn->atom[DATOM])) + printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + if(aj==&(moldyn->atom[DATOM])) + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), + f_c,b,f_a,f_r); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); + } #endif + /* virial */ + virial_calc(ai,&force,dist); + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ exchange->pre_dzeta=0.5*f_a*f_c*db; - /* energy contribution */ - energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism - moldyn->energy+=energy; - ai->e+=energy; + /* force contribution (drj derivative) */ + v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta); + v3_add(&(aj->f),&(aj->f),&force); + +#ifdef NDEBUG +if(aj->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag); +printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); +} +#endif + + /* virial */ + virial_calc(ai,&force,dist); + + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif /* reset k counter for second k loop */ - exchange->kcount=0; + exchange->kcnt=0; return 0; } -/* albe 3 body potential function (second k loop) */ -int albe_mult_3bp_k2(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { +/* second k loop within second j loop within first i loop */ +int albe_mult_i0_j2_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - int kcount; - t_3dvec dist_ik,dist_ij; - double d_ik2,d_ik,d_ij2,d_ij; - unsigned char brand; - double S2; - double g,dg,cos_theta; - double pre_dzeta; - double f_c_ik,df_c_ik; - double dijdik_inv,fcdg,dfcg; - t_3dvec dcosdrj,dcosdrk; - t_3dvec force,tmp; + + int j,k; + t_3dvec force; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - - if(kcount>ALBE_MAXN) - printf("FATAL: neighbours!\n"); - - /* d_ik2 */ - d_ik2=exchange->d_ik2[kcount]; - brand=ak->brand; - if(brand==ai->brand) - S2=params->S2[brand]; - else - S2=params->S2mixed; - - /* return if d_ik > S */ - if(d_ik2>S2) { - exchange->kcount++; + if(aj==ak) { + exchange->kcnt+=1; return 0; } - /* k!=j & check whether to run k */ - k=exchange->kcnt; - j=exchange->j2cnt; - if((k==j)|(exchange->skip[k])) { - exchange->kcnt+=1; - return 0; - } - - /* force contribution (drk derivative) */ - v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta); - v3_add(&(ak->f),&(ak->f),&force); + /* prefactor dzeta */ + pre_dzeta=exchange->pre_dzeta; - #ifdef NDEBUG - if(ak->tag==0) { - printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag); - printf(" t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z); - } + /* dist_ik, d_ik */ + dist_ik=exchange->dist_ik[kcount]; + d_ik=exchange->d_ik[kcount]; + + /* f_c_ik, df_c_ik */ + f_c_ik=exchange->f_c_ik[kcount]; + df_c_ik=exchange->df_c_ik[kcount]; + + /* dist_ij, d_ij2, d_ij */ + dist_ij=exchange->dist_ij; + d_ij2=exchange->d_ij2; + d_ij=exchange->d_ij; + + /* g, dg, cos_theta */ + g=exchange->g[kcount]; + dg=exchange->dg[kcount]; + cos_theta=exchange->cos_theta[kcount]; + + /* cos_theta derivatives wrt j,k */ + dijdik_inv=1.0/(d_ij*d_ik); + v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j + v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2); + v3_add(&dcosdrj,&dcosdrj,&tmp); + v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k + v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); + v3_add(&dcosdrk,&dcosdrk,&tmp); + + /* f_c_ik * dg, df_c_ik * g */ + fcdg=f_c_ik*dg; + dfcg=df_c_ik*g; + + /* derivative wrt j */ + v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); + + /* force contribution */ + v3_add(&(aj->f),&(aj->f),&force); + + #ifdef DEBUG + if(aj==&(moldyn->atom[DATOM])) { + printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } #endif /* virial */ - virial_calc(ai,&force,&(exchange->dist[k])); + virial_calc(ai,&force,&dist_ij); + /* force contribution to atom i */ v3_scale(&force,&force,-1.0); v3_add(&(ai->f),&(ai->f),&force); - #ifdef NDEBUG - if(ai->tag==0) { - printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k); - printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); - printf(" ## %f\n",exchange->d[k]); - } + /* derivative wrt k */ + v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik + v3_scale(&tmp,&dcosdrk,fcdg); + v3_add(&force,&force,&tmp); + v3_scale(&force,&force,pre_dzeta); + - /* force contribution */ - v3_add(&(ak->f),&(ak->f),&force); ++ v3_scale(&force,&force,-1.0); ++ v3_add(&(ai->f),&(ai->f),&force); + + #ifdef DEBUG + if(ak==&(moldyn->atom[DATOM])) { + printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } #endif + /* virial */ + virial_calc(ai,&force,&dist_ik); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; return 0; +} + +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->j2cnt+=1; + + return 0; } int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) { diff --combined potentials/albe.h index 39c4df8,89cdae2..9687655 --- a/potentials/albe.h +++ b/potentials/albe.h @@@ -8,36 -8,39 +8,36 @@@ #ifndef ALBE_H #define ALBE_H -#define ALBE_MAXN 16*27 +#define ALBE_MAXN (4*27) /* albe exchange type */ typedef struct s_albe_exchange { - t_3dvec dist_ij; - double d_ij2; - double d_ij; + t_3dvec dist[ALBE_MAXN]; + double d2[ALBE_MAXN]; + double d[ALBE_MAXN]; - t_3dvec dist_ik[ALBE_MAXN]; - double d_ik2[ALBE_MAXN]; - double d_ik[ALBE_MAXN]; + double f_c[ALBE_MAXN]; + double df_c[ALBE_MAXN]; - double f_c_ik[ALBE_MAXN]; - double df_c_ik[ALBE_MAXN]; + double zeta[ALBE_MAXN]; + t_3dvec dzeta[ALBE_MAXN][ALBE_MAXN]; - double g[ALBE_MAXN]; - double dg[ALBE_MAXN]; - double cos_theta[ALBE_MAXN]; + u8 skip[ALBE_MAXN]; - double *gamma_i; - double *c_i; - double *d_i; - double *h_i; + double *gamma_[ALBE_MAXN]; + double *c_[ALBE_MAXN]; + double *d_[ALBE_MAXN]; + double *c2_[ALBE_MAXN]; + double *d2_[ALBE_MAXN]; + double *c2d2_[ALBE_MAXN]; + double *h_[ALBE_MAXN]; - double ci2; - double di2; - double ci2di2; - - double zeta_ij; double pre_dzeta; - int kcount; + int jcnt; + int j2cnt; + int kcnt; } t_albe_exchange; /* albe mult (2!) potential parameters */ @@@ -62,33 -65,32 +62,34 @@@ typedef struct s_albe_mult_params double gamma[2]; double gamma_m; double c[2]; + double c2[2]; double c_mixed; + double c2[2]; double c2_mixed; double d[2]; + double d2[2]; double d_mixed; + double d2[2]; double d2_mixed; - double c2d2[2]; - double c2d2_m; double h[2]; double h_mixed; + double c2d2[2]; + double c2d2_m; t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */ } t_albe_mult_params; /* function prototypes */ int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2); - int albe_mult_i0(t_moldyn *moldyn,t_atom *ai); - int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int albe_mult_i0_j0_k0(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); - int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int albe_mult_i0_j2_k0(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); - int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int albe_mult_3bp_k1(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); + int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int albe_mult_3bp_k2(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc); + /* fast */ + int albe_potential_force_calc(t_moldyn *moldyn); /* albe potential parameter defines */ @@@ -102,35 -104,35 +103,35 @@@ #define ALBE_MU_SI (1.4761*sqrt(2.0/1.842)) #define ALBE_GAMMA_SI 0.114354 #define ALBE_C_SI 2.00494 - #define ALBE_D_SI 0.81472 + #define ALBE_D_SI 0.814719 #define ALBE_H_SI 0.259 #define ALBE_LC_SI 5.429 // carbon #define ALBE_R_C (2.00-0.15) #define ALBE_S_C (2.00+0.15) - #define ALBE_A_C (6.00*EV/1.167) - #define ALBE_B_C (-2.167*6.00*EV/1.167) - #define ALBE_R0_C 1.4276 - #define ALBE_LAMBDA_C (2.0099*sqrt(2.0*2.167)) - #define ALBE_MU_C (2.0099*sqrt(2.0/2.167)) - #define ALBE_GAMMA_C 0.11233 - #define ALBE_C_C 181.910 - #define ALBE_D_C 6.28433 - #define ALBE_H_C 0.5556 + #define ALBE_A_C (6.00*EV/1.1671419) + #define ALBE_B_C (-2.1671419*6.00*EV/1.1671419) + #define ALBE_R0_C 1.4276442 + #define ALBE_LAMBDA_C (2.0099457*sqrt(2.0*2.1671419)) + #define ALBE_MU_C (2.0099457*sqrt(2.0/2.1671419)) + #define ALBE_GAMMA_C 0.1123327 + #define ALBE_C_C 181.9100526 + #define ALBE_D_C 6.2843249 + #define ALBE_H_C 0.5556181 #define ALBE_LC_C 3.566 // mixed: silicon carbide #define ALBE_R_SIC (2.40-0.20) #define ALBE_S_SIC (2.40+0.20) - #define ALBE_A_SIC (4.36*EV/0.847) - #define ALBE_B_SIC (-1.847*4.36*EV/0.847) + #define ALBE_A_SIC (4.36*EV/0.8474739) + #define ALBE_B_SIC (-1.8474739*4.36*EV/0.8474739) #define ALBE_R0_SIC 1.79 - #define ALBE_LAMBDA_SIC (1.6991*sqrt(2.0*1.847)) - #define ALBE_MU_SIC (1.6991*sqrt(2.0/1.847)) - #define ALBE_GAMMA_SIC 0.011877 - #define ALBE_C_SIC 273987 - #define ALBE_D_SIC 180.314 + #define ALBE_LAMBDA_SIC (1.6990751*sqrt(2.0*1.8474739)) + #define ALBE_MU_SIC (1.6990751*sqrt(2.0/1.8474739)) + #define ALBE_GAMMA_SIC 0.0118769 + #define ALBE_C_SIC 273986.61 + #define ALBE_D_SIC 180.31411 #define ALBE_H_SIC 0.68 #define ALBE_LC_SIC 4.359