forked from lijiext/lammps
Added compute style temp/com
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1457 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
parent
2a2516693b
commit
eff5d1d72b
|
|
@ -466,7 +466,7 @@ optimized pair potentials for lj/cut, charmm/long, eam, morse : \
|
|||
James Fischer (High Performance Technologies), \
|
||||
David Richie and Vincent Natoli (Stone Ridge Technologies)
|
||||
fix wall/lj126 : Mark Stevens (Sandia)
|
||||
Stillinger-Weber and Tersoff potentials : Aidan Thompson (Sandia)
|
||||
Stillinger-Weber and Tersoff potentials : Aidan Thompson and Xiaowang Zhou (Sandia)
|
||||
region prism : Pieter in 't Veld (Sandia)
|
||||
LJ tail corrections for energy/pressure : Paul Crozier (Sandia)
|
||||
fix momentum and recenter : Naveen Michaud-Agrawal (Johns Hopkins U)
|
||||
|
|
|
|||
|
|
@ -119,6 +119,7 @@ available in LAMMPS:
|
|||
"sum"_compute_sum.html - sum per-atom quantities to a global value
|
||||
"temp"_compute_temp.html - temperature of group of atoms
|
||||
"temp/asphere"_compute_temp_asphere.html - temperature of aspherical particles
|
||||
"temp/com"_compute_temp_com.html - include center-of-mass momentum all directly-bonded neighbors
|
||||
"temp/deform"_compute_temp_deform.html - temperature excluding box deformation velocity
|
||||
"temp/dipole"_compute_temp_dipole.html - temperature of point dipolar particles
|
||||
"temp/partial"_compute_temp_partial.html - temperature excluding one or more dimensions of velocity
|
||||
|
|
|
|||
|
|
@ -0,0 +1,69 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
compute temp/com command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID temp/com :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command
|
||||
temp/com = style name of this compute command
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute newT oxygen temp/com :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
Define a compute to calculate the temperature of a group of atoms,
|
||||
where the momentum of each atom includes that of atoms to which
|
||||
it is directly bonded.
|
||||
A compute of this style can be used by any command that computes a temperature,
|
||||
e.g. "thermo_modify"_thermo_modify.html, "fix
|
||||
temp/rescale"_fix_temp_rescale.html, "fix npt"_fix_npt.html, etc.
|
||||
|
||||
The temperature is calculated by the formula KE = dim/2 N k T, where
|
||||
KE = total kinetic energy of the group of atoms
|
||||
The kinetic energy of each atom is based on the center of mass motion
|
||||
of itself and all the atoms directly bonded to it.
|
||||
N = number of atoms in the
|
||||
group, k = Boltzmann constant, and T = temperature.
|
||||
The dim parameter is adjusted to give the correct number of
|
||||
degrees of freedom.
|
||||
|
||||
The number of atoms contributing to the temperature is assumed to be
|
||||
constant for the duration of the run; use the {dynamic} option of the
|
||||
"compute_modify"_compute_modify.html command if this is not the case.
|
||||
|
||||
This command can also be used to compute a vector of two quantities. The
|
||||
first is the temperature described above; the second is the temperature
|
||||
based on the complementary part of the kinetic energy i.e. it uses the
|
||||
momentum of the atom and its directly-bonded neighbors, relative to their
|
||||
center of mass. It is scaled by the total number of atoms minus the total
|
||||
number of group atoms.
|
||||
|
||||
This compute should only be used in cases where none of the atoms in the group
|
||||
have any directly-bonded neighbors in common, including each other.
|
||||
If there are no directly-bonded
|
||||
neighbors, then this compute is equivalent to "compute temp"
|
||||
|
||||
[Output info:]
|
||||
|
||||
All values calculated by this compute are "intensive",
|
||||
meaning they are independent of the number of atoms in the simulation.
|
||||
|
||||
[Restrictions:] none
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"compute temp"_compute_temp.html, "compute
|
||||
temp/region"_compute_temp_region.html, "compute
|
||||
pressure"_compute_pressure.html
|
||||
|
||||
[Default:] none
|
||||
|
|
@ -5,11 +5,11 @@ SHELL = /bin/sh
|
|||
# System-specific settings
|
||||
|
||||
CC = c++
|
||||
CCFLAGS = -O -I../STUBS -I/Users/sjplimp/tools/fftw/include -DFFT_FFTW
|
||||
CCFLAGS = -O -I../STUBS -I/sw/include -DFFT_FFTW
|
||||
DEPFLAGS = -M
|
||||
LINK = c++
|
||||
LINKFLAGS = -O -L../STUBS -L/Users/sjplimp/tools/fftw/lib
|
||||
USRLIB = -lfftw -lmpi
|
||||
LINKFLAGS = -O -L/sw/lib
|
||||
USRLIB = -lfftw ../STUBS/mpi.o
|
||||
SYSLIB =
|
||||
SIZE = size
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,282 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "compute_temp_com.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "group.h"
|
||||
#include "error.h"
|
||||
#include "comm.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#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/com command");
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
size_vector = 2;
|
||||
extscalar = 0;
|
||||
extvector = 0;
|
||||
tempflag = 1;
|
||||
|
||||
vector = new double[2];
|
||||
|
||||
// set comm size needed by this Compute
|
||||
|
||||
comm_forward = 3;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
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();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeTempCoM::recount()
|
||||
{
|
||||
double natoms = group->count(igroup);
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *num_bond = atom->num_bond;
|
||||
|
||||
int atom1,m,i,nsum,nsumall;
|
||||
|
||||
nsum = 0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
nsum += num_bond[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&nsum,&nsumall,1,MPI_INT,MPI_SUM,world);
|
||||
|
||||
dof = 3 * natoms;
|
||||
dof -= extra_dof + fix_dof;
|
||||
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
|
||||
else tfactor = 0.0;
|
||||
|
||||
dof1 = 3 * nsumall;
|
||||
dof1 -= extra_dof + fix_dof;
|
||||
if (dof1 > 0) tfactor1 = force->mvv2e / (dof1 * force->boltz);
|
||||
else tfactor1 = 0.0;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double ComputeTempCoM::compute_scalar()
|
||||
{
|
||||
invoked |= INVOKED_SCALAR;
|
||||
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *num_bond = atom->num_bond;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int *tag = atom->tag;
|
||||
|
||||
double massone,t;
|
||||
double psum[3],msum;
|
||||
int atom1,m,i;
|
||||
|
||||
// communicate ghost particle velocities
|
||||
|
||||
comm->comm_compute(this);
|
||||
|
||||
t = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (mass) massone = mass[type[i]];
|
||||
else massone = rmass[i];
|
||||
psum[0] = v[i][0] * massone;
|
||||
psum[1] = v[i][1] * massone;
|
||||
psum[2] = v[i][2] * massone;
|
||||
msum = massone;
|
||||
for (m = 0; m < num_bond[i]; m++) {
|
||||
atom1 = atom->map(bond_atom[i][m]);
|
||||
if (atom1 == -1) {
|
||||
char str[128];
|
||||
sprintf(str,"Bond atoms %d %d missing",
|
||||
tag[i],bond_atom[i][m]);
|
||||
error->one(str);
|
||||
}
|
||||
if (mass) massone = mass[type[atom1]];
|
||||
else massone = rmass[atom1];
|
||||
psum[0] += v[atom1][0] * massone;
|
||||
psum[1] += v[atom1][1] * massone;
|
||||
psum[2] += v[atom1][2] * massone;
|
||||
msum += massone;
|
||||
}
|
||||
t += (psum[0]*psum[0] + psum[1]*psum[1] +
|
||||
psum[2]*psum[2]) / msum;
|
||||
|
||||
}
|
||||
|
||||
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
if (dynamic) recount();
|
||||
scalar *= tfactor;
|
||||
return scalar;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeTempCoM::compute_vector()
|
||||
{
|
||||
int i;
|
||||
|
||||
invoked |= INVOKED_VECTOR;
|
||||
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *num_bond = atom->num_bond;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int *tag = atom->tag;
|
||||
|
||||
double massone,t[2];
|
||||
double psum[3],msum,vav[3];
|
||||
int atom1,m;
|
||||
|
||||
// communicate ghost particle velocities
|
||||
|
||||
comm->comm_compute(this);
|
||||
|
||||
t[0] = 0.0;
|
||||
t[1] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (mass) massone = mass[type[i]];
|
||||
else massone = rmass[i];
|
||||
psum[0] = v[i][0] * massone;
|
||||
psum[1] = v[i][1] * massone;
|
||||
psum[2] = v[i][2] * massone;
|
||||
msum = massone;
|
||||
for (m = 0; m < num_bond[i]; m++) {
|
||||
atom1 = atom->map(bond_atom[i][m]);
|
||||
if (atom1 == -1) {
|
||||
char str[128];
|
||||
sprintf(str,"Bond atoms %d %d missing",
|
||||
tag[i],bond_atom[i][m]);
|
||||
error->one(str);
|
||||
}
|
||||
if (mass) massone = mass[type[atom1]];
|
||||
else massone = rmass[atom1];
|
||||
psum[0] += v[atom1][0] * massone;
|
||||
psum[1] += v[atom1][1] * massone;
|
||||
psum[2] += v[atom1][2] * massone;
|
||||
msum += massone;
|
||||
}
|
||||
t[0] += (psum[0]*psum[0] + psum[1]*psum[1] +
|
||||
psum[2]*psum[2]) / msum;
|
||||
|
||||
vav[0] = psum[0]/msum;
|
||||
vav[1] = psum[1]/msum;
|
||||
vav[2] = psum[2]/msum;
|
||||
|
||||
if (mass) massone = mass[type[i]];
|
||||
else massone = rmass[i];
|
||||
t[1] += (
|
||||
(v[i][0]-vav[0])*(v[i][0]-vav[0]) +
|
||||
(v[i][1]-vav[1])*(v[i][1]-vav[1]) +
|
||||
(v[i][2]-vav[2])*(v[i][2]-vav[2])
|
||||
) * massone;
|
||||
for (m = 0; m < num_bond[i]; m++) {
|
||||
atom1 = atom->map(bond_atom[i][m]);
|
||||
if (atom1 == -1) {
|
||||
char str[128];
|
||||
sprintf(str,"Bond atoms %d %d missing",
|
||||
tag[i],bond_atom[i][m]);
|
||||
error->one(str);
|
||||
}
|
||||
if (mass) massone = mass[type[atom1]];
|
||||
else massone = rmass[atom1];
|
||||
t[1] += (
|
||||
(v[atom1][0]-vav[0])*(v[atom1][0]-vav[0]) +
|
||||
(v[atom1][1]-vav[1])*(v[atom1][1]-vav[1]) +
|
||||
(v[atom1][2]-vav[2])*(v[atom1][2]-vav[2])
|
||||
) * massone;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
MPI_Allreduce(t,vector,2,MPI_DOUBLE,MPI_SUM,world);
|
||||
vector[0] *= tfactor;
|
||||
vector[1] *= tfactor1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeTempCoM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
double **v = atom->v;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = v[j][0];
|
||||
buf[m++] = v[j][1];
|
||||
buf[m++] = v[j][2];
|
||||
}
|
||||
return 3;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeTempCoM::unpack_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
double **v = atom->v;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
v[i][0] = buf[m++];
|
||||
v[i][1] = buf[m++];
|
||||
v[i][2] = buf[m++];
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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();
|
||||
int pack_comm(int, int *, double *, int, int *);
|
||||
void unpack_comm(int, int, double *);
|
||||
|
||||
private:
|
||||
int fix_dof;
|
||||
double tfactor;
|
||||
double dof1;
|
||||
double tfactor1;
|
||||
|
||||
void recount();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
@ -90,6 +90,7 @@ CommandStyle(write_restart,WriteRestart)
|
|||
#include "compute_temp.h"
|
||||
#include "compute_temp_deform.h"
|
||||
#include "compute_temp_partial.h"
|
||||
#include "compute_temp_com.h"
|
||||
#include "compute_temp_ramp.h"
|
||||
#include "compute_temp_region.h"
|
||||
#endif
|
||||
|
|
@ -109,6 +110,7 @@ ComputeStyle(sum,ComputeSum)
|
|||
ComputeStyle(temp,ComputeTemp)
|
||||
ComputeStyle(temp/deform,ComputeTempDeform)
|
||||
ComputeStyle(temp/partial,ComputeTempPartial)
|
||||
ComputeStyle(temp/com,ComputeTempCoM)
|
||||
ComputeStyle(temp/ramp,ComputeTempRamp)
|
||||
ComputeStyle(temp/region,ComputeTempRegion)
|
||||
#endif
|
||||
|
|
|
|||
Loading…
Reference in New Issue