basic p control added, virial still needed!

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 73becde..b6fa8d3 100644 (file)
--- 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",