X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=b6fa8d3df41dbf12e89c0596463a6cba4207255b;hb=88ea16029cc07a837977f737f4536fa3e881dcee;hp=73becdeed1a9b27220a714d55eed845d2ec71db2;hpb=72df64eacc634e315f2205fc7bd2406223f92bf3;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 73becde..b6fa8d3 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -78,6 +78,13 @@ int set_temperature(t_moldyn *moldyn,double t_ref) {
 	return 0;
 }
 
+int set_pressure(t_moldyn *moldyn,double p_ref) {
+
+	moldyn->p_ref=p_ref;
+
+	return 0;
+}
+
 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
 
 	moldyn->pt_scale=(ptype|ttype);
@@ -93,16 +100,19 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 	moldyn->dim.y=y;
 	moldyn->dim.z=z;
 
+	moldyn->volume=x*y*z;
+
 	if(visualize) {
 		moldyn->vis.dim.x=x;
 		moldyn->vis.dim.y=y;
 		moldyn->vis.dim.z=z;
 	}
 
-	printf("[moldyn] dimensions in A:\n");
+	printf("[moldyn] dimensions in A and A^2 respectively:\n");
 	printf("  x: %f\n",moldyn->dim.x);
 	printf("  y: %f\n",moldyn->dim.y);
 	printf("  z: %f\n",moldyn->dim.z);
+	printf("  volume: %f\n",moldyn->volume);
 	printf("  visualize simulation box: %s\n",visualize?"on":"off");
 
 	return 0;
@@ -390,7 +400,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 			count+=1;
 		}
 	}
-	if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN);
+	if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 	else return 0;	/* no atoms involved in scaling! */
 	
 	/* (temporary) hack for e,t = 0 */
@@ -423,6 +433,57 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 	return 0;
 }
 
+int scale_volume(t_moldyn *moldyn) {
+
+	t_atom *atom;
+	t_3dvec *dim,*vdim;
+	double virial,scale;
+	t_linkcell *lc;
+	int i;
+
+	atom=moldyn->atom;
+	dim=&(moldyn->dim);
+	vdim=&(moldyn->vis.dim);
+	lc=&(moldyn->lc);
+
+	for(i=0;i<moldyn->count;i++)
+		virial+=v3_norm(&(atom[i].virial));
+
+printf("%f\n",virial);
+	/* get pressure from virial */
+	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+	moldyn->p/=moldyn->volume;
+printf("%f\n",moldyn->p/(ATM));
+
+	/* scale factor */
+	if(moldyn->pt_scale&P_SCALE_BERENDSEN)
+		scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc);
+	else 
+		/* should actually never be used */
+		scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0);
+
+printf("scale = %f\n",scale);
+	/* actual scaling */
+	dim->x*=scale;
+	dim->y*=scale;
+	dim->z*=scale;
+	if(vdim->x) vdim->x=dim->x;
+	if(vdim->y) vdim->y=dim->y;
+	if(vdim->z) vdim->z=dim->z;
+	moldyn->volume*=(scale*scale*scale);
+
+	/* check whether we need a new linkcell init */
+	if((dim->x/moldyn->cutoff!=lc->nx)||
+	   (dim->y/moldyn->cutoff!=lc->ny)||
+	   (dim->z/moldyn->cutoff!=lc->nx)) {
+		link_cell_shutdown(moldyn);
+		link_cell_init(moldyn);
+	}
+
+	return 0;
+
+}
+
 double get_e_kin(t_moldyn *moldyn) {
 
 	int i;
@@ -724,6 +785,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		/* p/t scaling */
 		if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
 			scale_velocity(moldyn,FALSE);
+		if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
+{
+printf("going to do p scale ...\n");
+			scale_volume(moldyn);
+printf("done\n");
+}
 
 		/* check for log & visualization */
 		if(e) {
@@ -853,13 +920,16 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 	/* reset energy */
 	moldyn->energy=0.0;
-
+	
 	/* get energy and force of every atom */
 	for(i=0;i<count;i++) {
 
 		/* reset force */
 		v3_zero(&(itom[i].f));
 
+		/* reset viral of atom i */
+		v3_zero(&(itom[i].virial));
+
 		/* single particle potential/force */
 		if(itom[i].attr&ATOM_ATTR_1BP)
 			moldyn->func1b(moldyn,&(itom[i]));
@@ -1627,6 +1697,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
 		v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
 		                                  /* scaled with 0.5 ^ */
+
 	}
 
 	return 0;
@@ -1655,6 +1726,7 @@ x=dim->x/2;
 			printf("FATAL: atom %d: x: %.20f (%.20f)\n",
 			       i,atom[i].r.x,dim->x/2);
 			printf("diagnostic:\n");
+			printf("-----------\natom.r.x:\n");
 			for(j=0;j<8;j++) {
 				memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
 				for(k=0;k<8;k++)
@@ -1662,7 +1734,7 @@ x=dim->x/2;
 					((byte)&(1<<k))?1:0,
 					(k==7)?'\n':'|');
 			}
-			printf("---------------\n");
+			printf("---------------\nx=dim.x/2:\n");
 			for(j=0;j<8;j++) {
 				memcpy(&byte,(u8 *)(&x)+j,1);
 				for(k=0;k<8;k++)
@@ -1670,8 +1742,8 @@ x=dim->x/2;
 					((byte)&(1<<k))?1:0,
 					(k==7)?'\n':'|');
 			}
-			if(atom[i].r.x==x) printf("hier gleich!\n");
-			else printf("hier NICHT gleich!\n");
+			if(atom[i].r.x==x) printf("the same!\n");
+			else printf("different!\n");
 		}
 		if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
 			printf("FATAL: atom %d: y: %.20f (%.20f)\n",