X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=bd29afc9a3e6853bded38acae91dfafdf15a662a;hp=15dcedef5fd105b05ce69e2c9d49ad5a8c4c353e;hb=8524173a28f2c22a539ef1b0910a1136d9cb254b;hpb=f0e4c6422ec1aff0cf86597fef919335bba75c1b

diff --git a/moldyn.c b/moldyn.c
index 15dcede..bd29afc 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -81,6 +81,10 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 #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;
 }
@@ -2264,6 +2268,10 @@ 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;
@@ -2275,6 +2283,28 @@ int velocity_verlet(t_moldyn *moldyn) {
 		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 */
@@ -2288,6 +2318,30 @@ int velocity_verlet(t_moldyn *moldyn) {
 			delta.y=-delta.x;
 		}
 #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_add(&(atom[i].r),&(atom[i].r),&delta);
 #endif
@@ -2296,6 +2350,30 @@ int velocity_verlet(t_moldyn *moldyn) {
 		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]));
@@ -2328,6 +2406,28 @@ int velocity_verlet(t_moldyn *moldyn) {
 		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 */
@@ -3197,6 +3297,57 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	return 0;
 }
 
+int calculate_msd(t_moldyn *moldyn,double *msd) {
+
+	int i;
+	t_atom *atom;
+	t_3dvec dist;
+	t_3dvec final_r;
+	double d2;
+	int a_cnt;
+	int b_cnt;
+
+	atom=moldyn->atom;
+	msd[0]=0;
+	msd[1]=0;
+	msd[2]=0;
+	a_cnt=0;
+	b_cnt=0;
+
+	for(i=0;i<moldyn->count;i++) {
+
+		/* care for pb crossing */
+		if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) {
+			printf("[moldyn] msd pb crossings for atom %d\n",i);
+			printf("  x: %d y: %d z: %d\n",
+			       atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]);
+		}
+		final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+		final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+		final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+		/* calculate distance */
+		v3_sub(&dist,&final_r,&(atom[i].r_0));
+		d2=v3_absolute_square(&dist);
+
+		if(atom[i].brand) {
+			b_cnt+=1;
+			msd[1]+=d2;
+		}
+		else {
+			a_cnt+=1;
+			msd[0]+=d2;
+		}
+
+		msd[2]+=d2;
+	}
+
+	msd[0]/=a_cnt;
+	msd[1]/=b_cnt;
+	msd[2]/=moldyn->count;
+		
+	return 0;
+}
+
 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
 
 	return 0;