Commit ebe2ee9b authored by sjplimp's avatar sjplimp
Browse files 

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

src/ASPHERE/compute_erotate_asphere.cpp

+5 −5
Original line number Diff line number Diff line
@@ -25,7 +25,7 @@ using namespace LAMMPS_NS; 


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



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

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

  Compute(lmp, narg, arg)

{

  if (narg != 3) error->all("Illegal compute erotate/asphere command");
@@ -42,14 +42,14 @@ ComputeERotateASphere::ComputeERotateASphere(LAMMPS *lmp, int narg, char **arg) 


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



ComputeERotateASphere::~ComputeERotateASphere()

ComputeERotateAsphere::~ComputeERotateAsphere()

{

  memory->destroy_2d_double_array(inertia);

}



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



void ComputeERotateASphere::init()

void ComputeERotateAsphere::init()

{

  pfactor = 0.5 * force->mvv2e;
@@ -61,7 +61,7 @@ void ComputeERotateASphere::init() 


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



double ComputeERotateASphere::compute_scalar()

double ComputeERotateAsphere::compute_scalar()

{

  invoked |= INVOKED_SCALAR;
@@ -102,7 +102,7 @@ double ComputeERotateASphere::compute_scalar() 
   principal moments of inertia for ellipsoids

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



void ComputeERotateASphere::calculate_inertia()

void ComputeERotateAsphere::calculate_inertia()

{

  double *mass = atom->mass;

  double **shape = atom->shape;

src/ASPHERE/compute_erotate_asphere.h

+3 −3
Original line number Diff line number Diff line
@@ -18,10 +18,10 @@ 


namespace LAMMPS_NS {



class ComputeERotateASphere : public Compute {

class ComputeERotateAsphere : public Compute {

 public:

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

  ~ComputeERotateASphere();

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

  ~ComputeERotateAsphere();

  void init();

  double compute_scalar();

src/ASPHERE/compute_temp_asphere.cpp

+61 −0
Original line number Diff line number Diff line
@@ -70,6 +70,13 @@ void ComputeTempAsphere::init() 
    fix_dof += modify->fix[i]->dof(igroup);

  recount();



  if (id_bias) {

    tempbias = 1;

    int i = modify->find_compute(id_bias);

    if (i < 0) error->all("Could not find compute ID for temperature bias");

    tbias = modify->compute[i];

  } else tempbias = 0;



  calculate_inertia();

}
@@ -122,6 +129,12 @@ double ComputeTempAsphere::compute_scalar() 
{

  invoked |= INVOKED_SCALAR;



  if (tbias) {

    if (!(tbias->invoked & INVOKED_SCALAR))

      double tmp = tbias->compute_scalar();

    tbias->remove_bias_all();

  }



  double **v = atom->v;

  double **quat = atom->quat;

  double **angmom = atom->angmom;
@@ -163,6 +176,8 @@ double ComputeTempAsphere::compute_scalar() 
      }

    }



  if (tbias) tbias->restore_bias_all();



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

  if (dynamic) recount();

  scalar *= tfactor;
@@ -177,6 +192,11 @@ void ComputeTempAsphere::compute_vector() 


  invoked |= INVOKED_VECTOR;



  if (tbias) {

    if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector();

    tbias->remove_bias_all();

  }



  double **v = atom->v;

  double **quat = atom->quat;

  double **angmom = atom->angmom;
@@ -223,6 +243,8 @@ void ComputeTempAsphere::compute_vector() 
      t[5] += inertia[itype][2]*wbody[1]*wbody[2];

    }



  if (tbias) tbias->restore_bias_all();



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

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

}
@@ -245,3 +267,42 @@ void ComputeTempAsphere::calculate_inertia() 
      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;

  }

}



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

   remove velocity bias from atom I to leave thermal velocity

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



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

{

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

}



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

   remove velocity bias from all atoms to leave thermal velocity

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



void ComputeTempAsphere::remove_bias_all()

{

  if (tbias) tbias->remove_bias_all();

}



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

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

   assume remove_bias() was previously called

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



void ComputeTempAsphere::restore_bias(double *v)

{

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

}



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

   add back in velocity bias to all atoms removed by remove_bias_all()

   assume remove_bias_all() was previously called

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



void ComputeTempAsphere::restore_bias_all()

{

  if (tbias) tbias->restore_bias_all();

}

src/ASPHERE/compute_temp_asphere.h

+7 −0
Original line number Diff line number Diff line
@@ -26,11 +26,18 @@ class ComputeTempAsphere : public Compute { 
  double compute_scalar();

  void compute_vector();



  void remove_bias(int, double *);

  void remove_bias_all();

  void restore_bias(double *);

  void restore_bias_all();



 private:

  int fix_dof;

  double tfactor;

  double **inertia;



  Compute *tbias;     // ptr to additional bias compute



  void recount();

  void calculate_inertia();

};

src/ASPHERE/fix_npt_asphere.cpp

+8 −8
Original line number Diff line number Diff line
@@ -35,7 +35,7 @@ enum{NOBIAS,BIAS}; 


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



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

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

  FixNPT(lmp, narg, arg)

{

  if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag ||
@@ -46,7 +46,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : 


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



void FixNPTASphere::init()

void FixNPTAsphere::init()

{

  FixNPT::init();

  dtq = 0.5 * update->dt;
@@ -56,7 +56,7 @@ void FixNPTASphere::init() 
   1st half of Verlet update 

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



void FixNPTASphere::initial_integrate(int vflag)

void FixNPTAsphere::initial_integrate(int vflag)

{

  int i;

  double dtfm;
@@ -168,7 +168,7 @@ void FixNPTASphere::initial_integrate(int vflag) 
   2nd half of Verlet update 

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



void FixNPTASphere::final_integrate()

void FixNPTAsphere::final_integrate()

{

  int i;

  double dtfm;
@@ -251,7 +251,7 @@ void FixNPTASphere::final_integrate() 


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



void FixNPTASphere::reset_dt()

void FixNPTAsphere::reset_dt()

{

  FixNPT::reset_dt();

  dtq = 0.5 * update->dt;
@@ -261,7 +261,7 @@ void FixNPTASphere::reset_dt() 
   Richardson iteration to update quaternion accurately

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



void FixNPTASphere::richardson(double *q, double *m, double *moments)

void FixNPTAsphere::richardson(double *q, double *m, double *moments)

{

  // compute omega at 1/2 step from m at 1/2 step and q at 0
@@ -320,7 +320,7 @@ void FixNPTASphere::richardson(double *q, double *m, double *moments) 
     and divide by principal moments

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



void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia,

void FixNPTAsphere::omega_from_mq(double *q, double *m, double *inertia,

				  double *w)

{

  double rot[3][3];
@@ -339,7 +339,7 @@ void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia, 
   shape = x,y,z radii in body frame

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



void FixNPTASphere::calculate_inertia(double mass, double *shape,

void FixNPTAsphere::calculate_inertia(double mass, double *shape,

				      double *inertia)

{

  inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0;

CRA Git | Maintained and supported by SUSTech CRA and CCSE