Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88092995
modify.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Oct 16, 18:11
Size
23 KB
Mime Type
text/x-c
Expires
Fri, Oct 18, 18:11 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
21712754
Attached To
rLAMMPS lammps
modify.cpp
View Options
/* ----------------------------------------------------------------------
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 "stdio.h"
#include "string.h"
#include "modify.h"
#include "atom.h"
#include "comm.h"
#include "fix.h"
#include "compute.h"
#include "group.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#define ComputeInclude
#define FixInclude
#include "style.h"
#undef ComputeInclude
#undef FixInclude
using namespace LAMMPS_NS;
#define DELTA 4
// mask settings - same as in fix.cpp
#define INITIAL_INTEGRATE 1
#define PRE_EXCHANGE 2
#define PRE_NEIGHBOR 4
#define POST_FORCE 8
#define FINAL_INTEGRATE 16
#define END_OF_STEP 32
#define THERMO_ENERGY 64
#define INITIAL_INTEGRATE_RESPA 128
#define POST_FORCE_RESPA 256
#define FINAL_INTEGRATE_RESPA 512
#define MIN_POST_FORCE 1024
/* ---------------------------------------------------------------------- */
Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
{
nfix = maxfix = 0;
n_initial_integrate = 0;
n_pre_exchange = n_pre_neighbor = 0;
n_post_force = n_final_integrate = n_end_of_step = n_thermo_energy = 0;
n_initial_integrate_respa = n_post_force_respa = n_final_integrate_respa = 0;
n_min_post_force = 0;
fix = NULL;
fmask = NULL;
list_initial_integrate = NULL;
list_pre_exchange = list_pre_neighbor = NULL;
list_post_force = list_final_integrate = list_end_of_step = NULL;
list_thermo_energy = NULL;
list_initial_integrate_respa = list_post_force_respa = NULL;
list_final_integrate_respa = NULL;
list_min_post_force = NULL;
end_of_step_every = NULL;
nfix_restart_global = 0;
id_restart_global = style_restart_global = state_restart_global = NULL;
nfix_restart_peratom = 0;
id_restart_peratom = style_restart_peratom = NULL;
index_restart_peratom = NULL;
ncompute = maxcompute = 0;
compute = NULL;
}
/* ---------------------------------------------------------------------- */
Modify::~Modify()
{
// delete all fixes
// do it via delete_fix() so callbacks in Atom are also updated correctly
while (nfix) delete_fix(fix[0]->id);
memory->sfree(fix);
memory->sfree(fmask);
delete [] list_initial_integrate;
delete [] list_pre_exchange;
delete [] list_pre_neighbor;
delete [] list_post_force;
delete [] list_final_integrate;
delete [] list_end_of_step;
delete [] list_thermo_energy;
delete [] list_initial_integrate_respa;
delete [] list_post_force_respa;
delete [] list_final_integrate_respa;
delete [] list_min_post_force;
delete [] end_of_step_every;
restart_deallocate();
// delete all computes
for (int i = 0; i < ncompute; i++) delete compute[i];
memory->sfree(compute);
}
/* ----------------------------------------------------------------------
initialize all fixes and computes
------------------------------------------------------------------------- */
void Modify::init()
{
int i;
// delete storage of restart info since it is not valid after 1st run
restart_deallocate();
// init each fix
comm->maxforward_fix = comm->maxreverse_fix = 0;
for (i = 0; i < nfix; i++) fix[i]->init();
// create lists of fixes to call at each stage of run
list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate);
list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
list_init(POST_FORCE,n_post_force,list_post_force);
list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
list_init_thermo_energy(THERMO_ENERGY,n_thermo_energy,list_thermo_energy);
list_init(INITIAL_INTEGRATE_RESPA,
n_initial_integrate_respa,list_initial_integrate_respa);
list_init(POST_FORCE_RESPA,
n_post_force_respa,list_post_force_respa);
list_init(FINAL_INTEGRATE_RESPA,
n_final_integrate_respa,list_final_integrate_respa);
list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
// init each compute
for (i = 0; i < ncompute; i++) compute[i]->init();
}
/* ----------------------------------------------------------------------
setup for run, calls setup() of all fixes
------------------------------------------------------------------------- */
void Modify::setup()
{
if (update->whichflag == 0)
for (int i = 0; i < nfix; i++) fix[i]->setup();
else
for (int i = 0; i < nfix; i++) fix[i]->min_setup();
}
/* ----------------------------------------------------------------------
1st half of integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::initial_integrate()
{
for (int i = 0; i < n_initial_integrate; i++)
fix[list_initial_integrate[i]]->initial_integrate();
}
/* ----------------------------------------------------------------------
pre_exchange call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::pre_exchange()
{
for (int i = 0; i < n_pre_exchange; i++)
fix[list_pre_exchange[i]]->pre_exchange();
}
/* ----------------------------------------------------------------------
pre_neighbor call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::pre_neighbor()
{
for (int i = 0; i < n_pre_neighbor; i++)
fix[list_pre_neighbor[i]]->pre_neighbor();
}
/* ----------------------------------------------------------------------
force adjustment call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_force(int vflag)
{
for (int i = 0; i < n_post_force; i++)
fix[list_post_force[i]]->post_force(vflag);
}
/* ----------------------------------------------------------------------
2nd half of integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::final_integrate()
{
for (int i = 0; i < n_final_integrate; i++)
fix[list_final_integrate[i]]->final_integrate();
}
/* ----------------------------------------------------------------------
end-of-timestep call only for relevant fixes
only call fix->end_of_step() on timesteps that are multiples of nevery
------------------------------------------------------------------------- */
void Modify::end_of_step()
{
for (int i = 0; i < n_end_of_step; i++)
if (update->ntimestep % end_of_step_every[i] == 0)
fix[list_end_of_step[i]]->end_of_step();
}
/* ----------------------------------------------------------------------
thermo energy call only for relevant fixes
called by Thermo clas
arg to Fix thermo() is 0, so fix will return its energy contribution
------------------------------------------------------------------------- */
double Modify::thermo_energy()
{
double energy = 0.0;
for (int i = 0; i < n_thermo_energy; i++)
energy += fix[list_thermo_energy[i]]->thermo(0);
return energy;
}
/* ----------------------------------------------------------------------
1st half of rRESPA integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::initial_integrate_respa(int ilevel, int flag)
{
for (int i = 0; i < n_initial_integrate_respa; i++)
fix[list_initial_integrate_respa[i]]->initial_integrate_respa(ilevel,flag);
}
/* ----------------------------------------------------------------------
rRESPA force adjustment call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_force_respa(int vflag, int ilevel, int iloop)
{
for (int i = 0; i < n_post_force_respa; i++)
fix[list_post_force_respa[i]]->post_force_respa(vflag,ilevel,iloop);
}
/* ----------------------------------------------------------------------
2nd half of rRESPA integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::final_integrate_respa(int ilevel)
{
for (int i = 0; i < n_final_integrate_respa; i++)
fix[list_final_integrate_respa[i]]->final_integrate_respa(ilevel);
}
/* ----------------------------------------------------------------------
minimizer force adjustment call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::min_post_force(int vflag)
{
for (int i = 0; i < n_min_post_force; i++)
fix[list_min_post_force[i]]->min_post_force(vflag);
}
/* ----------------------------------------------------------------------
add a new fix or replace one with same ID
------------------------------------------------------------------------- */
void Modify::add_fix(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all("Fix command before simulation box is defined");
if (narg < 3) error->all("Illegal fix command");
// find group ID
int igroup = group->find(arg[1]);
if (igroup == -1) error->all("Could not find fix group ID");
// if fix ID exists:
// set newflag = 0 so create new fix in same location in fix list
// error if new style does not match old style
// since can't replace it (all when-to-invoke ptrs would be invalid)
// warn if new group != old group
// delete old fix
// set ptr to NULL in case new fix scans list of fixes
// if fix ID does not exist:
// set newflag = 1 so create new fix
// extend fix and fmask lists as necessary
int ifix,newflag;
for (ifix = 0; ifix < nfix; ifix++)
if (strcmp(arg[0],fix[ifix]->id) == 0) break;
if (ifix < nfix) {
newflag = 0;
if (strcmp(arg[2],fix[ifix]->style) != 0)
error->all("Replacing a fix, but new style != old style");
if (fix[ifix]->igroup != igroup && comm->me == 0)
error->warning("Replacing a fix, but new group != old group");
delete fix[ifix];
atom->update_callback(ifix);
fix[ifix] = NULL;
} else {
newflag = 1;
if (nfix == maxfix) {
maxfix += DELTA;
fix = (Fix **) memory->srealloc(fix,maxfix*sizeof(Fix *),"modify:fix");
fmask = (int *)
memory->srealloc(fmask,maxfix*sizeof(int),"modify:fmask");
}
}
// create the Fix
if (0) return; // dummy line to enable else-if macro expansion
#define FixClass
#define FixStyle(key,Class) \
else if (strcmp(arg[2],#key) == 0) fix[ifix] = new Class(lmp,narg,arg);
#include "style.h"
#undef FixClass
else error->all("Invalid fix style");
// if fix is new, set it's mask values and increment nfix
if (newflag) {
fmask[ifix] = fix[ifix]->setmask();
nfix++;
}
// check if Fix is in restart_global list
// if yes, pass state info to the Fix so it can reset itself
for (int i = 0; i < nfix_restart_global; i++)
if (strcmp(id_restart_global[i],fix[ifix]->id) == 0 &&
strcmp(style_restart_global[i],fix[ifix]->style) == 0) {
fix[ifix]->restart(state_restart_global[i]);
if (comm->me == 0) {
char *str = "Resetting global state of Fix %s Style %s "
"from restart file info\n";
if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style);
if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
}
}
// check if Fix is in restart_peratom list
// if yes, loop over atoms so they can extract info from atom->extra array
for (int i = 0; i < nfix_restart_peratom; i++)
if (strcmp(id_restart_peratom[i],fix[ifix]->id) == 0 &&
strcmp(style_restart_peratom[i],fix[ifix]->style) == 0) {
for (int j = 0; j < atom->nlocal; j++)
fix[ifix]->unpack_restart(j,index_restart_peratom[i]);
if (comm->me == 0) {
char *str = "Resetting per-atom state of Fix %s Style %s "
"from restart file info\n";
if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style);
if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
}
}
}
/* ----------------------------------------------------------------------
modify a Fix's parameters
------------------------------------------------------------------------- */
void Modify::modify_fix(int narg, char **arg)
{
if (narg < 2) error->all("Illegal fix_modify command");
// lookup Fix ID
int ifix;
for (ifix = 0; ifix < nfix; ifix++)
if (strcmp(arg[0],fix[ifix]->id) == 0) break;
if (ifix == nfix) error->all("Could not find fix_modify ID");
fix[ifix]->modify_params(narg-1,&arg[1]);
}
/* ----------------------------------------------------------------------
delete a Fix from list of Fixes
Atom class must update indices in its list of callbacks to fixes
------------------------------------------------------------------------- */
void Modify::delete_fix(char *id)
{
int ifix = find_fix(id);
if (ifix < 0) error->all("Could not find fix ID to delete");
delete fix[ifix];
atom->update_callback(ifix);
// move other Fixes and fmask down in list one slot
for (int i = ifix+1; i < nfix; i++) fix[i-1] = fix[i];
for (int i = ifix+1; i < nfix; i++) fmask[i-1] = fmask[i];
nfix--;
}
/* ----------------------------------------------------------------------
find a fix by ID
return index of fix or -1 if not found
------------------------------------------------------------------------- */
int Modify::find_fix(char *id)
{
int ifix;
for (ifix = 0; ifix < nfix; ifix++)
if (strcmp(id,fix[ifix]->id) == 0) break;
if (ifix == nfix) return -1;
return ifix;
}
/* ----------------------------------------------------------------------
add a new compute
------------------------------------------------------------------------- */
void Modify::add_compute(int narg, char **arg)
{
if (narg < 3) error->all("Illegal compute command");
// error checks
for (int icompute = 0; icompute < ncompute; icompute++)
if (strcmp(arg[0],compute[icompute]->id) == 0)
error->all("Reuse of compute ID");
int igroup = group->find(arg[1]);
if (igroup == -1) error->all("Could not find compute group ID");
// extend Compute list if necessary
if (ncompute == maxcompute) {
maxcompute += DELTA;
compute = (Compute **)
memory->srealloc(compute,maxcompute*sizeof(Compute *),"modify:compute");
}
// create the Compute
if (0) return; // dummy line to enable else-if macro expansion
#define ComputeClass
#define ComputeStyle(key,Class) \
else if (strcmp(arg[2],#key) == 0) \
compute[ncompute] = new Class(lmp,narg,arg);
#include "style.h"
else error->all("Invalid compute style");
ncompute++;
}
/* ----------------------------------------------------------------------
modify a Compute's parameters
------------------------------------------------------------------------- */
void Modify::modify_compute(int narg, char **arg)
{
if (narg < 2) error->all("Illegal compute_modify command");
// lookup Compute ID
int icompute;
for (icompute = 0; icompute < ncompute; icompute++)
if (strcmp(arg[0],compute[icompute]->id) == 0) break;
if (icompute == ncompute) error->all("Could not find compute_modify ID");
compute[icompute]->modify_params(narg-1,&arg[1]);
}
/* ----------------------------------------------------------------------
delete a Compute from list of Computes
------------------------------------------------------------------------- */
void Modify::delete_compute(char *id)
{
int icompute = find_compute(id);
if (icompute < 0) error->all("Could not find compute ID to delete");
delete compute[icompute];
// move other Computes down in list one slot
for (int i = icompute+1; i < ncompute; i++) compute[i-1] = compute[i];
ncompute--;
}
/* ----------------------------------------------------------------------
find a compute by ID
return index of compute or -1 if not found
------------------------------------------------------------------------- */
int Modify::find_compute(char *id)
{
int icompute;
for (icompute = 0; icompute < ncompute; icompute++)
if (strcmp(id,compute[icompute]->id) == 0) break;
if (icompute == ncompute) return -1;
return icompute;
}
/* ----------------------------------------------------------------------
write to restart file for all Fixes with restart info
(1) fixes that have global state
(2) fixes that store per-atom quantities
------------------------------------------------------------------------- */
void Modify::write_restart(FILE *fp)
{
int me = comm->me;
int count = 0;
for (int i = 0; i < nfix; i++)
if (fix[i]->restart_global) count++;
if (me == 0) fwrite(&count,sizeof(int),1,fp);
int n;
for (int i = 0; i < nfix; i++)
if (fix[i]->restart_global) {
if (me == 0) {
n = strlen(fix[i]->id) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(fix[i]->id,sizeof(char),n,fp);
n = strlen(fix[i]->style) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(fix[i]->style,sizeof(char),n,fp);
}
fix[i]->write_restart(fp);
}
count = 0;
for (int i = 0; i < nfix; i++)
if (fix[i]->restart_peratom) count++;
if (me == 0) fwrite(&count,sizeof(int),1,fp);
for (int i = 0; i < nfix; i++)
if (fix[i]->restart_peratom) {
if (me == 0) {
n = strlen(fix[i]->id) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(fix[i]->id,sizeof(char),n,fp);
n = strlen(fix[i]->style) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(fix[i]->style,sizeof(char),n,fp);
n = fix[i]->maxsize_restart();
fwrite(&n,sizeof(int),1,fp);
}
}
}
/* ----------------------------------------------------------------------
read in restart file data on all previously defined Fixes with restart info
(1) fixes that have global state
(2) fixes that store per-atom quantities
return maxsize of extra info that will be stored with any atom
------------------------------------------------------------------------- */
int Modify::read_restart(FILE *fp)
{
// nfix_restart_global = # of restart entries with global state info
int me = comm->me;
if (me == 0) fread(&nfix_restart_global,sizeof(int),1,fp);
MPI_Bcast(&nfix_restart_global,1,MPI_INT,0,world);
// allocate space for each entry
if (nfix_restart_global) {
id_restart_global = new char*[nfix_restart_global];
style_restart_global = new char*[nfix_restart_global];
state_restart_global = new char*[nfix_restart_global];
}
// read each entry and Bcast to all procs
// each entry has id string, style string, chunk of state data
int n;
for (int i = 0; i < nfix_restart_global; i++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
id_restart_global[i] = new char[n];
if (me == 0) fread(id_restart_global[i],sizeof(char),n,fp);
MPI_Bcast(id_restart_global[i],n,MPI_CHAR,0,world);
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
style_restart_global[i] = new char[n];
if (me == 0) fread(style_restart_global[i],sizeof(char),n,fp);
MPI_Bcast(style_restart_global[i],n,MPI_CHAR,0,world);
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
state_restart_global[i] = new char[n];
if (me == 0) fread(state_restart_global[i],sizeof(char),n,fp);
MPI_Bcast(state_restart_global[i],n,MPI_CHAR,0,world);
}
// nfix_restart_peratom = # of restart entries with peratom info
int maxsize = 0;
if (me == 0) fread(&nfix_restart_peratom,sizeof(int),1,fp);
MPI_Bcast(&nfix_restart_peratom,1,MPI_INT,0,world);
// allocate space for each entry
if (nfix_restart_peratom) {
id_restart_peratom = new char*[nfix_restart_peratom];
style_restart_peratom = new char*[nfix_restart_peratom];
index_restart_peratom = new int[nfix_restart_peratom];
}
// read each entry and Bcast to all procs
// each entry has id string, style string, maxsize of one atom's data
// set index = which set of extra data this fix represents
for (int i = 0; i < nfix_restart_peratom; i++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
id_restart_peratom[i] = new char[n];
if (me == 0) fread(id_restart_peratom[i],sizeof(char),n,fp);
MPI_Bcast(id_restart_peratom[i],n,MPI_CHAR,0,world);
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
style_restart_peratom[i] = new char[n];
if (me == 0) fread(style_restart_peratom[i],sizeof(char),n,fp);
MPI_Bcast(style_restart_peratom[i],n,MPI_CHAR,0,world);
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
maxsize += n;
index_restart_peratom[i] = i;
}
return maxsize;
}
/* ----------------------------------------------------------------------
delete all lists of restart file Fix info
------------------------------------------------------------------------- */
void Modify::restart_deallocate()
{
if (nfix_restart_global) {
for (int i = 0; i < nfix_restart_global; i++) {
delete [] id_restart_global[i];
delete [] style_restart_global[i];
delete [] state_restart_global[i];
}
delete [] id_restart_global;
delete [] style_restart_global;
delete [] state_restart_global;
}
if (nfix_restart_peratom) {
for (int i = 0; i < nfix_restart_peratom; i++) {
delete [] id_restart_peratom[i];
delete [] style_restart_peratom[i];
}
delete [] id_restart_peratom;
delete [] style_restart_peratom;
delete [] index_restart_peratom;
}
nfix_restart_global = nfix_restart_peratom = 0;
}
/* ----------------------------------------------------------------------
create list of fix indices for fixes which match mask
------------------------------------------------------------------------- */
void Modify::list_init(int mask, int &n, int *&list)
{
delete [] list;
n = 0;
for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
list = new int[n];
n = 0;
for (int i = 0; i < nfix; i++) if (fmask[i] & mask) list[n++] = i;
}
/* ----------------------------------------------------------------------
create list of fix indices for end_of_step fixes
also create end_of_step_every[]
------------------------------------------------------------------------- */
void Modify::list_init_end_of_step(int mask, int &n, int *&list)
{
delete [] list;
delete [] end_of_step_every;
n = 0;
for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
list = new int[n];
end_of_step_every = new int[n];
n = 0;
for (int i = 0; i < nfix; i++)
if (fmask[i] & mask) {
list[n] = i;
end_of_step_every[n++] = fix[i]->nevery;
}
}
/* ----------------------------------------------------------------------
create list of fix indices for thermo energy fixes
only added to list if its thermo_energy flag was set via fix_modify
------------------------------------------------------------------------- */
void Modify::list_init_thermo_energy(int mask, int &n, int *&list)
{
delete [] list;
n = 0;
for (int i = 0; i < nfix; i++)
if (fmask[i] & mask && fix[i]->thermo_energy) n++;
list = new int[n];
n = 0;
for (int i = 0; i < nfix; i++)
if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory from all fixes
------------------------------------------------------------------------- */
int Modify::memory_usage()
{
int bytes = 0;
for (int i = 0; i < nfix; i++) bytes += fix[i]->memory_usage();
for (int i = 0; i < ncompute; i++) bytes += compute[i]->memory_usage();
return bytes;
}
Event Timeline
Log In to Comment