no DEBUG by default

[physik/posic.git] / sic.c

diff --git a/sic.c b/sic.c
index e1103bd..406ad53 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -40,7 +40,7 @@ int insert_atoms(t_moldyn *moldyn) {
        int i,j;
        u8 run;
        t_3dvec r,v,dist;
-       double d;
+       double d,dmin;
 
        t_atom *atom;
 
@@ -51,38 +51,39 @@ int insert_atoms(t_moldyn *moldyn) {
        for(j=0;j<INS_ATOMS;j++) {
                run=1;
                while(run) {
+#ifdef INS_TETRA
                        // tetrahedral
-                       /*
                        r.x=0.0;
                        r.y=0.0;
                        r.z=0.0;
-                       */
+#endif
+#ifdef INS_HEXA
                        // hexagonal
-                       /*
                        r.x=-1.0/8.0*ALBE_LC_SI;
                        r.y=-1.0/8.0*ALBE_LC_SI;
                        r.z=1.0/8.0*ALBE_LC_SI;
-                       */
+#endif
+#ifdef INS_110DB
                        // 110 dumbbell
-                       /*
                        r.x=(-0.5+0.25+0.125)*ALBE_LC_SI;
                        r.y=(-0.5+0.25+0.125)*ALBE_LC_SI;
                        r.z=(-0.5+0.25)*ALBE_LC_SI;
                        md->atom[4372].r.x=(-0.5+0.125+0.125)*ALBE_LC_SI;
                        md->atom[4372].r.y=(-0.5+0.125+0.125)*ALBE_LC_SI;
-                       */
+#endif
+#ifdef INS_RAND
                        // random
-                       //
                        r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX;
                        r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY;
                        r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ;
-                       //
+#endif
                        // offset
                        r.x+=INS_OFFSET;
                        r.y+=INS_OFFSET;
                        r.z+=INS_OFFSET;
                        /* assume valid coordinates */
                        run=0;
+                       dmin=10000000000.0;             // for sure too high!
                        for(i=0;i<moldyn->count;i++) {
                                atom=&(moldyn->atom[i]);
                                v3_sub(&dist,&(atom->r),&r);
@@ -94,13 +95,16 @@ int insert_atoms(t_moldyn *moldyn) {
                                        run=1;
                                        break;
                                }
+                               if(d<dmin)
+                                       dmin=d;
                        }
                }
                add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
                         ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
-                        //ATOM_ATTR_HB|ATOM_ATTR_VB,
-                        ATOM_ATTR_HB,
+                        INS_ATTR,
                         &r,&v);
+               printf(" %02d: atom %d | %f %f %f | %f\n",
+                      j,moldyn->count-1,r.x,r.y,r.z,dmin);
        }
 
        return 0;
@@ -112,10 +116,14 @@ int sic_hook(void *moldyn,void *hook_params) {
        t_moldyn *md;
        int steps;
        double tau;
+       double dt;
 
        hp=hook_params;
        md=moldyn;
 
+       tau=1.0;
+       steps=0;
+
        /* switch on t scaling */
        if(md->schedule.count==0)
                set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
@@ -129,12 +137,15 @@ int sic_hook(void *moldyn,void *hook_params) {
        /* act according to state */
        switch(hp->state) {
                case STATE_INSERT:
+                       /* assigne values */
+                       steps=INS_RELAX;
+                       tau=INS_TAU;
                        /* check temperature */
-                       if(md->t_avg-md->t_ref>INS_DELTA_TC) {
-                               steps=INS_RELAX;
-                               tau=INS_TAU;
+                       dt=md->t_avg-md->t_ref;
+                       if(dt<0)
+                               dt=-dt;
+                       if(dt>INS_DELTA_TC)
                                break;
-                       }
                        /* insert atoms */
                        hp->insert_count+=1;
                        printf("   ### insert atoms (%d/%d) ###\n",
@@ -145,11 +156,15 @@ int sic_hook(void *moldyn,void *hook_params) {
                                hp->state=STATE_POSTRUN;
                        break;
                case STATE_POSTRUN:
-                       /* settings */
-                       if(md->t-md->t_ref>POST_DELTA_TC) {
-                               steps=POST_RELAX;
-                               tau=POST_TAU;
-                       }
+                       /* assigne values */
+                       steps=POST_RELAX;
+                       tau=POST_TAU;
+                       /* check temperature */
+                       dt=md->t_avg-md->t_ref;
+                       if(dt<0)
+                               dt=-dt;
+                       if(dt>POST_DELTA_TC)
+                               break;
                        /* decrease temperature */
                        hp->postrun_count+=1;
                        printf(" ### postrun (%d/%d) ###\n",
@@ -163,6 +178,9 @@ int sic_hook(void *moldyn,void *hook_params) {
                        break;
        }
 
+       /* reset the average counters */
+       average_reset(md);
+
        /* add schedule */
        moldyn_add_schedule(md,steps,tau);
 
@@ -301,34 +319,60 @@ int main(int argc,char **argv) {
 
        /* set (initial) dimensions of simulation volume */
 #ifdef ALBE
+ #ifdef INIT_SI
        set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE);
-       //set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
-       //set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
+ #endif
+ #ifdef INIT_C
+       set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
+ #endif
+ #ifdef INIT_3CSIC
+       set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
+ #endif
 #else
+ #ifdef INIT_SI
        set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE);
-       //set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
-       //set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
+ #endif
+ #ifdef INIT_C
+       set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
+ #endif
+ #ifdef INIT_3CSIC
+       set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
+ #endif
 #endif
 
        /* set periodic boundary conditions in all directions */
        set_pbc(&md,TRUE,TRUE,TRUE);
 
        /* create the lattice / place atoms */
-       //
+
+       // diamond
 #ifdef ALBE
+ #ifdef INIT_SI
        create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
-       //create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      0,LCNTX,LCNTY,LCNTZ,NULL);
+ #endif
+ #ifdef INIT_C
+       create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      1,LCNTX,LCNTY,LCNTZ,NULL);
+ #endif
 #else
+ #ifdef INIT_SI
        create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-#endif
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-       //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,LCNTX,LCNTY,LCNTZ,NULL);
-       //               1,LCNTX,LCNTY,LCNTZ,NULL);
+                      0,LCNTX,LCNTY,LCNTZ,NULL);
+ #endif
+ #ifdef INIT_C
+       create_lattice(&md,DIAMOND,LC_C,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      1,LCNTX,LCNTY,LCNTZ,NULL);
+ #endif
+#endif
 
-       /* create zinkblende structure */
-       /*
-#ifdef ALBE
+       // zinkblende 
+#ifdef INIT_3CSIC
+ #ifdef ALBE
        r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
        create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
@@ -337,7 +381,7 @@ int main(int argc,char **argv) {
        create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
                       1,LCNTX,LCNTY,LCNTZ,&r);
-#else
+ #else
        r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
        create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
@@ -346,8 +390,8 @@ int main(int argc,char **argv) {
        create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
                       1,LCNTX,LCNTY,LCNTZ,&r);
+ #endif
 #endif
-       */
 
        /* check for right atom placing */
        moldyn_bc_check(&md);
@@ -422,6 +466,9 @@ int main(int argc,char **argv) {
        moldyn_set_log(&md,SAVE_STEP,LOG_S);
        moldyn_set_log(&md,CREATE_REPORT,0);
 
+       /* next neighbour distance for critical checking */
+       set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0));
+
        /*
         * let's do the actual md algorithm now
         *