Commit 31d2b23f authored by alxvov's avatar alxvov
Browse files

rename cg2 -> cg

parent fabe611c

src/SPIN/min_spin_oso_cg.cpp

+266 −100
Original line number Diff line number Diff line
@@ -38,6 +38,7 @@ 
#include "modify.h"

#include "math_special.h"

#include "math_const.h"

#include "universe.h"



using namespace LAMMPS_NS;

using namespace MathConst;
@@ -66,7 +67,13 @@ MinSpinOSO_CG::MinSpinOSO_CG(LAMMPS *lmp) : 
{

  if (lmp->citeme) lmp->citeme->add(cite_minstyle_spin_oso_cg);

  nlocal_max = 0;

  alpha_damp = 1.0;



  // nreplica = number of partitions

  // ireplica = which world I am in universe



  nreplica = universe->nworlds;

  ireplica = universe->iworld;

  use_line_search = 1;

  discrete_factor = 10.0;

}
@@ -77,15 +84,20 @@ MinSpinOSO_CG::~MinSpinOSO_CG() 
  memory->destroy(g_old);

  memory->destroy(g_cur);

  memory->destroy(p_s);

  if (use_line_search)

    memory->destroy(sp_copy);

}



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



void MinSpinOSO_CG::init()

{



  local_iter = 0;

  der_e_cur = 0.0;

  der_e_pr = 0.0;



  Min::init();



  dts = dt = update->dt;

  last_negative = update->ntimestep;
@@ -95,6 +107,8 @@ void MinSpinOSO_CG::init() 
  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");

  if (use_line_search)

    memory->grow(sp_copy,nlocal_max,3,"min/spin/oso/cg:sp_copy");

}



/* ---------------------------------------------------------------------- */
@@ -117,9 +131,9 @@ void MinSpinOSO_CG::setup_style() 


int MinSpinOSO_CG::modify_param(int narg, char **arg)

{

  if (strcmp(arg[0],"alpha_damp") == 0) {

  if (strcmp(arg[0],"line_search") == 0) {

    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");

    alpha_damp = force->numeric(FLERR,arg[1]);

    use_line_search = force->numeric(FLERR,arg[1]);

    return 2;

  }

  if (strcmp(arg[0],"discrete_factor") == 0) {
@@ -160,15 +174,18 @@ int MinSpinOSO_CG::iterate(int maxiter) 
  bigint ntimestep;

  double fmdotfm;

  int flag, flagall;



  // grow tables if nlocal increased

  double **sp = atom->sp;

  double der_e_cur_tmp = 0.0;



  if (nlocal_max < nlocal) {

    nlocal_max = nlocal;

    local_iter = 0;

    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");

    if (use_line_search)

      memory->grow(sp_copy,nlocal_max,3,"min/spin/oso/cg:sp_copy");

  }



  for (int iter = 0; iter < maxiter; iter++) {
@@ -182,16 +199,51 @@ int MinSpinOSO_CG::iterate(int maxiter) 
    // optimize timestep accross processes / replicas

    // need a force calculation for timestep optimization



    if (local_iter == 0) energy_force(0);

    dts = evaluate_dt();

    if (use_line_search) {



      // here we need to do line search

      if (local_iter == 0)

        calc_gradient();



    calc_gradient(dts);

      calc_search_direction();

    advance_spins();

      der_e_cur = 0.0;

      for (int i = 0; i < 3 * nlocal; i++)

        der_e_cur += g_cur[i] * p_s[i];

      MPI_Allreduce(&der_e_cur,&der_e_cur_tmp,1,MPI_DOUBLE,MPI_SUM,world);

      der_e_cur = der_e_cur_tmp;

      if (update->multireplica == 1) {

        MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);

      }

      for (int i = 0; i < nlocal; i++)

        for (int j = 0; j < 3; j++)

      sp_copy[i][j] = sp[i][j];



      eprevious = ecurrent;

      der_e_pr = der_e_cur;

      calc_and_make_step(0.0, 1.0, 0);

    }

    else{



      // here we don't do line search

      // but use cutoff rotation angle

      // if gneb calc., nreplica > 1

      // then calculate gradients and advance spins

      // of intermediate replicas only



      if (nreplica > 1) {

      if(ireplica != 0 && ireplica != nreplica-1)

      calc_gradient();

      calc_search_direction();

      advance_spins();

      } else{

      calc_gradient();

      calc_search_direction();

      advance_spins();

      }

      eprevious = ecurrent;

      ecurrent = energy_force(0);

      neval++;

    }



    //// energy tolerance criterion

    //// only check after DELAYSTEP elapsed since velocties reset to 0
@@ -239,77 +291,28 @@ int MinSpinOSO_CG::iterate(int maxiter) 
  return MAXITER;

}



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

   evaluate max timestep

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



double MinSpinOSO_CG::evaluate_dt()

{

  double dtmax;

  double fmsq;

  double fmaxsqone,fmaxsqloc,fmaxsqall;

  int nlocal = atom->nlocal;

  double **fm = atom->fm;



  // finding max fm on this proc.



  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;

  for (int i = 0; i < nlocal; i++) {

    fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];

    fmaxsqone = MAX(fmaxsqone,fmsq);

  }



  // finding max fm on this replica



  fmaxsqloc = fmaxsqone;

  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);



  // finding max fm over all replicas, if necessary

  // this communicator would be invalid for multiprocess replicas



  fmaxsqall = fmaxsqloc;

  if (update->multireplica == 1) {

    fmaxsqall = fmaxsqloc;

    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);

  }



  if (fmaxsqall == 0.0)

    error->all(FLERR,"Incorrect fmaxsqall calculation");



  // define max timestep by dividing by the

  // inverse of max frequency by discrete_factor



  dtmax = MY_2PI/(discrete_factor*sqrt(fmaxsqall));



  return dtmax;

}



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

   calculate gradients

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



void MinSpinOSO_CG::calc_gradient(double dts)

void MinSpinOSO_CG::calc_gradient()

{

  int nlocal = atom->nlocal;

  double **sp = atom->sp;

  double **fm = atom->fm;

  double tdampx, tdampy, tdampz;

  double hbar = force->hplanck/MY_2PI;

  double factor;



  if (use_line_search)

    factor = hbar;

  else factor = evaluate_dt();



  // loop on all spins on proc.



  for (int i = 0; i < nlocal; i++) {

    

    // calc. damping torque

    

    tdampx = -alpha_damp*(fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);

    tdampy = -alpha_damp*(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);

    tdampz = -alpha_damp*(fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);



    // calculate gradients

    

    g_cur[3 * i + 0] = -tdampz * dts;

    g_cur[3 * i + 1] = tdampy * dts;

    g_cur[3 * i + 2] = -tdampx * dts;

    g_cur[3 * i + 0] = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]) * factor;

    g_cur[3 * i + 1] = -(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]) * factor;

    g_cur[3 * i + 2] = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]) * factor;

  }

}
@@ -329,6 +332,7 @@ void MinSpinOSO_CG::calc_search_direction() 


  double g2_global = 0.0;

  double g2old_global = 0.0;



  if (local_iter == 0 || local_iter % 5 == 0){ 	// steepest descent direction

    for (int i = 0; i < 3 * nlocal; i++) {

      p_s[i] = -g_cur[i];
@@ -341,7 +345,6 @@ void MinSpinOSO_CG::calc_search_direction() 
    }



    // now we need to collect/broadcast beta on this replica

    // different replica can have different beta for now.

    // need to check what is beta for GNEB



    MPI_Allreduce(&g2,&g2_global,1,MPI_DOUBLE,MPI_SUM,world);
@@ -354,12 +357,11 @@ void MinSpinOSO_CG::calc_search_direction() 
      MPI_Allreduce(&g2,&g2_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);

      MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);

    }



    if (fabs(g2_global) < 1.0e-60) beta = 0.0;

    else beta = g2_global / g2old_global;

    // calculate conjugate direction

    for (int i = 0; i < 3 * nlocal; i++) {

  	  p_s[i] = beta * p_s[i] - g_cur[i];

      p_s[i] = (beta * p_s[i] - g_cur[i]);

      g_old[i] = g_cur[i];

    }

  }
@@ -400,8 +402,12 @@ double MinSpinOSO_CG::max_torque() 
{

  double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall;

  int nlocal = atom->nlocal;

  double factor;

  double hbar = force->hplanck/MY_2PI;



  if (use_line_search) factor = 1.0;

  else factor = hbar;



  // finding max fm on this proc.



  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
@@ -426,7 +432,7 @@ double MinSpinOSO_CG::max_torque() 
    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);

  }



  return sqrt(fmaxsqall) * hbar;

  return sqrt(fmaxsqall) * factor;

}



/* ----------------------------------------------------------------------
@@ -453,8 +459,10 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out) 
  // if upp_tr is zero, return unity matrix

  for(int k = 0; k < 3; k++){

    for(int m = 0; m < 3; m++){

	if (m == k) out[3 * k + m] = 1.0;

	else out[3 * k + m] = 0.0;

      if (m == k)

        out[3 * k + m] = 1.0;

      else

        out[3 * k + m] = 0.0;

    }

  }

    return;
@@ -505,8 +513,166 @@ void MinSpinOSO_CG::vm3(const double *m, const double *v, double *out) 
{

  for(int i = 0; i < 3; i++){

    out[i] *= 0.0;

    for(int j = 0; j < 3; j++){

    for(int j = 0; j < 3; j++)

    out[i] += *(m + 3 * j + i) * v[j];

  }

}





void MinSpinOSO_CG::make_step(double c, double *energy_and_der)

{

  double p_scaled[3];

  int nlocal = atom->nlocal;

  double rot_mat[9]; // exponential of matrix made of search direction

  double s_new[3];

  double **sp = atom->sp;

  double der_e_cur_tmp = 0.0;;



  for (int i = 0; i < nlocal; i++) {



    // scale the search direction



    for (int j = 0; j < 3; j++) p_scaled[j] = c * p_s[3 * i + j];



    // calculate rotation matrix



    rodrigues_rotation(p_scaled, rot_mat);



    // rotate spins



    vm3(rot_mat, sp[i], s_new);

    for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];

  }



  ecurrent = energy_force(0);

  calc_gradient();

  neval++;

  der_e_cur = 0.0;

  for (int i = 0; i < 3 * nlocal; i++) {

    der_e_cur += g_cur[i] * p_s[i];

  }

  MPI_Allreduce(&der_e_cur,&der_e_cur_tmp,1,MPI_DOUBLE,MPI_SUM,world);

  der_e_cur = der_e_cur_tmp;

  if (update->multireplica == 1) {

    MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);

  }



  energy_and_der[0] = ecurrent;

  energy_and_der[1] = der_e_cur;

}



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

  Calculate step length which satisfies approximate Wolfe conditions

  using the cubic interpolation

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



int MinSpinOSO_CG::calc_and_make_step(double a, double b, int index)

{

  double e_and_d[2] = {0.0,0.0};

  double alpha,c1,c2,c3;

  double **sp = atom->sp;

  int nlocal = atom->nlocal;



  make_step(b,e_and_d);

  ecurrent = e_and_d[0];

  der_e_cur = e_and_d[1];

  index++;



  if (awc(der_e_pr,eprevious,e_and_d[1],e_and_d[0]) || index == 10){

    MPI_Bcast(&b,1,MPI_DOUBLE,0,world);

    for (int i = 0; i < 3 * nlocal; i++) {

      p_s[i] = b * p_s[i];

    }

    return 1;

  }

  else{

    double r,f0,f1,df0,df1;

    r = b - a;

    f0 = eprevious;

    f1 = ecurrent;

    df0 = der_e_pr;

    df1 = der_e_cur;



    c1 = -2.0*(f1-f0)/(r*r*r)+(df1+df0)/(r*r);

    c2 = 3.0*(f1-f0)/(r*r)-(df1+2.0*df0)/(r);

    c3 = df0;



    // f(x) = c1 x^3 + c2 x^2 + c3 x^1 + c4

    // has minimum at alpha below. We do not check boundaries.



    alpha = (-c2 + sqrt(c2*c2 - 3.0*c1*c3))/(3.0*c1);

    MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world);



    if (alpha < 0.0) alpha = r/2.0;



    for (int i = 0; i < nlocal; i++) {

      for (int j = 0; j < 3; j++) sp[i][j] = sp_copy[i][j];

    }

    calc_and_make_step(0.0, alpha, index);

   }



  return 0;

}



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

  Approximate Wolfe conditions:

  William W. Hager and Hongchao Zhang

  SIAM J. optim., 16(1), 170-192. (23 pages)

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



int MinSpinOSO_CG::awc(double der_phi_0, double phi_0, double der_phi_j, double phi_j){



  double eps = 1.0e-6;

  double delta = 0.1;

  double sigma = 0.9;



  if ((phi_j<=phi_0+eps*fabs(phi_0)) && ((2.0*delta-1.0) * der_phi_0>=der_phi_j>=sigma*der_phi_0))

    return 1;

  else

    return 0;

}



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

   evaluate max timestep

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



double MinSpinOSO_CG::evaluate_dt()

{

  double dtmax;

  double fmsq;

  double fmaxsqone,fmaxsqloc,fmaxsqall;

  int nlocal = atom->nlocal;

  double **fm = atom->fm;



  // finding max fm on this proc.



  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;

  for (int i = 0; i < nlocal; i++) {

    fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];

    fmaxsqone = MAX(fmaxsqone,fmsq);

  }



  // finding max fm on this replica



  fmaxsqloc = fmaxsqone;

  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);



  // finding max fm over all replicas, if necessary

  // this communicator would be invalid for multiprocess replicas



  fmaxsqall = fmaxsqloc;

  if (update->multireplica == 1) {

    fmaxsqall = fmaxsqloc;

    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);

  }



  if (fmaxsqall == 0.0)

    error->all(FLERR,"Incorrect fmaxsqall calculation");



  // define max timestep by dividing by the

  // inverse of max frequency by discrete_factor



  dtmax = MY_2PI/(discrete_factor*sqrt(fmaxsqall));



  return dtmax;

}
 
 No newline at end of file

src/SPIN/min_spin_oso_cg.h

+22 −21
Original line number Diff line number Diff line
@@ -33,33 +33,34 @@ class MinSpinOSO_CG : public Min { 
    int modify_param(int, char **);

    void reset_vectors();

    int iterate(int);



  private:

    double evaluate_dt();

    void advance_spins();

    double max_torque();

    void calc_gradient(double);

    void calc_search_direction();



    // global and spin timesteps



    double dt;

    double dts;

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



    double alpha_damp;		// damping for spin minimization

    double discrete_factor;	// factor for spin timestep evaluation



    double dt;        // global timestep

    double dts;       // spin timestep

    int ireplica,nreplica; // for neb

    double *spvec;		// variables for atomic dof, as 1d vector

    double *fmvec;		// variables for atomic dof, as 1d vector



    double *g_old;  		// gradient vector at previous iteration

    double *g_cur;  	// current gradient vector

    double *g_old;  	// gradient vector at previous step

    double *p_s;  		// search direction vector

    int local_iter;  // number of times we call search_direction

    double **sp_copy;   // copy of the spins

    int local_iter;     // for neb

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

    double discrete_factor;	// factor for spin timestep evaluation



    double evaluate_dt();

    void advance_spins();

    void calc_gradient();

    void calc_search_direction();

    double maximum_rotation(double *);

    void vm3(const double *, const double *, double *);

    void rodrigues_rotation(const double *, double *);

    int calc_and_make_step(double, double, int);

    int awc(double, double, double, double);

    void make_step(double, double *);

    double max_torque();

    double der_e_cur;   // current derivative along search dir.

    double der_e_pr;    // previous derivative along search dir.

    int use_line_search; // use line search or not.



    bigint last_negative;

};

src/SPIN/min_spin_oso_cg2.cpp

deleted100644 → 0
+0 −679

File deleted.

Preview size limit exceeded, changes collapsed.

src/SPIN/min_spin_oso_cg2.h

deleted100644 → 0
+0 −72
Original line number Diff line number Diff line 
/* -*- c++ -*- ----------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov



   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.



   See the README file in the top-level LAMMPS directory.

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



#ifdef MINIMIZE_CLASS



MinimizeStyle(spin/oso_cg2, MinSpinOSO_CG2)



#else



#ifndef LMP_MIN_SPIN_OSO_CG2_H

#define LMP_MIN_SPIN_OSO_CG2_H



#include "min.h"



namespace LAMMPS_NS {



class MinSpinOSO_CG2: public Min {

  public:

    MinSpinOSO_CG2(class LAMMPS *);

    virtual ~MinSpinOSO_CG2();

    void init();

    void setup_style();

    int modify_param(int, char **);

    void reset_vectors();

    int iterate(int);

  private:

    double dt;        // global timestep

    double dts;       // spin timestep

    int ireplica,nreplica; // for neb

    double *spvec;		// variables for atomic dof, as 1d vector

    double *fmvec;		// variables for atomic dof, as 1d vector

    double *g_cur;  	// current gradient vector

    double *g_old;  	// gradient vector at previous step

    double *p_s;  		// search direction vector

    double **sp_copy;   // copy of the spins

    int local_iter;     // for neb

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

    double discrete_factor;	// factor for spin timestep evaluation



    double evaluate_dt();

    void advance_spins();

    void calc_gradient();

    void calc_search_direction();

    double maximum_rotation(double *);

    void vm3(const double *, const double *, double *);

    void rodrigues_rotation(const double *, double *);

    int calc_and_make_step(double, double, int);

    int awc(double, double, double, double);

    void make_step(double, double *);

    double max_torque();

    double der_e_cur;   // current derivative along search dir.

    double der_e_pr;    // previous derivative along search dir.

    int use_line_search; // use line search or not.

    double maxepsrot;



    bigint last_negative;

};



}



#endif

#endif

CRA Git | Maintained and supported by SUSTech CRA and CCSE