- printf("placing silicon atoms\n");
- for(i=0;i<amount_si;i++) {
- si[i].x=RAND(LEN_X);
- si[i].y=RAND(LEN_Y);
- si[i].z=RAND(LEN_Z);
- si[i].vx=.0;
- si[i].vy=.0;
- si[i].vz=.0;
- si[i].fx=.0;
- si[i].fy=.0;
- si[i].fz=.0;
- }
-
- /* time */
- time=.0;
- tau=TAU;
- tau2=tau*tau;
-
- /* rasmol */
- printf("opening the rasmol file\n");
- fd=open("rasmol.xyz",O_WRONLY);
- if(fd<0) {
- perror("rasmol file open");
- return -1;
- }
-
- printf("starting velocity verlet: ");
- fflush(stdout);
-
- for(runs=0;runs<RUNS;runs++) {
-
- /*
- * velocity verlet
- *
- * r(t+h) = r(t) + h * dr/dt|t + h^2/2m * F(t)
- * dr/dt|(t+h) = dr/dt|t + h/2m * (F(t) + F(t+h))
+ moldyn_log_init(&md,&vis);
+ rand_init(&random,NULL,1);
+ random.status|=RAND_STAT_VERBOSE;
+
+ /* testing random numbers */
+ //for(a=0;a<1000000;a++)
+ // printf("%f %f\n",rand_get_gauss(&random),
+ // rand_get_gauss(&random));
+
+ a=LEN_X;
+ 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;
+
+ /* init lattice
+ printf("placing silicon atoms ... ");
+ md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
+ printf("(%d) ok!\n",md.count);
+ testing purpose */
+ md.count=2;
+ si=malloc(2*sizeof(t_atom));
+ si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
+ si[0].r.y=0;
+ si[0].r.z=0;
+ si[0].element=SI;
+ si[0].mass=M_SI;
+ si[1].r.x=-si[0].r.x;
+ si[1].r.y=0;
+ si[1].r.z=0;
+ si[1].element=SI;
+ si[1].mass=M_SI;
+ /* */
+
+ /* moldyn init (now si is a valid address) */
+ md.atom=si;
+ md.potential_force_function=lennard_jones;
+ //md.potential_force_function=harmonic_oscillator;
+ md.cutoff=R_CUTOFF*LC_SI;
+ md.pot_params=&lj;
+ //md.pot_params=&ho;
+ md.status=0;
+ md.visual=&vis;
+ /* dimensions of the simulation cell */
+ md.dim.x=a*LC_SI;
+ md.dim.y=b*LC_SI;
+ md.dim.z=c*LC_SI;
+
+ printf("setting thermal fluctuations (T=%f K)\n",md.t);
+ // thermal_init(&md,&random);
+ for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
+
+ /* check kinetic energy */
+
+ e=get_e_kin(si,md.count);
+ printf("kinetic energy: %.40f [J]\n",e);
+ printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
+ 1.5*md.count*K_BOLTZMANN*md.t,md.t);
+
+ /* check total momentum */
+ p=get_total_p(si,md.count);
+ printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
+
+ /* potential paramters */
+ lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
+ help=lj.sigma6*lj.sigma6;
+ lj.sigma6*=help;
+ lj.sigma12=lj.sigma6*lj.sigma6;
+ lj.epsilon4=4.0*LJ_EPSILON_SI;
+
+ ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
+ ho.spring_constant=1.0;
+
+ printf("estimated accurate time step: %.30f [s]\n",
+ estimate_time_step(&md,3.0,md.t));
+
+ /*
+ * let's do the actual md algorithm now