forked from lijiext/lammps
First version of the parallel algorithm
Performed by sectoring (1, 2, 4, or 8 chuncks) each process.
This commit is contained in:
julient31
parent
b934621651
commit
8a56b8ad3a
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -11,6 +11,7 @@
|
|||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <mpi.h>
|
||||
#include <string.h>
|
||||
#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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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++) {
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
|
|
|||
Loading…
Reference in New Issue