mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10368 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
28a774a70b
commit
2bbf3f1ef7
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
@ -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();
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue