Commit 8ec231a8 authored by julient31's avatar julient31
Browse files

Commit JT 040819 

- corrected bug in fix nve/spin (init tables with grow)
- start work for fix setforce/spin
parent 89886032

src/SPIN/fix_nve_spin.cpp

+6 −10
Original line number Diff line number Diff line
@@ -88,7 +88,6 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : 
  npairs = 0;

  npairspin = 0;





  // checking if map array or hash is defined



  if (atom->map_style == 0)
@@ -250,13 +249,10 @@ void FixNVESpin::init() 
  // init. size of stacking lists (sectoring)



  nlocal_max = atom->nlocal;

  stack_head = memory->grow(stack_head,nsectors,"NVE/spin:stack_head");

  stack_foot = memory->grow(stack_foot,nsectors,"NVE/spin:stack_foot");

  backward_stacks = memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");

  forward_stacks = memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");

  if (nlocal_max == 0)

    error->all(FLERR,"Incorrect value of nlocal_max");



  memory->grow(stack_head,nsectors,"NVE/spin:stack_head");

  memory->grow(stack_foot,nsectors,"NVE/spin:stack_foot");

  memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");

  memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");

}



/* ---------------------------------------------------------------------- */
@@ -389,8 +385,8 @@ void FixNVESpin::pre_neighbor() 


  if (nlocal_max < nlocal) {                    // grow linked lists if necessary

    nlocal_max = nlocal;

    backward_stacks = memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");

    forward_stacks = memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");

    memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");

    memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");

  }



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

src/SPIN/fix_setforce_spin.cpp

0 → 100644
+358 −0
Original line number Diff line number Diff line 
/* ----------------------------------------------------------------------

   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.

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



#include <cstring>

#include <cstdlib>

//#include "fix_setforce.h"

#include "fix_setforce_spin.h"

#include "atom.h"

#include "update.h"

#include "modify.h"

#include "domain.h"

#include "region.h"

#include "respa.h"

#include "input.h"

#include "variable.h"

#include "memory.h"

#include "error.h"

#include "force.h"



using namespace LAMMPS_NS;

using namespace FixConst;



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

}



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



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;



  // 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) {

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

      if (mask[i] & groupbit) {

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

        foriginal[0] += fm[i][0];

        foriginal[1] += fm[i][1];

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

        if (xstyle) fm[i][0] = xvalue;

        if (ystyle) fm[i][1] = yvalue;

        if (zstyle) fm[i][2] = zvalue;

      }



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



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

      if (mask[i] & groupbit) {

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

        foriginal[0] += fm[i][0];

        foriginal[1] += fm[i][1];

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

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

        else if (xstyle) fm[i][0] = xvalue;

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

        else if (ystyle) fm[i][1] = yvalue;

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

        else if (zstyle) fm[i][2] = zvalue;

      }

  }

}



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



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

{

  // set force to desired value on requested level, 0.0 on other levels



  if (ilevel == ilevel_respa) post_force(vflag);

  else {

    Region *region = NULL;

    if (iregion >= 0) {

      region = domain->regions[iregion];

      region->prematch();

    }



    double **x = atom->x;

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

    double **fm = atom->fm;

    int *mask = atom->mask;

    int nlocal = atom->nlocal;



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

      if (mask[i] & groupbit) {

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

        if (xstyle) fm[i][0] = 0.0;

        if (ystyle) fm[i][1] = 0.0;

        if (zstyle) fm[i][2] = 0.0;

      }

  }

}



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



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

0 → 100644
+90 −0
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 FIX_CLASS



FixStyle(setforce/spin,FixSetForceSpin)



#else



#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



namespace LAMMPS_NS {



//class FixSetForceSpin : public FixSetForce {

class FixSetForceSpin : public Fix {

 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;

};



}



#endif

#endif



/* ERROR/WARNING messages:



E: Illegal ... command



Self-explanatory.  Check the input script syntax and compare to the

documentation for the command.  You can use -echo screen as a

command-line option when running LAMMPS to see the offending line.



E: Region ID for fix setforce does not exist



Self-explanatory.



E: Variable name for fix setforce does not exist



Self-explanatory.



E: Variable for fix setforce is invalid style



Only equal-style variables can be used.



E: Cannot use non-zero forces in an energy minimization



Fix setforce cannot be used in this manner.  Use fix addforce

instead.



*/

CRA Git | Maintained and supported by SUSTech CRA and CCSE