Commit 3d7b3877 authored by julient31's avatar julient31
Browse files

Commit JT 040919 

- first working version of setforce
- modified fix nve/spin to account for it
parent 8ec231a8

src/SPIN/fix_nve_spin.cpp

+37 −18
Original line number Diff line number Diff line
@@ -31,9 +31,10 @@ 
#include "comm.h"

#include "domain.h"

#include "error.h"

#include "fix_precession_spin.h"

#include "fix_nve_spin.h"

#include "fix_precession_spin.h"

#include "fix_langevin_spin.h"

#include "fix_setforce_spin.h"

#include "force.h"

#include "math_vector.h"

#include "math_extra.h"
@@ -131,6 +132,7 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : 
  precession_spin_flag = 0;

  maglangevin_flag = 0;

  tdamp_flag = temp_flag = 0;

  setforce_spin_flag = 0;



}
@@ -237,6 +239,15 @@ void FixNVESpin::init() 
   if (locklangevinspin->temp_flag == 1) temp_flag = 1;

  }



  // ptrs FixSetForceSpin classes



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

    if (strstr(modify->fix[iforce]->style,"setforce/spin")) {

      setforce_spin_flag = 1;

      locksetforcespin = (FixSetForceSpin *) modify->fix[iforce];

    }

  }

  

  // setting the sector variables/lists



  nsectors = 0;
@@ -458,6 +469,14 @@ void FixNVESpin::ComputeInteractionsSpin(int i) 
    }

  }



  // update setforce of magnetic interactions



  if (setforce_spin_flag) {

    locksetforcespin->single_setforce_spin(i,fmi);

  }



  //printf("test after setforce: %g %g %g \n",fmi[0],fmi[1],fmi[2]);



  // replace the magnetic force fm[i] by its new value fmi



  fm[i][0] = fmi[0];
@@ -577,9 +596,9 @@ void FixNVESpin::AdvanceSingleSpin(int i) 
  g[1] += (fm[i][1]*energy - 0.5*sp[i][1]*fm2)*0.5*dts2;

  g[2] += (fm[i][2]*energy - 0.5*sp[i][2]*fm2)*0.5*dts2;



  g[0] /= (1+0.25*fm2*dts2);

  g[1] /= (1+0.25*fm2*dts2);

  g[2] /= (1+0.25*fm2*dts2);

  g[0] /= (1.0 + 0.25*fm2*dts2);

  g[1] /= (1.0 + 0.25*fm2*dts2);

  g[2] /= (1.0 + 0.25*fm2*dts2);



  sp[i][0] = g[0];

  sp[i][1] = g[1];
@@ -587,11 +606,11 @@ void FixNVESpin::AdvanceSingleSpin(int i) 


  // renormalization (check if necessary)



  msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];

  scale = 1.0/sqrt(msq);

  sp[i][0] *= scale;

  sp[i][1] *= scale;

  sp[i][2] *= scale;

  //msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];

  //scale = 1.0/sqrt(msq);

  //sp[i][0] *= scale;

  //sp[i][1] *= scale;

  //sp[i][2] *= scale;



  // comm. sp[i] to atoms with same tag (for serial algo)

src/SPIN/fix_nve_spin.h

+3 −1
Original line number Diff line number Diff line
@@ -61,11 +61,13 @@ friend class PairSpin; 
  int precession_spin_flag;		// magnetic precession flags

  int maglangevin_flag;			// magnetic langevin flags

  int tdamp_flag, temp_flag;

  int setforce_spin_flag;



  // pointers to magnetic fixes



  class FixPrecessionSpin *lockprecessionspin;

  class FixLangevinSpin *locklangevinspin;

  class FixSetForceSpin *locksetforcespin;



  // pointers to magnetic pair styles

src/SPIN/fix_setforce_spin.cpp

+90 −219
Original line number Diff line number Diff line
@@ -11,9 +11,18 @@ 
   See the README file in the top-level LAMMPS directory.

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



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

   Contributing authors: Julien Tranchida (SNL)

                         Aidan Thompson (SNL)



   Please cite the related publication:

   Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018).

   Massively parallel symplectic algorithm for coupled magnetic spin dynamics

   and molecular dynamics. Journal of Computational Physics.

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



#include <cstring>

#include <cstdlib>

//#include "fix_setforce.h"

#include "fix_setforce_spin.h"

#include "atom.h"

#include "update.h"
@@ -35,196 +44,13 @@ enum{NONE,CONSTANT,EQUAL,ATOM}; 
/* ---------------------------------------------------------------------- */



FixSetForceSpin::FixSetForceSpin(LAMMPS *lmp, int narg, char **arg) :

  //FixSetForce(lmp, narg, arg),

  Fix(lmp, narg, arg),

  xstr(NULL), ystr(NULL), zstr(NULL), idregion(NULL), sforce(NULL)

{

  if (narg < 6) error->all(FLERR,"Illegal fix setforce spin command");



  dynamic_group_allow = 1;

  vector_flag = 1;

  size_vector = 3;

  global_freq = 1;

  extvector = 1;

  respa_level_support = 1;

  ilevel_respa = nlevels_respa = 0;

  xstr = ystr = zstr = NULL;



  if (strstr(arg[3],"v_") == arg[3]) {

    int n = strlen(&arg[3][2]) + 1;

    xstr = new char[n];

    strcpy(xstr,&arg[3][2]);

  } else if (strcmp(arg[3],"NULL") == 0) {

    xstyle = NONE;

  } else {

    xvalue = force->numeric(FLERR,arg[3]);

    xstyle = CONSTANT;

  }

  if (strstr(arg[4],"v_") == arg[4]) {

    int n = strlen(&arg[4][2]) + 1;

    ystr = new char[n];

    strcpy(ystr,&arg[4][2]);

  } else if (strcmp(arg[4],"NULL") == 0) {

    ystyle = NONE;

  } else {

    yvalue = force->numeric(FLERR,arg[4]);

    ystyle = CONSTANT;

  }

  if (strstr(arg[5],"v_") == arg[5]) {

    int n = strlen(&arg[5][2]) + 1;

    zstr = new char[n];

    strcpy(zstr,&arg[5][2]);

  } else if (strcmp(arg[5],"NULL") == 0) {

    zstyle = NONE;

  } else {

    zvalue = force->numeric(FLERR,arg[5]);

    zstyle = CONSTANT;

  }



  // optional args



  iregion = -1;

  idregion = NULL;



  int iarg = 6;

  while (iarg < narg) {

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix setforce command");

      iregion = domain->find_region(arg[iarg+1]);

      if (iregion == -1)

        error->all(FLERR,"Region ID for fix setforce does not exist");

      int n = strlen(arg[iarg+1]) + 1;

      idregion = new char[n];

      strcpy(idregion,arg[iarg+1]);

      iarg += 2;

    } else error->all(FLERR,"Illegal fix setforce command");

  }



  force_flag = 0;

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



  maxatom = 1;

  memory->create(sforce,maxatom,3,"setforce:sforce");

}



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



FixSetForceSpin::~FixSetForceSpin()

{

  if (copymode) return;



  delete [] xstr;

  delete [] ystr;

  delete [] zstr;

  delete [] idregion;

  memory->destroy(sforce);

}



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



int FixSetForceSpin::setmask()

{

  int mask = 0;

  mask |= POST_FORCE;

  mask |= POST_FORCE_RESPA;

  mask |= MIN_POST_FORCE;

  return mask;

}



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



void FixSetForceSpin::init()

{

  // check variables



  if (xstr) {

    xvar = input->variable->find(xstr);

    if (xvar < 0)

      error->all(FLERR,"Variable name for fix setforce does not exist");

    if (input->variable->equalstyle(xvar)) xstyle = EQUAL;

    else if (input->variable->atomstyle(xvar)) xstyle = ATOM;

    else error->all(FLERR,"Variable for fix setforce is invalid style");

  }

  if (ystr) {

    yvar = input->variable->find(ystr);

    if (yvar < 0)

      error->all(FLERR,"Variable name for fix setforce does not exist");

    if (input->variable->equalstyle(yvar)) ystyle = EQUAL;

    else if (input->variable->atomstyle(yvar)) ystyle = ATOM;

    else error->all(FLERR,"Variable for fix setforce is invalid style");

  }

  if (zstr) {

    zvar = input->variable->find(zstr);

    if (zvar < 0)

      error->all(FLERR,"Variable name for fix setforce does not exist");

    if (input->variable->equalstyle(zvar)) zstyle = EQUAL;

    else if (input->variable->atomstyle(zvar)) zstyle = ATOM;

    else error->all(FLERR,"Variable for fix setforce is invalid style");

  }



  // set index and check validity of region



  if (iregion >= 0) {

    iregion = domain->find_region(idregion);

    if (iregion == -1)

      error->all(FLERR,"Region ID for fix setforce does not exist");

  }



  if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)

    varflag = ATOM;

  else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)

    varflag = EQUAL;

  else varflag = CONSTANT;



  if (strstr(update->integrate_style,"respa")) {

    nlevels_respa = ((Respa *) update->integrate)->nlevels;

    if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);

    else ilevel_respa = nlevels_respa-1;

  }



  // cannot use non-zero forces for a minimization since no energy is integrated

  // use fix addforce instead



  int flag = 0;

  if (update->whichflag == 2) {

    if (xstyle == EQUAL || xstyle == ATOM) flag = 1;

    if (ystyle == EQUAL || ystyle == ATOM) flag = 1;

    if (zstyle == EQUAL || zstyle == ATOM) flag = 1;

    if (xstyle == CONSTANT && xvalue != 0.0) flag = 1;

    if (ystyle == CONSTANT && yvalue != 0.0) flag = 1;

    if (zstyle == CONSTANT && zvalue != 0.0) flag = 1;

  }

  if (flag)

    error->all(FLERR,"Cannot use non-zero forces in an energy minimization");

}



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



void FixSetForceSpin::setup(int vflag)

{

  if (strstr(update->integrate_style,"verlet"))

    post_force(vflag);

  else

    for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {

      ((Respa *) update->integrate)->copy_flevel_f(ilevel);

      post_force_respa(vflag,ilevel,0);

      ((Respa *) update->integrate)->copy_f_flevel(ilevel);

    }

}



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



void FixSetForceSpin::min_setup(int vflag)

{

  post_force(vflag);

}

  FixSetForce(lmp, narg, arg) {}



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



void FixSetForceSpin::post_force(int /*vflag*/)

{

  double **x = atom->x;

  //double **f = atom->f;

  double **fm = atom->fm;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;
@@ -296,6 +122,85 @@ void FixSetForceSpin::post_force(int /*vflag*/) 


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



void FixSetForceSpin::single_setforce_spin(int i, double fmi[3])

{

  double **x = atom->x;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



  //printf("test inside before setforce: %g %g %g \n",fmi[0],fmi[1],fmi[2]);



  // update region if necessary



  Region *region = NULL;

  if (iregion >= 0) {

    region = domain->regions[iregion];

    region->prematch();

  }



  // reallocate sforce array if necessary



  if (varflag == ATOM && atom->nmax > maxatom) {

    maxatom = atom->nmax;

    memory->destroy(sforce);

    memory->create(sforce,maxatom,3,"setforce:sforce");

  }



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

  force_flag = 0;

  

  if (varflag == CONSTANT) {

    if (mask[i] & groupbit) {

      //if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;

      //if (region && region->match(x[i][0],x[i][1],x[i][2])){

      foriginal[0] += fmi[0];

      foriginal[1] += fmi[1];

      foriginal[2] += fmi[2];

      if (xstyle) fmi[0] = xvalue;

      if (ystyle) fmi[1] = yvalue;

      if (zstyle) fmi[2] = zvalue;

      //printf("test inside inter setforce: %g %g %g \n",fmi[0],fmi[1],fmi[2]);

      //}

    }



  // variable force, wrap with clear/add



  } else {



    modify->clearstep_compute();



    if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar);

    else if (xstyle == ATOM)

      input->variable->compute_atom(xvar,igroup,&sforce[0][0],3,0);

    if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar);

    else if (ystyle == ATOM)

      input->variable->compute_atom(yvar,igroup,&sforce[0][1],3,0);

    if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar);

    else if (zstyle == ATOM)

      input->variable->compute_atom(zvar,igroup,&sforce[0][2],3,0);



    modify->addstep_compute(update->ntimestep + 1);



    if (mask[i] & groupbit) {

      //if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;

      if (region && region->match(x[i][0],x[i][1],x[i][2])) {

        foriginal[0] += fmi[0];

        foriginal[1] += fmi[1];

        foriginal[2] += fmi[2];

        if (xstyle == ATOM) fmi[0] = sforce[i][0];

        else if (xstyle) fmi[0] = xvalue;

        if (ystyle == ATOM) fmi[1] = sforce[i][1];

        else if (ystyle) fmi[1] = yvalue;

        if (zstyle == ATOM) fmi[2] = sforce[i][2];

        else if (zstyle) fmi[2] = zvalue;

      }

    }

  }

  //printf("test inside after setforce: %g %g %g \n",fmi[0],fmi[1],fmi[2]);

}



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



void FixSetForceSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)

{

  // set force to desired value on requested level, 0.0 on other levels
@@ -309,7 +214,6 @@ void FixSetForceSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/) 
    }



    double **x = atom->x;

    //double **f = atom->f;

    double **fm = atom->fm;

    int *mask = atom->mask;

    int nlocal = atom->nlocal;
@@ -323,36 +227,3 @@ void FixSetForceSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/) 
      }

  }

}
 


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



void FixSetForceSpin::min_post_force(int vflag)

{

  post_force(vflag);

}



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

   return components of total force on fix group before force was changed

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



double FixSetForceSpin::compute_vector(int n)

{

  // only sum across procs one time



  if (force_flag == 0) {

    MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world);

    force_flag = 1;

  }

  return foriginal_all[n];

}



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

   memory usage of local atom-based array

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



double FixSetForceSpin::memory_usage()

{

  double bytes = 0.0;

  if (varflag == ATOM) bytes = maxatom*3 * sizeof(double);

  return bytes;

}

src/SPIN/fix_setforce_spin.h

+3 −28
Original line number Diff line number Diff line
@@ -20,41 +20,16 @@ FixStyle(setforce/spin,FixSetForceSpin) 
#ifndef LMP_FIX_SET_FORCE_SPIN_H

#define LMP_FIX_SET_FORCE_SPIN_H



#include "fix.h"

//#include "fix_setforce.h"



// only post_force and post_force_respa != from FixSetForce

#include "fix_setforce.h"



namespace LAMMPS_NS {



//class FixSetForceSpin : public FixSetForce {

class FixSetForceSpin : public Fix {

class FixSetForceSpin : public FixSetForce {

 public:

  FixSetForceSpin(class LAMMPS *, int, char **);

  virtual ~FixSetForceSpin();

  int setmask();

  virtual void init();

  void setup(int);

  void min_setup(int);

  virtual void post_force(int);

  void post_force_respa(int, int, int);

  void min_post_force(int);

  double compute_vector(int);



  double memory_usage();



 protected:

  double xvalue,yvalue,zvalue;

  int varflag,iregion;

  char *xstr,*ystr,*zstr;

  char *idregion;

  int xvar,yvar,zvar,xstyle,ystyle,zstyle;

  double foriginal[3],foriginal_all[3];

  int force_flag;

  int nlevels_respa,ilevel_respa;



  int maxatom;

  double **sforce;

  void single_setforce_spin(int, double *); 

};



}

CRA Git | Maintained and supported by SUSTech CRA and CCSE