ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
e12585148971a33f33cbd622616fcc5528119aa7
lammps/src/respa.cpp

631 lines
20 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
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: Mark Stevens (SNL), Paul Crozier (SNL)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "respa.h"
#include "neighbor.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "output.h"
#include "update.h"
#include "modify.h"
#include "fix_respa.h"
#include "memory.h"
#include "timer.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
{
if (narg < 1) error->all("Illegal run_style respa command");
nlevels = atoi(arg[0]);
if (nlevels < 1) error->all("Respa levels must be >= 1");
if (narg < nlevels) error->all("Illegal run_style respa command");
loop = new int[nlevels];
for (int iarg = 1; iarg < nlevels; iarg++) {
loop[iarg-1] = atoi(arg[iarg]);
if (loop[iarg-1] <= 0) error->all("Illegal run_style respa command");
}
loop[nlevels-1] = 1;
// set level at which each force is computed
// argument settings override defaults
level_bond = level_angle = level_dihedral = level_improper = -1;
level_pair = level_kspace = -1;
level_inner = level_middle = level_outer = -1;
int iarg = nlevels;
while (iarg < narg) {
if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_bond = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"angle") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_angle = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_dihedral = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"improper") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_improper = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_pair = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"inner") == 0) {
if (iarg+4 > narg) error->all("Illegal run_style respa command");
level_inner = atoi(arg[iarg+1]) - 1;
cutoff[0] = atof(arg[iarg+2]);
cutoff[1] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"middle") == 0) {
if (iarg+4 > narg) error->all("Illegal run_style respa command");
level_middle = atoi(arg[iarg+1]) - 1;
cutoff[2] = atof(arg[iarg+2]);
cutoff[3] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"outer") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_outer = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all("Illegal run_style respa command");
level_kspace = atoi(arg[iarg+1]) - 1;
iarg += 2;
} else error->all("Illegal run_style respa command");
}
// cannot specify both pair and inner/middle/outer
if (level_pair >= 0 &&
(level_inner >= 0 || level_middle >= 0 || level_outer >= 0))
error->all("Cannot set both respa pair and inner/middle/outer");
// if either inner and outer is specified, then both must be
if ((level_inner >= 0 && level_outer == -1) ||
(level_outer >= 0 && level_inner == -1))
error->all("Must set both respa inner and outer");
// middle cannot be set without inner/outer
if (level_middle >= 0 && level_inner == -1)
error->all("Cannot set respa middle without inner/outer");
// set defaults if user did not specify level
// bond to innermost level
// angle same as bond, dihedral same as angle, improper same as dihedral
// pair to outermost level if no inner/middle/outer
// inner/middle/outer have no defaults
// kspace same as pair or outer
if (level_bond == -1) level_bond = 0;
if (level_angle == -1) level_angle = level_bond;
if (level_dihedral == -1) level_dihedral = level_angle;
if (level_improper == -1) level_improper = level_dihedral;
if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1;
if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair;
if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer;
// print respa levels
if (comm->me == 0) {
if (screen) {
fprintf(screen,"Respa levels:\n");
for (int i = 0; i < nlevels; i++) {
fprintf(screen," %d =",i);
if (level_bond == i) fprintf(screen," bond");
if (level_angle == i) fprintf(screen," angle");
if (level_dihedral == i) fprintf(screen," dihedral");
if (level_improper == i) fprintf(screen," improper");
if (level_pair == i) fprintf(screen," pair");
if (level_inner == i) fprintf(screen," pair-inner");
if (level_middle == i) fprintf(screen," pair-middle");
if (level_outer == i) fprintf(screen," pair-outer");
if (level_kspace == i) fprintf(screen," kspace");
fprintf(screen,"\n");
}
}
if (logfile) {
fprintf(logfile,"Respa levels:\n");
for (int i = 0; i < nlevels; i++) {
fprintf(logfile," %d =",i);
if (level_bond == i) fprintf(logfile," bond");
if (level_angle == i) fprintf(logfile," angle");
if (level_dihedral == i) fprintf(logfile," dihedral");
if (level_improper == i) fprintf(logfile," improper");
if (level_pair == i) fprintf(logfile," pair");
if (level_inner == i) fprintf(logfile," pair-inner");
if (level_middle == i) fprintf(logfile," pair-middle");
if (level_outer == i) fprintf(logfile," pair-outer");
if (level_kspace == i) fprintf(logfile," kspace");
fprintf(logfile,"\n");
}
}
}
// check that levels are in correct order
if (level_angle < level_bond || level_dihedral < level_angle ||
level_improper < level_dihedral)
error->all("Invalid order of forces within respa levels");
if (level_pair >= 0) {
if (level_pair < level_improper || level_kspace < level_pair)
error->all("Invalid order of forces within respa levels");
}
if (level_pair == -1 && level_middle == -1) {
if (level_inner < level_improper || level_outer < level_inner ||
level_kspace != level_outer)
error->all("Invalid order of forces within respa levels");
}
if (level_pair == -1 && level_middle >= 0) {
if (level_inner < level_improper || level_middle < level_inner ||
level_outer < level_inner || level_kspace != level_outer)
error->all("Invalid order of forces within respa levels");
}
// warn if any levels are devoid of forces
int flag = 0;
for (int i = 0; i < nlevels; i++)
if (level_bond != i && level_angle != i && level_dihedral != i &&
level_improper != i && level_pair != i && level_inner != i &&
level_middle != i && level_outer != i && level_kspace != i) flag = 1;
if (flag && comm->me == 0)
error->warning("One or more respa levels compute no forces");
// check cutoff consistency if inner/middle/outer are enabled
if (level_inner >= 0 && cutoff[1] < cutoff[0])
error->all("Respa inner cutoffs are invalid");
if (level_middle >= 0 && (cutoff[3] < cutoff[2] || cutoff[2] < cutoff[1]))
error->all("Respa middle cutoffs are invalid");
// set outer pair of cutoffs to inner pair if middle is not enabled
if (level_inner >= 0 && level_middle < 0) {
cutoff[2] = cutoff[0];
cutoff[3] = cutoff[1];
}
// allocate other needed level arrays
newton = new int[nlevels];
step = new double[nlevels];
}
/* ---------------------------------------------------------------------- */
Respa::~Respa()
{
delete [] loop;
delete [] newton;
delete [] step;
}
/* ----------------------------------------------------------------------
initialization before run
------------------------------------------------------------------------- */
void Respa::init()
{
// warn if no fixes
if (modify->nfix == 0)
error->warning("No fixes defined, atoms won't move");
// respa not allowed with granular atom style
if (atom->check_style("granular"))
error->all("Respa not allowed with atom style granular");
// create fix needed for storing atom-based respa level forces
// will delete it at end of run
char **fixarg = new char*[4];
fixarg[0] = "RESPA";
fixarg[1] = "all";
fixarg[2] = "RESPA";
fixarg[3] = new char[8];
sprintf(fixarg[3],"%d",nlevels);
modify->add_fix(4,fixarg);
delete [] fixarg[3];
delete [] fixarg;
fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
// insure respa inner/middle/outer is using Pair class that supports it
if (level_inner >= 0)
if (force->pair && force->pair->respa_enable == 0)
error->all("Pair style does not support rRESPA inner/middle/outer");
// setup virial computations for timestepping
// virial_style = 1 (explicit) since never computed implicitly like Verlet
// virial_every = 1 if computed every timestep (NPT,NPH)
// fix arrays store info on fixes that need virial computed occasionally
virial_style = 1;
virial_every = 0;
nfix_virial = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->pressure_every == 1) virial_every = 1;
else if (modify->fix[i]->pressure_every > 1) nfix_virial++;
if (nfix_virial) {
delete [] fix_virial_every;
delete [] next_fix_virial;
fix_virial_every = new int[nfix_virial];
next_fix_virial = new int[nfix_virial];
nfix_virial = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->pressure_every > 1)
fix_virial_every[nfix_virial++] = modify->fix[i]->pressure_every;
}
// step[] = timestep for each level
step[nlevels-1] = update->dt;
for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
step[ilevel] = step[ilevel+1]/loop[ilevel];
// set newton flag for each level
for (int ilevel = 0; ilevel < nlevels; ilevel++) {
newton[ilevel] = 0;
if (force->newton_bond) {
if (level_bond == ilevel || level_angle == ilevel ||
level_dihedral == ilevel || level_improper == ilevel)
newton[ilevel] = 1;
}
if (force->newton_pair) {
if (level_pair == ilevel || level_inner == ilevel ||
level_middle == ilevel || level_outer == ilevel)
newton[ilevel] = 1;
}
}
}
/* ----------------------------------------------------------------------
setup before run
------------------------------------------------------------------------- */
void Respa::setup()
{
if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
neighbor->build();
neighbor->ncalls = 0;
// compute all forces
eflag = 1;
vflag = virial_style;
for (int ilevel = 0; ilevel < nlevels; ilevel++) {
force_clear(newton[ilevel]);
if (level_bond == ilevel && force->bond)
force->bond->compute(eflag,vflag);
if (level_angle == ilevel && force->angle)
force->angle->compute(eflag,vflag);
if (level_dihedral == ilevel && force->dihedral)
force->dihedral->compute(eflag,vflag);
if (level_improper == ilevel && force->improper)
force->improper->compute(eflag,vflag);
if (level_pair == ilevel && force->pair)
force->pair->compute(eflag,vflag);
if (level_inner == ilevel && force->pair)
force->pair->compute_inner();
if (level_middle == ilevel && force->pair)
force->pair->compute_middle();
if (level_outer == ilevel && force->pair)
force->pair->compute_outer(eflag,vflag);
if (level_kspace == ilevel && force->kspace) {
force->kspace->setup();
force->kspace->compute(eflag,vflag);
}
if (newton[ilevel]) comm->reverse_communicate();
copy_f_flevel(ilevel);
}
modify->setup();
sum_flevel_f();
output->setup(1);
// setup virial computations for timestepping
int ntimestep = update->ntimestep;
next_virial = 0;
if (virial_every) next_virial = ntimestep + 1;
else {
for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
next_fix_virial[ivirial] =
(ntimestep/fix_virial_every[ivirial])*fix_virial_every[ivirial] +
fix_virial_every[ivirial];
if (ivirial) next_virial = MIN(next_virial,next_fix_virial[ivirial]);
else next_virial = next_fix_virial[0];
}
}
}
/* ----------------------------------------------------------------------
iterate for n steps
------------------------------------------------------------------------- */
void Respa::iterate(int n)
{
int ntimestep;
for (int i = 0; i < n; i++) {
ntimestep = ++update->ntimestep;
// eflag/vflag = 0/1 for energy/virial computation
if (ntimestep == output->next_thermo) eflag = 1;
else eflag = 0;
if (ntimestep == output->next_thermo || ntimestep == next_virial) {
vflag = virial_style;
if (virial_every) next_virial++;
else next_virial = fix_virial(ntimestep);
} else vflag = 0;
recurse(nlevels-1);
if (modify->n_end_of_step) modify->end_of_step();
if (ntimestep == output->next) {
timer->stamp();
sum_flevel_f();
output->write(update->ntimestep);
timer->stamp(TIME_OUTPUT);
}
}
}
/* ----------------------------------------------------------------------
delete fix at end of run, so its atom arrays won't persist
------------------------------------------------------------------------- */
void Respa::cleanup()
{
modify->delete_fix("RESPA");
}
/* ---------------------------------------------------------------------- */
void Respa::recurse(int ilevel)
{
copy_flevel_f(ilevel);
for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
modify->initial_integrate_respa(ilevel,0);
if (ilevel) recurse(ilevel-1);
// at outermost level:
// call initial_integrate w/ flag = 1 so fix NPT,NPH perform 2nd dilate
// at correct symmetric position in rRESPA hierarchy
// check on rebuilding neighbor list
// at innermost level, communicate
// at middle levels, do nothing
if (ilevel == nlevels-1) {
modify->initial_integrate_respa(ilevel,1);
int nflag = neighbor->decide();
if (nflag) {
if (modify->n_pre_exchange) modify->pre_exchange();
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
comm->borders();
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
}
} else if (ilevel == 0) {
timer->stamp();
comm->communicate();
timer->stamp(TIME_COMM);
}
force_clear(newton[ilevel]);
timer->stamp();
if (level_bond == ilevel && force->bond) {
force->bond->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_angle == ilevel && force->angle) {
force->angle->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_dihedral == ilevel && force->dihedral) {
force->dihedral->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_improper == ilevel && force->improper) {
force->improper->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_pair == ilevel && force->pair) {
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (level_inner == ilevel && force->pair) {
force->pair->compute_inner();
timer->stamp(TIME_PAIR);
}
if (level_middle == ilevel && force->pair) {
force->pair->compute_middle();
timer->stamp(TIME_PAIR);
}
if (level_outer == ilevel && force->pair) {
force->pair->compute_outer(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (level_kspace == ilevel && force->kspace) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
if (newton[ilevel]) {
comm->reverse_communicate();
timer->stamp(TIME_COMM);
}
if (modify->n_post_force_respa) modify->post_force_respa(vflag,ilevel,iloop);
modify->final_integrate_respa(ilevel);
}
copy_f_flevel(ilevel);
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
------------------------------------------------------------------------- */
void Respa::force_clear(int newtonflag)
{
// clear global force array
// nall includes ghosts only if newton flag is set
int nall;
if (newtonflag) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
double **f = atom->f;
for (int i = 0; i < nall; i++) {
f[i][0] = 0.0;
f[i][1] = 0.0;
f[i][2] = 0.0;
}
}
/* ----------------------------------------------------------------------
copy force components from atom->f to FixRespa->f_level
------------------------------------------------------------------------- */
void Respa::copy_f_flevel(int ilevel)
{
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int i = 0; i < n; i++) {
f_level[i][ilevel][0] = f[i][0];
f_level[i][ilevel][1] = f[i][1];
f_level[i][ilevel][2] = f[i][2];
}
}
/* ----------------------------------------------------------------------
copy force components from FixRespa->f_level to atom->f
------------------------------------------------------------------------- */
void Respa::copy_flevel_f(int ilevel)
{
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int i = 0; i < n; i++) {
f[i][0] = f_level[i][ilevel][0];
f[i][1] = f_level[i][ilevel][1];
f[i][2] = f_level[i][ilevel][2];
}
}
/* ----------------------------------------------------------------------
sum all force components from FixRespa->f_level to create full atom->f
------------------------------------------------------------------------- */
void Respa::sum_flevel_f()
{
copy_flevel_f(0);
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int ilevel = 1; ilevel < nlevels; ilevel++) {
for (int i = 0; i < n; i++) {
f[i][0] += f_level[i][ilevel][0];
f[i][1] += f_level[i][ilevel][1];
f[i][2] += f_level[i][ilevel][2];
}
}
}
/* ----------------------------------------------------------------------
return next timestep virial should be computed
based on one or more fixes that need virial computed periodically
------------------------------------------------------------------------- */
int Respa::fix_virial(int ntimestep)
{
int next;
for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
if (ntimestep == next_fix_virial[ivirial])
next_fix_virial[ivirial] += fix_virial_every[ivirial];
if (ivirial) next = MIN(next,next_fix_virial[ivirial]);
else next = next_fix_virial[0];
}
return next;
}
Reference in New Issue View Git Blame Copy Permalink