Commit 75069ec5 authored by julient31's avatar julient31
Browse files

Commit JT-1 042618 

- new version DMI
parent 392816a8

src/SPIN/pair_spin_dmi.cpp

+46 −22
Changes for src/SPIN/pair_spin_dmi.cpp: 46 added lines, 22 removed lines.
Original line number Diff line number Diff line
@@ -49,6 +49,7 @@ PairSpinDmi::PairSpinDmi(LAMMPS *lmp) : PairSpin(lmp) 
{

  single_enable = 0;

  no_virial_fdotr_compute = 1;

  lattice_flag = 0;

}



/* ---------------------------------------------------------------------- */
@@ -165,6 +166,16 @@ void PairSpinDmi::init_style() 
  }

  if (ifix == modify->nfix)

    error->all(FLERR,"pair/spin style requires nve/spin");

  

  // get the lattice_flag from nve/spin



  for (int i = 0; i < modify->nfix; i++) {

    if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {

      lockfixnvespin = (FixNVESpin *) modify->fix[i];

      lattice_flag = lockfixnvespin->lattice_flag;

    }

  }



}



/* ----------------------------------------------------------------------
@@ -194,12 +205,13 @@ void *PairSpinDmi::extract(const char *str, int &dim) 


void PairSpinDmi::compute(int eflag, int vflag)

{

  int i,j,ii,jj,inum,jnum;

  int i,j,ii,jj,inum,jnum,itype,jtype;

  double evdwl, ecoul;

  double xi[3], rij[3];

  double xi[3], rij[3], eij[3];

  double spi[3], spj[3];

  double fi[3], fmi[3];

  double rsq;

  double local_cut2;

  double rsq, inorm;

  int *ilist,*jlist,*numneigh,**firstneigh;  



  evdwl = ecoul = 0.0;
@@ -224,6 +236,7 @@ void PairSpinDmi::compute(int eflag, int vflag) 


  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    itype = type[i];

    

    xi[0] = x[i][0];

    xi[1] = x[i][1];
@@ -242,6 +255,7 @@ void PairSpinDmi::compute(int eflag, int vflag) 
   

      j = jlist[jj];

      j &= NEIGHMASK;

      jtype = type[j];



      spj[0] = sp[j][0]; 

      spj[1] = sp[j][1];
@@ -252,16 +266,27 @@ void PairSpinDmi::compute(int eflag, int vflag) 
      fi[0] = fi[1] = fi[2] = 0.0;

      fmi[0] = fmi[1] = fmi[2] = 0.0;

      rij[0] = rij[1] = rij[2] = 0.0;

      eij[0] = eij[1] = eij[2] = 0.0;

     

      rij[0] = x[j][0] - xi[0];

      rij[1] = x[j][1] - xi[1];

      rij[2] = x[j][2] - xi[2];

      rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; 

      inorm = 1.0/sqrt(rsq);

      eij[0] = rij[0]*inorm;

      eij[1] = rij[1]*inorm;

      eij[2] = rij[2]*inorm;



      local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];



      // compute magnetic and mechanical components of soc_dmi

      

      compute_dmi(i,j,rsq,fmi,spi,spj);

      if (rsq <= local_cut2) {

	compute_dmi(i,j,rsq,eij,fmi,spi,spj);

	if (lattice_flag) {

	  compute_dmi_mech(i,j,fi,spi,spj);

	}

      }



      f[i][0] += fi[0];	 

      f[i][1] += fi[1];
@@ -300,15 +325,16 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) 
  int *type = atom->type;

  double **x = atom->x;

  double **sp = atom->sp;

  double local_cut2;



  double xi[3], rij[3];

  double xi[3], rij[3], eij[3];

  double spi[3], spj[3];



  int iexchange, idmi, ineel, ime;

  int i,j,jj,inum,jnum,itype,jtype;

  int *ilist,*jlist,*numneigh,**firstneigh;



  double rsq, rd, inorm;

  double rsq, inorm;



  inum = list->inum;

  ilist = list->ilist;
@@ -316,6 +342,7 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) 
  firstneigh = list->firstneigh;



  i = ilist[ii];

  itype = type[i];



  spi[0] = sp[i][0];

  spi[1] = sp[i][1];
@@ -332,7 +359,6 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) 


    j = jlist[jj];

    j &= NEIGHMASK;

    itype = type[ii];

    jtype = type[j];



    spj[0] = sp[j][0];
@@ -344,8 +370,15 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) 
    rij[2] = x[j][2] - xi[2];

    rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];

    inorm = 1.0/sqrt(rsq);

    eij[0] = rij[0]*inorm;

    eij[1] = rij[1]*inorm;

    eij[2] = rij[2]*inorm;



    local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];



    compute_dmi(i,j,rsq,fmi,spi,spj);

    if (rsq <= local_cut2) {

      compute_dmi(i,j,rsq,eij,fmi,spi,spj);

    }



  }
@@ -355,30 +388,21 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3]) 
   compute the dmi interaction between spin i and spin j

------------------------------------------------------------------------- */



void PairSpinDmi::compute_dmi(int i, int j, double rsq, double fmi[3], double spi[3], double spj[3])

void PairSpinDmi::compute_dmi(int i, int j, double rsq, double eij[3], double fmi[3], double spi[3], double spj[3])

{

  int *type = atom->type;  

  int itype, jtype;

  double dmix, dmiy, dmiz;	

  itype = type[i];

  jtype = type[j];



  double local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];



  if (rsq <= local_cut2) {

    double dmix, dmiy, dmiz;	



    //dmix = DM[itype][jtype] * v_dmx[itype][jtype];

    //dmiy = DM[itype][jtype] * v_dmy[itype][jtype];

    //dmiz = DM[itype][jtype] * v_dmz[itype][jtype];

    dmix = v_dmx[itype][jtype];

    dmiy = v_dmy[itype][jtype];

    dmiz = v_dmz[itype][jtype];

  dmix = eij[1]*v_dmz[itype][jtype] - eij[2]*v_dmy[itype][jtype];

  dmiy = eij[2]*v_dmx[itype][jtype] - eij[2]*v_dmz[itype][jtype];

  dmiz = eij[0]*v_dmy[itype][jtype] - eij[1]*v_dmx[itype][jtype];



  fmi[0] += (spj[1]*dmiz - spj[2]*dmiy);

  fmi[1] += (spj[2]*dmix - spj[0]*dmiz);

  fmi[2] += (spj[0]*dmiy - spj[1]*dmix);



  }

}



/* ----------------------------------------------------------------------

src/SPIN/pair_spin_dmi.h

+4 −1
Changes for src/SPIN/pair_spin_dmi.h: 4 added lines, 1 removed line.
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ class PairSpinDmi : public PairSpin { 
  void compute(int, int);

  void compute_single_pair(int, double *);



  void compute_dmi(int, int, double, double *, double *, double *);

  void compute_dmi(int, int, double, double *, double *, double *, double *);

  void compute_dmi_mech(int, int, double *, double *, double *);

  

  void write_restart(FILE *);
@@ -52,6 +52,9 @@ class PairSpinDmi : public PairSpin { 
  double **v_dmx, **v_dmy, **v_dmz;	// dmi direction

  double **cut_spin_dmi;      		// cutoff distance dmi



  int lattice_flag;                     // flag for mech force computation

  class FixNVESpin *lockfixnvespin;     // ptr to FixNVESpin for setups



  void allocate();

};

src/SPIN/pair_spin_exchange.cpp

+9 −22
Changes for src/SPIN/pair_spin_exchange.cpp: 9 added lines, 22 removed lines.
Original line number Diff line number Diff line
@@ -204,16 +204,15 @@ void PairSpinExchange::compute(int eflag, int vflag) 
  int i,j,ii,jj,inum,jnum,itype,jtype;  

  double evdwl,ecoul;

  double xi[3], rij[3], eij[3];

  double spi[3], spj[3];

  double fi[3], fmi[3];

  double cut_ex_2,cut_spin_exchange_global2;

  double local_cut2;

  double rsq, inorm;

  int *ilist,*jlist,*numneigh,**firstneigh;  

  double spi[3],spj[3];



  evdwl = ecoul = 0.0;

  if (eflag || vflag) ev_setup(eflag,vflag);

  else evflag = vflag_fdotr = 0;

  cut_spin_exchange_global2 = cut_spin_exchange_global*cut_spin_exchange_global;



  double **x = atom->x;

  double **f = atom->f;
@@ -268,10 +267,11 @@ void PairSpinExchange::compute(int eflag, int vflag) 
      eij[1] = inorm*rij[1]; 

      eij[2] = inorm*rij[2]; 



      cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];

      local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];



      // compute exchange interaction

      if (rsq <= cut_ex_2) {

      

      if (rsq <= local_cut2) {

	compute_exchange(i,j,rsq,fmi,spi,spj);

        if (lattice_flag) {

	  compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
@@ -285,7 +285,6 @@ void PairSpinExchange::compute(int eflag, int vflag) 
      fm[i][1] += fmi[1];	  	  

      fm[i][2] += fmi[2];



      // check newton pair  =>  see if needs correction

      if (newton_pair || j < nlocal) {

	f[j][0] -= fi[0];	 

        f[j][1] -= fi[1];
@@ -293,7 +292,7 @@ void PairSpinExchange::compute(int eflag, int vflag) 
      }



      if (eflag) {

	if (rsq <= cut_ex_2) {

	if (rsq <= local_cut2) {

	  evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);

	  evdwl *= hbar;

	}
@@ -375,21 +374,17 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3]) 
}



/* ----------------------------------------------------------------------

   compute the magnetic exchange interaction between spin i and spin j

   compute exchange interaction between spins i and j

------------------------------------------------------------------------- */



void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3], double spi[3], double spj[3])

{

  int *type = atom->type;  

  int itype, jtype;

  double Jex, ra;

  itype = type[i];

  jtype = type[j];



  double local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];

  

  if (rsq <= local_cut2) {

    double Jex, ra;



  ra = rsq/J3[itype][jtype]/J3[itype][jtype]; 

  Jex = 4.0*J1_mag[itype][jtype]*ra;

  Jex *= (1.0-J2[itype][jtype]*ra);
@@ -398,8 +393,6 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3], 
  fmi[0] += Jex*spj[0];

  fmi[1] += Jex*spj[1];

  fmi[2] += Jex*spj[2];



  }

}



/* ----------------------------------------------------------------------
@@ -410,13 +403,10 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ri 
{

  int *type = atom->type;  

  int itype, jtype;

  double Jex, Jex_mech, ra, rr, iJ3;

  itype = type[i];

  jtype = type[j];

    

  double local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];



  if (rsq <= local_cut2) {

    double Jex, Jex_mech, ra, rr, iJ3;

  Jex = J1_mech[itype][jtype];

  iJ3 = 1.0/(J3[itype][jtype]*J3[itype][jtype]);
@@ -430,9 +420,6 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ri 
  fi[0] -= Jex_mech*rij[0];

  fi[1] -= Jex_mech*rij[1];

  fi[2] -= Jex_mech*rij[2];



  }

          

}



/* ----------------------------------------------------------------------

src/SPIN/pair_spin_me.cpp

+32 −9
Changes for src/SPIN/pair_spin_me.cpp: 32 added lines, 9 removed lines.
Original line number Diff line number Diff line
@@ -49,6 +49,7 @@ PairSpinMe::PairSpinMe(LAMMPS *lmp) : PairSpin(lmp) 
{

  single_enable = 0;

  no_virial_fdotr_compute = 1;

  lattice_flag = 0;

}



/* ---------------------------------------------------------------------- */
@@ -163,6 +164,16 @@ void PairSpinMe::init_style() 
  }

  if (ifix == modify->nfix)

    error->all(FLERR,"pair/spin style requires nve/spin");

  

  // get the lattice_flag from nve/spin



  for (int i = 0; i < modify->nfix; i++) {

    if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {

      lockfixnvespin = (FixNVESpin *) modify->fix[i];

      lattice_flag = lockfixnvespin->lattice_flag;

    }

  }



}



/* ----------------------------------------------------------------------
@@ -195,11 +206,12 @@ void PairSpinMe::compute(int eflag, int vflag) 
{

  int i,j,ii,jj,inum,jnum,itype,jtype;  

  double evdwl, ecoul;

  double xi[3], rij[3];

  double xi[3], rij[3], eij[3];

  double spi[3], spj[3];

  double fi[3], fj[3];

  double fmi[3], fmj[3];

  double rsq, rd, inorm;

  double local_cut2;

  double rsq, inorm;

  int *ilist,*jlist,*numneigh,**firstneigh;  



  evdwl = ecoul = 0.0;
@@ -224,6 +236,8 @@ void PairSpinMe::compute(int eflag, int vflag) 


  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    itype = type[i];



    xi[0] = x[i][0];

    xi[1] = x[i][1];

    xi[2] = x[i][2];
@@ -256,17 +270,22 @@ void PairSpinMe::compute(int eflag, int vflag) 
      rij[2] = x[j][2] - xi[2];

      rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; 

      inorm = 1.0/sqrt(rsq);

      rij[0] *= inorm;

      rij[1] *= inorm;

      rij[2] *= inorm;

      eij[0] *= inorm;

      eij[1] *= inorm;

      eij[2] *= inorm;



      itype = type[i];

      jtype = type[j];



      local_cut2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];



      // compute me interaction



      compute_me(i,j,rsq,rij,fmi,spi,spj);

      if (rsq <= local_cut2) {

	compute_me(i,j,rsq,eij,fmi,spi,spj);

	if (lattice_flag) {

	  compute_me_mech(i,j,fi,spi,spj);

	}

      }



      f[i][0] += fi[0];	 

      f[i][1] += fi[1];
@@ -305,6 +324,7 @@ void PairSpinMe::compute_single_pair(int ii, double fmi[3]) 
  int *type = atom->type;

  double **x = atom->x;

  double **sp = atom->sp;

  double local_cut2;



  double xi[3], rij[3], eij[3];

  double spi[3], spj[3];
@@ -313,7 +333,7 @@ void PairSpinMe::compute_single_pair(int ii, double fmi[3]) 
  int i,j,jj,inum,jnum,itype,jtype;

  int *ilist,*jlist,*numneigh,**firstneigh;



  double rsq, rd, inorm;

  double rsq, inorm;



  inum = list->inum;

  ilist = list->ilist;
@@ -321,6 +341,7 @@ void PairSpinMe::compute_single_pair(int ii, double fmi[3]) 
  firstneigh = list->firstneigh;



  i = ilist[ii];

  itype = type[i];



  spi[0] = sp[i][0];

  spi[1] = sp[i][1];
@@ -339,8 +360,8 @@ void PairSpinMe::compute_single_pair(int ii, double fmi[3]) 


    j = jlist[jj];

    j &= NEIGHMASK;

    itype = type[ii];

    jtype = type[j];

    local_cut2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];



    spj[0] = sp[j][0];

    spj[1] = sp[j][1];
@@ -355,8 +376,10 @@ void PairSpinMe::compute_single_pair(int ii, double fmi[3]) 
    eij[1] = inorm*rij[1];

    eij[2] = inorm*rij[2];



    if (rsq <= local_cut2) {

      compute_me(i,j,rsq,eij,fmi,spi,spj);

    }

  }



}

src/SPIN/pair_spin_me.h

+3 −0
Changes for src/SPIN/pair_spin_me.h: 3 added lines, 0 removed lines.
Original line number Diff line number Diff line
@@ -52,6 +52,9 @@ class PairSpinMe : public PairSpin { 
  double **v_mex, **v_mey, **v_mez;	// me direction

  double **cut_spin_me;			// me cutoff distance 



  int lattice_flag;                     // flag for mech force computation

  class FixNVESpin *lockfixnvespin;     // ptr to FixNVESpin for setups



  void allocate();

};

CRA Git | Maintained and supported by SUSTech CRA and CCSE