diff --git a/examples/SPIN/pppm_spin/in.spin.spin_dipolar_cut b/examples/SPIN/pppm_spin/in.spin.spin_dipolar_cut new file mode 100644 index 0000000000..838181e6d8 --- /dev/null +++ b/examples/SPIN/pppm_spin/in.spin.spin_dipolar_cut @@ -0,0 +1,62 @@ +# 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 8.0 0.0 8.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 spin/dipolar/cut 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/dipolar/cut long 8.0 +#pair_coeff * * spin/long long 8.0 + +neighbor 0.1 bin +neigh_modify every 10 check yes delay 20 + +#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 +run 10 diff --git a/src/SPIN/pair_spin_dipolar_cut.cpp b/src/SPIN/pair_spin_dipolar_cut.cpp new file mode 100644 index 0000000000..f686d12926 --- /dev/null +++ b/src/SPIN/pair_spin_dipolar_cut.cpp @@ -0,0 +1,528 @@ +/* ---------------------------------------------------------------------- + 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_dipolar_cut.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 + +/* ---------------------------------------------------------------------- */ + +PairSpinDipolarCut::PairSpinDipolarCut(LAMMPS *lmp) : PairSpin(lmp), +lockfixnvespin(NULL) +{ + single_enable = 0; + 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 +------------------------------------------------------------------------- */ + +PairSpinDipolarCut::~PairSpinDipolarCut() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cut_spin_long); + memory->destroy(cutsq); + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairSpinDipolarCut::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 PairSpinDipolarCut::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 PairSpinDipolarCut::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; + } + } + +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairSpinDipolarCut::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 *PairSpinDipolarCut::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 PairSpinDipolarCut::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double rinv,r2inv,r3inv,rsq; + double evdwl,ecoul; + double xi[3],rij[3]; + double spi[4],spj[4],fi[3],fmi[3]; + double local_cut2; + 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 **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; + + // 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]; + + fi[0] = fi[1] = fi[2] = 0.0; + fmi[0] = fmi[1] = fmi[2] = 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; + + compute_dipolar(i,j,rij,fmi,spi,spj,r3inv); + if (lattice_flag) compute_dipolar_mech(i,j,rij,fmi,spi,spj,r2inv); + } + + // force accumulation + + f[i][0] += fi[0]; + f[i][1] += fi[1]; + f[i][2] += fi[2]; + fm[i][0] += fmi[0]; + fm[i][1] += fmi[1]; + fm[i][2] += fmi[2]; + + 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 PairSpinDipolarCut::compute_single_pair(int ii, double fmi[3]) +{ + int i,j,jj,jnum,itype,jtype; + double rsq,rinv,r2inv,r3inv; + double xi[3],rij[3]; + double spi[4],spj[4]; + double local_cut2; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **sp = atom->sp; + int *type = atom->type; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // 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]; + + 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]; + + 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; + + // compute dipolar interaction + + compute_dipolar(i,j,rij,fmi,spi,spj,r3inv); + } + } + + //fmi[0] *= mub2mu0hbinv; + //fmi[1] *= mub2mu0hbinv; + //fmi[2] *= mub2mu0hbinv; +} + +/* ---------------------------------------------------------------------- + compute dipolar interaction between spins i and j +------------------------------------------------------------------------- */ + +void PairSpinDipolarCut::compute_dipolar(int i, int j, double rij[3], + double fmi[3], double spi[4], double spj[4], double r3inv) +{ + double sjdotr; + double gigjri2,pre; + + sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2]; + gigjri2 = (spi[3] * spj[3])*r3inv; + pre = mub2mu0hbinv * gigjri2 / 4.0 / MY_PI; + + fmi[0] += pre * gigjri2 * (3.0 * sjdotr *rij[0] - spj[0]); + fmi[1] += pre * gigjri2 * (3.0 * sjdotr *rij[1] - spj[1]); + fmi[2] += pre * gigjri2 * (3.0 * sjdotr *rij[2] - spj[2]); +} + +/* ---------------------------------------------------------------------- + compute the mechanical force due to the dipolar interaction between + atom i and atom j +------------------------------------------------------------------------- */ + +void PairSpinDipolarCut::compute_dipolar_mech(int i, int j, double rij[3], + double fi[3], double spi[3], double spj[3], double r2inv) +{ + double sdots,sidotr,sjdotr,b2,b3; + double gigjri4,bij,pre; + + gigjri4 = (spi[3] * spj[3])/r2inv/r2inv; + 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]; + + bij = sdots - 5.0 * sidotr*sjdotr; + pre = mub2mu0 * bij / 4.0 / MY_PI; + fi[0] += pre * (rij[0] * bij + (sjdotr*spi[0] + sidotr*spj[0])); + fi[1] += pre * (rij[1] * bij + (sjdotr*spi[1] + sidotr*spj[1])); + fi[2] += pre * (rij[2] * bij + (sjdotr*spi[2] + sidotr*spj[2])); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairSpinDipolarCut::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 PairSpinDipolarCut::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 PairSpinDipolarCut::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 PairSpinDipolarCut::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 PairSpinDipolarCut::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_dipolar_cut.h b/src/SPIN/pair_spin_dipolar_cut.h new file mode 100644 index 0000000000..ac17ac2120 --- /dev/null +++ b/src/SPIN/pair_spin_dipolar_cut.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/dipolar/cut,PairSpinDipolarCut) + +#else + +#ifndef LMP_PAIR_SPIN_DIPOLAR_CUT_H +#define LMP_PAIR_SPIN_DIPOLAR_CUT_H + +#include "pair_spin.h" + +namespace LAMMPS_NS { + +class PairSpinDipolarCut : public PairSpin { + public: + double cut_coul; + double **sigma; + + PairSpinDipolarCut(class LAMMPS *); + ~PairSpinDipolarCut(); + 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_dipolar(int, int, double *, double *, double *, + double *, double); + void compute_dipolar_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. + +*/