Commit a62fdded authored by sjplimp's avatar sjplimp
Browse files

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1643 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent d96de318
Loading
Loading
Loading
Loading
+4 −0
Original line number Diff line number Diff line
@@ -5,6 +5,7 @@ if ($1 == 1) then
  cp style_asphere.h ..

  cp atom_vec_ellipsoid.cpp ..
  cp compute_erotate_asphere.cpp ..
  cp compute_temp_asphere.cpp ..
  cp fix_npt_asphere.cpp ..
  cp fix_nve_asphere.cpp ..
@@ -13,6 +14,7 @@ if ($1 == 1) then
  cp pair_resquared.cpp ..

  cp atom_vec_ellipsoid.h ..
  cp compute_erotate_asphere.h ..
  cp compute_temp_asphere.h ..
  cp fix_npt_asphere.h ..
  cp fix_nve_asphere.h ..
@@ -26,6 +28,7 @@ else if ($1 == 0) then
  touch ../style_asphere.h

  rm ../atom_vec_ellipsoid.cpp
  rm ../compute_erotate_asphere.cpp
  rm ../compute_temp_asphere.cpp
  rm ../fix_npt_asphere.cpp
  rm ../fix_nve_asphere.cpp
@@ -34,6 +37,7 @@ else if ($1 == 0) then
  rm ../pair_resquared.cpp

  rm ../atom_vec_ellipsoid.h
  rm ../compute_erotate_asphere.h
  rm ../compute_temp_asphere.h
  rm ../fix_npt_asphere.h
  rm ../fix_nve_asphere.h
+118 −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 "mpi.h"
#include "compute_erotate_asphere.h"
#include "math_extra.h"
#include "atom.h"
#include "force.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define INVOKED_SCALAR 1

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

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

  if (!atom->angmom_flag || !atom->quat_flag)
    error->all("Compute erotate/asphere requires atom attributes angmom, quat");

  scalar_flag = 1;
  extscalar = 1;

  inertia = 
    memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia");
}

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

ComputeERotateASphere::~ComputeERotateASphere()
{
  memory->destroy_2d_double_array(inertia);
}

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

void ComputeERotateASphere::init()
{
  pfactor = 0.5 * force->mvv2e;

  if (!atom->shape)
    error->all("Compute erotate/asphere requires atom attribute shape");

  calculate_inertia();
}

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

double ComputeERotateASphere::compute_scalar()
{
  invoked |= INVOKED_SCALAR;

  double **quat = atom->quat;
  double **angmom = atom->angmom;
  int *mask = atom->mask;
  int *type = atom->type;
  int nlocal = atom->nlocal;

  int itype;
  double wbody[3];
  double rot[3][3];
  double erotate = 0.0;

  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      itype = type[i];

      // wbody = angular velocity in body frame

      MathExtra::quat_to_mat(quat[i],rot);
      MathExtra::transpose_times_column3(rot,angmom[i],wbody);
      wbody[0] /= inertia[itype][0];
      wbody[1] /= inertia[itype][1];
      wbody[2] /= inertia[itype][2];
      
      erotate += inertia[itype][0]*wbody[0]*wbody[0]+
	inertia[itype][1]*wbody[1]*wbody[1]+
	inertia[itype][2]*wbody[2]*wbody[2];
    }

  MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
  scalar *= pfactor;
  return scalar;
}

/* ----------------------------------------------------------------------
   principal moments of inertia for ellipsoids
------------------------------------------------------------------------- */

void ComputeERotateASphere::calculate_inertia()
{
  double *mass = atom->mass;
  double **shape = atom->shape;

  for (int i = 1; i <= atom->ntypes; i++) {
    inertia[i][0] = mass[i] * 
      (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]) / 5.0;
    inertia[i][1] = mass[i] * 
      (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]) / 5.0;
    inertia[i][2] = mass[i] * 
      (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0;
  }
}
+37 −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.
------------------------------------------------------------------------- */

#ifndef COMPUTE_EROTATE_ASPHERE_H
#define COMPUTE_EROTATE_ASPHERE_H

#include "compute.h"

namespace LAMMPS_NS {

class ComputeERotateASphere : public Compute {
 public:
  ComputeERotateASphere(class LAMMPS *, int, char **);
  ~ComputeERotateASphere();
  void init();
  double compute_scalar();

 private:
  double pfactor;
  double **inertia;

  void calculate_inertia();
};

}

#endif
+2 −0
Original line number Diff line number Diff line
@@ -20,10 +20,12 @@ AtomStyle(ellipsoid,AtomVecEllipsoid)
# endif

#ifdef ComputeInclude
#include "compute_erotate_asphere.h"
#include "compute_temp_asphere.h"
#endif

#ifdef ComputeClass
ComputeStyle(erotate/asphere,ComputeERotateASphere)
ComputeStyle(temp/asphere,ComputeTempAsphere)
#endif