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) */
_______________________________________________
ESPResSo mailing list
address@hidden
https://fias.uni-frankfurt.de/mailman/listinfo/espresso