Commit f80c527b authored by julient31's avatar julient31
Browse files

Commit JT 100619

- modified precession and Langevin/spin
parent 816546d0
Loading
Loading
Loading
Loading
+2 −1
Original line number Diff line number Diff line
@@ -104,7 +104,8 @@ void ComputeSpin::compute_vector()
        mag[0] += sp[i][0];
        mag[1] += sp[i][1];
        mag[2] += sp[i][2];
        magenergy -= 2.0*(sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
        // magenergy -= 2.0*(sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
        magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
        tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
        ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
        tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
+14 −12
Original line number Diff line number Diff line
@@ -30,7 +30,8 @@
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "random_park.h"
// #include "random_park.h"
#include "random_mars.h"
#include "respa.h"
#include "update.h"
#include "utils.h"
@@ -74,7 +75,8 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :

  // initialize Marsaglia RNG with processor-unique seed

  random = new RanPark(lmp,seed + comm->me);
  // random = new RanPark(lmp,seed + comm->me);
  random = new RanMars(lmp,seed + comm->me);

}

@@ -82,8 +84,6 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :

FixLangevinSpin::~FixLangevinSpin()
{
  memory->destroy(spi);
  memory->destroy(fmi);
  delete random;
}

@@ -113,15 +113,14 @@ void FixLangevinSpin::init()
  }
  if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin has to come after all other spin fixes");

  memory->create(spi,3,"langevin:spi");
  memory->create(fmi,3,"langevin:fmi");

  gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
  dts = update->dt;
  dts = 0.25 * update->dt;

  double hbar = force->hplanck/MY_2PI;  // eV/(rad.THz)
  double kb = force->boltz;             // eV/K
  D = (MY_2PI*alpha_t*gil_factor*kb*temp);
  // D = (MY_2PI*alpha_t*gil_factor*kb*temp);
  D = (1.0/MY_2PI)*(MY_2PI*alpha_t*gil_factor*kb*temp);
  // D = (12.0/MY_2PI)*(MY_2PI*alpha_t*gil_factor*kb*temp);
  D /= (hbar*dts);
  sigma = sqrt(2.0*D);
}
@@ -157,9 +156,12 @@ void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3])
void FixLangevinSpin::add_temperature(double fmi[3])
{

  double rx = sigma*(2.0*random->uniform() - 1.0);
  double ry = sigma*(2.0*random->uniform() - 1.0);
  double rz = sigma*(2.0*random->uniform() - 1.0);
  // double rx = sigma*(2.0*random->uniform() - 1.0);
  // double ry = sigma*(2.0*random->uniform() - 1.0);
  // double rz = sigma*(2.0*random->uniform() - 1.0);
  double rx = sigma*random->gaussian();
  double ry = sigma*random->gaussian();
  double rz = sigma*random->gaussian();

  // adding the random field

+5 −5
Original line number Diff line number Diff line
@@ -32,12 +32,11 @@ class FixLangevinSpin : public Fix {
  void init();
  void setup(int);
  void post_force_respa(int, int, int);
  void add_tdamping(double spi[3], double fmi[3]);      // add transverse damping
  void add_temperature(double fmi[3]);                  // add temperature
  void add_tdamping(double *, double *);                 // add transverse damping
  void add_temperature(double *);                        // add temperature
  int tdamp_flag, ldamp_flag, temp_flag;                 // damping and temperature flags

 protected:
  double *spi, *fmi;
  double alpha_t;               // transverse mag. damping
  double dts;                   // magnetic timestep
  double temp;                  // spin bath temperature
@@ -48,7 +47,8 @@ class FixLangevinSpin : public Fix {
  class Compute *temperature;

  int nlevels_respa;
  class RanPark *random;
  // class RanPark *random;
  class RanMars *random;
  int seed;

};
+8 −4
Original line number Diff line number Diff line
@@ -257,7 +257,8 @@ void FixPrecessionSpin::post_force(int /* vflag */)

      if (zeeman_flag) {          // compute Zeeman interaction
        compute_zeeman(i,fmi);
        epreci -= 2.0*hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
        // epreci -= 2.0*hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
        epreci -= hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
      }

      if (aniso_flag) {           // compute magnetic anisotropy
@@ -295,9 +296,12 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f
void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
{
  double **sp = atom->sp;
  fmi[0] += 0.5*sp[i][3]*hx;
  fmi[1] += 0.5*sp[i][3]*hy;
  fmi[2] += 0.5*sp[i][3]*hz;
  // fmi[0] += 0.5*sp[i][3]*hx;
  // fmi[1] += 0.5*sp[i][3]*hy;
  // fmi[2] += 0.5*sp[i][3]*hz;
  fmi[0] += sp[i][3]*hx;
  fmi[1] += sp[i][3]*hy;
  fmi[2] += sp[i][3]*hz;
}

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