scaling to fs, angstrom, amu done, probs with velocity scaling!

diff --git a/Makefile b/Makefile
index 742a7a5..6f3e39e 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,5 @@
 CC=gcc-3.4
+#CC=gcc
 CFLAGS=-Wall -O3 #-DSIMPLE_TESTING
 LDFLAGS=-lm
 

diff --git a/moldyn.c b/moldyn.c
index 8f6f262..9d9901c 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -385,12 +385,17 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
        else
                if(moldyn->pt_scale&T_SCALE_BERENDSEN)
                        scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc;
+printf("scale=%f\n",scale);
        scale=sqrt(scale);
+printf("debug: %f %f %f %f \n",scale,moldyn->t_ref,moldyn->t,moldyn->tau);
 
        /* velocity scaling */
-       for(i=0;i<moldyn->count;i++)
+       for(i=0;i<moldyn->count;i++) {
+printf("vorher: %f\n",atom[i].v.x);
                if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
                        v3_scale(&(atom[i].v),&(atom[i].v),scale);
+printf("nachher: %f\n",atom[i].v.x);
+       }
 
        return 0;
 }
@@ -447,12 +452,6 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 
        /* nn_dist is the nearest neighbour distance */
 
-       if(moldyn->t==5.0) {
-               printf("[moldyn] i do not estimate timesteps below %f K!\n",
-                      MOLDYN_CRITICAL_EST_TEMP);
-               return 23.42;
-       }
-
        tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 
        return tau;     
@@ -813,6 +812,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 /* generic potential and force calculation */
 
 int potential_force_calc(t_moldyn *moldyn) {
+printf("start pot force calc\n");
 
        int i,j,k,count;
        t_atom *itom,*jtom,*ktom;
@@ -833,6 +833,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
        /* get energy and force of every atom */
        for(i=0;i<count;i++) {
+printf("atom %d: %f\n",i,itom[i].r.x);
+printf("atom %d: %f\n",i,itom[i].v.x);
+printf("atom %d: %f\n",i,itom[i].f.x);
 
                /* reset force */
                v3_zero(&(itom[i].f));
@@ -927,9 +930,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                        }
                }
-//printf("DEBUG: %.15f \n",itom[i].f.x);
        }
 
+printf("end pot force calc\n");
        return 0;
 }
 
@@ -1050,8 +1053,8 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
        p->mu_m=0.5*(p->mu[0]+p->mu[1]);
 
        printf("[moldyn] tersoff mult parameter info:\n");
-       printf("  S (m)  | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
-       printf("  R (m)  | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
+       printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
+       printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
        printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
        printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
        printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],

diff --git a/moldyn.h b/moldyn.h
index 534287b..a2223c0 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -247,15 +247,25 @@ typedef struct s_tersoff_mult_params {
  *
  */
 
-/* default values */
+/*
+ * 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 )
+ */
+
+#define METER                          1e10                    /* A */
+#define SECOND                         1e15                    /* fs */
+#define AMU                            1.6605388628e-27        /* kg */
+#define KILOGRAM                       (1.0/AMU)               /* amu */
+#define NEWTON (METER*KILOGRAM/(SECOND*SECOND))        /* A amu / fs^2 */
 
 #define MOLDYN_TEMP                    273.0
-#define MOLDYN_TAU                     1.0e-15
-#define MOLDYN_CUTOFF                  1.0e-9
+#define MOLDYN_TAU                     1.0
+#define MOLDYN_CUTOFF                  10.0
 #define MOLDYN_RUNS                    1000000
 
-#define MOLDYN_CRITICAL_EST_TEMP       5.0
-
 #define MOLDYN_INTEGRATE_VERLET                0x00
 #define MOLDYN_INTEGRATE_DEFAULT       MOLDYN_INTEGRATE_VERLET
 
@@ -271,56 +281,57 @@ typedef struct s_tersoff_mult_params {
 #define TRUE                           1
 #define FALSE                          0
 
-#define ACCEPTABLE_ERROR               1e-15
-
 /*
  *
  * phsical values / constants
  *
  */
 
-#define K_BOLTZMANN            1.3807e-27                      /* Nm/K */
-#define AMU                    1.660540e-27                    /* kg */
-#define EV                     1.60217733e-19                  /* Nm */
-
-#define FCC                    0x01
-#define DIAMOND                        0x02
+#define K_BOLTZMANN            (1.380650524e-23*METER*NEWTON)  /* NA/K */
+#define EV                     (1.6021765314e-19*METER*NEWTON) /* NA */
 
 #define C                      0x06
-#define M_C                    (12.011*AMU)
+#define M_C                    12.011                          /* amu */
 
 #define SI                     0x0e
-#define LC_SI                  0.543105e-9                             /* m */
-#define M_SI                   (28.085*AMU)                            /* kg */
-#define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* m */
-#define LJ_EPSILON_SI          (2.1678*EV)                             /* Nm */
-
-#define TM_R_SI                        2.7e-10                         /* m */
-#define TM_S_SI                        3.0e-10                         /* m */
-#define TM_A_SI                        (1830.8*EV)                     /* Nm */
-#define TM_B_SI                        (471.18*EV)                     /* Nm */
-#define TM_LAMBDA_SI           2.4799e10                       /* 1/m */
-#define TM_MU_SI               1.7322e10                       /* 1/m */
+#define LC_SI                  (0.543105e-9*METER)             /* A */
+#define M_SI                   28.08553                        /* amu */
+#define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
+#define LJ_EPSILON_SI          (2.1678*EV)                             /* NA */
+
+#define TM_R_SI                        (2.7e-10*METER)                 /* A */
+#define TM_S_SI                        (3.0e-10*METER)                 /* A */
+#define TM_A_SI                        (1830.8*EV)                     /* NA */
+#define TM_B_SI                        (471.18*EV)                     /* NA */
+#define TM_LAMBDA_SI           (2.4799e10/METER)               /* 1/A */
+#define TM_MU_SI               (1.7322e10/METER)               /* 1/A */
 #define TM_BETA_SI             1.1000e-6
 #define TM_N_SI                        0.78734
 #define TM_C_SI                        1.0039e5
 #define TM_D_SI                        16.217
-#define TM_H_SI                        (-0.59825)
-
-#define TM_R_C                 1.8e-10                         /* m */
-#define TM_S_C                 2.1e-10                         /* m */
-#define TM_A_C                 (1393.6*EV)                     /* Nm */
-#define TM_B_C                 (346.7*EV)                      /* Nm */
-#define TM_LAMBDA_C            3.4879e10                       /* 1/m */
-#define TM_MU_C                        2.2119e10                       /* 1/m */
+#define TM_H_SI                        -0.59825
+
+#define TM_R_C                 (1.8e-10*METER)                 /* A */
+#define TM_S_C                 (2.1e-10*METER)                 /* A */
+#define TM_A_C                 (1393.6*EV)                     /* NA */
+#define TM_B_C                 (346.7*EV)                      /* NA */
+#define TM_LAMBDA_C            (3.4879e10/METER)               /* 1/A */
+#define TM_MU_C                        (2.2119e10/METER)               /* 1/A */
 #define TM_BETA_C              1.5724e-7
 #define TM_N_C                 0.72751
 #define TM_C_C                 3.8049e4
 #define TM_D_C                 4.384
-#define TM_H_C                 (-0.57058)
+#define TM_H_C                 -0.57058
 
 #define TM_CHI_SIC             0.9776
 
+/*
+ * lattice constants
+ */
+
+#define FCC                    0x01
+#define DIAMOND                        0x02
+
 
 /*
  *

diff --git a/sic.c b/sic.c
index 0d3b798..8c0e648 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -107,20 +107,20 @@ int main(int argc,char **argv) {
        set_pbc(&md,TRUE,TRUE,TRUE);
 
        /* create the lattice / place atoms */
-       printf("[sic] creating atoms\n");
-       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);
+       //printf("[sic] creating atoms\n");
+       //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);
 
        /* testing configuration */
-       //r.x=-0.55*0.25*sqrt(3.0)*LC_SI;       v.x=0;
-       //r.y=0;                        v.y=0;
-       //r.z=0;                        v.z=0;
-       //add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
-       //r.x=+0.55*0.25*sqrt(3.0)*LC_SI;       v.x=0;
-       //r.y=0;                        v.y=0;
-       //r.z=0;                        v.z=0;
-       //add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
+       r.x=-0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
+       r.y=0;                  v.y=0;
+       r.z=0;                  v.z=0;
+       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
+       r.x=+0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
+       r.y=0;                  v.y=0;
+       r.z=0;                  v.z=0;
+       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -140,7 +140,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,100,1.0e-15);
+       moldyn_add_schedule(&md,100,1.0);
 
        /* activate logging */
        printf("[sic] activate logging\n");

diff --git a/visual/visual.c b/visual/visual.c
index 687c561..d03903b 100644 (file)
--- a/visual/visual.c
+++ b/visual/visual.c
@@ -74,11 +74,11 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        char file[128+64];
        t_3dvec dim;
 
-       dim.x=10e9*v->dim.x;
-       dim.y=10e9*v->dim.y;
-       dim.z=10e9*v->dim.z;
+       dim.x=v->dim.x;
+       dim.y=v->dim.y;
+       dim.z=v->dim.z;
 
-       sprintf(file,"%s-%.15f.xyz",v->fb,time);
+       sprintf(file,"%s-%07f.xyz",v->fb,time);
        fd=open(file,O_WRONLY|O_CREAT|O_TRUNC);
        if(fd<0) {
                perror("open visual save file fd");
@@ -93,7 +93,7 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dprintf(v->fd,"set ambient 20\n");
        dprintf(v->fd,"set specular on\n");
        dprintf(v->fd,"set boundbox on\n");
-       dprintf(v->fd,"label\n");
+       //dprintf(v->fd,"label\n");
        sprintf(file,"%s-%.15f.ppm",v->fb,time);
        dprintf(v->fd,"write ppm %s\n",file);
        dprintf(v->fd,"zap\n");
@@ -103,9 +103,9 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dprintf(fd,"atoms at time %.15f\n",time);
        for(i=0;i<n;i++)
                dprintf(fd,"%s %f %f %f\n",pse[atom[i].element],
-                                          10e9*atom[i].r.x,
-                                          10e9*atom[i].r.y,
-                                          10e9*atom[i].r.z);
+                                          atom[i].r.x,
+                                          atom[i].r.y,
+                                          atom[i].r.z);
        if(dim.x) {
                dprintf(fd,"%s %f %f %f\n",pse[0],dim.x/2,dim.y/2,dim.z/2);
                dprintf(fd,"%s %f %f %f\n",pse[0],-dim.x/2,dim.y/2,dim.z/2);