diff --git a/doc/Section_intro.txt b/doc/Section_intro.txt index feece7de40..b9774250b4 100644 --- a/doc/Section_intro.txt +++ b/doc/Section_intro.txt @@ -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) diff --git a/doc/compute.txt b/doc/compute.txt index d89bd9e928..e52f05f860 100644 --- a/doc/compute.txt +++ b/doc/compute.txt @@ -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 diff --git a/doc/compute_temp_com.txt b/doc/compute_temp_com.txt new file mode 100644 index 0000000000..62d03cd78b --- /dev/null +++ b/doc/compute_temp_com.txt @@ -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 diff --git a/src/MAKE/Makefile.mac b/src/MAKE/Makefile.mac index e5d695535e..827e80d0ce 100755 --- a/src/MAKE/Makefile.mac +++ b/src/MAKE/Makefile.mac @@ -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 diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp new file mode 100644 index 0000000000..565bed9410 --- /dev/null +++ b/src/compute_temp_com.cpp @@ -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++]; + } +} diff --git a/src/compute_temp_com.h b/src/compute_temp_com.h new file mode 100644 index 0000000000..de81e3c8ca --- /dev/null +++ b/src/compute_temp_com.h @@ -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 diff --git a/src/style.h b/src/style.h index 4be8738257..097c855f83 100644 --- a/src/style.h +++ b/src/style.h @@ -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