X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=33224a69f82fd62c175ea6e300f757f70b8851bb;hb=98df942a2412e6ace72ae221f2e967f24bba7108;hp=bd29afc9a3e6853bded38acae91dfafdf15a662a;hpb=cba0d347201db1c99dc23736658ed78ed12ad5bd;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index bd29afc..33224a6 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -55,6 +55,11 @@ pthread_mutex_t *amutex;
 pthread_mutex_t emutex;
 #endif
 
+/* fully constrained relaxation technique - global pointers */
+u8 crtt;
+u8 *constraints;
+double *trafo_angle;
+
 /*
  * the moldyn functions
  */
@@ -78,14 +83,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;
 }
 
@@ -2260,6 +2257,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) {
@@ -2268,113 +2302,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));
@@ -2401,38 +2374,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);