basic p control added, virial still needed!

diff --git a/moldyn.c b/moldyn.c
index 73becdeed1a9b27220a714d55eed845d2ec71db2..b6fa8d3df41dbf12e89c0596463a6cba4207255b 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",

diff --git a/moldyn.h b/moldyn.h
index 72ff30bdfb00ab876cce6261360b827298c77043..4203c79e4cd92ffc7cf1da846a1c7d589cd02a46 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -27,6 +27,7 @@ typedef struct s_atom {
        t_3dvec r;              /* position */
        t_3dvec v;              /* velocity */
        t_3dvec f;              /* force */
+       t_3dvec virial;         /* virial */
        int element;            /* number of element in pse */
        double mass;            /* atom mass */
        u8 bnum;                /* brand number */
@@ -67,6 +68,7 @@ typedef struct s_moldyn {
        t_atom *atom;           /* pointer to the atoms */
 
        t_3dvec dim;            /* dimensions of the simulation volume */
+       double volume;          /* volume of sim cell (dim.x*dim.y*dim.z) */
 
        /* potential force function and parameter pointers */
        int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
@@ -92,7 +94,7 @@ typedef struct s_moldyn {
        double p;               /* actual pressure */
 
        /* pressure and temperature control (velocity/volume scaling) */
-       /* (in units of tau!) */
+       /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */
        unsigned char pt_scale; /* type of p and t scaling */
        double t_tc;            /* t berendsen control time constant */
        double p_tc;            /* p berendsen control time constant */
@@ -129,13 +131,15 @@ typedef struct s_moldyn {
        int debug;              /* debugging stuff, ignore */
 } t_moldyn;
 
-#define MOLDYN_STAT_PBX                        0x08    /* periodic boudaries in x */
-#define MOLDYN_STAT_PBY                        0x10    /* y */
-#define MOLDYN_STAT_PBZ                        0x20    /* and z direction */
+#define MOLDYN_STAT_PBX                        0x01    /* periodic boudaries in x */
+#define MOLDYN_STAT_PBY                        0x02    /* y */
+#define MOLDYN_STAT_PBZ                        0x04    /* and z direction */
 
-#define MOLDYN_1BP                     0x00    /* care about single */
-#define MOLDYN_2BP                     0x01    /* 2 body */
-#define MOLDYN_3BP                     0x02    /* and 3 body particle pots */
+#define MOLDYN_PSCALE                  0x08    /* size controlled by piston */
+
+#define MOLDYN_1BP                     0x10    /* care about single */
+#define MOLDYN_2BP                     0x20    /* 2 body */
+#define MOLDYN_3BP                     0x40    /* and 3 body particle pots */
 
 #define T_SCALE_BERENDSEN              0x01    /* berendsen t control */
 #define T_SCALE_DIRECT                 0x02    /* direct t control */
@@ -246,12 +250,17 @@ typedef struct s_tersoff_mult_params {
  *
  */
 
+#define ONE_THIRD      (1.0/3.0)
+
 /*
  * default values
  *
  * - length unit: 1 A (1 A = 1e-10 m)
  * - time unit: 1 fs (1 fs = 1e-15 s)
  * - mass unit: 1 amu (1 amu = 1.6605388628e-27 kg )
+ *
+ * fyi: in the following 1 N = (amu*A)/(fs*fs)
+ *
  */
 
 #define METER                          1e10                    /* A */
@@ -259,6 +268,8 @@ typedef struct s_tersoff_mult_params {
 #define AMU                            1.6605388628e-27        /* kg */
 #define KILOGRAM                       (1.0/AMU)               /* amu */
 #define NEWTON (METER*KILOGRAM/(SECOND*SECOND))        /* A amu / fs^2 */
+#define PASCAL (NEWTON/(METER*METER))                  /* N / A^2 */
+#define ATM    (1.0133e5*PASCAL)                       /* N / A^2 */
 
 #define MOLDYN_TEMP                    273.0
 #define MOLDYN_TAU                     1.0
@@ -284,6 +295,7 @@ typedef struct s_tersoff_mult_params {
  *
  * phsical values / constants
  *
+ *
  */
 
 #define K_BOLTZMANN            (1.380650524e-23*METER*NEWTON)  /* NA/K */
@@ -349,6 +361,7 @@ int moldyn_shutdown(t_moldyn *moldyn);
 int set_int_alg(t_moldyn *moldyn,u8 algo);
 int set_cutoff(t_moldyn *moldyn,double cutoff);
 int set_temperature(t_moldyn *moldyn,double t_ref);
+int set_pressure(t_moldyn *moldyn,double p_ref);
 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
 int set_nn_dist(t_moldyn *moldyn,double dist);
@@ -370,6 +383,7 @@ int destroy_atoms(t_moldyn *moldyn);
 
 int thermal_init(t_moldyn *moldyn,u8 equi_init);
 int scale_velocity(t_moldyn *moldyn,u8 equi_init);
+int scale_volume(t_moldyn *moldyn);
 
 double get_e_kin(t_moldyn *moldyn);
 double get_e_pot(t_moldyn *moldyn);

diff --git a/sic.c b/sic.c
index bc10e314fe3b11121fe3b6b1b86a475cb13af6fe..516ee01ac59fffcf34482fafc97e408ec52de94a 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -111,6 +111,7 @@ int main(int argc,char **argv) {
        create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
                       0,5,5,5);
+       moldyn_bc_check(&md);
 
        /* testing configuration */
        //r.x=2.45/2;   v.x=0;
@@ -131,13 +132,20 @@ int main(int argc,char **argv) {
 
        /* set temperature */
        printf("[sic] setting temperature\n");
+       set_temperature(&md,273.0+1410.0);
        //set_temperature(&md,273.0+450.0);
-       set_temperature(&md,1.0);
        //set_temperature(&md,273.0);
+       //set_temperature(&md,1.0);
+       //set_temperature(&md,0.0);
+
+       /* set pressure */
+       printf("[sic] setting pressure\n");
+       set_pressure(&md,ATM);
 
        /* set p/t scaling */
        printf("[sic] set p/t scaling\n");
-       set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
+       //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
+       //                 T_SCALE_BERENDSEN,100.0);
        
        /* initial thermal fluctuations of particles (in equilibrium) */
        printf("[sic] thermal init\n");
@@ -145,13 +153,13 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,1000,1.0);
+       moldyn_add_schedule(&md,100001,1.0);
 
        /* activate logging */
        printf("[sic] activate logging\n");
        moldyn_set_log_dir(&md,"saves/test");
-       moldyn_set_log(&md,LOG_TOTAL_ENERGY,10);
-       moldyn_set_log(&md,VISUAL_STEP,10);
+       moldyn_set_log(&md,LOG_TOTAL_ENERGY,200);
+       moldyn_set_log(&md,VISUAL_STEP,200);
 
        /*
         * let's do the actual md algorithm now