diff --git a/src/MISC/compute_msd_nongauss.cpp b/src/MISC/compute_msd_nongauss.cpp new file mode 100644 index 0000000000..eb038cc356 --- /dev/null +++ b/src/MISC/compute_msd_nongauss.cpp @@ -0,0 +1,106 @@ +/* ---------------------------------------------------------------------- + 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 authors: Rob Hoy +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_msd_nongauss.h" +#include "atom.h" +#include "update.h" +#include "group.h" +#include "domain.h" +#include "fix_store.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeMSDNonGauss::ComputeMSDNonGauss(LAMMPS *lmp, int narg, char **arg) : + ComputeMSD(lmp, narg, arg) +{ + size_vector = 3; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMSDNonGauss::compute_vector() +{ + invoked_vector = update->ntimestep; + + // cm = current center of mass + + double cm[3]; + if (comflag) group->xcm(igroup,masstotal,cm); + else cm[0] = cm[1] = cm[2] = 0.0; + + // dx,dy,dz = displacement of atom from original position + // original unwrapped position is stored by fix + // relative to center of mass if comflag is set + // for triclinic, need to unwrap current atom coord via h matrix + + double **xoriginal = fix->astore; + + double **x = atom->x; + int *mask = atom->mask; + tagint *image = atom->image; + int nlocal = atom->nlocal; + + double *h = domain->h; + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + double dx,dy,dz; + int xbox,ybox,zbox; + + double msd[2]; + msd[0] = msd[1] = 0.0; + + if (domain->triclinic == 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + dx = x[i][0] + xbox*xprd - cm[0] - xoriginal[i][0]; + dy = x[i][1] + ybox*yprd - cm[1] - xoriginal[i][1]; + dz = x[i][2] + zbox*zprd - cm[2] - xoriginal[i][2]; + msd[0] += dx*dx + dy*dy + dz*dz; + msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz); + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - + cm[0] - xoriginal[i][0]; + dy = x[i][1] + h[1]*ybox + h[3]*zbox - cm[1] - xoriginal[i][1]; + dz = x[i][2] + h[2]*zbox - cm[2] - xoriginal[i][2]; + msd[0] += dx*dx + dy*dy + dz*dz; + msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz); + } + } + + MPI_Allreduce(msd,vector,2,MPI_DOUBLE,MPI_SUM,world); + if (nmsd) { + vector[0] /= nmsd; + vector[1] /= nmsd; + vector[2] = (3*vector[1])/(5*vector[0]*vector[0]) - 1; + } +} diff --git a/src/MISC/compute_msd_nongauss.h b/src/MISC/compute_msd_nongauss.h new file mode 100644 index 0000000000..86bb653c56 --- /dev/null +++ b/src/MISC/compute_msd_nongauss.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + 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 COMPUTE_CLASS + +ComputeStyle(msd/nongauss,ComputeMSDNonGauss) + +#else + +#ifndef LMP_COMPUTE_MSD_NONGAUSS_H +#define LMP_COMPUTE_MSD_NONGAUSS_H + +#include "compute_msd.h" + +namespace LAMMPS_NS { + +class ComputeMSDNonGauss : public ComputeMSD { + public: + ComputeMSDNonGauss(class LAMMPS *, int, char **); + ~ComputeMSDNonGauss() {} + void compute_vector(); +}; + +} + +#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: Could not find compute msd fix ID + +Self-explanatory. + +*/ diff --git a/src/MISC/compute_ti.cpp b/src/MISC/compute_ti.cpp new file mode 100644 index 0000000000..e47bbda556 --- /dev/null +++ b/src/MISC/compute_ti.cpp @@ -0,0 +1,241 @@ +/* ---------------------------------------------------------------------- + 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: Sai Jayaraman (University of Notre Dame) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "atom.h" +#include "string.h" +#include "compute_ti.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "pair.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{PAIR,TAIL,KSPACE}; + +/* ---------------------------------------------------------------------- */ + +ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all(FLERR,"Illegal compute ti command"); + + peflag = 1; + peratom_flag = 1; + peatomflag = 1; + scalar_flag = 1; + extscalar = 1; + timeflag = 1; + + // terms come in triplets + // changed to quadruplets to include atom type + + nterms = (narg-3) / 4; + if (narg != 4*nterms + 3) error->all(FLERR,"Illegal compute ti command"); + + which = new int[nterms]; + ivar1 = new int[nterms]; + ivar2 = new int[nterms]; + ilo = new int[nterms]; + ihi = new int[nterms]; + var1 = new char*[nterms]; + var2 = new char*[nterms]; + pptr = new Pair*[nterms]; + pstyle = new char*[nterms]; + + for (int m = 0; m < nterms; m++) pstyle[m] = NULL; + + // parse keywords + + nterms = 0; + + int iarg = 3; + while (iarg < narg) { + if (iarg+4 > narg) error->all(FLERR,"Illegal compute ti command"); + if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE; + else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL; + else which[nterms] = PAIR; + + int n = strlen(arg[iarg]) + 1; + pstyle[nterms] = new char[n]; + strcpy(pstyle[nterms],arg[iarg]); + force->bounds(arg[iarg+1],atom->ntypes,ilo[nterms],ihi[nterms]); + iarg += 1; + + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + int n = strlen(&arg[iarg+1][2]) + 1; + var1[nterms] = new char[n]; + strcpy(var1[nterms],&arg[iarg+1][2]); + } else error->all(FLERR,"Illegal compute ti command"); + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + var2[nterms] = new char[n]; + strcpy(var2[nterms],&arg[iarg+2][2]); + } else error->all(FLERR,"Illegal compute ti command"); + + nterms++; + iarg += 3; + } +} + +/* --------------------------------------------------------------------- */ + +ComputeTI::~ComputeTI() +{ + for (int m = 0; m < nterms; m++) { + delete [] var1[m]; + delete [] var2[m]; + delete [] pstyle[m]; + } + delete [] which; + delete [] ivar1; + delete [] ivar2; + delete [] var1; + delete [] var2; + delete [] ilo; + delete [] ihi; + delete [] pptr; + delete [] pstyle; +} + +/* --------------------------------------------------------------------- */ + +void ComputeTI::init() +{ + // setup and error checks + + for (int m = 0; m < nterms; m++) { + ivar1[m] = input->variable->find(var1[m]); + ivar2[m] = input->variable->find(var2[m]); + if (ivar1[m] < 0 || ivar2 < 0) + error->all(FLERR,"Variable name for compute ti does not exist"); + if (!input->variable->equalstyle(ivar1[m]) || + !input->variable->equalstyle(ivar2[m])) + error->all(FLERR,"Variable for compute ti is invalid style"); + + if (which[m] == PAIR) { + pptr[m] = force->pair_match(pstyle[m],1); + if (pptr[m] == NULL) + error->all(FLERR,"Compute ti pair style does not exist"); + + } else if (which[m] == TAIL) { + if (force->pair == NULL || force->pair->tail_flag == 0) + error->all(FLERR,"Compute ti tail when pair style does not " + "compute tail corrections"); + + } else if (which[m] == KSPACE) { + if (force->kspace == NULL) + error->all(FLERR,"Compute ti kspace style does not exist"); + } + } +} + +/* --------------------------------------------------------------------- */ + +double ComputeTI::compute_scalar() +{ + double eng,engall,value1,value2; + + invoked_scalar = update->ntimestep; + if (update->eflag_global != invoked_scalar) + error->all(FLERR,"Energy was not tallied on needed timestep"); + + double dUdl = 0.0; + + for (int m = 0; m < nterms; m++) { + int total_flag = 0; + if ((ihi[m]-ilo[m])==atom->ntypes) total_flag = 1; + eng = 0.0; + value1 = input->variable->compute_equal(ivar1[m]); + value2 = input->variable->compute_equal(ivar2[m]); + if (value1 == 0.0) continue; + + if (which[m] == PAIR) { + int ntypes = atom->ntypes; + int *mask = atom->mask; + if (total_flag) { + eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + } + else { + int nlocal = atom->nlocal; + int npair = nlocal; + int *mask = atom->mask; + int *type = atom->type; + + double *eatom = pptr[m]->eatom; + + if (force->newton) npair += atom->nghost; + for (int i = 0; i < npair; i++) + if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) eng += eatom[i]; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + } + dUdl += engall/value1 * value2; + + } else if (which[m] == TAIL) { + double vol = domain->xprd*domain->yprd*domain->zprd; + if (total_flag) + eng = force->pair->etail / vol; + else { + eng = 0; + for (int it = 1; it <= atom->ntypes; it++) { + int jt; + if ((it >= ilo[m])&&(it <=ihi[m])) jt = it; + else jt = ilo[m]; + for (; jt <=ihi[m];jt++) { + if ((force->pair->tail_flag)&&(force->pair->setflag[it][jt])) { + double cut = force->pair->init_one(it,jt); + eng += force->pair->etail_ij; + } + if (it !=jt) eng += force->pair->etail_ij; + } + } + eng /= vol; + } + dUdl += eng/value1 * value2; + + } else if (which[m] == KSPACE) { + int ntypes = atom->ntypes; + int *mask = atom->mask; + if (total_flag) + eng = force->kspace->energy; + else { + int nlocal = atom->nlocal; + int npair = nlocal; + int *mask = atom->mask; + int *type = atom->type; + double *eatom = force->kspace->eatom; + eng = 0; + for(int i = 0; i < nlocal; i++) + if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) + eng += eatom[i]; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + eng = engall; + } + dUdl += eng/value1 * value2; + } + } + + scalar = dUdl; + return scalar; +} diff --git a/src/MISC/compute_ti.h b/src/MISC/compute_ti.h new file mode 100644 index 0000000000..23d3e3a4e9 --- /dev/null +++ b/src/MISC/compute_ti.h @@ -0,0 +1,83 @@ +/* ---------------------------------------------------------------------- + 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 COMPUTE_CLASS + +ComputeStyle(ti,ComputeTI) + +#else + +#ifndef COMPUTE_TI_H +#define COMPUTE_TI_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTI : public Compute { + public: + ComputeTI(class LAMMPS *, int, char **); + ~ComputeTI(); + void init(); + double compute_scalar(); + + private: + int nterms; + int *which; + int *ivar1,*ivar2; + int *ilo, *ihi; + char **var1,**var2; + class Pair **pptr; + char **pstyle; +}; + +} + +#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: Variable name for compute ti does not exist + +Self-explanatory. + +E: Variable for compute ti is invalid style + +Self-explanatory. + +E: Compute ti pair style does not exist + +Self-explanatory. + +E: Compute ti tail when pair style does not compute tail corrections + +Self-explanatory. + +E: Compute ti kspace style does not exist + +Self-explanatory. + +E: Energy was not tallied on needed timestep + +You are using a thermo keyword that requires potentials to +have tallied energy, but they didn't on this timestep. See the +variable doc page for ideas on how to make this work. + +*/ diff --git a/src/MISC/fix_adapt.cpp b/src/MISC/fix_adapt.cpp new file mode 100644 index 0000000000..2f1793dfba --- /dev/null +++ b/src/MISC/fix_adapt.cpp @@ -0,0 +1,389 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "string.h" +#include "stdlib.h" +#include "fix_adapt.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "force.h" +#include "pair.h" +#include "pair_hybrid.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +enum{PAIR,KSPACE,ATOM}; +enum{DIAMETER,CHARGE}; + +/* ---------------------------------------------------------------------- */ + +FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); + nevery = force->inumeric(FLERR,arg[3]); + if (nevery < 0) error->all(FLERR,"Illegal fix adapt command"); + + // count # of adaptations + + nadapt = 0; + + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"kspace") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 2; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 3; + } else break; + } + + if (nadapt == 0) error->all(FLERR,"Illegal fix adapt command"); + adapt = new Adapt[nadapt]; + + // parse keywords + + nadapt = 0; + diamflag = 0; + chgflag = 0; + + iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); + adapt[nadapt].which = PAIR; + int n = strlen(arg[iarg+1]) + 1; + adapt[nadapt].pstyle = new char[n]; + strcpy(adapt[nadapt].pstyle,arg[iarg+1]); + n = strlen(arg[iarg+2]) + 1; + adapt[nadapt].pparam = new char[n]; + strcpy(adapt[nadapt].pparam,arg[iarg+2]); + force->bounds(arg[iarg+3],atom->ntypes, + adapt[nadapt].ilo,adapt[nadapt].ihi); + force->bounds(arg[iarg+4],atom->ntypes, + adapt[nadapt].jlo,adapt[nadapt].jhi); + if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) { + n = strlen(&arg[iarg+5][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+5][2]); + } else error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"kspace") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + adapt[nadapt].which = KSPACE; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + int n = strlen(&arg[iarg+1][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+1][2]); + } else error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 2; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); + adapt[nadapt].which = ATOM; + if (strcmp(arg[iarg+1],"diameter") == 0) { + adapt[nadapt].aparam = DIAMETER; + diamflag = 1; + } else if (strcmp(arg[iarg+1],"charge") == 0) { + adapt[nadapt].aparam = CHARGE; + chgflag = 1; + } else error->all(FLERR,"Illegal fix adapt command"); + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+2][2]); + } else error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 3; + } else break; + } + + // optional keywords + + resetflag = 0; + scaleflag = 0; + + while (iarg < narg) { + if (strcmp(arg[iarg],"reset") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1; + else error->all(FLERR,"Illegal fix adapt command"); + iarg += 2; + } else if (strcmp(arg[iarg],"scale") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix adapt command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix adapt command"); + } + + // allocate pair style arrays + + int n = atom->ntypes; + for (int m = 0; m < nadapt; m++) + if (adapt[m].which == PAIR) + memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); +} + +/* ---------------------------------------------------------------------- */ + +FixAdapt::~FixAdapt() +{ + for (int m = 0; m < nadapt; m++) { + delete [] adapt[m].var; + if (adapt[m].which == PAIR) { + delete [] adapt[m].pstyle; + delete [] adapt[m].pparam; + memory->destroy(adapt[m].array_orig); + } + } + delete [] adapt; +} + +/* ---------------------------------------------------------------------- */ + +int FixAdapt::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + mask |= POST_RUN; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::init() +{ + int i,j; + + // setup and error checks + + anypair = 0; + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + + ad->ivar = input->variable->find(ad->var); + if (ad->ivar < 0) + error->all(FLERR,"Variable name for fix adapt does not exist"); + if (!input->variable->equalstyle(ad->ivar)) + error->all(FLERR,"Variable for fix adapt is invalid style"); + + if (ad->which == PAIR) { + anypair = 1; + + Pair *pair = force->pair_match(ad->pstyle,1); + if (pair == NULL) error->all(FLERR,"Fix adapt pair style does not exist"); + void *ptr = pair->extract(ad->pparam,ad->pdim); + if (ptr == NULL) + error->all(FLERR,"Fix adapt pair style param not supported"); + + ad->pdim = 2; + if (ad->pdim == 0) ad->scalar = (double *) ptr; + if (ad->pdim == 2) ad->array = (double **) ptr; + + // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style + + if (ad->pdim == 2 && (strcmp(force->pair_style,"hybrid") == 0 || + strcmp(force->pair_style,"hybrid/overlay") == 0)) { + PairHybrid *pair = (PairHybrid *) force->pair; + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + if (!pair->check_ijtype(i,j,ad->pstyle)) + error->all(FLERR,"Fix adapt type pair range is not valid for " + "pair hybrid sub-style"); + } + + } else if (ad->which == KSPACE) { + if (force->kspace == NULL) + error->all(FLERR,"Fix adapt kspace style does not exist"); + kspace_scale = (double *) force->kspace->extract("scale"); + + } else if (ad->which == ATOM) { + if (ad->aparam == DIAMETER) { + if (!atom->radius_flag) + error->all(FLERR,"Fix adapt requires atom attribute diameter"); + } + if (ad->aparam == CHARGE) { + if (!atom->q_flag) + error->all(FLERR,"Fix adapt requires atom attribute charge"); + } + } + } + + // make copy of original pair array values + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + if (ad->which == PAIR && ad->pdim == 2) { + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array_orig[i][j] = ad->array[i][j]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::setup_pre_force(int vflag) +{ + change_settings(); +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::pre_force(int vflag) +{ + if (nevery == 0) return; + if (update->ntimestep % nevery) return; + change_settings(); +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::post_run() +{ + if (resetflag) restore_settings(); +} + +/* ---------------------------------------------------------------------- + change pair,kspace,atom parameters based on variable evaluation +------------------------------------------------------------------------- */ + +void FixAdapt::change_settings() +{ + int i,j; + + // variable evaluation may invoke computes so wrap with clear/add + + modify->clearstep_compute(); + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + double value = input->variable->compute_equal(ad->ivar); + + // set global scalar or type pair array values + + if (ad->which == PAIR) { + if (ad->pdim == 0) { + if (scaleflag) *ad->scalar = value * ad->scalar_orig; + else *ad->scalar = value; + } else if (ad->pdim == 2) { + if (scaleflag) + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = value*ad->array_orig[i][j]; + else + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = value; + } + + // set kspace scale factor + + } else if (ad->which == KSPACE) { + *kspace_scale = value; + + } else if (ad->which == ATOM) { + + // set radius from diameter + // also scale rmass to new value + + if (ad->aparam == DIAMETER) { + int mflag = 0; + if (atom->rmass_flag) mflag = 1; + double density; + + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (mflag == 0) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + radius[i] = 0.5*value; + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + density = rmass[i] / (4.0*MY_PI/3.0 * + radius[i]*radius[i]*radius[i]); + radius[i] = 0.5*value; + rmass[i] = 4.0*MY_PI/3.0 * + radius[i]*radius[i]*radius[i] * density; + } + } + } else if (ad->aparam == CHARGE) { + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) q[i] = value; + } + } + } + + modify->addstep_compute(update->ntimestep + nevery); + + // re-initialize pair styles if any PAIR settings were changed + // this resets other coeffs that may depend on changed values, + // and also offset and tail corrections + + if (anypair) force->pair->reinit(); +} + +/* ---------------------------------------------------------------------- + restore pair,kspace.atom parameters to original values +------------------------------------------------------------------------- */ + +void FixAdapt::restore_settings() +{ + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + if (ad->which == PAIR) { + if (ad->pdim == 0) *ad->scalar = ad->scalar_orig; + else if (ad->pdim == 2) { + for (int i = ad->ilo; i <= ad->ihi; i++) + for (int j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = ad->array_orig[i][j]; + } + + } else if (ad->which == KSPACE) { + *kspace_scale = 1.0; + + } else if (ad->which == ATOM) { + + } + } + + if (anypair) force->pair->reinit(); +} diff --git a/src/MISC/fix_adapt.h b/src/MISC/fix_adapt.h new file mode 100644 index 0000000000..88d94d169a --- /dev/null +++ b/src/MISC/fix_adapt.h @@ -0,0 +1,107 @@ +/* ---------------------------------------------------------------------- + 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(adapt,FixAdapt) + +#else + +#ifndef LMP_FIX_ADAPT_H +#define LMP_FIX_ADAPT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAdapt : public Fix { + public: + int diamflag; // 1 if atom diameters will vary, for AtomVecGranular + int chgflag; + + FixAdapt(class LAMMPS *, int, char **); + ~FixAdapt(); + int setmask(); + void init(); + void setup_pre_force(int); + void pre_force(int); + void post_run(); + + private: + int nadapt,resetflag,scaleflag; + int anypair; + + struct Adapt { + int which,ivar; + char *var; + char *pstyle,*pparam; + int ilo,ihi,jlo,jhi; + int pdim; + double *scalar,scalar_orig; + double **array,**array_orig; + int aparam; + }; + + Adapt *adapt; + double *kspace_scale; + + void change_settings(); + void restore_settings(); +}; + +} + +#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: Variable name for fix adapt does not exist + +Self-explanatory. + +E: Variable for fix adapt is invalid style + +Only equal-style variables can be used. + +E: Fix adapt pair style does not exist + +Self-explanatory + +E: Fix adapt pair style param not supported + +The pair style does not know about the parameter you specified. + +E: Fix adapt type pair range is not valid for pair hybrid sub-style + +Self-explanatory. + +E: Fix adapt kspace style does not exist + +Self-explanatory. + +E: Fix adapt requires atom attribute diameter + +The atom style being used does not specify an atom diameter. + +E: Fix adapt requires atom attribute charge + +The atom style being used does not specify an atom charge. + +*/