From: hackbard Date: Thu, 5 Jul 2007 14:58:10 +0000 (+0000) Subject: added fluctuation calc code X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=25c22fe95e80365056c6a7fadc548119360ca8ce;p=physik%2Fposic.git added fluctuation calc code --- diff --git a/Makefile b/Makefile index bb520c4..ce40484 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,4 @@ -CC=gcc-3.4 -#CC=gcc +CC=gcc CFLAGS=-Wall CFLAGS+=-O3 CFLAGS+=-ffloat-store diff --git a/fluctuation_calc.c b/fluctuation_calc.c new file mode 100644 index 0000000..210f382 --- /dev/null +++ b/fluctuation_calc.c @@ -0,0 +1,121 @@ +/* + * calcultae fluctuations and related values + * + * author: frank.zirkelbach@physik.uni-augsburg.de + * + */ + +#include +#include +#include +#include +#include +#include +#include + +#define AMU 1.6605388628e-27 +#define METER 1e10 +#define SECOND 1e15 +#define KILOGRAM (1.0/AMU) +#define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) +#define JOULE (NEWTON*METER) +#define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) +#define K_B2 (K_BOLTZMANN*K_BOLTZMANN) +#define EV (1.6021765314e-19*METER*NEWTON) + +int get_line(int fd,char *line,int max) { + + int count,ret; + + count=0; + + while(1) { + if(count==max) return count; + ret=read(fd,line+count,1); + if(ret<=0) return ret; + if(line[count]=='\n') { + line[count]='\0'; + return count+1; + } + count+=1; + } +} + +int main(int argc,char **argv) { + + double from,to; + int fd; + char file[256]; + char buf[256]; + double e,e2,de2,val,time; + char *ptr; + int i,count,lnr; + double np,m,t; + + if(argc<5) { + printf("usage:\n"); + printf("%s \n",argv[0]); + printf("\n"); + printf("optional add: "); + printf("\n"); + return -1; + } + + strncpy(file,argv[1],255); + from=atof(argv[2]); + to=atof(argv[3]); + lnr=atoi(argv[4]); + np=0; + m=0.0; + t=0.0; + + if(argc==8) { + np=atoi(argv[5]); + m=atof(argv[6])*AMU*np; + t=atof(argv[7])+273.0; + } + + fd=open(file,O_RDONLY); + if(fd<2) { + perror("fd open"); + return -1; + } + + count=0; + e=0.0; + e2=0.0; + + while(1) { + if(get_line(fd,buf,255)<=0) + break; + ptr=strtok(buf," "); + time=atof(ptr); + if((time>=from)&(time<=to)) { + count+=1; + for(i=1;i fluctuation [eV/atom]: %.20f\n",de2); + if(argc==8) { + de2*=(np*np)*(EV*EV); + printf(" specific heat capacity (T=%f K) [J/(kg K)]:\n", + t); + printf(" nve: %f\n", + (1.5*np*K_BOLTZMANN/(1.0-de2/(1.5*np*K_B2*t*t))/ + (m*JOULE))); + printf(" nvt: %f\n", + (de2/(K_BOLTZMANN*t*t)+1.5*np*K_BOLTZMANN)/ + (m*JOULE)); + } + + close(fd); + + return 0; +} + diff --git a/moldyn.c b/moldyn.c index cd2a803..16c811e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -218,6 +218,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) { return 0; } +int set_mean_skip(t_moldyn *moldyn,int skip) { + + printf("[moldyn] skip %d steps before starting average calc\n",skip); + moldyn->mean_skip=skip; + + return 0; +} + int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { strncpy(moldyn->vlsdir,dir,127); @@ -728,8 +736,12 @@ double temperature_calc(t_moldyn *moldyn) { /* assume up to date kinetic energy, which is 3/2 N k_B T */ moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); + + if(moldyn->total_stepsmean_skip) + return 0; + moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/moldyn->total_steps; + moldyn->mean_t=moldyn->t_sum/(moldyn->total_steps+1-moldyn->mean_skip); return moldyn->t; } @@ -826,41 +838,54 @@ double pressure_calc(t_moldyn *moldyn) { /* virial sum and mean virial */ moldyn->virial_sum+=v; - moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps; + if(moldyn->total_steps>=moldyn->mean_skip) + moldyn->mean_v=moldyn->virial_sum/ + (moldyn->total_steps+1-moldyn->mean_skip); /* assume up to date kinetic energy */ moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; moldyn->p/=(3.0*moldyn->volume); - moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; + if(moldyn->total_steps>=moldyn->mean_skip) { + moldyn->p_sum+=moldyn->p; + moldyn->mean_p=moldyn->p_sum/ + (moldyn->total_steps+1-moldyn->mean_skip); + } /* pressure from 'absolute coordinates' virial */ virial=&(moldyn->virial); v=virial->xx+virial->yy+virial->zz; moldyn->gp=2.0*moldyn->ekin+v; moldyn->gp/=(3.0*moldyn->volume); - moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps; + if(moldyn->total_steps>=moldyn->mean_skip) { + moldyn->gp_sum+=moldyn->gp; + moldyn->mean_gp=moldyn->gp_sum/ + (moldyn->total_steps+1-moldyn->mean_skip); + } return moldyn->p; } int energy_fluctuation_calc(t_moldyn *moldyn) { + if(moldyn->total_stepsmean_skip) + return 0; + /* assume up to date energies */ /* kinetic energy */ moldyn->k_sum+=moldyn->ekin; moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin); - moldyn->k_mean=moldyn->k_sum/moldyn->total_steps; - moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps; + moldyn->k_mean=moldyn->k_sum/(moldyn->total_steps+1-moldyn->mean_skip); + moldyn->k2_mean=moldyn->k2_sum/ + (moldyn->total_steps+1-moldyn->mean_skip); moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean); /* potential energy */ moldyn->v_sum+=moldyn->energy; moldyn->v2_sum+=(moldyn->energy*moldyn->energy); - moldyn->v_mean=moldyn->v_sum/moldyn->total_steps; - moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps; + moldyn->v_mean=moldyn->v_sum/(moldyn->total_steps+1-moldyn->mean_skip); + moldyn->v2_mean=moldyn->v2_sum/ + (moldyn->total_steps+1-moldyn->mean_skip); moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean); return 0; @@ -870,6 +895,10 @@ int get_heat_capacity(t_moldyn *moldyn) { double temp2,ighc; + /* averages needed for heat capacity calc */ + if(moldyn->total_stepsmean_skip) + return 0; + /* (temperature average)^2 */ temp2=moldyn->mean_t*moldyn->mean_t; printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", @@ -890,6 +919,7 @@ int get_heat_capacity(t_moldyn *moldyn) { printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE); printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE); +printf(" --> sim: %f experimental: %f\n",moldyn->dv2_mean,1.5*moldyn->count*K_B2*moldyn->mean_t*moldyn->mean_t*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); return 0; } diff --git a/moldyn.h b/moldyn.h index 2d1a438..868db62 100644 --- a/moldyn.h +++ b/moldyn.h @@ -102,6 +102,8 @@ typedef struct s_moldyn { t_linkcell lc; /* linked cell list interface */ + int mean_skip; /* amount of steps without average calc */ + double t_ref; /* reference temperature */ double t; /* actual temperature */ double t_sum; /* sum over all t */ @@ -309,7 +311,7 @@ typedef struct s_moldyn { #define TM_CHI_SIC 0.9776 -#define TM_LC_3C_SIC (0.432e-9*METER) /* A */ +#define TM_LC_SIC 4.32 /* A */ #define ALBE_R_SI (2.82-0.14) #define ALBE_S_SI (2.82+0.14) @@ -323,7 +325,7 @@ typedef struct s_moldyn { #define ALBE_D_SI 0.81472 #define ALBE_H_SI 0.259 -#define LC_SI_ALBE 5.429 +#define ALBE_LC_SI 5.429 #define ALBE_R_C (2.00-0.15) #define ALBE_S_C (2.00+0.15) @@ -337,7 +339,7 @@ typedef struct s_moldyn { #define ALBE_D_C 6.28433 #define ALBE_H_C 0.5556 -#define LC_C_ALBE 3.566 +#define ALBE_LC_C 3.566 #define ALBE_R_SIC (2.40-0.20) #define ALBE_S_SIC (2.40+0.10) @@ -351,7 +353,7 @@ typedef struct s_moldyn { #define ALBE_D_SIC 180.314 #define ALBE_H_SIC 0.68 -#define LC_SIC_ALBE 4.359 +#define ALBE_LC_SIC 4.359 /* @@ -393,6 +395,8 @@ int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func); int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func); int set_potential_params(t_moldyn *moldyn,void *params); +int set_mean_skip(t_moldyn *moldyn,int skip); + int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); diff --git a/random/random.c b/random/random.c index df77980..902cdef 100644 --- a/random/random.c +++ b/random/random.c @@ -36,6 +36,7 @@ int rand_init(t_random *random,char *randomfile,int logfd) { return -1; } } + random->buffer=NULL; random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int)); if(random->buffer==NULL) { perror("malloc random buffer"); diff --git a/sic.c b/sic.c index 426998a..30de328 100644 --- a/sic.c +++ b/sic.c @@ -210,12 +210,13 @@ int main(int argc,char **argv) { /* set (initial) dimensions of simulation volume */ #ifdef ALBLE - set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE); - //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE); - //set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE); + //set_dim(&md,8*ALBE_LC_SI,8*ALBE_LC_SI,8*ALBE_LC_SI,TRUE); + //set_dim(&md,8*ALBE_LC_C,8*ALBE_LC_C,8*ALBE_LC_C,TRUE); + set_dim(&md,8*ALBE_LC_SIC,8*ALBE_LC_SIC,8*ALBE_LC_SIC,TRUE); #else - set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE); - //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE); + //set_dim(&md,8*LC_SI,8*LC_SI,8*LC_SI,TRUE); + //set_dim(&md,8*LC_C,8*LC_C,8*LC_C,TRUE); + set_dim(&md,8*TM_LC_SIC,8*TM_LC_SIC,8*TM_LC_SIC,TRUE); #endif /* set periodic boundary conditions in all directions */ @@ -223,28 +224,40 @@ int main(int argc,char **argv) { /* create the lattice / place atoms */ #ifdef ALBLE - create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI, + //create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI, + //create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C, #else - create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, - //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C, + //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, #endif - ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, // ATOM_ATTR_2BP|ATOM_ATTR_HB, - 0,6,6,6,NULL); - // 1,6,6,6,NULL); + // 0,8,8,8,NULL); + // 1,8,8,8,NULL); - /* create centered zinc blende lattice */ - /* - r.x=0.5*0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,LC_SIC_ALBE,SI,M_SI, + /* create zinkblende structure */ + /**/ +#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, + 0,8,8,8,&r); + r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, - 0,6,6,6,&r); - r.x+=0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,LC_SIC_ALBE,C,M_C, + 1,8,8,8,&r); +#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, + 0,8,8,8,&r); + r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,TM_LC_SIC,C,M_C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, - 1,6,6,6,&r); - */ + 1,8,8,8,&r); +#endif + /**/ + /* check for right atom placing */ moldyn_bc_check(&md); /* testing configuration */ @@ -283,6 +296,9 @@ int main(int argc,char **argv) { set_temperature(&md,atof(argv[2])+273.0); set_pressure(&md,BAR); + /* set amount of steps to skip before average calc */ + set_mean_skip(&md,500); + /* set p/t scaling */ //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001, @@ -295,7 +311,7 @@ int main(int argc,char **argv) { /* create the simulation schedule */ /* initial configuration */ - moldyn_add_schedule(&md,3000,1.0); + moldyn_add_schedule(&md,1500,1.0); //moldyn_add_schedule(&md,1000,1.0); //moldyn_add_schedule(&md,1000,1.0); //moldyn_add_schedule(&md,1000,1.0);