jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Commit JT 082318 

- corrected memory errors in pppm_dipole and pppm_dipole_spin
- created fm_long in atom_vec_spin
- fm_long added to fm in initial_integrate (in ComputeInteractionsSpin)
This commit is contained in:
julient31 2018-08-23 15:18:30 -06:00
parent 8d79db03d3
commit cf1d421e10
17 changed files with 288 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy Symbolic link
View File

@ -0,0 +1 @@
../cobalt_fcc/Co_PurjaPun_2012.eam.alloy

32
examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_fit.py Normal file
View File

@ -0,0 +1,32 @@
#Program fitting the exchange interaction
#Model curve: Bethe-Slater function
import numpy as np, pylab, tkinter
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from decimal import *
print("Loop begin")
#Definition of the Bethe-Slater function
def func(x,a,b,c):
return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)
#Exchange coeff table (data to fit)
rdata, Jdata = np.loadtxt('exchange_hcp_co.dat', usecols=(0,1), unpack=True)
plt.plot(rdata, Jdata, 'b-', label='data')
#Perform the fit
popt, pcov = curve_fit(func, rdata, Jdata, bounds=(0, [500.,5.,5.]))
plt.plot(rdata, func(rdata, *popt), 'r--', label='fit')
#Print the fitted params
print("Parameters: a={:.10} (in meV), b={:.10} (adim), c={:.10} (in Ang)".format(*popt))
#Ploting the result
plt.xlabel('r_ij')
pylab.xlim([0,6.5])
plt.ylabel('J_ij')
plt.legend()
plt.show()
print("Loop end")

9
examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_hcp_co.dat Normal file
View File

@ -0,0 +1,9 @@
2.25569176882662 73.37931034482759
2.3817863397548162 47.99999999999999
2.4518388791593697 34.39080459770115
2.507880910683012 31.816091954022987
2.5359019264448337 28.137931034482747
2.5779334500875657 25.011494252873554
2.6339754816112086 19.126436781609186
2.760070052539404 13.241379310344826
3.5446584938704033 6.068965517241367

55
examples/SPIN/pppm_spin/in.dipole.pppm_dipole Normal file
View File

@ -0,0 +1,55 @@
# 3d Lennard-Jones melt
units lj
#atom_style charge
atom_style hybrid sphere dipole
processors * 1 1
lattice fcc 0.8442
#region box block 0 10 0 10 0 10
region box block 0 5 0 5 0 5
create_box 3 box
create_atoms 1 box
mass * 1.0
region long block 3 6 0 10 0 10
set region long type 2
set group all dipole/random 98934 0.75
#set type 1:2 charge 0.0
velocity all create 1.0 87287
#pair_style lj/long/coul/long long off 2.5
#pair_coeff * * 1.0 1.0 2.5
#pair_coeff * 2 1.0 1.0 5.0
pair_style lj/cut/dipole/long 3.0
pair_coeff * * 1.0 1.0
#kspace_style pppm/disp 1.0e-4
kspace_style pppm/dipole 1.0e-4
kspace_modify gewald/disp 0.1
neighbor 0.3 bin
neigh_modify every 2 delay 4 check yes
group fast type 1
group slow type 2
fix 0 all balance 20 1.0 shift x 5 1.0 &
weight group 2 fast 1.0 slow 2.0 weight time 0.66
fix 1 all nve
#dump id all atom 50 dump.melt
#dump 2 all image 25 image.*.jpg type type &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3
#dump 3 all movie 25 movie.mpg type type &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3
#thermo 50
thermo 1
#run 500
run 5

61
examples/SPIN/pppm_spin/in.spin.pppm_spin Normal file
View File

@ -0,0 +1,61 @@
# hcp cobalt in a 3d periodic box
clear
units metal
atom_style spin
dimension 3
boundary p p p
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice hcp 2.5071
region box block 0.0 20.0 0.0 20.0 0.0 8.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for hcp cobalt
mass 1 58.93
#set group all spin/random 31 1.72
set group all spin 1.72 0.0 0.0 1.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/long 8.0
#pair_style hybrid/overlay eam/alloy spin/exchange 4.0
pair_coeff * * eam/alloy ../examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy Co
pair_coeff * * spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567
pair_coeff * * spin/long long 8.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
kspace_style pppm/dipole/spin 1.0e-4
#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
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
timestep 0.0001
compute out_mag all compute/spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo_style custom step time v_magnorm v_emag temp etotal
thermo 10
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_cobalt_hcp.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
run 20000

6
src/KSPACE/pppm.cpp
View File

@ -100,6 +100,9 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg),
nyhi_in = nylo_in = nyhi_out = nylo_out = 0;
nzhi_in = nzlo_in = nzhi_out = nzlo_out = 0;
// test
nlower = nupper = 0;
density_brick = vdx_brick = vdy_brick = vdz_brick = NULL;
density_fft = NULL;
u_brick = NULL;
@ -1428,12 +1431,13 @@ void PPPM::set_grid_local()
double zprd = prd[2];
double zprd_slab = zprd*slab_volfactor;
double dist[3];
double dist[3] = {0.0,0.0,0.0};
double cuthalf = 0.5*neighbor->skin + qdist;
if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf;
else kspacebbox(cuthalf,&dist[0]);
int nlo,nhi;
nlo = nhi = 0;
nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) *
nx_pppm/xprd + shift) - OFFSET;

93
src/KSPACE/pppm_dipole.cpp
View File

@ -116,9 +116,10 @@ void PPPMDipole::init()
// error check
dipoleflag = atom->mu?1:0;
qsum_qsq(0);
qsum_qsq(0); // q[i] might not be declared ?
if (dipoleflag && q2)
error->all(FLERR,"Cannot (yet) uses charges with Kspace style PPPMDipole");
error->all(FLERR,"Cannot (yet) use charges with Kspace style PPPMDipole");
triclinic_check();
@ -168,7 +169,9 @@ void PPPMDipole::init()
cutoff = *p_cutoff;
// kspace TIP4P not yet supported
// qdist = offset only for TIP4P fictitious charge
qdist = 0.0;
if (tip4pflag)
error->all(FLERR,"Cannot yet use TIP4P with PPPMDipole");
@ -206,7 +209,7 @@ void PPPMDipole::init()
"beyond nearest neighbor processor");
compute_gf_denom();
set_grid_global(mu2);
set_grid_global();
set_grid_local();
if (overlap_allowed) break;
@ -235,7 +238,7 @@ void PPPMDipole::init()
// calculate the final accuracy
double estimated_accuracy = final_accuracy_dipole(mu2);
double estimated_accuracy = final_accuracy_dipole();
// print stats
@ -607,7 +610,7 @@ void PPPMDipole::allocate()
// summation coeffs
order_allocated = order;
memory->create(gf_b,order,"pppm_dipole:gf_b");
if (!gf_b) memory->create(gf_b,order,"pppm_dipole:gf_b");
memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm_dipole:rho1d");
memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm_dipole:drho1d");
memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm_dipole:rho_coeff");
@ -791,7 +794,8 @@ void PPPMDipole::deallocate_peratom()
used for charge accumulation, FFTs, and electric field interpolation
------------------------------------------------------------------------- */
void PPPMDipole::set_grid_global(double dipole2)
//void PPPMDipole::set_grid_global(double dipole2)
void PPPMDipole::set_grid_global()
{
// use xprd,yprd,zprd
// adjust z dimension for 2d slab PPPMDipole
@ -813,14 +817,11 @@ void PPPMDipole::set_grid_global(double dipole2)
if (!gewaldflag) {
if (accuracy <= 0.0)
error->all(FLERR,"KSpace accuracy must be > 0");
//if (mu2 == 0.0)
if (dipole2 == 0.0)
if (mu2 == 0.0)
error->all(FLERR,"Must use kspace_modify gewald for systems with no dipoles");
g_ewald = (1.35 - 0.15*log(accuracy))/cutoff;
//Try Newton Solver
double g_ewald_new =
find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,dipole2);
//find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
if (g_ewald_new > 0.0) g_ewald = g_ewald_new;
else error->warning(FLERR,"PPPMDipole dipole Newton solver failed, "
"using old method to estimate g_ewald");
@ -837,6 +838,7 @@ void PPPMDipole::set_grid_global(double dipole2)
while (1) {
// set grid dimension
nx_pppm = static_cast<int> (xprd/h_x);
ny_pppm = static_cast<int> (yprd/h_y);
nz_pppm = static_cast<int> (zprd_slab/h_z);
@ -845,7 +847,8 @@ void PPPMDipole::set_grid_global(double dipole2)
if (ny_pppm <= 1) ny_pppm = 2;
if (nz_pppm <= 1) nz_pppm = 2;
//set local grid dimension
// set local grid dimension
int npey_fft,npez_fft;
if (nz_pppm >= nprocs) {
npey_fft = 1;
@ -862,7 +865,7 @@ void PPPMDipole::set_grid_global(double dipole2)
nzlo_fft = me_z*nz_pppm/npez_fft;
nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1;
double df_kspace = compute_df_kspace_dipole(dipole2);
double df_kspace = compute_df_kspace_dipole();
count++;
@ -895,7 +898,7 @@ void PPPMDipole::set_grid_global(double dipole2)
compute estimated kspace force error for dipoles
------------------------------------------------------------------------- */
double PPPMDipole::compute_df_kspace_dipole(double dipole2)
double PPPMDipole::compute_df_kspace_dipole()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
@ -903,8 +906,7 @@ double PPPMDipole::compute_df_kspace_dipole(double dipole2)
double zprd_slab = zprd*slab_volfactor;
bigint natoms = atom->natoms;
double qopt = compute_qopt_dipole();
//double df_kspace = sqrt(qopt/natoms)*mu2/(3.0*xprd*yprd*zprd_slab);
double df_kspace = sqrt(qopt/natoms)*dipole2/(3.0*xprd*yprd*zprd_slab);
double df_kspace = sqrt(qopt/natoms)*mu2/(3.0*xprd*yprd*zprd_slab);
return df_kspace;
}
@ -1101,27 +1103,27 @@ void PPPMDipole::compute_gf_dipole()
calculate f(x) for use in Newton-Raphson solver
------------------------------------------------------------------------- */
//double PPPMDipole::newton_raphson_f()
//{
// double xprd = domain->xprd;
// double yprd = domain->yprd;
// double zprd = domain->zprd;
// bigint natoms = atom->natoms;
//
// double df_rspace,df_kspace;
// double vol = xprd*yprd*zprd;
// double a = cutoff*g_ewald;
// double rg2 = a*a;
// double rg4 = rg2*rg2;
// double rg6 = rg4*rg2;
// double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
// double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
// df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
// sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2));
// df_kspace = compute_df_kspace_dipole();
//
// return df_rspace - df_kspace;
//}
double PPPMDipole::newton_raphson_f()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
bigint natoms = atom->natoms;
double df_rspace,df_kspace;
double vol = xprd*yprd*zprd;
double a = cutoff*g_ewald;
double rg2 = a*a;
double rg4 = rg2*rg2;
double rg6 = rg4*rg2;
double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2));
df_kspace = compute_df_kspace_dipole();
return df_rspace - df_kspace;
}
/* ----------------------------------------------------------------------
find g_ewald parameter for dipoles based on desired accuracy
@ -1184,7 +1186,7 @@ double PPPMDipole::derivf_dipole(double x, double Rc,
calculate the final estimate of the accuracy
------------------------------------------------------------------------- */
double PPPMDipole::final_accuracy_dipole(double dipole2)
double PPPMDipole::final_accuracy_dipole()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
@ -1193,7 +1195,7 @@ double PPPMDipole::final_accuracy_dipole(double dipole2)
bigint natoms = atom->natoms;
if (natoms == 0) natoms = 1; // avoid division by zero
double df_kspace = compute_df_kspace_dipole(mu2);
double df_kspace = compute_df_kspace_dipole();
double a = cutoff*g_ewald;
double rg2 = a*a;
@ -1201,10 +1203,7 @@ double PPPMDipole::final_accuracy_dipole(double dipole2)
double rg6 = rg4*rg2;
double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
//double df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
// sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) *
// exp(-rg2));
double df_rspace = (dipole2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
double df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) *
exp(-rg2));
@ -2521,11 +2520,11 @@ double PPPMDipole::memory_usage()
double bytes = nmax*3 * sizeof(double);
int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
(nzhi_out-nzlo_out+1);
bytes += 6 * nfft_both * sizeof(double); // vg
bytes += nfft_both * sizeof(double); // greensfn
bytes += 6 * nfft_both * sizeof(double); // vg
bytes += nfft_both * sizeof(double); // greensfn
bytes += nfft_both*5 * sizeof(FFT_SCALAR); // work*2*2
bytes += 9 * nbrick * sizeof(FFT_SCALAR); // ubrick*3 + vdbrick*6
bytes += nfft_both*7 * sizeof(FFT_SCALAR); //density_ffx*3 + work*2*2
bytes += 9 * nbrick * sizeof(FFT_SCALAR); // ubrick*3 + vdbrick*6
bytes += nfft_both*7 * sizeof(FFT_SCALAR); // density_ffx*3 + work*2*2
if (peratom_allocate_flag)
bytes += 21 * nbrick * sizeof(FFT_SCALAR);

8
src/KSPACE/pppm_dipole.h
View File

@ -37,8 +37,8 @@ class PPPMDipole : public PPPM {
double memory_usage();
protected:
void set_grid_global(double);
//double newton_raphson_f();
void set_grid_global();
double newton_raphson_f();
void allocate();
void allocate_peratom();
@ -76,7 +76,7 @@ class PPPMDipole : public PPPM {
double find_gewald_dipole(double, double, bigint, double, double);
double newton_raphson_f_dipole(double, double, bigint, double, double);
double derivf_dipole(double, double, bigint, double, double);
double compute_df_kspace_dipole(double);
double compute_df_kspace_dipole();
double compute_qopt_dipole();
void compute_gf_dipole();
void make_rho_dipole();
@ -85,7 +85,7 @@ class PPPMDipole : public PPPM {
void poisson_peratom_dipole();
void fieldforce_ik_dipole();
void fieldforce_peratom_dipole();
double final_accuracy_dipole(double dipole2);
double final_accuracy_dipole();
void musum_musq();
};

84
src/KSPACE/pppm_dipole_spin.cpp
View File

@ -50,8 +50,8 @@ using namespace MathSpecial;
#define SMALL 0.00001
#define EPS_HOC 1.0e-7
enum{REVERSE_SP};
enum{FORWARD_SP,FORWARD_SP_PERATOM};
enum{REVERSE_MU};
enum{FORWARD_MU,FORWARD_MU_PERATOM};
#ifdef FFT_SINGLE
#define ZEROF 0.0f
@ -68,10 +68,12 @@ PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp, int narg, char **arg) :
{
dipoleflag = 0;
spinflag = 1;
group_group_enable = 0;
cg_dipole = NULL;
cg_peratom_dipole = NULL;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
mub = 5.78901e-5; // in eV/T
mu_0 = 1.2566370614e-6; // in T.m/A
mub2mu0 = mub * mub * mu_0; // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
}
/* ----------------------------------------------------------------------
@ -104,7 +106,10 @@ void PPPMDipoleSpin::init()
// error check
spinflag = atom->sp?1:0;
//qsum_qsq(0); // q[i] is probably not declared ?
//if (spinflag && q2)
// error->all(FLERR,"Cannot use charges with Kspace style PPPMDipoleSpin");
triclinic_check();
if (triclinic != domain->triclinic)
@ -118,8 +123,8 @@ void PPPMDipoleSpin::init()
if (!atom->sp) error->all(FLERR,"Kspace style requires atom attribute sp");
if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff"
" ad with spins");
if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use"
" kspace_modify diff ad with spins");
if (spinflag && strcmp(update->unit_style,"metal") != 0)
error->all(FLERR,"'metal' units have to be used with spins");
@ -148,21 +153,19 @@ void PPPMDipoleSpin::init()
int itmp = 0;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
// probably not the correct extract here
if (p_cutoff == NULL)
error->all(FLERR,"KSpace style is incompatible with Pair style");
cutoff = *p_cutoff;
// kspace TIP4P not yet supported
// qdist = offset only for TIP4P fictitious charge
qdist = 0.0;
if (tip4pflag)
error->all(FLERR,"Cannot yet use TIP4P with PPPMDipoleSpin");
scale = 1.0;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
mub = 5.78901e-5; // in eV/T
mu_0 = 1.2566370614e-6; // in T.m/A
mub2mu0 = mub * mub * mu_0; // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
spsum_spsq();
natoms_original = atom->natoms;
@ -188,14 +191,14 @@ void PPPMDipoleSpin::init()
GridComm *cgtmp = NULL;
int iteration = 0;
while (order >= minorder) {
if (iteration && me == 0)
error->warning(FLERR,"Reducing PPPMDipoleSpin order b/c stencil extends "
"beyond nearest neighbor processor");
compute_gf_denom();
set_grid_global(sp2);
set_grid_global();
set_grid_local();
if (overlap_allowed) break;
@ -224,7 +227,7 @@ void PPPMDipoleSpin::init()
// calculate the final accuracy
double estimated_accuracy = final_accuracy_dipole(sp2);
double estimated_accuracy = final_accuracy_dipole();
// print stats
@ -310,7 +313,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
// return if there are no spins
if (spsqsum == 0.0) return;
if (musqsum == 0.0) return;
// convert atoms from box to lamda coords
@ -334,7 +337,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
// to fully sum contribution in their 3d bricks
// remap from 3d decomposition to FFT decomposition
cg_dipole->reverse_comm(this,REVERSE_SP);
cg_dipole->reverse_comm(this,REVERSE_MU);
brick2fft_dipole();
// compute potential gradient on my FFT grid and
@ -347,12 +350,12 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
// all procs communicate E-field values
// to fill ghost cells surrounding their 3d bricks
cg_dipole->forward_comm(this,FORWARD_SP);
cg_dipole->forward_comm(this,FORWARD_MU);
// extra per-atom energy/virial communication
if (evflag_atom) {
cg_peratom_dipole->forward_comm(this,FORWARD_SP_PERATOM);
cg_peratom_dipole->forward_comm(this,FORWARD_MU_PERATOM);
}
// calculate the force on my particles
@ -374,7 +377,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
energy = energy_all;
energy *= 0.5*volume;
energy -= spsqsum*2.0*g3/3.0/MY_PIS;
energy -= musqsum*2.0*g3/3.0/MY_PIS;
energy *= spscale;
}
@ -510,7 +513,7 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
double spx,spy,spz;
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **fm_long = atom->fm_long;
int nlocal = atom->nlocal;
@ -548,7 +551,7 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
}
}
// convert M-field to mech. and mag. forces
// convert M-field and store mech. forces
const double spfactor = mub2mu0 * scale;
spx = sp[i][0]*sp[i][3];
@ -558,13 +561,12 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
f[i][1] += spfactor*(vxy*spx + vyy*spy + vyz*spz);
f[i][2] += spfactor*(vxz*spx + vyz*spy + vzz*spz);
// store long-range mag. precessions
const double spfactorh = mub2mu0hbinv * scale;
fm[i][0] += spfactorh*ex;
fm[i][1] += spfactorh*ey;
fm[i][2] += spfactorh*ez;
// create a new vector (in atom_spin style ?) to store long-range fm tables
fm_long[i][0] += spfactorh*ex;
fm_long[i][1] += spfactorh*ey;
fm_long[i][2] += spfactorh*ez;
}
}
@ -708,9 +710,11 @@ void PPPMDipoleSpin::slabcorr()
// add on mag. force corrections
double ffact = spscale * (-4.0*MY_PI/volume);
double **fm = atom->fm;
//double **fm = atom->fm;
double **fm_long = atom->fm_long;
for (int i = 0; i < nlocal; i++) {
fm[i][2] += ffact * spin_all;
//fm[i][2] += ffact * spin_all;
fm_long[i][2] += ffact * spin_all;
}
}
@ -723,13 +727,13 @@ void PPPMDipoleSpin::spsum_spsq()
{
const int nlocal = atom->nlocal;
spsum = spsqsum = sp2 = 0.0;
musum = musqsum = mu2 = 0.0;
if (atom->sp_flag) {
double **sp = atom->sp;
double spx, spy, spz;
double spsum_local(0.0), spsqsum_local(0.0);
// not exactly the good loop: need to add norm of spins
// sum (direction x norm) of all spins
for (int i = 0; i < nlocal; i++) {
spx = sp[i][0]*sp[i][3];
@ -739,12 +743,14 @@ void PPPMDipoleSpin::spsum_spsq()
spsqsum_local += spx*spx + spy*spy + spz*spz;
}
MPI_Allreduce(&spsum_local,&spsum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&spsqsum_local,&spsqsum,1,MPI_DOUBLE,MPI_SUM,world);
// store results into pppm_dipole quantities
sp2 = spsqsum * mub2mu0;
MPI_Allreduce(&spsum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&spsqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world);
mu2 = musqsum * mub2mu0;
}
if (sp2 == 0 && comm->me == 0)
if (mu2 == 0 && comm->me == 0)
error->all(FLERR,"Using kspace solver PPPMDipoleSpin on system with no spins");
}

1
src/KSPACE/pppm_dipole_spin.h
View File

@ -42,7 +42,6 @@ class PPPMDipoleSpin : public PPPMDipole {
// spin
double spsum,spsqsum,sp2;
void make_rho_spin();
void fieldforce_ik_spin();
void fieldforce_peratom_spin();

1
src/SPIN/atom_vec_spin.cpp
View File

@ -106,7 +106,6 @@ void AtomVecSpin::grow_reset()
sp = atom->sp; fm = atom->fm; fm_long = atom->fm_long;
}
/* ----------------------------------------------------------------------
copy atom I info to atom J
------------------------------------------------------------------------- */

18
src/SPIN/fix_nve_spin.cpp
View File

@ -126,6 +126,7 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
// initialize the magnetic interaction flags
pair_spin_flag = 0;
long_spin_flag = 0;
precession_spin_flag = 0;
maglangevin_flag = 0;
tdamp_flag = temp_flag = 0;
@ -209,8 +210,16 @@ void FixNVESpin::init()
if (count != npairspin)
error->all(FLERR,"Incorrect number of spin pairs");
// set pair/spin and long/spin flags
if (npairspin >= 1) pair_spin_flag = 1;
for (int i = 0; i<npairs; i++) {
if (force->pair_match("spin/long",0,i)) {
long_spin_flag = 1;
}
}
// ptrs FixPrecessionSpin classes
int iforce;
@ -425,6 +434,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
double **sp = atom->sp;
double **fm = atom->fm;
double **fm_long = atom->fm_long;
// force computation for spin i
@ -442,6 +452,14 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
}
}
// update magnetic long-range components
if (long_spin_flag) {
fmi[0] += fm_long[i][0];
fmi[1] += fm_long[i][1];
fmi[2] += fm_long[i][2];
}
// update magnetic precession interactions
if (precession_spin_flag) {

3
src/SPIN/fix_nve_spin.h
View File

@ -58,7 +58,8 @@ friend class PairSpin;
int nlocal_max; // max value of nlocal (for lists size)
int pair_spin_flag; // magnetic pair flags
int precession_spin_flag; // magnetic precession flags
int long_spin_flag; // magnetic long-range flag
int precession_spin_flag; // magnetic precession flags
int maglangevin_flag; // magnetic langevin flags
int tdamp_flag, temp_flag;

11
src/SPIN/pair_spin_long.cpp
View File

@ -81,6 +81,7 @@ PairSpinLong::~PairSpinLong()
if (allocated) {
memory->destroy(setflag);
memory->destroy(cut_spin_long);
memory->destroy(cutsq);
}
}
@ -97,7 +98,6 @@ void PairSpinLong::settings(int narg, char **arg)
error->all(FLERR,"Spin simulations require metal unit style");
cut_spin_long_global = force->numeric(FLERR,arg[0]);
//cut_spin = force->numeric(FLERR,arg[0]);
// reset cutoffs that have been explicitly set
@ -126,7 +126,7 @@ void PairSpinLong::coeff(int narg, char **arg)
if (strcmp(arg[2],"long") != 0)
error->all(FLERR,"Incorrect args in pair_style command");
if (narg != 3)
if (narg < 1 || narg > 4)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
@ -197,9 +197,9 @@ void PairSpinLong::init_style()
double PairSpinLong::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cut_spin_long[j][i] = cut_spin_long[i][j];
return cut_spin_long_global;
}
@ -505,7 +505,8 @@ void PairSpinLong::allocate()
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cut_spin_long,n+1,n+1,"pair:cut_spin_long");
memory->create(cut_spin_long,n+1,n+1,"pair/spin/long:cut_spin_long");
memory->create(cutsq,n+1,n+1,"pair/spin/long:cutsq");
}
/* ----------------------------------------------------------------------

3
src/atom.cpp
View File

@ -98,7 +98,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// SPIN package
sp = fm = NULL;
sp = fm = fm_long = NULL;
// USER-DPD
@ -277,6 +277,7 @@ Atom::~Atom()
memory->destroy(sp);
memory->destroy(fm);
memory->destroy(fm_long);
memory->destroy(vfrac);
memory->destroy(s0);

1
src/atom.h
View File

@ -65,6 +65,7 @@ class Atom : protected Pointers {
double **sp;
double **fm;
double **fm_long;
// PERI package

2
src/pair.h
View File

@ -61,7 +61,7 @@ class Pair : protected Pointers {
int dispersionflag; // 1 if compatible with LJ/dispersion solver
int tip4pflag; // 1 if compatible with TIP4P solver
int dipoleflag; // 1 if compatible with dipole solver
int spinflag; // 1 if compatible with spin long solver
int spinflag; // 1 if compatible with spin solver
int reinitflag; // 1 if compatible with fix adapt and alike
int tail_flag; // pair_modify flag for LJ tail correction