From 82a5346ab1f1ddfa8ad114a7dceff47d0cdbe0ee Mon Sep 17 00:00:00 2001 From: julient31 Date: Fri, 14 Sep 2018 15:09:59 -0600 Subject: [PATCH] Commit JT 091418 - created pair_spin_long_qsymp - modified ewald_dipole --- examples/SPIN/pppm_spin/in.spin.pppm_spin | 7 +- src/KSPACE/ewald_dipole.cpp | 100 ++-- src/KSPACE/ewald_dipole.h | 6 +- src/SPIN/fix_nve_spin.cpp | 17 +- src/SPIN/fix_nve_spin.h | 1 - src/SPIN/pair_spin_long.cpp | 24 +- src/SPIN/pair_spin_long_qsymp.cpp | 655 ++++++++++++++++++++++ src/SPIN/pair_spin_long_qsymp.h | 100 ++++ 8 files changed, 828 insertions(+), 82 deletions(-) create mode 100644 src/SPIN/pair_spin_long_qsymp.cpp create mode 100644 src/SPIN/pair_spin_long_qsymp.h diff --git a/examples/SPIN/pppm_spin/in.spin.pppm_spin b/examples/SPIN/pppm_spin/in.spin.pppm_spin index f7e462c343..9e57797f55 100644 --- a/examples/SPIN/pppm_spin/in.spin.pppm_spin +++ b/examples/SPIN/pppm_spin/in.spin.pppm_spin @@ -23,17 +23,20 @@ mass 1 58.93 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 spin/long 8.0 +pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/long/qsymp 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 +pair_coeff * * spin/long/qsymp long 8.0 +#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 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 fix 2 all langevin/spin 0.0 0.0 21 diff --git a/src/KSPACE/ewald_dipole.cpp b/src/KSPACE/ewald_dipole.cpp index d03d93a7c5..3b3e3b93db 100644 --- a/src/KSPACE/ewald_dipole.cpp +++ b/src/KSPACE/ewald_dipole.cpp @@ -79,9 +79,9 @@ void EwaldDipole::init() triclinic_check(); - // set triclinic to 0 for now (no triclinic calc.) - triclinic = 0; - + // no triclinic ewald dipole (yet) + + triclinic = domain->triclinic; if (triclinic) error->all(FLERR,"Cannot (yet) use EwaldDipole with triclinic box"); @@ -100,9 +100,6 @@ void EwaldDipole::init() if (domain->xperiodic != 1 || domain->yperiodic != 1 || domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) error->all(FLERR,"Incorrect boundaries with slab EwaldDipole"); - //if (domain->triclinic) - // error->all(FLERR,"Cannot (yet) use EwaldDipole with triclinic box " - // "and slab correction"); } // extract short-range Coulombic cutoff from pair style @@ -278,23 +275,6 @@ void EwaldDipole::setup() kymax_orig = kymax; kzmax_orig = kzmax; - // scale lattice vectors for triclinic skew - - //if (triclinic) { - // double tmp[3]; - // tmp[0] = kxmax/xprd; - // tmp[1] = kymax/yprd; - // tmp[2] = kzmax/zprd; - // lamda2xT(&tmp[0],&tmp[0]); - // kxmax = MAX(1,static_cast(tmp[0])); - // kymax = MAX(1,static_cast(tmp[1])); - // kzmax = MAX(1,static_cast(tmp[2])); - - // kmax = MAX(kxmax,kymax); - // kmax = MAX(kmax,kzmax); - // kmax3d = 4*kmax*kmax*kmax + 6*kmax*kmax + 3*kmax; - //} - } else { kxmax = kx_ewald; @@ -340,10 +320,6 @@ void EwaldDipole::setup() // pre-compute EwaldDipole coefficients coeffs(); - //if (triclinic == 0) - // coeffs(); - //else - // coeffs_triclinic(); } /* ---------------------------------------------------------------------- @@ -380,7 +356,6 @@ void EwaldDipole::compute(int eflag, int vflag) // if atom count has changed, update qsum and qsqsum if (atom->natoms != natoms_original) { - //qsum_qsq(); musum_musq(); natoms_original = atom->natoms; } @@ -407,10 +382,6 @@ void EwaldDipole::compute(int eflag, int vflag) // partial structure factors on each processor // total structure factor by summing over procs - //if (triclinic == 0) - // eik_dot_r(); - //else - // eik_dot_r_triclinic(); eik_dot_r(); MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world); @@ -421,7 +392,8 @@ void EwaldDipole::compute(int eflag, int vflag) // perform per-atom calculations if needed double **f = atom->f; - double *q = atom->q; + //double *q = atom->q; + double **mu = atom->mu; int nlocal = atom->nlocal; int kx,ky,kz; @@ -439,21 +411,34 @@ void EwaldDipole::compute(int eflag, int vflag) kz = kzvecs[k]; for (i = 0; i < nlocal; i++) { + + // calculating exp(i*k*ri) + cypz = cs[ky][1][i]*cs[kz][2][i] - sn[ky][1][i]*sn[kz][2][i]; sypz = sn[ky][1][i]*cs[kz][2][i] + cs[ky][1][i]*sn[kz][2][i]; exprl = cs[kx][0][i]*cypz - sn[kx][0][i]*sypz; expim = sn[kx][0][i]*cypz + cs[kx][0][i]*sypz; + + // taking im-part of struct_fact x exp(i*k*ri) (for force calc.) + partial = expim*sfacrl_all[k] - exprl*sfacim_all[k]; - ek[i][0] += partial*eg[k][0]; - ek[i][1] += partial*eg[k][1]; - ek[i][2] += partial*eg[k][2]; + //ek[i][0] += partial*eg[k][0]; + //ek[i][1] += partial*eg[k][1]; + //ek[i][2] += partial*eg[k][2]; + ek[i][0] += kx*partial*eg[k][0]; + ek[i][1] += ky*partial*eg[k][1]; + ek[i][2] += kz*partial*eg[k][2]; if (evflag_atom) { + + // taking re-part of struct_fact x exp(i*k*ri) (for energy calc.) + partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k]; //if (eflag_atom) eatom[i] += q[i]*ug[k]*partial_peratom; if (eflag_atom) eatom[i] += muk[k][i]*ug[k]*partial_peratom; if (vflag_atom) for (j = 0; j < 6; j++) + // to be done vatom[i][j] += ug[k]*vg[k][j]*partial_peratom; } } @@ -461,30 +446,33 @@ void EwaldDipole::compute(int eflag, int vflag) // convert E-field to force - const double qscale = qqrd2e * scale; + //const double qscale = qqrd2e * scale; const double muscale = qqrd2e * scale; for (i = 0; i < nlocal; i++) { - for (k = 0; k < kcount; k++) { - //f[i][0] += qscale * q[i]*ek[i][0]; - //f[i][1] += qscale * q[i]*ek[i][1]; - //if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2]; - f[i][0] += muscale * muk[k][i] * ek[i][0]; - f[i][1] += muscale * muk[k][i] * ek[i][1]; - if (slabflag != 2) f[i][2] += muscale * muk[k][i] * ek[i][2]; - } + //f[i][0] += qscale * q[i]*ek[i][0]; + //f[i][1] += qscale * q[i]*ek[i][1]; + //if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2]; + f[i][0] += muscale * mu[i][0] * ek[i][0]; + f[i][1] += muscale * mu[i][1] * ek[i][1]; + if (slabflag != 2) f[i][2] += muscale * mu[i][2] * ek[i][2]; } // sum global energy across Kspace vevs and add in volume-dependent term if (eflag_global) { - for (k = 0; k < kcount; k++) + for (k = 0; k < kcount; k++) { + + // taking the re-part of struct_fact_i x struct_fact_j + energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); + } - energy -= g_ewald*qsqsum/MY_PIS + - MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); - energy *= qscale; + // substracting self energy and scaling + + energy -= musqsum*2.0*g3/3.0/MY_PIS; + energy *= muscale; } // global virial @@ -495,7 +483,8 @@ void EwaldDipole::compute(int eflag, int vflag) uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); for (j = 0; j < 6; j++) virial[j] += uk*vg[k][j]; } - for (j = 0; j < 6; j++) virial[j] *= qscale; + //for (j = 0; j < 6; j++) virial[j] *= qscale; + for (j = 0; j < 6; j++) virial[j] *= muscale; } // per-atom energy/virial @@ -504,9 +493,9 @@ void EwaldDipole::compute(int eflag, int vflag) if (evflag_atom) { if (eflag_atom) { for (i = 0; i < nlocal; i++) { - eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum / - (g_ewald*g_ewald*volume); - eatom[i] *= qscale; + eatom[i] -= (mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]) + *2.0*g3/3.0/MY_PIS; + eatom[i] *= muscale; } } @@ -520,7 +509,9 @@ void EwaldDipole::compute(int eflag, int vflag) if (slabflag == 1) slabcorr(); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + compute the +------------------------------------------------------------------------- */ void EwaldDipole::eik_dot_r() { @@ -530,7 +521,6 @@ void EwaldDipole::eik_dot_r() double mux, muy, muz; double **x = atom->x; - //double *q = atom->q; double **mu = atom->mu; int nlocal = atom->nlocal; diff --git a/src/KSPACE/ewald_dipole.h b/src/KSPACE/ewald_dipole.h index 5cd969bbee..401742ed3a 100644 --- a/src/KSPACE/ewald_dipole.h +++ b/src/KSPACE/ewald_dipole.h @@ -34,17 +34,13 @@ class EwaldDipole : public Ewald { protected: double musum,musqsum,mu2; - double **muk; // mu_i dot k + double **muk; // store mu_i dot k void musum_musq(); double rms_dipole(int, double, bigint); virtual void eik_dot_r(); void slabcorr(); - // triclinic - - //void eik_dot_r_triclinic(); - }; } diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index 996bd3c2da..5e972cd14d 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -295,6 +295,13 @@ void FixNVESpin::initial_integrate(int vflag) } } + // update fm_kspace if long-range + // remove short-range comp. of fm_kspace + + if (long_spin_flag) { + + } + // update half s for all atoms if (sector_flag) { // sectoring seq. update @@ -434,7 +441,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i) double **sp = atom->sp; double **fm = atom->fm; - double **fm_long = atom->fm_long; + //double **fm_long = atom->fm_long; // force computation for spin i @@ -452,14 +459,6 @@ 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) { diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 565de13e92..9fcbfb3803 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -48,7 +48,6 @@ friend class PairSpin; int lattice_flag; // lattice_flag = 0 if spins only // lattice_flag = 1 if spin-lattice calc. - protected: int sector_flag; // sector_flag = 0 if serial algorithm // sector_flag = 1 if parallel algorithm diff --git a/src/SPIN/pair_spin_long.cpp b/src/SPIN/pair_spin_long.cpp index efedea3247..d7ecdf5edd 100644 --- a/src/SPIN/pair_spin_long.cpp +++ b/src/SPIN/pair_spin_long.cpp @@ -362,13 +362,15 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3]) int i,j,jj,jnum,itype,jtype; double r,rinv,r2inv,rsq; double grij,expm2,t,erfc; - double bij[4],xi[3],rij[3],spi[4],spj[4]; + double bij[4],xi[3],rij[3]; + double spi[4],spj[4]; double local_cut2; double pre1,pre2,pre3; int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double **sp = atom->sp; + double **fm_long = atom->fm_long; int *type = atom->type; ilist = list->ilist; @@ -433,15 +435,16 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3]) } } - // force accumulation + // adding the kspace components to fm + + fmi[0] += fm_long[i][0]; + fmi[1] += fm_long[i][1]; + fmi[2] += fm_long[i][2]; - fmi[0] *= mub2mu0hbinv; - fmi[1] *= mub2mu0hbinv; - fmi[2] *= mub2mu0hbinv; } /* ---------------------------------------------------------------------- - compute exchange interaction between spins i and j + compute dipolar interaction between spins i and j ------------------------------------------------------------------------- */ void PairSpinLong::compute_long(int i, int j, double rij[3], @@ -456,13 +459,14 @@ void PairSpinLong::compute_long(int i, int j, double rij[3], b1 = bij[1]; b2 = bij[2]; - fmi[0] += gigj * (b2 * sjdotr *rij[0] - b1 * spj[0]); - fmi[1] += gigj * (b2 * sjdotr *rij[1] - b1 * spj[1]); - fmi[2] += gigj * (b2 * sjdotr *rij[2] - b1 * spj[2]); + fmi[0] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[0] - b1 * spj[0]); + fmi[1] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[1] - b1 * spj[1]); + fmi[2] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[2] - b1 * spj[2]); } /* ---------------------------------------------------------------------- - compute the mechanical force due to the exchange interaction between atom i and atom j + compute the mechanical force due to the dipolar interaction between + atom i and atom j ------------------------------------------------------------------------- */ void PairSpinLong::compute_long_mech(int i, int j, double rij[3], diff --git a/src/SPIN/pair_spin_long_qsymp.cpp b/src/SPIN/pair_spin_long_qsymp.cpp new file mode 100644 index 0000000000..3b499d0ef7 --- /dev/null +++ b/src/SPIN/pair_spin_long_qsymp.cpp @@ -0,0 +1,655 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ------------------------------------------------------------------------ + Contributing authors: Julien Tranchida (SNL) + Aidan Thompson (SNL) + + Please cite the related publication: + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics + and molecular dynamics. Journal of Computational Physics. +------------------------------------------------------------------------- */ + +#include +#include +#include +#include + +#include "pair_spin_long_qsymp.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "fix_nve_spin.h" +#include "force.h" +#include "kspace.h" +#include "math_const.h" +#include "memory.h" +#include "modify.h" +#include "error.h" +#include "update.h" + + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairSpinLongQsymp::PairSpinLongQsymp(LAMMPS *lmp) : PairSpin(lmp), +lockfixnvespin(NULL) +{ + single_enable = 0; + ewaldflag = pppmflag = spinflag = 1; + respa_enable = 0; + no_virial_fdotr_compute = 1; + lattice_flag = 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 + +} + +/* ---------------------------------------------------------------------- + free all arrays +------------------------------------------------------------------------- */ + +PairSpinLongQsymp::~PairSpinLongQsymp() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cut_spin_long); + memory->destroy(cutsq); + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::settings(int narg, char **arg) +{ + if (narg < 1 || narg > 2) + error->all(FLERR,"Incorrect args in pair_style command"); + + if (strcmp(update->unit_style,"metal") != 0) + error->all(FLERR,"Spin simulations require metal unit style"); + + cut_spin_long_global = force->numeric(FLERR,arg[0]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) { + for (j = i+1; j <= atom->ntypes; j++) { + if (setflag[i][j]) { + cut_spin_long[i][j] = cut_spin_long_global; + } + } + } + } + +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + + // check if args correct + + if (strcmp(arg[2],"long") != 0) + error->all(FLERR,"Incorrect args in pair_style command"); + if (narg < 1 || narg > 4) + error->all(FLERR,"Incorrect args in pair_style command"); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double spin_long_cut_one = force->numeric(FLERR,arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + setflag[i][j] = 1; + cut_spin_long[i][j] = spin_long_cut_one; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::init_style() +{ + if (!atom->sp_flag) + error->all(FLERR,"Pair spin requires atom/spin style"); + + // need a full neighbor list + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + + // checking if nve/spin is a listed fix + + int ifix = 0; + while (ifix < modify->nfix) { + if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; + ifix++; + } + if (ifix == modify->nfix) + error->all(FLERR,"pair/spin style requires nve/spin"); + + // get the lattice_flag from nve/spin + + for (int i = 0; i < modify->nfix; i++) { + if (strcmp(modify->fix[i]->style,"nve/spin") == 0) { + lockfixnvespin = (FixNVESpin *) modify->fix[i]; + lattice_flag = lockfixnvespin->lattice_flag; + } + } + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + + g_ewald = force->kspace->g_ewald; + +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairSpinLongQsymp::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; +} + +/* ---------------------------------------------------------------------- + extract the larger cutoff if "cut" or "cut_coul" +------------------------------------------------------------------------- */ + +void *PairSpinLongQsymp::extract(const char *str, int &dim) +{ + if (strcmp(str,"cut") == 0) { + dim = 0; + return (void *) &cut_spin_long_global; + } else if (strcmp(str,"cut_coul") == 0) { + dim = 0; + return (void *) &cut_spin_long_global; + } else if (strcmp(str,"ewald_order") == 0) { + ewald_order = 0; + ewald_order |= 1<<1; + ewald_order |= 1<<3; + dim = 0; + return (void *) &ewald_order; + } else if (strcmp(str,"ewald_mix") == 0) { + dim = 0; + return (void *) &mix_flag; + } + return NULL; +} + +/* ---------------------------------------------------------------------- */ + +void PairSpinLongQsymp::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double r,rinv,r2inv,rsq; + double grij,expm2,t,erfc,erf; + double sjdotr,gigj; + double bij[4]; + double cij[4]; + double evdwl,ecoul; + double xi[3],rij[3]; + double spi[4],spj[4],fi[3],fmi[3]; + double fmx_erf_s,fmy_erf_s,fmz_erf_s; + double local_cut2; + double pre1,pre2,pre3; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double **fm = atom->fm; + double **fm_long = atom->fm_long; + double **sp = atom->sp; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + pre1 = 2.0 * g_ewald / MY_PIS; + pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS; + pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS; + + // computation of the exchange interaction + // loop over atoms and their neighbors + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xi[0] = x[i][0]; + xi[1] = x[i][1]; + xi[2] = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + spi[0] = sp[i][0]; + spi[1] = sp[i][1]; + spi[2] = sp[i][2]; + spi[3] = sp[i][3]; + itype = type[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = type[j]; + + spj[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; + spj[3] = sp[j][3]; + + evdwl = 0.0; + + fi[0] = fi[1] = fi[2] = 0.0; + fmi[0] = fmi[1] = fmi[2] = 0.0; + fmx_erf_s = fmy_erf_s = fmz_erf_s = 0.0; + bij[0] = bij[1] = bij[2] = bij[3] = 0.0; + cij[0] = cij[1] = cij[2] = cij[3] = 0.0; + + rij[0] = x[j][0] - xi[0]; + rij[1] = x[j][1] - xi[1]; + rij[2] = x[j][2] - xi[2]; + rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; + + local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype]; + + if (rsq < local_cut2) { + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + erf = 1.0 - t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + // evaluating erfc for mech. force and energy calc. + + bij[0] = erfc * rinv; + bij[1] = (bij[0] + pre1*expm2) * r2inv; + bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; + bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; + + compute_long(i,j,rij,bij,fmi,spi,spj); + compute_long_mech(i,j,rij,bij,fmi,spi,spj); + + // evaluating erf comp. for fm_kspace correction + + cij[0] = erf * rinv; + cij[1] = (bij[0] + pre1*expm2) * r2inv; + cij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; + //cij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; + + gigj = spi[3] * spj[3]; + sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2]; + + // evaluating short-range correction to the kspace part on [0,rc] + + fmx_erf_s += gigj * (cij[2] * sjdotr * rij[0] - cij[1] * spj[0]); + fmy_erf_s += gigj * (cij[2] * sjdotr * rij[1] - cij[1] * spj[1]); + fmz_erf_s += gigj * (cij[2] * sjdotr * rij[2] - cij[1] * spj[2]); + + } + + // force accumulation + + f[i][0] += fi[0] * mub2mu0; + f[i][1] += fi[1] * mub2mu0; + f[i][2] += fi[2] * mub2mu0; + fm[i][0] += fmi[0] * mub2mu0hbinv; + fm[i][1] += fmi[1] * mub2mu0hbinv; + fm[i][2] += fmi[2] * mub2mu0hbinv; + + // correction of the fm_kspace + + fm_long[i][0] -= mub2mu0hbinv * fmx_erf_s; + fm_long[i][1] -= mub2mu0hbinv * fmy_erf_s; + fm_long[i][2] -= mub2mu0hbinv * fmz_erf_s; + + if (newton_pair || j < nlocal) { + f[j][0] -= fi[0]; + f[j][1] -= fi[1]; + f[j][2] -= fi[2]; + } + + if (eflag) { + if (rsq <= local_cut2) { + evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] + + spi[2]*fmi[2]; + evdwl *= hbar; + } + } else evdwl = 0.0; + + + if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, + evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); + + } + } +} + +/* ---------------------------------------------------------------------- + update the pair interaction fmi acting on the spin ii + adding 1/r (for r in [0,rc]) contribution to the pair + removing erf(r)/r (for r in [0,rc]) from the kspace force +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::compute_single_pair(int ii, double fmi[3]) +{ + int i,j,jj,jnum,itype,jtype; + double r,rinv,r2inv,r3inv,rsq; + double grij,expm2,t,erf; + double sjdotr,sjdotrr3inv; + double b1,b2,gigj; + double bij[4],xi[3],rij[3]; + double spi[4],spj[4]; + double local_cut2; + double pre1,pre2,pre3; + int *ilist,*jlist,*numneigh,**firstneigh; + //double fmx_erf_s,fmy_erf_s,fmz_erf_s; + double fmx_s,fmy_s,fmz_s; + //double fmx_long,fmy_long,fmz_long; + + double **x = atom->x; + double **sp = atom->sp; + double **fm_long = atom->fm_long; + int *type = atom->type; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + pre1 = 2.0 * g_ewald / MY_PIS; + pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS; + pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS; + + // computation of the exchange interaction + // loop over neighbors of atom i + + i = ilist[ii]; + xi[0] = x[i][0]; + xi[1] = x[i][1]; + xi[2] = x[i][2]; + spi[0] = sp[i][0]; + spi[1] = sp[i][1]; + spi[2] = sp[i][2]; + spi[3] = sp[i][3]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + itype = type[i]; + + //fmx_long = fmy_long = fmz_long = 0.0; + //fmx_erf_s = fmy_erf_s = fmz_erf_s = 0.0; + fmx_s = fmy_s = fmz_s = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = type[j]; + + spj[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; + spj[3] = sp[j][3]; + + bij[0] = bij[1] = bij[2] = bij[3] = 0.0; + + rij[0] = x[j][0] - xi[0]; + rij[1] = x[j][1] - xi[1]; + rij[2] = x[j][2] - xi[2]; + rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; + + local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype]; + + if (rsq < local_cut2) { + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + r3inv = r2inv*rinv; + + r = sqrt(rsq); + //grij = g_ewald * r; + //expm2 = exp(-grij*grij); + //t = 1.0 / (1.0 + EWALD_P*grij); + + // evaluating erf instead of erfc + + //erf = 1.0 - t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + //bij[0] = erf * rinv; + //bij[1] = (bij[0] + pre1*expm2) * r2inv; + //bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; + //bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; + + //gigj = spi[3] * spj[3]; + //sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2]; + + //b1 = bij[1]; + //b2 = bij[2]; + + // evaluating short-range correction to the kspace part on [0,rc] + + //fmx_erf_s += mub2mu0hbinv * gigj * (b2 * sjdotr * rij[0] - b1 * spj[0]); + //fmy_erf_s += mub2mu0hbinv * gigj * (b2 * sjdotr * rij[1] - b1 * spj[1]); + //fmz_erf_s += mub2mu0hbinv * gigj * (b2 * sjdotr * rij[2] - b1 * spj[2]); + + // evaluating real dipolar interaction on [0,rc] + + sjdotrr3inv = 3.0 * sjdotr * r3inv; + fmx_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[0] - sp[i][0]/r3inv); + fmy_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[1] - sp[i][1]/r3inv); + fmz_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[2] - sp[i][2]/r3inv); + + } + } + + // removing short-range erf function from kspace force + + //fmx_long = fm_long[i][0] - fmx_erf_s; + //fmy_long = fm_long[i][1] - fmy_erf_s; + //fmz_long = fm_long[i][2] - fmz_erf_s; + + // adding truncated kspace force and short-range full force + + //fmi[0] += fmx_s + fmx_long; + //fmi[1] += fmy_s + fmy_long; + //fmi[2] += fmz_s + fmz_long; + fmi[0] += (fmx_s + fm_long[i][0]); + fmi[1] += (fmy_s + fm_long[i][1]); + fmi[2] += (fmz_s + fm_long[i][2]); +} + +/* ---------------------------------------------------------------------- + compute dipolar interaction between spins i and j +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::compute_long(int i, int j, double rij[3], + double bij[4], double fmi[3], double spi[4], double spj[4]) +{ + double sjdotr; + double b1,b2,gigj; + + gigj = spi[3] * spj[3]; + sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2]; + + b1 = bij[1]; + b2 = bij[2]; + + fmi[0] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[0] - b1 * spj[0]); + fmi[1] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[1] - b1 * spj[1]); + fmi[2] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[2] - b1 * spj[2]); +} + +/* ---------------------------------------------------------------------- + compute the mechanical force due to the dipolar interaction between + atom i and atom j +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::compute_long_mech(int i, int j, double rij[3], + double bij[4], double fi[3], double spi[3], double spj[3]) +{ + double sdots,sidotr,sjdotr,b2,b3; + double g1,g2,g1b2_g2b3,gigj; + + gigj = spi[3] * spj[3]; + sdots = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2]; + sidotr = spi[0]*rij[0] + spi[1]*rij[1] + spi[2]*rij[2]; + sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2]; + + b2 = bij[2]; + b3 = bij[3]; + g1 = sdots; + g2 = -sidotr*sjdotr; + g1b2_g2b3 = g1*b2 + g2*b3; + + fi[0] += gigj * (rij[0] * g1b2_g2b3 + + b2 * (sjdotr*spi[0] + sidotr*spj[0])); + fi[1] += gigj * (rij[1] * g1b2_g2b3 + + b2 * (sjdotr*spi[1] + sidotr*spj[1])); + fi[2] += gigj * (rij[2] * g1b2_g2b3 + + b2 * (sjdotr*spi[2] + sidotr*spj[2])); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + 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"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) { + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&cut_spin_long[i][j],sizeof(int),1,fp); + } + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) { + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&cut_spin_long[i][j],sizeof(int),1,fp); + } + MPI_Bcast(&cut_spin_long[i][j],1,MPI_INT,0,world); + } + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::write_restart_settings(FILE *fp) +{ + fwrite(&cut_spin_long_global,sizeof(double),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairSpinLongQsymp::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_spin_long_global,sizeof(double),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_spin_long_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} diff --git a/src/SPIN/pair_spin_long_qsymp.h b/src/SPIN/pair_spin_long_qsymp.h new file mode 100644 index 0000000000..ae8c5a3864 --- /dev/null +++ b/src/SPIN/pair_spin_long_qsymp.h @@ -0,0 +1,100 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(spin/long/qsymp,PairSpinLongQsymp) + +#else + +#ifndef LMP_PAIR_SPIN_LONG_QSYMP_H +#define LMP_PAIR_SPIN_LONG_QSYMP_H + +#include "pair_spin.h" + +namespace LAMMPS_NS { + +class PairSpinLongQsymp : public PairSpin { + public: + double cut_coul; + double **sigma; + + PairSpinLongQsymp(class LAMMPS *); + ~PairSpinLongQsymp(); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + void *extract(const char *, int &); + + void compute(int, int); + void compute_single_pair(int, double *); + + void compute_long(int, int, double *, double *, double *, + double *, double *); + void compute_long_mech(int, int, double *, double *, double *, + double *, double *); + + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + + double cut_spin_long_global; // global long cutoff distance + + protected: + double hbar; // reduced Planck's constant + double mub; // Bohr's magneton + double mu_0; // vacuum permeability + double mub2mu0; // prefactor for mech force + double mub2mu0hbinv; // prefactor for mag force + + double **cut_spin_long; // cutoff distance long + + double g_ewald; + int ewald_order; + + int lattice_flag; // flag for mech force computation + class FixNVESpin *lockfixnvespin; // ptr for setups + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args in pair_style command + +Self-explanatory. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair dipole/long requires atom attributes q, mu, torque + +The atom style defined does not have these attributes. + +E: Cannot (yet) use 'electron' units with dipoles + +This feature is not yet supported. + +E: Pair style requires a KSpace style + +No kspace style is defined. + +*/