v+=(virial->xx+virial->yy+virial->zz);
}
+ /* virial sum and mean virial */
+ moldyn->virial_sum+=v;
+ moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps;
+
/* assume up to date kinetic energy */
- moldyn->p=2.0*moldyn->ekin+v;
+ 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;
double thermodynamic_pressure_calc(t_moldyn *moldyn) {
t_3dvec dim,*tp;
- double u,p;
- double scale,dv;
+ double u_up,u_down,dv;
+ double scale,p;
t_atom *store;
/*
*
* => p = - dU/dV
*
- * dV: dx,y,z = 0.001 x,y,z
*/
- scale=1.00000000000001;
-printf("\n\nP-DEBUG:\n");
+ scale=0.00001;
+ dv=8*scale*scale*scale*moldyn->volume;
- tp=&(moldyn->tp);
store=malloc(moldyn->count*sizeof(t_atom));
if(store==NULL) {
printf("[moldyn] allocating store mem failed\n");
}
/* save unscaled potential energy + atom/dim configuration */
- u=moldyn->energy;
memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
dim=moldyn->dim;
- /* derivative with respect to x direction */
- scale_dim(moldyn,scale,TRUE,0,0);
- scale_atoms(moldyn,scale,TRUE,0,0);
- dv=0.00000000000001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z;
- link_cell_shutdown(moldyn);
- link_cell_init(moldyn,QUIET);
- potential_force_calc(moldyn);
- tp->x=(moldyn->energy-u)/dv;
- p=tp->x*tp->x;
-
- /* restore atomic configuration + dim */
- memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
- moldyn->dim=dim;
-
- /* derivative with respect to y direction */
- scale_dim(moldyn,scale,0,TRUE,0);
- scale_atoms(moldyn,scale,0,TRUE,0);
- dv=0.00000000000001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z;
+ /* scale up dimension and atom positions */
+ scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
+ scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
link_cell_shutdown(moldyn);
link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
- tp->y=(moldyn->energy-u)/dv;
- p+=tp->y*tp->y;
+ u_up=moldyn->energy;
/* restore atomic configuration + dim */
memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
moldyn->dim=dim;
- /* derivative with respect to z direction */
- scale_dim(moldyn,scale,0,0,TRUE);
- scale_atoms(moldyn,scale,0,0,TRUE);
- dv=0.00000000000001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y;
+ /* scale down dimension and atom positions */
+ scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
+ scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
link_cell_shutdown(moldyn);
link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
- tp->z=(moldyn->energy-u)/dv;
- p+=tp->z*tp->z;
+ u_down=moldyn->energy;
+
+ /* calculate pressure */
+ p=-(u_up-u_down)/dv;
+printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
/* restore atomic configuration + dim */
memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
moldyn->dim=dim;
/* restore energy */
- moldyn->energy=u;
+ potential_force_calc(moldyn);
link_cell_shutdown(moldyn);
link_cell_init(moldyn,QUIET);
- return sqrt(p);
+ return p;
}
double get_pressure(t_moldyn *moldyn) {
}
-int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
+int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
t_3dvec *dim;
dim=&(moldyn->dim);
+ if(dir==SCALE_UP)
+ scale=1.0+scale;
+
+ if(dir==SCALE_DOWN)
+ scale=1.0-scale;
+
if(x) dim->x*=scale;
if(y) dim->y*=scale;
if(z) dim->z*=scale;
return 0;
}
-int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
+int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
int i;
t_3dvec *r;
+ if(dir==SCALE_UP)
+ scale=1.0+scale;
+
+ if(dir==SCALE_DOWN)
+ scale=1.0-scale;
+
for(i=0;i<moldyn->count;i++) {
r=&(moldyn->atom[i].r);
if(x) r->x*=scale;
moldyn->debug=scale;
/* scale the atoms and dimensions */
- scale_atoms(moldyn,scale,TRUE,TRUE,TRUE);
- scale_dim(moldyn,scale,TRUE,TRUE,TRUE);
+ scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+ scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
/* visualize dimensions */
if(vdim->x!=0) {
temperature_calc(moldyn);
pressure_calc(moldyn);
//tp=thermodynamic_pressure_calc(moldyn);
-//printf("thermodynamic p: %f %f %f - %f\n",moldyn->tp.x/BAR,moldyn->tp.y/BAR,moldyn->tp.z/BAR,tp/BAR);
+//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR);
/* p/t scaling */
if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
/* single particle potential/force */
if(itom[i].attr&ATOM_ATTR_1BP)
- moldyn->func1b(moldyn,&(itom[i]));
+ if(moldyn->func1b)
+ moldyn->func1b(moldyn,&(itom[i]));
if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
continue;
double gp_sum; /* sum over all gp */
double mean_gp; /* mean value of gp */
+ double mean_v; /* mean of virial */
+ double virial_sum; /* sum over all calculated virials */
+
double p_ref; /* reference pressure */
double p; /* actual pressure (computed by virial) */
double p_sum; /* sum over all p */
#define VERBOSE 1
#define QUIET 0
+#define SCALE_UP 'u'
+#define SCALE_DOWN 'd'
+#define SCALE_DIRECT 'D'
+
/*
* potential related phsical values / constants
*
#define ONE_THIRD (1.0/3.0)
#define C 0x06
-//#define LC_C 3.567 /* A */
-#define LC_C 3.560 /* A */
+#define LC_C 3.567 /* A */
#define M_C 12.011 /* amu */
#define SI 0x0e
#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* A */
//#define LJ_SIGMA_SI (LC_SI/1.122462) /* A */
-//#define LJ_SIGMA_SI (0.5*sqrt(2.0)*LC_SI/1.122462) /* A */
+//#define LJ_SIGMA_SI (0.5*sqrt(2.0)*LC_SI/1.122462) /* A */
#define LJ_EPSILON_SI (2.1678*EV) /* NA */
#define TM_R_SI 2.7 /* A */
#define ALBE_D_SI 0.81472
#define ALBE_H_SI 0.259
+#define LC_SI_ALBE 5.429
+
#define ALBE_R_C (2.00-0.15)
#define ALBE_S_C (2.00+0.15)
#define ALBE_A_C (6.00*EV/1.167)
#define ALBE_D_C 6.28433
#define ALBE_H_C 0.5556
+#define LC_C_ALBE 3.566
+
#define ALBE_R_SIC (2.40-0.20)
#define ALBE_S_SIC (2.40+0.10)
#define ALBE_A_SIC (4.36*EV/0.847)
#define ALBE_D_SIC 180.314
#define ALBE_H_SIC 0.68
-#define ALBE_CHI_SIC 1.0
+#define LC_SIC_ALBE 4.359
+
/*
- * lattice constants
+ * lattice types
*/
#define CUBIC 0x01
double thermodynamic_pressure_calc(t_moldyn *moldyn);
double get_pressure(t_moldyn *moldyn);
int scale_volume(t_moldyn *moldyn);
-int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z);
-int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z);
+int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z);
+int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z);
double e_kin_calc(t_moldyn *moldyn);
double get_total_energy(t_moldyn *moldyn);
return 0;
}
-/* albe 1 body part */
-int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
-
- int brand;
- t_albe_mult_params *params;
- t_albe_exchange *exchange;
-
- brand=ai->brand;
- params=moldyn->pot_params;
- exchange=&(params->exchange);
-
- /*
- * simple: point constant parameters only depending on atom i to
- * their right values
- */
-
- exchange->gamma_i=&(params->gamma[brand]);
- exchange->c_i=&(params->c[brand]);
- exchange->d_i=&(params->d[brand]);
- exchange->h_i=&(params->h[brand]);
-
- exchange->ci2=params->c[brand]*params->c[brand];
- exchange->di2=params->d[brand]*params->d[brand];
- exchange->ci2di2=exchange->ci2/exchange->di2;
-
- return 0;
-}
-
/* albe 3 body potential function (first ij loop) */
int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
*/
brand=ai->brand;
-
- if(brand==aj->brand)
+ if(brand==aj->brand) {
S2=params->S2[brand];
- else
+ }
+ else {
S2=params->S2mixed;
+ }
/* dist_ij, d_ij2 */
v3_sub(&dist_ij,&(aj->r),&(ai->r));
R=params->R[brand];
S=params->S[brand];
S2=params->S2[brand];
+ /* albe needs i,k depending c,d,h and gamma values */
+ exchange->gamma_i=&(params->gamma[brand]);
+ exchange->c_i=&(params->c[brand]);
+ exchange->d_i=&(params->d[brand]);
+ exchange->h_i=&(params->h[brand]);
}
else {
R=params->Rmixed;
S=params->Smixed;
S2=params->S2mixed;
+ /* albe needs i,k depending c,d,h and gamma values */
+ exchange->gamma_i=&(params->gamma_m);
+ exchange->c_i=&(params->c_mixed);
+ exchange->d_i=&(params->d_mixed);
+ exchange->h_i=&(params->h_mixed);
}
+ exchange->ci2=*(exchange->c_i)**(exchange->c_i);
+ exchange->di2=*(exchange->d_i)**(exchange->d_i);
+ exchange->ci2di2=exchange->ci2/exchange->di2;
/* dist_ik, d_ik2 */
v3_sub(&dist_ik,&(ak->r),&(ai->r));
cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
/* g_ijk */
- h_cos=*(exchange->h_i)-cos_theta;
+ h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
d2_h_cos2=exchange->di2+(h_cos*h_cos);
frac=exchange->ci2/d2_h_cos2;
g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
- dg=-2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2;
+ dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
/* zeta sum += f_c_ik * g_ijk */
if(d_ik<=R) {
exchange->zeta_ij+=f_c_ik*g;
}
+#ifdef DEBUG
+ if((ai==&(moldyn->atom[0]))|
+ (aj==&(moldyn->atom[864]))|
+ (ak==&(moldyn->atom[1003]))) {
+ printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
+ }
+#endif
+
/* store even more data for second k loop */
exchange->g[kcount]=g;
exchange->dg[kcount]=dg;
}
/* force contribution */
- scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
+ scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*f_a));
v3_scale(&force,&(exchange->dist_ij),scale);
v3_add(&(ai->f),&(ai->f),&force);
v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
f_c,b,f_a,f_r);
+ printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
}
#endif
/* virial */
- if(aj<ai)
- virial_calc(ai,&force,&(exchange->dist_ij));
+ virial_calc(ai,&force,&(exchange->dist_ij));
/* dzeta prefactor = - 0.5 f_c f_a db */
exchange->pre_dzeta=-0.5*f_a*f_c*db;
/* virial */
//v3_scale(&force,&force,-1.0);
- if(aj<ai)
- virial_calc(ai,&force,&dist_ij);
+ //virial_calc(ai,&force,&dist_ij);
/* derivative wrt k */
v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
/* virial */
//v3_scale(&force,&force,-1.0);
- if(aj<ai)
- virial_calc(ai,&force,&dist_ik);
+ //virial_calc(ai,&force,&dist_ik);
/* increase k counter */
exchange->kcount++;
double mu_m; /* mixed mu */
double gamma[2];
+ double gamma_m;
double c[2];
+ double c_mixed;
double d[2];
+ double d_mixed;
double h[2];
+ double h_mixed;
t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */
} t_albe_mult_params;
/* function prototypes */
int albe_mult_complete_params(t_albe_mult_params *p);
-int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai);
-int albe_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
int albe_mult_3bp_k1(t_moldyn *moldyn,
t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
exchange->zeta_ij+=f_c_ik*g;
}
+#ifdef DEBUG
+ if((ai==&(moldyn->atom[0]))|
+ (aj==&(moldyn->atom[864]))|
+ (ak==&(moldyn->atom[1003]))) {
+ printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
+ }
+#endif
+
/* store even more data for second k loop */
exchange->g[kcount]=g;
exchange->dg[kcount]=dg;
printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
if(aj==&(moldyn->atom[0]))
printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+ printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
+ f_c,b,f_a,f_r);
+ printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
}
#endif
t_moldyn md;
/* potential parameters */
- t_lj_params lj;
- t_ho_params ho;
t_tersoff_mult_params tp;
t_albe_mult_params ap;
set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
/* choose potential */
- set_potential1b(&md,tersoff_mult_1bp);
-#ifdef TERSOFF_ORIG
- set_potential3b_j1(&md,tersoff_mult_2bp);
- set_potential3b_k1(&md,tersoff_mult_3bp);
- set_potential3b_j2(&md,tersoff_mult_post_2bp);
-#elif ALBE
- set_potential1b(&md,albe_mult_1bp);
+#ifdef ALBE
set_potential3b_j1(&md,albe_mult_3bp_j1);
set_potential3b_k1(&md,albe_mult_3bp_k1);
set_potential3b_j2(&md,albe_mult_3bp_j2);
set_potential3b_k2(&md,albe_mult_3bp_k2);
#else
+ set_potential1b(&md,tersoff_mult_1bp);
set_potential3b_j1(&md,tersoff_mult_3bp_j1);
set_potential3b_k1(&md,tersoff_mult_3bp_k1);
set_potential3b_j2(&md,tersoff_mult_3bp_j2);
set_potential3b_k2(&md,tersoff_mult_3bp_k2);
#endif
- //set_potential2b(&md,lennard_jones);
- //set_potential2b(&md,harmonic_oscillator);
#ifdef ALBE
set_potential_params(&md,&ap);
#else
set_potential_params(&md,&tp);
#endif
- //set_potential_params(&md,&lj);
- //set_potential_params(&md,&ho);
/* cutoff radius */
+#ifdef ALBE
set_cutoff(&md,ALBE_S_SI);
- //set_cutoff(&md,TM_S_SI);
- //set_cutoff(&md,LC_SI*sqrt(3.0));
- //set_cutoff(&md,2.0*LC_SI);
+#else
+ set_cutoff(&md,TM_S_SI);
+#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;
- 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.equilibrium_distance=LC_SI;
- ho.spring_constant=LJ_EPSILON_SI;
-
/*
* tersoff mult potential parameters for SiC
*/
ap.r0_mixed=ALBE_R0_SIC;
ap.lambda_m=ALBE_LAMBDA_SIC;
ap.mu_m=ALBE_MU_SIC;
+ ap.gamma_m=ALBE_GAMMA_SIC;
+ ap.c_mixed=ALBE_C_SIC;
+ ap.d_mixed=ALBE_D_SIC;
+ ap.h_mixed=ALBE_H_SIC;
albe_mult_complete_params(&ap);
/* set (initial) dimensions of simulation volume */
- set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+ //set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE);
+ //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+ //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE);
//set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
- //set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE);
+ set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE);
/* set periodic boundary conditions in all directions */
set_pbc(&md,TRUE,TRUE,TRUE);
/* create the lattice / place atoms */
- //create_lattice(&md,CUBIC,LC_SI,SI,M_SI,
- //create_lattice(&md,FCC,LC_SI,SI,M_SI,
- create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
- //create_lattice(&md,DIAMOND,LC_C,C,M_C,
- ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+ //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_ALBE,SI,M_SI,
+ // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
// ATOM_ATTR_2BP|ATOM_ATTR_HB,
- 0,6,6,6,NULL);
+ // 0,6,6,6,NULL);
// 1,6,6,6,NULL);
/* create centered zinc blende lattice */
- //r.x=0.5*0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
- //create_lattice(&md,FCC,TM_LC_3C_SIC,SI,M_SI,
- // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
- // 0,6,6,6,&r);
- //r.x+=0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
- //create_lattice(&md,FCC,TM_LC_3C_SIC,C,M_C,
- // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
- // 1,6,6,6,&r);
+ /**/
+ 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,
+ 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,
+ ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+ 1,6,6,6,&r);
+ /**/
moldyn_bc_check(&md);
/* create the simulation schedule */
/* initial configuration */
- moldyn_add_schedule(&md,100,1.0);
+ moldyn_add_schedule(&md,10000,1.0);
/* adding atoms */
//for(inject=0;inject<INJECT;inject++) {
// /* injecting atom and run with enabled t scaling */
moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
moldyn_set_log(&md,LOG_TEMPERATURE,1);
moldyn_set_log(&md,LOG_PRESSURE,1);
- moldyn_set_log(&md,VISUAL_STEP,1);
- moldyn_set_log(&md,SAVE_STEP,1);
+ moldyn_set_log(&md,VISUAL_STEP,100);
+ moldyn_set_log(&md,SAVE_STEP,100);
moldyn_set_log(&md,CREATE_REPORT,0);
/*
projects / physik / posic.git / commitdiff