From d9e7f195bbb219ad4c2de0e5f54d023ef9e669fb Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 11 Dec 2006 22:50:51 +0000 Subject: [PATCH] scaling to fs, angstrom, amu done, probs with velocity scaling! --- Makefile | 1 + moldyn.c | 23 +++++++------- moldyn.h | 79 ++++++++++++++++++++++++++++--------------------- sic.c | 26 ++++++++-------- visual/visual.c | 16 +++++----- 5 files changed, 80 insertions(+), 65 deletions(-) diff --git a/Makefile b/Makefile index 742a7a5..6f3e39e 100644 --- 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 --- 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;icount;i++) + for(i=0;icount;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;imu_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 --- 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 --- 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 --- 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