From: hackbard Date: Mon, 20 Feb 2012 15:55:56 +0000 (+0100) Subject: Merge branch 'leadoff' X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;p=physik%2Fposic.git Merge branch 'leadoff' Conflicts: Makefile moldyn.c potentials/albe.c potentials/albe.h runmd ... hopefully fixed! Just a test to merge leadoff back to master! :) --- 6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9 diff --cc moldyn.c index 010b8ed,33224a6..1d6e0b7 --- a/moldyn.c +++ b/moldyn.c @@@ -216,28 -263,18 +264,28 @@@ int set_potential(t_moldyn *moldyn,u8 t switch(type) { case MOLDYN_POTENTIAL_TM: - moldyn->func_i0=tersoff_mult_1bp; - moldyn->func_j1=tersoff_mult_3bp_j1; - moldyn->func_j1_k0=tersoff_mult_3bp_k1; - moldyn->func_j1c=tersoff_mult_3bp_j2; - moldyn->func_j1_k1=tersoff_mult_3bp_k2; + //moldyn->func1b=tersoff_mult_1bp; + moldyn->func3b_j1=tersoff_mult_3bp_j1; + moldyn->func3b_k1=tersoff_mult_3bp_k1; + moldyn->func3b_j2=tersoff_mult_3bp_j2; + moldyn->func3b_k2=tersoff_mult_3bp_k2; moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; + case MOLDYN_POTENTIAL_AO: + moldyn->func_j1=albe_orig_mult_3bp_j1; + moldyn->func_j1_k0=albe_orig_mult_3bp_k1; + moldyn->func_j1c=albe_orig_mult_3bp_j2; + moldyn->func_j1_k1=albe_orig_mult_3bp_k2; + moldyn->check_2b_bond=albe_orig_mult_check_2b_bond; + break; case MOLDYN_POTENTIAL_AM: - moldyn->func3b_j1=albe_mult_3bp_j1; - moldyn->func3b_k1=albe_mult_3bp_k1; - moldyn->func3b_j2=albe_mult_3bp_j2; - moldyn->func3b_k2=albe_mult_3bp_k2; + moldyn->func_i0=albe_mult_i0; + moldyn->func_j0=albe_mult_i0_j0; + moldyn->func_j0_k0=albe_mult_i0_j0_k0; + moldyn->func_j0e=albe_mult_i0_j1; + moldyn->func_j1=albe_mult_i0_j2; + moldyn->func_j1_k0=albe_mult_i0_j2_k0; + moldyn->func_j1c=albe_mult_i0_j3; moldyn->check_2b_bond=albe_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_HO: @@@ -1908,86 -2536,17 +2553,74 @@@ int potential_force_calc(t_moldyn *mold if(jtom==&(itom[i])) continue; + /* reset 3bp run */ + moldyn->run3bp=1; + if((jtom->attr&ATOM_ATTR_2BP)& (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + moldyn->func_j0(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + + /* 3 body potential/force */ + + /* in j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) + continue; + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; + + if(moldyn->func_j0_k0==NULL) + continue; + + /* first loop over atoms k in first j loop */ + for(k=0;k<27;k++) { + + bc_ik=(kstart==NULL) + continue; + + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + //if(ktom==jtom) + // continue; + + if(ktom==&(itom[i])) + continue; + + moldyn->func_j0_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS } - #else - } while(list_next_f(that)!=\ - L_NO_NEXT_ELEMENT); - #endif - - } - - /* finish of first j loop after first k loop */ - if(moldyn->func_j0e) - moldyn->func_j0e(moldyn, - &(itom[i]), - jtom, - bc_ij); - #ifdef STATIC_LISTS } + #elif LOWMEM_LISTS + } #else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif diff --cc potentials/albe.c index 9d8336e,1ec3938..9f98282 --- a/potentials/albe.c +++ b/potentials/albe.c @@@ -271,88 -239,52 +280,66 @@@ int albe_mult_i0_j0_k0(t_moldyn *moldyn return 0; } - /* d_ik */ - d_ik=sqrt(d_ik2); - - /* dist_ij, d_ij */ - dist_ij=exchange->dist_ij; - d_ij=exchange->d_ij; + /* distances */ + distj=exchange->dist[j]; + distk=exchange->dist[k]; + dj=exchange->d[j]; + dk=exchange->d[k]; + djdk_inv=1.0/(dj*dk); /* cos theta */ - cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - - /* g_ijk */ - h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism - d2_h_cos2=exchange->di2+(h_cos*h_cos); - frac=exchange->ci2/d2_h_cos2; - g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac); - dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. - - /* zeta sum += f_c_ik * g_ijk */ - if(d_ik<=R) { - exchange->zeta_ij+=g; - f_c_ik=1.0; - df_c_ik=0.0; + cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv; + + /* g(cos(theta)) ij and ik values */ + h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism + d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j); + frac_j=*exchange->c2_[j]/d2_h_cos2_j; + gj=1.0+*exchange->c2d2_[j]-frac_j; + gj*=*(exchange->gamma_[j]); + dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe + if(ak->brand==aj->brand) { + gk=gj; + dgk=dgj; } else { - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); - exchange->zeta_ij+=f_c_ik*g; + h_cos_k=*(exchange->h_[k])+cos_theta; + d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k); + frac_k=*exchange->c2_[k]/d2_h_cos2_k; + gk=1.0+*exchange->c2d2_[k]-frac_k; + gk*=*(exchange->gamma_[k]); + dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k; } - /* zeta - for albe: ik depending g function */ - //if(ai->tag==0) { - // printf("------> %.15f %.15f\n",dj,dk); - // printf("------> %.15f %.15f\n",dj,dk); - //} - - exchange->zeta[j]+=(exchange->f_c[k]*gk); - exchange->zeta[k]+=(exchange->f_c[j]*gj); - - /* cos theta derivatives */ - v3_scale(&dcosdrj,&distk,djdk_inv); // j - v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]); - v3_add(&dcosdrj,&dcosdrj,&tmp); - v3_scale(&dcosdrk,&distj,djdk_inv); // k - v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]); - v3_add(&dcosdrk,&dcosdrk,&tmp); + #ifdef DEBUG + if(ai==&(moldyn->atom[DATOM])) + printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik); + #endif - /* zeta derivatives */ - dzjj=&(exchange->dzeta[j][j]); - dzkk=&(exchange->dzeta[k][k]); - dzjk=&(exchange->dzeta[j][k]); - dzkj=&(exchange->dzeta[k][j]); - v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk); - v3_add(dzjj,dzjj,&tmp); // j j - v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj); - v3_add(dzkk,dzkk,&tmp); // k k - v3_scale(&tmp,&distk,-exchange->df_c[k]*gk); // dk rik = - di rik - v3_add(dzjk,dzjk,&tmp); - v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk); - v3_add(dzjk,dzjk,&tmp); // j k - v3_scale(&tmp,&distj,-exchange->df_c[j]*gj); // dj rij = - di rij - v3_add(dzkj,dzkj,&tmp); - v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj); - v3_add(dzkj,dzkj,&tmp); // k j + /* store even more data for second k loop */ + exchange->g[kcount]=g; + exchange->dg[kcount]=dg; + exchange->d_ik[kcount]=d_ik; + exchange->cos_theta[kcount]=cos_theta; + exchange->f_c_ik[kcount]=f_c_ik; + exchange->df_c_ik[kcount]=df_c_ik; /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; + + return 0; +} + +/* first j loop within first i loop */ +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->jcnt+=1; return 0; } @@@ -440,16 -368,23 +427,26 @@@ printf(" t: %.15f %.15f %.15f\n",ai- v3_scale(&force,&force,-1.0); // dri rij = - drj rij v3_add(&(aj->f),&(aj->f),&force); - #ifdef NDEBUG - if(aj->tag==0) { - printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]); - printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); - } + /* virial */ + virial_calc(ai,&force,&(exchange->dist_ij)); + + #ifdef DEBUG + if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) { + printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + if(ai==&(moldyn->atom[DATOM])) + printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + if(aj==&(moldyn->atom[DATOM])) + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), + f_c,b,f_a,f_r); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); + } #endif + /* virial */ + virial_calc(ai,&force,dist); + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ exchange->pre_dzeta=0.5*f_a*f_c*db; @@@ -526,32 -496,37 +555,50 @@@ int albe_mult_i0_j2_k0(t_moldyn *moldyn v3_scale(&force,&force,-1.0); v3_add(&(ai->f),&(ai->f),&force); - #ifdef NDEBUG - if(ai->tag==0) { - printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k); - printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); - printf(" ## %f\n",exchange->d[k]); - } + /* derivative wrt k */ + v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik + v3_scale(&tmp,&dcosdrk,fcdg); + v3_add(&force,&force,&tmp); + v3_scale(&force,&force,pre_dzeta); + - /* force contribution */ - v3_add(&(ak->f),&(ak->f),&force); ++ v3_scale(&force,&force,-1.0); ++ v3_add(&(ai->f),&(ai->f),&force); + + #ifdef DEBUG + if(ak==&(moldyn->atom[DATOM])) { + printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } #endif + /* virial */ + virial_calc(ai,&force,&dist_ik); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; return 0; +} + +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->j2cnt+=1; + + return 0; } int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) { diff --cc potentials/albe.h index 39c4df8,89cdae2..9687655 --- a/potentials/albe.h +++ b/potentials/albe.h @@@ -62,17 -65,17 +62,19 @@@ typedef struct s_albe_mult_params double gamma[2]; double gamma_m; double c[2]; + double c2[2]; double c_mixed; + double c2[2]; double c2_mixed; double d[2]; + double d2[2]; double d_mixed; + double d2[2]; double d2_mixed; - double c2d2[2]; - double c2d2_m; double h[2]; double h_mixed; + double c2d2[2]; + double c2d2_m; t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */ } t_albe_mult_params;