X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1d6e0b7dec9e3032df43f33baa692ba078bb5f7a;hp=ccc0b56f99159f0bab42b18051116241355eb0f8;hb=HEAD;hpb=bc17c58e6e7c5adde9e7cf89efd6dfe37147446d

diff --git a/moldyn.c b/moldyn.c
index ccc0b56..1d6e0b7 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -34,6 +34,7 @@
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
 #include "potentials/albe.h"
+#include "potentials/albe_orig.h"
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
 #else
@@ -55,12 +56,10 @@ pthread_mutex_t *amutex;
 pthread_mutex_t emutex;
 #endif
 
-#ifdef CRT
 /* fully constrained relaxation technique - global pointers */
-u8 crt;
+u8 crtt;
 u8 *constraints;
 double *trafo_angle;
-#endif
 
 /*
  * the moldyn functions
@@ -85,14 +84,6 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 	pthread_mutex_init(&emutex,NULL);
 #endif
 
-#ifdef CONSTRAINT_110_5832
-	printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
-	printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
-#endif
-#ifdef CONSTRAINT_11X_5832
-	printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
-	printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
-#endif
 	return 0;
 }
 
@@ -280,19 +271,29 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 			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:
-			moldyn->func2b=harmonic_oscillator;
+			moldyn->func_j0=harmonic_oscillator;
 			moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
 			break;
 		case MOLDYN_POTENTIAL_LJ:
-			moldyn->func2b=lennard_jones;
+			moldyn->func_j0=lennard_jones;
 			moldyn->check_2b_bond=lennard_jones_check_2b_bond;
 			break;
 		default:
@@ -2267,6 +2268,43 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 	return 0;
 }
 
+/* basis trafo */
+
+#define FORWARD		0
+#define BACKWARD	1
+
+int basis_trafo(t_3dvec *r,u8 dir,double z,double x) {
+
+	t_3dvec tmp;
+
+	if(dir==FORWARD) {
+		if(z!=0.0) {
+			v3_copy(&tmp,r);
+			r->x=cos(z)*tmp.x-sin(z)*tmp.y;
+			r->y=sin(z)*tmp.x+cos(z)*tmp.y;
+		}
+		if(x!=0.0) {
+			v3_copy(&tmp,r);
+			r->y=cos(x)*tmp.y-sin(x)*tmp.z;
+			r->z=sin(x)*tmp.y+cos(x)*tmp.z;
+		}
+	}
+	else {
+		if(x!=0.0) {
+			v3_copy(&tmp,r);
+			r->y=cos(-x)*tmp.y-sin(-x)*tmp.z;
+			r->z=sin(-x)*tmp.y+cos(-x)*tmp.z;
+		}
+		if(z!=0.0) {
+			v3_copy(&tmp,r);
+			r->x=cos(-z)*tmp.x-sin(-z)*tmp.y;
+			r->y=sin(-z)*tmp.x+cos(-z)*tmp.y;
+		}
+	}
+
+	return 0;
+}
+
 /* velocity verlet */
 
 int velocity_verlet(t_moldyn *moldyn) {
@@ -2275,113 +2313,52 @@ int velocity_verlet(t_moldyn *moldyn) {
 	double tau,tau_square,h;
 	t_3dvec delta;
 	t_atom *atom;
-#ifdef CONSTRAINT_11X_5832
-	double xt,yt,zt;
-	double xtt,ytt,ztt;
-#endif
 
 	atom=moldyn->atom;
 	count=moldyn->count;
 	tau=moldyn->tau;
 	tau_square=moldyn->tau_square;
 
-#ifdef CONSTRAINT_110_5832
-	if(count==5833) {
-		atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
-		atom[5832].f.y=-atom[5832].f.x;
-	}
-#endif
-#ifdef CONSTRAINT_11X_5832
-	if(count==5833) {
-		// second trafo
-		xt=atom[5832].f.x;
-		yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
-		zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
-		// first trafo
-		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
-		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
-		ztt=zt;
-		// apply constraints
-		ytt=0.0;
-		// first trafo backwards
-		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		zt=ztt;
-		// second trafo backwards
-		atom[5832].f.x=xt;
-		atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-		atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-	}
-#endif
 	for(i=0;i<count;i++) {
+
 		/* check whether fixed atom */
 		if(atom[i].attr&ATOM_ATTR_FP)
 			continue;
+
 		/* new positions */
 		h=0.5/atom[i].mass;
-		v3_scale(&delta,&(atom[i].v),tau);
-#ifdef CONSTRAINT_110_5832
-		if(i==5832) {
-			delta.y=-delta.x;
+
+		/* constraint relaxation */
+		if(crtt) {
+			// forces
+			basis_trafo(&(atom[i].f),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].f.x=0;
+			if(constraints[3*i+1])
+				atom[i].f.y=0;
+			if(constraints[3*i+2])
+				atom[i].f.z=0;
+			basis_trafo(&(atom[i].f),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			// velocities
+			basis_trafo(&(atom[i].v),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].v.x=0;
+			if(constraints[3*i+1])
+				atom[i].v.y=0;
+			if(constraints[3*i+2])
+				atom[i].v.z=0;
+			basis_trafo(&(atom[i].v),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
 		}
-#endif
-#ifdef JHDVHJDV
-#ifdef CONSTRAINT_11X_5832
-	if(i==5833) {
-		// second trafo
-		xt=delta.x;
-		yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
-		zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
-		// first trafo
-		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
-		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
-		ztt=zt;
-		// apply constraints
-		ytt=0.0;
-		// first trafo backwards
-		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		zt=ztt;
-		// second trafo backwards
-		delta.x=xt;
-		delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-		delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-	}
-#endif
-#endif
+
 #ifndef QUENCH
+		v3_scale(&delta,&(atom[i].v),tau);
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 #endif
 		v3_scale(&delta,&(atom[i].f),h*tau_square);
-#ifdef CONSTRAINT_110_5832
-		if(i==5832) {
-			delta.y=-delta.x;
-		}
-#endif
-#ifdef GHDJHBSJBSD
-#ifdef CONSTRAINT_11X_5832
-	if(i==5833) {
-		// second trafo
-		xt=delta.x;
-		yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
-		zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
-		// first trafo
-		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
-		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
-		ztt=zt;
-		// apply constraints
-		ytt=0.0;
-		// first trafo backwards
-		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		zt=ztt;
-		// second trafo backwards
-		delta.x=xt;
-		delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-		delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-	}
-#endif
-#endif
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 		//check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
 		check_per_bound(moldyn,&(atom[i].r));
@@ -2408,38 +2385,37 @@ int velocity_verlet(t_moldyn *moldyn) {
 	albe_potential_force_calc(moldyn);
 #endif
 
-#ifdef CONSTRAINT_110_5832
-	if(count==5833) {
-		atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
-		atom[5832].f.y=-atom[5832].f.x;
-	}
-#endif
-#ifdef CONSTRAINT_11X_5832
-	if(count==5833) {
-		// second trafo
-		xt=atom[5832].f.x;
-		yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
-		zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
-		// first trafo
-		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
-		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
-		ztt=zt;
-		// apply constraints
-		ytt=0.0;
-		// first trafo backwards
-		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
-		zt=ztt;
-		// second trafo backwards
-		atom[5832].f.x=xt;
-		atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-		atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-	}
-#endif
 	for(i=0;i<count;i++) {
 		/* check whether fixed atom */
 		if(atom[i].attr&ATOM_ATTR_FP)
 			continue;
+
+		/* constraint relaxation */
+		if(crtt) {
+			// forces
+			basis_trafo(&(atom[i].f),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].f.x=0;
+			if(constraints[3*i+1])
+				atom[i].f.y=0;
+			if(constraints[3*i+2])
+				atom[i].f.z=0;
+			basis_trafo(&(atom[i].f),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			// velocities
+			basis_trafo(&(atom[i].v),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].v.x=0;
+			if(constraints[3*i+1])
+				atom[i].v.y=0;
+			if(constraints[3*i+2])
+				atom[i].v.z=0;
+			basis_trafo(&(atom[i].v),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+		}
+
 		/* again velocities [actually v(t+tau)] */
 		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -2522,8 +2498,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		/* single particle potential/force */
 		if(itom[i].attr&ATOM_ATTR_1BP)
-			if(moldyn->func1b)
-				moldyn->func1b(moldyn,&(itom[i]));
+			if(moldyn->func_i0)
+				moldyn->func_i0(moldyn,&(itom[i]));
 
 		if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
 			continue;
@@ -2538,8 +2514,14 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		dnlc=lc->dnlc;
 
+#ifndef STATIC_LISTS
+		/* check for later 3 body interaction */
+		if(itom[i].attr&ATOM_ATTR_3BP)
+			memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
+
 		/* first loop over atoms j */
-		if(moldyn->func2b) {
+		if(moldyn->func_j0) {
 			for(j=0;j<27;j++) {
 
 				bc_ij=(j<dnlc)?0:1;
@@ -2571,12 +2553,69 @@ 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,
-						               &(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
 					}
 #ifdef STATIC_LISTS
 				}
@@ -2589,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 			}
 		}
 
-		/* 3 body potential/force */
+		/* continued 3 body potential/force */
 
 		if(!(itom[i].attr&ATOM_ATTR_3BP))
 			continue;
@@ -2642,18 +2681,18 @@ int potential_force_calc(t_moldyn *moldyn) {
 				/* reset 3bp run */
 				moldyn->run3bp=1;
 
-				if(moldyn->func3b_j1)
-					moldyn->func3b_j1(moldyn,
-					                  &(itom[i]),
-					                  jtom,
-					                  bc_ij);
+				if(moldyn->func_j1)
+					moldyn->func_j1(moldyn,
+					                &(itom[i]),
+					                jtom,
+					                bc_ij);
 
-				/* in first j loop, 3bp run can be skipped */
+				/* in j loop, 3bp run can be skipped */
 				if(!(moldyn->run3bp))
 					continue;
 			
-				/* first loop over atoms k */
-				if(moldyn->func3b_k1) {
+				/* first loop over atoms k in second j loop */
+				if(moldyn->func_j1_k0) {
 
 				for(k=0;k<27;k++) {
 
@@ -2686,8 +2725,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
 
-						if(ktom==jtom)
-							continue;
+						//if(ktom==jtom)
+						//	continue;
 
 						if(ktom==&(itom[i]))
 							continue;
@@ -2711,14 +2750,15 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				}
 
-				if(moldyn->func3b_j2)
-					moldyn->func3b_j2(moldyn,
-					                  &(itom[i]),
-					                  jtom,
-					                  bc_ij);
+				/* continued j loop after first k loop */
+				if(moldyn->func_j1c)
+					moldyn->func_j1c(moldyn,
+					                 &(itom[i]),
+					                 jtom,
+					                 bc_ij);
 
 				/* second loop over atoms k */
-				if(moldyn->func3b_k2) {
+				if(moldyn->func_j1_k1) {
 
 				for(k=0;k<27;k++) {
 
@@ -2751,17 +2791,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
 
-						if(ktom==jtom)
-							continue;
+						//if(ktom==jtom)
+						//	continue;
 
 						if(ktom==&(itom[i]))
 							continue;
 
-						moldyn->func3b_k2(moldyn,
-						                  &(itom[i]),
-						                  jtom,
-						                  ktom,
-						                  bc_ik|bc_ij);
+						moldyn->func_j1_k1(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
 
 #ifdef STATIC_LISTS
 					}
@@ -2776,11 +2816,11 @@ int potential_force_calc(t_moldyn *moldyn) {
 				
 				}
 
-				/* 2bp post function */
-				if(moldyn->func3b_j3) {
-					moldyn->func3b_j3(moldyn,
-					                  &(itom[i]),
-					                  jtom,bc_ij);
+				/* finish of j loop after second k loop */
+				if(moldyn->func_j1e) {
+					moldyn->func_j1e(moldyn,
+					                 &(itom[i]),
+					                 jtom,bc_ij);
 				}
 #ifdef STATIC_LISTS
 			}
@@ -2791,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
 		
 		}
-		
+
 #ifdef DEBUG
 	//printf("\n\n");
 #endif