From 4615a859c3cfb19b113dbdff222e7842270ee2e9 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 14 Jun 2016 22:36:31 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15174 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_spring_chunk.cpp | 249 +++++++++++++++++++++++++++++++++++++++ src/fix_spring_chunk.h | 78 ++++++++++++ 2 files changed, 327 insertions(+) create mode 100644 src/fix_spring_chunk.cpp create mode 100644 src/fix_spring_chunk.h diff --git a/src/fix_spring_chunk.cpp b/src/fix_spring_chunk.cpp new file mode 100644 index 0000000000..d6dc36a73c --- /dev/null +++ b/src/fix_spring_chunk.cpp @@ -0,0 +1,249 @@ +/* ---------------------------------------------------------------------- + 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 +#include +#include +#include "fix_spring_chunk.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "respa.h" +#include "domain.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "compute_com_chunk.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define SMALL 1.0e-10 + +/* ---------------------------------------------------------------------- */ + +FixSpringChunk::FixSpringChunk(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 6) error->all(FLERR,"Illegal fix spring/chunk command"); + + scalar_flag = 1; + global_freq = 1; + extscalar = 1; + + k_spring = force->numeric(FLERR,arg[3]); + + int n = strlen(arg[4]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[4]); + + n = strlen(arg[5]) + 1; + idcom = new char[n]; + strcpy(idcom,arg[5]); + + esprings = 0.0; + com0 = fcom = NULL; + nchunk = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixSpringChunk::~FixSpringChunk() +{ + memory->destroy(com0); + memory->destroy(fcom); + + // decrement lock counter in compute chunk/atom, it if still exists + + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->unlock(this); + cchunk->lockcount--; + } + + delete [] idchunk; + delete [] idcom; +} + +/* ---------------------------------------------------------------------- */ + +int FixSpringChunk::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::init() +{ + // current indices for idchunk and idcom + + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for fix spring/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Fix spring/chunk does not use chunk/atom compute"); + + icompute = modify->find_compute(idcom); + if (icompute < 0) + error->all(FLERR,"Com/chunk compute does not exist for fix spring/chunk"); + ccom = (ComputeCOMChunk *) modify->compute[icompute]; + if (strcmp(ccom->style,"com/chunk") != 0) + error->all(FLERR,"Fix spring/chunk does not use com/chunk compute"); + + // check that idchunk is consistent with ccom->idchunk + + if (strcmp(idchunk,ccom->idchunk) != 0) + error->all(FLERR,"Fix spring chunk chunkID not same as comID chunkID"); + + if (strstr(update->integrate_style,"respa")) + nlevels_respa = ((Respa *) update->integrate)->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::min_setup(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::post_force(int vflag) +{ + int i,m; + double dx,dy,dz,r; + + // check if first time cchunk will be queried via ccom + // if so, lock idchunk for as long as this fix is in place + // will be unlocked in destructor + // necessary b/c this fix stores original COM + + if (com0 == NULL) cchunk->lock(this,update->ntimestep,-1); + + // calculate current centers of mass for each chunk + // extract pointers from idchunk and idcom + + ccom->compute_array(); + + nchunk = cchunk->nchunk; + int *ichunk = cchunk->ichunk; + double *masstotal = ccom->masstotal; + double **com = ccom->array; + + // check if first time cchunk was queried via ccom + // if so, allocate com0,fcom and store initial COM + + if (com0 == NULL) { + memory->create(com0,nchunk,3,"spring/chunk:com0"); + memory->create(fcom,nchunk,3,"spring/chunk:fcom"); + + for (m = 0; m < nchunk; m++) { + com0[m][0] = com[m][0]; + com0[m][1] = com[m][1]; + com0[m][2] = com[m][2]; + } + } + + // calculate fcom = force on each COM, divided by masstotal + + esprings = 0.0; + for (m = 0; m < nchunk; m++) { + dx = com[m][0] - com0[m][0]; + dy = com[m][1] - com0[m][1]; + dz = com[m][2] - com0[m][2]; + r = sqrt(dx*dx + dy*dy + dz*dz); + r = MAX(r,SMALL); + + if (masstotal[m]) { + fcom[m][0] = k_spring*dx/r / masstotal[m]; + fcom[m][1] = k_spring*dy/r / masstotal[m]; + fcom[m][2] = k_spring*dz/r / masstotal[m]; + esprings += 0.5*k_spring*r*r; + } + } + + // apply restoring force to atoms in each chunk + + double **f = atom->f; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + double massone; + + if (rmass) { + for (i = 0; i < nlocal; i++) { + m = ichunk[i]-1; + if (m < 0) continue; + massone = rmass[i]; + f[i][0] -= fcom[m][0]*massone; + f[i][1] -= fcom[m][1]*massone; + f[i][2] -= fcom[m][2]*massone; + } + } else { + for (i = 0; i < nlocal; i++) { + m = ichunk[i]-1; + if (m < 0) continue; + massone = mass[type[i]]; + f[i][0] -= fcom[m][0]*massone; + f[i][1] -= fcom[m][1]*massone; + f[i][2] -= fcom[m][2]*massone; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixSpringChunk::min_post_force(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- + energy of stretched spring +------------------------------------------------------------------------- */ + +double FixSpringChunk::compute_scalar() +{ + return esprings; +} diff --git a/src/fix_spring_chunk.h b/src/fix_spring_chunk.h new file mode 100644 index 0000000000..c95d9063ca --- /dev/null +++ b/src/fix_spring_chunk.h @@ -0,0 +1,78 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(spring/chunk,FixSpringChunk) + +#else + +#ifndef LMP_FIX_SPRING_CHUNK_H +#define LMP_FIX_SPRING_CHUNK_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixSpringChunk : public Fix { + public: + FixSpringChunk(class LAMMPS *, int, char **); + ~FixSpringChunk(); + int setmask(); + void init(); + void setup(int); + void min_setup(int); + void post_force(int); + void post_force_respa(int, int, int); + void min_post_force(int); + double compute_scalar(); + + private: + int nlevels_respa; + double k_spring; + double esprings; + char *idchunk,*idcom; + + int nchunk; + double **com0,**fcom; + + class ComputeChunkAtom *cchunk; + class ComputeCOMChunk *ccom; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: R0 < 0 for fix spring command + +Equilibrium spring length is invalid. + +E: Fix spring couple group ID does not exist + +Self-explanatory. + +E: Two groups cannot be the same in fix spring couple + +Self-explanatory. + +*/