Commit 6c6d45d7 authored by sjplimp's avatar sjplimp
Browse files 

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

src/compute_temp_com.cpp

0 → 100644
+159 −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 "stdlib.h"

#include "string.h"

#include "compute_temp_com.h"

#include "atom.h"

#include "force.h"

#include "group.h"

#include "modify.h"

#include "fix.h"

#include "domain.h"

#include "lattice.h"

#include "error.h"



using namespace LAMMPS_NS;



#define MIN(A,B) ((A) < (B)) ? (A) : (B)

#define MAX(A,B) ((A) > (B)) ? (A) : (B)



#define INVOKED_SCALAR 1

#define INVOKED_VECTOR 2



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



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

  Compute(lmp, narg, arg)

{

  if (narg != 3) error->all("Illegal compute temp command");



  scalar_flag = vector_flag = 1;

  size_vector = 6;

  extscalar = 0;

  extvector = 1;

  tempflag = 1;



  vector = new double[6];

}



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



ComputeTempCOM::~ComputeTempCOM()

{

  delete [] vector;

}



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



void ComputeTempCOM::init()

{

  fix_dof = 0;

  for (int i = 0; i < modify->nfix; i++)

    fix_dof += modify->fix[i]->dof(igroup);

  recount();

  masstotal = group->mass(igroup);

}



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



void ComputeTempCOM::recount()

{

  double natoms = group->count(igroup);

  dof = domain->dimension * natoms;

  dof -= extra_dof + fix_dof;

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

  else tfactor = 0.0;

}



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



double ComputeTempCOM::compute_scalar()

{

  double vcm[3],vthermal[3];



  invoked |= INVOKED_SCALAR;



  if (dynamic) masstotal = group->mass(igroup);

  group->vcm(igroup,masstotal,vcm);



  double **v = atom->v;

  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



  double t = 0.0;

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

    if (mask[i] & groupbit) {

      vthermal[0] = v[i][0] - vcm[0];

      vthermal[1] = v[i][1] - vcm[1];

      vthermal[2] = v[i][2] - vcm[2];

      if (mass)

	t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + 

	      vthermal[2]*vthermal[2]) * mass[type[i]];

      else

	t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + 

	      vthermal[2]*vthermal[2]) * rmass[i];

    }



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

  if (dynamic) recount();

  scalar *= tfactor;

  return scalar;

}



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



void ComputeTempCOM::compute_vector()

{

  int i;

  double vcm[3],vthermal[3];



  invoked |= INVOKED_VECTOR;



  if (dynamic) masstotal = group->mass(igroup);

  group->vcm(igroup,masstotal,vcm);



  double **x = atom->x;

  double **v = atom->v;

  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;



  double massone,t[6];

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



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

    if (mask[i] & groupbit) {

      vthermal[0] = v[i][0] - vcm[0];

      vthermal[1] = v[i][1] - vcm[1];

      vthermal[2] = v[i][2] - vcm[2];



      if (mass) massone = mass[type[i]];

      else massone = rmass[i];

      t[0] += massone * vthermal[0]*vthermal[0];

      t[1] += massone * vthermal[1]*vthermal[1];

      t[2] += massone * vthermal[2]*vthermal[2];

      t[3] += massone * vthermal[0]*vthermal[1];

      t[4] += massone * vthermal[0]*vthermal[2];

      t[5] += massone * vthermal[1]*vthermal[2];

    }



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

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

}

src/compute_temp_com.h

0 → 100644
+38 −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_TEMP_COM_H

#define COMPUTE_TEMP_COM_H



#include "compute.h"



namespace LAMMPS_NS {



class ComputeTempCOM : public Compute {

 public:

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

  ~ComputeTempCOM();

  void init();

  double compute_scalar();

  void compute_vector();



 private:

  int fix_dof;

  double tfactor,masstotal;



  void recount();

};



}



#endif

src/style.h

+2 −0
Original line number Diff line number Diff line
@@ -89,6 +89,7 @@ CommandStyle(write_restart,WriteRestart) 
#include "compute_rotate_gran.h"

#include "compute_stress_atom.h"

#include "compute_temp.h"

#include "compute_temp_com.h"

#include "compute_temp_deform.h"

#include "compute_temp_partial.h"

#include "compute_temp_ramp.h"
@@ -109,6 +110,7 @@ ComputeStyle(rotate/dipole,ComputeRotateDipole) 
ComputeStyle(rotate/gran,ComputeRotateGran)

ComputeStyle(stress/atom,ComputeStressAtom)

ComputeStyle(temp,ComputeTemp)

ComputeStyle(temp/com,ComputeTempCOM)

ComputeStyle(temp/deform,ComputeTempDeform)

ComputeStyle(temp/partial,ComputeTempPartial)

ComputeStyle(temp/ramp,ComputeTempRamp)

CRA Git | Maintained and supported by SUSTech CRA and CCSE