diff --git a/in.cobalt b/in.cobalt index 06f7aca33a..8a9de336a7 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 @@ -20,30 +20,35 @@ atom_style spin atom_modify map array +#newton off for pair spin in SEQNEI +#newton off off + ########################### #######Create atoms######## ########################### #Lattice constant of fcc Cobalt -#lattice fcc 3.54 -lattice sc 2.50 +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) +#variable a equal sqrt(3.0)/8.0 +#variable b equal 3.0*sqrt(3.0)/8.0 +#variable c equal sqrt(3.0)/4.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 +#lattice custom 2.5 a1 1.0 0.0 0.0 & +# a2 0.0 1.0 0.0 & +# a3 0.0 0.0 1.0 & +# basis 0.0 $a 0.0 & +# basis 0.25 $a 0.0 & +# basis 0.375 0.0 0.0 & +# basis 0.25 $b 0.0 & +# basis 0.5 $b 0.0 & +# basis 0.625 $c 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 5.0 0.0 5.0 0.0 1.0 +region box block 0.0 2.0 0.0 2.0 0.0 2.0 #Creating a simulation box based on the specified region. Entries: number of atom types and box ref. create_box 1 box @@ -65,27 +70,18 @@ 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/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 +#type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange) +pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 +#type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME) +#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. +#Define a skin distance, update neigh list every #neighbor 1.0 bin #neigh_modify every 10 check yes delay 20 neighbor 0.0 bin @@ -93,13 +89,15 @@ neigh_modify every 1 check no delay 0 #Magnetic field fix #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 zeeman 10.0 0.0 0.0 1.0 #fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0 +#fix 1 all force/spin anisotropy -0.1 0.0 0.0 1.0 +#fix 1 all force/spin anisotropy 0.1 0.0 1.0 0.0 #Fix Langevin spins (merging damping and temperature) #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 +#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 @@ -107,13 +105,13 @@ fix 3 all nve/spin #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 +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_seqnei_test.dat #Defining a computation that will be performed on a group of atoms. #Entries: ID(user assigned), group-ID(group of atoms to peform the sim on), style(temp, pe, ...), args #Setting the timestep for the simulation -timestep 0.00005 +timestep 0.00001 ################## #######run######## @@ -123,6 +121,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 60000 -#run 10 +#run 200000 +run 1 diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index 3807ac0e56..d31f4691d1 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -46,6 +46,7 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) forceclearflag = 1; atom->mumag_flag = atom->sp_flag = 1; + } diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp index df927858ae..27408456d6 100644 --- a/src/SPIN/fix_force_spin.cpp +++ b/src/SPIN/fix_force_spin.cpp @@ -26,6 +26,7 @@ #include "math_const.h" #include "error.h" #include "force.h" +#include "memory.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -59,10 +60,13 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a Hx = Hy = Hz = 0.0; Ka = 0.0; Kax = Kay = Kaz = 0.0; - + + zeeman_flag = aniso_flag = 0; + if (strcmp(arg[3],"zeeman") == 0) { if (narg != 8) error->all(FLERR,"Illegal fix zeeman command"); style = ZEEMAN; + zeeman_flag = 1; H_field = force->numeric(FLERR,arg[4]); Hx = force->numeric(FLERR,arg[5]); Hy = force->numeric(FLERR,arg[6]); @@ -71,6 +75,7 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a } else if (strcmp(arg[3],"anisotropy") == 0) { if (narg != 8) error->all(FLERR,"Illegal fix anisotropy command"); style = ANISOTROPY; + aniso_flag = 1; Ka = force->numeric(FLERR,arg[4]); Kax = force->numeric(FLERR,arg[5]); Kay = force->numeric(FLERR,arg[6]); @@ -88,6 +93,8 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a FixForceSpin::~FixForceSpin() { + memory->destroy(spi); + memory->destroy(fmi); delete [] magstr; } @@ -133,7 +140,10 @@ void FixForceSpin::init() // set magnetic field components once and for all if (varflag == CONSTANT) set_magneticforce(); - + + memory->create(spi,3,"forcespin:spi"); + memory->create(fmi,3,"forcespin:fmi"); + } /* ---------------------------------------------------------------------- */ @@ -160,7 +170,7 @@ void FixForceSpin::post_force(int vflag) set_magneticforce(); //Update value of the mag. field if time-dependent } - double **x = atom->x; +// double **x = atom->x; double **sp = atom->sp; double *mumag = atom->mumag; double **fm = atom->fm; @@ -169,22 +179,42 @@ void FixForceSpin::post_force(int vflag) eflag = 0; emag = 0.0; - - if (style == ZEEMAN) { - for (int i = 0; i < nlocal; i++) { - fm[i][0] += mumag[i]*xmag; - fm[i][1] += mumag[i]*ymag; - fm[i][2] += mumag[i]*zmag; - } - } - if (style == ANISOTROPY) { - for (int i = 0; i < nlocal; i++) { - scalar = Kax*sp[i][0] + Kay*sp[i][1] + Kaz*sp[i][2]; - fm[i][0] += scalar*xmag; - fm[i][1] += scalar*ymag; - fm[i][2] += scalar*zmag; - } - } + + for (int i = 0; i < nlocal; i++) { + fmi[0] = fmi[1] = fmi[2] = 0.0; + //if (style == ZEEMAN) { + if (zeeman_flag) { + compute_zeeman(i,fmi); + } + //if (style == ANISOTROPY) { + if (aniso_flag) { + spi[0] = sp[i][0]; + spi[1] = sp[i][1]; + spi[2] = sp[i][2]; + compute_anisotropy(i,spi,fmi); + } + fm[i][0] += fmi[0]; + fm[i][1] += fmi[1]; + fm[i][2] += fmi[2]; + } +} + + +/* ---------------------------------------------------------------------- */ +void FixForceSpin::compute_zeeman(int i, double *fmi) +{ +double *mumag = atom->mumag; + fmi[0] += mumag[i]*xmag; + fmi[1] += mumag[i]*ymag; + fmi[2] += mumag[i]*zmag; +} + +void FixForceSpin::compute_anisotropy(int i, double * spi, double *fmi) +{ + double scalar = Kax*spi[0] + Kay*spi[1] + Kaz*spi[2]; + fmi[0] += scalar*xmag; + fmi[1] += scalar*ymag; + fmi[2] += scalar*zmag; } /* ---------------------------------------------------------------------- */ diff --git a/src/SPIN/fix_force_spin.h b/src/SPIN/fix_force_spin.h index 75805a7734..a422abad2c 100644 --- a/src/SPIN/fix_force_spin.h +++ b/src/SPIN/fix_force_spin.h @@ -37,6 +37,10 @@ class FixForceSpin : public Fix { virtual void post_force_respa(int, int, int); double compute_scalar(); + int zeeman_flag, aniso_flag; + void compute_zeeman(int, double *); + void compute_anisotropy(int, double *, double *); + protected: int style; @@ -58,8 +62,10 @@ class FixForceSpin : public Fix { double Ka; //Magnetic anisotropy intensity and direction double Kax, Kay, Kaz; + double *spi, *fmi; //Temp var. in compute + void set_magneticforce(); - + }; } diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index f09d35c854..ab76051f08 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -145,25 +145,25 @@ void FixLangevinSpin::post_force(int vflag) // 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]; - sy = sp[i][1]; - sz = sp[i][2]; - - fmx = fm[i][0]; - fmy = fm[i][1]; - fmz = fm[i][2]; + sx = sp[i][0]; + sy = sp[i][1]; + sz = sp[i][2]; + + fmx = fm[i][0]; + fmy = fm[i][1]; + fmz = fm[i][2]; - cpx = fmy*sz - fmz*sy;//Computing cross product - cpy = fmz*sx - fmx*sz; - cpz = fmx*sy - fmy*sx; + cpx = fmy*sz - fmz*sy;//Computing cross product + cpy = fmz*sx - fmx*sz; + cpz = fmx*sy - fmy*sx; - fmx -= alpha_t*cpx;//Taking the damping value away - fmy -= alpha_t*cpy; - fmz -= alpha_t*cpz; + fmx -= alpha_t*cpx;//Taking the damping value away + fmy -= alpha_t*cpy; + fmz -= alpha_t*cpz; - fm[i][0] = fmx; - fm[i][1] = fmy; - fm[i][2] = fmz; + fm[i][0] = fmx; + fm[i][1] = fmy; + fm[i][2] = fmz; } //apply thermal effects adding random fields to fm diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index 9e0f3f934c..d2b1ce328b 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -26,8 +26,15 @@ #include "math_const.h" #include "modify.h" +//Add headers (see delete later) #include "pair.h" #include "timer.h" +#include "integrate.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "pair_spin.h" +#include "memory.h" +#include "fix_force_spin.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -57,7 +64,24 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : // error checks if (extra == SPIN && !atom->mumag_flag) error->all(FLERR,"Fix nve/spin requires spin attribute mumag"); - + +#if defined SEQNEI + lockpairspin = NULL; + lockforcespin = NULL; + exch_flag = dmi_flag = me_flag = 0; + zeeman_flag = aniso_flag = 0; +#endif +} + +/* ---------------------------------------------------------------------- */ +FixNVESpin::~FixNVESpin(){ +#if defined SEQNEI + //delete lockpairspin; + //delete lockforcespin; + memory->destroy(spi); + memory->destroy(fmi); + memory->destroy(fmj); +#endif } /* ---------------------------------------------------------------------- */ @@ -68,6 +92,26 @@ void FixNVESpin::init() dts = update->dt; + #if defined SEQNEI + lockpairspin = (PairSpin *) force->pair; + + memory->create(spi,3,"nves:spi"); + memory->create(fmi,3,"nves:fmi"); + memory->create(fmj,3,"nves:fmj"); + + int iforce; + for (iforce = 0; iforce < modify->nfix; iforce++) + if (strstr(modify->fix[iforce]->style,"force/spin")) break; + lockforcespin = (FixForceSpin *) modify->fix[iforce]; + + exch_flag = lockpairspin->exch_flag; + dmi_flag = lockpairspin->dmi_flag; + me_flag = lockpairspin->me_flag; + + zeeman_flag = lockforcespin->zeeman_flag; + aniso_flag = lockforcespin->aniso_flag; + #endif + /*int idamp; for (idamp = 0; idamp < modify->nfix; idamp++) if (strstr(modify->fix[idamp]->style,"damping/spin")) break; @@ -101,9 +145,18 @@ void FixNVESpin::initial_integrate(int vflag) // 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); - } +#if defined SEQNEI + for (int i = 0; i < nlocal; i++){ + ComputeSpinInteractionsNei(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } +#endif +#if defined SEQ + for (int i = 0; i < nlocal; i++){ + AdvanceSingleSpin(i,0.5*dts,sp,fm); + ComputeSpinInteractions(); + } +#endif } // update half v all particles @@ -125,18 +178,293 @@ void FixNVESpin::initial_integrate(int vflag) x[i][2] += 0.5 * dtv * v[i][2]; } } +} -//#define FORCE_PRINT -#if defined FORCE_PRINT - FILE* file_force=NULL; - file_force=fopen("spin_force_Lammps.dat","a"); - fprintf(file_force,"---------------------------------- \n"); - for(int i=0;inlocal; + + int n_post_integrate = modify->n_post_integrate; + int n_pre_exchange = modify->n_pre_exchange; + int n_pre_neighbor = modify->n_pre_neighbor; + int n_pre_force = modify->n_pre_force; + int n_pre_reverse = modify->n_pre_reverse; + int n_post_force = modify->n_post_force; + int n_end_of_step = modify->n_end_of_step; + + bigint ntimestep; + ntimestep = update->ntimestep; + + //Force compute quantities + int i,j,jj,inum,jnum,itype,jtype; + int *ilist,*jlist,*numneigh,**firstneigh; + double **x = atom->x; + double **sp = atom->sp; + double **fm = atom->fm; + int *type = atom->type; + int newton_pair = force->newton_pair; + + inum = lockpairspin->list->inum; + ilist = lockpairspin->list->ilist; + numneigh = lockpairspin->list->numneigh; + firstneigh = lockpairspin->list->firstneigh; + + double xtmp,ytmp,ztmp; + double rsq,rd,delx,dely,delz; + double cut_ex_2, cut_dmi_2, cut_me_2; + cut_ex_2 = cut_dmi_2 = cut_me_2 = 0.0; + + int eflag = 1; + int vflag = 0; + int pair_compute_flag = 1; + + if (atom->sortfreq > 0) sortflag = 1; + else sortflag = 0; + if (n_post_integrate) modify->post_integrate(); + timer->stamp(Timer::MODIFY); + + // regular communication vs neighbor list rebuild + nflag = neighbor->decide(); + if (nflag == 0) { + timer->stamp(); + comm->forward_comm(); + timer->stamp(Timer::COMM); + } else { + if (n_pre_exchange) { + timer->stamp(); + modify->pre_exchange(); + timer->stamp(Timer::MODIFY); + } + //if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + if (domain->box_change) { + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + } + timer->stamp(); + comm->exchange(); + if (sortflag && ntimestep >= atom->nextsort) atom->sort(); + comm->borders(); + //if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + timer->stamp(Timer::COMM); + if (n_pre_neighbor) { + modify->pre_neighbor(); + timer->stamp(Timer::MODIFY); + } + neighbor->build(); + timer->stamp(Timer::NEIGH); + } + + ///////Force computation for spin i///////////// + i = ilist[ii]; + //Clear atom i + fm[i][0] = fm[i][1] = fm[i][2] = 0.0; + + timer->stamp(); + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + fmi[0] = fmi[1] = fmi[2] = 0.0; + fmj[0] = fmj[1] = fmj[2] = 0.0; + jlist = firstneigh[i]; + jnum = numneigh[i]; + +// printf("Test inum: %g \n",inum); +/* + //Pair interaction + for (int jj = 0; jj < inum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + itype = type[ii]; + jtype = type[j]; + + if (exch_flag) { + cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]); + if (rsq <= cut_ex_2) { + lockpairspin->compute_exchange(i,j,rsq,fmi,fmj); + } + } + if (dmi_flag) { + cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]); + if (rsq <= cut_dmi_2) { + lockpairspin->compute_dmi(i,j,fmi,fmj); + } + } + if (me_flag) { + cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]); + if (rsq <= cut_me_2) { + lockpairspin->compute_me(i,j,fmi,fmj); + } + } + } +*/ + + //Pair interaction + int natom = nlocal + atom->nghost; + for (int k = 0; k < natom; k++) { + delx = xtmp - x[k][0]; + dely = ytmp - x[k][1]; + delz = ztmp - x[k][2]; + rsq = delx*delx + dely*dely + delz*delz; + itype = type[ii]; + jtype = type[k]; + + if (exch_flag) { + cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]); + if (rsq <= cut_ex_2) { + lockpairspin->compute_exchange(i,k,rsq,fmi,fmj); + } + } + if (dmi_flag) { + cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]); + if (rsq <= cut_dmi_2) { + lockpairspin->compute_dmi(i,k,fmi,fmj); + } + } + if (me_flag) { + cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]); + if (rsq <= cut_me_2) { + lockpairspin->compute_me(i,k,fmi,fmj); + } + } + } + + + //post force + if (zeeman_flag) { + lockforcespin->compute_zeeman(i,fmi); + } + if (aniso_flag) { + spi[0] = sp[i][0]; + spi[1] = sp[i][1]; + spi[2] = sp[i][2]; + lockforcespin->compute_anisotropy(i,spi,fmi); + } + + //Replace the force by its new value + fm[i][0] = fmi[0]; + fm[i][1] = fmi[1]; + fm[i][2] = fmi[2]; + +} #endif +/* ---------------------------------------------------------------------- */ +void FixNVESpin::ComputeSpinInteractions() +{ + int nflag,sortflag; + int nlocal = atom->nlocal; + + int n_post_integrate = modify->n_post_integrate; + int n_pre_exchange = modify->n_pre_exchange; + int n_pre_neighbor = modify->n_pre_neighbor; + int n_pre_force = modify->n_pre_force; + int n_pre_reverse = modify->n_pre_reverse; + int n_post_force = modify->n_post_force; + int n_end_of_step = modify->n_end_of_step; + + bigint ntimestep; + ntimestep = update->ntimestep; + + //int eflag = update->integrate->eflag; + //int vflag = update->integrate->vflag; + int eflag = 1; + int vflag = 0; + int pair_compute_flag = 1; + + if (atom->sortfreq > 0) sortflag = 1; + else sortflag = 0; + if (n_post_integrate) modify->post_integrate(); + timer->stamp(Timer::MODIFY); + + // regular communication vs neighbor list rebuild + nflag = neighbor->decide(); + if (nflag == 0) { + timer->stamp(); + comm->forward_comm(); + timer->stamp(Timer::COMM); + } else { + if (n_pre_exchange) { + timer->stamp(); + modify->pre_exchange(); + timer->stamp(Timer::MODIFY); + } + //if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + if (domain->box_change) { + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + } + timer->stamp(); + comm->exchange(); + if (sortflag && ntimestep >= atom->nextsort) atom->sort(); + comm->borders(); + //if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + timer->stamp(Timer::COMM); + if (n_pre_neighbor) { + modify->pre_neighbor(); + timer->stamp(Timer::MODIFY); + } + neighbor->build(); + timer->stamp(Timer::NEIGH); + } + + // force computations + // important for pair to come before bonded contributions + // since some bonded potentials tally pairwise energy/virial + // and Pair:ev_tally() needs to be called before any tallying + + size_t nbytes; + nbytes = sizeof(double) * nlocal; + if (force->newton) nbytes += sizeof(double) * atom->nghost; + + atom->avec->force_clear(0,nbytes); + + timer->stamp(); + + if (n_pre_force) { + modify->pre_force(vflag); + timer->stamp(Timer::MODIFY); + } + + if (pair_compute_flag) { + force->pair->compute(eflag,vflag); + timer->stamp(Timer::PAIR); + } + + /*if (kspace_compute_flag) { + force->kspace->compute(eflag,vflag); + timer->stamp(Timer::KSPACE); + }*/ + + if (n_pre_reverse) { + modify->pre_reverse(eflag,vflag); + timer->stamp(Timer::MODIFY); + } + + // reverse communication of forces + + if (force->newton) { + comm->reverse_comm(); + timer->stamp(Timer::COMM); + } + + // force modifications + + if (n_post_force) modify->post_force(vflag); + timer->stamp(Timer::MODIFY); + } /* ---------------------------------------------------------------------- */ @@ -218,8 +546,17 @@ void FixNVESpin::final_integrate() // 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); - } +#if defined SEQNEI + for (int i = nlocal-1; i >= 0; i--){ + ComputeSpinInteractionsNei(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } +#endif +#if defined SEQ + for (int i = nlocal-1; i >= 0; i--){ + AdvanceSingleSpin(i,0.5*dts,sp,fm); + ComputeSpinInteractions(); + } +#endif } } diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 31ef30de3e..8f516b2409 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -25,24 +25,34 @@ FixStyle(nve/spin,FixNVESpin) namespace LAMMPS_NS { class FixNVESpin : public FixNVE { - friend class FixSpinDamping; public: FixNVESpin(class LAMMPS *, int, char **); - virtual ~FixNVESpin() {} + virtual ~FixNVESpin(); void init(); virtual void initial_integrate(int); void AdvanceSingleSpin(int, double, double **, double **); virtual void final_integrate(); -//Sorting atoms/spins routine -void SortSpins(); +//#define SEQ +#define SEQNEI + void ComputeSpinInteractions(); + void ComputeSpinInteractionsNei(int); protected: int extra; double dts; //double alpha_t; - + +#if defined SEQNEI + private: + int exch_flag, dmi_flag, me_flag; + int zeeman_flag, aniso_flag; + class PairSpin *lockpairspin; + double *spi, *fmi, *fmj; //Temp var. for compute + class FixForceSpin *lockforcespin; +#endif + //private: //class FixSpinDamping *lockspindamping; }; diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index f09b4d05fd..cdd9ca8e0d 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -25,6 +25,9 @@ #include "update.h" #include +//Add. lib. for full neighb. list +#include "neigh_request.h" + using namespace LAMMPS_NS; using namespace MathConst; @@ -34,7 +37,9 @@ PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; exch_flag = 0; - dmi_flag = 0; + dmi_flag = 0; + me_flag = 0; + } /* ---------------------------------------------------------------------- */ @@ -76,7 +81,7 @@ void PairSpin::compute(int eflag, int vflag) double evdwl,ecoul; double xtmp,ytmp,ztmp; double fmix,fmiy,fmiz,fmjx,fmjy,fmjz; - double cut_ex,cut_dmi,cut_me; + double cut_ex_2,cut_dmi_2,cut_me_2; double rsq,rd,delx,dely,delz; int *ilist,*jlist,*numneigh,**firstneigh; @@ -120,28 +125,27 @@ void PairSpin::compute(int eflag, int vflag) 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 itype = type[i]; jtype = type[j]; //Exchange interaction if (exch_flag) { - cut_ex = cut_spin_exchange[itype][jtype]; - if (rd <= cut_ex) { + cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype]; + if (rsq <= cut_ex_2) { compute_exchange(i,j,rsq,fmi,fmj); } } //DM interaction if (dmi_flag){ - cut_dmi = cut_spin_dmi[itype][jtype]; - if (rd <= cut_dmi){ + cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype]; + if (rsq <= cut_dmi_2){ compute_dmi(i,j,fmi,fmj); } } //ME interaction if (me_flag){ - cut_me = cut_spin_me[itype][jtype]; - if (rd <= cut_me){ + cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype]; + if (rsq <= cut_me_2){ compute_me(i,j,fmi,fmj); } } @@ -162,7 +166,7 @@ void PairSpin::compute(int eflag, int vflag) } /* ---------------------------------------------------------------------- */ -inline void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj) +void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj) { int *type = atom->type; int itype, jtype; @@ -184,13 +188,12 @@ inline void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, d 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) +void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj) { - int *type = atom->type; int itype, jtype; double **sp = atom->sp; @@ -212,7 +215,7 @@ inline void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj) } /* ---------------------------------------------------------------------- */ -inline void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) +void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) { int *type = atom->type; int itype, jtype; @@ -247,8 +250,6 @@ inline void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) 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]); - } /* ---------------------------------------------------------------------- @@ -312,6 +313,7 @@ void PairSpin::settings(int narg, char **arg) if (setflag[i][j]) { cut_spin_exchange[i][j] = cut_spin_pair_global; cut_spin_dmi[i][j] = cut_spin_pair_global; + cut_spin_me[i][j] = cut_spin_pair_global; } } @@ -438,6 +440,13 @@ void PairSpin::init_style() error->all(FLERR,"Pair spin requires atom attributes sp, mumag"); neighbor->request(this,instance_me); + +#define FULLNEI +#if defined FULLNEI + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +#endif } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index 8b3e174de6..8df192b69f 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -39,19 +39,22 @@ class PairSpin : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); - 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 *); + void compute_exchange(int, int, double, double *, double *); + void compute_dmi(int, int, double *, double *); + void compute_me(int, int, double *, double *); - protected: + //Test for seq. integ. + //protected: int exch_flag,dmi_flag,me_flag; double cut_spin_pair_global; double cut_spin_dipolar_global; double **cut_spin_exchange; //cutting distance exchange - double **cut_spin_dmi; //cutting distance dmi - double **cut_spin_me; //cutting distance me - + double **cut_spin_dmi; //cutting distance dmi + double **cut_spin_me; //cutting distance me + + //Test for seq. integ. + protected: double **J_1, **J_2, **J_3; //exchange coeffs Jij //J1 in eV, J2 adim and J3 in Ang double **DM; diff --git a/src/modify.cpp b/src/modify.cpp index 64970f2cf9..86cfdda8b6 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -1,4 +1,5 @@ /* ---------------------------------------------------------------------- +eoundary p f f LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov