From a9fbc66448c52bc4138176739b33d17ba86b7eae Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 1 Dec 2006 14:59:08 +0000 Subject: [PATCH] ready for tersoff 3bp debugging --- moldyn.c | 170 ++++++++++++++++++++++++++++++++++--------------------- moldyn.h | 24 ++++++-- sic.c | 16 +++++- 3 files changed, 136 insertions(+), 74 deletions(-) diff --git a/moldyn.c b/moldyn.c index 80f4fa1..af0132d 100644 --- a/moldyn.c +++ b/moldyn.c @@ -71,9 +71,18 @@ int set_cutoff(t_moldyn *moldyn,double cutoff) { return 0; } -int set_temperature(t_moldyn *moldyn,double t) { +int set_temperature(t_moldyn *moldyn,double t_ref) { - moldyn->t=t; + moldyn->t_ref=t_ref; + + return 0; +} + +int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { + + moldyn->pt_scale=(ptype|ttype); + moldyn->t_tc=ttc; + moldyn->p_tc=ptc; return 0; } @@ -278,7 +287,7 @@ int destroy_atoms(t_moldyn *moldyn) { return 0; } -int thermal_init(t_moldyn *moldyn) { +int thermal_init(t_moldyn *moldyn,u8 equi_init) { /* * - gaussian distribution of velocities @@ -298,7 +307,7 @@ int thermal_init(t_moldyn *moldyn) { /* gaussian distribution of velocities */ v3_zero(&p_total); for(i=0;icount;i++) { - sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass); + sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass); /* x direction */ v=sigma*rand_get_gauss(random); atom[i].v.x=v; @@ -321,12 +330,12 @@ int thermal_init(t_moldyn *moldyn) { } /* velocity scaling */ - scale_velocity(moldyn,VSCALE_INIT_EQUI); + scale_velocity(moldyn,equi_init); return 0; } -int scale_velocity(t_moldyn *moldyn,u8 type) { +int scale_velocity(t_moldyn *moldyn,u8 equi_init) { int i; double e,scale; @@ -339,29 +348,40 @@ int scale_velocity(t_moldyn *moldyn,u8 type) { * - velocity scaling (E = 3/2 N k T), E: kinetic energy */ + /* get kinetic energy / temperature & count involved atoms */ e=0.0; count=0; for(i=0;icount;i++) { - if(atom[i].attr&ATOM_ATTR_HB) { + if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) { e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); count+=1; } } - - /* temporary hack for e,t = 0 */ + if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN); + else return 0; /* no atoms involved in scaling! */ + + /* (temporary) hack for e,t = 0 */ if(e==0.0) { - if(moldyn->t!=0.0) - thermal_init(moldyn); + moldyn->t=0.0; + if(moldyn->t_ref!=0.0) + thermal_init(moldyn,equi_init); else - return 0; + return 0; /* no scaling needed */ } - /* direct scaling */ - scale=(1.5*count*K_BOLTZMANN*moldyn->t)/e; - if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */ + + /* get scaling factor */ + scale=moldyn->t_ref/moldyn->t; + if(equi_init&TRUE) + scale*=2.0; + else + if(moldyn->pt_scale&T_SCALE_BERENDSEN) + scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc; scale=sqrt(scale); + + /* velocity scaling */ for(i=0;icount;i++) - if(atom[i].attr&ATOM_ATTR_HB) + if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) v3_scale(&(atom[i].v),&(atom[i].v),scale); return 0; @@ -666,6 +686,10 @@ int moldyn_integrate(t_moldyn *moldyn) { /* integration step */ moldyn->integrate(moldyn); + /* p/t scaling */ + if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) + scale_velocity(moldyn,FALSE); + /* increase absolute time */ moldyn->time+=moldyn->tau; @@ -777,7 +801,9 @@ int potential_force_calc(t_moldyn *moldyn) { int i,j,k,count; t_atom *itom,*jtom,*ktom; t_linkcell *lc; - t_list neighbour_i[27],neighbour_j[27]; + t_list neighbour_i[27]; + t_list neighbour_i2[27]; + //t_list neighbour_j[27]; t_list *this,*that; u8 bc_ij,bc_ijk; int countn,dnlc; @@ -839,47 +865,61 @@ int potential_force_calc(t_moldyn *moldyn) { !(jtom->attr&ATOM_ATTR_3BP)) continue; - link_cell_neighbour_index(moldyn, - (jtom->r.x+moldyn->dim.x/2)/lc->x, - (jtom->r.y+moldyn->dim.y/2)/lc->y, - (jtom->r.z+moldyn->dim.z/2)/lc->z, - neighbour_j); - - /* neighbours of j */ - for(k=0;kcountn;k++) { - - that=&(neighbour_j[k]); - list_reset(that); - - if(that->start==NULL) - continue; - - bc_ijk=(kdnlc)?0:1; - - do { - - ktom=that->current->data; - - if(!(ktom->attr&ATOM_ATTR_3BP)) - continue; - - if(ktom==jtom) - continue; - - if(ktom==&(itom[i])) - continue; - - moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); - - } while(list_next(that)!=\ - L_NO_NEXT_ELEMENT); - - } - - /* neighbours of i */ + /* + * according to mr. nordlund, we dont need to take the + * sum over all atoms now, as 'this is centered' around + * atom i ... + * i am not quite sure though! there is a not vanishing + * part even if f_c_ik is zero ... + * this analytical potentials suck! + * switching from mc to md to dft soon! + */ + + // link_cell_neighbour_index(moldyn, + // (jtom->r.x+moldyn->dim.x/2)/lc->x, + // (jtom->r.y+moldyn->dim.y/2)/lc->y, + // (jtom->r.z+moldyn->dim.z/2)/lc->z, + // neighbour_j); + +// /* neighbours of j */ +// for(k=0;kcountn;k++) { +// +// that=&(neighbour_j[k]); +// list_reset(that); +// +// if(that->start==NULL) +// continue; +// +// bc_ijk=(kdnlc)?0:1; +// +// do { +// +// ktom=that->current->data; +// +// if(!(ktom->attr&ATOM_ATTR_3BP)) +// continue; +// +// if(ktom==jtom) +// continue; +// +// if(ktom==&(itom[i])) +// continue; +// +// moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); +// +/* } while(list_next(that)!=\ */ +// L_NO_NEXT_ELEMENT); +// +// } + + /* copy the neighbour lists */ + memcpy(neighbour_i2,neighbour_i, + 27*sizeof(t_list)); + + /* get neighbours of i */ for(k=0;kstart==NULL) @@ -900,7 +940,9 @@ int potential_force_calc(t_moldyn *moldyn) { if(ktom==&(itom[i])) continue; +printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); +printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); } while(list_next(that)!=\ L_NO_NEXT_ELEMENT); @@ -1094,8 +1136,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { if(bc) check_per_bound(moldyn,&dist_ij); - /* save for use in 3bp */ /* REALLY ?!?!?! */ - exchange->dist_ij=dist_ij; + d_ij=v3_norm(&dist_ij); + + /* save for use in 3bp */ + exchange->dist_ij=dist_ij; /* <- needed ? */ + exchange->d_ij=d_ij; /* constants */ if(num==aj->bnum) { @@ -1121,11 +1166,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { params->exchange.chi=params->chi; } - d_ij=v3_norm(&dist_ij); - - /* save for use in 3bp */ - exchange->d_ij=d_ij; - if(d_ij>S) return 0; @@ -1153,8 +1193,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* add forces */ v3_add(&(ai->f),&(ai->f),&force); - /* energy is 0.5 f_r f_c, but we will sum it up twice ... */ - moldyn->energy+=(0.25*f_r*f_c); + /* energy is 0.5 f_r f_c ... */ + moldyn->energy+=(0.5*f_r*f_c); /* save for use in 3bp */ exchange->f_c=f_c; diff --git a/moldyn.h b/moldyn.h index 2d4b563..8b40d4c 100644 --- a/moldyn.h +++ b/moldyn.h @@ -85,7 +85,16 @@ typedef struct s_moldyn { t_linkcell lc; /* linked cell list interface */ - double t; /* temperature */ + double t_ref; /* reference temperature */ + double t; /* actual temperature */ + + double p_ref; /* reference pressure */ + double p; /* actual pressure */ + + /* pressure and temperature control (velocity/volume scaling) */ + unsigned char pt_scale; /* type of p and t scaling */ + double t_tc; /* t berendsen control time constant */ + double p_tc; /* p berendsen control time constant */ /* simulation schedule */ t_moldyn_schedule schedule; @@ -127,8 +136,10 @@ typedef struct s_moldyn { #define MOLDYN_2BP 0x01 /* 2 body */ #define MOLDYN_3BP 0x02 /* and 3 body particle pots */ -#define VSCALE_INIT_EQUI 0x00 /* initial, eq positions */ -#define VSCALE_BERENDSEN 0x01 /* berendsen control */ +#define T_SCALE_BERENDSEN 0x01 /* berendsen t control */ +#define T_SCALE_DIRECT 0x02 /* direct t control */ +#define P_SCALE_BERENDSEN 0x04 /* berendsen p control */ +#define P_SCALE_DIRECT 0x08 /* direct p control */ /* @@ -310,7 +321,8 @@ int moldyn_shutdown(t_moldyn *moldyn); int set_int_alg(t_moldyn *moldyn,u8 algo); int set_cutoff(t_moldyn *moldyn,double cutoff); -int set_temperature(t_moldyn *moldyn,double t); +int set_temperature(t_moldyn *moldyn,double t_ref); +int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_nn_dist(t_moldyn *moldyn,double dist); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); @@ -327,8 +339,8 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr, t_3dvec *r,t_3dvec *v); int destroy_atoms(t_moldyn *moldyn); -int thermal_init(t_moldyn *moldyn); -int scale_velocity(t_moldyn *moldyn,u8 type); +int thermal_init(t_moldyn *moldyn,u8 equi_init); +int scale_velocity(t_moldyn *moldyn,u8 equi_init); double get_e_kin(t_moldyn *moldyn); double get_e_pot(t_moldyn *moldyn); diff --git a/sic.c b/sic.c index fa1cb66..b10c681 100644 --- a/sic.c +++ b/sic.c @@ -23,6 +23,12 @@ int main(int argc,char **argv) { t_ho_params ho; t_tersoff_mult_params tp; + /* misc parameters */ + double tau; + + /* values */ + tau=1.0e-15; /* delta t = 1 fs */ + /* initialize moldyn */ printf("[sic] moldyn init\n"); moldyn_init(&md,argc,argv); @@ -106,11 +112,15 @@ int main(int argc,char **argv) { /* set temperature */ printf("[sic] setting temperature\n"); - set_temperature(&md,273.0); + set_temperature(&md,0.0); + + /* set p/t scaling */ + printf("[sic] set p/t scaling\n"); + set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100*tau); - /* initial thermal fluctuations of particles */ + /* initial thermal fluctuations of particles (in equilibrium) */ printf("[sic] thermal init\n"); - thermal_init(&md); + thermal_init(&md,TRUE); /* create the simulation schedule */ printf("[sic] adding schedule\n"); -- 2.20.1