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

whitespace cleanup

Browse Source
This commit is contained in:
Axel Kohlmeyer
2018-06-18 17:54:04 -04:00
parent 83ae0ad26f
commit da5931d65d
20 changed files with 450 additions and 449 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
src/SPIN/atom_vec_spin.cpp
View File

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

4
src/SPIN/atom_vec_spin.h
View File

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

42
src/SPIN/compute_spin.cpp
View File

@ -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;
}
/* ----------------------------------------------------------------------

34
src/SPIN/fix_langevin_spin.cpp
View File

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

2
src/SPIN/fix_langevin_spin.h
View File

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

116
src/SPIN/fix_nve_spin.cpp
View File

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

26
src/SPIN/fix_nve_spin.h
View File

@ -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.
*/

46
src/SPIN/fix_precession_spin.cpp
View File

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

22
src/SPIN/fix_precession_spin.h
View File

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

6
src/SPIN/pair_spin.cpp
View File

@ -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.
------------------------------------------------------------------------- */

2
src/SPIN/pair_spin.h
View File

@ -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() {}
};

112
src/SPIN/pair_spin_dmi.cpp
View File

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

6
src/SPIN/pair_spin_dmi.h
View File

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

118
src/SPIN/pair_spin_exchange.cpp
View File

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

6
src/SPIN/pair_spin_exchange.h
View File

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

148
src/SPIN/pair_spin_magelec.cpp
View File

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

10
src/SPIN/pair_spin_magelec.h
View File

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

152
src/SPIN/pair_spin_neel.cpp
View File

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

6
src/SPIN/pair_spin_neel.h
View File

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

1
src/atom.cpp
View File

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