X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=d33de02efa776fa52dd3a10bbd5ca00d93ba7b4c;hb=15b4727e1137600f8f46af027aefd2b5c7a56420;hp=716a087fa90f1b89b073411af22e2d55995aba8a;hpb=15dd7ee084a5db24a3e1a15428cd47e8d61fda5b;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 716a087..d33de02 100644 --- a/moldyn.c +++ b/moldyn.c @@ -17,20 +17,8 @@ #include "moldyn.h" -#include "math/math.h" -#include "init/init.h" -#include "random/random.h" -#include "visual/visual.h" -#include "list/list.h" - - int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { - //int ret; - - //ret=moldyn_parse_argv(moldyn,argc,argv); - //if(ret<0) return ret; - memset(moldyn,0,sizeof(t_moldyn)); rand_init(&(moldyn->random),NULL,1); @@ -78,6 +66,13 @@ int set_temperature(t_moldyn *moldyn,double t_ref) { return 0; } +int set_pressure(t_moldyn *moldyn,double p_ref) { + + moldyn->p_ref=p_ref; + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); @@ -93,12 +88,21 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->dim.y=y; moldyn->dim.z=z; + moldyn->volume=x*y*z; + if(visualize) { moldyn->vis.dim.x=x; moldyn->vis.dim.y=y; moldyn->vis.dim.z=z; } + printf("[moldyn] dimensions in A and A^3 respectively:\n"); + printf(" x: %f\n",moldyn->dim.x); + printf(" y: %f\n",moldyn->dim.y); + printf(" z: %f\n",moldyn->dim.z); + printf(" volume: %f\n",moldyn->volume); + printf(" visualize simulation box: %s\n",visualize?"on":"off"); + return 0; } @@ -141,7 +145,7 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) { - moldyn->func2b=func; + moldyn->func2b_post=func; moldyn->pot2b_params=params; return 0; @@ -155,38 +159,56 @@ int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { return 0; } -int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { +int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { + + strncpy(moldyn->vlsdir,dir,127); + + return 0; +} + +int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { + + char filename[128]; + int ret; switch(type) { case LOG_TOTAL_ENERGY: moldyn->ewrite=timer; - moldyn->efd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + snprintf(filename,127,"%s/energy",moldyn->vlsdir); + moldyn->efd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); if(moldyn->efd<0) { - perror("[moldyn] efd open"); + perror("[moldyn] energy log fd open"); return moldyn->efd; } dprintf(moldyn->efd,"# total energy log file\n"); break; case LOG_TOTAL_MOMENTUM: moldyn->mwrite=timer; - moldyn->mfd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + snprintf(filename,127,"%s/momentum",moldyn->vlsdir); + moldyn->mfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); if(moldyn->mfd<0) { - perror("[moldyn] mfd open"); + perror("[moldyn] momentum log fd open"); return moldyn->mfd; } dprintf(moldyn->efd,"# total momentum log file\n"); break; case SAVE_STEP: moldyn->swrite=timer; - strncpy(moldyn->sfb,fb,63); break; case VISUAL_STEP: moldyn->vwrite=timer; - strncpy(moldyn->vfb,fb,63); - visual_init(&(moldyn->vis),fb); + ret=visual_init(&(moldyn->vis),moldyn->vlsdir); + if(ret<0) { + printf("[moldyn] visual init failure\n"); + return ret; + } break; default: - printf("unknown log mechanism: %02x\n",type); + printf("[moldyn] unknown log mechanism: %02x\n",type); return -1; } @@ -203,33 +225,43 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { return 0; } +/* + * creating lattice functions + */ + int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, - u8 attr,u8 bnum,int a,int b,int c) { + u8 attr,u8 brand,int a,int b,int c) { - int count; + int new,count; int ret; t_3dvec origin; + void *ptr; + t_atom *atom; - count=a*b*c; - - if(type==FCC) count*=4; + new=a*b*c; + count=moldyn->count; - if(type==DIAMOND) count*=8; + /* how many atoms do we expect */ + if(type==FCC) new*=4; + if(type==DIAMOND) new*=8; - moldyn->atom=malloc(count*sizeof(t_atom)); - if(moldyn->atom==NULL) { - perror("malloc (atoms)"); + /* allocate space for atoms */ + ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (create lattice)"); return -1; } - + moldyn->atom=ptr; + atom=&(moldyn->atom[count]); + v3_zero(&origin); switch(type) { case FCC: - ret=fcc_init(a,b,c,lc,moldyn->atom,&origin); + ret=fcc_init(a,b,c,lc,atom,&origin); break; case DIAMOND: - ret=diamond_init(a,b,c,lc,moldyn->atom,&origin); + ret=diamond_init(a,b,c,lc,atom,&origin); break; default: printf("unknown lattice type (%02x)\n",type); @@ -237,30 +269,114 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } /* debug */ - if(ret!=count) { - printf("ok, there is something wrong ...\n"); - printf("calculated -> %d atoms\n",count); - printf("created -> %d atoms\n",ret); + if(ret!=new) { + printf("[moldyn] creating lattice failed\n"); + printf(" amount of atoms\n"); + printf(" - expected: %d\n",new); + printf(" - created: %d\n",ret); return -1; } - moldyn->count=count; - printf("[moldyn] created lattice with %d atoms\n",count); + moldyn->count+=new; + printf("[moldyn] created lattice with %d atoms\n",new); - while(count) { - count-=1; - moldyn->atom[count].element=element; - moldyn->atom[count].mass=mass; - moldyn->atom[count].attr=attr; - moldyn->atom[count].bnum=bnum; - check_per_bound(moldyn,&(moldyn->atom[count].r)); + for(ret=0;retatom; - count=++(moldyn->count); + count=(moldyn->count)++; - ptr=realloc(atom,count*sizeof(t_atom)); + ptr=realloc(atom,(count+1)*sizeof(t_atom)); if(!ptr) { perror("[moldyn] realloc (add atom)"); return -1; @@ -278,12 +394,13 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr, moldyn->atom=ptr; atom=moldyn->atom; - atom[count-1].r=*r; - atom[count-1].v=*v; - atom[count-1].element=element; - atom[count-1].mass=mass; - atom[count-1].bnum=bnum; - atom[count-1].attr=attr; + atom[count].r=*r; + atom[count].v=*v; + atom[count].element=element; + atom[count].mass=mass; + atom[count].brand=brand; + atom[count].tag=count; + atom[count].attr=attr; return 0; } @@ -343,6 +460,28 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) { return 0; } +double temperature_calc(t_moldyn *moldyn) { + + double double_ekin; + int i; + t_atom *atom; + + atom=moldyn->atom; + + for(i=0;icount;i++) + double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v)); + + /* kinetic energy = 3/2 N k_B T */ + moldyn->t=double_ekin/(3.0*K_BOLTZMANN*moldyn->count); + + return moldyn->t; +} + +double get_temperature(t_moldyn *moldyn) { + + return moldyn->t; +} + int scale_velocity(t_moldyn *moldyn,u8 equi_init) { int i; @@ -361,18 +500,21 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { count=0; for(i=0;icount;i++) { if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) { - e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + e+=atom[i].mass*v3_absolute_square(&(atom[i].v)); count+=1; } } - if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN); + e*=0.5; + if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN); else return 0; /* no atoms involved in scaling! */ /* (temporary) hack for e,t = 0 */ if(e==0.0) { moldyn->t=0.0; - if(moldyn->t_ref!=0.0) + if(moldyn->t_ref!=0.0) { thermal_init(moldyn,equi_init); + return 0; + } else return 0; /* no scaling needed */ } @@ -384,17 +526,109 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { scale*=2.0; else if(moldyn->pt_scale&T_SCALE_BERENDSEN) - scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc; + scale=1.0+(scale-1.0)/moldyn->t_tc; scale=sqrt(scale); /* velocity scaling */ - for(i=0;icount;i++) + for(i=0;icount;i++) { if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) v3_scale(&(atom[i].v),&(atom[i].v),scale); + } return 0; } +double pressure_calc(t_moldyn *moldyn) { + + int i; + t_atom *atom; + double p1,p2,p=0; + + for(i=0;icount;i++) { + + + } + + p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1); + p1/=moldyn->volume; + + p2=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt2); + p2/=moldyn->volume; + + printf("compare pressures: %f %f\n",p1/ATM,p2/ATM); + + return moldyn->p; +} + +double get_pressure(t_moldyn *moldyn) { + + return moldyn->p; + +} + +int scale_volume(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim,*vdim; + double scale,v; + t_virial virial; + t_linkcell *lc; + int i; + + atom=moldyn->atom; + dim=&(moldyn->dim); + vdim=&(moldyn->vis.dim); + lc=&(moldyn->lc); + + memset(&virial,0,sizeof(t_virial)); + + for(i=0;icount;i++) { + virial.xx+=atom[i].virial.xx; + virial.yy+=atom[i].virial.yy; + virial.zz+=atom[i].virial.zz; + virial.xy+=atom[i].virial.xy; + virial.xz+=atom[i].virial.xz; + virial.yz+=atom[i].virial.yz; + } + + /* just a guess so far ... */ + v=virial.xx+virial.yy+virial.zz; + +printf("%f\n",v); + /* get pressure from virial */ + moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v; + moldyn->p/=moldyn->volume; +printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM); + + /* scale factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) + scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc); + else + /* should actually never be used */ + scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0); + +printf("scale = %f\n",scale); + /* actual scaling */ + dim->x*=scale; + dim->y*=scale; + dim->z*=scale; + if(vdim->x) vdim->x=dim->x; + if(vdim->y) vdim->y=dim->y; + if(vdim->z) vdim->z=dim->z; + moldyn->volume*=(scale*scale*scale); + + /* check whether we need a new linkcell init */ + if((dim->x/moldyn->cutoff!=lc->nx)|| + (dim->y/moldyn->cutoff!=lc->ny)|| + (dim->z/moldyn->cutoff!=lc->nx)) { + link_cell_shutdown(moldyn); + link_cell_init(moldyn); + } + + return 0; + +} + double get_e_kin(t_moldyn *moldyn) { int i; @@ -409,11 +643,6 @@ double get_e_kin(t_moldyn *moldyn) { return moldyn->ekin; } -double get_e_pot(t_moldyn *moldyn) { - - return moldyn->energy; -} - double update_e_kin(t_moldyn *moldyn) { return(get_e_kin(moldyn)); @@ -447,12 +676,6 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { /* nn_dist is the nearest neighbour distance */ - if(moldyn->t==5.0) { - printf("[moldyn] i do not estimate timesteps below %f K!\n", - MOLDYN_CRITICAL_EST_TEMP); - return 23.42; - } - tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t); return tau; @@ -468,9 +691,6 @@ int link_cell_init(t_moldyn *moldyn) { t_linkcell *lc; int i; - int fd; - - fd=open("/dev/null",O_WRONLY); lc=&(moldyn->lc); @@ -485,11 +705,13 @@ int link_cell_init(t_moldyn *moldyn) { lc->cells=lc->nx*lc->ny*lc->nz; lc->subcell=malloc(lc->cells*sizeof(t_list)); + if(lc->cells<27) + printf("[moldyn] FATAL: less then 27 subcells!\n"); + printf("[moldyn] initializing linked cells (%d)\n",lc->cells); for(i=0;icells;i++) - //list_init(&(lc->subcell[i]),1); - list_init(&(lc->subcell[i]),fd); + list_init_f(&(lc->subcell[i])); link_cell_update(moldyn); @@ -499,26 +721,30 @@ int link_cell_init(t_moldyn *moldyn) { int link_cell_update(t_moldyn *moldyn) { int count,i,j,k; - int nx,ny,nz; + int nx,ny; t_atom *atom; t_linkcell *lc; + double x,y,z; atom=moldyn->atom; lc=&(moldyn->lc); nx=lc->nx; ny=lc->ny; - nz=lc->nz; + + x=moldyn->dim.x/2; + y=moldyn->dim.y/2; + z=moldyn->dim.z/2; for(i=0;icells;i++) - list_destroy(&(moldyn->lc.subcell[i])); + list_destroy_f(&(lc->subcell[i])); for(count=0;countcount;count++) { - i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x; - j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y; - k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z; - list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]), - &(atom[count])); + i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x); + j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y); + k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z); + list_add_immediate_f(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]), + &(atom[count])); } return 0; @@ -576,9 +802,8 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } lc->dnlc=count1; - lc->countn=27; - return count2; + return count1; } int link_cell_shutdown(t_moldyn *moldyn) { @@ -589,7 +814,9 @@ int link_cell_shutdown(t_moldyn *moldyn) { lc=&(moldyn->lc); for(i=0;inx*lc->ny*lc->nz;i++) - list_shutdown(&(moldyn->lc.subcell[i])); + list_destroy_f(&(moldyn->lc.subcell[i])); + + free(lc->subcell); return 0; } @@ -601,23 +828,23 @@ int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) { t_moldyn_schedule *schedule; schedule=&(moldyn->schedule); - count=++(schedule->content_count); + count=++(schedule->total_sched); - ptr=realloc(moldyn->schedule.runs,count*sizeof(int)); + ptr=realloc(schedule->runs,count*sizeof(int)); if(!ptr) { perror("[moldyn] realloc (runs)"); return -1; } - moldyn->schedule.runs=ptr; - moldyn->schedule.runs[count-1]=runs; + schedule->runs=ptr; + schedule->runs[count-1]=runs; ptr=realloc(schedule->tau,count*sizeof(double)); if(!ptr) { perror("[moldyn] realloc (tau)"); return -1; } - moldyn->schedule.tau=ptr; - moldyn->schedule.tau[count-1]=tau; + schedule->tau=ptr; + schedule->tau[count-1]=tau; return 0; } @@ -640,16 +867,16 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { - int i,sched; + int i; unsigned int e,m,s,v; t_3dvec p; - t_moldyn_schedule *schedule; + t_moldyn_schedule *sched; t_atom *atom; int fd; - char fb[128]; + char dir[128]; double ds; - schedule=&(moldyn->schedule); + sched=&(moldyn->schedule); atom=moldyn->atom; /* initialize linked cell method */ @@ -664,10 +891,11 @@ int moldyn_integrate(t_moldyn *moldyn) { /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; + /* calculate initial forces */ potential_force_calc(moldyn); - /* do some checks before we actually start calculating bullshit */ + /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) printf("[moldyn] warning: cutoff > 0.5 x dim.x\n"); if(moldyn->cutoff>0.5*moldyn->dim.y) @@ -680,12 +908,17 @@ int moldyn_integrate(t_moldyn *moldyn) { /* zero absolute time */ moldyn->time=0.0; - for(sched=0;schedschedule.content_count;sched++) { + + /* debugging, ignore */ + moldyn->debug=0; + + /* executing the schedule */ + for(sched->count=0;sched->counttotal_sched;sched->count++) { /* setting amount of runs and finite time step size */ - moldyn->tau=schedule->tau[sched]; + moldyn->tau=sched->tau[sched->count]; moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->time_steps=schedule->runs[sched]; + moldyn->time_steps=sched->runs[sched->count]; /* integration according to schedule */ @@ -697,32 +930,38 @@ int moldyn_integrate(t_moldyn *moldyn) { /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); - - /* increase absolute time */ - moldyn->time+=moldyn->tau; + if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) + scale_volume(moldyn); /* check for log & visualization */ +//double ax; +//double ao; +//double av; if(e) { if(!(i%e)) +//ao=sqrt(0.1/M_SI); +//ax=((0.28-0.25)*sqrt(3)*LC_SI/2)*cos(ao*i); +//av=ao*(0.28-0.25)*sqrt(3)*LC_SI/2*sin(ao*i); + update_e_kin(moldyn); dprintf(moldyn->efd, - "%.15f %.45f %.45f %.45f\n", - moldyn->time,update_e_kin(moldyn), + "%f %f %f %f\n", + moldyn->time,moldyn->ekin, moldyn->energy, get_total_energy(moldyn)); +//moldyn->atom[0].r.x,ax,av*av*M_SI,0.1*ax*ax,av*av*M_SI+0.1*ax*ax); } if(m) { if(!(i%m)) { p=get_total_p(moldyn); dprintf(moldyn->mfd, - "%.15f %.45f\n",moldyn->time, - v3_norm(&p)); + "%f %f\n",moldyn->time,v3_norm(&p)); } } if(s) { if(!(i%s)) { - snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb, - moldyn->t,i*moldyn->tau); - fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT); + snprintf(dir,128,"%s/s-%07.f.save", + moldyn->vlsdir,moldyn->time); + fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT); if(fd<0) perror("[moldyn] save fd open"); else { write(fd,moldyn,sizeof(t_moldyn)); @@ -736,16 +975,23 @@ int moldyn_integrate(t_moldyn *moldyn) { if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d",sched,i); + printf("\rsched: %d, steps: %d, debug: %d", + sched->count,i,moldyn->debug); fflush(stdout); } } + /* increase absolute time */ + moldyn->time+=moldyn->tau; + } /* check for hooks */ - if(schedule->hook) - schedule->hook(moldyn,schedule->hook_params); + if(sched->hook) + sched->hook(moldyn,sched->hook_params); + + /* get a new info line */ + printf("\n"); } @@ -757,7 +1003,7 @@ int moldyn_integrate(t_moldyn *moldyn) { int velocity_verlet(t_moldyn *moldyn) { int i,count; - double tau,tau_square; + double tau,tau_square,h; t_3dvec delta; t_atom *atom; @@ -768,14 +1014,15 @@ int velocity_verlet(t_moldyn *moldyn) { for(i=0;ipotential_force_function(moldyn); for(i=0;icount; itom=moldyn->atom; @@ -823,151 +1069,155 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; + moldyn->vt2=0.0; + + /* get energy and force of every atom */ for(i=0;ixx=0.0; + virial->yy=0.0; + virial->zz=0.0; + virial->xy=0.0; + virial->xz=0.0; + virial->yz=0.0; + moldyn->vt1=0.0; + + /* reset site energy */ + itom[i].e=0.0; + /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) moldyn->func1b(moldyn,&(itom[i])); + if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) + continue; + /* 2 body pair potential/force */ - if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) { - link_cell_neighbour_index(moldyn, - (itom[i].r.x+moldyn->dim.x/2)/lc->x, - (itom[i].r.y+moldyn->dim.y/2)/lc->y, - (itom[i].r.z+moldyn->dim.z/2)/lc->z, - neighbour_i); + link_cell_neighbour_index(moldyn, + (itom[i].r.x+moldyn->dim.x/2)/lc->x, + (itom[i].r.y+moldyn->dim.y/2)/lc->y, + (itom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour_i); + + dnlc=lc->dnlc; + + for(j=0;j<27;j++) { + + this=&(neighbour_i[j]); + list_reset_f(this); - countn=lc->countn; - dnlc=lc->dnlc; + if(this->start==NULL) + continue; - for(j=0;jcurrent->data; - if(this->start==NULL) + if(jtom==&(itom[i])) continue; - bc_ij=(jattr&ATOM_ATTR_2BP)& + (itom[i].attr&ATOM_ATTR_2BP)) { + moldyn->func2b(moldyn, + &(itom[i]), + jtom, + bc_ij); + } - do { - jtom=this->current->data; + /* 3 body potential/force */ - if(jtom==&(itom[i])) - continue; + if(!(itom[i].attr&ATOM_ATTR_3BP)|| + !(jtom->attr&ATOM_ATTR_3BP)) + continue; - if((jtom->attr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + /* copy the neighbour lists */ + memcpy(neighbour_i2,neighbour_i, + 27*sizeof(t_list)); - /* 3 body potential/force */ + /* get neighbours of i */ + for(k=0;k<27;k++) { - if(!(itom[i].attr&ATOM_ATTR_3BP)|| - !(jtom->attr&ATOM_ATTR_3BP)) + that=&(neighbour_i2[k]); + list_reset_f(that); + + if(that->start==NULL) continue; - /* - * according to mr. nordlund, we dont need to take the - * sum over all atoms now, as 'this is centered' around - * atom i ... - * i am not quite sure though! there is a not vanishing - * part even if f_c_ik is zero ... - * this analytical potentials suck! - * switching from mc to md to dft soon! - */ - - // link_cell_neighbour_index(moldyn, - // (jtom->r.x+moldyn->dim.x/2)/lc->x, - // (jtom->r.y+moldyn->dim.y/2)/lc->y, - // (jtom->r.z+moldyn->dim.z/2)/lc->z, - // neighbour_j); - -// /* neighbours of j */ -// for(k=0;kcountn;k++) { -// -// that=&(neighbour_j[k]); -// list_reset(that); -// -// if(that->start==NULL) -// continue; -// -// bc_ijk=(kdnlc)?0:1; -// -// do { -// -// ktom=that->current->data; -// -// if(!(ktom->attr&ATOM_ATTR_3BP)) -// continue; -// -// if(ktom==jtom) -// continue; -// -// if(ktom==&(itom[i])) -// continue; -// -// moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); -// -/* } while(list_next(that)!=\ */ -// L_NO_NEXT_ELEMENT); -// -// } - - /* copy the neighbour lists */ - memcpy(neighbour_i2,neighbour_i, - 27*sizeof(t_list)); - - /* get neighbours of i */ - for(k=0;kstart==NULL) - continue; + bc_ik=(kcurrent->data; - ktom=that->current->data; + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; - if(!(ktom->attr&ATOM_ATTR_3BP)) - continue; + if(ktom==jtom) + continue; - if(ktom==jtom) - continue; + if(ktom==&(itom[i])) + continue; - if(ktom==&(itom[i])) - continue; + moldyn->func3b(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); -//printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); - moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); -//printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); - } while(list_next(that)!=\ - L_NO_NEXT_ELEMENT); + } - } - - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - /* 2bp post function */ - if(moldyn->func2b_post) + if(moldyn->func2b_post) { moldyn->func2b_post(moldyn, &(itom[i]), jtom,bc_ij); - - } + } + + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); + } + } +#ifdef DEBUG +printf("\n\n"); +#endif +#ifdef VDEBUG +printf("\n\n"); +#endif + + moldyn->vt2=0.0; + for(i=0;ivt2-=v3_scalar_product(&(itom[i].r),&(itom[i].f)); + +printf("compare: vt1: %f vt2: %f\n",moldyn->vt1,moldyn->vt2); + +pressure_calc(moldyn); + + return 0; +} + +/* + * virial calculation + */ + +inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { + + a->virial.xx-=f->x*d->x; + a->virial.yy-=f->y*d->y; + a->virial.zz-=f->z*d->z; + a->virial.xy-=f->x*d->y; + a->virial.xz-=f->x*d->z; + a->virial.yz-=f->y*d->z; return 0; } @@ -976,16 +1226,16 @@ int potential_force_calc(t_moldyn *moldyn) { * periodic boundayr checking */ -int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { +inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { double x,y,z; t_3dvec *dim; dim=&(moldyn->dim); - x=0.5*dim->x; - y=0.5*dim->y; - z=0.5*dim->z; + x=dim->x/2; + y=dim->y/2; + z=dim->z/2; if(moldyn->status&MOLDYN_STAT_PBX) { if(a->x>=x) a->x-=dim->x; @@ -1014,22 +1264,29 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_ho_params *params; t_3dvec force,distance; - double d; + double d,f; double sc,equi_dist; params=moldyn->pot2b_params; sc=params->spring_constant; equi_dist=params->equilibrium_distance; - v3_sub(&distance,&(ai->r),&(aj->r)); + if(air),&(ai->r)); if(bc) check_per_bound(moldyn,&distance); d=v3_norm(&distance); if(d<=moldyn->cutoff) { - /* energy is 1/2 (d-d0)^2, but we will add this twice ... */ - moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist)); - v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d))); + moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); + /* f = -grad E; grad r_ij = -1 1/r_ij distance */ + f=sc*(1.0-equi_dist/d); + v3_scale(&force,&distance,f); v3_add(&(ai->f),&(ai->f),&force); + virial_calc(ai,&force,&distance); + virial_calc(aj,&force,&distance); /* f and d signe switched */ + v3_scale(&force,&distance,-f); + v3_add(&(aj->f),&(aj->f),&force); } return 0; @@ -1049,7 +1306,9 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { sig6=params->sigma6; sig12=params->sigma12; - v3_sub(&distance,&(ai->r),&(aj->r)); + if(air),&(ai->r)); if(bc) check_per_bound(moldyn,&distance); d=v3_absolute_square(&distance); /* 1/r^2 */ if(d<=moldyn->cutoff_square) { @@ -1057,8 +1316,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { h2=d*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ - /* energy is eps*..., but we will add this twice ... */ - moldyn->energy+=0.5*eps*(sig12*h1-sig6*h2); + moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); h2*=d; /* 1/r^8 */ h1*=d; /* 1/r^14 */ h2*=6*sig6; @@ -1066,7 +1324,12 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { d=+h1-h2; d*=eps; v3_scale(&force,&distance,d); + v3_add(&(aj->f),&(aj->f),&force); + v3_scale(&force,&distance,-1.0*d); /* f = - grad E */ v3_add(&(ai->f),&(ai->f),&force); + virial_calc(ai,&force,&distance); + virial_calc(aj,&force,&distance); /* f and d signe switched */ + moldyn->vt1-=v3_scalar_product(&force,&distance); } return 0; @@ -1080,7 +1343,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int tersoff_mult_complete_params(t_tersoff_mult_params *p) { printf("[moldyn] tersoff parameter completion\n"); + p->S2[0]=p->S[0]*p->S[0]; + p->S2[1]=p->S[1]*p->S[1]; p->Smixed=sqrt(p->S[0]*p->S[1]); + p->S2mixed=p->Smixed*p->Smixed; p->Rmixed=sqrt(p->R[0]*p->R[1]); p->Amixed=sqrt(p->A[0]*p->A[1]); p->Bmixed=sqrt(p->B[0]*p->B[1]); @@ -1088,8 +1354,8 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { p->mu_m=0.5*(p->mu[0]+p->mu[1]); printf("[moldyn] tersoff mult parameter info:\n"); - printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed); - printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed); + printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed); + printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed); printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV); printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], @@ -1108,11 +1374,11 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { /* tersoff 1 body part */ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { - int num; + int brand; t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - num=ai->bnum; + brand=ai->brand; params=moldyn->pot1b_params; exchange=&(params->exchange); @@ -1121,17 +1387,16 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { * their right values */ - exchange->beta=&(params->beta[num]); - exchange->n=&(params->n[num]); - exchange->c=&(params->c[num]); - exchange->d=&(params->d[num]); - exchange->h=&(params->h[num]); + exchange->beta_i=&(params->beta[brand]); + exchange->n_i=&(params->n[brand]); + exchange->c_i=&(params->c[brand]); + exchange->d_i=&(params->d[brand]); + exchange->h_i=&(params->h[brand]); - exchange->betan=pow(*(exchange->beta),*(exchange->n)); - exchange->n_betan=*(exchange->n)*exchange->betan; - exchange->c2=params->c[num]*params->c[num]; - exchange->d2=params->d[num]*params->d[num]; - exchange->c2d2=exchange->c2/exchange->d2; + exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i)); + exchange->ci2=params->c[brand]*params->c[brand]; + exchange->di2=params->d[brand]*params->d[brand]; + exchange->ci2di2=exchange->ci2/exchange->di2; return 0; } @@ -1142,50 +1407,52 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; t_3dvec dist_ij,force; - double d_ij; - double A,B,R,S,lambda,mu; + double d_ij,d_ij2; + double A,B,R,S,S2,lambda,mu; double f_r,df_r; double f_c,df_c; - int num; + int brand; double s_r; double arg; - double scale; params=moldyn->pot2b_params; - num=ai->bnum; + brand=aj->brand; exchange=&(params->exchange); + /* clear 3bp and 2bp post run */ exchange->run3bp=0; + exchange->run2bp_post=0; + + /* reset S > r > R mark */ + exchange->d_ij_between_rs=0; /* - * we need: f_c, df_c, f_r, df_r + * calc of 2bp contribution of V_ij and dV_ij/ji + * + * for Vij and dV_ij we need: + * - f_c_ij, df_c_ij + * - f_r_ij, df_r_ij + * + * for dV_ji we need: + * - f_c_ji = f_c_ij, df_c_ji = df_c_ij + * - f_r_ji = f_r_ij; df_r_ji = df_r_ij * - * therefore we need: R, S, A, lambda */ - v3_sub(&dist_ij,&(ai->r),&(aj->r)); - - if(bc) check_per_bound(moldyn,&dist_ij); - - d_ij=v3_norm(&dist_ij); - - /* save for use in 3bp */ - exchange->d_ij=d_ij; - exchange->dist_ij=dist_ij; - exchange->d_ij2=d_ij*d_ij; - /* constants */ - if(num==aj->bnum) { - S=params->S[num]; - R=params->R[num]; - A=params->A[num]; - B=params->B[num]; - lambda=params->lambda[num]; - mu=params->mu[num]; - params->exchange.chi=1.0; + if(brand==ai->brand) { + S=params->S[brand]; + S2=params->S2[brand]; + R=params->R[brand]; + A=params->A[brand]; + B=params->B[brand]; + lambda=params->lambda[brand]; + mu=params->mu[brand]; + exchange->chi=1.0; } else { S=params->Smixed; + S2=params->S2mixed; R=params->Rmixed; A=params->Amixed; B=params->Bmixed; @@ -1194,47 +1461,115 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { params->exchange.chi=params->chi; } - if(d_ij>S) + /* dist_ij, d_ij */ + v3_sub(&dist_ij,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist_ij); + d_ij2=v3_absolute_square(&dist_ij); + + /* if d_ij2 > S2 => no force & potential energy contribution */ + if(d_ij2>S2) return 0; + /* now we will need the distance */ + //d_ij=v3_norm(&dist_ij); + d_ij=sqrt(d_ij2); + + /* save for use in 3bp */ + exchange->d_ij=d_ij; + exchange->d_ij2=d_ij2; + exchange->dist_ij=dist_ij; + + /* more constants */ + exchange->beta_j=&(params->beta[brand]); + exchange->n_j=&(params->n[brand]); + exchange->c_j=&(params->c[brand]); + exchange->d_j=&(params->d[brand]); + exchange->h_j=&(params->h[brand]); + if(brand==ai->brand) { + exchange->betajnj=exchange->betaini; + exchange->cj2=exchange->ci2; + exchange->dj2=exchange->di2; + exchange->cj2dj2=exchange->ci2di2; + } + else { + exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j)); + exchange->cj2=params->c[brand]*params->c[brand]; + exchange->dj2=params->d[brand]*params->d[brand]; + exchange->cj2dj2=exchange->cj2/exchange->dj2; + } + + /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */ f_r=A*exp(-lambda*d_ij); - df_r=-lambda*f_r/d_ij; + df_r=lambda*f_r/d_ij; - /* f_a, df_a calc + save for 3bp use */ + /* f_a, df_a calc (again, same for ij and ji) | save for later use! */ exchange->f_a=-B*exp(-mu*d_ij); - exchange->df_a=-mu*exchange->f_a/d_ij; + exchange->df_a=mu*exchange->f_a/d_ij; + /* f_c, df_c calc (again, same for ij and ji) */ if(d_ij r > R */ + exchange->d_ij_between_rs=1; } - /* add forces */ + /* add forces of 2bp (ij, ji) contribution + * dVij = dVji and we sum up both: no 1/2) */ v3_add(&(ai->f),&(ai->f),&force); - /* energy is 0.5 f_r f_c ... */ + + /* virial */ + ai->virial.xx-=force.x*dist_ij.x; + ai->virial.yy-=force.y*dist_ij.y; + ai->virial.zz-=force.z*dist_ij.z; + ai->virial.xy-=force.x*dist_ij.y; + ai->virial.xz-=force.x*dist_ij.z; + ai->virial.yz-=force.y*dist_ij.z; + +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij, dVji (2bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij, dVji (2bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz); +} +#endif + + /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */ moldyn->energy+=(0.5*f_r*f_c); /* save for use in 3bp */ exchange->f_c=f_c; exchange->df_c=df_c; - /* enable the run of 3bp function */ + /* enable the run of 3bp function and 2bp post processing */ exchange->run3bp=1; + exchange->run2bp_post=1; /* reset 3bp sums */ - exchange->sum1_3bp=0.0; - exchange->sum2_3bp=0.0; - v3_zero(&(exchange->db_ij)); + exchange->zeta_ij=0.0; + exchange->zeta_ji=0.0; + v3_zero(&(exchange->dzeta_ij)); + v3_zero(&(exchange->dzeta_ji)); return 0; } @@ -1243,51 +1578,150 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - /* here we have to allow for the 3bp sums */ + /* + * here we have to allow for the 3bp sums + * + * that is: + * - zeta_ij, dzeta_ij + * - zeta_ji, dzeta_ji + * + * to compute the 3bp contribution to: + * - Vij, dVij + * - dVji + * + */ t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - t_3dvec force,temp,*db_ij,*dist_ij; - double db_ij_scale1,db_ij_scale2; - double b_ij; + t_3dvec force,temp; + t_3dvec *dist_ij; + double b,db,tmp; double f_c,df_c,f_a,df_a; - double chi,betan; - double help; - double n; + double chi,ni,betaini,nj,betajnj; + double zeta; params=moldyn->pot2b_params; exchange=&(params->exchange); - db_ij=&(exchange->db_ij); + /* we do not run if f_c_ij was detected to be 0! */ + if(!(exchange->run2bp_post)) + return 0; + f_c=exchange->f_c; df_c=exchange->df_c; f_a=exchange->f_a; df_a=exchange->df_a; - betan=exchange->betan; - n=*(exchange->n); + betaini=exchange->betaini; + betajnj=exchange->betajnj; + ni=*(exchange->n_i); + nj=*(exchange->n_j); + chi=exchange->chi; dist_ij=&(exchange->dist_ij); + + /* Vij and dVij */ + zeta=exchange->zeta_ij; + if(zeta==0.0) { + moldyn->debug++; /* just for debugging ... */ + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ + db=chi*pow(b,-1.0/(2*ni)-1); /* x(...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ij),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,f_c); + } + v3_scale(&temp,dist_ij,df_c*b*f_a); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,-0.5); - db_ij_scale1=(1+betan*exchange->sum1_3bp); - db_ij_scale2=(exchange->n_betan*exchange->sum2_3bp); - help=pow(db_ij_scale1,-1.0/(2*n)-1); - b_ij=chi*db_ij_scale1*help; - db_ij_scale1=-chi/(2*n)*help; - - v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2)); - v3_scale(db_ij,db_ij,f_a); + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); - v3_scale(&temp,dist_ij,b_ij*df_a); + /* virial */ + ai->virial.xx-=force.x*dist_ij->x; + ai->virial.yy-=force.y*dist_ij->y; + ai->virial.zz-=force.z*dist_ij->z; + ai->virial.xy-=force.x*dist_ij->y; + ai->virial.xz-=force.x*dist_ij->z; + ai->virial.yz-=force.y*dist_ij->z; + +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij (3bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij (3bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); +} +#endif - v3_add(&force,&temp,db_ij); - v3_scale(&force,&force,f_c); + /* add energy of 3bp sum */ + moldyn->energy+=(0.5*f_c*b*f_a); + + /* dVji */ + zeta=exchange->zeta_ji; + if(zeta==0.0) { + moldyn->debug++; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ + db=chi*pow(b,-1.0/(2*nj)-1); /* x(...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ji),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,f_c); + } + v3_scale(&temp,dist_ij,df_c*b*f_a); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,-0.5); - v3_scale(&temp,dist_ij,f_a*b_ij*df_c); + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); - /* add energy of 3bp sum */ - moldyn->energy+=(0.5*f_c*b_ij*f_a); - /* add force of 3bp calculation */ - v3_add(&(ai->f),&temp,&force); + /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ +// TEST ... with a minus instead + ai->virial.xx-=force.x*dist_ij->x; + ai->virial.yy-=force.y*dist_ij->y; + ai->virial.zz-=force.z*dist_ij->z; + ai->virial.xy-=force.x*dist_ij->y; + ai->virial.xz-=force.x*dist_ij->z; + ai->virial.yz-=force.y*dist_ij->z; + +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVji (3bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVji (3bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); +} +#endif return 0; } @@ -1299,151 +1733,347 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; t_3dvec dist_ij,dist_ik,dist_jk; - t_3dvec temp,force; - double R,S,s_r; - double d_ij,d_ij2,d_ik,d_jk; - double f_c,df_c,f_a,df_a; + t_3dvec temp1,temp2; + t_3dvec *dzeta; + double R,S,S2,s_r; + double B,mu; + double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2; + double rr,dd; + double f_c,df_c; double f_c_ik,df_c_ik,arg; + double f_c_jk; double n,c,d,h; double c2,d2,c2d2; - double numer,denom; - double theta,cos_theta,sin_theta; - double d_theta,d_theta1,d_theta2; + double cos_theta,d_costheta1,d_costheta2; double h_cos,d2_h_cos2; - double frac,bracket,bracket_n_1,bracket_n; - double g; - int num; + double frac,g,zeta,chi; + double tmp; + int brand; params=moldyn->pot3b_params; - num=ai->bnum; exchange=&(params->exchange); if(!(exchange->run3bp)) return 0; /* - * we need: f_c, d_fc, b_ij, db_ij, f_a, df_a + * calc of 3bp contribution of V_ij and dV_ij/ji/jk & + * 2bp contribution of dV_jk + * + * for Vij and dV_ij we still need: + * - b_ij, db_ij (zeta_ij) + * - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk * - * we got f_c, df_c, f_a, df_a from 2bp calculation + * for dV_ji we still need: + * - b_ji, db_ji (zeta_ji) + * - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik + * + * for dV_jk we need: + * - f_c_jk + * - f_a_jk + * - db_jk (zeta_jk) + * - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki + * + */ + + /* + * get exchange data */ + /* dist_ij, d_ij - this is < S_ij ! */ + dist_ij=exchange->dist_ij; d_ij=exchange->d_ij; d_ij2=exchange->d_ij2; - dist_ij=exchange->dist_ij; - - f_a=params->exchange.f_a; - df_a=params->exchange.df_a; + /* f_c_ij, df_c_ij (same for ji) */ f_c=exchange->f_c; df_c=exchange->df_c; - - /* d_ij is <= S, as we didn't return so far! */ /* - * calc of b_ij (scalar) and db_ij (vector) - * - * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta - * - * - for db_ij: d_theta, sin_theta, cos_theta, f_c_ik, df_c_ik, - * w_ik, - * + * calculate unknown values now ... */ - v3_sub(&dist_ik,&(ai->r),&(ak->r)); - if(bc) check_per_bound(moldyn,&dist_ik); - d_ik=v3_norm(&dist_ik); + /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */ - /* constants for f_c_ik calc */ - if(num==ak->bnum) { - R=params->R[num]; - S=params->S[num]; + /* dist_ik, d_ik */ + v3_sub(&dist_ik,&(ak->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist_ik); + d_ik2=v3_absolute_square(&dist_ik); + + /* ik constants */ + brand=ai->brand; + if(brand==ak->brand) { + R=params->R[brand]; + S=params->S[brand]; + S2=params->S2[brand]; } else { R=params->Rmixed; S=params->Smixed; + S2=params->S2mixed; } - /* calc of f_c_ik */ - if(d_ik>S) { - f_c_ik=0.0; - df_c_ik=0.0; - } - else if(d_ikn_i); + c=*(exchange->c_i); + d=*(exchange->d_i); + h=*(exchange->h_i); + c2=exchange->ci2; + d2=exchange->di2; + c2d2=exchange->ci2di2; + + /* cosine of theta_ijk by scalaproduct */ + rr=v3_scalar_product(&dist_ij,&dist_ik); + dd=d_ij*d_ik; + cos_theta=rr/dd; + + /* d_costheta */ + tmp=1.0/dd; + d_costheta1=cos_theta/d_ij2-tmp; + d_costheta2=cos_theta/d_ik2-tmp; + + /* some usefull values */ + h_cos=(h-cos_theta); + d2_h_cos2=d2+(h_cos*h_cos); + frac=c2/(d2_h_cos2); + + /* g(cos_theta) */ + g=1.0+c2d2-frac; + + /* d_costheta_ij and dg(cos_theta) - needed in any case! */ + v3_scale(&temp1,&dist_ij,d_costheta1); + v3_scale(&temp2,&dist_ik,d_costheta2); + v3_add(&temp1,&temp1,&temp2); + v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ + + /* f_c_ik & df_c_ik + {d,}zeta contribution */ + dzeta=&(exchange->dzeta_ij); + if(d_ik f_c_ik=1.0; + // => df_c_ik=0.0; of course we do not set this! + + /* zeta_ij */ + exchange->zeta_ij+=g; + + /* dzeta_ij */ + v3_add(dzeta,dzeta,&temp1); + } + else { + /* {d,}f_c_ik */ + 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)); + + /* zeta_ij */ + exchange->zeta_ij+=f_c_ik*g; + + /* dzeta_ij */ + v3_scale(&temp1,&temp1,f_c_ik); + v3_scale(&temp2,&dist_ik,g*df_c_ik); + v3_add(&temp1,&temp1,&temp2); + v3_add(dzeta,dzeta,&temp1); + } } - - v3_sub(&dist_jk,&(aj->r),&(ak->r)); + + /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */ + + /* dist_jk, d_jk */ + v3_sub(&dist_jk,&(ak->r),&(aj->r)); if(bc) check_per_bound(moldyn,&dist_jk); - d_jk=v3_norm(&dist_jk); - - /* get exchange data */ - n=*(exchange->n); - c=*(exchange->c); - d=*(exchange->d); - h=*(exchange->h); - c2=exchange->c2; - d2=exchange->d2; - c2d2=exchange->c2d2; - - numer=d_ij2+d_ik*d_ik-d_jk*d_jk; - denom=2*d_ij*d_ik; - cos_theta=numer/denom; - /* prefere law of cosines, dot product -> nan (often) */ - //cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - sin_theta=sqrt(1.0-(cos_theta*cos_theta)); - theta=acos(cos_theta); - d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom); - d_theta1=2*denom-numer*2*d_ik/d_ij; - d_theta2=2*denom-numer*2*d_ij/d_ik; - d_theta1*=d_theta; - d_theta2*=d_theta; - - h_cos=(h-cos_theta); - d2_h_cos2=d2+(h_cos*h_cos); - - frac=c2/(d2_h_cos2); - g=1.0+c2d2-frac; - - if(f_c_ik==0.0) { - bracket=0.0; - bracket_n_1=0.0; - bracket_n=0.0; - printf("Foo -> 0: "); + d_jk2=v3_absolute_square(&dist_jk); + + /* jk constants */ + brand=aj->brand; + if(brand==ak->brand) { + R=params->R[brand]; + S=params->S[brand]; + S2=params->S2[brand]; + B=params->B[brand]; + mu=params->mu[brand]; + chi=1.0; } else { - bracket=f_c_ik*g; - bracket_n_1=pow(bracket,n-1.0); - bracket_n=bracket_n_1*bracket; - printf("Foo -> 1: "); + R=params->Rmixed; + S=params->Smixed; + S2=params->S2mixed; + B=params->Bmixed; + mu=params->mu_m; + chi=params->chi; } -//printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket); - /* calc of db_ij and the 2 sums */ - exchange->sum1_3bp+=bracket_n; - exchange->sum2_3bp+=bracket_n_1; - - /* derivation of theta */ - v3_scale(&force,&dist_ij,d_theta1); - v3_scale(&temp,&dist_ik,d_theta2); - v3_add(&force,&force,&temp); + /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */ + if(d_jk2n_j); + c=*(exchange->c_j); + d=*(exchange->d_j); + h=*(exchange->h_j); + c2=exchange->cj2; + d2=exchange->dj2; + c2d2=exchange->cj2dj2; + + /* cosine of theta_jik by scalaproduct */ + rr=-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji */ + dd=d_ij*d_jk; + cos_theta=rr/dd; + + /* d_costheta */ + d_costheta1=1.0/dd; + d_costheta2=cos_theta/d_ij2; + + /* some usefull values */ + h_cos=(h-cos_theta); + d2_h_cos2=d2+(h_cos*h_cos); + frac=c2/(d2_h_cos2); + + /* g(cos_theta) */ + g=1.0+c2d2-frac; + + /* d_costheta_jik and dg(cos_theta) - needed in any case! */ + v3_scale(&temp1,&dist_jk,d_costheta1); + v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */ + //v3_add(&temp1,&temp1,&temp2); + v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */ + v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ + + /* store dg in temp2 and use it for dVjk later */ + v3_copy(&temp2,&temp1); + + /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ + dzeta=&(exchange->dzeta_ji); + if(d_jkzeta_ji+=g; + + /* dzeta_ji */ + v3_add(dzeta,dzeta,&temp1); + } + else { + /* f_c_jk */ + s_r=S-R; + arg=M_PI*(d_jk-R)/s_r; + f_c_jk=0.5+0.5*cos(arg); + + /* zeta_ji */ + exchange->zeta_ji+=f_c_jk*g; + + /* dzeta_ji */ + v3_scale(&temp1,&temp1,f_c_jk); + v3_add(dzeta,dzeta,&temp1); + } - /* part 1 of db_ij */ - v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2); + /* dV_jk stuff | add force contribution on atom i immediately */ + if(exchange->d_ij_between_rs) { + zeta=f_c*g; + v3_scale(&temp1,&temp2,f_c); + v3_scale(&temp2,&dist_ij,df_c*g); + v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */ + } + else { + zeta=g; + // dzeta_jk is simply dg, which is stored in temp2 + } + /* betajnj * zeta_jk ^ nj-1 */ + tmp=exchange->betajnj*pow(zeta,(n-1.0)); + tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp; + v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); + v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ + /* scaled with 0.5 ^ */ + + /* virial */ + ai->virial.xx-=temp2.x*dist_jk.x; + ai->virial.yy-=temp2.y*dist_jk.y; + ai->virial.zz-=temp2.z*dist_jk.z; + ai->virial.xy-=temp2.x*dist_jk.y; + ai->virial.xz-=temp2.x*dist_jk.z; + ai->virial.yz-=temp2.y*dist_jk.z; + +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVjk (3bp) contrib:\n"); + printf("%f | %f\n",temp2.x,ai->f.x); + printf("%f | %f\n",temp2.y,ai->f.y); + printf("%f | %f\n",temp2.z,ai->f.z); +} +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVjk (3bp) contrib:\n"); + printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx); + printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy); + printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz); +} +#endif - /* part 2 of db_ij */ - v3_scale(&temp,&dist_ik,df_c_ik*g); + } - /* sum up and add to db_ij */ - v3_add(&temp,&temp,&force); - v3_add(&(exchange->db_ij),&(exchange->db_ij),&temp); - return 0; } + +/* + * debugging / critical check functions + */ + +int moldyn_bc_check(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim; + int i; + double x; + u8 byte; + int j,k; + + atom=moldyn->atom; + dim=&(moldyn->dim); + x=dim->x/2; + + for(i=0;icount;i++) { + if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) { + printf("FATAL: atom %d: x: %.20f (%.20f)\n", + i,atom[i].r.x,dim->x/2); + printf("diagnostic:\n"); + printf("-----------\natom.r.x:\n"); + for(j=0;j<8;j++) { + memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1); + for(k=0;k<8;k++) + printf("%d%c", + ((byte)&(1<=dim->y/2||-atom[i].r.y>dim->y/2) + printf("FATAL: atom %d: y: %.20f (%.20f)\n", + i,atom[i].r.y,dim->y/2); + if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2) + printf("FATAL: atom %d: z: %.20f (%.20f)\n", + i,atom[i].r.z,dim->z/2); + } + + return 0; +}