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

diff --git a/moldyn.c b/moldyn.c
index 056047c..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,6 +56,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,10 +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
 	return 0;
 }
 
@@ -269,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:
@@ -520,7 +532,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   t_part_params *p_params,t_defect_params *d_params) {
+                   t_part_params *p_params,t_defect_params *d_params,
+                   t_offset_params *o_params) {
 
 	int new,count;
 	int ret;
@@ -592,6 +605,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 
 	switch(type) {
 		case CUBIC:
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,lc);
 			ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"cubic");
@@ -599,6 +614,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 		case FCC:
 			if(!origin)
 				v3_scale(&orig,&orig,0.5);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 			ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"fcc");
@@ -606,6 +623,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 		case DIAMOND:
 			if(!origin)
 				v3_scale(&orig,&orig,0.25);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 			ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"diamond");
@@ -806,8 +825,6 @@ int del_atom(t_moldyn *moldyn,int tag) {
 			case DEFECT_STYPE_DB_Z:\
 				d_o.z=d_params->od;\
 				d_d.z=d_params->dd;\
-d_d.x=0.9;\
-d_d.y=0.9;\
 				break;\
 			case DEFECT_STYPE_DB_R:\
 				break;\
@@ -2251,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) {
@@ -2265,33 +2319,46 @@ int velocity_verlet(t_moldyn *moldyn) {
 	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
 	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
+
 #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
 		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));
@@ -2318,16 +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
 	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);
@@ -2410,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;
@@ -2426,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;
@@ -2459,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
 				}
@@ -2477,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;
@@ -2530,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++) {
 
@@ -2574,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;
@@ -2599,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++) {
 
@@ -2639,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
 					}
@@ -2664,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
 			}
@@ -2679,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
 		
 		}
-		
+
 #ifdef DEBUG
 	//printf("\n\n");
 #endif
@@ -3192,6 +3344,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;