ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Commit2 JT 092418

Browse Source 
- initialized g_ewald before Newton solver
- mu2 is now adim in ewald_dipole_spin
This commit is contained in:
julient31
2018-09-24 16:40:59 -06:00
parent 53a779067e
commit 6b4303c405
4 changed files with 23 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
examples/SPIN/pppm_spin/in.dipole.pppm_dipole
View File

@ -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

4
examples/SPIN/pppm_spin/in.spin.ewald_spin
View File

@ -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

9
src/KSPACE/ewald_dipole.cpp
View File

@ -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;

12
src/KSPACE/ewald_dipole_spin.cpp
View File

@ -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)