X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=def5079dbbd949fa55e4f65b9aa4a77e165077db;hb=4c2140b0f76fb191bdd9b9c2a329877eb0aae531;hp=84cf7003c13f70cd6e7180a7a23b53ddda6722da;hpb=2346d076f3d0955fb3a7bfb47f2b17b9eae6dd5b;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 84cf700..def5079 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -20,6 +20,17 @@
 #include "moldyn.h"
 #include "report/report.h"
 
+/* potential includes */
+#include "potentials/harmonic_oscillator.h"
+#include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
+#include "potentials/tersoff.h"
+#endif
+
+
 /*
  * global variables, pse and atom colors (only needed here)
  */ 
@@ -228,58 +239,29 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 	return 0;
 }
 
-int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
-
-	moldyn->func1b=func;
-
-	return 0;
-}
-
-int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func2b=func;
-
-	return 0;
-}
-
-int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func3b_j1=func;
-
-	return 0;
-}
-
-int set_potential3b_j2(t_moldyn *moldyn,pf_func2b 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;
-}
-
-int set_potential_params(t_moldyn *moldyn,void *params) {
+int set_potential(t_moldyn *moldyn,u8 type) {
 
-	moldyn->pot_params=params;
+	switch(type) {
+		case MOLDYN_POTENTIAL_TM:
+			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;
+			// missing: check 2b bond func
+			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->check_2b_bond=albe_mult_check_2b_bond;
+			break;
+		default:
+			printf("[moldyn] set potential: unknown type %02x\n",
+			       type);
+			return -1;
+	}
 
 	return 0;
 }
@@ -2508,8 +2490,8 @@ int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
 	if(bc) check_per_bound(moldyn,&dist);
 	d=v3_absolute_square(&dist);
 
-	/* ignore if greater or equal 2 times cutoff */
-	if(d>=4.0*moldyn->cutoff_square)
+	/* ignore if greater cutoff */
+	if(d>moldyn->cutoff_square)
 		return 0;
 
 	/* fill the slots */
@@ -2548,7 +2530,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 	int i;
 
 	pcc.dr=dr;
-	pcc.o1=2.0*moldyn->cutoff/dr;
+	pcc.o1=moldyn->cutoff/dr;
 	pcc.o2=2*pcc.o1;
 
 	if(pcc.o1*dr<=moldyn->cutoff)
@@ -2612,9 +2594,15 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	d=v3_absolute_square(&dist);
 
 	/* ignore if greater or equal cutoff */
-	if(d>=moldyn->cutoff_square)
+	if(d>moldyn->cutoff_square)
 		return 0;
 
+	/* check for potential bond */
+	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+		return 0;
+
+	d=sqrt(d);
+
 	/* now count this bonding ... */
 	ba=data;
 
@@ -2641,6 +2629,7 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	// by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
 
 	int qcnt;
+	int ccnt,cset;
 	t_ba ba;
 	int i;
 	t_atom *atom;
@@ -2660,8 +2649,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
 	ba.tcnt=0;
-
 	qcnt=0;
+	ccnt=0;
+	cset=0;
 
 	atom=moldyn->atom;
 
@@ -2673,17 +2663,25 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 				qcnt+=4;
 		}
 		else {
-			if((ba.acnt[i]==4)&(ba.bcnt[i]==0))
+			if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
 				qcnt+=4;
+				ccnt+=1;
+			}
+			cset+=1;
 		}
 	}
 
-printf("%d %d\n",qcnt,ba.tcnt);
-	if(quality)
-		*quality=1.0*qcnt/ba.tcnt;
-	else
-		printf("[moldyn] bond analyze: quality = %f\n",
-		       1.0*qcnt/ba.tcnt);
+	printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
+	printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
+
+	if(quality) {
+		quality[0]=1.0*ccnt/cset;
+		quality[1]=1.0*qcnt/ba.tcnt;
+	}
+	else {
+		printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
+		printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
+	}
 
 	return 0;
 }