added insertion type controled by config.h

[physik/posic.git] / sic.c

diff --git a/sic.c b/sic.c
index e1103bd..cea4c39 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);
 
@@ -422,6 +440,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
         *