+int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ t_stage *stage;
+ t_insert_mixed_atoms_params *imp;
+ t_atom *atom;
+ double x,x0,y,y0,z,z0;
+ double dmin,d,c;
+ u8 retry;
+ t_3dvec r,v,dist;
+ int i;
+
+
+ stage=mdrun->stage.current->data;
+ imp=stage->params;
+
+ x=moldyn->dim.x;
+ x0=-x/2.0;
+ y=moldyn->dim.y;
+ y0=-y/2.0;
+ z=moldyn->dim.z;
+ z0=-z/2.0;
+
+ v.x=0.0;
+ v.y=0.0;
+ v.z=0.0;
+
+ c=imp->cr*imp->cr;
+
+ while(imp->amount1|imp->amount2) {
+ if(imp->amount1) {
+ retry=1;
+ while(retry) {
+ retry=0;
+ r.x=rand_get_double(&(moldyn->random))*x;
+ r.y=rand_get_double(&(moldyn->random))*y;
+ r.z=rand_get_double(&(moldyn->random))*z;
+ r.x+=x0;
+ r.y+=y0;
+ r.z+=z0;
+ dmin=1000.0; // for sure too high!
+ for(i=0;i<moldyn->count;i++) {
+ atom=&(moldyn->atom[i]);
+ v3_sub(&dist,&(atom->r),&r);
+ check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
+ if(d<c) {
+ retry=1;
+ break;
+ }
+ if(d<dmin)
+ dmin=d;
+
+ }
+ }
+ add_atom(moldyn,imp->element1,pse_mass[imp->element1],
+ imp->brand1,imp->attr1,&r,&v);
+ printf("%s (mixed) atom inserted (%d): %f %f %f\n",
+ ME,imp->amount1,r.x,r.y,r.z);
+ printf(" -> d2 = %f/%f\n",dmin,c);
+ imp->amount1-=1;
+ }
+ if(imp->amount2) {
+ retry=1;
+ while(retry) {
+ retry=0;
+ r.x=rand_get_double(&(moldyn->random))*x;
+ r.y=rand_get_double(&(moldyn->random))*y;
+ r.z=rand_get_double(&(moldyn->random))*z;
+ r.x+=x0;
+ r.y+=y0;
+ r.z+=z0;
+ dmin=1000.0; // for sure too high!
+ for(i=0;i<moldyn->count;i++) {
+ atom=&(moldyn->atom[i]);
+ v3_sub(&dist,&(atom->r),&r);
+ check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
+ if(d<c) {
+ retry=1;
+ break;
+ }
+ if(d<dmin)
+ dmin=d;
+
+ }
+ }
+ add_atom(moldyn,imp->element2,pse_mass[imp->element2],
+ imp->brand2,imp->attr2,&r,&v);
+ printf("%s (mixed) atom inserted (%d): %f %f %f\n",
+ ME,imp->amount2,r.x,r.y,r.z);
+ printf(" -> d2 = %f/%f\n",dmin,c);
+ imp->amount2-=1;
+ }
+ }
+
+ return 0;
+}
+