diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index 17b9d3eb22..873a21ea66 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -67,6 +67,9 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : H_field = 0.0; nhx = nhy = nhz = 0.0; hx = hy = hz = 0.0; + stt_field = 0.0; + nsttx = nstty = nsttz = 0.0; + sttx = stty = sttz = 0.0; Ka = 0.0; nax = nay = naz = 0.0; Kax = Kay = Kaz = 0.0; @@ -74,8 +77,11 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nc1x = nc1y = nc1z = 0.0; nc2x = nc2y = nc2z = 0.0; nc3x = nc3y = nc3z = 0.0; + K6 = 0.0; + n6x = n6y = n6z = 0.0; + m6x = m6y = m6z = 0.0; - zeeman_flag = aniso_flag = cubic_flag = 0; + zeeman_flag = stt_flag = aniso_flag = cubic_flag = hexaniso_flag = 0; int iarg = 3; while (iarg < narg) { @@ -87,6 +93,14 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nhy = utils::numeric(FLERR,arg[iarg+3],false,lmp); nhz = utils::numeric(FLERR,arg[iarg+4],false,lmp); iarg += 5; + } else if (strcmp(arg[iarg],"stt") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command"); + stt_flag = 1; + stt_field = utils::numeric(FLERR,arg[iarg+1],false,lmp); + nsttx = utils::numeric(FLERR,arg[iarg+2],false,lmp); + nstty = utils::numeric(FLERR,arg[iarg+3],false,lmp); + nsttz = utils::numeric(FLERR,arg[iarg+4],false,lmp); + iarg += 5; } else if (strcmp(arg[iarg],"anisotropy") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command"); aniso_flag = 1; @@ -110,6 +124,17 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nc3y = utils::numeric(FLERR,arg[iarg+10],false,lmp); nc3z = utils::numeric(FLERR,arg[iarg+11],false,lmp); iarg += 12; + } else if (strcmp(arg[iarg],"hexaniso") == 0) { + if (iarg+7 > narg) error->all(FLERR,"Illegal fix precession/spin command"); + hexaniso_flag = 1; + K6 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + n6x = utils::numeric(FLERR,arg[iarg+2],false,lmp); + n6y = utils::numeric(FLERR,arg[iarg+3],false,lmp); + n6z = utils::numeric(FLERR,arg[iarg+4],false,lmp); + m6x = utils::numeric(FLERR,arg[iarg+5],false,lmp); + m6y = utils::numeric(FLERR,arg[iarg+6],false,lmp); + m6z = utils::numeric(FLERR,arg[iarg+7],false,lmp); + iarg += 8; } else error->all(FLERR,"Illegal precession/spin command"); } @@ -123,6 +148,13 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nhz *= inorm; } + if (stt_flag) { + inorm = 1.0/sqrt(nsttx*nsttx + nstty*nstty + nsttz*nsttz); + nsttx *= inorm; + nstty *= inorm; + nsttz *= inorm; + } + if (aniso_flag) { inorm = 1.0/sqrt(nax*nax + nay*nay + naz*naz); nax *= inorm; @@ -144,6 +176,27 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nc3y *= inorm; nc3z *= inorm; } + + if (hexaniso_flag) { + inorm = 1.0/sqrt(n6x*n6x + n6y*n6y + n6z*n6z); + n6x *= inorm; + n6y *= inorm; + n6z *= inorm; + inorm = 1.0/sqrt(m6x*m6x + m6y*m6y + m6z*m6z); + m6x *= inorm; + m6y *= inorm; + m6z *= inorm; + l6x = (n6z*m6y-n6y*m6z); + l6y = (n6x*m6z-n6z*m6x); + l6z = (n6y*m6x-n6x*m6y); + inorm = 1.0/sqrt(l6x*l6x + l6y*l6y + l6z*l6z); + l6x *= inorm; + l6y *= inorm; + l6z *= inorm; + m6x = (l6z*n6y-l6y*n6z); + m6y = (l6x*n6z-l6z*n6x); + m6z = (l6y*n6x-l6x*n6y); + } degree2rad = MY_PI/180.0; time_origin = update->ntimestep; @@ -185,8 +238,9 @@ void FixPrecessionSpin::init() Kah = Ka/hbar; k1ch = k1c/hbar; k2ch = k2c/hbar; + K6h = K6/hbar; - if (utils::strmatch(update->integrate_style,"^respa")) { + if (strstr(update->integrate_style,"respa")) { ilevel_respa = ((Respa *) update->integrate)->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); } @@ -201,12 +255,11 @@ void FixPrecessionSpin::init() // check that fix precession/spin is only declared once - int iprec = 0; - for (int iforce = 0; iforce < modify->nfix; iforce++) - if (strstr(modify->fix[iforce]->style,"precession/spin")) iprec++; - if (iprec > 1) - error->all(FLERR,"precession/spin command can only be declared once"); - + // int iprec = 0; + // for (int iforce = 0; iforce < modify->nfix; iforce++) + // if (strstr(modify->fix[iforce]->style,"precession/spin")) iprec++; + // if (iprec > 1) + // error->all(FLERR,"precession/spin command can only be declared once"); varflag = CONSTANT; if (magfieldstyle != CONSTANT) varflag = EQUAL; @@ -225,7 +278,7 @@ void FixPrecessionSpin::init() void FixPrecessionSpin::setup(int vflag) { - if (utils::strmatch(update->integrate_style,"^verlet")) + if (strstr(update->integrate_style,"verlet")) post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa); @@ -283,16 +336,26 @@ void FixPrecessionSpin::post_force(int /* vflag */) epreci -= compute_zeeman_energy(spi); } + if (stt_flag) { // compute Spin Transfer Torque + compute_stt(spi,fmi); + epreci -= compute_stt_energy(spi); + } + if (aniso_flag) { // compute magnetic anisotropy compute_anisotropy(spi,fmi); epreci -= compute_anisotropy_energy(spi); } - if (cubic_flag) { // compute cubic anisotropy + if (cubic_flag) { // compute cubic anisotropy compute_cubic(spi,fmi); epreci -= compute_cubic_energy(spi); } + if (hexaniso_flag) { // compute hexagonal anisotropy + compute_hexaniso(spi,fmi); + epreci -= compute_hexaniso_energy(spi); + } + emag[i] += epreci; eprec += epreci; fm[i][0] += fmi[0]; @@ -309,12 +372,16 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f int *mask = atom->mask; if (mask[i] & groupbit) { if (zeeman_flag) compute_zeeman(i,fmi); + if (stt_flag) compute_stt(spi,fmi); if (aniso_flag) compute_anisotropy(spi,fmi); if (cubic_flag) compute_cubic(spi,fmi); + if (hexaniso_flag) compute_hexaniso(spi,fmi); } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Zeeman +------------------------------------------------------------------------- */ void FixPrecessionSpin::compute_zeeman(int i, double fmi[3]) { @@ -334,14 +401,39 @@ double FixPrecessionSpin::compute_zeeman_energy(double spi[4]) return energy; } +/* ---------------------------------------------------------------------- + STT +------------------------------------------------------------------------- */ + +void FixPrecessionSpin::compute_stt(double spi[3], double fmi[3]) +{ + double sx = spi[0]; + double sy = spi[1]; + double sz = spi[2]; + fmi[0] += 1.0*stt_field*( sy*nsttz-sz*nstty); + fmi[1] += 1.0*stt_field*(-sx*nsttz+sz*nsttx); + fmi[2] += 1.0*stt_field*( sx*nstty-sy*nsttx); +} + /* ---------------------------------------------------------------------- */ +double FixPrecessionSpin::compute_stt_energy(double spi[3]) +{ + double energy = 0.0; // Non-conservative force + return energy; +} + +/* ---------------------------------------------------------------------- + compute uniaxial anisotropy interaction for spin i +------------------------------------------------------------------------- */ + void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3]) { double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2]; fmi[0] += scalar*Kax; fmi[1] += scalar*Kay; fmi[2] += scalar*Kaz; + // printf("fm pres1: %g %g %g \n",fmi[0],fmi[1],fmi[2]); } /* ---------------------------------------------------------------------- */ @@ -393,9 +485,7 @@ void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3]) fmi[2] += (fourz + sixz); } -/* ---------------------------------------------------------------------- - compute cubic aniso energy of spin i -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ double FixPrecessionSpin::compute_cubic_energy(double spi[3]) { @@ -412,6 +502,62 @@ double FixPrecessionSpin::compute_cubic_energy(double spi[3]) return energy; } +/* ---------------------------------------------------------------------- + compute hexagonal anisotropy interaction for spin i +------------------------------------------------------------------------- */ + +void FixPrecessionSpin::compute_hexaniso(double spi[3], double fmi[3]) +{ + double s_x,s_y,s_z; + double pf, phi, ssint2; + + // changing to the axes' frame + + s_x = l6x*spi[0]+l6y*spi[1]+l6z*spi[2]; + s_y = m6x*spi[0]+m6y*spi[1]+m6z*spi[2]; + s_z = n6x*spi[0]+n6y*spi[1]+n6z*spi[2]; + + // hexagonal anisotropy in the axes' frame + + phi = atan2(s_y,s_x); + ssint2 = s_x*s_x + s_y*s_y; // s^2sin^2(theta) + pf = 6.0 * K6h * ssint2*ssint2*sqrt(ssint2); // 6*K_6*s^5*sin^5(theta) + double fm_x = pf*cos(5*phi); + double fm_y = -pf*sin(5*phi); + double fm_z = 0; + + // back to the lab's frame + + fmi[0] += fm_x*l6x+fm_y*m6x+fm_z*n6x; + fmi[1] += fm_x*l6y+fm_y*m6y+fm_z*n6y; + fmi[2] += fm_x*l6z+fm_y*m6z+fm_z*n6z; +} + +/* ---------------------------------------------------------------------- + compute hexagonal aniso energy of spin i +------------------------------------------------------------------------- */ + +double FixPrecessionSpin::compute_hexaniso_energy(double spi[3]) +{ + double energy = 0.0; + double s_x,s_y,s_z, phi,ssint2; + + // changing to the axes' frame + + s_x = l6x*spi[0]+l6y*spi[1]+l6z*spi[2]; + s_y = m6x*spi[0]+m6y*spi[1]+m6z*spi[2]; + s_z = n6x*spi[0]+n6y*spi[1]+n6z*spi[2]; + + // hexagonal anisotropy in the axes' frame + + phi = atan2(s_y,s_z); + ssint2 = s_x*s_x + s_y*s_y; + + energy = K6 * ssint2*ssint2*ssint2*cos(6*phi); + + return 2.0*energy; +} + /* ---------------------------------------------------------------------- */ void FixPrecessionSpin::set_magneticprecession() @@ -421,6 +567,13 @@ void FixPrecessionSpin::set_magneticprecession() hy = H_field*nhy; hz = H_field*nhz; } + + if (stt_flag) { + sttx = stt_field*nsttx; + stty = stt_field*nstty; + sttz = stt_field*nsttz; + } + if (aniso_flag) { Kax = 2.0*Kah*nax; Kay = 2.0*Kah*nay; diff --git a/src/SPIN/fix_precession_spin.h b/src/SPIN/fix_precession_spin.h index 9c3c616077..d06a04f8f1 100644 --- a/src/SPIN/fix_precession_spin.h +++ b/src/SPIN/fix_precession_spin.h @@ -39,7 +39,7 @@ class FixPrecessionSpin : public Fix { void min_post_force(int); double compute_scalar(); - int zeeman_flag, aniso_flag, cubic_flag; + int zeeman_flag, stt_flag, aniso_flag, cubic_flag, hexaniso_flag; void compute_single_precession(int, double *, double *); // zeeman calculations @@ -47,6 +47,11 @@ class FixPrecessionSpin : public Fix { void compute_zeeman(int, double *); double compute_zeeman_energy(double *); + // stt calculations + + void compute_stt(double *, double *); + double compute_stt_energy(double *); + // uniaxial aniso calculations void compute_anisotropy(double *, double *); @@ -57,6 +62,11 @@ class FixPrecessionSpin : public Fix { void compute_cubic(double *, double *); double compute_cubic_energy(double *); + // hexagonal aniso calculations + + void compute_hexaniso(double *, double *); + double compute_hexaniso_energy(double *); + // storing magnetic energies int nlocal_max; // max nlocal (for list size) @@ -83,6 +93,12 @@ class FixPrecessionSpin : public Fix { double nhx, nhy, nhz; double hx, hy, hz; // temp. force variables + // STT intensity and direction + + double stt_field; + double nsttx, nstty, nsttz; + double sttx, stty, sttz; + // magnetic anisotropy intensity and direction double Ka; // aniso const. in eV @@ -98,6 +114,13 @@ class FixPrecessionSpin : public Fix { double nc2x,nc2y,nc2z; double nc3x,nc3y,nc3z; + // hexagonal anisotropy + double K6; // hexagonal aniso const. in eV + double K6h; // hexagonal aniso const. in rad.THz + double n6x,n6y,n6z; // main axis + double m6x,m6y,m6z; // secondary (perpendicular) axis + double l6x,l6y,l6z; // =(m x n) + void set_magneticprecession(); };