X-Git-Url: https://www.hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=sic.c;h=375c6d54b38756aef74f12431075dd8fd6b446fc;hp=7c40a7e6c8f603f8ff375b08b40c30e8f37ce9ea;hb=HEAD;hpb=180ff86ae35ca0ea6ee857e92e56120e65315d17 diff --git a/sic.c b/sic.c index 7c40a7e..375c6d5 100644 --- a/sic.c +++ b/sic.c @@ -8,112 +8,460 @@ #include #include "moldyn.h" -#include "math/math.h" -#include "init/init.h" -#include "visual/visual.h" -#include "posic.h" +/* potential */ +#include "potentials/harmonic_oscillator.h" +#include "potentials/lennard_jones.h" +#include "potentials/albe.h" +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else +#include "potentials/tersoff.h" +#endif + +typedef struct s_hp { + int prerun_count; /* prerun count */ + int insert_count; /* insert count */ + int postrun_count; /* post run count */ + unsigned char state; /* current state */ + int argc; /* arg count */ + char **argv; /* args */ +} t_hp; + +#define STATE_PRERUN 0x00 +#define STATE_INSERT 0x01 +#define STATE_POSTRUN 0x02 + +/* include the config file */ +#include "config.h" + +int insert_atoms(t_moldyn *moldyn) { + + int i,j; + u8 run; + t_3dvec r,v,dist; + double d,dmin; + + t_atom *atom; + + atom=moldyn->atom; + + v.x=0; v.y=0; v.z=0; + + for(j=0;jatom[4372].r.x=(-0.5+0.125+0.125)*ALBE_LC_SI; + moldyn->atom[4372].r.y=(-0.5+0.125+0.125)*ALBE_LC_SI; +#endif +#ifdef INS_001DB + // 001 dumbbell + r.x=(-0.5+0.25)*ALBE_LC_SI; + r.y=(-0.5+0.25)*ALBE_LC_SI; + r.z=(-0.1)*ALBE_LC_SI; + moldyn->atom[4372].r.z=(-0.4)*ALBE_LC_SI; +#endif +#ifdef INS_USER + // 001 dumbbell + r.x=INS_UX*ALBE_LC_SI; + r.y=INS_UY*ALBE_LC_SI; + r.z=INS_UZ*ALBE_LC_SI; +#endif +#ifdef INS_RAND + // random +#ifdef INS_DYNAMIC_LEN + r.x=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.x; + r.y=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.y; + r.z=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.z; +#else + r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX; + r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY; + r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ; +#endif +#endif + // offset + r.x+=INS_OFFSET; + r.y+=INS_OFFSET; + r.z+=INS_OFFSET; + /* assume valid coordinates */ + run=0; + dmin=10000000000.0; // for sure too high! + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); + v3_sub(&dist,&(atom->r),&r); + check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + /* reject coordinates */ + if(dcount-1,r.x,r.y,r.z,dmin); + } + + return 0; +} + +int sic_hook(void *moldyn,void *hook_params) { + + t_hp *hp; + t_moldyn *md; + int steps; + double tau; + double dt; + double dp; + + hp=hook_params; + md=moldyn; + + tau=1.0; + steps=0; + + /* switch on t scaling */ + if(md->schedule.count==0) + set_pt_scale(md,P_SCALE_BERENDSEN,P_SCALE_TAU, + T_SCALE_BERENDSEN,T_SCALE_TAU); + + /* my lousy state machine ! */ + + /* switch to insert state immediately */ + if(hp->state==STATE_PRERUN) + hp->state=STATE_INSERT; + + switch(hp->state) { + case STATE_INSERT: + goto insert; + break; + case STATE_POSTRUN: + goto postrun; + break; + default: + printf("[sic hook] unknown state\n"); + return -1; + } + + /* act according to state */ + +insert: + + /* assigne values */ + steps=INS_RELAX; + tau=INS_TAU; + + /* check temperature */ + dt=md->t_avg-md->t_ref; + dp=md->p_avg-md->p_ref; + if(dt<0) + dt=-dt; + if(dp<0) + dp=-dp; + if((dt>INS_DELTA_TC)|(dp>INS_DELTA_PC)) + goto addsched; + + /* immediately go on if no job is to be done */ + if(hp->insert_count==INS_RUNS) { + printf(" --- leaving insert state ---\n"); + hp->state=STATE_POSTRUN; + goto postrun; + } + + /* else -> insert atoms */ + hp->insert_count+=1; + printf(" ### insert atoms (%d/%d) ###\n", + hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS); + insert_atoms(md); + goto addsched; + +postrun: + + /* assigne values */ + steps=POST_RELAX; + tau=POST_TAU; + + /* check temperature */ + dt=md->t_avg-md->t_ref; + dp=md->p_avg-md->p_ref; + if(dt<0) + dt=-dt; + if(dp<0) + dp=-dp; + if((dt>POST_DELTA_TC)|(dp>POST_DELTA_PC)) + goto addsched; + + /* immediately return if no job is to be done */ + if(hp->postrun_count==POST_RUNS) { + printf(" --- leaving post run state ---\n"); + return 0; + } + + /* postrun action */ + hp->postrun_count+=1; + printf(" ### postrun (%d/%d) ###\n", + hp->postrun_count,POST_RUNS); + set_temperature(md,md->t_ref-POST_DT); + +addsched: + + /* reset the average counters */ + average_reset(md); + + /* add schedule */ + moldyn_add_schedule(md,steps,tau); + + return 0; +} int main(int argc,char **argv) { /* main moldyn structure */ t_moldyn md; - /* potential parameters */ - t_lj_params lj; - t_ho_params ho; - t_tersoff_mult_params tp; + /* hook parameter structure */ + t_hp hookparam; - /* misc variables, mainly to initialize stuff */ + /* testing location & velocity vector */ t_3dvec r,v; - - /* temperature */ - double t; + memset(&r,0,sizeof(t_3dvec)); + memset(&v,0,sizeof(t_3dvec)); /* initialize moldyn */ - printf("[sic] moldyn init\n"); moldyn_init(&md,argc,argv); /* choose integration algorithm */ - printf("[sic] setting integration algorithm\n"); set_int_alg(&md,MOLDYN_INTEGRATE_VERLET); /* choose potential */ - printf("[sic] selecting potential\n"); - //set_potential1b(&md,tersoff_mult_1bp,&tp); - //set_potential2b(&md,tersoff_mult_2bp,&tp); - //set_potential3b(&md,tersoff_mult_3bp,&tp); - set_potential2b(&md,lennard_jones,&lj); +#ifdef ALBE + if(set_potential(&md,MOLDYN_POTENTIAL_AM)<0) + return -1; +#else + if(set_potential(&md,MOLDYN_POTENTIAL_TM)<0) + return -1; +#endif + + /* cutoff radius & bondlen */ +#ifdef ALBE + set_cutoff(&md,ALBE_S_SI); + //set_cutoff(&md,ALBE_S_C); +#else + set_cutoff(&md,TM_S_SI); + //set_cutoff(&md,TM_S_C); +#endif /* * potential parameters */ - /* lennard jones */ - lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI; - lj.sigma6*=lj.sigma6; - lj.sigma12=lj.sigma6*lj.sigma6; - lj.epsilon4=4.0*LJ_EPSILON_SI; - - /* harmonic oscillator */ - ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; - ho.spring_constant=1; - - /* cutoff radius */ - printf("[sic] setting cutoff radius\n"); - set_cutoff(&md,1*LC_SI); +#ifndef ALBE + /* + * tersoff mult potential parameters for SiC + */ + tersoff_mult_set_params(&md,SI,C); +#else + /* + * albe mult potential parameters for SiC + */ + albe_mult_set_params(&md,SI,C); +#endif /* set (initial) dimensions of simulation volume */ - printf("[sic] setting dimensions\n"); - set_dim(&md,4*LC_SI,4*LC_SI,4*LC_SI,TRUE); +#ifdef ALBE + #ifdef INIT_SI + set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE); + #endif + #ifdef INIT_C + set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE); + #endif + #ifdef INIT_3CSIC + set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE); + #endif +#else + #ifdef INIT_SI + set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE); + #endif + #ifdef INIT_C + set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE); + #endif + #ifdef INIT_3CSIC + set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE); + #endif +#endif /* set periodic boundary conditions in all directions */ - printf("[sic] setting periodic boundary conditions\n"); set_pbc(&md,TRUE,TRUE,TRUE); /* create the lattice / place atoms */ - printf("[sic] creating atoms\n"); - //memset(&v,0,sizeof(t_3dvec)); - //r.y=0; - //r.z=0; - //r.x=0.23*sqrt(3.0)*LC_SI/2.0; - //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v); - //r.x=-r.x; - //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v); - create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,ATOM_ATTR_2BP,0,4,4,4); - - /* setting a nearest neighbour distance for the moldyn checks */ - set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */ - - /* set temperature */ - printf("[sic] setting temperature\n"); - set_temperature(&md,273.0); + + // diamond +#ifdef ALBE + #ifdef INIT_SI + create_lattice(&md,DIAMOND,ALBE_LC_SI,SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); + #endif + #ifdef INIT_C + create_lattice(&md,DIAMOND,ALBE_LC_C,C, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); + #endif +#else + #ifdef INIT_SI + create_lattice(&md,DIAMOND,LC_SI,SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); + #endif + #ifdef INIT_C + create_lattice(&md,DIAMOND,LC_C,SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); + #endif +#endif + + // zinkblende +#ifdef INIT_3CSIC + #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, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 0,LCNTX,LCNTY,LCNTZ,&r,0,NULL); + r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,ALBE_LC_SIC,C, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB, + 1,LCNTX,LCNTY,LCNTZ,&r,0,NULL); + #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, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 0,LCNTX,LCNTY,LCNTZ,&r,0,NULL); + r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x; + create_lattice(&md,FCC,TM_LC_SIC,C, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + 1,LCNTX,LCNTY,LCNTZ,&r,0,NULL); + #endif +#endif + + /* check for right atom placing */ + moldyn_bc_check(&md); + + /* testing configuration */ + //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0; + //r.x=(TM_S_SI+TM_R_SI)/4.0; 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, + // ATOM_ATTR_2BP|ATOM_ATTR_HB, + // &r,&v); + //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_HB, + // &r,&v); + //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0; + //r.x=(TM_S_SI+TM_R_SI)/4.0; v.x=0; + //r.y=0; v.y=0; + //r.x=0; v.x=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_HB, + // &r,&v); + //r.z=-r.z; v.z=-v.z; + //r.y=0; v.y=0; + //r.x=0; v.x=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_HB, + // &r,&v); + + /* set temperature & pressure */ + set_temperature(&md,atof(argv[2])+273.0); + set_pressure(&md,0.0); + + /* set amount of steps to skip before average calc */ + set_avg_skip(&md,AVG_SKIP); + + /* set p/t scaling */ + //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,P_SCALE_BERENDSEN,0.01/(100*GPA), + // T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); + //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0); - /* initial thermal fluctuations of particles */ - printf("[sic] thermal init\n"); - thermal_init(&md); + /* initial thermal fluctuations of particles (in equilibrium) */ + thermal_init(&md,TRUE); /* create the simulation schedule */ - printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,10000,1.0e-15); + moldyn_add_schedule(&md,PRERUN,PRE_TAU); + + /* schedule hook function */ + memset(&hookparam,0,sizeof(t_hp)); + hookparam.argc=argc; + hookparam.argv=argv; + moldyn_set_schedule_hook(&md,&sic_hook,&hookparam); + //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam); + //moldyn_add_schedule(&md,POSTRUN,1.0); /* activate logging */ - printf("[sic] activate logging\n"); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",100); - moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",100); + moldyn_set_log_dir(&md,argv[1]); + moldyn_set_report(&md,"Frank Zirkelbach",R_TITLE); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E); + moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T); + moldyn_set_log(&md,LOG_PRESSURE,LOG_P); + moldyn_set_log(&md,LOG_VOLUME,LOG_V); + moldyn_set_log(&md,VISUAL_STEP,LOG_A); + moldyn_set_log(&md,SAVE_STEP,LOG_S); + moldyn_set_log(&md,CREATE_REPORT,0); + + /* next neighbour distance for critical checking */ + set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0)); /* * let's do the actual md algorithm now * * integration of newtons equations */ - - printf("[sic] integration start, go get a coffee ...\n"); moldyn_integrate(&md); +#ifdef dEBUG +return 0; +#endif - /* close */ + /* + * post processing the data + */ - printf("[sic] shutdown\n"); + /* close */ moldyn_shutdown(&md); return 0;