From 4d375e72f0bfb5c51f2c809e5953a9c02292fbc9 Mon Sep 17 00:00:00 2001 From: julient31 Date: Mon, 5 Jun 2017 17:13:20 -0600 Subject: [PATCH] Changes: - DMI and ME interactions - Computation optimisations - lot of removed prints Next work: - Sequential algo implemetation - temperature simulations (check) - Work on parallelization --- in.cobalt | 59 ++++-- src/SPIN/compute_spin.cpp | 1 + src/SPIN/fix_force_spin.cpp | 6 - src/SPIN/fix_langevin_spin.cpp | 31 ++- src/SPIN/fix_nve_spin.cpp | 138 ++---------- src/SPIN/fix_nve_spin.h | 3 + src/SPIN/pair_spin.cpp | 376 +++++++++++++++++++++++++-------- src/SPIN/pair_spin.h | 38 ++-- 8 files changed, 396 insertions(+), 256 deletions(-) diff --git a/in.cobalt b/in.cobalt index 0d2975d17e..06f7aca33a 100644 --- a/in.cobalt +++ b/in.cobalt @@ -10,8 +10,8 @@ units metal dimension 3 #setting boundary conditions. (p for periodic, f for fixed, ...) -boundary p p p -#boundary f f f +#boundary p p p +boundary f f f #setting atom_style, defines what can of atoms to use in simulation (atomic, molecule, angle, dipole, ...) atom_style spin @@ -25,13 +25,25 @@ atom_modify map array ########################### #Lattice constant of fcc Cobalt -lattice fcc 3.54 -#lattice sc 2.50 -#lattice bcc 3.54 +#lattice fcc 3.54 +lattice sc 2.50 +#Test Kagome +#variable a equal sqrt(3.0) +#variable d equal 1.0/(2.0*sqrt(3.0)+1) +#variable p_1 equal $d*sqrt(3.0) +#variable p_2 equal $d*2.2*sqrt(3.0) + +#lattice custom 1.0 a1 2.0 0.0 0.0 & +# a2 1.0 $a 0.0 & +# a3 0.0 0.0 1.0 & +# basis ${p_1} ${p_2} 0.0 & +# basis ${p_2} ${p_2} 0.0 & +# basis 0.5 0.0 0.0 + #Defining a geometric region of space. Sets ID(user's choice), style(block, sphere, ...), then, args depends on the style chosen #(for block, one has x0, xf, y0, yf, z0, zf, in distance units) -region box block 0.0 2.0 0.0 2.0 0.0 2.0 +region box block 0.0 5.0 0.0 5.0 0.0 1.0 #Creating a simulation box based on the specified region. Entries: number of atom types and box ref. create_box 1 box @@ -53,13 +65,24 @@ create_atoms 1 box mass 1 1.0 #set group all mass 1.0 #Setting spins orientation and moment -#set group all spin/random 11 1.72 -set group all spin 1.72 1.0 0.0 0.0 +set group all spin/random 11 1.72 +#set group all spin 1.72 1.0 0.0 0.0 #Magnetic exchange interaction coefficient for bulk fcc Cobalt -pair_style pair/spin 4.0 -pair_coeff * * 0.0446928 0.003496 1.4885 +#pair_style pair/spin/exchange 4.0 +#type i and j | J1 | J2 | J3 +#pair_coeff * * 0.0446928 0.003496 1.4885 #pair_coeff * * 0.0 0.003496 1.4885 +#pair_style pair/spin/dmi 4.0 +# type i and j | DM (in eV) | Directions +#pair_coeff * * 0.001 0.0 0.0 1.0 +#pair_style hybrid/overlay pair/spin/exchange 4.0 pair/spin/dmi 4.0 +pair_style pair/spin 4.0 +#type i and j | interaction type | cutoff | J1 | J2 | J3 +pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 +pair_coeff * * dmi 2.6 0.01 1.0 0.0 0.0 +#pair_coeff * * me 2.6 0.01 1.0 1.0 1.0 + #Fix Langevin spins (merging damping and temperature) #Defines a cutof distance for the neighbors. @@ -67,20 +90,22 @@ pair_coeff * * 0.0446928 0.003496 1.4885 #neigh_modify every 10 check yes delay 20 neighbor 0.0 bin neigh_modify every 1 check no delay 0 + #Magnetic field fix -fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0 -#fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0 +#Type | Intensity (T or eV) | Direction +fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0 #fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0 #Fix Langevin spins (merging damping and temperature) -#fix 2 all langevin/spin 0.0 0.0 0.1 21 +#Temp | Alpha_trans | Alpha_long | Seed +fix 2 all langevin/spin 0.0 0.1 0.0 21 #fix 2 all langevin/spin 0.0 0.0 0.0 21 #Magnetic integration fix fix 3 all nve/spin -#compute total magnetization and magnetic energy -#Ite Time Mx My Mz |M| Em Tm +#compute real time, total magnetization, magnetic energy, and spin temperature +#Iteration | Time | Mx | My | Mz | |M| | Em | Tm compute mag all compute/spin fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag.dat @@ -98,6 +123,6 @@ timestep 0.00005 dump 1 all custom 50 dump_spin.lammpstrj type x y z spx spy spz #Running the simulations for N timesteps -run 300000 -#run 1 +run 60000 +#run 10 diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index bd1b09fb44..370a2848a7 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -87,6 +87,7 @@ void ComputeSpin::compute_vector() int nlocal = atom->nlocal; // compute total magnetization and magnetic energy + // compute spin temperature; See Nurdin et al., Phys. Rev. E 61, 2000 for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (atom->mumag_flag && atom->sp_flag) { diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp index e5abe3a2ca..df927858ae 100644 --- a/src/SPIN/fix_force_spin.cpp +++ b/src/SPIN/fix_force_spin.cpp @@ -185,8 +185,6 @@ void FixForceSpin::post_force(int vflag) fm[i][2] += scalar*zmag; } } - //printf("test force. 1;i=0, fx=%g, fy=%g, fz=%g, mumag=%g \n",fm[0][0],fm[0][1],fm[0][2],mumag[0]); - //printf("Field force compute, fm[0][2]=%g \n",fm[0][2]); } /* ---------------------------------------------------------------------- */ @@ -198,8 +196,6 @@ void FixForceSpin::post_force_respa(int vflag, int ilevel, int iloop) /* ---------------------------------------------------------------------- */ //No acceleration for magnetic EOM, only a "magnetic force" -//(keeping set_magneticforce in case of time--dependent mag. field implementation) - void FixForceSpin::set_magneticforce() { if (style == ZEEMAN) { @@ -212,10 +208,8 @@ void FixForceSpin::set_magneticforce() ymag = 2.0*Ka*Kay; zmag = 2.0*Ka*Kaz; } - } - /* ---------------------------------------------------------------------- potential energy in magnetic field ------------------------------------------------------------------------- */ diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index c22f224f59..f09d35c854 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -142,11 +142,10 @@ void FixLangevinSpin::post_force(int vflag) double cpx, cpy, cpz; double rx, ry, rz; - // apply transverse magnetic damping to spins - // add the damping to the effective field + // add the damping to the effective field of each spin for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - sx = sp[i][0];//Getting Mag. and Mag. force components + sx = sp[i][0]; sy = sp[i][1]; sz = sp[i][2]; @@ -164,22 +163,23 @@ void FixLangevinSpin::post_force(int vflag) fm[i][0] = fmx; fm[i][1] = fmy; - fm[i][2] = fmz; + fm[i][2] = fmz; } - //printf("test damping. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]); - //apply thermal effects - //add random field to fm + //apply thermal effects adding random fields to fm for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { + #define GAUSSIAN_R + #if defined GAUSSIAN_R rx = sigma*random->gaussian();//Drawing random distributions ry = sigma*random->gaussian(); rz = sigma*random->gaussian(); - - //rx = sigma*(random->uniform() - 0.5); - //ry = sigma*(random->uniform() - 0.5); - //rz = sigma*(random->uniform() - 0.5); - + #else + rx = sigma*(random->uniform() - 0.5); + ry = sigma*(random->uniform() - 0.5); + rz = sigma*(random->uniform() - 0.5); + #endif + fm[i][0] += rx;//Adding random field fm[i][1] += ry; fm[i][2] += rz; @@ -189,13 +189,6 @@ void FixLangevinSpin::post_force(int vflag) fm[i][2] *= Gil_factor; } - - //printf("test langevin 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]); - - //printf("test rand var: %g, sigma=%g \n",(random->uniform()-0.5),sigma); - //printf("test rand var: %g, sigma=%g \n",random->gaussian(),sigma); - //printf("test random 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]); - //printf("test dt: %g, sigma=%g \n",dts,random->gaussian(),sigma); } /* ---------------------------------------------------------------------- */ diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index 20aaf4afef..9e0f3f934c 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -89,14 +89,22 @@ 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 **sp = atom->sp; double **fm = atom->fm; + double *rmass = atom->rmass; + double *mass = atom->mass; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; int *type = atom->type; int *mask = atom->mask; + + // Advance half spins all particles + //See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011 + if (extra == SPIN) { + for (int i = 0; i < nlocal; i++){ + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + } // update half v all particles for (int i = 0; i < nlocal; i++) { @@ -109,7 +117,7 @@ void FixNVESpin::initial_integrate(int vflag) } } - // update half x for all particles + // update x for all particles for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { x[i][0] += 0.5 * dtv * v[i][0]; @@ -118,82 +126,7 @@ void FixNVESpin::initial_integrate(int vflag) } } - - size_t nbytes; - nbytes = sizeof(double) * nlocal; - int eflag = 3; - // update sp for all particles - if (extra == SPIN) { - // Advance spins - //See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011 - -#define CONC -#if defined CONC - for (int i = 0; i < nlocal; i++){ - AdvanceSingleSpin(i,dts,sp,fm); - } -#endif - -//#define SEQ -#if defined SEQ - //advance N-1 spins to half - for (int i = 0; i < nlocal-1; i++){ - //Recomp field - atom->avec->force_clear(0,nbytes); - timer->stamp(); - modify->pre_force(vflag); - timer->stamp(Timer::PAIR); - force->pair->compute(eflag,vflag); - timer->stamp(Timer::PAIR); - modify->pre_reverse(eflag,vflag); - timer->stamp(Timer::MODIFY); - comm->reverse_comm(); - timer->stamp(Timer::COMM); - modify->post_force(vflag); - - AdvanceSingleSpin(i,0.5*dts,sp,fm); - } - - //advance N spin - //Recomp field - atom->avec->force_clear(0,nbytes); - timer->stamp(); - modify->pre_force(vflag); - timer->stamp(Timer::PAIR); - force->pair->compute(eflag,vflag); - timer->stamp(Timer::PAIR); - modify->pre_reverse(eflag,vflag); - timer->stamp(Timer::MODIFY); - comm->reverse_comm(); - timer->stamp(Timer::COMM); - modify->post_force(vflag); - - AdvanceSingleSpin(nlocal-1,dts,sp,fm); - - - //advance N-1 spins to half - for (int i = nlocal-2; i >= 0; i--){ - //Recomp field - atom->avec->force_clear(0,nbytes); - timer->stamp(); - modify->pre_force(vflag); - timer->stamp(Timer::PAIR); - force->pair->compute(eflag,vflag); - timer->stamp(Timer::PAIR); - modify->pre_reverse(eflag,vflag); - timer->stamp(Timer::MODIFY); - comm->reverse_comm(); - timer->stamp(Timer::COMM); - modify->post_force(vflag); - - AdvanceSingleSpin(i,0.5*dts,sp,fm); - } - -#endif - - } - -#define FORCE_PRINT +//#define FORCE_PRINT #if defined FORCE_PRINT FILE* file_force=NULL; file_force=fopen("spin_force_Lammps.dat","a"); @@ -206,8 +139,6 @@ void FixNVESpin::initial_integrate(int vflag) } - - /* ---------------------------------------------------------------------- */ void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm) @@ -230,30 +161,6 @@ void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm) cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2]; cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0]; -//#define ALG_MA //Ma algo -#if defined ALG_MA - double A[3]; - A[0]= A[1] = A[2] = 0.0; - double xi=fmsq*dts; - double zeta=0.0; // this is because omega contains already the transverse damping - double chi=energy/fmsq; - double expo=exp(2.0*zeta); - double K1=1.0+expo+chi*(1.0-expo); - double K=2.0*exp(zeta)/K1; - double Ktrigo=1.0+0.25*xi*xi; - double cosinus=(1.0-0.25*xi*xi)/Ktrigo; //cos(xi) - double sinus=xi/Ktrigo; //sin(xi) - A[0]=K*cosinus; - A[1]=K*sinus; - A[2]=(1.0-expo+chi*(1.0+expo-2.0*exp(zeta)*cosinus))/K1; - - g[0] = A[0]*sp[i][0]+A[1]*cp[0]/fmsq+A[2]*fm[i][0]/fmsq; - g[1] = A[0]*sp[i][1]+A[1]*cp[1]/fmsq+A[2]*fm[i][0]/fmsq; - g[2] = A[0]*sp[i][2]+A[1]*cp[2]/fmsq+A[2]*fm[i][0]/fmsq; -#endif - -#define ALG_OM //Omelyan algo -#if defined ALG_OM g[0] = sp[i][0]+cp[0]*dts; g[1] = sp[i][1]+cp[1]*dts; g[2] = sp[i][2]+cp[2]*dts; @@ -265,7 +172,6 @@ void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm) g[0] /= (1+0.25*fm2*dts*dts); g[1] /= (1+0.25*fm2*dts*dts); g[2] /= (1+0.25*fm2*dts*dts); -#endif sp[i][0] = g[0]; sp[i][1] = g[1]; @@ -277,7 +183,6 @@ void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm) sp[i][0] *= scale; sp[i][1] *= scale; sp[i][2] *= scale; - } /* ---------------------------------------------------------------------- */ @@ -290,6 +195,8 @@ void FixNVESpin::final_integrate() double **x = atom->x; double **v = atom->v; double **f = atom->f; + double **sp = atom->sp; + double **fm = atom->fm; double *rmass = atom->rmass; double *mass = atom->mass; int nlocal = atom->nlocal; @@ -297,15 +204,6 @@ void FixNVESpin::final_integrate() int *type = atom->type; int *mask = atom->mask; - // update half x for all particles - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += 0.5 * dtv * v[i][0]; - x[i][1] += 0.5 * dtv * v[i][1]; - x[i][2] += 0.5 * dtv * v[i][2]; - } - } - // update half v for all particles for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -316,4 +214,12 @@ void FixNVESpin::final_integrate() v[i][2] += dtfm * f[i][2]; } } + + // Advance half spins all particles + //See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011 + if (extra == SPIN) { + for (int i = 0; i < nlocal; i++){ + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + } } diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 82e9822bf1..31ef30de3e 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -35,6 +35,9 @@ class FixNVESpin : public FixNVE { void AdvanceSingleSpin(int, double, double **, double **); virtual void final_integrate(); +//Sorting atoms/spins routine +void SortSpins(); + protected: int extra; double dts; diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index b87c14b7b9..f09b4d05fd 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -33,6 +33,8 @@ using namespace MathConst; PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; + exch_flag = 0; + dmi_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -43,10 +45,25 @@ PairSpin::~PairSpin() memory->destroy(setflag); memory->destroy(cut_spin_exchange); - memory->destroy(cut_spin_dipolar); memory->destroy(J_1); memory->destroy(J_2); - memory->destroy(J_2); + memory->destroy(J_2); + + memory->destroy(cut_spin_dmi); + memory->destroy(DM); + memory->destroy(v_dmx); + memory->destroy(v_dmy); + memory->destroy(v_dmz); + + memory->destroy(cut_spin_me); + memory->destroy(ME); + memory->destroy(v_mex); + memory->destroy(v_mey); + memory->destroy(v_mez); + + memory->destroy(fmi); + memory->destroy(fmj); + memory->destroy(cutsq); } } @@ -57,8 +74,9 @@ void PairSpin::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double evdwl,ecoul; - double xtmp,ytmp,ztmp,fmix,fmiy,fmiz,fmjx,fmjy,fmjz,omx,omy,omz; - double cut, Jex, ra; + double xtmp,ytmp,ztmp; + double fmix,fmiy,fmiz,fmjx,fmjy,fmjz; + double cut_ex,cut_dmi,cut_me; double rsq,rd,delx,dely,delz; int *ilist,*jlist,*numneigh,**firstneigh; @@ -90,47 +108,52 @@ void PairSpin::compute(int eflag, int vflag) jlist = firstneigh[i]; jnum = numneigh[i]; - //Exchange interaction + //Loop on Neighbors for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - fmix = fmiy = fmiz = 0.0; - fmjx = fmjy = fmjz = 0.0; + fmi[0] = fmi[1] = fmi[2] = 0.0; + fmj[0] = fmj[1] = fmj[2] = 0.0; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; //square or inter-atomic distance rd = sqrt(rsq); //Inter-atomic distance - cut = cut_spin_exchange_global; - - if (rd <= cut) { - itype = type[i]; - jtype = type[j]; + itype = type[i]; + jtype = type[j]; - ra = (rd/J_3[itype][jtype])*(rd/J_3[itype][jtype]); - Jex = 4.0*J_1[itype][jtype]*ra; - Jex *= (1.0-J_2[itype][jtype]*ra); - Jex *= exp(-ra); - - fmix = Jex*sp[j][0]; - fmiy = Jex*sp[j][1]; - fmiz = Jex*sp[j][2]; - - fmjx = Jex*sp[i][0]; - fmjy = Jex*sp[i][1]; - fmjz = Jex*sp[i][2]; - } - - fm[i][0] += fmix; - fm[i][1] += fmiy; - fm[i][2] += fmiz; + //Exchange interaction + if (exch_flag) { + cut_ex = cut_spin_exchange[itype][jtype]; + if (rd <= cut_ex) { + compute_exchange(i,j,rsq,fmi,fmj); + } + } + //DM interaction + if (dmi_flag){ + cut_dmi = cut_spin_dmi[itype][jtype]; + if (rd <= cut_dmi){ + compute_dmi(i,j,fmi,fmj); + } + } + //ME interaction + if (me_flag){ + cut_me = cut_spin_me[itype][jtype]; + if (rd <= cut_me){ + compute_me(i,j,fmi,fmj); + } + } + + fm[i][0] += fmi[0]; + fm[i][1] += fmi[1]; + fm[i][2] += fmi[2]; if (newton_pair || j < nlocal) { - fm[j][0] += fmjx; - fm[j][1] += fmjy; - fm[j][2] += fmjz; + fm[j][0] += fmj[0]; + fm[j][1] += fmj[1]; + fm[j][2] += fmj[2]; } } } @@ -138,6 +161,95 @@ void PairSpin::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } +/* ---------------------------------------------------------------------- */ +inline void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj) +{ + int *type = atom->type; + int itype, jtype; + double **sp = atom->sp; + double dmix,dmiy,dmiz; + double Jex, ra; + itype = type[i]; + jtype = type[j]; + + ra = rsq/J_3[itype][jtype]/J_3[itype][jtype]; + Jex = 4.0*J_1[itype][jtype]*ra; + Jex *= (1.0-J_2[itype][jtype]*ra); + Jex *= exp(-ra); + + fmi[0] += Jex*sp[j][0]; + fmi[1] += Jex*sp[j][1]; + fmi[2] += Jex*sp[j][2]; + + fmj[0] += Jex*sp[i][0]; + fmj[1] += Jex*sp[i][1]; + fmj[2] += Jex*sp[i][2]; + +} + +/* ---------------------------------------------------------------------- */ +inline void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj) +{ + + int *type = atom->type; + int itype, jtype; + double **sp = atom->sp; + double dmix,dmiy,dmiz; + itype = type[i]; + jtype = type[j]; + + dmix = DM[itype][jtype]*v_dmx[itype][jtype]; + dmiy = DM[itype][jtype]*v_dmy[itype][jtype]; + dmiz = DM[itype][jtype]*v_dmz[itype][jtype]; + + fmi[0] += sp[j][1]*dmiz-sp[j][2]*dmiy; + fmi[1] += sp[j][2]*dmix-sp[j][0]*dmiz; + fmi[2] += sp[j][0]*dmiy-sp[j][1]*dmix; + + fmj[0] -= sp[i][1]*dmiz-sp[i][2]*dmiy; + fmj[1] -= sp[i][2]*dmix-sp[i][0]*dmiz; + fmj[2] -= sp[i][0]*dmiy-sp[i][1]*dmix; +} + +/* ---------------------------------------------------------------------- */ +inline void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) +{ + int *type = atom->type; + int itype, jtype; + itype = type[i]; + jtype = type[j]; + double **sp = atom->sp; + double **x = atom->x; + double meix,meiy,meiz; + double rx, ry, rz, inorm; + + rx = x[j][0] - x[i][0]; + ry = x[j][1] - x[i][1]; + rz = x[j][2] - x[i][2]; + inorm = 1.0/sqrt(rx*rx+ry*ry+rz*rz); + rx *= inorm; + ry *= inorm; + rz *= inorm; + + meix = v_mey[itype][jtype]*rz - v_mez[itype][jtype]*ry; + meiy = v_mez[itype][jtype]*rx - v_mex[itype][jtype]*rz; + meiz = v_mex[itype][jtype]*ry - v_mey[itype][jtype]*rx; + + meix *= ME[itype][jtype]; + meiy *= ME[itype][jtype]; + meiz *= ME[itype][jtype]; + + fmi[0] += sp[j][1]*meiz - sp[j][2]*meiy; + fmi[1] += sp[j][2]*meix - sp[j][0]*meiz; + fmi[2] += sp[j][0]*meiy - sp[j][1]*meix; + + fmj[0] -= sp[i][1]*meiz - sp[i][2]*meiy; + fmj[1] -= sp[i][2]*meix - sp[i][0]*meiz; + fmj[2] -= sp[i][0]*meiy - sp[i][1]*meix; + + // printf("test val fmi=%g, fmj=%g \n",fmi[2],fmj[2]); + +} /* ---------------------------------------------------------------------- allocate all arrays @@ -154,11 +266,25 @@ void PairSpin::allocate() setflag[i][j] = 0; memory->create(cut_spin_exchange,n+1,n+1,"pair:cut_spin_exchange"); - memory->create(cut_spin_dipolar,n+1,n+1,"pair:cut_spin_dipolar"); memory->create(J_1,n+1,n+1,"pair:J_1"); memory->create(J_2,n+1,n+1,"pair:J_2"); memory->create(J_3,n+1,n+1,"pair:J_3"); - + + 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(cut_spin_me,n+1,n+1,"pair:cut_spin_me"); + memory->create(ME,n+1,n+1,"pair:ME"); + memory->create(v_mex,n+1,n+1,"pair:ME_vector_x"); + memory->create(v_mey,n+1,n+1,"pair:ME_vector_y"); + memory->create(v_mez,n+1,n+1,"pair:ME_vector_z"); + + memory->create(fmi,3,"pair:fmi"); + memory->create(fmj,3,"pair:fmj"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); } @@ -172,14 +298,11 @@ void PairSpin::settings(int narg, char **arg) if (narg < 1 || narg > 2) error->all(FLERR,"Incorrect number of args in pair_style pair/spin command"); - if (strcmp(update->unit_style,"electron") == 0) - error->all(FLERR,"Cannot (yet) use 'electron' units with spins"); + 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]); + cut_spin_pair_global = force->numeric(FLERR,arg[0]); - if (narg == 1) cut_spin_dipolar_global = cut_spin_exchange_global; - else cut_spin_dipolar_global = force->numeric(FLERR,arg[1]); - // reset cutoffs that have been explicitly set if (allocated) { @@ -187,8 +310,8 @@ void PairSpin::settings(int narg, char **arg) for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) { - cut_spin_exchange[i][j] = cut_spin_exchange_global; - cut_spin_dipolar[i][j] = cut_spin_dipolar_global; + cut_spin_exchange[i][j] = cut_spin_pair_global; + cut_spin_dmi[i][j] = cut_spin_pair_global; } } @@ -201,38 +324,108 @@ void PairSpin::settings(int narg, char **arg) void PairSpin::coeff(int narg, char **arg) { - if (narg != 5) - error->all(FLERR,"Incorrect number of args for pair spin coefficients"); - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - double J1 = force->numeric(FLERR,arg[2]); - double J2 = force->numeric(FLERR,arg[3]); - double J3 = force->numeric(FLERR,arg[4]); - - double hbar = force->hplanck/MY_2PI; - J1 /= hbar; - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - J_1[i][j] = J1; - J_2[i][j] = J2; - J_3[i][j] = J3; - setflag[i][j] = 1; - count++; - } - } - if (count == 0) error->all(FLERR,"Incorrect args for spinpair coefficients"); - - //Simple (Anti)Ferromagnetic exchange for now. - //Check if Jex [][] still works for Ferrimagnetic exchange - -} + if (!allocated) allocate(); + if (strcmp(arg[2],"exchange")==0){ + if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command"); + exch_flag = 1; + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double rij = force->numeric(FLERR,arg[3]); + double J1 = force->numeric(FLERR,arg[4]); + double J2 = force->numeric(FLERR,arg[5]); + double J3 = force->numeric(FLERR,arg[6]); + + double hbar = force->hplanck/MY_2PI; + J1 /= hbar; + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_spin_exchange[i][j] = rij; + J_1[i][j] = J1; + J_2[i][j] = J2; + J_3[i][j] = J3; + setflag[i][j] = 1; + count++; + } + } + if (count == 0) error->all(FLERR,"Incorrect args in pair_style command"); + } else if (strcmp(arg[2],"dmi")==0) { + if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command"); + dmi_flag = 1; + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double rij = force->numeric(FLERR,arg[3]); + double dm = force->numeric(FLERR,arg[4]); + 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; + + double hbar = force->hplanck/MY_2PI; + dm /= hbar; + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_spin_dmi[i][j] = rij; + DM[i][j] = dm; + v_dmx[i][j] = dmx; + v_dmy[i][j] = dmy; + v_dmz[i][j] = dmz; + setflag[i][j] = 1; + count++; + } + } + if (count == 0) error->all(FLERR,"Incorrect args in pair_style command"); + } else if (strcmp(arg[2],"me")==0) { + if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command"); + me_flag = 1; + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double rij = force->numeric(FLERR,arg[3]); + double me = 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 inorm = 1.0/(mex*mex+mey*mey+mez*mez); + mex *= inorm; + mey *= inorm; + mez *= inorm; + + double hbar = force->hplanck/MY_2PI; + me /= hbar; + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_spin_me[i][j] = rij; + DM[i][j] = me; + v_mex[i][j] = mex; + v_mey[i][j] = mey; + v_mez[i][j] = mez; + setflag[i][j] = 1; + count++; + } + } + if (count == 0) error->all(FLERR,"Incorrect args in pair_style command"); + } else error->all(FLERR,"Incorrect args in pair_style command"); + + //Check if Jex [][] still works for Ferrimagnetic exchange +} /* ---------------------------------------------------------------------- @@ -256,7 +449,7 @@ double PairSpin::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; + return cut_spin_pair_global; } /* ---------------------------------------------------------------------- @@ -272,10 +465,26 @@ void PairSpin::write_restart(FILE *fp) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { - fwrite(&J_1[i][j],sizeof(double),1,fp); - fwrite(&J_2[i][j],sizeof(double),1,fp); - fwrite(&J_3[i][j],sizeof(double),1,fp); - fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp); + if (exch_flag){ + fwrite(&J_1[i][j],sizeof(double),1,fp); + fwrite(&J_2[i][j],sizeof(double),1,fp); + fwrite(&J_3[i][j],sizeof(double),1,fp); + fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp); + } + if (dmi_flag) { + fwrite(&DM[i][j],sizeof(double),1,fp); + fwrite(&v_dmx[i][j],sizeof(double),1,fp); + fwrite(&v_dmy[i][j],sizeof(double),1,fp); + fwrite(&v_dmz[i][j],sizeof(double),1,fp); + fwrite(&cut_spin_dmi[i][j],sizeof(double),1,fp); + } + if (me_flag) { + fwrite(&ME[i][j],sizeof(double),1,fp); + fwrite(&v_mex[i][j],sizeof(double),1,fp); + fwrite(&v_mey[i][j],sizeof(double),1,fp); + fwrite(&v_mez[i][j],sizeof(double),1,fp); + fwrite(&cut_spin_me[i][j],sizeof(double),1,fp); + } } } } @@ -310,16 +519,15 @@ void PairSpin::read_restart(FILE *fp) } } } - + /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairSpin::write_restart_settings(FILE *fp) { - fwrite(&cut_spin_exchange_global,sizeof(double),1,fp); - fwrite(&cut_spin_dipolar_global,sizeof(double),1,fp); + fwrite(&cut_spin_pair_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } @@ -331,13 +539,11 @@ void PairSpin::write_restart_settings(FILE *fp) void PairSpin::read_restart_settings(FILE *fp) { if (comm->me == 0) { - fread(&cut_spin_exchange_global,sizeof(double),1,fp); - fread(&cut_spin_dipolar_global,sizeof(double),1,fp); + fread(&cut_spin_pair_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } - MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_spin_dipolar_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_spin_pair_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index fd5822c163..8b3e174de6 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -33,24 +33,36 @@ class PairSpin : public Pair { void coeff(int, char **); void init_style(); double init_one(int, int); + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - - //Test transf. force -//#define TRANS_FORCE -#if defined TRANS_FORCE - void transferfm(double **); -#endif + inline void compute_exchange(int, int, double, double *, double *); + inline void compute_dmi(int, int, double *, double *); + inline void compute_me(int, int, double *, double *); + protected: - double cut_spin_exchange_global, cut_spin_dipolar_global; //Global cutting distance - double **cut_spin_exchange; //cutting distance for each exchange interaction - double **cut_spin_dipolar; //cutting distance for the dipolar interaction + int exch_flag,dmi_flag,me_flag; + double cut_spin_pair_global; + double cut_spin_dipolar_global; - double **J_1, **J_2, **J_3; //coefficients for computing the exchange interaction Jij - //J1 is an energy (in eV), J2 is adim and J3 is a distance (in Ang) + double **cut_spin_exchange; //cutting distance exchange + double **cut_spin_dmi; //cutting distance dmi + double **cut_spin_me; //cutting distance me + + double **J_1, **J_2, **J_3; //exchange coeffs Jij + //J1 in eV, J2 adim and J3 in Ang + double **DM; + double **v_dmx, **v_dmy, **v_dmz;//DMI coeffs + //DM int. in eV, v direction + + double **ME; + double **v_mex, **v_mey, **v_mez;//ME coeffs + //ME in eV, v direction + + double *fmi, *fmj; //Temp var. in compute void allocate(); }; @@ -66,9 +78,9 @@ E: Incorrect args in pair_style command Self-explanatory. -E: Cannot (yet) use 'electron' units with spins +E: Spin simulations require metal unit style -This feature is not yet supported. +Self-explanatory. E: Incorrect args for pair coefficients