diff --git a/examples/SPIN/pppm_spin/in.dipole.pppm_dipole b/examples/SPIN/pppm_spin/in.dipole.pppm_dipole index c2c49e3caf..86ac5198b0 100644 --- a/examples/SPIN/pppm_spin/in.dipole.pppm_dipole +++ b/examples/SPIN/pppm_spin/in.dipole.pppm_dipole @@ -28,7 +28,7 @@ pair_coeff * * 0.0 0.0 #kspace_style pppm/disp 1.0e-4 #kspace_style pppm/dipole 1.0e-4 kspace_style ewald/dipole 1.0e-4 -kspace_modify gewald 0.1 +#kspace_modify gewald 0.1 neighbor 0.3 bin neigh_modify every 2 delay 4 check yes diff --git a/examples/SPIN/pppm_spin/in.spin.ewald_spin b/examples/SPIN/pppm_spin/in.spin.ewald_spin index d9ca46b27a..f2192676c7 100644 --- a/examples/SPIN/pppm_spin/in.spin.ewald_spin +++ b/examples/SPIN/pppm_spin/in.spin.ewald_spin @@ -36,8 +36,8 @@ neigh_modify every 10 check yes delay 20 #kspace_style pppm/dipole/spin 1.0e-4 #kspace_style ewald/dipole/spin 1.0e-4 -kspace_style ewald/dipole/spin 1.0e-2 -kspace_modify mesh 32 32 32 +kspace_style ewald/dipole/spin 1.0e-4 +#kspace_modify mesh 32 32 32 #fix 1 all precession/spin zeeman 1.0 0.0 0.0 1.0 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 diff --git a/src/KSPACE/ewald_dipole.cpp b/src/KSPACE/ewald_dipole.cpp index 42d850dcc2..c3a3818013 100644 --- a/src/KSPACE/ewald_dipole.cpp +++ b/src/KSPACE/ewald_dipole.cpp @@ -163,6 +163,15 @@ void EwaldDipole::init() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); + + // initial guess with old method + + g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*mu2); + if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; + else g_ewald = sqrt(-log(g_ewald)) / cutoff; + + // try Newton solver + double g_ewald_new = NewtonSolve(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2); if (g_ewald_new > 0.0) g_ewald = g_ewald_new; diff --git a/src/KSPACE/ewald_dipole_spin.cpp b/src/KSPACE/ewald_dipole_spin.cpp index 5522b18e03..43b9b32c76 100644 --- a/src/KSPACE/ewald_dipole_spin.cpp +++ b/src/KSPACE/ewald_dipole_spin.cpp @@ -164,6 +164,15 @@ void EwaldDipoleSpin::init() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); + + // initial guess with old method + + g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*mu2); + if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; + else g_ewald = sqrt(-log(g_ewald)) / cutoff; + + // try Newton solver + double g_ewald_new = NewtonSolve(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2); if (g_ewald_new > 0.0) g_ewald = g_ewald_new; @@ -887,7 +896,8 @@ void EwaldDipoleSpin::spsum_musq() MPI_Allreduce(&musum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(&musqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world); - mu2 = musqsum * mub2mu0; + //mu2 = musqsum * mub2mu0; + mu2 = musqsum; } if (mu2 == 0 && comm->me == 0)