diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index d31f4691d1..61ead88129 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -32,10 +32,9 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) molecular = 0; mass_type = 1; //check why - //comm_x_only = 0; - comm_x_only = 1; - //comm_f_only = 1; - comm_f_only = 0; + comm_x_only = 0; + comm_f_only = 1; + size_forward = 7; size_reverse = 6; size_border = 11; diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index 370a2848a7..9516bcde90 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include #include #include "compute_spin.h" #include "atom.h" @@ -112,7 +113,7 @@ void ComputeSpin::compute_vector() MPI_Allreduce(&magenergy,&magenergytot,1,MPI_DOUBLE,MPI_SUM,world); 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_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&countsp,&countsptot,1,MPI_INT,MPI_SUM,world); double scale = 1.0/countsptot; magtot[0] *= scale; diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp index 27408456d6..44948f8ea3 100644 --- a/src/SPIN/fix_force_spin.cpp +++ b/src/SPIN/fix_force_spin.cpp @@ -114,9 +114,9 @@ int FixForceSpin::setmask() void FixForceSpin::init() { - double hbar = force->hplanck/MY_2PI; //eV/(rad.THz) - double mub = 5.78901e-5; //in eV/T - double gyro = mub/hbar; //in rad.THz/T + 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 H_field *= gyro; //in rad.THz Ka /= hbar; //in rad.THz @@ -174,7 +174,7 @@ void FixForceSpin::post_force(int vflag) double **sp = atom->sp; double *mumag = atom->mumag; double **fm = atom->fm; - int nlocal = atom->nlocal; + const int nlocal = atom->nlocal; double scalar; eflag = 0; diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index ab76051f08..8079000d59 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -63,9 +63,29 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) : alpha_l = force->numeric(FLERR,arg[5]); seed = force->inumeric(FLERR,arg[6]); - if (alpha_t < 0.0) error->all(FLERR,"Fix langevin/spin transverse damping must be >= 0.0"); - if (alpha_l < 0.0) error->all(FLERR,"Fix langevin/spin transverse damping must be >= 0.0"); - if (seed <= 0) error->all(FLERR,"Illegal fix langevin/spin seed must be > 0"); + if (alpha_t < 0.0) { + error->all(FLERR,"Illegal fix/langevin/spin command"); + } else if (alpha_t == 0.0) { + tdamp_flag = 0; + } else { + tdamp_flag = 1; + } + + if (alpha_l < 0.0) { + error->all(FLERR,"Illegal fix/langevin/spin command"); + } else if (alpha_l == 0.0) { + ldamp_flag = 0; + } else { + ldamp_flag = 1; + } + + if (temp < 0.0) { + error->all(FLERR,"Illegal fix/langevin/spin command"); + } else if (temp == 0.0) { + temp_flag = 0; + } else { + temp_flag = 1; + } // initialize Marsaglia RNG with processor-unique seed //random = new RanMars(lmp,seed + comm->me); @@ -77,6 +97,8 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) : FixLangevinSpin::~FixLangevinSpin() { + memory->destroy(spi); + memory->destroy(fmi); delete random; } @@ -106,6 +128,9 @@ void FixLangevinSpin::init() } if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin should come after all other spin fixes"); + memory->create(spi,3,"pair:spi"); + memory->create(fmi,3,"pair:fmi"); + dts = update->dt; Gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t)); @@ -135,7 +160,7 @@ void FixLangevinSpin::post_force(int vflag) double **sp = atom->sp; double **fm = atom->fm; int *mask = atom->mask; - int nlocal = atom->nlocal; + const int nlocal = atom->nlocal; double sx, sy, sz; double fmx, fmy, fmz; @@ -143,54 +168,69 @@ void FixLangevinSpin::post_force(int vflag) double rx, ry, rz; // add the damping to the effective field of each spin - for (int i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - sx = sp[i][0]; - sy = sp[i][1]; - sz = sp[i][2]; + spi[0] = sp[i][0]; + spi[1] = sp[i][1]; + spi[2] = 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; - - 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; - } + fmi[0] = fm[i][0]; + fmi[1] = fm[i][1]; + fmi[2] = fm[i][2]; - //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(); - #else - rx = sigma*(random->uniform() - 0.5); - ry = sigma*(random->uniform() - 0.5); - rz = sigma*(random->uniform() - 0.5); - #endif + if (tdamp_flag) { + add_tdamping(spi,fmi); + } - fm[i][0] += rx;//Adding random field - fm[i][1] += ry; - fm[i][2] += rz; - - fm[i][0] *= Gil_factor;//Multiplying by Gilbert's prefactor - fm[i][1] *= Gil_factor; - fm[i][2] *= Gil_factor; + if (temp_flag) { + add_temperature(fmi); + } - } + fm[i][0] = fmi[0]; + fm[i][1] = fmi[1]; + fm[i][2] = fmi[2]; + } + } + } +/* ---------------------------------------------------------------------- */ +void FixLangevinSpin::add_tdamping(double *spi, double *fmi) +{ + double cpx = fmi[1]*spi[2] - fmi[2]*spi[1]; + double cpy = fmi[2]*spi[0] - fmi[0]*spi[2]; + double cpz = fmi[0]*spi[1] - fmi[1]*spi[0]; + + fmi[0] -= alpha_t*cpx;//Taking the damping value away + fmi[1] -= alpha_t*cpy; + fmi[2] -= alpha_t*cpz; +} + +/* ---------------------------------------------------------------------- */ +void FixLangevinSpin::add_temperature(double *fmi) +{ +//#define GAUSSIAN_R +#if defined GAUSSIAN_R + double rx = sigma*random->gaussian();//Drawing random dist + double ry = sigma*random->gaussian(); + double rz = sigma*random->gaussian(); +#else + double rx = sigma*(random->uniform() - 0.5); + double ry = sigma*(random->uniform() - 0.5); + double rz = sigma*(random->uniform() - 0.5); +#endif + + fmi[0] += rx;//Adding random field + fmi[1] += ry; + fmi[2] += rz; + + fmi[0] *= Gil_factor;//Multiplying by Gilbert's prefactor + fmi[1] *= Gil_factor; + fmi[2] *= Gil_factor; + +} + + /* ---------------------------------------------------------------------- */ void FixLangevinSpin::post_force_respa(int vflag, int ilevel, int iloop) diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h index ffd94745e2..a6b9bc5df3 100644 --- a/src/SPIN/fix_langevin_spin.h +++ b/src/SPIN/fix_langevin_spin.h @@ -33,10 +33,12 @@ class FixLangevinSpin : public Fix { void setup(int); virtual void post_force(int); void post_force_respa(int, int, int); + void add_tdamping(double *, double *); + void add_temperature(double *); + int tdamp_flag, ldamp_flag, temp_flag; protected: - //First mag. quantities - int transv_damp_flag, long_damp_flag; //Flags for transverse or longitudinal mag. dampings + double *spi, *fmi; double alpha_t, alpha_l; //Transverse and long. damping value double dts,temp,D,sigma;//timestep, temp, noise double Gil_factor;//Gilbert's prefactor diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index d2b1ce328b..11ec1ccf11 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -27,14 +27,13 @@ #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.h" #include "pair_spin.h" #include "memory.h" #include "fix_force_spin.h" +#include "fix_langevin_spin.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -65,23 +64,32 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : 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 + tdamp_flag = temp_flag = 0; + + lockpairspin = NULL; + lockforcespin = NULL; + locklangevinspin = NULL; } /* ---------------------------------------------------------------------- */ FixNVESpin::~FixNVESpin(){ -#if defined SEQNEI //delete lockpairspin; //delete lockforcespin; + memory->destroy(xi); +#if defined SECTORING + memory->destroy(sec); + memory->destroy(rsec); + memory->destroy(seci); +#endif +#if defined SECTOR_PRINT + fclose(file_sect); +#endif memory->destroy(spi); + memory->destroy(spj); memory->destroy(fmi); memory->destroy(fmj); -#endif } /* ---------------------------------------------------------------------- */ @@ -91,36 +99,56 @@ void FixNVESpin::init() FixNVE::init(); dts = update->dt; - - #if defined SEQNEI - lockpairspin = (PairSpin *) force->pair; - + memory->create(xi,3,"nves:xi"); +#if defined SECTORING + memory->create(sec,3,"nves:sec"); + memory->create(rsec,3,"nves:rsec"); + memory->create(seci,3,"nves:seci"); +#endif memory->create(spi,3,"nves:spi"); + memory->create(spj,3,"nves:spj"); memory->create(fmi,3,"nves:fmi"); memory->create(fmj,3,"nves:fmj"); + lockpairspin = (PairSpin *) force->pair; + int iforce; for (iforce = 0; iforce < modify->nfix; iforce++) if (strstr(modify->fix[iforce]->style,"force/spin")) break; lockforcespin = (FixForceSpin *) modify->fix[iforce]; + for (iforce = 0; iforce < modify->nfix; iforce++) + if (strstr(modify->fix[iforce]->style,"langevin/spin")) break; + locklangevinspin = (FixLangevinSpin *) 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; - if (idamp == modify->nfix) - error->all(FLERR,"Integration of spin systems requires use of fix damping (set damping to 0.0 for NVE)"); - - lockspindamping = (FixSpinDamping *) modify->fix[idamp]; - alpha_t = lockspindamping->get_damping(0); - */ + + tdamp_flag = locklangevinspin->tdamp_flag; + temp_flag = locklangevinspin->temp_flag; + + +#if defined SECTORING + sectoring(); +#endif + +#if defined SECTOR_PRINT + file_sect=fopen("sectoring.lammpstrj", "w"); + fprintf(file_sect,"ITEM: TIMESTEP\n"); + fprintf(file_sect,"%g\n",0.0); + fprintf(file_sect,"ITEM: NUMBER OF ATOMS\n"); + //int natoms = atom->natoms; + int natoms = atom->nlocal; + fprintf(file_sect,"%d\n",natoms); + fprintf(file_sect,"ITEM: BOX BOUNDS\n"); + for(int d=0; d<3; d++) fprintf(file_sect,"%lf %lf\n",domain->boxlo[d],domain->boxhi[d]); + fprintf(file_sect,"ITEM: ATOMS type x y z vx vy vz\n"); +#endif + } /* ---------------------------------------------------------------------- */ @@ -142,22 +170,6 @@ void FixNVESpin::initial_integrate(int vflag) 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) { -#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 for (int i = 0; i < nlocal; i++) { @@ -170,6 +182,59 @@ void FixNVESpin::initial_integrate(int vflag) } } + // 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]; + } + } + +#if defined SECTORING + int nseci; + // Seq. update spins for all particles + if (extra == SPIN) { + for (int j = 0; j < nsectors; j++) { + comm->forward_comm(); + for (int i = 0; i < nlocal; i++) { + xi[0] = x[i][0]; + xi[1] = x[i][1]; + xi[2] = x[i][2]; + nseci = coords2sector(xi); + if (j != nseci) continue; + ComputeSpinInteractionsNeigh(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + } + for (int j = nsectors-1; j >= 0; j--) { + comm->forward_comm(); + for (int i = nlocal-1; i >= 0; i--) { + xi[0] = x[i][0]; + xi[1] = x[i][1]; + xi[2] = x[i][2]; + nseci = coords2sector(xi); + if (j != nseci) continue; + ComputeSpinInteractionsNeigh(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + } + } +#else + // Seq. update spins for all particles + if (extra == SPIN) { + for (int i = 0; i < nlocal; i++){ + ComputeSpinInteractionsNeigh(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + + for (int i = nlocal-1; i >= 0; i--){ + ComputeSpinInteractionsNeigh(i); + AdvanceSingleSpin(i,0.5*dts,sp,fm); + } + } +#endif + // update x for all particles for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -177,26 +242,32 @@ void FixNVESpin::initial_integrate(int vflag) x[i][1] += 0.5 * dtv * v[i][1]; x[i][2] += 0.5 * dtv * v[i][2]; } - } + } + + +#if defined SECTOR_PRINT + int my_rank; + MPI_Comm_rank(world, &my_rank); + if (my_rank == 0) { + for (int j = 0; j < nsectors; j++) { + for (int i = 0; i < nlocal; i++) { + xi[0] = x[i][0]; + xi[1] = x[i][1]; + xi[2] = x[i][2]; + nseci = coords2sector(xi); + if (j != nseci) continue; + fprintf(file_sect,"%d %lf %lf %lf %lf %lf %lf\n",j,xi[0],xi[1],xi[2],0.0,0.0,1.0); + } + } + } +#endif + } -#if defined SEQNEI /* ---------------------------------------------------------------------- */ -void FixNVESpin::ComputeSpinInteractionsNei(int ii) +void FixNVESpin::ComputeSpinInteractionsNeigh(int ii) { - 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; + const int nlocal = atom->nlocal; //Force compute quantities int i,j,jj,inum,jnum,itype,jtype; @@ -205,7 +276,7 @@ void FixNVESpin::ComputeSpinInteractionsNei(int ii) double **sp = atom->sp; double **fm = atom->fm; int *type = atom->type; - int newton_pair = force->newton_pair; + const int newton_pair = force->newton_pair; inum = lockpairspin->list->inum; ilist = lockpairspin->list->ilist; @@ -221,69 +292,37 @@ void FixNVESpin::ComputeSpinInteractionsNei(int ii) 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); - } +// comm->forward_comm(); ///////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]; + + 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]; 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++) { + for (int 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]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; + delx = xi[0] - x[j][0]; + dely = xi[1] - x[j][1]; + delz = xi[2] - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; itype = type[ii]; jtype = type[j]; @@ -291,191 +330,144 @@ void FixNVESpin::ComputeSpinInteractionsNei(int ii) 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); + lockpairspin->compute_exchange(i,j,rsq,fmi,fmj,spi,spj); } } + 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); + lockpairspin->compute_dmi(i,j,fmi,fmj,spi,spj); } } + 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); + lockpairspin->compute_me(i,j,fmi,fmj,spi,spj); } } } -*/ - - //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); } - + + if (tdamp_flag) { + locklangevinspin->add_tdamping(spi,fmi); + } + + if (temp_flag) { + locklangevinspin->add_temperature(fmi); + } + //Replace the force by its new value fm[i][0] = fmi[0]; fm[i][1] = fmi[1]; fm[i][2] = fmi[2]; } -#endif +#if defined SECTORING /* ---------------------------------------------------------------------- */ -void FixNVESpin::ComputeSpinInteractions() +void FixNVESpin::sectoring() { - 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); + double sublo[3],subhi[3]; + double* sublotmp = domain->sublo; + double* subhitmp = domain->subhi; + for (int dim = 0 ; dim<3 ; dim++) { + sublo[dim]=sublotmp[dim]; + subhi[dim]=subhitmp[dim]; } - // 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 + const double rsx = subhi[0] - sublo[0]; + const double rsy = subhi[1] - sublo[1]; + const double rsz = subhi[2] - sublo[2]; - size_t nbytes; - nbytes = sizeof(double) * nlocal; - if (force->newton) nbytes += sizeof(double) * atom->nghost; + const double rv = lockpairspin->cut_spin_pair_global; - atom->avec->force_clear(0,nbytes); + double rax = rsx/rv; + double ray = rsy/rv; + double raz = rsz/rv; + + sec[0] = 1; + sec[1] = 1; + sec[2] = 1; + if (rax >= 2.0) sec[0] = 2; + if (ray >= 2.0) sec[1] = 2; + if (raz >= 2.0) sec[2] = 2; - timer->stamp(); + nsectors = sec[0]*sec[1]*sec[2]; - if (n_pre_force) { - modify->pre_force(vflag); - timer->stamp(Timer::MODIFY); - } + rsec[0] = rsx; + rsec[1] = rsy; + rsec[2] = rsz; + if (sec[0] == 2) rsec[0] = rsx/2.0; + if (sec[1] == 2) rsec[1] = rsy/2.0; + if (sec[2] == 2) rsec[2] = rsz/2.0; - 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); + if (2.0 * rv >= rsx && sec[0] >= 2) + error->all(FLERR,"Illegal number of sectors"); + + if (2.0 * rv >= rsy && sec[1] >= 2) + error->all(FLERR,"Illegal number of sectors"); + + if (2.0 * rv >= rsz && sec[2] >= 2) + error->all(FLERR,"Illegal number of sectors"); } +/* ---------------------------------------------------------------------- */ +int FixNVESpin::coords2sector(double *xi) +{ + int nseci; + double sublo[3]; + double* sublotmp = domain->sublo; + for (int dim = 0 ; dim<3 ; dim++) { + sublo[dim]=sublotmp[dim]; + } + + double rix = (xi[0] - sublo[0])/rsec[0]; + double riy = (xi[1] - sublo[1])/rsec[1]; + double riz = (xi[2] - sublo[2])/rsec[2]; + + seci[0] = (int)rix; + seci[1] = (int)riy; + seci[2] = (int)riz; + + if (nsectors == 1) { + nseci == 0; + } else if (nsectors == 2) { + nseci = seci[0] + seci[1] + seci[2]; + } else if (nsectors == 4) { + if (sec[1]*sec[2] == 4) { //plane normal to x + nseci = (seci[1] + 2*seci[2]); + } else if (sec[0]*sec[2] == 4) { //plane normal to y + nseci = (seci[0] + 2*seci[2]); + } else if (sec[0]*sec[1] == 4) { //plane normal to z + nseci = (seci[0] + 2*seci[1]); + } + } else if (nsectors == 8) { + nseci = (seci[0] + 2*seci[1] + 4*seci[2]); + } + + return nseci; +} + +#endif + /* ---------------------------------------------------------------------- */ void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm) { - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - int *type = atom->type; - int *mask = atom->mask; - double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy; double cp[3],g[3]; @@ -523,8 +515,6 @@ 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; @@ -543,20 +533,4 @@ 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) { -#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 8f516b2409..d77ea4e37a 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -33,28 +33,41 @@ class FixNVESpin : public FixNVE { virtual void initial_integrate(int); void AdvanceSingleSpin(int, double, double **, double **); virtual void final_integrate(); - -//#define SEQ -#define SEQNEI void ComputeSpinInteractions(); - void ComputeSpinInteractionsNei(int); + void ComputeSpinInteractionsNeigh(int); + +#define SECTORING +#if defined SECTORING + void sectoring(); + int coords2sector(double *); +#endif protected: int extra; double dts; - //double alpha_t; - -#if defined SEQNEI - private: + int exch_flag, dmi_flag, me_flag; int zeeman_flag, aniso_flag; + int tdamp_flag, temp_flag; + class PairSpin *lockpairspin; - double *spi, *fmi, *fmj; //Temp var. for compute class FixForceSpin *lockforcespin; + class FixLangevinSpin *locklangevinspin; + + double *spi, *spj, *fmi, *fmj; //Temp var. for compute + double *xi; + +#if defined SECTORING + int nsectors; + int *sec; + int *seci; + double *rsec; #endif - //private: - //class FixSpinDamping *lockspindamping; +//#define SECTOR_PRINT +#if defined SECTOR_PRINT + FILE* file_sect=NULL; +#endif }; } diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index cdd9ca8e0d..f2c9cb8d36 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -66,6 +66,8 @@ PairSpin::~PairSpin() memory->destroy(v_mey); memory->destroy(v_mez); + memory->destroy(spi); + memory->destroy(spj); memory->destroy(fmi); memory->destroy(fmj); @@ -111,12 +113,18 @@ void PairSpin::compute(int eflag, int vflag) ytmp = x[i][1]; ztmp = x[i][2]; jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; + 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]; fmi[0] = fmi[1] = fmi[2] = 0.0; fmj[0] = fmj[1] = fmj[2] = 0.0; @@ -132,21 +140,21 @@ void PairSpin::compute(int eflag, int vflag) if (exch_flag) { 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); + compute_exchange(i,j,rsq,fmi,fmj,spi,spj); } } //DM interaction if (dmi_flag){ cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype]; if (rsq <= cut_dmi_2){ - compute_dmi(i,j,fmi,fmj); + compute_dmi(i,j,fmi,fmj,spi,spj); } } //ME interaction if (me_flag){ cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype]; if (rsq <= cut_me_2){ - compute_me(i,j,fmi,fmj); + compute_me(i,j,fmi,fmj,spi,spj); } } @@ -166,11 +174,10 @@ void PairSpin::compute(int eflag, int vflag) } /* ---------------------------------------------------------------------- */ -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, double *spi, double *spj) { int *type = atom->type; int itype, jtype; - double **sp = atom->sp; double dmix,dmiy,dmiz; double Jex, ra; itype = type[i]; @@ -181,22 +188,21 @@ void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double * 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]; + fmi[0] += Jex*spj[0]; + fmi[1] += Jex*spj[1]; + fmi[2] += Jex*spj[2]; - fmj[0] += Jex*sp[i][0]; - fmj[1] += Jex*sp[i][1]; - fmj[2] += Jex*sp[i][2]; + fmj[0] += Jex*spi[0]; + fmj[1] += Jex*spi[1]; + fmj[2] += Jex*spi[2]; } /* ---------------------------------------------------------------------- */ -void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj) +void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj, double *spi, double *spj) { int *type = atom->type; int itype, jtype; - double **sp = atom->sp; double dmix,dmiy,dmiz; itype = type[i]; jtype = type[j]; @@ -205,17 +211,17 @@ void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj) 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; + fmi[0] += spj[1]*dmiz-spj[2]*dmiy; + fmi[1] += spj[2]*dmix-spj[0]*dmiz; + fmi[2] += spj[0]*dmiy-spj[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; + fmj[0] -= spi[1]*dmiz-spi[2]*dmiy; + fmj[1] -= spi[2]*dmix-spi[0]*dmiz; + fmj[2] -= spi[0]*dmiy-spi[1]*dmix; } /* ---------------------------------------------------------------------- */ -void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) +void PairSpin::compute_me(int i, int j, double *fmi, double *fmj, double *spi, double *spj) { int *type = atom->type; int itype, jtype; @@ -242,13 +248,13 @@ void PairSpin::compute_me(int i, int j, double *fmi, double *fmj) 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; + fmi[0] += spj[1]*meiz - spj[2]*meiy; + fmi[1] += spj[2]*meix - spj[0]*meiz; + fmi[2] += spj[0]*meiy - spj[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; + fmj[0] -= spi[1]*meiz - spi[2]*meiy; + fmj[1] -= spi[2]*meix - spi[0]*meiz; + fmj[2] -= spi[0]*meiy - spi[1]*meix; } @@ -283,6 +289,8 @@ void PairSpin::allocate() 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(spi,3,"pair:spi"); + memory->create(spj,3,"pair:spj"); memory->create(fmi,3,"pair:fmi"); memory->create(fmj,3,"pair:fmj"); @@ -325,6 +333,7 @@ void PairSpin::settings(int narg, char **arg) void PairSpin::coeff(int narg, char **arg) { + const double hbar = force->hplanck/MY_2PI; if (!allocated) allocate(); @@ -336,14 +345,11 @@ void PairSpin::coeff(int narg, char **arg) 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]); + const double rij = force->numeric(FLERR,arg[3]); + const double J1 = (force->numeric(FLERR,arg[4]))/hbar; + const double J2 = force->numeric(FLERR,arg[5]); + const 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++) { @@ -363,8 +369,8 @@ void PairSpin::coeff(int narg, char **arg) 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]); + 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]); @@ -374,9 +380,6 @@ void PairSpin::coeff(int narg, char **arg) 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++) { @@ -397,8 +400,8 @@ void PairSpin::coeff(int narg, char **arg) 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]); + const double rij = force->numeric(FLERR,arg[3]); + const double me = (force->numeric(FLERR,arg[4]))/hbar; double mex = force->numeric(FLERR,arg[5]); double mey = force->numeric(FLERR,arg[6]); double mez = force->numeric(FLERR,arg[7]); @@ -408,9 +411,6 @@ void PairSpin::coeff(int narg, char **arg) 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++) { diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index 8df192b69f..2c65c62264 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -39,9 +39,9 @@ class PairSpin : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); - void compute_exchange(int, int, double, double *, double *); - void compute_dmi(int, int, double *, double *); - void compute_me(int, int, double *, double *); + void compute_exchange(int, int, double, double *, double *,double *, double *); + void compute_dmi(int, int, double *, double *, double *, double *); + void compute_me(int, int, double *, double *, double *, double *); //Test for seq. integ. //protected: @@ -61,10 +61,10 @@ class PairSpin : public Pair { 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 **ME; //ME in eV + double **v_mex, **v_mey, **v_mez;//ME direction + double *spi, *spj; double *fmi, *fmj; //Temp var. in compute void allocate();