From 15b4727e1137600f8f46af027aefd2b5c7a56420 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 12 Jan 2007 19:35:12 +0000 Subject: [PATCH] still full of bugs ... --- Makefile | 3 +- moldyn.c | 149 +++++++++++++++++++++++++++++++++++++++++++++---------- moldyn.h | 11 +++- sic.c | 51 ++++++++++--------- 4 files changed, 162 insertions(+), 52 deletions(-) diff --git a/Makefile b/Makefile index 39457f7..c51f8bd 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,8 @@ CFLAGS=-Wall CFLAGS+=-O3 CFLAGS+=-ffloat-store CFLAGS+=-g -CFLAGS+=-DDEBUG +#CFLAGS+=-DDEBUG +#CFLAGS+=-DVDEBUG LDFLAGS=-lm OBJS=visual/visual.o random/random.o diff --git a/moldyn.c b/moldyn.c index c21f59f..d33de02 100644 --- a/moldyn.c +++ b/moldyn.c @@ -460,6 +460,28 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) { return 0; } +double temperature_calc(t_moldyn *moldyn) { + + double double_ekin; + int i; + t_atom *atom; + + atom=moldyn->atom; + + for(i=0;icount;i++) + double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v)); + + /* kinetic energy = 3/2 N k_B T */ + moldyn->t=double_ekin/(3.0*K_BOLTZMANN*moldyn->count); + + return moldyn->t; +} + +double get_temperature(t_moldyn *moldyn) { + + return moldyn->t; +} + int scale_velocity(t_moldyn *moldyn,u8 equi_init) { int i; @@ -478,10 +500,11 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { count=0; for(i=0;icount;i++) { if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) { - e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + e+=atom[i].mass*v3_absolute_square(&(atom[i].v)); count+=1; } } + e*=0.5; if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN); else return 0; /* no atoms involved in scaling! */ @@ -515,6 +538,34 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { return 0; } +double pressure_calc(t_moldyn *moldyn) { + + int i; + t_atom *atom; + double p1,p2,p=0; + + for(i=0;icount;i++) { + + + } + + p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1); + p1/=moldyn->volume; + + p2=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt2); + p2/=moldyn->volume; + + printf("compare pressures: %f %f\n",p1/ATM,p2/ATM); + + return moldyn->p; +} + +double get_pressure(t_moldyn *moldyn) { + + return moldyn->p; + +} + int scale_volume(t_moldyn *moldyn) { t_atom *atom; @@ -592,11 +643,6 @@ double get_e_kin(t_moldyn *moldyn) { return moldyn->ekin; } -double get_e_pot(t_moldyn *moldyn) { - - return moldyn->energy; -} - double update_e_kin(t_moldyn *moldyn) { return(get_e_kin(moldyn)); @@ -659,6 +705,9 @@ int link_cell_init(t_moldyn *moldyn) { lc->cells=lc->nx*lc->ny*lc->nz; lc->subcell=malloc(lc->cells*sizeof(t_list)); + if(lc->cells<27) + printf("[moldyn] FATAL: less then 27 subcells!\n"); + printf("[moldyn] initializing linked cells (%d)\n",lc->cells); for(i=0;icells;i++) @@ -885,13 +934,21 @@ int moldyn_integrate(t_moldyn *moldyn) { scale_volume(moldyn); /* check for log & visualization */ +//double ax; +//double ao; +//double av; if(e) { if(!(i%e)) +//ao=sqrt(0.1/M_SI); +//ax=((0.28-0.25)*sqrt(3)*LC_SI/2)*cos(ao*i); +//av=ao*(0.28-0.25)*sqrt(3)*LC_SI/2*sin(ao*i); + update_e_kin(moldyn); dprintf(moldyn->efd, "%f %f %f %f\n", - moldyn->time,update_e_kin(moldyn), + moldyn->time,moldyn->ekin, moldyn->energy, get_total_energy(moldyn)); +//moldyn->atom[0].r.x,ax,av*av*M_SI,0.1*ax*ax,av*av*M_SI+0.1*ax*ax); } if(m) { if(!(i%m)) { @@ -946,7 +1003,7 @@ int moldyn_integrate(t_moldyn *moldyn) { int velocity_verlet(t_moldyn *moldyn) { int i,count; - double tau,tau_square; + double tau,tau_square,h; t_3dvec delta; t_atom *atom; @@ -957,14 +1014,15 @@ int velocity_verlet(t_moldyn *moldyn) { for(i=0;ienergy=0.0; + + moldyn->vt2=0.0; /* get energy and force of every atom */ for(i=0;ixy=0.0; virial->xz=0.0; virial->yz=0.0; + moldyn->vt1=0.0; /* reset site energy */ itom[i].e=0.0; @@ -1134,6 +1195,30 @@ printf("\n\n"); printf("\n\n"); #endif + moldyn->vt2=0.0; + for(i=0;ivt2-=v3_scalar_product(&(itom[i].r),&(itom[i].f)); + +printf("compare: vt1: %f vt2: %f\n",moldyn->vt1,moldyn->vt2); + +pressure_calc(moldyn); + + return 0; +} + +/* + * virial calculation + */ + +inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { + + a->virial.xx-=f->x*d->x; + a->virial.yy-=f->y*d->y; + a->virial.zz-=f->z*d->z; + a->virial.xy-=f->x*d->y; + a->virial.xz-=f->x*d->z; + a->virial.yz-=f->y*d->z; + return 0; } @@ -1179,23 +1264,29 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_ho_params *params; t_3dvec force,distance; - double d; + double d,f; double sc,equi_dist; params=moldyn->pot2b_params; sc=params->spring_constant; equi_dist=params->equilibrium_distance; + if(air),&(ai->r)); if(bc) check_per_bound(moldyn,&distance); d=v3_norm(&distance); if(d<=moldyn->cutoff) { - /* energy is 1/2 (d-d0)^2, but we will add this twice ... */ - moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist)); + moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); /* f = -grad E; grad r_ij = -1 1/r_ij distance */ - v3_scale(&force,&distance,sc*(1.0-(equi_dist/d))); + f=sc*(1.0-equi_dist/d); + v3_scale(&force,&distance,f); v3_add(&(ai->f),&(ai->f),&force); + virial_calc(ai,&force,&distance); + virial_calc(aj,&force,&distance); /* f and d signe switched */ + v3_scale(&force,&distance,-f); + v3_add(&(aj->f),&(aj->f),&force); } return 0; @@ -1215,6 +1306,8 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { sig6=params->sigma6; sig12=params->sigma12; + if(air),&(ai->r)); if(bc) check_per_bound(moldyn,&distance); d=v3_absolute_square(&distance); /* 1/r^2 */ @@ -1223,16 +1316,20 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { h2=d*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ - /* energy is eps*..., but we will add this twice ... */ - moldyn->energy+=0.5*eps*(sig12*h1-sig6*h2); + moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc); h2*=d; /* 1/r^8 */ h1*=d; /* 1/r^14 */ h2*=6*sig6; h1*=12*sig12; d=+h1-h2; d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(aj->f),&(aj->f),&force); v3_scale(&force,&distance,-1.0*d); /* f = - grad E */ v3_add(&(ai->f),&(ai->f),&force); + virial_calc(ai,&force,&distance); + virial_calc(aj,&force,&distance); /* f and d signe switched */ + moldyn->vt1-=v3_scalar_product(&force,&distance); } return 0; @@ -1526,7 +1623,6 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { zeta=exchange->zeta_ij; if(zeta==0.0) { moldyn->debug++; /* just for debugging ... */ - db=0.0; b=chi; v3_scale(&force,dist_ij,df_a*b*f_c); } @@ -1602,12 +1698,13 @@ if(ai==&(moldyn->atom[0])) { v3_add(&(ai->f),&(ai->f),&force); /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ - ai->virial.xx+=force.x*dist_ij->x; - ai->virial.yy+=force.y*dist_ij->y; - ai->virial.zz+=force.z*dist_ij->z; - ai->virial.xy+=force.x*dist_ij->y; - ai->virial.xz+=force.x*dist_ij->z; - ai->virial.yz+=force.y*dist_ij->z; +// TEST ... with a minus instead + ai->virial.xx-=force.x*dist_ij->x; + ai->virial.yy-=force.y*dist_ij->y; + ai->virial.zz-=force.z*dist_ij->z; + ai->virial.xy-=force.x*dist_ij->y; + ai->virial.xz-=force.x*dist_ij->z; + ai->virial.yz-=force.y*dist_ij->z; #ifdef DEBUG if(ai==&(moldyn->atom[0])) { diff --git a/moldyn.h b/moldyn.h index 23b5100..c5cc064 100644 --- a/moldyn.h +++ b/moldyn.h @@ -82,6 +82,7 @@ typedef struct s_moldyn { t_3dvec dim; /* dimensions of the simulation volume */ double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ + double vt1,vt2; /* potential force function and parameter pointers */ int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); @@ -183,6 +184,7 @@ typedef struct s_lj_params { double sigma6; double sigma12; double epsilon4; + double uc; } t_lj_params; /* @@ -401,11 +403,15 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, int destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); +double temperature_calc(t_moldyn *moldyn); +double get_temperature(t_moldyn *moldyn); int scale_velocity(t_moldyn *moldyn,u8 equi_init); +double pressure_calc(t_moldyn *moldyn); +double get_pressure(t_moldyn *moldyn); int scale_volume(t_moldyn *moldyn); double get_e_kin(t_moldyn *moldyn); -double get_e_pot(t_moldyn *moldyn); +double update_e_kin(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); t_3dvec get_total_p(t_moldyn *moldyn); @@ -423,9 +429,10 @@ int moldyn_integrate(t_moldyn *moldyn); int velocity_verlet(t_moldyn *moldyn); int potential_force_calc(t_moldyn *moldyn); +inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) + __attribute__((always_inline)); inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) __attribute__((always_inline)); -int check_per_bound(t_moldyn *moldyn,t_3dvec *a); int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_complete_params(t_tersoff_mult_params *p); diff --git a/sic.c b/sic.c index aabb79c..e2e8626 100644 --- a/sic.c +++ b/sic.c @@ -46,11 +46,17 @@ int main(int argc,char **argv) { /* choose potential */ printf("[sic] selecting potential\n"); - set_potential1b(&md,tersoff_mult_1bp,&tp); - set_potential2b(&md,tersoff_mult_2bp,&tp); - set_potential2b_post(&md,tersoff_mult_post_2bp,&tp); - set_potential3b(&md,tersoff_mult_3bp,&tp); - //set_potential2b(&md,lennard_jones,&lj); + //set_potential1b(&md,tersoff_mult_1bp,&tp); + //set_potential2b(&md,tersoff_mult_2bp,&tp); + //set_potential2b_post(&md,tersoff_mult_post_2bp,&tp); + //set_potential3b(&md,tersoff_mult_3bp,&tp); + set_potential2b(&md,lennard_jones,&lj); + //set_potential2b(&md,harmonic_oscillator,&ho); + + /* cutoff radius */ + printf("[sic] setting cutoff radius\n"); + //set_cutoff(&md,TM_S_SI); + set_cutoff(&md,3*LC_SI); /* * potential parameters @@ -61,10 +67,11 @@ int main(int argc,char **argv) { lj.sigma6*=lj.sigma6; lj.sigma12=lj.sigma6*lj.sigma6; lj.epsilon4=4.0*LJ_EPSILON_SI; + lj.uc=lj.epsilon4*(lj.sigma12/pow(md.cutoff,12.0)-lj.sigma6/pow(md.cutoff,6)); /* harmonic oscillator */ ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; - ho.spring_constant=1; + ho.spring_constant=.1; /* * tersoff mult potential parameters for SiC @@ -97,14 +104,9 @@ int main(int argc,char **argv) { tersoff_mult_complete_params(&tp); - /* cutoff radius */ - printf("[sic] setting cutoff radius\n"); - set_cutoff(&md,TM_S_SI); - //set_cutoff(&md,2*LC_SI); - /* set (initial) dimensions of simulation volume */ printf("[sic] setting dimensions\n"); - set_dim(&md,5*LC_SI,5*LC_SI,5*LC_SI,TRUE); + set_dim(&md,10*LC_SI,10*LC_SI,10*LC_SI,TRUE); /* set periodic boundary conditions in all directions */ printf("[sic] setting periodic boundary conditions\n"); @@ -113,24 +115,26 @@ int main(int argc,char **argv) { /* 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); + // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + ATOM_ATTR_2BP|ATOM_ATTR_HB, + 0,10,10,10); moldyn_bc_check(&md); /* testing configuration */ - //r.x=2.8/2; v.x=0; + //r.x=0.28*sqrt(3)*LC_SI/2; v.x=0; + //r.x=1.75*LC_SI; v.x=-0.01; //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, - // ATOM_ATTR_2BP, + // ATOM_ATTR_2BP|ATOM_ATTR_HB, // &r,&v); - //r.x=-2.8/2; v.x=0; + //r.x=-r.x; v.x=-v.x; //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, - // ATOM_ATTR_2BP, + // ATOM_ATTR_2BP|ATOM_ATTR_HB, // &r,&v); /* setting a nearest neighbour distance for the moldyn checks */ @@ -153,21 +157,22 @@ int main(int argc,char **argv) { printf("[sic] set p/t scaling\n"); //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, // T_SCALE_BERENDSEN,100.0); - set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0); /* initial thermal fluctuations of particles (in equilibrium) */ printf("[sic] thermal init\n"); - //thermal_init(&md,TRUE); + thermal_init(&md,TRUE); /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,100,1.0); + moldyn_add_schedule(&md,10001,1.0); /* activate logging */ printf("[sic] activate logging\n"); moldyn_set_log_dir(&md,argv[1]); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,1); - moldyn_set_log(&md,VISUAL_STEP,1); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,10); + moldyn_set_log(&md,VISUAL_STEP,100); /* * let's do the actual md algorithm now -- 2.20.1