Commit cc698f0b authored by alxvov's avatar alxvov
Browse files

Merge branch 'OSO' into OSO_LBFGS

parents 61b12a09 5c8e8124
Loading
Loading
Loading
Loading
+1 −1
Changes for src/SPIN/fix_nve_spin.h: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -54,7 +54,7 @@ friend class PairSpin;

  double dtv, dtf, dts;			// velocity, force, and spin timesteps

  int nlocal_max;			// max value of nlocal (for lists size)
  int nlocal_max;			// max value of nlocal (for size of lists)

  int pair_spin_flag;			// magnetic pair flags
  int long_spin_flag;			// magnetic long-range flag
+54 −32
Changes for src/SPIN/min_spin_oso_cg.cpp: 54 added lines, 32 removed lines.
Original line number Diff line number Diff line
@@ -34,6 +34,7 @@
#include "output.h"
#include "timer.h"
#include "error.h"
#include "memory.h"
#include "modify.h"
#include "math_special.h"
#include "math_const.h"
@@ -60,8 +61,20 @@ static const char cite_minstyle_spin_oso_cg[] =

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

MinSpinOSO_CG::MinSpinOSO_CG(LAMMPS *lmp) : Min(lmp) {
MinSpinOSO_CG::MinSpinOSO_CG(LAMMPS *lmp) : 
  Min(lmp), g_old(NULL), g_cur(NULL), p_s(NULL)
{
  if (lmp->citeme) lmp->citeme->add(cite_minstyle_spin_oso_cg);
  nlocal_max = 0;
}

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

MinSpinOSO_CG::~MinSpinOSO_CG()
{
  memory->destroy(g_old);
  memory->destroy(g_cur);
  memory->destroy(p_s);
}

/* ---------------------------------------------------------------------- */
@@ -75,6 +88,13 @@ void MinSpinOSO_CG::init()

  dts = dt = update->dt;
  last_negative = update->ntimestep;
  
  // allocate tables

  nlocal_max = atom->nlocal;
  memory->grow(g_old,3*nlocal_max,"min/spin/oso/cg:g_old");
  memory->grow(g_cur,3*nlocal_max,"min/spin/oso/cg:g_cur");
  memory->grow(p_s,3*nlocal_max,"min/spin/oso/cg:p_s");
}

/* ---------------------------------------------------------------------- */
@@ -134,17 +154,21 @@ void MinSpinOSO_CG::reset_vectors()
   minimization via damped spin dynamics
------------------------------------------------------------------------- */

// g_old g_cur p_s

int MinSpinOSO_CG::iterate(int maxiter)
{
  int nlocal = atom->nlocal;
  bigint ntimestep;
  double fmdotfm;
  int flag, flagall;

  // not sure it is best place to allocate memory
  int nlocal = atom->nlocal;
  g = (double *) calloc(3*nlocal, sizeof(double));
  p = (double *) calloc(3*nlocal, sizeof(double));
  g_old = (double *) calloc(3*nlocal, sizeof(double));
  if (nlocal_max < nlocal) {
    nlocal_max = nlocal;
    memory->grow(g_old,3*nlocal_max,"min/spin/oso/cg:g_old");
    memory->grow(g_cur,3*nlocal_max,"min/spin/oso/cg:g_cur");
    memory->grow(p_s,3*nlocal_max,"min/spin/oso/cg:p_s");
  }

  for (int iter = 0; iter < maxiter; iter++) {
    
@@ -211,10 +235,6 @@ int MinSpinOSO_CG::iterate(int maxiter)
    }
  }

  free(p);
  free(g);
  free(g_old);

  return MAXITER;
}

@@ -286,9 +306,9 @@ void MinSpinOSO_CG::calc_gradient(double dts)

    // calculate gradients
    
    g[3 * i + 0] = -tdampz * dts;
    g[3 * i + 1] = tdampy * dts;
    g[3 * i + 2] = -tdampx * dts;
    g_cur[3 * i + 0] = -tdampz * dts;
    g_cur[3 * i + 1] = tdampy * dts;
    g_cur[3 * i + 2] = -tdampx * dts;
  }
}

@@ -312,15 +332,15 @@ void MinSpinOSO_CG::calc_search_direction(int iter)
  if (iter <= 2 || iter % 5 == 0){ 	// steepest descent direction
    for (int i = 0; i < nlocal; i++) {
      for (int j = 0; j < 3; j++){
          p[3 * i + j] = -g[3 * i + j];
          g_old[3 * i + j] = g[3 * i + j];
          p_s[3 * i + j] = -g_cur[3 * i + j];
          g_old[3 * i + j] = g_cur[3 * i + j];
      }
    }
  } else { 				// conjugate direction
    for (int i = 0; i < nlocal; i++) {
      for (int j = 0; j < 3; j++){
	g2old += g_old[3 * i + j] * g_old[3 * i + j];
	g2 += g[3 * i + j] * g[3 * i + j];
	g2 += g_cur[3 * i + j] * g_cur[3 * i + j];
      }
    }

@@ -337,8 +357,8 @@ void MinSpinOSO_CG::calc_search_direction(int iter)
    
    for (int i = 0; i < nlocal; i++) {
      for (int j = 0; j < 3; j++){
	p[3 * i + j] = beta * p[3 * i + j] - g[3 * i + j];
	g_old[3 * i + j] = g[3 * i + j];
	p_s[3 * i + j] = beta * p_s[3 * i + j] - g_cur[3 * i + j];
	g_old[3 * i + j] = g_cur[3 * i + j];
      }
    }
  }
@@ -360,7 +380,7 @@ void MinSpinOSO_CG::advance_spins()
  // loop on all spins on proc.

  for (int i = 0; i < nlocal; i++) {
    rodrigues_rotation(p + 3 * i, rot_mat);
    rodrigues_rotation(p_s + 3 * i, rot_mat);
    
    // rotate spins
    
@@ -418,6 +438,8 @@ double MinSpinOSO_CG::fmnorm_sqr()

void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out)
{
  double theta,A,B,D,x,y,z;
  double s1,s2,s3,a1,a2,a3;

  if (fabs(upp_tr[0]) < 1.0e-40 &&
      fabs(upp_tr[1]) < 1.0e-40 &&
@@ -433,16 +455,16 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out)
    return;
  }

  double theta = sqrt(upp_tr[0] * upp_tr[0] +
  theta = sqrt(upp_tr[0] * upp_tr[0] +
               upp_tr[1] * upp_tr[1] +
               upp_tr[2] * upp_tr[2]);

  double A = cos(theta);
  double B = sin(theta);
  double D = 1 - A;
  double x = upp_tr[0]/theta;
  double y = upp_tr[1]/theta;
  double z = upp_tr[2]/theta;
  A = cos(theta);
  B = sin(theta);
  D = 1 - A;
  x = upp_tr[0]/theta;
  y = upp_tr[1]/theta;
  z = upp_tr[2]/theta;

  // diagonal elements of U
  
@@ -452,13 +474,13 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out)

  // off diagonal of U
  
  double s1 = -y * z *D;
  double s2 = x * z * D;
  double s3 = -x * y * D;
  s1 = -y * z *D;
  s2 = x * z * D;
  s3 = -x * y * D;

  double a1 = x * B;
  double a2 = y * B;
  double a3 = z * B;
  a1 = x * B;
  a2 = y * B;
  a3 = z * B;

  out[1] = s1 + a1;
  out[3] = s1 - a1;
+6 −5
Changes for src/SPIN/min_spin_oso_cg.h: 6 added lines, 5 removed lines.
Original line number Diff line number Diff line
@@ -28,7 +28,7 @@ class MinSpinOSO_CG : public Min {

public:
    MinSpinOSO_CG(class LAMMPS *);
    ~MinSpinOSO_CG() {}
    virtual ~MinSpinOSO_CG();
    void init();
    void setup_style();
    int modify_param(int, char **);
@@ -43,6 +43,7 @@ public:
private:
    // global and spin timesteps

    int nlocal_max;		// max value of nlocal (for size of lists)
    double dt;
    double dts;

@@ -53,11 +54,11 @@ private:
    double *fmvec;		// variables for atomic dof, as 1d vector

    double *g_old;  		// gradient vector at previous iteration
    double *g;  		// gradient vector
    double *p;  		// search direction vector
    double *g_cur;  		// current gradient vector
    double *p_s;  		// search direction vector

    void vm3(const double *m, const double *v, double *out);
    void rodrigues_rotation(const double *upp_tr, double *out);
    void vm3(const double *, const double *, double *);
    void rodrigues_rotation(const double *, double *);

    bigint last_negative;
};