diff --git a/src/fix_heat.cpp b/src/fix_heat.cpp new file mode 100644 index 0000000000..4813643689 --- /dev/null +++ b/src/fix_heat.cpp @@ -0,0 +1,91 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_heat.h" +#include "atom.h" +#include "domain.h" +#include "group.h" +#include "force.h" +#include "update.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixHeat::FixHeat(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal fix heat command"); + + nevery = atoi(arg[3]); + if (nevery <= 0) error->all("Illegal fix heat command"); + + heat = atof(arg[4]); + heat *= nevery*update->dt*force->ftm2v; + + // cannot have 0 atoms in group + + if (group->count(igroup) == 0.0) error->all("Fix heat group has no atoms"); +} + +/* ---------------------------------------------------------------------- */ + +int FixHeat::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixHeat::init() +{ + masstotal = group->mass(igroup); +} + +/* ---------------------------------------------------------------------- */ + +void FixHeat::end_of_step() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double vsub[3],vcm[3]; + + double ke = group->ke(igroup); + group->vcm(igroup,masstotal,vcm); + double vcmsq = vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]; + double escale = (ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal); + if (escale < 0.0) error->all("Illegal fix heat attempt"); + double r = sqrt(escale); + + vsub[0] = (r-1.0) * vcm[0]; + vsub[1] = (r-1.0) * vcm[1]; + vsub[2] = (r-1.0) * vcm[2]; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] = r*v[i][0] - vsub[0]; + v[i][1] = r*v[i][1] - vsub[1]; + v[i][2] = r*v[i][2] - vsub[2]; + } +} diff --git a/src/fix_heat.h b/src/fix_heat.h new file mode 100644 index 0000000000..f05e3f3bdc --- /dev/null +++ b/src/fix_heat.h @@ -0,0 +1,35 @@ +/* ---------------------------------------------------------------------- + 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 FIX_HEAT_H +#define FIX_HEAT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixHeat : public Fix { + public: + FixHeat(class LAMMPS *, int, char **); + int setmask(); + void init(); + void end_of_step(); + + private: + double heat; + double masstotal; +}; + +} + +#endif diff --git a/src/group.cpp b/src/group.cpp index 9c369101ad..a59d2e6821 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -20,6 +20,7 @@ #include "domain.h" #include "atom.h" #include "atom_vec_hybrid.h" +#include "force.h" #include "region.h" #include "memory.h" #include "error.h" @@ -467,19 +468,19 @@ double Group::mass(int igroup) double *rmass = atom->rmass; int nlocal = atom->nlocal; - double m = 0.0; + double one = 0.0; if (mass) { for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) m += mass[type[i]]; + if (mask[i] & groupbit) one += mass[type[i]]; } else { for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) m += rmass[i]; + if (mask[i] & groupbit) one += rmass[i]; } - double mall; - MPI_Allreduce(&m,&mall,1,MPI_DOUBLE,MPI_SUM,world); - return mall; + double all; + MPI_Allreduce(&one,&all,1,MPI_DOUBLE,MPI_SUM,world); + return all; } /* ---------------------------------------------------------------------- @@ -671,6 +672,41 @@ void Group::fcm(int igroup, double *cm) MPI_Allreduce(flocal,cm,3,MPI_DOUBLE,MPI_SUM,world); } +/* ---------------------------------------------------------------------- + compute the total kinetic energy of group and return it +------------------------------------------------------------------------- */ + +double Group::ke(int igroup) +{ + int groupbit = bitmask[igroup]; + + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + double one = 0.0; + + if (mass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + rmass[i]; + } + + double all; + MPI_Allreduce(&one,&all,1,MPI_DOUBLE,MPI_SUM,world); + all *= 0.5 * force->mvv2e; + return all; +} + /* ---------------------------------------------------------------------- compute the radius-of-gyration of group around center-of-mass cm must unwrap atoms to compute Rg correctly diff --git a/src/group.h b/src/group.h index 7d22d0c38c..34789e5c75 100644 --- a/src/group.h +++ b/src/group.h @@ -42,6 +42,7 @@ class Group : protected Pointers { void xcm(int, double, double *); // center-of-mass coords of group void vcm(int, double, double *); // center-of-mass velocity of group void fcm(int, double *); // total force on group + double ke(int); // kinetic energy of group double gyration(int, double, double *); // radius-of-gyration of group void angmom(int, double *, double *); // angular momentum of group void inertia(int, double *, double [3][3]); // inertia tensor diff --git a/src/style.h b/src/style.h index f2e17b0b7f..1ff730a2a8 100644 --- a/src/style.h +++ b/src/style.h @@ -137,6 +137,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_enforce2d.h" #include "fix_gravity.h" #include "fix_gyration.h" +#include "fix_heat.h" #include "fix_indent.h" #include "fix_langevin.h" #include "fix_line_force.h" @@ -184,6 +185,7 @@ FixStyle(efield,FixEfield) FixStyle(enforce2d,FixEnforce2D) FixStyle(gravity,FixGravity) FixStyle(gyration,FixGyration) +FixStyle(heat,FixHeat) FixStyle(indent,FixIndent) FixStyle(langevin,FixLangevin) FixStyle(lineforce,FixLineForce)