From: hackbard Date: Fri, 27 Apr 2007 16:39:13 +0000 (+0000) Subject: add another way of calculating tersoff + small moldyn mods X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=0b96eb313c9bfec6272b1f8de0d99c4ce26d1686;p=physik%2Fposic.git add another way of calculating tersoff + small moldyn mods --- diff --git a/Makefile b/Makefile index 3e31459..70fcb8b 100644 --- a/Makefile +++ b/Makefile @@ -8,10 +8,11 @@ CFLAGS+=-g #CFLAGS+=-DVDEBUG LDFLAGS=-lm -SOURCE=moldyn.c visual/visual.c random/random.c +SOURCE=moldyn.c visual/visual.c random/random.c list/list.c -POT_SRC=potentials/tersoff.c potentials/lennard_jones.c -POT_SRC+= potentials/harmonic_oscillator.c +POT_SRC=potentials/tersoff.c +#POT_SRC=potentials/tersoff_orig.c +POT_SRC+= potentials/lennard_jones.c potentials/harmonic_oscillator.c all: sic diff --git a/moldyn.c b/moldyn.c index 9a62b9f..4d73617 100644 --- a/moldyn.c +++ b/moldyn.c @@ -82,7 +82,7 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { moldyn->p_ref=p_ref; - printf("[moldyn] pressure [atm]: %f\n",moldyn->p_ref/ATM); + printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR); return 0; } @@ -176,16 +176,37 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func) { return 0; } -int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func) { +int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { - moldyn->func2b_post=func; + moldyn->func3b_j1=func; return 0; } -int set_potential3b(t_moldyn *moldyn,pf_func3b func) { +int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) { - moldyn->func3b=func; + moldyn->func3b_j2=func; + + return 0; +} + +int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) { + + moldyn->func3b_j3=func; + + return 0; +} + +int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) { + + moldyn->func3b_k1=func; + + return 0; +} + +int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) { + + moldyn->func3b_k2=func; return 0; } @@ -294,6 +315,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { perror("[moldyn] report fd open"); return moldyn->rfd; } + printf("report -> "); if(moldyn->efd) { snprintf(filename,127,"%s/e_plot.scr", moldyn->vlsdir); @@ -306,6 +328,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { } dprintf(moldyn->epfd,e_plot_script); close(moldyn->epfd); + printf("energy "); } if(moldyn->pfd) { snprintf(filename,127,"%s/pressure_plot.scr", @@ -319,6 +342,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { } dprintf(moldyn->ppfd,pressure_plot_script); close(moldyn->ppfd); + printf("pressure "); } if(moldyn->tfd) { snprintf(filename,127,"%s/temperature_plot.scr", @@ -332,9 +356,11 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { } dprintf(moldyn->tpfd,temperature_plot_script); close(moldyn->tpfd); + printf("temperature "); } dprintf(moldyn->rfd,report_start, moldyn->rauthor,moldyn->rtitle); + printf("\n"); break; default: printf("unknown log type: %02x\n",type); @@ -672,6 +698,8 @@ double temperature_calc(t_moldyn *moldyn) { /* 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); + moldyn->t_sum+=moldyn->t; + moldyn->mean_t=moldyn->t_sum/moldyn->total_steps; return moldyn->t; } @@ -769,6 +797,16 @@ double pressure_calc(t_moldyn *moldyn) { /* assume up to date kinetic energy */ moldyn->p=2.0*moldyn->ekin+v; moldyn->p/=(3.0*moldyn->volume); + moldyn->p_sum+=moldyn->p; + moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; + + /* pressure from 'absolute coordinates' virial */ + virial=&(moldyn->virial); + v=virial->xx+virial->yy+virial->zz; + moldyn->gp=2.0*moldyn->ekin+v; + moldyn->gp/=(3.0*moldyn->volume); + moldyn->gp_sum+=moldyn->gp; + moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps; return moldyn->p; } @@ -788,7 +826,7 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * dV: dx,y,z = 0.001 x,y,z */ - scale=1.00001; + scale=1.00000000000001; printf("\n\nP-DEBUG:\n"); tp=&(moldyn->tp); @@ -806,13 +844,12 @@ printf("\n\nP-DEBUG:\n"); /* derivative with respect to x direction */ scale_dim(moldyn,scale,TRUE,0,0); scale_atoms(moldyn,scale,TRUE,0,0); - dv=0.00001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z; + dv=0.00000000000001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); tp->x=(moldyn->energy-u)/dv; p=tp->x*tp->x; -printf("e: %f eV de: %f eV dV: %f A^3\n",moldyn->energy/moldyn->count/EV,(moldyn->energy-u)/moldyn->count/EV,dv); /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); @@ -821,7 +858,7 @@ printf("e: %f eV de: %f eV dV: %f A^3\n",moldyn->energy/moldyn->count/EV,(moldyn /* derivative with respect to y direction */ scale_dim(moldyn,scale,0,TRUE,0); scale_atoms(moldyn,scale,0,TRUE,0); - dv=0.00001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z; + dv=0.00000000000001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); @@ -835,7 +872,7 @@ printf("e: %f eV de: %f eV dV: %f A^3\n",moldyn->energy/moldyn->count/EV,(moldyn /* derivative with respect to z direction */ scale_dim(moldyn,scale,0,0,TRUE); scale_atoms(moldyn,scale,0,0,TRUE); - dv=0.00001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y; + dv=0.00000000000001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); @@ -939,7 +976,7 @@ moldyn->debug=scale; } -double get_e_kin(t_moldyn *moldyn) { +double e_kin_calc(t_moldyn *moldyn) { int i; t_atom *atom; @@ -953,11 +990,6 @@ double get_e_kin(t_moldyn *moldyn) { return moldyn->ekin; } -double update_e_kin(t_moldyn *moldyn) { - - return(get_e_kin(moldyn)); -} - double get_total_energy(t_moldyn *moldyn) { return(moldyn->ekin+moldyn->energy); @@ -1018,7 +1050,12 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { if(lc->cells<27) printf("[moldyn] FATAL: less then 27 subcells!\n"); - if(vol) printf("[moldyn] initializing linked cells (%d)\n",lc->cells); + if(vol) { + printf("[moldyn] initializing linked cells (%d)\n",lc->cells); + printf(" x: %d x %f A\n",lc->nx,lc->x); + printf(" y: %d x %f A\n",lc->ny,lc->y); + printf(" z: %d x %f A\n",lc->nz,lc->z); + } for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); @@ -1053,7 +1090,7 @@ int link_cell_update(t_moldyn *moldyn) { 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]), + list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]), &(atom[count])); } @@ -1229,6 +1266,7 @@ int moldyn_integrate(t_moldyn *moldyn) { /* zero absolute time */ moldyn->time=0.0; + moldyn->total_steps=0; /* debugging, ignore */ moldyn->debug=0; @@ -1252,10 +1290,11 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->integrate(moldyn); /* calculate kinetic energy, temperature and pressure */ - update_e_kin(moldyn); + e_kin_calc(moldyn); temperature_calc(moldyn); pressure_calc(moldyn); //tp=thermodynamic_pressure_calc(moldyn); +//printf("thermodynamic p: %f %f %f - %f\n",moldyn->tp.x/BAR,moldyn->tp.y/BAR,moldyn->tp.z/BAR,tp/BAR); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1284,13 +1323,16 @@ int moldyn_integrate(t_moldyn *moldyn) { if(p) { if(!(i%p)) { dprintf(moldyn->pfd, - "%f %f\n",moldyn->time,moldyn->p/ATM); + "%f %f %f %f %f\n",moldyn->time, + moldyn->p/BAR,moldyn->mean_p/BAR, + moldyn->gp/BAR,moldyn->mean_gp/BAR); } } if(t) { if(!(i%t)) { dprintf(moldyn->tfd, - "%f %f\n",moldyn->time,moldyn->t); + "%f %f %f\n", + moldyn->time,moldyn->t,moldyn->mean_t); } } if(s) { @@ -1316,14 +1358,18 @@ int moldyn_integrate(t_moldyn *moldyn) { /* display progress */ if(!(i%10)) { - printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f", + printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f", sched->count,i, - moldyn->t,moldyn->p/ATM,moldyn->volume); + moldyn->mean_t, + moldyn->mean_p/BAR, + moldyn->mean_gp/BAR, + moldyn->volume); fflush(stdout); } /* increase absolute time */ moldyn->time+=moldyn->tau; + moldyn->total_steps+=1; } @@ -1410,6 +1456,9 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; + /* reset global virial */ + memset(&(moldyn->virial),0,sizeof(t_virial)); + /* reset force, site energy and virial of every atom */ for(i=0;idnlc; + /* first loop over atoms j */ + if(moldyn->func2b) { + for(j=0;j<27;j++) { + + this=&(neighbour_i[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + bc_ij=(jcurrent->data; + + if(jtom==&(itom[i])) + continue; + + if((jtom->attr&ATOM_ATTR_2BP)& + (itom[i].attr&ATOM_ATTR_2BP)) { + moldyn->func2b(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); + + } + } + + /* 3 body potential/force */ + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + /* copy the neighbour lists */ + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); + + /* second loop over atoms j */ for(j=0;j<27;j++) { this=&(neighbour_i[j]); @@ -1466,25 +1556,70 @@ int potential_force_calc(t_moldyn *moldyn) { if(jtom==&(itom[i])) continue; - if((jtom->attr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; - /* 3 body potential/force */ + /* reset 3bp run */ + moldyn->run3bp=1; - if(!(itom[i].attr&ATOM_ATTR_3BP)|| - !(jtom->attr&ATOM_ATTR_3BP)) + if(moldyn->func3b_j1) + moldyn->func3b_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); + + /* in first j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) continue; + + /* first loop over atoms k */ + if(moldyn->func3b_k1) { + + for(k=0;k<27;k++) { + + that=&(neighbour_i2[k]); + list_reset_f(that); + + if(that->start==NULL) + continue; + + bc_ik=(kcurrent->data; + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + if(ktom==jtom) + continue; + + if(ktom==&(itom[i])) + continue; + + moldyn->func3b_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); + + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); + + } + + } + + if(moldyn->func3b_j2) + moldyn->func3b_j2(moldyn, + &(itom[i]), + jtom, + bc_ij); + + /* second loop over atoms k */ + if(moldyn->func3b_k2) { - /* get neighbours of i */ for(k=0;k<27;k++) { that=&(neighbour_i2[k]); @@ -1508,37 +1643,53 @@ int potential_force_calc(t_moldyn *moldyn) { if(ktom==&(itom[i])) continue; - moldyn->func3b(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func3b_k2(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); } + + } /* 2bp post function */ - if(moldyn->func2b_post) { - moldyn->func2b_post(moldyn, - &(itom[i]), - jtom,bc_ij); + if(moldyn->func3b_j3) { + moldyn->func3b_j3(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 } #ifdef DEBUG -printf("\n\n"); -#endif -#ifdef VDEBUG -printf("\n\n"); + printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z); #endif + /* calculate global virial */ + for(i=0;ivirial.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; + moldyn->virial.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; + moldyn->virial.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; + moldyn->virial.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; + moldyn->virial.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; + moldyn->virial.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; + } + return 0; } diff --git a/moldyn.h b/moldyn.h index 85a4950..5b12048 100644 --- a/moldyn.h +++ b/moldyn.h @@ -85,11 +85,15 @@ typedef struct s_moldyn { /* 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 (*func2b_post)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak, - u8 bck); + 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); void *pot_params; - //int (*potential_force_function)(struct s_moldyn *moldyn); + unsigned char run3bp; double cutoff; /* cutoff radius */ double cutoff_square; /* square of the cutoff radius */ @@ -99,9 +103,18 @@ typedef struct s_moldyn { double t_ref; /* reference temperature */ double t; /* actual temperature */ + double t_sum; /* sum over all t */ + double mean_t; /* mean value of t */ + + t_virial virial; /* global virial (absolute coordinates) */ + double gp; /* pressure computed from global virial */ + double gp_sum; /* sum over all gp */ + double mean_gp; /* mean value of gp */ double p_ref; /* reference pressure */ double p; /* actual pressure (computed by virial) */ + double p_sum; /* sum over all p */ + double mean_p; /* mean value of p */ t_3dvec tp; /* thermodynamic pressure dU/dV */ double dv; /* dV for thermodynamic pressure calc */ @@ -121,7 +134,7 @@ typedef struct s_moldyn { double tau; /* delta t */ double time; /* absolute time */ double tau_square; /* delta t squared */ - double elapsed; /* total elapsed time */ + int total_steps; /* total steps */ double energy; /* potential energy */ double ekin; /* kinetic energy */ @@ -175,7 +188,7 @@ typedef struct s_moldyn { #define P_SCALE_DIRECT 0x08 /* direct p control */ /* - * default values + * default values & units * * - length unit: 1 A (1 A = 1e-10 m) * - time unit: 1 fs (1 fs = 1e-15 s) @@ -191,7 +204,9 @@ typedef struct s_moldyn { #define KILOGRAM (1.0/AMU) /* amu */ #define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) /* A amu / fs^2 */ #define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */ -#define ATM ((1.0133e5*PASCAL)) /* N / A^2 */ +#define BAR ((1.0e5*PASCAL)) /* N / A^2 */ +#define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ +#define EV (1.6021765314e-19*METER*NEWTON) /* NA */ #define MOLDYN_TEMP 273.0 #define MOLDYN_TAU 1.0 @@ -220,17 +235,12 @@ typedef struct s_moldyn { #define QUIET 0 /* - * - * phsical values / constants - * + * potential related phsical values / constants * */ #define ONE_THIRD (1.0/3.0) -#define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ -#define EV (1.6021765314e-19*METER*NEWTON) /* NA */ - #define C 0x06 #define M_C 12.011 /* amu */ @@ -284,10 +294,9 @@ typedef struct s_moldyn { * */ -typedef int (*pf_func1b)(t_moldyn *,t_atom *ai); -typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8 bc); -typedef int (*pf_func2b_post)(t_moldyn *,t_atom *,t_atom *,u8 bc); -typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8 bc); +typedef int (*pf_func1b)(t_moldyn *,t_atom *); +typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8); +typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8); int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); @@ -302,8 +311,11 @@ int set_nn_dist(t_moldyn *moldyn,double dist); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); int set_potential1b(t_moldyn *moldyn,pf_func1b func); int set_potential2b(t_moldyn *moldyn,pf_func2b func); -int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func); -int set_potential3b(t_moldyn *moldyn,pf_func3b func); +int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func); +int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func); +int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func); +int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func); +int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func); int set_potential_params(t_moldyn *moldyn,void *params); int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); @@ -331,8 +343,7 @@ int scale_volume(t_moldyn *moldyn); int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); -double get_e_kin(t_moldyn *moldyn); -double update_e_kin(t_moldyn *moldyn); +double e_kin_calc(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); t_3dvec get_total_p(t_moldyn *moldyn); diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 16f771c..531d292 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -85,10 +85,9 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { 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,d_ij2; - double A,B,R,S,S2,lambda,mu; + double A,R,S,S2,lambda; double f_r,df_r; double f_c,df_c; int brand; @@ -96,31 +95,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double arg; /* use newtons third law */ - //if(aipot_params; 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; - - /* - * 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 - * - */ /* determine cutoff square */ if(brand==ai->brand) @@ -128,71 +106,41 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { else S2=params->S2mixed; - /* dist_ij, d_ij */ + /* 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); /* if d_ij2 > S2 => no force & potential energy contribution */ - if(d_ij2>S2) + 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) { S=params->S[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; - exchange->betajnj=exchange->betaini; - exchange->cj2=exchange->ci2; - exchange->dj2=exchange->di2; - exchange->cj2dj2=exchange->ci2di2; } else { S=params->Smixed; R=params->Rmixed; A=params->Amixed; - B=params->Bmixed; lambda=params->lambda_m; - mu=params->mu_m; - params->exchange.chi=params->chi; - 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_ij, df_r_ij */ f_r=A*exp(-lambda*d_ij); df_r=lambda*f_r/d_ij; - /* 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; - - /* f_c, df_c calc (again, same for ij and ji) */ + /* f_c, df_c */ if(d_ij r > R */ - exchange->d_ij_between_rs=1; } - /* add forces of 2bp (ij, ji) contribution - * dVij = dVji and we sum up both: no 1/2) */ + /* add forces */ v3_add(&(ai->f),&(ai->f),&force); + v3_sub(&(aj->f),&(aj->f),&force); + +#ifdef DEBUG + if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { + printf("force 2bp: [%d %d]\n",ai->tag,aj->tag); + printf("adding %f %f %f\n",force.x,force.y,force.z); + if(ai==&(moldyn->atom[0])) + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + if(aj==&(moldyn->atom[0])) + printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + } +#endif /* virial */ - virial_calc(ai,&force,&dist_ij); + //virial_calc(ai,&force,&dist_ij); + //virial_calc(aj,&force,&dist_ij); //ai->virial.xx-=force.x*dist_ij.x; //ai->virial.yy-=force.y*dist_ij.y; //ai->virial.zz-=force.z*dist_ij.z; @@ -219,270 +176,86 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { //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 and 2bp post processing */ - exchange->run3bp=1; - exchange->run2bp_post=1; - - /* reset 3bp sums */ - exchange->zeta_ij=0.0; - exchange->zeta_ji=0.0; - v3_zero(&(exchange->dzeta_ij)); - v3_zero(&(exchange->dzeta_ji)); + /* energy 2bp contribution */ + moldyn->energy+=f_r*f_c; return 0; } -/* tersoff 2 body post part */ - -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 - * - * that is: - * - zeta_ij, dzeta_ij - * - zeta_ji, dzeta_ji - * - * to compute the 3bp contribution to: - * - Vij, dVij - * - dVji - * - */ +/* tersoff 3 body potential function (first ij loop) */ +int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - - t_3dvec force,temp; - t_3dvec *dist_ij; - double b,db,tmp; - double f_c,df_c,f_a,df_a; - double chi,ni,betaini,nj,betajnj; - double zeta; + unsigned char brand; + double S2; + t_3dvec dist_ij; + double d_ij2,d_ij; params=moldyn->pot_params; exchange=&(params->exchange); - /* 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; - 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); - - /* add force */ - v3_add(&(ai->f),&(ai->f),&force); + /* reset zeta sum */ + exchange->zeta_ij=0.0; - /* virial */ - virial_calc(ai,&force,dist_ij); - //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; + /* + * set ij depending values + */ -#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 + brand=ai->brand; + + if(brand==aj->brand) + S2=params->S2[brand]; + else + S2=params->S2mixed; - /* add energy of 3bp sum */ - moldyn->energy+=(0.5*f_c*b*f_a); + /* 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); - /* 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); + /* if d_ij2 > S2 => no force & potential energy contribution */ + if(d_ij2>S2) { + moldyn->run3bp=0; + return 0; } - v3_scale(&temp,dist_ij,df_c*b*f_a); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,-0.5); - /* add force */ - v3_add(&(ai->f),&(ai->f),&force); - - /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ -// TEST ... with a minus instead - virial_calc(ai,&force,dist_ij); - //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; + /* d_ij */ + d_ij=sqrt(d_ij2); -#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 + /* store values */ + exchange->dist_ij=dist_ij; + exchange->d_ij2=d_ij2; + exchange->d_ij=d_ij; + /* reset k counter for first k loop */ + exchange->kcount=0; + return 0; } -/* tersoff 3 body part */ - -int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { +/* tersoff 3 body potential function (first k loop) */ +int tersoff_mult_3bp_k1(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 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 cos_theta,d_costheta1,d_costheta2; - double h_cos,d2_h_cos2; - double frac,g,zeta,chi; - double tmp; - int brand; + 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; params=moldyn->pot_params; exchange=&(params->exchange); + kcount=exchange->kcount; - if(!(exchange->run3bp)) - return 0; - - /* - * 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 - * - * 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; - - /* f_c_ij, df_c_ij (same for ji) */ - f_c=exchange->f_c; - df_c=exchange->df_c; - - /* - * calculate unknown values now ... - */ - - /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */ - - /* 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); + if(kcount>TERSOFF_MAXN) { + printf("FATAL: neighbours = %d\n",kcount); + printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag); + } /* ik constants */ brand=ai->brand; @@ -497,215 +270,298 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { S2=params->S2mixed; } - /* zeta_ij/dzeta_ij contribution only for d_ik < S */ - if(d_ik2n_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); - } + /* 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++; + return 0; + } + + /* d_ik */ + d_ik=sqrt(d_ik2); + + /* dist_ij, d_ij */ + dist_ij=exchange->dist_ij; + d_ij=exchange->d_ij; + + /* cos theta */ + cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); + + /* g_ijk */ + h_cos=*(exchange->h_i)-cos_theta; + d2_h_cos2=exchange->di2+(h_cos*h_cos); + frac=exchange->ci2/d2_h_cos2; + g=1.0+exchange->ci2di2-frac; + dg=-2.0*frac*h_cos/d2_h_cos2; + + /* 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; + } + 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; } - /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */ + /* 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; - /* dist_jk, d_jk */ - v3_sub(&dist_jk,&(ak->r),&(aj->r)); - if(bc) check_per_bound(moldyn,&dist_jk); - d_jk2=v3_absolute_square(&dist_jk); + /* increase k counter */ + exchange->kcount++; + + return 0; +} + +int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; + t_3dvec force; + double f_a,df_a,b,db,f_c,df_c; + double mu,B,chi; + double d_ij; + unsigned char brand; + double ni,tmp; + double S,R,s_r,arg; + + params=moldyn->pot_params; + exchange=&(params->exchange); - /* jk constants */ brand=aj->brand; - if(brand==ak->brand) { - R=params->R[brand]; + if(brand==ai->brand) { S=params->S[brand]; - S2=params->S2[brand]; + R=params->R[brand]; B=params->B[brand]; mu=params->mu[brand]; chi=1.0; } else { - R=params->Rmixed; S=params->Smixed; - S2=params->S2mixed; + R=params->Rmixed; B=params->Bmixed; mu=params->mu_m; chi=params->chi; } - /* 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); - } - - /* 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; + d_ij=exchange->d_ij; + + /* f_c, df_c */ + if(d_ijzeta_ij==0.0) { + b=chi; + db=0.0; + } + else { + ni=*(exchange->n_i); + tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0); + b=(1.0+exchange->zeta_ij*tmp); + db=chi*pow(b,-1.0/(2*ni)-1.0); + b=db*b; + db*=-0.5*tmp; + } + + /* force contribution */ + v3_scale(&force,&(exchange->dist_ij),df_a*f_c+f_a*df_c); + v3_scale(&force,&force,-0.5*b); + v3_add(&(ai->f),&(ai->f),&force); + v3_sub(&(aj->f),&(aj->f),&force); #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); -} + if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { + 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[0])) + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + if(aj==&(moldyn->atom[0])) + printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->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); + + /* virial */ + //virial_calc(ai,&force,&(exchange->dist_ij)); + //virial_calc(aj,&force,&(exchange->dist_ij)); + + /* dzeta prefactor = - 0.5 f_c f_a db */ + exchange->pre_dzeta=-0.5*f_a*f_c*db; + + /* energy contribution */ + moldyn->energy+=0.5*f_c*b*f_a; + + /* reset k counter for second k loop */ + exchange->kcount=0; + + return 0; } + +/* tersoff 3 body potential function (second k loop) */ +int tersoff_mult_3bp_k2(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + + t_tersoff_mult_params *params; + t_tersoff_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 dcosdri,dcosdrj,dcosdrk; + t_3dvec force,tmp; + + params=moldyn->pot_params; + exchange=&(params->exchange); + kcount=exchange->kcount; + + if(kcount>TERSOFF_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++; + return 0; + } + + /* prefactor dzeta */ + pre_dzeta=exchange->pre_dzeta; + + /* 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 i,j,k */ + dijdik_inv=1.0/(d_ij*d_ik); + v3_scale(&dcosdrj,&dist_ik,dijdik_inv); + v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2); + v3_add(&dcosdrj,&dcosdrj,&tmp); + v3_scale(&dcosdrk,&dist_ij,dijdik_inv); + v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); + v3_add(&dcosdrk,&dcosdrk,&tmp); + v3_add(&dcosdri,&dcosdrj,&dcosdrk); + v3_scale(&dcosdri,&dcosdri,-1.0); + + /* f_c_ik * dg, df_c_ik * g */ + fcdg=f_c_ik*dg; + dfcg=df_c_ik*g; + + /* derivative wrt i */ + v3_scale(&force,&dist_ik,dfcg); + v3_scale(&tmp,&dcosdri,fcdg); + v3_add(&force,&force,&tmp); + v3_scale(&force,&force,pre_dzeta); + + /* force contribution */ + v3_add(&(ai->f),&(ai->f),&force); + +#ifdef DEBUG + if(ai==&(moldyn->atom[0])) { + 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 i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + } +#endif + + /* virial */ + //virial_calc(ai,&force,&dist_ij); + + /* derivatice wrt j */ + v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); + + /* force contribution */ + v3_add(&(aj->f),&(aj->f),&force); + +#ifdef DEBUG + if(aj==&(moldyn->atom[0])) { + 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); + } #endif + /* virial */ + //virial_calc(aj,&force,&dist_ij); + + /* derivative wrt k */ + v3_scale(&force,&dist_ik,dfcg); + 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); + +#ifdef DEBUG + if(ak==&(moldyn->atom[0])) { + 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); } +#endif + + /* virial */ + virial_calc(ak,&force,&dist_ik); + + /* increase k counter */ + exchange->kcount++; return 0; -} +} diff --git a/potentials/tersoff.h b/potentials/tersoff.h index a04ed02..52e848f 100644 --- a/potentials/tersoff.h +++ b/potentials/tersoff.h @@ -8,45 +8,41 @@ #ifndef TERSOFF_H #define TERSOFF_H +#define TERSOFF_MAXN 16*27 + /* tersoff exchange type */ typedef struct s_tersoff_echange { - double f_c,df_c; - double f_a,df_a; t_3dvec dist_ij; double d_ij2; double d_ij; - double chi; + t_3dvec dist_ik[TERSOFF_MAXN]; + double d_ik2[TERSOFF_MAXN]; + double d_ik[TERSOFF_MAXN]; + + double f_c_ik[TERSOFF_MAXN]; + double df_c_ik[TERSOFF_MAXN]; + + double g[TERSOFF_MAXN]; + double dg[TERSOFF_MAXN]; + double cos_theta[TERSOFF_MAXN]; double *beta_i; - double *beta_j; double *n_i; - double *n_j; double *c_i; - double *c_j; double *d_i; - double *d_j; double *h_i; - double *h_j; double ci2; - double cj2; double di2; - double dj2; double ci2di2; - double cj2dj2; double betaini; - double betajnj; - - u8 run3bp; - u8 run2bp_post; - u8 d_ij_between_rs; double zeta_ij; - double zeta_ji; - t_3dvec dzeta_ij; - t_3dvec dzeta_ji; + double pre_dzeta; + + int kcount; } t_tersoff_exchange; /* tersoff mult (2!) potential parameters */ @@ -81,7 +77,11 @@ typedef struct s_tersoff_mult_params { int tersoff_mult_complete_params(t_tersoff_mult_params *p); int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai); 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); -int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int tersoff_mult_3bp_k1(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int tersoff_mult_3bp_k2(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); #endif diff --git a/potentials/tersoff_orig.c b/potentials/tersoff_orig.c new file mode 100644 index 0000000..90121f2 --- /dev/null +++ b/potentials/tersoff_orig.c @@ -0,0 +1,707 @@ +/* + * tersoff_orig.c - tersoff potential + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../moldyn.h" +#include "../math/math.h" +#include "tersoff_orig.h" + +/* create mixed terms from parameters and set them */ +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]); + p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); + p->mu_m=0.5*(p->mu[0]+p->mu[1]); + + printf("[moldyn] tersoff mult parameter info:\n"); + 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], + p->lambda_m); + printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); + printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]); + printf(" n | %f | %f\n",p->n[0],p->n[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(" chi | %f \n",p->chi); + + return 0; +} + +/* tersoff 1 body part */ +int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { + + int brand; + t_tersoff_mult_params *params; + t_tersoff_exchange *exchange; + + brand=ai->brand; + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* + * simple: point constant parameters only depending on atom i to + * their right values + */ + + 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->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; +} + +/* tersoff 2 body part */ +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,d_ij2; + double A,B,R,S,S2,lambda,mu; + double f_r,df_r; + double f_c,df_c; + int brand; + double s_r; + double arg; + + /* use newtons third law */ + //if(aipot_params; + 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; + + /* + * 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 + * + */ + + /* determine cutoff square */ + if(brand==ai->brand) + S2=params->S2[brand]; + else + S2=params->S2mixed; + + /* 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) { + S=params->S[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; + exchange->betajnj=exchange->betaini; + exchange->cj2=exchange->ci2; + exchange->dj2=exchange->di2; + exchange->cj2dj2=exchange->ci2di2; + } + else { + S=params->Smixed; + R=params->Rmixed; + A=params->Amixed; + B=params->Bmixed; + lambda=params->lambda_m; + mu=params->mu_m; + exchange->chi=params->chi; + 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; + + /* 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; + + /* f_c, df_c calc (again, same for ij and ji) */ + if(d_ij r > R */ + exchange->d_ij_between_rs=1; + } + + /* add forces of 2bp (ij, ji) contribution + * dVij = dVji and we sum up both: no 1/2) */ + v3_add(&(ai->f),&(ai->f),&force); + + /* virial */ + virial_calc(ai,&force,&dist_ij); + //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: [%d %d]\n",ai->tag,aj->tag); + printf("adding %f %f %f\n",force.x,force.y,force.z); + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,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 and 2bp post processing */ + exchange->run3bp=1; + exchange->run2bp_post=1; + + /* reset 3bp sums */ + exchange->zeta_ij=0.0; + exchange->zeta_ji=0.0; + v3_zero(&(exchange->dzeta_ij)); + v3_zero(&(exchange->dzeta_ji)); + + return 0; +} + +/* tersoff 2 body post part */ + +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 + * + * 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; + t_3dvec *dist_ij; + double b,db,tmp; + double f_c,df_c,f_a,df_a; + double chi,ni,betaini,nj,betajnj; + double zeta; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* 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; + 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); + + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); + + /* virial */ + virial_calc(ai,&force,dist_ij); + //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: [%d %d sum]\n",ai->tag,aj->tag); + printf("adding %f %f %f\n",force.x,force.y,force.z); + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,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 + + /* 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); + + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); + + /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ +// TEST ... with a minus instead + virial_calc(ai,&force,dist_ij); + //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: [%d %d sum]\n",ai->tag,aj->tag); + printf("adding %f %f %f\n",force.x,force.y,force.z); + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,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; +} + +/* tersoff 3 body part */ + +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 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 cos_theta,d_costheta1,d_costheta2; + double h_cos,d2_h_cos2; + double frac,g,zeta,chi; + double tmp; + int brand; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + if(!(exchange->run3bp)) + return 0; + + /* + * 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 + * + * 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; + + /* f_c_ij, df_c_ij (same for ji) */ + f_c=exchange->f_c; + df_c=exchange->df_c; + + /* + * calculate unknown values now ... + */ + + /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */ + + /* 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; + } + + /* zeta_ij/dzeta_ij contribution only for d_ik < S */ + if(d_ik2n_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); + } + } + + /* 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_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 { + R=params->Rmixed; + S=params->Smixed; + S2=params->S2mixed; + B=params->Bmixed; + mu=params->mu_m; + chi=params->chi; + } + + /* 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); + } + + /* 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: [%d %d %d]\n",ai->tag,aj->tag,ak->tag); + printf("adding %f %f %f\n",temp2.x,temp2.y,temp2.z); + printf("total i: %f %f %f\n",ai->f.x,ai->f.y,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 + + } + + return 0; +} + diff --git a/potentials/tersoff_orig.h b/potentials/tersoff_orig.h new file mode 100644 index 0000000..d5d4a66 --- /dev/null +++ b/potentials/tersoff_orig.h @@ -0,0 +1,87 @@ +/* + * tersoff_orig.h - tersoff potential header file + * + * author: Frank Zirkelbach + * + */ + +#ifndef TERSOFF_H +#define TERSOFF_H + +/* tersoff exchange type */ +typedef struct s_tersoff_echange { + double f_c,df_c; + double f_a,df_a; + + t_3dvec dist_ij; + double d_ij2; + double d_ij; + + double chi; + + double *beta_i; + double *beta_j; + double *n_i; + double *n_j; + double *c_i; + double *c_j; + double *d_i; + double *d_j; + double *h_i; + double *h_j; + + double ci2; + double cj2; + double di2; + double dj2; + double ci2di2; + double cj2dj2; + double betaini; + double betajnj; + + u8 run3bp; + u8 run2bp_post; + u8 d_ij_between_rs; + + double zeta_ij; + double zeta_ji; + t_3dvec dzeta_ij; + t_3dvec dzeta_ji; +} t_tersoff_exchange; + +/* tersoff mult (2!) potential parameters */ +typedef struct s_tersoff_mult_params { + double S[2]; /* tersoff cutoff radii */ + double S2[2]; /* tersoff cutoff radii squared */ + double R[2]; /* tersoff cutoff radii */ + double Smixed; /* mixed S radius */ + double S2mixed; /* mixed S radius squared */ + double Rmixed; /* mixed R radius */ + double A[2]; /* factor of tersoff attractive part */ + double B[2]; /* factor of tersoff repulsive part */ + double Amixed; /* mixed A factor */ + double Bmixed; /* mixed B factor */ + double lambda[2]; /* tersoff lambda */ + double lambda_m; /* mixed lambda */ + double mu[2]; /* tersoff mu */ + double mu_m; /* mixed mu */ + + double chi; + + double beta[2]; + double n[2]; + double c[2]; + double d[2]; + double h[2]; + + t_tersoff_exchange exchange; /* exchange between 2bp and 3bp calc */ +} t_tersoff_mult_params; + +/* function prototypes */ +int tersoff_mult_complete_params(t_tersoff_mult_params *p); +int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai); +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); +int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); + +#endif diff --git a/report/report.h b/report/report.h index 4fe42f7..bd4ce88 100644 --- a/report/report.h +++ b/report/report.h @@ -87,10 +87,10 @@ set xtic auto \n\ set ytic auto \n\ set title 'Pressure vs. time' \n\ set xlabel 'Time [fs]' \n\ -set ylabel 'Pressure [atm]' \n\ +set ylabel 'Pressure [bar]' \n\ set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\ set output 'pressure.eps' \n\ -plot \"pressure\" using 1:2 title 'Pressure' with lines \ +plot \"pressure\" using 1:2 title 'P' with lines , \"pressure\" using 1:3 title '

' with lines , \"pressure\" using 1:4 title 'P (global virial)' with lines , \"pressure\" using 1:5 title '

' with lines \ "; static char temperature_plot_script[]="\ @@ -104,7 +104,7 @@ set xlabel 'Time [fs]' \n\ set ylabel 'Temperature [K]' \n\ set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\ set output 'temperature.eps' \n\ -plot \"temperature\" using 1:2 title 'Temperature' with lines \ +plot \"temperature\" using 1:2 title 'T' with lines , \"temperature\" using 1:3 title '' with lines \ "; #endif diff --git a/sic.c b/sic.c index 682d89e..a78f3a0 100644 --- a/sic.c +++ b/sic.c @@ -14,9 +14,10 @@ #include "potentials/harmonic_oscillator.h" #include "potentials/lennard_jones.h" #include "potentials/tersoff.h" +//#include "potentials/tersoff_orig.h" -#define INJECT 20 -#define NR_ATOMS 20 +#define INJECT 1 +#define NR_ATOMS 4 int hook(void *moldyn,void *hook_params) { @@ -97,8 +98,13 @@ int main(int argc,char **argv) { /* choose potential */ set_potential1b(&md,tersoff_mult_1bp); set_potential2b(&md,tersoff_mult_2bp); - set_potential2b_post(&md,tersoff_mult_post_2bp); - set_potential3b(&md,tersoff_mult_3bp); + //set_potential3b_j1(&md,tersoff_mult_2bp); + //set_potential3b_k1(&md,tersoff_mult_3bp); + //set_potential3b_j3(&md,tersoff_mult_post_2bp); + set_potential3b_j1(&md,tersoff_mult_3bp_j1); + set_potential3b_k1(&md,tersoff_mult_3bp_k1); + set_potential3b_j2(&md,tersoff_mult_3bp_j2); + set_potential3b_k2(&md,tersoff_mult_3bp_k2); //set_potential2b(&md,lennard_jones); //set_potential2b(&md,harmonic_oscillator); set_potential_params(&md,&tp); @@ -198,7 +204,7 @@ int main(int argc,char **argv) { /* set temperature & pressure */ set_temperature(&md,atof(argv[2])+273.0); - set_pressure(&md,ATM); + set_pressure(&md,BAR); /* set p/t scaling */ //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001, @@ -211,14 +217,14 @@ int main(int argc,char **argv) { /* create the simulation schedule */ /* initial configuration */ - moldyn_add_schedule(&md,500,1.0); + moldyn_add_schedule(&md,10000,1.0); /* adding atoms */ - for(inject=0;inject