X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=a6f5c4557efc68b4ff01e6996792b782b5d2bcea;hb=ff44bcba9df0d021568abcb553d1d490c442f696;hp=bc4fdffdee29be7341c67692f1a3282bc049927b;hpb=dfb47bde6bf974d1b231048ec8bcda712404ccc1;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index bc4fdff..a6f5c45 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,6 +83,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 	pthread_mutex_init(&emutex,NULL);
 #endif
 
+	if(crtt)
+		printf("USING CRT\n");
+
 	return 0;
 }
 
@@ -258,7 +266,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
 	switch(type) {
 		case MOLDYN_POTENTIAL_TM:
-			moldyn->func1b=tersoff_mult_1bp;
+			//moldyn->func1b=tersoff_mult_1bp;
 			moldyn->func3b_j1=tersoff_mult_3bp_j1;
 			moldyn->func3b_k1=tersoff_mult_3bp_k1;
 			moldyn->func3b_j2=tersoff_mult_3bp_j2;
@@ -516,7 +524,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;
@@ -588,6 +597,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");
@@ -595,6 +606,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");
@@ -602,6 +615,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");
@@ -736,6 +751,9 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
 	t_atom *new,*old;
 	int cnt;
+#if defined LOWMEM_LISTS || defined PTHREADS
+	void *ptr;
+#endif
 
 	old=moldyn->atom;
 
@@ -758,6 +776,28 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
 	free(old);
 
+#ifdef LOWMEM_LISTS
+	ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+	if(!ptr) {
+		perror("[moldyn] list realloc (del atom)");
+		return -1;
+	}
+	moldyn->lc.subcell->list=ptr;
+	// update lists
+	link_cell_update(moldyn);
+#endif
+
+#ifdef PTHREADS
+	ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t));
+	if(!ptr) {
+		perror("[moldyn] mutex realloc (add atom)");
+		return -1;
+	}
+	amutex=ptr;
+	pthread_mutex_destroy(&(amutex[moldyn->count+1]));
+#endif
+
+
 	return 0;
 }
 
@@ -2046,11 +2086,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
-		/*
+#ifdef PDEBUG	
 		thermodynamic_pressure_calc(moldyn);
 		printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
 		       moldyn->tp/BAR);
-		*/
+#endif
 
 		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
@@ -2065,10 +2105,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(e) {
 			if(!(moldyn->total_steps%e))
 				dprintf(moldyn->efd,
-				        "%f %f %f %f\n",
+				        "%f %f %f %f %f %f\n",
 				        moldyn->time,moldyn->ekin/energy_scale,
 				        moldyn->energy/energy_scale,
-				        get_total_energy(moldyn)/energy_scale);
+				        get_total_energy(moldyn)/energy_scale,
+				        moldyn->ekin/EV,moldyn->energy/EV);
 		}
 		if(m) {
 			if(!(moldyn->total_steps%m)) {
@@ -2156,15 +2197,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 
 		/* display progress */
+#ifndef PDEBUG
 		if(!(i%10)) {
+#endif
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
 printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
        moldyn->p/BAR,moldyn->p_avg/BAR,
-       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+       moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2172,7 +2218,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
 			/* copy over time */
 			t1=t2;
+#ifndef PDEBUG
 		}
+#endif
 
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
@@ -2193,6 +2241,59 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
 	}
 
+	/* writing a final save file! */
+	if(s) {
+		snprintf(dir,128,"%s/s-final.save",moldyn->vlsdir);
+		fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,S_IRUSR|S_IWUSR);
+		if(fd<0) perror("[moldyn] save fd open");
+		else {
+			write(fd,moldyn,sizeof(t_moldyn));
+			write(fd,moldyn->atom,
+			      moldyn->count*sizeof(t_atom));
+		}
+		close(fd);
+	}
+	/* writing a final visual file! */
+	if(a)
+		visual_atoms(moldyn);
+
+	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;
 }
 
@@ -2211,15 +2312,35 @@ int velocity_verlet(t_moldyn *moldyn) {
 	tau_square=moldyn->tau_square;
 
 	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;
+
+		/* constraint relaxation */
+		if(crtt) {
+			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]);
+		}
+
+#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);
 		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));
 
 		/* velocities [actually v(t+tau/2)] */
@@ -2235,7 +2356,11 @@ int velocity_verlet(t_moldyn *moldyn) {
 
 	/* forces depending on chosen potential */
 #ifndef ALBE_FAST
-	potential_force_calc(moldyn);
+	// if albe, use albe force calc routine
+	//if(moldyn->func3b_j1==albe_mult_3bp_j1)
+	//	albe_potential_force_calc(moldyn);
+	//else
+		potential_force_calc(moldyn);
 #else
 	albe_potential_force_calc(moldyn);
 #endif
@@ -2244,6 +2369,21 @@ int velocity_verlet(t_moldyn *moldyn) {
 		/* check whether fixed atom */
 		if(atom[i].attr&ATOM_ATTR_FP)
 			continue;
+
+		/* constraint relaxation */
+		if(crtt) {
+			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]);
+		}
+
 		/* 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);
@@ -2687,6 +2827,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 	return 0;
 }
         
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+	
+	double x,y,z;
+	t_3dvec *dim;
+
+	dim=&(moldyn->dim);
+
+	x=dim->x/2;
+	y=dim->y/2;
+	z=dim->z/2;
+
+	if(moldyn->status&MOLDYN_STAT_PBX) {
+		if(a->r.x>=x) {
+			a->pbc[0]+=1;
+			a->r.x-=dim->x;
+		}
+		else if(-a->r.x>x) {
+			a->pbc[0]-=1;
+			a->r.x+=dim->x;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBY) {
+		if(a->r.y>=y) {
+			a->pbc[1]+=1;
+			a->r.y-=dim->y;
+		}
+		else if(-a->r.y>y) {
+			a->pbc[1]-=1;
+			a->r.y+=dim->y;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBZ) {
+		if(a->r.z>=z) {
+			a->pbc[2]+=1;
+			a->r.z-=dim->z;
+		}
+		else if(-a->r.z>z) {
+			a->pbc[2]-=1;
+			a->r.z+=dim->z;
+		}
+	}
+
+	return 0;
+}
+        
 /*
  * debugging / critical check functions
  */
@@ -2777,7 +2962,7 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 	if(fsize!=sizeof(t_moldyn)+size) {
 		corr=fsize-sizeof(t_moldyn)-size;
 		printf("[moldyn] WARNING: lsf (illegal file size)\n");
-		printf("  moifying offset:\n");
+		printf("  modifying offset:\n");
 		printf("  - current pos: %d\n",sizeof(t_moldyn));
 		printf("  - atom size: %d\n",size);
 		printf("  - file size: %d\n",fsize);
@@ -3022,6 +3207,7 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	int i;
 	t_atom *atom;
 	t_3dvec dist;
+	t_3dvec final_r;
 	double d2;
 	int a_cnt;
 	int b_cnt;
@@ -3035,8 +3221,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
 	for(i=0;i<moldyn->count;i++) {
 
-		v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
-		check_per_bound(moldyn,&dist);
+		/* care for pb crossing */
+		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) {
@@ -3058,6 +3248,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;
@@ -3203,9 +3444,8 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	ba=data;
 
 	/* increase total bond counter
-	 * ... double counting!
 	 */
-	ba->tcnt+=2;
+	ba->tcnt+=1;
 
 	if(itom->brand==0)
 		ba->acnt[jtom->tag]+=1;
@@ -3222,10 +3462,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 
 int bond_analyze(t_moldyn *moldyn,double *quality) {
 
-	// by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
-
-	int qcnt;
-	int ccnt,cset;
+	int qcnt4;
+	int qcnt3;
+	int ccnt4;
+	int ccnt3;
+	int bcnt;
 	t_ba ba;
 	int i;
 	t_atom *atom;
@@ -3245,9 +3486,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
 	ba.tcnt=0;
-	qcnt=0;
-	ccnt=0;
-	cset=0;
+	qcnt4=0; qcnt3=0;
+	ccnt4=0; ccnt3=0;
+	bcnt=0;
 
 	atom=moldyn->atom;
 
@@ -3256,27 +3497,30 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	for(i=0;i<moldyn->count;i++) {
 		if(atom[i].brand==0) {
 			if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
-				qcnt+=4;
+				qcnt4+=4;
+			if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
+				qcnt3+=3;
 		}
 		else {
 			if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
-				qcnt+=4;
-				ccnt+=1;
+				qcnt4+=4;
+				ccnt4+=1;
+			}
+			if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
+				qcnt3+=4;
+				ccnt3+=1;
 			}
-			cset+=1;
+			bcnt+=1;
 		}
 	}
 
-	printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
-	printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
-
 	if(quality) {
-		quality[0]=1.0*ccnt/cset;
-		quality[1]=1.0*qcnt/ba.tcnt;
+		quality[0]=1.0*ccnt4/bcnt;
+		quality[1]=1.0*ccnt3/bcnt;
 	}
 	else {
-		printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
-		printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
+		printf("[moldyn] bond analyze: %f %f\n",
+		       1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
 	}
 
 	return 0;