Commit 3db8850f authored by sjplimp's avatar sjplimp
Browse files 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14598 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent e98647ef

src/BODY/compute_temp_body.cpp

0 → 100755
+394 −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.

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



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

   Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)

   based on ComputeTempAsphere

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



#include <mpi.h>

#include <string.h>

#include "compute_temp_body.h"

#include "math_extra.h"

#include "atom.h"

#include "atom_vec_body.h"

#include "update.h"

#include "force.h"

#include "domain.h"

#include "modify.h"

#include "group.h"

#include "memory.h"

#include "error.h"



using namespace LAMMPS_NS;



enum{ROTATE,ALL};



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



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

  Compute(lmp, narg, arg)

{

  if (narg < 3) error->all(FLERR,"Illegal compute temp/body command");



  scalar_flag = vector_flag = 1;

  size_vector = 6;

  extscalar = 0;

  extvector = 1;

  tempflag = 1;



  tempbias = 0;

  id_bias = NULL;

  mode = ALL;



  int iarg = 3;

  while (iarg < narg) {

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

      if (iarg+2 > narg)

        error->all(FLERR,"Illegal compute temp/body command");

      tempbias = 1;

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

      id_bias = new char[n];

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

      iarg += 2;

    } else if (strcmp(arg[iarg],"dof") == 0) {

      if (iarg+2 > narg)

        error->all(FLERR,"Illegal compute temp/body command");

      if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;

      else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;

      else error->all(FLERR,"Illegal compute temp/body command");

      iarg += 2;

    } else error->all(FLERR,"Illegal compute temp/body command");

  }



  vector = new double[6];



}



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



ComputeTempBody::~ComputeTempBody()

{

  delete [] id_bias;

  delete [] vector;

}



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



void ComputeTempBody::init()

{

  // error check



  avec = (AtomVecBody *) atom->style_match("body");

  if (!avec)

    error->all(FLERR,"Compute temp/body requires atom style body");



  // check that all particles are finite-size, no point particles allowed



  int *body = atom->body;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



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

    if (mask[i] & groupbit)

      if (body[i] < 0)

        error->one(FLERR,"Compute temp/body requires bodies");



  if (tempbias) {

    int i = modify->find_compute(id_bias);

    if (i < 0)

      error->all(FLERR,"Could not find compute ID for temperature bias");

    tbias = modify->compute[i];

    if (tbias->tempflag == 0)

      error->all(FLERR,"Bias compute does not calculate temperature");

    if (tbias->tempbias == 0)

      error->all(FLERR,"Bias compute does not calculate a velocity bias");

    if (tbias->igroup != igroup)

      error->all(FLERR,"Bias compute group does not match compute group");

    if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2;

    else tempbias = 1;



    // init and setup bias compute because

    // this compute's setup()->dof_compute() may be called first



    tbias->init();

    tbias->setup();

  }

}



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



void ComputeTempBody::setup()

{

  dynamic = 0;

  if (dynamic_user || group->dynamic[igroup]) dynamic = 1;

  dof_compute();

}



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



void ComputeTempBody::dof_compute()

{

  adjust_dof_fix();



  // 6 dof for 3d, 3 dof for 2d

  // which dof are included also depends on mode

  // assume full rotation of extended particles

  // user should correct this via compute_modify if needed



  natoms_temp = group->count(igroup);

  int nper;

  if (domain->dimension == 3) {

    if (mode == ALL) nper = 6;

    else nper = 3;

  } else {

    if (mode == ALL) nper = 3;

    else nper = 1;

  }

  dof = nper*natoms_temp;



  // additional adjustments to dof



  if (tempbias == 1) {

    if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms_temp;



  } else if (tempbias == 2) {

    int *mask = atom->mask;

    int nlocal = atom->nlocal;



    tbias->dof_remove_pre();



    int count = 0;

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

      if (mask[i] & groupbit)

        if (tbias->dof_remove(i)) count++;

    int count_all;

    MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world);

    dof -= nper*count_all;

  }



  dof -= extra_dof + fix_dof;

  if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);

  else tfactor = 0.0;

}



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



double ComputeTempBody::compute_scalar()

{

  invoked_scalar = update->ntimestep;



  if (tempbias) {

    if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();

    tbias->remove_bias_all();

  }



  AtomVecBody::Bonus *bonus = avec->bonus;

  double **v = atom->v;

  double **angmom = atom->angmom;

  double *rmass = atom->rmass;

  int *body = atom->body;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



  double *inertia,*quat;

  double wbody[3];

  double rot[3][3];



  // sum translational and rotational energy for each particle



  double t = 0.0;



  if (mode == ALL) {

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

      if (mask[i] & groupbit) {

        t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];



        // principal moments of inertia



        inertia = bonus[body[i]].inertia;

        quat = bonus[body[i]].quat;



        // wbody = angular velocity in body frame



        MathExtra::quat_to_mat(quat,rot);

        MathExtra::transpose_matvec(rot,angmom[i],wbody);

        if (inertia[0] == 0.0) wbody[0] = 0.0;

        else wbody[0] /= inertia[0];

        if (inertia[1] == 0.0) wbody[1] = 0.0;

        else wbody[1] /= inertia[1];

        if (inertia[2] == 0.0) wbody[2] = 0.0;

        else wbody[2] /= inertia[2];



        t += inertia[0]*wbody[0]*wbody[0] +

          inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];

      }



  } else {

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

      if (mask[i] & groupbit) {



        // principal moments of inertia



        inertia = bonus[body[i]].inertia;

        quat = bonus[body[i]].quat;



        // wbody = angular velocity in body frame



        MathExtra::quat_to_mat(quat,rot);

        MathExtra::transpose_matvec(rot,angmom[i],wbody);

        if (inertia[0] == 0.0) wbody[0] = 0.0;

        else wbody[0] /= inertia[0];

        if (inertia[1] == 0.0) wbody[1] = 0.0;

        else wbody[1] /= inertia[1];

        if (inertia[2] == 0.0) wbody[2] = 0.0;

        else wbody[2] /= inertia[2];



        t += inertia[0]*wbody[0]*wbody[0] +

          inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];

      }

  }



  if (tempbias) tbias->restore_bias_all();



  MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);

  if (dynamic || tempbias == 2) dof_compute();

  if (dof < 0.0 && natoms_temp > 0.0)

    error->all(FLERR,"Temperature compute degrees of freedom < 0");

  scalar *= tfactor;

  return scalar;

}



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



void ComputeTempBody::compute_vector()

{

  int i;



  invoked_vector = update->ntimestep;



  if (tempbias) {

    if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector();

    tbias->remove_bias_all();

  }



  AtomVecBody::Bonus *bonus = avec->bonus;

  double **v = atom->v;

  double **angmom = atom->angmom;

  double *rmass = atom->rmass;

  int *body = atom->body;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



  double *inertia,*quat;

  double wbody[3],t[6];

  double rot[3][3];

  double massone;



  // sum translational and rotational energy for each particle



  for (i = 0; i < 6; i++) t[i] = 0.0;



  if (mode == ALL) {

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

      if (mask[i] & groupbit) {

        massone = rmass[i];

        t[0] += massone * v[i][0]*v[i][0];

        t[1] += massone * v[i][1]*v[i][1];

        t[2] += massone * v[i][2]*v[i][2];

        t[3] += massone * v[i][0]*v[i][1];

        t[4] += massone * v[i][0]*v[i][2];

        t[5] += massone * v[i][1]*v[i][2];



        // principal moments of inertia



        inertia = bonus[body[i]].inertia;

        quat = bonus[body[i]].quat;



        // wbody = angular velocity in body frame



        MathExtra::quat_to_mat(quat,rot);

        MathExtra::transpose_matvec(rot,angmom[i],wbody);

        if (inertia[0] == 0.0) wbody[0] = 0.0;

        else wbody[0] /= inertia[0];

        if (inertia[1] == 0.0) wbody[1] = 0.0;

        else wbody[1] /= inertia[1];

        if (inertia[2] == 0.0) wbody[2] = 0.0;

        else wbody[2] /= inertia[2];



        // rotational kinetic energy



        t[0] += inertia[0]*wbody[0]*wbody[0];

        t[1] += inertia[1]*wbody[1]*wbody[1];

        t[2] += inertia[2]*wbody[2]*wbody[2];

        t[3] += inertia[0]*wbody[0]*wbody[1];

        t[4] += inertia[1]*wbody[0]*wbody[2];

        t[5] += inertia[2]*wbody[1]*wbody[2];

      }



  } else {

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

      if (mask[i] & groupbit) {



        // principal moments of inertia



        inertia = bonus[body[i]].inertia;

        quat = bonus[body[i]].quat;

        massone = rmass[i];



        // wbody = angular velocity in body frame



        MathExtra::quat_to_mat(quat,rot);

        MathExtra::transpose_matvec(rot,angmom[i],wbody);

        if (inertia[0] == 0.0) wbody[0] = 0.0;

        else wbody[0] /= inertia[0];

        if (inertia[1] == 0.0) wbody[1] = 0.0;

        else wbody[1] /= inertia[1];

        if (inertia[2] == 0.0) wbody[2] = 0.0;

        else wbody[2] /= inertia[2];



        // rotational kinetic energy



        t[0] += inertia[0]*wbody[0]*wbody[0];

        t[1] += inertia[1]*wbody[1]*wbody[1];

        t[2] += inertia[2]*wbody[2]*wbody[2];

        t[3] += inertia[0]*wbody[0]*wbody[1];

        t[4] += inertia[1]*wbody[0]*wbody[2];

        t[5] += inertia[2]*wbody[1]*wbody[2];

      }

  }



  if (tempbias) tbias->restore_bias_all();



  MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);

  for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;

}



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

   remove velocity bias from atom I to leave thermal velocity

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



void ComputeTempBody::remove_bias(int i, double *v)

{

  if (tbias) tbias->remove_bias(i,v);

}



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

   add back in velocity bias to atom I removed by remove_bias()

   assume remove_bias() was previously called

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



void ComputeTempBody::restore_bias(int i, double *v)

{

  if (tbias) tbias->restore_bias(i,v);

}

src/BODY/compute_temp_body.h

0 → 100755
+96 −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 COMPUTE_CLASS



ComputeStyle(temp/body,ComputeTempBody)



#else



#ifndef LMP_COMPUTE_TEMP_BODY_H

#define LMP_COMPUTE_TEMP_BODY_H



#include "compute.h"



namespace LAMMPS_NS {



class ComputeTempBody : public Compute {

 public:

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

  ~ComputeTempBody();

  void init();

  void setup();

  double compute_scalar();

  void compute_vector();



  void remove_bias(int, double *);

  void restore_bias(int, double *);



 private:

  int mode;

  double tfactor;

  char *id_bias;

  class Compute *tbias;              // ptr to additional bias compute

  class AtomVecBody *avec;



  void dof_compute();

};



}



#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: Compute temp/body requires atom style body



Self-explanatory.



E: Compute temp/body requires bodies



This compute can only be applied to body particles.



E: Could not find compute ID for temperature bias



Self-explanatory.



E: Bias compute does not calculate temperature



The specified compute must compute temperature.



E: Bias compute does not calculate a velocity bias



The specified compute must compute a bias for temperature.



E: Bias compute group does not match compute group



The specified compute must operate on the same group as the parent

compute.



E: Temperature compute degrees of freedom < 0



This should not happen if you are calculating the temperature

on a valid set of atoms.



U: Compute temp/asphere requires atom style ellipsoid



Self-explanatory.



*/

src/BODY/fix_nh_body.cpp

0 → 100644
+154 −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.

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



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

   Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)

   based on FixNHAsphere

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



#include <string.h>

#include <stdlib.h>

#include <math.h>

#include "math_extra.h"

#include "fix_nh_body.h"

#include "atom.h"

#include "atom_vec_body.h"

#include "group.h"

#include "memory.h"

#include "error.h"



using namespace LAMMPS_NS;

using namespace FixConst;



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



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

  FixNH(lmp, narg, arg)

{

}



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



void FixNHBody::init()

{

  avec = (AtomVecBody *) atom->style_match("body");

  if (!avec)

    error->all(FLERR,

               "Compute nvt/nph/npt body requires atom style body");



  // check that all particles are finite-size

  // no point particles allowed, spherical is OK



  int *body = atom->body;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



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

    if (mask[i] & groupbit)

      if (body[i] < 0)

        error->one(FLERR,"Fix nvt/nph/npt body requires bodies");



  FixNH::init();

}



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

   perform half-step update of angular momentum

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



void FixNHBody::nve_v()

{

  // standard nve_v velocity update



  FixNH::nve_v();



  double **angmom = atom->angmom;

  double **torque = atom->torque;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  if (igroup == atom->firstgroup) nlocal = atom->nfirst;



  // update angular momentum by 1/2 step for all particles



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

    if (mask[i] & groupbit) {

      angmom[i][0] += dtf*torque[i][0];

      angmom[i][1] += dtf*torque[i][1];

      angmom[i][2] += dtf*torque[i][2];

    }

  }

}



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

   perform full-step update of orientation

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



void FixNHBody::nve_x()

{

  double omega[3];

  double *quat,*inertia;



  // standard nve_x position update



  FixNH::nve_x();



  AtomVecBody::Bonus *bonus = avec->bonus;

  int *body = atom->body;

  double **angmom = atom->angmom;

  double *rmass = atom->rmass;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  if (igroup == atom->firstgroup) nlocal = atom->nfirst;



  // set timestep here since dt may have changed or come via rRESPA



  dtq = 0.5 * dtv;



  // update quaternion a full step via Richardson iteration



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

    if (mask[i] & groupbit) {



      // compute omega at 1/2 step from angmom at 1/2 step and current q

      // update quaternion a full step via Richardson iteration

      // returns new normalized quaternion



      inertia = bonus[body[i]].inertia;

      quat = bonus[body[i]].quat;

      MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);

      MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);

    }

}



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

   perform half-step temperature scaling of angular momentum

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



void FixNHBody::nh_v_temp()

{

  // standard nh_v_temp scaling



  FixNH::nh_v_temp();



  double **angmom = atom->angmom;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  if (igroup == atom->firstgroup) nlocal = atom->nfirst;



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

    if (mask[i] & groupbit) {

      angmom[i][0] *= factor_eta;

      angmom[i][1] *= factor_eta;

      angmom[i][2] *= factor_eta;

    }

  }

}

src/BODY/fix_nh_body.h

0 → 100644
+50 −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.

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



#ifndef LMP_FIX_NH_BODY_H

#define LMP_FIX_NH_BODY_H



#include "fix_nh.h"



namespace LAMMPS_NS {



class FixNHBody : public FixNH {

 public:

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

  virtual ~FixNHBody() {}

  void init();



 protected:

  double dtq;

  class AtomVecBody *avec;



  void nve_v();

  void nve_x();

  void nh_v_temp();

};



}



#endif



/* ERROR/WARNING messages:



E: Compute nvt/nph/npt body requires atom style body



Self-explanatory.



E: Fix nvt/nph/npt body requires bodies



Self-explanatory.



*/

src/BODY/fix_nph_body.cpp

0 → 100644
+72 −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.

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



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

   Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)

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



#include <string.h>

#include "fix_nph_body.h"

#include "modify.h"

#include "error.h"



using namespace LAMMPS_NS;

using namespace FixConst;



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



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

  FixNHBody(lmp, narg, arg)

{

  if (tstat_flag)

    error->all(FLERR,"Temperature control can not be used with fix nph/body");

  if (!pstat_flag)

    error->all(FLERR,"Pressure control must be used with fix nph/body");



  // create a new compute temp style

  // id = fix-ID + temp

  // compute group = all since pressure is always global (group all)

  // and thus its KE/temperature contribution should use group all



  int n = strlen(id) + 6;

  id_temp = new char[n];

  strcpy(id_temp,id);

  strcat(id_temp,"_temp");



  char **newarg = new char*[3];

  newarg[0] = id_temp;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "temp/body";



  modify->add_compute(3,newarg);

  delete [] newarg;

  tflag = 1;



  // create a new compute pressure style

  // id = fix-ID + press, compute group = all

  // pass id_temp as 4th arg to pressure constructor



  n = strlen(id) + 7;

  id_press = new char[n];

  strcpy(id_press,id);

  strcat(id_press,"_press");



  newarg = new char*[4];

  newarg[0] = id_press;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "pressure";

  newarg[3] = id_temp;

  modify->add_compute(4,newarg);

  delete [] newarg;

  pflag = 1;

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE