From: hackbard Date: Mon, 5 Mar 2007 09:01:53 +0000 (+0000) Subject: new src file layout (warning: doesnt compile by now!) X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=caa3bc828974c35df2462fde737c31c0a618ee4e;p=physik%2Fposic.git new src file layout (warning: doesnt compile by now!) --- diff --git a/Makefile b/Makefile index c51f8bd..37fafcc 100644 --- a/Makefile +++ b/Makefile @@ -8,13 +8,16 @@ CFLAGS+=-g #CFLAGS+=-DVDEBUG LDFLAGS=-lm -OBJS=visual/visual.o random/random.o -OBJS+=moldyn.o +SOURCE=moldyn.c visual/visual.c random/random.c -all: clean sic +POT_SRC=potentials/tersoff.c potentials/lennard_jones.c +POT_SRC+= potentials/harmonic_oscillator.c -sic: $(OBJS) +all: sic + +sic: + $(CC) $(CFLAGS) $(LDFLAGS) $(SOURCE) $(POT_SRC) sic.c -o sic .PHONY:clean clean: - rm -f *.o sic */*.o + rm -f sic diff --git a/moldyn.c b/moldyn.c index f6c3c81..1ab88eb 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1531,740 +1531,6 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { return 0; } -/* lennard jones potential & force for one sort of atoms */ - -int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - - t_lj_params *params; - t_3dvec force,distance; - double d,h1,h2; - double eps,sig6,sig12; - - params=moldyn->pot2b_params; - eps=params->epsilon4; - sig6=params->sigma6; - sig12=params->sigma12; - - if(air),&(ai->r)); - if(bc) check_per_bound(moldyn,&distance); - d=v3_absolute_square(&distance); /* 1/r^2 */ - if(d<=moldyn->cutoff_square) { - d=1.0/d; /* 1/r^2 */ - h2=d*d; /* 1/r^4 */ - h2*=d; /* 1/r^6 */ - h1=h2*h2; /* 1/r^12 */ - moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); - h2*=d; /* 1/r^8 */ - h1*=d; /* 1/r^14 */ - h2*=6*sig6; - h1*=12*sig12; - d=+h1-h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(aj->f),&(aj->f),&force); - v3_scale(&force,&force,-1.0); /* f = - grad E */ - v3_add(&(ai->f),&(ai->f),&force); - virial_calc(ai,&force,&distance); -if(force.x*distance.x<=0) printf("virial xx: %.15f -> %f %f %f\n",force.x*distance.x,distance.x,distance.y,distance.z); - virial_calc(aj,&force,&distance); /* f and d signe switched */ - } - - return 0; -} - -/* - * tersoff potential & force for 2 sorts of atoms - */ - -/* 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->pot1b_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; - - params=moldyn->pot2b_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 - * - */ - - /* constants */ - 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; - lambda=params->lambda_m; - mu=params->mu_m; - params->exchange.chi=params->chi; - } - - /* 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; - - /* 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 */ - 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 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->pot2b_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 */ - 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 - - /* 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 - 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; -} - -/* 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->pot3b_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:\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 - - } - - return 0; -} - - /* * debugging / critical check functions */ diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c new file mode 100644 index 0000000..0d933b6 --- /dev/null +++ b/potentials/harmonic_oscillator.c @@ -0,0 +1,53 @@ +/* + * harmonic_oscillator.c - harmonic oscillator potential + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../moldyn.h" +#include "../math/math.h" +//#include "harmonic_oscillator.h" + +int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_ho_params *params; + t_3dvec force,distance; + double d,f; + double sc,equi_dist; + + params=moldyn->pot2b_params; + sc=params->spring_constant; + equi_dist=params->equilibrium_distance; + + if(air),&(ai->r)); + + if(bc) check_per_bound(moldyn,&distance); + d=v3_norm(&distance); + if(d<=moldyn->cutoff) { + 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; +} + diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c new file mode 100644 index 0000000..7d45a17 --- /dev/null +++ b/potentials/lennard_jones.c @@ -0,0 +1,60 @@ +/* + * lennard_jones.c - lennard jones potential + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../moldyn.h" +#inlcude "../math/math.h" +//#include "lennard_jones.h" + +int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_lj_params *params; + t_3dvec force,distance; + double d,h1,h2; + double eps,sig6,sig12; + + params=moldyn->pot2b_params; + eps=params->epsilon4; + sig6=params->sigma6; + sig12=params->sigma12; + + if(air),&(ai->r)); + if(bc) check_per_bound(moldyn,&distance); + d=v3_absolute_square(&distance); /* 1/r^2 */ + if(d<=moldyn->cutoff_square) { + d=1.0/d; /* 1/r^2 */ + h2=d*d; /* 1/r^4 */ + h1=h2*h2; /* 1/r^12 */ + moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); + h2*=d; /* 1/r^8 */ + h1*=d; /* 1/r^14 */ + h2*=6*sig6; + h1*=12*sig12; + d=+h1-h2; + d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(aj->f),&(aj->f),&force); + v3_scale(&force,&force,-1.0); /* 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 */ + } + + return 0; +} + diff --git a/potentials/tersoff.c b/potentials/tersoff.c new file mode 100644 index 0000000..87bb187 --- /dev/null +++ b/potentials/tersoff.c @@ -0,0 +1,707 @@ +/* + * tersoff.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.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->pot1b_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; + + params=moldyn->pot2b_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 + * + */ + + /* constants */ + 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; + lambda=params->lambda_m; + mu=params->mu_m; + params->exchange.chi=params->chi; + } + + /* 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; + + /* 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 */ + 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 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->pot2b_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 */ + 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 + + /* 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 + 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; +} + +/* 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->pot3b_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:\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 + + } + + return 0; +} +