From: hackbard Date: Wed, 13 Dec 2006 17:47:42 +0000 (+0000) Subject: basic p control added, virial still needed! X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=88ea16029cc07a837977f737f4536fa3e881dcee;p=physik%2Fposic.git basic p control added, virial still needed! --- diff --git a/moldyn.c b/moldyn.c index 73becde..b6fa8d3 100644 --- a/moldyn.c +++ b/moldyn.c @@ -78,6 +78,13 @@ int set_temperature(t_moldyn *moldyn,double t_ref) { return 0; } +int set_pressure(t_moldyn *moldyn,double p_ref) { + + moldyn->p_ref=p_ref; + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); @@ -93,16 +100,19 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->dim.y=y; moldyn->dim.z=z; + moldyn->volume=x*y*z; + if(visualize) { moldyn->vis.dim.x=x; moldyn->vis.dim.y=y; moldyn->vis.dim.z=z; } - printf("[moldyn] dimensions in A:\n"); + printf("[moldyn] dimensions in A and A^2 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); printf(" y: %f\n",moldyn->dim.y); printf(" z: %f\n",moldyn->dim.z); + printf(" volume: %f\n",moldyn->volume); printf(" visualize simulation box: %s\n",visualize?"on":"off"); return 0; @@ -390,7 +400,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { count+=1; } } - if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN); + if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN); else return 0; /* no atoms involved in scaling! */ /* (temporary) hack for e,t = 0 */ @@ -423,6 +433,57 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { return 0; } +int scale_volume(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim,*vdim; + double virial,scale; + t_linkcell *lc; + int i; + + atom=moldyn->atom; + dim=&(moldyn->dim); + vdim=&(moldyn->vis.dim); + lc=&(moldyn->lc); + + for(i=0;icount;i++) + virial+=v3_norm(&(atom[i].virial)); + +printf("%f\n",virial); + /* get pressure from virial */ + moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial; + moldyn->p/=moldyn->volume; +printf("%f\n",moldyn->p/(ATM)); + + /* scale factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) + scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc); + else + /* should actually never be used */ + scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0); + +printf("scale = %f\n",scale); + /* actual scaling */ + dim->x*=scale; + dim->y*=scale; + dim->z*=scale; + if(vdim->x) vdim->x=dim->x; + if(vdim->y) vdim->y=dim->y; + if(vdim->z) vdim->z=dim->z; + moldyn->volume*=(scale*scale*scale); + + /* check whether we need a new linkcell init */ + if((dim->x/moldyn->cutoff!=lc->nx)|| + (dim->y/moldyn->cutoff!=lc->ny)|| + (dim->z/moldyn->cutoff!=lc->nx)) { + link_cell_shutdown(moldyn); + link_cell_init(moldyn); + } + + return 0; + +} + double get_e_kin(t_moldyn *moldyn) { int i; @@ -724,6 +785,12 @@ int moldyn_integrate(t_moldyn *moldyn) { /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); + if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) +{ +printf("going to do p scale ...\n"); + scale_volume(moldyn); +printf("done\n"); +} /* check for log & visualization */ if(e) { @@ -853,13 +920,16 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; - + /* get energy and force of every atom */ for(i=0;ifunc1b(moldyn,&(itom[i])); @@ -1627,6 +1697,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ /* scaled with 0.5 ^ */ + } return 0; @@ -1655,6 +1726,7 @@ x=dim->x/2; printf("FATAL: atom %d: x: %.20f (%.20f)\n", i,atom[i].r.x,dim->x/2); printf("diagnostic:\n"); + printf("-----------\natom.r.x:\n"); for(j=0;j<8;j++) { memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1); for(k=0;k<8;k++) @@ -1662,7 +1734,7 @@ x=dim->x/2; ((byte)&(1<x/2; ((byte)&(1<=dim->y/2||-atom[i].r.y>dim->y/2) printf("FATAL: atom %d: y: %.20f (%.20f)\n", diff --git a/moldyn.h b/moldyn.h index 72ff30b..4203c79 100644 --- a/moldyn.h +++ b/moldyn.h @@ -27,6 +27,7 @@ typedef struct s_atom { t_3dvec r; /* position */ t_3dvec v; /* velocity */ t_3dvec f; /* force */ + t_3dvec virial; /* virial */ int element; /* number of element in pse */ double mass; /* atom mass */ u8 bnum; /* brand number */ @@ -67,6 +68,7 @@ typedef struct s_moldyn { t_atom *atom; /* pointer to the atoms */ t_3dvec dim; /* dimensions of the simulation volume */ + double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ /* potential force function and parameter pointers */ int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); @@ -92,7 +94,7 @@ typedef struct s_moldyn { double p; /* actual pressure */ /* pressure and temperature control (velocity/volume scaling) */ - /* (in units of tau!) */ + /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */ 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 */ @@ -129,13 +131,15 @@ typedef struct s_moldyn { int debug; /* debugging stuff, ignore */ } t_moldyn; -#define MOLDYN_STAT_PBX 0x08 /* periodic boudaries in x */ -#define MOLDYN_STAT_PBY 0x10 /* y */ -#define MOLDYN_STAT_PBZ 0x20 /* and z direction */ +#define MOLDYN_STAT_PBX 0x01 /* periodic boudaries in x */ +#define MOLDYN_STAT_PBY 0x02 /* y */ +#define MOLDYN_STAT_PBZ 0x04 /* and z direction */ -#define MOLDYN_1BP 0x00 /* care about single */ -#define MOLDYN_2BP 0x01 /* 2 body */ -#define MOLDYN_3BP 0x02 /* and 3 body particle pots */ +#define MOLDYN_PSCALE 0x08 /* size controlled by piston */ + +#define MOLDYN_1BP 0x10 /* care about single */ +#define MOLDYN_2BP 0x20 /* 2 body */ +#define MOLDYN_3BP 0x40 /* and 3 body particle pots */ #define T_SCALE_BERENDSEN 0x01 /* berendsen t control */ #define T_SCALE_DIRECT 0x02 /* direct t control */ @@ -246,12 +250,17 @@ typedef struct s_tersoff_mult_params { * */ +#define ONE_THIRD (1.0/3.0) + /* * default values * * - length unit: 1 A (1 A = 1e-10 m) * - time unit: 1 fs (1 fs = 1e-15 s) * - mass unit: 1 amu (1 amu = 1.6605388628e-27 kg ) + * + * fyi: in the following 1 N = (amu*A)/(fs*fs) + * */ #define METER 1e10 /* A */ @@ -259,6 +268,8 @@ typedef struct s_tersoff_mult_params { #define AMU 1.6605388628e-27 /* kg */ #define KILOGRAM (1.0/AMU) /* amu */ #define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) /* A amu / fs^2 */ +#define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */ +#define ATM (1.0133e5*PASCAL) /* N / A^2 */ #define MOLDYN_TEMP 273.0 #define MOLDYN_TAU 1.0 @@ -284,6 +295,7 @@ typedef struct s_tersoff_mult_params { * * phsical values / constants * + * */ #define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ @@ -349,6 +361,7 @@ 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_ref); +int set_pressure(t_moldyn *moldyn,double p_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); @@ -370,6 +383,7 @@ int destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); int scale_velocity(t_moldyn *moldyn,u8 equi_init); +int scale_volume(t_moldyn *moldyn); double get_e_kin(t_moldyn *moldyn); double get_e_pot(t_moldyn *moldyn); diff --git a/sic.c b/sic.c index bc10e31..516ee01 100644 --- a/sic.c +++ b/sic.c @@ -111,6 +111,7 @@ int main(int argc,char **argv) { 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); + moldyn_bc_check(&md); /* testing configuration */ //r.x=2.45/2; v.x=0; @@ -131,13 +132,20 @@ int main(int argc,char **argv) { /* set temperature */ printf("[sic] setting temperature\n"); + set_temperature(&md,273.0+1410.0); //set_temperature(&md,273.0+450.0); - set_temperature(&md,1.0); //set_temperature(&md,273.0); + //set_temperature(&md,1.0); + //set_temperature(&md,0.0); + + /* set pressure */ + printf("[sic] setting pressure\n"); + set_pressure(&md,ATM); /* set p/t scaling */ printf("[sic] set p/t scaling\n"); - set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, + // T_SCALE_BERENDSEN,100.0); /* initial thermal fluctuations of particles (in equilibrium) */ printf("[sic] thermal init\n"); @@ -145,13 +153,13 @@ int main(int argc,char **argv) { /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,1000,1.0); + moldyn_add_schedule(&md,100001,1.0); /* activate logging */ printf("[sic] activate logging\n"); moldyn_set_log_dir(&md,"saves/test"); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,10); - moldyn_set_log(&md,VISUAL_STEP,10); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,200); + moldyn_set_log(&md,VISUAL_STEP,200); /* * let's do the actual md algorithm now