t_atom *atom;
t_3dvec *dim,*vdim;
- double virial,scale;
+ double scale,v;
+ t_virial virial;
t_linkcell *lc;
int i;
vdim=&(moldyn->vis.dim);
lc=&(moldyn->lc);
- for(i=0;i<moldyn->count;i++)
- virial+=v3_norm(&(atom[i].virial));
+ memset(&virial,0,sizeof(t_virial));
+
+ for(i=0;i<moldyn->count;i++) {
+ virial.xx+=atom[i].virial.xx;
+ virial.yy+=atom[i].virial.yy;
+ virial.zz+=atom[i].virial.zz;
+ virial.xy+=atom[i].virial.xy;
+ virial.xz+=atom[i].virial.xz;
+ virial.yz+=atom[i].virial.yz;
+ }
+
+ /* just a guess so far ... */
+ v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz);
-printf("%f\n",virial);
+printf("%f\n",v);
/* get pressure from virial */
- moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+ moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v;
moldyn->p/=moldyn->volume;
printf("%f\n",moldyn->p/(ATM));
int i,j,k,count;
t_atom *itom,*jtom,*ktom;
+ t_virial *virial;
t_linkcell *lc;
t_list neighbour_i[27];
t_list neighbour_i2[27];
v3_zero(&(itom[i].f));
/* reset viral of atom i */
- v3_zero(&(itom[i].virial));
+ virial=&(itom[i].virial);
+ virial->xx=0.0;
+ virial->yy=0.0;
+ virial->zz=0.0;
+ virial->xy=0.0;
+ virial->xz=0.0;
+ virial->yz=0.0;
/* reset site energy */
itom[i].e=0.0;
int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
printf("[moldyn] tersoff parameter completion\n");
+ p->S2[0]=p->S[0]*p->S[0];
+ p->S2[1]=p->S[1]*p->S[1];
p->Smixed=sqrt(p->S[0]*p->S[1]);
+ p->S2mixed=p->Smixed*p->Smixed;
p->Rmixed=sqrt(p->R[0]*p->R[1]);
p->Amixed=sqrt(p->A[0]*p->A[1]);
p->Bmixed=sqrt(p->B[0]*p->B[1]);
t_tersoff_mult_params *params;
t_tersoff_exchange *exchange;
t_3dvec dist_ij,force;
- double d_ij;
- double A,B,R,S,lambda,mu;
+ double d_ij,d_ij2;
+ double A,B,R,S,S2,lambda,mu;
double f_r,df_r;
double f_c,df_c;
int brand;
*
*/
- /* dist_ij, d_ij */
- v3_sub(&dist_ij,&(aj->r),&(ai->r));
- if(bc) check_per_bound(moldyn,&dist_ij);
- d_ij=v3_norm(&dist_ij);
-
- /* save for use in 3bp */
- exchange->d_ij=d_ij;
- exchange->dist_ij=dist_ij;
-
/* constants */
if(brand==ai->brand) {
S=params->S[brand];
+ S2=params->S2[brand];
R=params->R[brand];
A=params->A[brand];
B=params->B[brand];
}
else {
S=params->Smixed;
+ S2=params->S2mixed;
R=params->Rmixed;
A=params->Amixed;
B=params->Bmixed;
params->exchange.chi=params->chi;
}
- /* if d_ij > S => no force & potential energy contribution */
- if(d_ij>S)
+ /* dist_ij, d_ij */
+ v3_sub(&dist_ij,&(aj->r),&(ai->r));
+ if(bc) check_per_bound(moldyn,&dist_ij);
+ d_ij2=v3_absolute_square(&dist_ij);
+
+ /* if d_ij2 > S2 => no force & potential energy contribution */
+ if(d_ij2>S2)
return 0;
+ /* now we will need the distance */
+ //d_ij=v3_norm(&dist_ij);
+ d_ij=sqrt(d_ij2);
+
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
+ exchange->dist_ij=dist_ij;
+
/* more constants */
exchange->beta_j=&(params->beta[brand]);
exchange->n_j=&(params->n[brand]);
s_r=S-R;
arg=M_PI*(d_ij-R)/s_r;
f_c=0.5+0.5*cos(arg);
- //df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* MARK! */
df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
/* two body contribution (ij, ji) */
v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
s_r=S-R;
arg=M_PI*(d_ik-R)/s_r;
f_c_ik=0.5+0.5*cos(arg);
- //df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
/* zeta_ij */
/* general */
typedef unsigned char u8;
+/* virial */
+typedef struct s_virial {
+ double xx; /* | xx xy xz | */
+ double yy; /* V = | yx yy yz | */
+ double zz; /* | zx zy zz | */
+ double xy; /* */
+ double xz; /* with: xy=yx, xz=zx, yz=zy */
+ double yz; /* */
+} t_virial;
+
/* the atom of the md simulation */
typedef struct s_atom {
t_3dvec r; /* position */
t_3dvec v; /* velocity */
t_3dvec f; /* force */
- t_3dvec virial; /* virial (v_xx, v_yy, v_zz) */
+ t_virial virial; /* virial */
double e; /* site energy */
int element; /* number of element in pse */
double mass; /* atom mass */
/* tersoff multi (2!) potential parameters */
typedef struct s_tersoff_mult_params {
double S[2]; /* tersoff cutoff radii */
+ double S2[2]; /* tersoff cutoff radii squared */
double R[2]; /* tersoff cutoff radii */
double Smixed; /* mixed S radius */
+ double S2mixed; /* mixed S radius squared */
double Rmixed; /* mixed R radius */
double A[2]; /* factor of tersoff attractive part */
double B[2]; /* factor of tersoff repulsive part */
#define SI 0x0e
#define LC_SI (0.543105e-9*METER) /* A */
#define M_SI 28.08553 /* amu */
+
#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* A */
#define LJ_EPSILON_SI (2.1678*EV) /* NA */
/* 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);
- moldyn_bc_check(&md);
+ //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.95/2; 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,
- // &r,&v);
- //r.x=-2.95/2; 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,
- // &r,&v);
+ r.x=2.8/2; 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,
+ &r,&v);
+ r.x=-2.8/2; 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,
+ &r,&v);
/* setting a nearest neighbour distance for the moldyn checks */
set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
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);
/* 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,20000,.1);
- moldyn_add_schedule(&md,10000,.2);
- moldyn_add_schedule(&md,6667,.3);
- moldyn_add_schedule(&md,5000,.4);
- moldyn_add_schedule(&md,4001,.5);
+ moldyn_add_schedule(&md,10000,.1);
/* activate logging */
printf("[sic] activate logging\n");
moldyn_set_log_dir(&md,argv[1]);
- moldyn_set_log(&md,LOG_TOTAL_ENERGY,100);
- moldyn_set_log(&md,VISUAL_STEP,100);
+ moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
+ moldyn_set_log(&md,VISUAL_STEP,20);
/*
* let's do the actual md algorithm now
projects / physik / posic.git / commitdiff