X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=posic.c;h=7c1c41848aff58df3068542f0e6b7b9541bce733;hb=792f14f95b47989f7f12df0ea70b54619be016ee;hp=4eb2d043e4be7972848ef059b4f69ad39204ce6d;hpb=eceebe3ee412aa8cea3e6a7f0038883707f78460;p=physik%2Fposic.git diff --git a/posic.c b/posic.c index 4eb2d04..7c1c418 100644 --- a/posic.c +++ b/posic.c @@ -1,9 +1,11 @@ /* * posic.c - precipitation process of silicon carbide in silicon * - * author: Frank Zirkelbach + * author: Frank Zirkelbach * */ + +#include #include "moldyn.h" #include "math/math.h" @@ -14,6 +16,8 @@ int main(int argc,char **argv) { + t_moldyn md; + t_atom *si; t_visual vis; @@ -21,11 +25,14 @@ int main(int argc,char **argv) { t_random random; int a,b,c; - double t,e; + double t,e,u; + double help; t_3dvec p; int count; - char fb[32]="saves/fcc_test"; + t_lj_params lj; + + char fb[32]="saves/lj_test"; /* init */ @@ -43,28 +50,81 @@ int main(int argc,char **argv) { b=LEN_Y; c=LEN_Z; + /* set for 'bounding atoms' */ + //vis.dim.x=a*LC_SI; + //vis.dim.y=b*LC_SI; + //vis.dim.z=c*LC_SI; + t=TEMPERATURE; printf("placing silicon atoms ... "); - count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si); + count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si); printf("(%d) ok!\n",count); + /* testing purpose + count=2; + si=malloc(2*sizeof(t_atom)); + si[0].r.x=0.16e-9; + si[0].r.y=0; + si[0].r.z=0; + si[0].element=SI; + si[0].mass=M_SI; + si[1].r.x=-0.16e-9; + si[1].r.y=0; + si[1].r.z=0; + si[1].element=SI; + si[1].mass=M_SI; + */ + printf("setting thermal fluctuations\n"); thermal_init(si,&random,count,t); - - /* visualize */ - - visual_atoms(&vis,0.0,si,count); + //v3_zero(&(si[0].v)); + //v3_zero(&(si[1].v)); /* check kinetic energy */ e=get_e_kin(si,count); - printf("kinetic energy: %f\n",e); - printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t); + printf("kinetic energy: %.40f [J]\n",e); + printf("3/2 N k T = %.40f [J]\n",1.5*count*K_BOLTZMANN*t); /* check total momentum */ p=get_total_p(si,count); - printf("total momentum: %f\n",v3_norm(&p)); + printf("total momentum: %.30f [Ns]\n",v3_norm(&p)); + + /* check potential energy */ + md.count=count; + md.atom=si; + md.potential=potential_lennard_jones; + md.force=force_lennard_jones; + md.cutoff=R_CUTOFF; + md.cutoff_square=(R_CUTOFF*R_CUTOFF); + md.pot_params=&lj; + md.integrate=velocity_verlet; + md.time_steps=RUNS; + md.tau=TAU; + md.status=0; + md.visual=&vis; + md.write=WRITE_FILE; + + lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI; + help=lj.sigma6*lj.sigma6; + lj.sigma6*=help; + lj.sigma12=lj.sigma6*lj.sigma6; + lj.epsilon=LJ_EPSILON_SI; + + u=get_e_pot(&md); + + printf("potential energy: %.40f [J]\n",u); + printf("total energy (1): %.40f [J]\n",e+u); + printf("total energy (2): %.40f [J]\n",get_total_energy(&md)); + + md.dim.x=a*LC_SI; + md.dim.y=b*LC_SI; + md.dim.z=c*LC_SI; + + printf("estimated accurate time step: %.30f [s]\n", + estimate_time_step(&md,3.0,t)); + /* * let's do the actual md algorithm now @@ -72,139 +132,17 @@ int main(int argc,char **argv) { * integration of newtons equations */ + moldyn_integrate(&md); + + printf("total energy (after integration): %.40f [J]\n", + get_total_energy(&md)); + /* close */ visual_tini(&vis); rand_close(&random); - - //printf("starting velocity verlet: "); - //fflush(stdout); - - //for(runs=0;runsLX) si[i].x-=LEN_X; -// else if(si[i].x<-LX) si[i].x+=LEN_X; -// si[i].y+=(tau2*si[i].fy/m2); -// if(si[i].y>LY) si[i].y-=LEN_Y; -// else if(si[i].y<-LY) si[i].y+=LEN_Y; -// si[i].z+=(tau2*si[i].fz/m2); -// if(si[i].z>LZ) si[i].z-=LEN_Z; -// else if(si[i].z<-LZ) si[i].z+=LEN_Z; -// /* calculation of velocities v(t+h/2) */ -// si[i].vx+=(tau*si[i].fx/m2); -// si[i].vy+=(tau*si[i].fy/m2); -// si[i].vz+=(tau*si[i].fz/m2); -// /* reset of forces */ -// si[i].fx=.0; -// si[i].fy=.0; -// si[i].fz=.0; -// } -// for(i=0;iLX) deltax-=LEN_X; -// else if(-deltax>LX) deltax+=LEN_X; -// deltax2=deltax*deltax; -// deltay=si[i].y-si[j].y; -// if(deltay>LY) deltay-=LEN_Y; -// else if(-deltay>LY) deltay+=LEN_Y; -// deltay2=deltay*deltay; -// deltaz=si[i].z-si[j].z; -// if(deltaz>LZ) deltaz-=LEN_Z; -// else if(-deltaz>LZ) deltaz+=LEN_Z; -// deltaz2=deltaz*deltaz; -// distance=deltax2+deltay2+deltaz2; -// if(distance<=R2_CUTOFF) { -// tmp=1.0/distance; // 1/r^2 -// lj1=tmp; // 1/r^2 -// tmp*=tmp; // 1/r^4 -// lj1*=tmp; // 1/r^6 -// tmp*=tmp; // 1/r^8 -// lj2=tmp; // 1/r^8 -// lj1*=tmp; // 1/r^14 -// lj1*=LJ_SIGMA_12; -// lj2*=LJ_SIGMA_06; -// lj=-2*lj1+lj2; -// si[i].fx-=lj*deltax; -// si[i].fy-=lj*deltay; -// si[i].fz-=lj*deltaz; -// si[j].fx+=lj*deltax; -// si[j].fy+=lj*deltay; -// si[j].fz+=lj*deltaz; -// } -// } -// } -// for(i=0;i