X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=bf7b9f0cde519d988f962b6576a6c58d947142c8;hb=959801961cfdd51e080e6500f51647b3397d336c;hp=82ae543c073457a47372253f6c238f76d1bfe6b5;hpb=0d1dfb1e5027d215fced8ca306dd658f692a2a44;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 82ae543..bf7b9f0 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -30,101 +30,6 @@
 #include "potentials/tersoff.h"
 #endif
 
-
-/*
- * global variables, pse and atom colors (only needed here)
- */ 
-
-static char *pse_name[]={
-	"*",
-	"H",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"C",
-	"N",
-	"O",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Si",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-static char *pse_col[]={
-	"*",
-	"White",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"Gray",
-	"N",
-	"Blue",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Yellow",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-/*
-static double pse_mass[]={
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	M_C,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	M_SI,
-	0,
-	0,
-	0,
-	0,
-};
-
-static double pse_lc[]={
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	LC_C,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	LC_SI,
-	0,
-	0,
-	0,
-	0,
-};
-*/
-
 /*
  * the moldyn functions
  */
@@ -208,7 +113,7 @@ int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
 	moldyn->pt_scale|=ptype;
 	moldyn->p_tc=ptc;
 
-	printf("[moldyn] p/t scaling:\n");
+	printf("[moldyn] p scaling:\n");
 
 	printf("  p: %s",ptype?"yes":"no ");
 	if(ptype)
@@ -224,7 +129,7 @@ int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
 	moldyn->pt_scale|=ttype;
 	moldyn->t_tc=ttc;
 
-	printf("[moldyn] p/t scaling:\n");
+	printf("[moldyn] t scaling:\n");
 
 	printf("  t: %s",ttype?"yes":"no ");
 	if(ttype)
@@ -315,8 +220,9 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 			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
+			moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
 			break;
+		/*
 		case MOLDYN_POTENTIAL_AM:
 			moldyn->func3b_j1=albe_mult_3bp_j1;
 			moldyn->func3b_k1=albe_mult_3bp_k1;
@@ -324,6 +230,17 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 			moldyn->func3b_k2=albe_mult_3bp_k2;
 			moldyn->check_2b_bond=albe_mult_check_2b_bond;
 			break;
+		*/
+		case MOLDYN_POTENTIAL_AM:
+			moldyn->func1b=albe_mult_i0;
+			moldyn->func2b=albe_mult_i0_j0;
+			moldyn->func3b_0=albe_mult_i0_j0_k0;
+			moldyn->func3b_1=albe_mult_i0_j1;
+			moldyn->func3b_j1=albe_mult_i0_j2;
+			moldyn->func3b_k2=albe_mult_i0_j2_k0;
+			moldyn->func3b_j2=albe_mult_i0_j3;
+			moldyn->check_2b_bond=albe_mult_check_2b_bond;
+			break;
 		case MOLDYN_POTENTIAL_HO:
 			moldyn->func2b=harmonic_oscillator;
 			moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
@@ -1674,8 +1591,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		printf("[moldyn] WARNING: forces too high / tau too small!\n");
 
 	/* zero absolute time */
-	moldyn->time=0.0;
-	moldyn->total_steps=0;
+	// should have right values!
+	//moldyn->time=0.0;
+	//moldyn->total_steps=0;
 
 	/* debugging, ignore */
 	moldyn->debug=0;
@@ -1995,6 +1913,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 					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,
@@ -2002,6 +1923,78 @@ int potential_force_calc(t_moldyn *moldyn) {
 						               jtom,
 					        	       bc_ij);
 					}
+
+// REWRITE ONCE WORKING!!!
+					/* 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->func3b_0==NULL)
+						continue;
+
+				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->func3b_0(moldyn,
+						                 &(itom[i]),
+						                 jtom,
+						                 ktom,
+						                 bc_ik|bc_ij);
+#ifdef STATIC_LISTS
+					}
+#else
+					} while(list_next_f(that)!=\
+					        L_NO_NEXT_ELEMENT);
+#endif
+
+				}
+
+
+				if(moldyn->func3b_1)
+					moldyn->func3b_1(moldyn,
+					                 &(itom[i]),
+					                 jtom,
+					                 bc_ij);
+
+// UNTIL HERE + change function names!
+
+
+
 				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 
@@ -2681,8 +2674,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
 		return 0;
 
-	d=sqrt(d);
-
 	/* now count this bonding ... */
 	ba=data;