-CC=gcc-3.4
-#CC=gcc
+CC=gcc
CFLAGS=-Wall
CFLAGS+=-O3
CFLAGS+=-ffloat-store
--- /dev/null
+/*
+ * calcultae fluctuations and related values
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#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 <file> <from> <to> <linenumber>\n",argv[0]);
+ printf("\n");
+ printf("optional add: <nr particles> <mass> <temperature>");
+ 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<lnr;i++)
+ ptr=strtok(NULL," ");
+ val=atof(ptr);
+ e+=val;
+ e2+=(val*val);
+ }
+ }
+
+ de2=e2/count-e*e/(count*count);
+ printf("--> 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;
+}
+
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);
/* 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_steps<moldyn->mean_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;
}
/* 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_steps<moldyn->mean_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;
double temp2,ighc;
+ /* averages needed for heat capacity calc */
+ if(moldyn->total_steps<moldyn->mean_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",
printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+printf(" --> <dV2> 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;
}
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 */
#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)
#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)
#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)
#define ALBE_D_SIC 180.314
#define ALBE_H_SIC 0.68
-#define LC_SIC_ALBE 4.359
+#define ALBE_LC_SIC 4.359
/*
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);
return -1;
}
}
+ random->buffer=NULL;
random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int));
if(random->buffer==NULL) {
perror("malloc random buffer");
/* 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 */
/* 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 */
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,
/* 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);
projects / physik / posic.git / commitdiff