[ESPResSo] magnetic dipolar binary mixtures

[Valera, Manuel]

Dear users,

Hello to all,

I am trying to simulate a system with binary mixtures of magnetic dipolar 
particles having dopoles mu1 and mu2
However, setting the field dipm for the particle does not seem to have any 
effect.
I have attached the function add_p3m_dipolar_pair_force, the values dipm are 
only used if NPT is enable, lines 349-354 below.
Is this a feature? Is there a way to simulate binary dipolar systems?
Thanks in advanced.

Manuel


00286 MDINLINE double add_p3m_dipolar_pair_force(Particle *p1, Particle *p2,
00287                                            double *d,double dist2,double 
dist,double force[3])
00288 {
00289   int j;
00290 #ifdef NPT
00291   double fac1;
00292 #endif
00293   double adist, erfc_part_ri, coeff, exp_adist2, dist2i;
00294   double mimj, mir, mjr;
00295   double B_r, C_r, D_r;
00296   double alpsq = p3m.Dalpha * p3m.Dalpha;
00297   double mixmj[3], mixr[3], mjxr[3];
00298 
00299   if(dist < p3m.Dr_cut && dist > 0) {
00300     adist = p3m.Dalpha * dist;
00301     #if USE_ERFC_APPROXIMATION
00302        erfc_part_ri = AS_erfc_part(adist) / dist;
00303     #else
00304        erfc_part_ri = erfc(adist) / dist;
00305     #endif
00306 
00307   //Calculate scalar multiplications for vectors mi, mj, rij
00308   mimj = p1->r.dip[0]*p2->r.dip[0] + p1->r.dip[1]*p2->r.dip[1] + 
p1->r.dip[2]*p2->r.dip[2];
00309   mir = p1->r.dip[0]*d[0] + p1->r.dip[1]*d[1] + p1->r.dip[2]*d[2];
00310   mjr = p2->r.dip[0]*d[0] + p2->r.dip[1]*d[1] + p2->r.dip[2]*d[2];
00311 
00312   coeff = 2.0*p3m.Dalpha*wupii;
00313   dist2i = 1 / dist2;
00314   exp_adist2 = exp(-adist*adist);
00315 
00316   if(p3m.Daccuracy > 5e-06)
00317     B_r = (erfc_part_ri + coeff) * exp_adist2 * dist2i;
00318   else
00319     B_r = (erfc(adist)/dist + coeff * exp_adist2) * dist2i;
00320     
00321   C_r = (3*B_r + 2*alpsq*coeff*exp_adist2) * dist2i;
00322   D_r = (5*C_r + 4*coeff*alpsq*alpsq*exp_adist2) * dist2i;
00323 
00324   // Calculate real-space forces
00325   for(j=0;j<3;j++)
00326     force[j] += coulomb.Dprefactor *((mimj*d[j] + p1->r.dip[j]*mjr + 
p2->r.dip[j]*mir) * C_r - mir*mjr*D_r*d[j]) ;
00327 
00328   //Calculate vector multiplications for vectors mi, mj, rij
00329 
00330   mixmj[0] = p1->r.dip[1]*p2->r.dip[2] - p1->r.dip[2]*p2->r.dip[1];
00331   mixmj[1] = p1->r.dip[2]*p2->r.dip[0] - p1->r.dip[0]*p2->r.dip[2];
00332   mixmj[2] = p1->r.dip[0]*p2->r.dip[1] - p1->r.dip[1]*p2->r.dip[0];
00333 
00334   mixr[0] = p1->r.dip[1]*d[2] - p1->r.dip[2]*d[1];
00335   mixr[1] = p1->r.dip[2]*d[0] - p1->r.dip[0]*d[2];
00336   mixr[2] = p1->r.dip[0]*d[1] - p1->r.dip[1]*d[0];
00337 
00338   mjxr[0] = p2->r.dip[1]*d[2] - p2->r.dip[2]*d[1];
00339   mjxr[1] = p2->r.dip[2]*d[0] - p2->r.dip[0]*d[2];
00340   mjxr[2] = p2->r.dip[0]*d[1] - p2->r.dip[1]*d[0];
00341 
00342   // Calculate real-space torques
00343 #ifdef ROTATION
00344   for(j=0;j<3;j++){
00345     p1->f.torque[j] += coulomb.Dprefactor *(-mixmj[j]*B_r + 
mixr[j]*mjr*C_r);
00346     p2->f.torque[j] += coulomb.Dprefactor *( mixmj[j]*B_r + 
mjxr[j]*mir*C_r);
00347   }
00348 #endif
00349 #ifdef NPT
00350 #if USE_ERFC_APPROXIMATION
00351   fac1 = coulomb.Dprefactor * p1->p.dipm*p2->p.dipm * exp(-adist*adist);
00352 #else
00353   fac1 = coulomb.Dprefactor * p1->p.dipm*p2->p.dipm;
00354 #endif
00355   return fac1 * ( mimj*B_r - mir*mjr * C_r );
00356 #endif
00357   }
00358   return 0.0;
00359 }
00360 
00362 MDINLINE double p3m_dipolar_pair_energy(Particle *p1, Particle *p2,
00363                                       double *d,double dist2,double dist)
00364 {
00365   double /* fac1,*/ adist, erfc_part_ri, coeff, exp_adist2, dist2i;
00366   double mimj, mir, mjr;
00367   double B_r, C_r;
00368   double alpsq = p3m.Dalpha * p3m.Dalpha;
00369  
00370   if(dist < p3m.Dr_cut && dist > 0) {
00371     adist = p3m.Dalpha * dist;
00372     /*fac1 = coulomb.Dprefactor;*/
00373 
00374 #if USE_ERFC_APPROXIMATION
00375     erfc_part_ri = AS_erfc_part(adist) / dist;
00376     /*  fac1 = coulomb.Dprefactor * p1->p.dipm*p2->p.dipm; IT WAS WRONG 
*/ /* *exp(-adist*adist); */ 
00377 #else
00378     erfc_part_ri = erfc(adist) / dist;
00379     /* fac1 = coulomb.Dprefactor * p1->p.dipm*p2->p.dipm;  IT WAS WRONG*/
00380 #endif
00381 
00382     //Calculate scalar multiplications for vectors mi, mj, rij
00383     mimj = p1->r.dip[0]*p2->r.dip[0] + p1->r.dip[1]*p2->r.dip[1] + 
p1->r.dip[2]*p2->r.dip[2];
00384     mir = p1->r.dip[0]*d[0] + p1->r.dip[1]*d[1] + p1->r.dip[2]*d[2];
00385     mjr = p2->r.dip[0]*d[0] + p2->r.dip[1]*d[1] + p2->r.dip[2]*d[2];
00386 
00387     coeff = 2.0*p3m.Dalpha*wupii;
00388     dist2i = 1 / dist2;
00389     exp_adist2 = exp(-adist*adist);
00390 
00391     if(p3m.Daccuracy > 5e-06)
00392       B_r = (erfc_part_ri + coeff) * exp_adist2 * dist2i;
00393     else
00394       B_r = (erfc(adist)/dist + coeff * exp_adist2) * dist2i;
00395   
00396     C_r = (3*B_r + 2*alpsq*coeff*exp_adist2) * dist2i;
00397 
00398     /*
00399       printf("(%4i %4i) pair energy = %f (B_r=%15.12f 
C_r=%15.12f)\n",p1->p.identity,p2->p.identity,fac1*(mimj*B_r-mir*mjr*C_r),B_r,C_r);
00400     */
00401   
00402     /* old line return fac1 * ( mimj*B_r - mir*mjr * C_r );*/
00403     return coulomb.Dprefactor * ( mimj*B_r - mir*mjr * C_r );
00404   }
00405   return 0.0;
00406 }
00407 
00409 #endif /* of defined(MAGNETOSTATICS) */