From: hackbard <hackbard@hackdaworld.org>
Date: Mon, 20 Feb 2012 15:55:56 +0000 (+0100)
Subject: Merge branch 'leadoff'
X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9

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=(k<dnlc)?0:1;
 +#ifdef STATIC_LISTS
 +					q=0;
 +
 +					while(neighbour_i[j][q]!=0) {
 +
 +						ktom=&(atom[neighbour_i[k][q]]);
 +						q++;
 +#else
 +					that=&(neighbour_i2[k]);
 +					list_reset_f(that);
 +					
 +					if(that->start==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;