+ /* shift partition values */
+ if(p_params->type) {
+ p.x=p_params->p.x+(a*lc)/2.0;
+ p.y=p_params->p.y+(b*lc)/2.0;
+ p.z=p_params->p.z+(c*lc)/2.0;
+ }
+
+ r.x=o.x;
+ for(i=0;i<a;i++) {
+ r.y=o.y;
+ for(j=0;j<b;j++) {
+ r.z=o.z;
+ for(k=0;k<c;k++) {
+ switch(p_params->type) {
+ case PART_INSIDE_R:
+ v3_sub(&dist,&r,&p);
+ if(v3_absolute_square(&dist)<
+ (p_params->r*p_params->r)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_OUTSIDE_R:
+ v3_sub(&dist,&r,&p);
+ if(v3_absolute_square(&dist)>=
+ (p_params->r*p_params->r)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_INSIDE_D:
+ v3_sub(&dist,&r,&p);
+ if((fabs(dist.x)<p_params->d.x)&&
+ (fabs(dist.y)<p_params->d.y)&&
+ (fabs(dist.z)<p_params->d.z)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_OUTSIDE_D:
+ v3_sub(&dist,&r,&p);
+ if((fabs(dist.x)>=p_params->d.x)||
+ (fabs(dist.y)>=p_params->d.y)||
+ (fabs(dist.z)>=p_params->d.z)) {
+ set_atom_positions(r);
+ }
+ break;
+ default:
+ set_atom_positions(r);
+ break;
+ }
+ r.z+=lc;
+ }
+ r.y+=lc;
+ }
+ r.x+=lc;
+ }
+
+ for(i=0;i<count;i++) {
+ atom[i].r.x-=(a*lc)/2.0;
+ atom[i].r.y-=(b*lc)/2.0;
+ atom[i].r.z-=(c*lc)/2.0;
+ }
+
+ return count;
+}
+
+/* fcc lattice init */
+int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+ t_part_params *p_params,t_defect_params *d_params) {
+
+ int count;
+ int i,j,k,l;
+ t_3dvec o,r,n;
+ t_3dvec basis[3];
+ t_3dvec dist;
+ t_3dvec p;
+ t_3dvec d_d,d_o,dr;
+
+ p.x=0; p.y=0; p.z=0;
+
+ count=0;
+ if(origin)
+ v3_copy(&o,origin);
+ else
+ v3_zero(&o);
+
+ /* construct the basis */
+ memset(basis,0,3*sizeof(t_3dvec));
+ basis[0].x=0.5*lc;
+ basis[0].y=0.5*lc;
+ basis[1].x=0.5*lc;
+ basis[1].z=0.5*lc;
+ basis[2].y=0.5*lc;
+ basis[2].z=0.5*lc;
+
+ /* shift partition values */
+ if(p_params->type) {
+ p.x=p_params->p.x+(a*lc)/2.0;
+ p.y=p_params->p.y+(b*lc)/2.0;
+ p.z=p_params->p.z+(c*lc)/2.0;
+ }
+
+ /* fill up the room */
+ r.x=o.x;
+ for(i=0;i<a;i++) {
+ r.y=o.y;
+ for(j=0;j<b;j++) {
+ r.z=o.z;
+ for(k=0;k<c;k++) {
+ /* first atom */
+ switch(p_params->type) {
+ case PART_INSIDE_R:
+ v3_sub(&dist,&r,&p);
+ if(v3_absolute_square(&dist)<
+ (p_params->r*p_params->r)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_OUTSIDE_R:
+ v3_sub(&dist,&r,&p);
+ if(v3_absolute_square(&dist)>=
+ (p_params->r*p_params->r)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_INSIDE_D:
+ v3_sub(&dist,&r,&p);
+ if((fabs(dist.x)<p_params->d.x)&&
+ (fabs(dist.y)<p_params->d.y)&&
+ (fabs(dist.z)<p_params->d.z)) {
+ set_atom_positions(r);
+ }
+ break;
+ case PART_OUTSIDE_D:
+ v3_sub(&dist,&r,&p);
+ if((fabs(dist.x)>=p_params->d.x)||
+ (fabs(dist.y)>=p_params->d.y)||
+ (fabs(dist.z)>=p_params->d.z)) {
+ set_atom_positions(r);
+ }
+ break;
+ default:
+ set_atom_positions(r);
+ break;
+ }
+ /* the three face centered atoms */
+ for(l=0;l<3;l++) {
+ v3_add(&n,&r,&basis[l]);
+ switch(p_params->type) {
+ case PART_INSIDE_R:
+ v3_sub(&dist,&n,&p);
+ if(v3_absolute_square(&dist)<
+ (p_params->r*p_params->r)) {
+ set_atom_positions(n);
+ }
+ break;
+ case PART_OUTSIDE_R:
+ v3_sub(&dist,&n,&p);
+ if(v3_absolute_square(&dist)>=
+ (p_params->r*p_params->r)) {
+ set_atom_positions(n);
+ }
+ break;
+ case PART_INSIDE_D:
+ v3_sub(&dist,&n,&p);
+ if((fabs(dist.x)<p_params->d.x)&&
+ (fabs(dist.y)<p_params->d.y)&&
+ (fabs(dist.z)<p_params->d.z)) {
+ set_atom_positions(n);
+ }
+ break;
+ case PART_OUTSIDE_D:
+ v3_sub(&dist,&n,&p);
+ if((fabs(dist.x)>=p_params->d.x)||
+ (fabs(dist.y)>=p_params->d.y)||
+ (fabs(dist.z)>=p_params->d.z)) {
+ set_atom_positions(n);
+ }
+ break;
+ default:
+ set_atom_positions(n);
+ break;
+ }
+ }
+ r.z+=lc;
+ }
+ r.y+=lc;
+ }
+ r.x+=lc;
+ }
+
+ /* coordinate transformation */
+ for(i=0;i<count;i++) {
+ atom[i].r.x-=(a*lc)/2.0;
+ atom[i].r.y-=(b*lc)/2.0;
+ atom[i].r.z-=(c*lc)/2.0;
+ }
+
+ return count;
+}
+
+int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+ t_part_params *p_params,t_defect_params *d_params) {
+
+ int count;
+ t_3dvec o;
+
+ count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params);
+
+ o.x=0.25*lc;
+ o.y=0.25*lc;
+ o.z=0.25*lc;
+
+ if(origin) v3_add(&o,&o,origin);
+
+ count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params);
+
+ return count;
+}
+
+int destroy_atoms(t_moldyn *moldyn) {
+
+ if(moldyn->atom) free(moldyn->atom);
+
+ return 0;
+}
+
+int thermal_init(t_moldyn *moldyn,u8 equi_init) {
+
+ /*
+ * - gaussian distribution of velocities
+ * - zero total momentum
+ * - velocity scaling (E = 3/2 N k T), E: kinetic energy
+ */
+
+ int i;
+ double v,sigma;
+ t_3dvec p_total,delta;
+ t_atom *atom;
+ t_random *random;
+
+ atom=moldyn->atom;
+ random=&(moldyn->random);
+
+ printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
+
+ /* gaussian distribution of velocities */
+ v3_zero(&p_total);
+ for(i=0;i<moldyn->count;i++) {
+ sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
+ /* x direction */
+ v=sigma*rand_get_gauss(random);
+ atom[i].v.x=v;
+ p_total.x+=atom[i].mass*v;
+ /* y direction */
+ v=sigma*rand_get_gauss(random);
+ atom[i].v.y=v;
+ p_total.y+=atom[i].mass*v;
+ /* z direction */
+ v=sigma*rand_get_gauss(random);
+ atom[i].v.z=v;
+ p_total.z+=atom[i].mass*v;
+ }
+
+ /* zero total momentum */
+ v3_scale(&p_total,&p_total,1.0/moldyn->count);
+ for(i=0;i<moldyn->count;i++) {
+ v3_scale(&delta,&p_total,1.0/atom[i].mass);
+ v3_sub(&(atom[i].v),&(atom[i].v),&delta);
+ }
+
+ /* velocity scaling */
+ scale_velocity(moldyn,equi_init);
+
+ return 0;
+}
+
+double total_mass_calc(t_moldyn *moldyn) {
+
+ int i;
+
+ moldyn->mass=0.0;
+
+ for(i=0;i<moldyn->count;i++)
+ moldyn->mass+=moldyn->atom[i].mass;
+
+ return moldyn->mass;
+}
+
+double temperature_calc(t_moldyn *moldyn) {
+
+ /* assume up to date kinetic energy, which is 3/2 N k_B T */
+
+ if(moldyn->count)
+ moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+ else moldyn->t=0.0;
+
+ return moldyn->t;
+}
+
+double get_temperature(t_moldyn *moldyn) {
+
+ return moldyn->t;
+}