diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index d548300179..20888290a0 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -16,10 +16,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -102,7 +102,7 @@ void AtomVecSpin::grow_reset() { tag = atom->tag; type = atom->type; mask = atom->mask; image = atom->image; - x = atom->x; v = atom->v; f = atom->f; + x = atom->x; v = atom->v; f = atom->f; sp = atom->sp; fm = atom->fm; } @@ -361,7 +361,7 @@ void AtomVecSpin::unpack_reverse(int n, int *list, double *buf) fm[j][0] += buf[m++]; fm[j][1] += buf[m++]; fm[j][2] += buf[m++]; - } + } } /* ---------------------------------------------------------------------- */ @@ -522,7 +522,7 @@ int AtomVecSpin::pack_border_hybrid(int n, int *list, double *buf) buf[m++] = sp[j][2]; buf[m++] = sp[j][3]; } - + return m; } @@ -552,7 +552,7 @@ void AtomVecSpin::unpack_border(int n, int first, double *buf) for (int iextra = 0; iextra < atom->nextra_border; iextra++) m += modify->fix[atom->extra_border[iextra]]-> unpack_border(n,first,&buf[m]); - + } /* ---------------------------------------------------------------------- */ @@ -584,7 +584,7 @@ void AtomVecSpin::unpack_border_vel(int n, int first, double *buf) for (int iextra = 0; iextra < atom->nextra_border; iextra++) m += modify->fix[atom->extra_border[iextra]]-> unpack_border(n,first,&buf[m]); - + } /* ---------------------------------------------------------------------- */ @@ -601,7 +601,7 @@ int AtomVecSpin::unpack_border_hybrid(int n, int first, double *buf) sp[i][2] = buf[m++]; sp[i][3] = buf[m++]; } - + return m; } @@ -667,7 +667,7 @@ int AtomVecSpin::unpack_exchange(double *buf) unpack_exchange(nlocal,&buf[m]); atom->nlocal++; - + return m; } @@ -784,7 +784,7 @@ void AtomVecSpin::create_atom(int itype, double *coord) mask[nlocal] = 1; image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; - v[nlocal][0] = 0.0; + v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; @@ -843,7 +843,7 @@ void AtomVecSpin::data_atom(double *coord, imageint imagetmp, char **values) int AtomVecSpin::data_atom_hybrid(int nlocal, char **values) { - + sp[nlocal][0] = atof(values[0]); sp[nlocal][1] = atof(values[1]); sp[nlocal][2] = atof(values[2]); @@ -854,7 +854,7 @@ int AtomVecSpin::data_atom_hybrid(int nlocal, char **values) sp[nlocal][1] *= inorm; sp[nlocal][2] *= inorm; sp[nlocal][3] = atof(values[3]); - + return 4; } @@ -878,7 +878,7 @@ void AtomVecSpin::pack_data(double **buf) buf[i][9] = ubuf((image[i] & IMGMASK) - IMGMAX).d; buf[i][10] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; buf[i][11] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; - } + } } /* ---------------------------------------------------------------------- @@ -886,7 +886,7 @@ void AtomVecSpin::pack_data(double **buf) ------------------------------------------------------------------------- */ int AtomVecSpin::pack_data_hybrid(int i, double *buf) -{ +{ buf[0] = sp[i][0]; buf[1] = sp[i][1]; buf[2] = sp[i][2]; @@ -899,7 +899,7 @@ int AtomVecSpin::pack_data_hybrid(int i, double *buf) ------------------------------------------------------------------------- */ void AtomVecSpin::write_data(FILE *fp, int n, double **buf) -{ +{ for (int i = 0; i < n; i++) fprintf(fp,TAGINT_FORMAT \ " %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e " @@ -908,7 +908,7 @@ void AtomVecSpin::write_data(FILE *fp, int n, double **buf) buf[i][2],buf[i][3],buf[i][4], buf[i][5],buf[i][6],buf[i][7],buf[i][8], (int) ubuf(buf[i][10]).i,(int) ubuf(buf[i][11]).i, - (int) ubuf(buf[i][12]).i); + (int) ubuf(buf[i][12]).i); } /* ---------------------------------------------------------------------- @@ -936,7 +936,7 @@ bigint AtomVecSpin::memory_usage() if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3); if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3); if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3); - + if (atom->memcheck("sp")) bytes += memory->usage(sp,nmax,4); if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3); @@ -945,8 +945,8 @@ bigint AtomVecSpin::memory_usage() void AtomVecSpin::force_clear(int n, size_t nbytes) { - memset(&atom->f[0][0],0,3*nbytes); - memset(&atom->fm[0][0],0,3*nbytes); + memset(&atom->f[0][0],0,3*nbytes); + memset(&atom->fm[0][0],0,3*nbytes); } diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h index 99d4a86189..34bc55b99f 100644 --- a/src/SPIN/atom_vec_spin.h +++ b/src/SPIN/atom_vec_spin.h @@ -55,12 +55,12 @@ class AtomVecSpin : public AtomVec { void pack_data(double **); int pack_data_hybrid(int, double *); void write_data(FILE *, int, double **); - int write_data_hybrid(FILE *, double *); + int write_data_hybrid(FILE *, double *); bigint memory_usage(); // clear magnetic and mechanic forces - void force_clear(int, size_t); + void force_clear(int, size_t); private: diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index ab3bab2487..54818a9b70 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -40,7 +40,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) : +ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command"); @@ -52,7 +52,7 @@ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) : init(); allocate(); - + } /* ---------------------------------------------------------------------- */ @@ -79,27 +79,27 @@ void ComputeSpin::compute_vector() double mag[4], magtot[4]; double magenergy, magenergytot; double tempnum, tempnumtot; - double tempdenom, tempdenomtot; + double tempdenom, tempdenomtot; double spintemperature; - + invoked_vector = update->ntimestep; - + countsp = countsptot = 0.0; - mag[0] = mag[1] = mag[2] = mag[3] = 0.0; - magtot[0] = magtot[1] = magtot[2] = magtot[3] = 0.0; - magenergy = magenergytot = 0.0; + mag[0] = mag[1] = mag[2] = mag[3] = 0.0; + magtot[0] = magtot[1] = magtot[2] = magtot[3] = 0.0; + magenergy = magenergytot = 0.0; tempnum = tempnumtot = 0.0; - tempdenom = tempdenomtot = 0.0; - spintemperature = 0.0; + tempdenom = tempdenomtot = 0.0; + spintemperature = 0.0; int *mask = atom->mask; - double **sp = atom->sp; + double **sp = atom->sp; double **fm = atom->fm; double tx,ty,tz; int nlocal = atom->nlocal; - // compute total magnetization and magnetic energy + // compute total magnetization and magnetic energy // compute spin temperature (Nurdin et al., Phys. Rev. E 61, 2000) for (i = 0; i < nlocal; i++) { @@ -108,7 +108,7 @@ void ComputeSpin::compute_vector() mag[0] += sp[i][0]; mag[1] += sp[i][1]; mag[2] += sp[i][2]; - magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]); + magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]); tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1]; ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2]; tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0]; @@ -125,22 +125,22 @@ void ComputeSpin::compute_vector() MPI_Allreduce(&tempnum,&tempnumtot,1,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(&tempdenom,&tempdenomtot,1,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(&countsp,&countsptot,1,MPI_INT,MPI_SUM,world); - + double scale = 1.0/countsptot; magtot[0] *= scale; magtot[1] *= scale; magtot[2] *= scale; magtot[3] = sqrt((magtot[0]*magtot[0])+(magtot[1]*magtot[1])+(magtot[2]*magtot[2])); - spintemperature = hbar*tempnumtot; + spintemperature = hbar*tempnumtot; spintemperature /= (kb*tempdenomtot); - + vector[0] = magtot[0]; vector[1] = magtot[1]; vector[2] = magtot[2]; vector[3] = magtot[3]; - vector[4] = magenergytot*hbar; + vector[4] = magenergytot*hbar; vector[5] = spintemperature; - + } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index d12ec11237..97b33197ce 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -85,7 +85,7 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) : } // initialize Marsaglia RNG with processor-unique seed - + random = new RanPark(lmp,seed + comm->me); } @@ -116,21 +116,21 @@ int FixLangevinSpin::setmask() void FixLangevinSpin::init() { // fix_langevin_spin has to be the last defined fix - + int flag_force = 0; int flag_lang = 0; - for (int i = 0; i < modify->nfix; i++) { + for (int i = 0; i < modify->nfix; i++) { if (strcmp("precession/spin",modify->fix[i]->style)==0) flag_force = MAX(flag_force,i); if (strcmp("langevin/spin",modify->fix[i]->style)==0) flag_lang = i; } - if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin has to come after all other spin fixes"); + if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin has to come after all other spin fixes"); memory->create(spi,3,"langevin:spi"); memory->create(fmi,3,"langevin:fmi"); gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t)); - dts = update->dt; - + dts = update->dt; + double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) double kb = force->boltz; // eV/K D = (MY_2PI*alpha_t*gil_factor*kb*temp); @@ -168,24 +168,24 @@ void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3]) /* ---------------------------------------------------------------------- */ -void FixLangevinSpin::add_temperature(double fmi[3]) +void FixLangevinSpin::add_temperature(double fmi[3]) { double rx = sigma*(2.0*random->uniform() - 1.0); double ry = sigma*(2.0*random->uniform() - 1.0); double rz = sigma*(2.0*random->uniform() - 1.0); - // adding the random field + // adding the random field - fmi[0] += rx; + fmi[0] += rx; fmi[1] += ry; fmi[2] += rz; - - // adding gilbert's prefactor - fmi[0] *= gil_factor; - fmi[1] *= gil_factor; - fmi[2] *= gil_factor; + // adding gilbert's prefactor + + fmi[0] *= gil_factor; + fmi[1] *= gil_factor; + fmi[2] *= gil_factor; } diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h index 0f90a77c14..5358438396 100644 --- a/src/SPIN/fix_langevin_spin.h +++ b/src/SPIN/fix_langevin_spin.h @@ -43,7 +43,7 @@ class FixLangevinSpin : public Fix { double temp; // spin bath temperature double D,sigma; // bath intensity var. double gil_factor; // gilbert's prefactor - + char *id_temp; class Compute *temperature; diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index 699426de9c..d91636af58 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -39,7 +39,7 @@ #include "math_extra.h" #include "math_const.h" #include "memory.h" -#include "modify.h" +#include "modify.h" #include "neighbor.h" #include "neigh_list.h" #include "pair.h" @@ -68,15 +68,15 @@ enum{NONE}; /* ---------------------------------------------------------------------- */ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), + Fix(lmp, narg, arg), pair(NULL), spin_pairs(NULL), - rsec(NULL), stack_head(NULL), stack_foot(NULL), + rsec(NULL), stack_head(NULL), stack_foot(NULL), backward_stacks(NULL), forward_stacks(NULL) { if (lmp->citeme) lmp->citeme->add(cite_fix_nve_spin); - if (narg < 4) error->all(FLERR,"Illegal fix/NVE/spin command"); - + if (narg < 4) error->all(FLERR,"Illegal fix/NVE/spin command"); + time_integrate = 1; sector_flag = NONE; lattice_flag = 1; @@ -90,7 +90,7 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : if (atom->map_style == 0) error->all(FLERR,"Fix NVE/spin requires an atom map, see atom_modify"); - // defining sector_flag + // defining sector_flag int nprocs_tmp = comm->nprocs; if (nprocs_tmp == 1) { @@ -99,10 +99,10 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : sector_flag = 1; } else error->all(FLERR,"Illegal fix/NVE/spin command"); - // defining lattice_flag + // defining lattice_flag int iarg = 3; - while (iarg < narg) { + while (iarg < narg) { if (strcmp(arg[iarg],"lattice") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix/NVE/spin command"); if (strcmp(arg[iarg+1],"no") == 0) lattice_flag = 0; @@ -151,7 +151,7 @@ int FixNVESpin::setmask() mask |= INITIAL_INTEGRATE; mask |= PRE_NEIGHBOR; mask |= FINAL_INTEGRATE; - return mask; + return mask; } /* ---------------------------------------------------------------------- */ @@ -182,7 +182,7 @@ void FixNVESpin::init() npairspin ++; } } - } + } // init length of vector of ptrs to Pair/Spin styles @@ -208,7 +208,7 @@ void FixNVESpin::init() if (count != npairspin) error->all(FLERR,"Incorrect number of spin pairs"); - if (npairspin >= 1) pair_spin_flag = 1; + if (npairspin >= 1) pair_spin_flag = 1; // ptrs FixPrecessionSpin classes @@ -216,29 +216,29 @@ void FixNVESpin::init() for (iforce = 0; iforce < modify->nfix; iforce++) { if (strstr(modify->fix[iforce]->style,"precession/spin")) { precession_spin_flag = 1; - lockprecessionspin = (FixPrecessionSpin *) modify->fix[iforce]; + lockprecessionspin = (FixPrecessionSpin *) modify->fix[iforce]; } } // ptrs on the FixLangevinSpin class - + for (iforce = 0; iforce < modify->nfix; iforce++) { if (strstr(modify->fix[iforce]->style,"langevin/spin")) { maglangevin_flag = 1; - locklangevinspin = (FixLangevinSpin *) modify->fix[iforce]; - } + locklangevinspin = (FixLangevinSpin *) modify->fix[iforce]; + } } if (maglangevin_flag) { if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1; if (locklangevinspin->temp_flag == 1) temp_flag = 1; } - + // setting the sector variables/lists nsectors = 0; memory->create(rsec,3,"NVE/spin:rsec"); - + // perform the sectoring operation if (sector_flag) sectoring(); @@ -264,20 +264,20 @@ void FixNVESpin::initial_integrate(int vflag) double **x = atom->x; double **v = atom->v; double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; + double *rmass = atom->rmass; + double *mass = atom->mass; int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; int *type = atom->type; - int *mask = atom->mask; + int *mask = atom->mask; // update half v for all atoms - + if (lattice_flag) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (rmass) dtfm = dtf / rmass[i]; - else dtfm = dtf / mass[type[i]]; + else dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -297,7 +297,7 @@ void FixNVESpin::initial_integrate(int vflag) i = forward_stacks[i]; } } - for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal + for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal comm->forward_comm(); int i = stack_head[j]; while (i >= 0) { @@ -342,7 +342,7 @@ void FixNVESpin::initial_integrate(int vflag) i = forward_stacks[i]; } } - for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal + for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal comm->forward_comm(); int i = stack_head[j]; while (i >= 0) { @@ -365,7 +365,7 @@ void FixNVESpin::initial_integrate(int vflag) } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- setup pre_neighbor() ---------------------------------------------------------------------- */ @@ -374,7 +374,7 @@ void FixNVESpin::setup_pre_neighbor() pre_neighbor(); } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- store in two linked lists the advance order of the spins (sectoring) ---------------------------------------------------------------------- */ @@ -414,7 +414,7 @@ void FixNVESpin::pre_neighbor() } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- compute the magnetic torque for the spin ii ---------------------------------------------------------------------- */ @@ -430,7 +430,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i) spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + fmi[0] = fmi[1] = fmi[2] = 0.0; // update magnetic pair interactions @@ -440,10 +440,10 @@ void FixNVESpin::ComputeInteractionsSpin(int i) spin_pairs[k]->compute_single_pair(i,fmi); } } - + // update magnetic precession interactions - if (precession_spin_flag) { + if (precession_spin_flag) { lockprecessionspin->compute_single_precession(i,spi,fmi); } @@ -451,11 +451,11 @@ void FixNVESpin::ComputeInteractionsSpin(int i) if (maglangevin_flag) { // mag. langevin if (tdamp_flag) { // transverse damping - locklangevinspin->add_tdamping(spi,fmi); + locklangevinspin->add_tdamping(spi,fmi); } if (temp_flag) { // spin temperature locklangevinspin->add_temperature(fmi); - } + } } // replace the magnetic force fm[i] by its new value fmi @@ -466,7 +466,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i) } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- divide each domain into 8 sectors ---------------------------------------------------------------------- */ @@ -481,9 +481,9 @@ void FixNVESpin::sectoring() subhi[dim]=subhitmp[dim]; } - const double rsx = subhi[0] - sublo[0]; - const double rsy = subhi[1] - sublo[1]; - const double rsz = subhi[2] - sublo[2]; + const double rsx = subhi[0] - sublo[0]; + const double rsy = subhi[1] - sublo[1]; + const double rsz = subhi[2] - sublo[2]; // extract larger cutoff from PairSpin styles @@ -498,10 +498,10 @@ void FixNVESpin::sectoring() if (rv == 0.0) error->all(FLERR,"Illegal sectoring operation"); - double rax = rsx/rv; - double ray = rsy/rv; - double raz = rsz/rv; - + double rax = rsx/rv; + double ray = rsy/rv; + double raz = rsz/rv; + sec[0] = 1; sec[1] = 1; sec[2] = 1; @@ -511,7 +511,7 @@ void FixNVESpin::sectoring() nsectors = sec[0]*sec[1]*sec[2]; - if (sector_flag == 1 && nsectors != 8) + if (sector_flag == 1 && nsectors != 8) error->all(FLERR,"Illegal sectoring operation"); rsec[0] = rsx; @@ -523,7 +523,7 @@ void FixNVESpin::sectoring() } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- define sector for an atom at a position x[i] ---------------------------------------------------------------------- */ @@ -541,12 +541,12 @@ int FixNVESpin::coords2sector(double *x) seci[1] = x[1] > (sublo[1] + rsec[1]); seci[2] = x[2] > (sublo[2] + rsec[2]); - nseci = (seci[0] + 2*seci[1] + 4*seci[2]); + nseci = (seci[0] + 2*seci[1] + 4*seci[2]); return nseci; } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- advance the spin i of a timestep dts ---------------------------------------------------------------------- */ @@ -557,7 +557,7 @@ void FixNVESpin::AdvanceSingleSpin(int i) double **sp = atom->sp; double **fm = atom->fm; double msq,scale,fm2,energy,dts2; - double cp[3],g[3]; + double cp[3],g[3]; cp[0] = cp[1] = cp[2] = 0.0; g[0] = g[1] = g[2] = 0.0; @@ -572,18 +572,18 @@ void FixNVESpin::AdvanceSingleSpin(int i) g[0] = sp[i][0]+cp[0]*dts; g[1] = sp[i][1]+cp[1]*dts; g[2] = sp[i][2]+cp[2]*dts; - + g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2; g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2; g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2; - + g[0] /= (1+0.25*fm2*dts2); g[1] /= (1+0.25*fm2*dts2); g[2] /= (1+0.25*fm2*dts2); - + sp[i][0] = g[0]; sp[i][1] = g[1]; - sp[i][2] = g[2]; + sp[i][2] = g[2]; // renormalization (check if necessary) @@ -618,11 +618,11 @@ void FixNVESpin::final_integrate() double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; - double *mass = atom->mass; + double *mass = atom->mass; int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; int *type = atom->type; - int *mask = atom->mask; + int *mask = atom->mask; // update half v for all particles @@ -630,10 +630,10 @@ void FixNVESpin::final_integrate() for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (rmass) dtfm = dtf / rmass[i]; - else dtfm = dtf / mass[type[i]]; + else dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; + v[i][2] += dtfm * f[i][2]; } } } diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 70781c6d8c..afc1db14d6 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -36,17 +36,17 @@ friend class PairSpin; virtual void initial_integrate(int); virtual void final_integrate(); - void ComputeInteractionsSpin(int); // compute and advance single spin functions + void ComputeInteractionsSpin(int); // compute and advance single spin functions void AdvanceSingleSpin(int); - void sectoring(); // sectoring operation functions + void sectoring(); // sectoring operation functions int coords2sector(double *); void setup_pre_neighbor(); void pre_neighbor(); int lattice_flag; // lattice_flag = 0 if spins only - // lattice_flag = 1 if spin-lattice calc. + // lattice_flag = 1 if spin-lattice calc. protected: @@ -54,7 +54,7 @@ friend class PairSpin; // sector_flag = 1 if parallel algorithm double dtv, dtf, dts; // velocity, force, and spin timesteps - + int nlocal_max; // max value of nlocal (for lists size) int pair_spin_flag; // magnetic pair flags @@ -63,24 +63,24 @@ friend class PairSpin; int tdamp_flag, temp_flag; // pointers to magnetic fixes - + class FixPrecessionSpin *lockprecessionspin; - class FixLangevinSpin *locklangevinspin; + class FixLangevinSpin *locklangevinspin; // pointers to magnetic pair styles - + int npairs, npairspin; // # of pairs, and # of spin pairs class Pair *pair; class PairSpin **spin_pairs; // vector of spin pairs - + // sectoring variables - + int nsectors; double *rsec; // stacking variables for sectoring algorithm - - int *stack_head; // index of first atom in backward_stacks + + int *stack_head; // index of first atom in backward_stacks int *stack_foot; // index of first atom in forward_stacks int *backward_stacks; // index of next atom in backward stack int *forward_stacks; // index of next atom in forward stack @@ -96,7 +96,7 @@ friend class PairSpin; E: Illegal fix NVE/spin command -Self-explanatory. Check the input script syntax and compare to the +Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. @@ -106,7 +106,7 @@ An atom/spin style with this attribute is needed. E: Illegal sectoring operation -The number of processes does not match the size of the system. +The number of processes does not match the size of the system. See the documentation of the sectoring method. */ diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index 412b9d903b..9b33fac551 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -50,7 +50,7 @@ enum{CONSTANT,EQUAL}; FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 7) error->all(FLERR,"Illegal precession/spin command"); - + // magnetic interactions coded for cartesian coordinates hbar = force->hplanck/MY_2PI; @@ -61,10 +61,10 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm extscalar = 1; respa_level_support = 1; ilevel_respa = 0; - + magstr = NULL; magfieldstyle = CONSTANT; - + H_field = 0.0; nhx = nhy = nhz = 0.0; hx = hy = hz = 0.0; @@ -99,7 +99,7 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm time_origin = update->ntimestep; eflag = 0; - emag = 0.0; + emag = 0.0; } /* ---------------------------------------------------------------------- */ @@ -126,8 +126,8 @@ int FixPrecessionSpin::setmask() void FixPrecessionSpin::init() { const double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) - const double mub = 5.78901e-5; // in eV/T - const double gyro = mub/hbar; // in rad.THz/T + const double mub = 5.78901e-5; // in eV/T + const double gyro = mub/hbar; // in rad.THz/T H_field *= gyro; // in rad.THz Ka /= hbar; // in rad.THz @@ -139,19 +139,19 @@ void FixPrecessionSpin::init() if (magstr) { magvar = input->variable->find(magstr); - if (magvar < 0) + if (magvar < 0) error->all(FLERR,"Illegal precession/spin command"); if (!input->variable->equalstyle(magvar)) error->all(FLERR,"Illegal precession/spin command"); } - + varflag = CONSTANT; if (magfieldstyle != CONSTANT) varflag = EQUAL; - + // set magnetic field components if (varflag == CONSTANT) set_magneticprecession(); - + } /* ---------------------------------------------------------------------- */ @@ -179,10 +179,10 @@ void FixPrecessionSpin::post_force(int vflag) set_magneticprecession(); // update mag. field if time-dep. } - double **sp = atom->sp; + double **sp = atom->sp; double **fm = atom->fm; - double spi[3], fmi[3]; - const int nlocal = atom->nlocal; + double spi[3], fmi[3]; + const int nlocal = atom->nlocal; eflag = 0; emag = 0.0; @@ -192,13 +192,13 @@ void FixPrecessionSpin::post_force(int vflag) spi[1] = sp[i][1]; spi[2] = sp[i][2]; fmi[0] = fmi[1] = fmi[2] = 0.0; - + if (zeeman_flag) { // compute Zeeman interaction compute_zeeman(i,fmi); emag -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); emag *= hbar; - } - + } + if (aniso_flag) { // compute magnetic anisotropy compute_anisotropy(spi,fmi); emag -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); @@ -228,7 +228,7 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f void FixPrecessionSpin::compute_zeeman(int i, double fmi[3]) { - double **sp = atom->sp; + double **sp = atom->sp; fmi[0] -= sp[i][3]*hx; fmi[1] -= sp[i][3]*hy; fmi[2] -= sp[i][3]*hz; @@ -258,12 +258,12 @@ void FixPrecessionSpin::set_magneticprecession() if (zeeman_flag) { hx = H_field*nhx; hy = H_field*nhy; - hz = H_field*nhz; + hz = H_field*nhz; } if (aniso_flag) { Kax = 2.0*Ka*nax; Kay = 2.0*Ka*nay; - Kaz = 2.0*Ka*naz; + Kaz = 2.0*Ka*naz; } } @@ -274,7 +274,7 @@ void FixPrecessionSpin::set_magneticprecession() double FixPrecessionSpin::compute_scalar() { // only sum across procs one time - + if (eflag == 0) { MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world); eflag = 1; diff --git a/src/SPIN/fix_precession_spin.h b/src/SPIN/fix_precession_spin.h index 9e7bd1a830..53ae4ba124 100644 --- a/src/SPIN/fix_precession_spin.h +++ b/src/SPIN/fix_precession_spin.h @@ -38,13 +38,13 @@ class FixPrecessionSpin : public Fix { double compute_scalar(); int zeeman_flag, aniso_flag; - void compute_single_precession(int, double *, double *); + void compute_single_precession(int, double *, double *); void compute_zeeman(int, double *); void compute_anisotropy(double *, double *); protected: int style; // style of the magnetic precession - + double degree2rad; double hbar; int ilevel_respa; @@ -56,21 +56,21 @@ class FixPrecessionSpin : public Fix { int magfieldstyle; int magvar; char *magstr; - + // zeeman field intensity and direction - double H_field; + double H_field; double nhx, nhy, nhz; double hx, hy, hz; // temp. force variables - + // magnetic anisotropy intensity and direction - - double Ka; + + double Ka; double nax, nay, naz; double Kax, Kay, Kaz; // temp. force variables void set_magneticprecession(); - + }; } @@ -86,7 +86,7 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. -precession/spin fix command has 7 arguments: -fix ID group precession/spin magnitude (T or eV) style (zeeman or anisotropy) -direction (3 cartesian coordinates) +precession/spin fix command has 7 arguments: +fix ID group precession/spin magnitude (T or eV) style (zeeman or anisotropy) +direction (3 cartesian coordinates) */ diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index e8a4c126da..acb7ffe548 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index 94f8433108..100eec1732 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -39,7 +39,7 @@ friend class FixNVESpin; virtual void compute_single_pair(int, double *) {} protected: - double hbar; // Planck constant (eV.ps.rad-1) + double hbar; // Planck constant (eV.ps.rad-1) virtual void allocate() {} }; diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp index 5f735cd1a8..07ae684939 100755 --- a/src/SPIN/pair_spin_dmi.cpp +++ b/src/SPIN/pair_spin_dmi.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -80,9 +80,9 @@ void PairSpinDmi::settings(int narg, char **arg) if (strcmp(update->unit_style,"metal") != 0) error->all(FLERR,"Spin simulations require metal unit style"); - + cut_spin_dmi_global = force->numeric(FLERR,arg[0]); - + // reset cutoffs that have been explicitly set if (allocated) { @@ -95,7 +95,7 @@ void PairSpinDmi::settings(int narg, char **arg) } } } - + } /* ---------------------------------------------------------------------- @@ -106,28 +106,28 @@ void PairSpinDmi::coeff(int narg, char **arg) { if (!allocated) allocate(); - // check if args correct - + // check if args correct + if (strcmp(arg[2],"dmi") != 0) error->all(FLERR,"Incorrect args in pair_style command"); if (narg != 8) 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); - + const double rij = force->numeric(FLERR,arg[3]); const double dm = (force->numeric(FLERR,arg[4]))/hbar; - double dmx = force->numeric(FLERR,arg[5]); - double dmy = force->numeric(FLERR,arg[6]); - double dmz = force->numeric(FLERR,arg[7]); + double dmx = force->numeric(FLERR,arg[5]); + double dmy = force->numeric(FLERR,arg[6]); + double dmz = force->numeric(FLERR,arg[7]); double inorm = 1.0/(dmx*dmx+dmy*dmy+dmz*dmz); - dmx *= inorm; - dmy *= inorm; - dmz *= inorm; - + dmx *= inorm; + dmy *= inorm; + dmz *= inorm; + int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { @@ -140,8 +140,8 @@ void PairSpinDmi::coeff(int narg, char **arg) count++; } } - if (count == 0) - error->all(FLERR,"Incorrect args in pair_style command"); + if (count == 0) + error->all(FLERR,"Incorrect args in pair_style command"); } @@ -169,7 +169,7 @@ void PairSpinDmi::init_style() } 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++) { @@ -187,7 +187,7 @@ void PairSpinDmi::init_style() double PairSpinDmi::init_one(int i, int j) { - + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); return cut_spin_dmi_global; @@ -215,20 +215,20 @@ void PairSpinDmi::compute(int eflag, int vflag) double fi[3], fmi[3]; double local_cut2; double rsq, inorm; - int *ilist,*jlist,*numneigh,**firstneigh; + 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 *type = atom->type; + int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -240,28 +240,28 @@ void PairSpinDmi::compute(int eflag, int vflag) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; itype = type[i]; - + jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; 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[0] = sp[i][0]; + spi[1] = sp[i][1]; spi[2] = sp[i][2]; - - + + // loop on neighbors 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[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; evdwl = 0.0; @@ -269,11 +269,11 @@ void PairSpinDmi::compute(int eflag, int vflag) fmi[0] = fmi[1] = fmi[2] = 0.0; rij[0] = rij[1] = rij[2] = 0.0; eij[0] = eij[1] = eij[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]; + rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; inorm = 1.0/sqrt(rsq); eij[0] = rij[0]*inorm; eij[1] = rij[1]*inorm; @@ -282,7 +282,7 @@ void PairSpinDmi::compute(int eflag, int vflag) local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype]; // compute dmi interaction - + if (rsq <= local_cut2) { compute_dmi(i,j,eij,fmi,spj); if (lattice_flag) { @@ -290,16 +290,16 @@ void PairSpinDmi::compute(int eflag, int vflag) } } - f[i][0] += fi[0]; - f[i][1] += fi[1]; + 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][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][0] -= fi[0]; + f[j][1] -= fi[1]; f[j][2] -= fi[2]; } @@ -311,7 +311,7 @@ void PairSpinDmi::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); } - } + } if (vflag_fdotr) virial_fdotr_compute(); @@ -340,11 +340,11 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) i = ilist[ii]; itype = type[i]; - + xi[0] = x[i][0]; xi[1] = x[i][1]; xi[2] = x[i][2]; - + jlist = firstneigh[i]; jnum = numneigh[i]; @@ -383,7 +383,7 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; double dmix, dmiy, dmiz; itype = type[i]; @@ -422,15 +422,15 @@ void PairSpinDmi::allocate() for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; - + memory->create(cut_spin_dmi,n+1,n+1,"pair:cut_spin_dmi"); memory->create(DM,n+1,n+1,"pair:DM"); memory->create(v_dmx,n+1,n+1,"pair:DM_vector_x"); memory->create(v_dmy,n+1,n+1,"pair:DM_vector_y"); memory->create(v_dmz,n+1,n+1,"pair:DM_vector_z"); - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + } @@ -441,7 +441,7 @@ void PairSpinDmi::allocate() void PairSpinDmi::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++) { @@ -490,7 +490,7 @@ void PairSpinDmi::read_restart(FILE *fp) } } - + /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -515,5 +515,5 @@ void PairSpinDmi::read_restart_settings(FILE *fp) } MPI_Bcast(&cut_spin_dmi_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/SPIN/pair_spin_dmi.h b/src/SPIN/pair_spin_dmi.h index cbd234c192..a309f0c8d5 100755 --- a/src/SPIN/pair_spin_dmi.h +++ b/src/SPIN/pair_spin_dmi.h @@ -39,16 +39,16 @@ class PairSpinDmi : public PairSpin { void compute_dmi(int, int, double *, double *, double *); void compute_dmi_mech(double *); - + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - + double cut_spin_dmi_global; // short range pair cutoff protected: - double **DM; // dmi coeff in eV + double **DM; // dmi coeff in eV double **v_dmx, **v_dmy, **v_dmz; // dmi direction double **cut_spin_dmi; // cutoff distance dmi diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index 39c5823d03..7b05d7337e 100755 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -64,7 +64,7 @@ PairSpinExchange::~PairSpinExchange() memory->destroy(J1_mag); memory->destroy(J1_mech); memory->destroy(J2); - memory->destroy(J3); + memory->destroy(J3); memory->destroy(cutsq); // to be deleted } } @@ -72,7 +72,7 @@ PairSpinExchange::~PairSpinExchange() /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ - + void PairSpinExchange::settings(int narg, char **arg) { if (narg < 1 || narg > 2) @@ -80,9 +80,9 @@ void PairSpinExchange::settings(int narg, char **arg) if (strcmp(update->unit_style,"metal") != 0) error->all(FLERR,"Spin simulations require metal unit style"); - + cut_spin_exchange_global = force->numeric(FLERR,arg[0]); - + // reset cutoffs that have been explicitly set if (allocated) { @@ -93,7 +93,7 @@ void PairSpinExchange::settings(int narg, char **arg) cut_spin_exchange[i][j] = cut_spin_exchange_global; } } - + } /* ---------------------------------------------------------------------- @@ -103,29 +103,29 @@ void PairSpinExchange::settings(int narg, char **arg) void PairSpinExchange::coeff(int narg, char **arg) { if (!allocated) allocate(); - + // check if args correct if (strcmp(arg[2],"exchange") != 0) error->all(FLERR,"Incorrect args in pair_style command"); - if (narg != 7) + if (narg != 7) 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); - + // get exchange arguments from input command const double rc = force->numeric(FLERR,arg[3]); const double j1 = force->numeric(FLERR,arg[4]); - const double j2 = force->numeric(FLERR,arg[5]); - const double j3 = force->numeric(FLERR,arg[6]); - + const double j2 = force->numeric(FLERR,arg[5]); + const double j3 = force->numeric(FLERR,arg[6]); + int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - cut_spin_exchange[i][j] = rc; + cut_spin_exchange[i][j] = rc; J1_mag[i][j] = j1/hbar; J1_mech[i][j] = j1; J2[i][j] = j2; @@ -134,12 +134,12 @@ void PairSpinExchange::coeff(int narg, char **arg) count++; } } - if (count == 0) - error->all(FLERR,"Incorrect args in pair_style command"); + if (count == 0) + error->all(FLERR,"Incorrect args in pair_style command"); } /* ---------------------------------------------------------------------- - init specific to this pair style + init specific to this pair style ------------------------------------------------------------------------- */ void PairSpinExchange::init_style() @@ -180,7 +180,7 @@ void PairSpinExchange::init_style() double PairSpinExchange::init_one(int i, int j) { - + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); return cut_spin_exchange_global; @@ -201,14 +201,14 @@ void *PairSpinExchange::extract(const char *str, int &dim) void PairSpinExchange::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double evdwl, ecoul; double xi[3], rij[3], eij[3]; double spi[3], spj[3]; double fi[3], fmi[3]; double local_cut2; double rsq, inorm; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist,*jlist,*numneigh,**firstneigh; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -218,10 +218,10 @@ void PairSpinExchange::compute(int eflag, int vflag) double **f = atom->f; double **fm = atom->fm; double **sp = atom->sp; - int *type = atom->type; - int nlocal = atom->nlocal; + int *type = atom->type; + int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -233,16 +233,16 @@ void PairSpinExchange::compute(int eflag, int vflag) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; itype = type[i]; - + jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; 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[0] = sp[i][0]; + spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + // loop on neighbors for (jj = 0; jj < jnum; jj++) { @@ -250,28 +250,28 @@ void PairSpinExchange::compute(int eflag, int vflag) j &= NEIGHMASK; jtype = type[j]; - spj[0] = sp[j][0]; - spj[1] = sp[j][1]; - spj[2] = sp[j][2]; + spj[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; evdwl = 0.0; 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]; inorm = 1.0/sqrt(rsq); - eij[0] = inorm*rij[0]; - eij[1] = inorm*rij[1]; - eij[2] = inorm*rij[2]; + eij[0] = inorm*rij[0]; + eij[1] = inorm*rij[1]; + eij[2] = inorm*rij[2]; local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype]; // compute exchange interaction - + if (rsq <= local_cut2) { compute_exchange(i,j,rsq,fmi,spj); if (lattice_flag) { @@ -279,16 +279,16 @@ void PairSpinExchange::compute(int eflag, int vflag) } } - f[i][0] += fi[0]; - f[i][1] += fi[1]; + 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][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][0] -= fi[0]; + f[j][1] -= fi[1]; f[j][2] -= fi[2]; } @@ -300,17 +300,17 @@ void PairSpinExchange::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); } - } + } if (vflag_fdotr) virial_fdotr_compute(); - + } /* ---------------------------------------------------------------------- update the pair interactions fmi acting on the spin ii ------------------------------------------------------------------------- */ -void PairSpinExchange::compute_single_pair(int ii, double fmi[3]) +void PairSpinExchange::compute_single_pair(int ii, double fmi[3]) { int *type = atom->type; @@ -336,7 +336,7 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3]) xi[0] = x[i][0]; xi[1] = x[i][1]; xi[2] = x[i][2]; - + jlist = firstneigh[i]; jnum = numneigh[i]; @@ -370,13 +370,13 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3]) void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; double Jex, ra; itype = type[i]; jtype = type[j]; - ra = rsq/J3[itype][jtype]/J3[itype][jtype]; + ra = rsq/J3[itype][jtype]/J3[itype][jtype]; Jex = 4.0*J1_mag[itype][jtype]*ra; Jex *= (1.0-J2[itype][jtype]*ra); Jex *= exp(-ra); @@ -392,21 +392,21 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3], void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double rij[3], double fi[3], double spi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; double Jex, Jex_mech, ra, rr, iJ3; itype = type[i]; jtype = type[j]; - + Jex = J1_mech[itype][jtype]; iJ3 = 1.0/(J3[itype][jtype]*J3[itype][jtype]); - ra = rsq*iJ3; + ra = rsq*iJ3; rr = sqrt(rsq)*iJ3; Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra); Jex_mech *= 8.0*Jex*rr*exp(-ra); - Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]); + Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]); fi[0] -= Jex_mech*rij[0]; fi[1] -= Jex_mech*rij[1]; @@ -426,13 +426,13 @@ void PairSpinExchange::allocate() for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; - + memory->create(cut_spin_exchange,n+1,n+1,"pair/spin/exchange:cut_spin_exchange"); memory->create(J1_mag,n+1,n+1,"pair/spin/exchange:J1_mag"); memory->create(J1_mech,n+1,n+1,"pair/spin/exchange:J1_mech"); - memory->create(J2,n+1,n+1,"pair/spin/exchange:J2"); + memory->create(J2,n+1,n+1,"pair/spin/exchange:J2"); memory->create(J3,n+1,n+1,"pair/spin/exchange:J3"); - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); } @@ -444,7 +444,7 @@ void PairSpinExchange::allocate() void PairSpinExchange::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++) { @@ -494,7 +494,7 @@ void PairSpinExchange::read_restart(FILE *fp) } } - + /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -519,6 +519,6 @@ void PairSpinExchange::read_restart_settings(FILE *fp) } MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/SPIN/pair_spin_exchange.h b/src/SPIN/pair_spin_exchange.h index 89b35c1bd8..b524a513eb 100755 --- a/src/SPIN/pair_spin_exchange.h +++ b/src/SPIN/pair_spin_exchange.h @@ -34,12 +34,12 @@ class PairSpinExchange : public PairSpin { double init_one(int, int); void *extract(const char *, int &); - void compute(int, int); + void compute(int, int); void compute_single_pair(int, double *); void compute_exchange(int, int, double, double *, double *); void compute_exchange_mech(int, int, double, double *, double *, double *, double *); - + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); @@ -49,7 +49,7 @@ class PairSpinExchange : public PairSpin { protected: double **J1_mag; // exchange coeffs in eV - double **J1_mech; // mech exchange coeffs in + double **J1_mech; // mech exchange coeffs in double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang double **cut_spin_exchange; // cutoff distance exchange diff --git a/src/SPIN/pair_spin_magelec.cpp b/src/SPIN/pair_spin_magelec.cpp index 7b6fd32333..b8ef9db609 100755 --- a/src/SPIN/pair_spin_magelec.cpp +++ b/src/SPIN/pair_spin_magelec.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -47,7 +47,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ PairSpinMagelec::PairSpinMagelec(LAMMPS *lmp) : PairSpin(lmp), -lockfixnvespin(NULL) +lockfixnvespin(NULL) { single_enable = 0; no_virial_fdotr_compute = 1; @@ -81,9 +81,9 @@ void PairSpinMagelec::settings(int narg, char **arg) if (strcmp(update->unit_style,"metal") != 0) error->all(FLERR,"Spin simulations require metal unit style"); - + cut_spin_magelec_global = force->numeric(FLERR,arg[0]); - + // reset cutoffs that have been explicitly set if (allocated) { @@ -94,7 +94,7 @@ void PairSpinMagelec::settings(int narg, char **arg) cut_spin_magelec[i][j] = cut_spin_magelec_global; } } - + } /* ---------------------------------------------------------------------- @@ -105,27 +105,27 @@ void PairSpinMagelec::coeff(int narg, char **arg) { if (!allocated) allocate(); - // check if args correct - + // check if args correct + if (strcmp(arg[2],"magelec") != 0) error->all(FLERR,"Incorrect args in pair_style command"); if (narg != 8) 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); - + const double rij = force->numeric(FLERR,arg[3]); const double magelec = (force->numeric(FLERR,arg[4])); - double mex = force->numeric(FLERR,arg[5]); - double mey = force->numeric(FLERR,arg[6]); - double mez = force->numeric(FLERR,arg[7]); + double mex = force->numeric(FLERR,arg[5]); + double mey = force->numeric(FLERR,arg[6]); + double mez = force->numeric(FLERR,arg[7]); double inorm = 1.0/(mex*mex+mey*mey+mez*mez); - mex *= inorm; - mey *= inorm; - mez *= inorm; + mex *= inorm; + mey *= inorm; + mez *= inorm; int count = 0; for (int i = ilo; i <= ihi; i++) { @@ -140,7 +140,7 @@ void PairSpinMagelec::coeff(int narg, char **arg) count++; } } - if (count == 0) + if (count == 0) error->all(FLERR,"Incorrect args in pair_style command"); } @@ -155,7 +155,7 @@ void PairSpinMagelec::init_style() 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; @@ -169,7 +169,7 @@ void PairSpinMagelec::init_style() } 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++) { @@ -187,7 +187,7 @@ void PairSpinMagelec::init_style() double PairSpinMagelec::init_one(int i, int j) { - + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); return cut_spin_magelec_global; @@ -209,27 +209,27 @@ void *PairSpinMagelec::extract(const char *str, int &dim) void PairSpinMagelec::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double evdwl, ecoul; double xi[3], rij[3], eij[3]; double spi[3], spj[3]; double fi[3], fmi[3]; double local_cut2; double rsq, inorm; - int *ilist,*jlist,*numneigh,**firstneigh; + 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 *type = atom->type; + int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -243,14 +243,14 @@ void PairSpinMagelec::compute(int eflag, int vflag) itype = type[i]; jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; 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[0] = sp[i][0]; + spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + // loop on neighbors for (jj = 0; jj < jnum; jj++) { @@ -258,19 +258,19 @@ void PairSpinMagelec::compute(int eflag, int vflag) j &= NEIGHMASK; jtype = type[j]; - spj[0] = sp[j][0]; - spj[1] = sp[j][1]; - spj[2] = sp[j][2]; + spj[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; evdwl = 0.0; 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]; + rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; inorm = 1.0/sqrt(rsq); eij[0] = inorm*rij[0]; eij[1] = inorm*rij[1]; @@ -287,19 +287,19 @@ void PairSpinMagelec::compute(int eflag, int vflag) } } - f[i][0] += fi[0]; - f[i][1] += fi[1]; + 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][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][0] -= fi[0]; + f[j][1] -= fi[1]; f[j][2] -= fi[2]; } - + if (eflag) { evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl *= hbar; @@ -308,15 +308,15 @@ void PairSpinMagelec::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); } - } + } if (vflag_fdotr) virial_fdotr_compute(); - + } /* ---------------------------------------------------------------------- */ -void PairSpinMagelec::compute_single_pair(int ii, double fmi[3]) +void PairSpinMagelec::compute_single_pair(int ii, double fmi[3]) { int *type = atom->type; double **x = atom->x; @@ -343,7 +343,7 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3]) xi[2] = x[i][2]; eij[0] = eij[1] = eij[2] = 0.0; - + jlist = firstneigh[i]; jnum = numneigh[i]; @@ -376,36 +376,36 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3]) /* ---------------------------------------------------------------------- */ -void PairSpinMagelec::compute_magelec(int i, int j, double rsq, double eij[3], double fmi[3], double spj[3]) +void PairSpinMagelec::compute_magelec(int i, int j, double rsq, double eij[3], double fmi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; itype = type[i]; jtype = type[j]; - - double local_cut2 = cut_spin_magelec[itype][jtype]*cut_spin_magelec[itype][jtype]; - + + double local_cut2 = cut_spin_magelec[itype][jtype]*cut_spin_magelec[itype][jtype]; + if (rsq <= local_cut2) { double meix,meiy,meiz; double rx, ry, rz; double vx, vy, vz; - + rx = eij[0]; - ry = eij[1]; - rz = eij[2]; - + ry = eij[1]; + rz = eij[2]; + vx = v_mex[itype][jtype]; vy = v_mey[itype][jtype]; vz = v_mez[itype][jtype]; - - meix = vy*rz - vz*ry; - meiy = vz*rx - vx*rz; - meiz = vx*ry - vy*rx; - - meix *= ME[itype][jtype]; - meiy *= ME[itype][jtype]; - meiz *= ME[itype][jtype]; - + + meix = vy*rz - vz*ry; + meiy = vz*rx - vx*rz; + meiz = vx*ry - vy*rx; + + meix *= ME[itype][jtype]; + meiy *= ME[itype][jtype]; + meiz *= ME[itype][jtype]; + fmi[0] += spj[1]*meiz - spj[2]*meiy; fmi[1] += spj[2]*meix - spj[0]*meiz; fmi[2] += spj[0]*meiy - spj[1]*meix; @@ -416,7 +416,7 @@ void PairSpinMagelec::compute_magelec(int i, int j, double rsq, double eij[3], d void PairSpinMagelec::compute_magelec_mech(int i, int j, double fi[3], double spi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; itype = type[i]; jtype = type[j]; @@ -432,16 +432,16 @@ void PairSpinMagelec::compute_magelec_mech(int i, int j, double fi[3], double sp meiy = spi[2]*spi[0] - spi[0]*spj[2]; meiz = spi[0]*spi[1] - spi[1]*spj[0]; - meix *= ME_mech[itype][jtype]; - meiy *= ME_mech[itype][jtype]; - meiz *= ME_mech[itype][jtype]; + meix *= ME_mech[itype][jtype]; + meiy *= ME_mech[itype][jtype]; + meiz *= ME_mech[itype][jtype]; fi[0] = meiy*vz - meiz*vy; fi[1] = meiz*vx - meix*vz; fi[2] = meix*vy - meiy*vx; } - + /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ @@ -455,14 +455,14 @@ void PairSpinMagelec::allocate() for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; - + memory->create(cut_spin_magelec,n+1,n+1,"pair/spin/me:cut_spin_magelec"); memory->create(ME,n+1,n+1,"pair/spin/me:ME"); memory->create(ME_mech,n+1,n+1,"pair/spin/me:ME_mech"); memory->create(v_mex,n+1,n+1,"pair/spin/me:ME_vector_x"); memory->create(v_mey,n+1,n+1,"pair/spin/me:ME_vector_y"); memory->create(v_mez,n+1,n+1,"pair/spin/me:ME_vector_z"); - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); } /* ---------------------------------------------------------------------- @@ -472,7 +472,7 @@ void PairSpinMagelec::allocate() void PairSpinMagelec::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++) { @@ -545,5 +545,5 @@ void PairSpinMagelec::read_restart_settings(FILE *fp) } MPI_Bcast(&cut_spin_magelec_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/SPIN/pair_spin_magelec.h b/src/SPIN/pair_spin_magelec.h index bc1e3a6296..ce13476271 100755 --- a/src/SPIN/pair_spin_magelec.h +++ b/src/SPIN/pair_spin_magelec.h @@ -37,20 +37,20 @@ class PairSpinMagelec : public PairSpin { void compute(int, int); void compute_single_pair(int, double *); - void compute_magelec(int, int, double, double *, double *, double *); - void compute_magelec_mech(int, int, double *, double *, double *); - + void compute_magelec(int, int, double, double *, double *, double *); + void compute_magelec_mech(int, int, double *, double *, double *); + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - + double cut_spin_magelec_global; // global me cutoff protected: double **ME, **ME_mech; // magelec coeff in eV double **v_mex, **v_mey, **v_mez; // magelec direction - double **cut_spin_magelec; // magelec cutoff distance + double **cut_spin_magelec; // magelec cutoff distance int lattice_flag; // flag for mech force computation class FixNVESpin *lockfixnvespin; // ptr to FixNVESpin for setups diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp index 36deab81f1..0160d8a69b 100755 --- a/src/SPIN/pair_spin_neel.cpp +++ b/src/SPIN/pair_spin_neel.cpp @@ -14,10 +14,10 @@ /* ------------------------------------------------------------------------ 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 + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. arXiv preprint arXiv:1801.10233. ------------------------------------------------------------------------- */ @@ -65,11 +65,11 @@ PairSpinNeel::~PairSpinNeel() memory->destroy(g1); memory->destroy(g1_mech); memory->destroy(g2); - memory->destroy(g3); + memory->destroy(g3); memory->destroy(q1); memory->destroy(q1_mech); memory->destroy(q2); - memory->destroy(q3); + memory->destroy(q3); memory->destroy(cutsq); // to be deleted } } @@ -85,9 +85,9 @@ void PairSpinNeel::settings(int narg, char **arg) if (strcmp(update->unit_style,"metal") != 0) error->all(FLERR,"Spin simulations require metal unit style"); - + cut_spin_neel_global = force->numeric(FLERR,arg[0]); - + // reset cutoffs that have been explicitly set if (allocated) { @@ -100,7 +100,7 @@ void PairSpinNeel::settings(int narg, char **arg) } } } - + } /* ---------------------------------------------------------------------- @@ -110,30 +110,30 @@ void PairSpinNeel::settings(int narg, char **arg) void PairSpinNeel::coeff(int narg, char **arg) { if (!allocated) allocate(); - - // check if args correct + + // check if args correct if (strcmp(arg[2],"neel") != 0) error->all(FLERR,"Incorrect args in pair_style command"); if (narg != 10) 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); - + const double rij = force->numeric(FLERR,arg[3]); const double k1 = force->numeric(FLERR,arg[4]); - const double k2 = force->numeric(FLERR,arg[5]); - const double k3 = force->numeric(FLERR,arg[6]); + const double k2 = force->numeric(FLERR,arg[5]); + const double k3 = force->numeric(FLERR,arg[6]); const double l1 = force->numeric(FLERR,arg[7]); - const double l2 = force->numeric(FLERR,arg[8]); - const double l3 = force->numeric(FLERR,arg[9]); - + const double l2 = force->numeric(FLERR,arg[8]); + const double l3 = force->numeric(FLERR,arg[9]); + int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - cut_spin_neel[i][j] = rij; + cut_spin_neel[i][j] = rij; g1[i][j] = k1/hbar; q1[i][j] = l1/hbar; g1_mech[i][j] = k1; @@ -146,8 +146,8 @@ void PairSpinNeel::coeff(int narg, char **arg) count++; } } - if (count == 0) - error->all(FLERR,"Incorrect args in pair_style command"); + if (count == 0) + error->all(FLERR,"Incorrect args in pair_style command"); } @@ -160,8 +160,8 @@ void PairSpinNeel::init_style() if (!atom->sp_flag) error->all(FLERR,"Pair spin requires atom/spin style"); - // need a full neighbor list - + // need a full neighbor list + int irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; @@ -175,7 +175,7 @@ void PairSpinNeel::init_style() } 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++) { @@ -193,7 +193,7 @@ void PairSpinNeel::init_style() double PairSpinNeel::init_one(int i, int j) { - + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); return cut_spin_neel_global; @@ -214,27 +214,27 @@ void *PairSpinNeel::extract(const char *str, int &dim) void PairSpinNeel::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double evdwl,ecoul; double xi[3], rij[3], eij[3]; double spi[3], spj[3]; double fi[3], fmi[3]; double local_cut2; double rsq, inorm; - int *ilist,*jlist,*numneigh,**firstneigh; + 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 *type = atom->type; + int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -248,34 +248,34 @@ void PairSpinNeel::compute(int eflag, int vflag) itype = type[i]; jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; 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[0] = sp[i][0]; + spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + // loop on neighbors for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - spj[0] = sp[j][0]; - spj[1] = sp[j][1]; - spj[2] = sp[j][2]; + spj[0] = sp[j][0]; + spj[1] = sp[j][1]; + spj[2] = sp[j][2]; evdwl = 0.0; fi[0] = fi[1] = fi[2] = 0.0; fmi[0] = fmi[1] = fmi[2] = 0.0; rij[0] = rij[1] = rij[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]; + rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; inorm = 1.0/sqrt(rsq); eij[0] = rij[0]*inorm; eij[1] = rij[1]*inorm; @@ -294,17 +294,17 @@ void PairSpinNeel::compute(int eflag, int vflag) compute_neel_mech(i,j,rsq,eij,fi,spi,spj); } } - - f[i][0] += fi[0]; - f[i][1] += fi[1]; + + 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][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][0] -= fi[0]; + f[j][1] -= fi[1]; f[j][2] -= fi[2]; } @@ -316,15 +316,15 @@ void PairSpinNeel::compute(int eflag, int vflag) if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); } - } + } if (vflag_fdotr) virial_fdotr_compute(); - + } /* ---------------------------------------------------------------------- */ -void PairSpinNeel::compute_single_pair(int ii, double fmi[3]) +void PairSpinNeel::compute_single_pair(int ii, double fmi[3]) { int *type = atom->type; double **x = atom->x; @@ -349,13 +349,13 @@ void PairSpinNeel::compute_single_pair(int ii, double fmi[3]) spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + xi[0] = x[i][0]; xi[1] = x[i][1]; xi[2] = x[i][2]; - + eij[0] = eij[1] = eij[2] = 0.0; - + jlist = firstneigh[i]; jnum = numneigh[i]; @@ -391,7 +391,7 @@ void PairSpinNeel::compute_single_pair(int ii, double fmi[3]) void PairSpinNeel::compute_neel(int i, int j, double rsq, double eij[3], double fmi[3], double spi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; itype = type[i]; jtype = type[j]; @@ -402,8 +402,8 @@ void PairSpinNeel::compute_neel(int i, int j, double rsq, double eij[3], double double pq2x, pq2y, pq2z; // pseudo-dipolar component - - ra = rsq/g3[itype][jtype]/g3[itype][jtype]; + + ra = rsq/g3[itype][jtype]/g3[itype][jtype]; gij = 4.0*g1[itype][jtype]*ra; gij *= (1.0-g2[itype][jtype]*ra); gij *= exp(-ra); @@ -412,7 +412,7 @@ void PairSpinNeel::compute_neel(int i, int j, double rsq, double eij[3], double double scalar_eij_sj = eij[0]*spj[0] + eij[1]*spj[1] + eij[2]*spj[2]; double scalar_si_sj = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2]; - double gij_eij_sj = gij*scalar_eij_sj; + double gij_eij_sj = gij*scalar_eij_sj; double gij_3 = gij/3.0; pdx = gij_eij_sj*eij[0] - gij_3*spj[0]; pdy = gij_eij_sj*eij[1] - gij_3*spj[1]; @@ -460,11 +460,11 @@ void PairSpinNeel::compute_neel(int i, int j, double rsq, double eij[3], double void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], double fi[3], double spi[3], double spj[3]) { - int *type = atom->type; + int *type = atom->type; int itype, jtype; itype = type[i]; jtype = type[j]; - + double g_mech, gij, dgij; double q_mech, q1ij, dq1ij; double q2ij, dq2ij; @@ -480,11 +480,11 @@ void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], do double scalar_eij_sj = eij[0]*spj[0]+eij[1]*spj[1]+eij[2]*spj[2]; // pseudo-dipolar component - - g_mech = g1_mech[itype][jtype]; + + g_mech = g1_mech[itype][jtype]; ig3 = 1.0/(g3[itype][jtype]*g3[itype][jtype]); - ra = rsq*ig3; + ra = rsq*ig3; rr = drij*ig3; gij = 4.0*g_mech*ra; @@ -503,18 +503,18 @@ void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], do pdz = -(pdt1*eij[2] + pdt2*spi[2] + pdt3*spj[2]); // pseudo-quadrupolar component - - q_mech = q1_mech[itype][jtype]; + + q_mech = q1_mech[itype][jtype]; iq3 = 1.0/(q3[itype][jtype]*q3[itype][jtype]); - ra = rsq*iq3; + ra = rsq*iq3; rr = drij*iq3; q1ij = 4.0*q_mech*ra; q1ij *= (1.0-q2[itype][jtype]*ra); q1ij *= exp(-ra); q2ij = -2.0*q1ij/9.0; - + dq1ij = 1.0-ra-q2[itype][jtype]*ra*(2.0-ra); dq1ij *= 8.0*q_mech*rr*exp(-ra); dq2ij = -2.0*dq1ij/9.0; @@ -525,13 +525,13 @@ void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], do double pqt2 = scalar_eij_sj_2 - scalar_si_sj/3.0; pq1x = dq1ij * pqt1 * pqt2 * eij[0]; pq1y = dq1ij * pqt1 * pqt2 * eij[1]; - pq1z = dq1ij * pqt1 * pqt2 * eij[2]; + pq1z = dq1ij * pqt1 * pqt2 * eij[2]; double scalar_eij_si_3 = scalar_eij_si*scalar_eij_si*scalar_eij_si; double scalar_eij_sj_3 = scalar_eij_sj*scalar_eij_sj*scalar_eij_sj; double scalar_si_sj_2 = scalar_si_sj*scalar_si_sj; -/* double pqt3 = 2.0*scalar_eij_si*scalar_eij_sj_3/drij; - double pqt4 = 2.0*scalar_eij_sj*scalar_eij_si_3/drij; +/* double pqt3 = 2.0*scalar_eij_si*scalar_eij_sj_3/drij; + double pqt4 = 2.0*scalar_eij_sj*scalar_eij_si_3/drij; double pqt5 = -2.0*scalar_si_sj*scalar_eij_si/(3.0*drij); double pqt6 = -2.0*scalar_si_sj*scalar_eij_sj/(3.0*drij); // pq1x += q1ij*(spi[0]*(pqt3+pqt6) + spj[0]*(pqt4+)); @@ -600,21 +600,21 @@ void PairSpinNeel::allocate() for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; - + memory->create(cut_spin_neel,n+1,n+1,"pair/spin/soc/neel:cut_spin_neel"); memory->create(g1,n+1,n+1,"pair/spin/soc/neel:g1"); memory->create(g1_mech,n+1,n+1,"pair/spin/soc/neel:g1_mech"); - memory->create(g2,n+1,n+1,"pair/spin/soc/neel:g2"); + memory->create(g2,n+1,n+1,"pair/spin/soc/neel:g2"); memory->create(g3,n+1,n+1,"pair/spin/soc/neel:g3"); memory->create(q1,n+1,n+1,"pair/spin/soc/neel:q1"); memory->create(q1_mech,n+1,n+1,"pair/spin/soc/neel:q1_mech"); - memory->create(q2,n+1,n+1,"pair/spin/soc/neel:q2"); + memory->create(q2,n+1,n+1,"pair/spin/soc/neel:q2"); memory->create(q3,n+1,n+1,"pair/spin/soc/neel:q3"); - - memory->create(cutsq,n+1,n+1,"pair/spin/soc/neel:cutsq"); - + + memory->create(cutsq,n+1,n+1,"pair/spin/soc/neel:cutsq"); + } @@ -625,7 +625,7 @@ void PairSpinNeel::allocate() void PairSpinNeel::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++) { @@ -710,5 +710,5 @@ void PairSpinNeel::read_restart_settings(FILE *fp) } MPI_Bcast(&cut_spin_neel_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/SPIN/pair_spin_neel.h b/src/SPIN/pair_spin_neel.h index e1973d35b6..934d4a93ad 100755 --- a/src/SPIN/pair_spin_neel.h +++ b/src/SPIN/pair_spin_neel.h @@ -34,12 +34,12 @@ class PairSpinNeel : public PairSpin { double init_one(int, int); void *extract(const char *, int &); - void compute(int, int); + void compute(int, int); void compute_single_pair(int, double *); void compute_neel(int, int, double, double *, double *, double *, double *); void compute_neel_mech(int, int, double, double *, double *, double *, double *); - + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); @@ -50,7 +50,7 @@ class PairSpinNeel : public PairSpin { protected: // pseudo-dipolar and pseudo-quadrupolar coeff. - + double **g1, **g1_mech; // exchange coeffs gij double **g2, **g3; // g1 in eV, g2 adim, g3 in Ang double **q1, **q1_mech; // exchange coeffs qij diff --git a/src/atom.cpp b/src/atom.cpp index a8e0d1556d..cf4d20a71e 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -507,6 +507,7 @@ AtomVec *Atom::new_avec(const char *style, int trysuffix, int &sflag) AtomVecCreator avec_creator = (*avec_map)[style]; return avec_creator(lmp); } + error->all(FLERR,"Unknown atom style"); return NULL; }