git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3251 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
@ -14,7 +14,7 @@ OBJ = $(SRC:.cpp=.o)
|
|||||||
# Package variables
|
# Package variables
|
||||||
|
|
||||||
PACKAGE = asphere class2 colloid dipole dpd gpu granular \
|
PACKAGE = asphere class2 colloid dipole dpd gpu granular \
|
||||||
kspace manybody meam molecule opt peri poems reax xtc
|
kspace manybody meam molecule opt peri poems prd reax xtc
|
||||||
|
|
||||||
PACKUSER = user-ackland user-atc user-cg-cmm user-ewaldn user-smd
|
PACKUSER = user-ackland user-atc user-cg-cmm user-ewaldn user-smd
|
||||||
|
|
||||||
|
|||||||
@ -1,9 +1,9 @@
|
|||||||
# Settings for libraries used by specific LAMMPS packages
|
# Settings for libraries used by specific LAMMPS packages
|
||||||
# this file is auto-edited when those packages are included/excluded
|
# this file is auto-edited when those packages are included/excluded
|
||||||
|
|
||||||
PKG_INC = -I../../lib/atc -I../../lib/reax -I../../lib/poems -I../../lib/meam
|
PKG_INC = -I../../lib/reax -I../../lib/poems -I../../lib/meam
|
||||||
PKG_PATH = -L../../lib/atc -L../../lib/reax -L../../lib/poems -L../../lib/meam -L../../lib/gpu
|
PKG_PATH = -L../../lib/reax -L../../lib/poems -L../../lib/meam
|
||||||
PKG_LIB = -latc -lreax -lpoems -lmeam -lgpu
|
PKG_LIB = -lreax -lpoems -lmeam
|
||||||
|
|
||||||
PKG_SYSPATH = $(user-atc_SYSPATH) $(reax_SYSPATH) $(meam_SYSPATH) $(gpu_SYSPATH)
|
PKG_SYSPATH = $(reax_SYSPATH) $(meam_SYSPATH)
|
||||||
PKG_SYSLIB = $(user-atc_SYSLIB) $(reax_SYSLIB) $(meam_SYSLIB) $(gpu_SYSLIB)
|
PKG_SYSLIB = $(reax_SYSLIB) $(meam_SYSLIB)
|
||||||
|
|||||||
28
src/PRD/Install.csh
Normal file
28
src/PRD/Install.csh
Normal file
@ -0,0 +1,28 @@
|
|||||||
|
# Install/unInstall package classes in LAMMPS
|
||||||
|
|
||||||
|
if ($1 == 1) then
|
||||||
|
|
||||||
|
cp style_prd.h ..
|
||||||
|
|
||||||
|
cp compute_event_displace.cpp ..
|
||||||
|
cp fix_event.cpp ..
|
||||||
|
cp prd.cpp ..
|
||||||
|
|
||||||
|
cp compute_event_displace.h ..
|
||||||
|
cp fix_event.h ..
|
||||||
|
cp prd.h ..
|
||||||
|
|
||||||
|
else if ($1 == 0) then
|
||||||
|
|
||||||
|
rm ../style_prd.h
|
||||||
|
touch ../style_prd.h
|
||||||
|
|
||||||
|
rm ../compute_event_displace.cpp
|
||||||
|
rm ../fix_event.cpp
|
||||||
|
rm ../prd.cpp
|
||||||
|
|
||||||
|
rm ../compute_event_displace.h
|
||||||
|
rm ../fix_event.h
|
||||||
|
rm ../prd.h
|
||||||
|
|
||||||
|
endif
|
||||||
150
src/PRD/compute_event_displace.cpp
Normal file
150
src/PRD/compute_event_displace.cpp
Normal file
@ -0,0 +1,150 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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: Mike Brown (SNL)
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "math.h"
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "compute_event_displace.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "domain.h"
|
||||||
|
#include "modify.h"
|
||||||
|
#include "fix.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
#include "update.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
|
||||||
|
#define INVOKED_SCALAR 1
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
ComputeEventDisplace::ComputeEventDisplace(LAMMPS *lmp, int narg, char **arg) :
|
||||||
|
Compute(lmp, narg, arg)
|
||||||
|
{
|
||||||
|
if (narg != 4) error->all("Illegal compute event/displace command");
|
||||||
|
|
||||||
|
scalar_flag = 1;
|
||||||
|
extscalar = 0;
|
||||||
|
|
||||||
|
double displace_dist = atof(arg[3]);
|
||||||
|
if (displace_dist <= 0.0)
|
||||||
|
error->all("Distnace must be > 0 for compute event/displace");
|
||||||
|
displace_distsq = displace_dist * displace_dist;
|
||||||
|
|
||||||
|
// fix event ID will be set later by PRD
|
||||||
|
|
||||||
|
id_event = NULL;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
ComputeEventDisplace::~ComputeEventDisplace()
|
||||||
|
{
|
||||||
|
delete [] id_event;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void ComputeEventDisplace::init()
|
||||||
|
{
|
||||||
|
// set fix which stores original atom coords
|
||||||
|
// check if is correct style
|
||||||
|
|
||||||
|
if (id_event == NULL)
|
||||||
|
error->all("Compute event/displace has not had fix event assigned");
|
||||||
|
|
||||||
|
int ifix = modify->find_fix(id_event);
|
||||||
|
if (ifix < 0) error->all("Could not find compute event/displace fix ID");
|
||||||
|
fix = modify->fix[ifix];
|
||||||
|
|
||||||
|
if (strcmp(fix->style,"EVENT") != 0)
|
||||||
|
error->all("Compute event/displace has invalid fix event assigned");
|
||||||
|
|
||||||
|
triclinic = domain->triclinic;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
return non-zero if an atom has moved > displace_dist since last event
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double ComputeEventDisplace::compute_scalar()
|
||||||
|
{
|
||||||
|
invoked_scalar = update->ntimestep;
|
||||||
|
|
||||||
|
double event = 0.0;
|
||||||
|
double **xevent = fix->vector_atom;
|
||||||
|
|
||||||
|
double **x = atom->x;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int *image = atom->image;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
double *h = domain->h;
|
||||||
|
double xprd = domain->xprd;
|
||||||
|
double yprd = domain->yprd;
|
||||||
|
double zprd = domain->zprd;
|
||||||
|
int xbox,ybox,zbox;
|
||||||
|
double dx,dy,dz,rsq;
|
||||||
|
|
||||||
|
if (triclinic == 0) {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
xbox = (image[i] & 1023) - 512;
|
||||||
|
ybox = (image[i] >> 10 & 1023) - 512;
|
||||||
|
zbox = (image[i] >> 20) - 512;
|
||||||
|
dx = x[i][0] + xbox*xprd - xevent[i][0];
|
||||||
|
dy = x[i][1] + ybox*yprd - xevent[i][1];
|
||||||
|
dz = x[i][2] + zbox*zprd - xevent[i][2];
|
||||||
|
rsq = dx*dx + dy*dy + dz*dz;
|
||||||
|
if (rsq >= displace_distsq) {
|
||||||
|
event = 1.0;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
} else {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
xbox = (image[i] & 1023) - 512;
|
||||||
|
ybox = (image[i] >> 10 & 1023) - 512;
|
||||||
|
zbox = (image[i] >> 20) - 512;
|
||||||
|
dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xevent[i][0];
|
||||||
|
dy = x[i][1] + h[1]*ybox + h[3]*zbox - xevent[i][1];
|
||||||
|
dz = x[i][2] + h[2]*zbox - xevent[i][2];
|
||||||
|
rsq = dx*dx + dy*dy + dz*dz;
|
||||||
|
if (rsq >= displace_distsq) {
|
||||||
|
event = 1.0;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
MPI_Allreduce(&event,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
|
||||||
|
return scalar;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void ComputeEventDisplace::reset_extra_compute_fix(char *id_new)
|
||||||
|
{
|
||||||
|
delete [] id_event;
|
||||||
|
int n = strlen(id_new) + 1;
|
||||||
|
id_event = new char[n];
|
||||||
|
strcpy(id_event,id_new);
|
||||||
|
}
|
||||||
38
src/PRD/compute_event_displace.h
Normal file
38
src/PRD/compute_event_displace.h
Normal file
@ -0,0 +1,38 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifndef COMPUTE_EVENT_DISPLACE_H
|
||||||
|
#define COMPUTE_EVENT_DISPLACE_H
|
||||||
|
|
||||||
|
#include "compute.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class ComputeEventDisplace : public Compute {
|
||||||
|
public:
|
||||||
|
ComputeEventDisplace(class LAMMPS *, int, char **);
|
||||||
|
~ComputeEventDisplace();
|
||||||
|
void init();
|
||||||
|
double compute_scalar();
|
||||||
|
void reset_extra_compute_fix(char *);
|
||||||
|
|
||||||
|
private:
|
||||||
|
int triclinic;
|
||||||
|
double displace_distsq;
|
||||||
|
char *id_event;
|
||||||
|
class Fix *fix;
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
249
src/PRD/fix_event.cpp
Normal file
249
src/PRD/fix_event.cpp
Normal file
@ -0,0 +1,249 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||||
|
http://lammps.sandia.gov, Sandia National Laboratories
|
||||||
|
Steve Plimpton, sjplimp@sandia.gov
|
||||||
|
|
||||||
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||||
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||||
|
certain rights in this software. This software is distributed under
|
||||||
|
the GNU General Public License.
|
||||||
|
|
||||||
|
See the README file in the top-level LAMMPS directory.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Contributing author: Mike Brown (SNL)
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "fix_event.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "update.h"
|
||||||
|
#include "domain.h"
|
||||||
|
#include "neighbor.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "universe.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) :
|
||||||
|
Fix(lmp, narg, arg)
|
||||||
|
{
|
||||||
|
if (narg != 3) error->all("Illegal fix event command");
|
||||||
|
|
||||||
|
restart_global = 1;
|
||||||
|
|
||||||
|
// perform initial allocation of atom-based array
|
||||||
|
// register with Atom class
|
||||||
|
|
||||||
|
xevent = NULL;
|
||||||
|
xold = NULL;
|
||||||
|
imageold = NULL;
|
||||||
|
grow_arrays(atom->nmax);
|
||||||
|
atom->add_callback(0);
|
||||||
|
|
||||||
|
event_number = 0;
|
||||||
|
event_timestep = 0;
|
||||||
|
clock = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
FixEvent::~FixEvent()
|
||||||
|
{
|
||||||
|
// unregister callbacks to this fix from Atom class
|
||||||
|
|
||||||
|
atom->delete_callback(id,0);
|
||||||
|
|
||||||
|
// delete locally stored array
|
||||||
|
|
||||||
|
memory->destroy_2d_double_array(xevent);
|
||||||
|
memory->destroy_2d_double_array(xold);
|
||||||
|
memory->sfree(imageold);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixEvent::setmask()
|
||||||
|
{
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
save current atom coords as an event
|
||||||
|
called when an event occurs in some replica
|
||||||
|
set event_timestep = when event occurred in a particular replica
|
||||||
|
update clock = elapsed time since last event, across all replicas
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::store_event(int timestep, int delta_clock)
|
||||||
|
{
|
||||||
|
double **x = atom->x;
|
||||||
|
int *image = atom->image;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
domain->unmap(x[i],image[i],xevent[i]);
|
||||||
|
|
||||||
|
event_timestep = timestep;
|
||||||
|
clock += delta_clock;
|
||||||
|
event_number++;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
store state of all atoms
|
||||||
|
called before quench and subsequent check for event
|
||||||
|
so can later restore pre-quench state if no event occurs
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::store_state()
|
||||||
|
{
|
||||||
|
double **x = atom->x;
|
||||||
|
double **f = atom->f;
|
||||||
|
int *image = atom->image;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
xold[i][0] = x[i][0];
|
||||||
|
xold[i][1] = x[i][1];
|
||||||
|
xold[i][2] = x[i][2];
|
||||||
|
imageold[i] = image[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
restore state of all atoms to pre-quench state
|
||||||
|
called after no event detected so can continue
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::restore_state()
|
||||||
|
{
|
||||||
|
double **x = atom->x;
|
||||||
|
int *image = atom->image;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
x[i][0] = xold[i][0];
|
||||||
|
x[i][1] = xold[i][1];
|
||||||
|
x[i][2] = xold[i][2];
|
||||||
|
image[i] = imageold[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
memory usage of local atom-based array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double FixEvent::memory_usage()
|
||||||
|
{
|
||||||
|
double bytes = 6*atom->nmax * sizeof(double);
|
||||||
|
bytes += atom->nmax*sizeof(int);
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
allocate atom-based array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::grow_arrays(int nmax)
|
||||||
|
{
|
||||||
|
xevent = memory->grow_2d_double_array(xevent,nmax,3,"event:xevent");
|
||||||
|
xold = memory->grow_2d_double_array(xold,nmax,3,"event:xold");
|
||||||
|
imageold = (int *)
|
||||||
|
memory->srealloc(imageold,nmax*sizeof(int),"event:imageold");
|
||||||
|
|
||||||
|
// allow compute event to access stored event coords
|
||||||
|
|
||||||
|
vector_atom = xevent;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
copy values within local atom-based array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::copy_arrays(int i, int j)
|
||||||
|
{
|
||||||
|
xevent[j][0] = xevent[i][0];
|
||||||
|
xevent[j][1] = xevent[i][1];
|
||||||
|
xevent[j][2] = xevent[i][2];
|
||||||
|
xold[j][0] = xold[i][0];
|
||||||
|
xold[j][1] = xold[i][1];
|
||||||
|
xold[j][2] = xold[i][2];
|
||||||
|
imageold[j] = imageold[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
pack values in local atom-based array for exchange with another proc
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixEvent::pack_exchange(int i, double *buf)
|
||||||
|
{
|
||||||
|
buf[0] = xevent[i][0];
|
||||||
|
buf[1] = xevent[i][1];
|
||||||
|
buf[2] = xevent[i][2];
|
||||||
|
buf[3] = xold[i][0];
|
||||||
|
buf[4] = xold[i][1];
|
||||||
|
buf[5] = xold[i][2];
|
||||||
|
buf[6] = imageold[i];
|
||||||
|
|
||||||
|
return 7;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
unpack values in local atom-based array from exchange with another proc
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixEvent::unpack_exchange(int nlocal, double *buf)
|
||||||
|
{
|
||||||
|
xevent[nlocal][0] = buf[0];
|
||||||
|
xevent[nlocal][1] = buf[1];
|
||||||
|
xevent[nlocal][2] = buf[2];
|
||||||
|
xold[nlocal][0] = buf[3];
|
||||||
|
xold[nlocal][1] = buf[4];
|
||||||
|
xold[nlocal][2] = buf[5];
|
||||||
|
imageold[nlocal] = static_cast<int>(buf[6]);
|
||||||
|
|
||||||
|
return 7;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
pack entire state of Fix into one write
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::write_restart(FILE *fp)
|
||||||
|
{
|
||||||
|
int n = 0;
|
||||||
|
double list[5];
|
||||||
|
list[n++] = event_number;
|
||||||
|
list[n++] = event_timestep;
|
||||||
|
list[n++] = clock;
|
||||||
|
list[n++] = replica_number;
|
||||||
|
list[n++] = correlated_event;
|
||||||
|
|
||||||
|
if (comm->me == 0) {
|
||||||
|
int size = n * sizeof(double);
|
||||||
|
fwrite(&size,sizeof(int),1,fp);
|
||||||
|
fwrite(&list,sizeof(double),n,fp);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
use state info from restart file to restart the Fix
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixEvent::restart(char *buf)
|
||||||
|
{
|
||||||
|
int n = 0;
|
||||||
|
double *list = (double *) buf;
|
||||||
|
|
||||||
|
event_number = static_cast<int> (list[n++]);
|
||||||
|
event_timestep = static_cast<int> (list[n++]);
|
||||||
|
clock = static_cast<int> (list[n++]);
|
||||||
|
replica_number = static_cast<int> (list[n++]);
|
||||||
|
correlated_event = static_cast<int> (list[n++]);
|
||||||
|
}
|
||||||
55
src/PRD/fix_event.h
Normal file
55
src/PRD/fix_event.h
Normal file
@ -0,0 +1,55 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifndef FIX_EVENT_H
|
||||||
|
#define FIX_EVENT_H
|
||||||
|
|
||||||
|
#include "fix.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class FixEvent : public Fix {
|
||||||
|
public:
|
||||||
|
int event_number; // event counter
|
||||||
|
int event_timestep; // timestep of last event on any replica
|
||||||
|
int clock; // total elapsed timesteps across all replicas
|
||||||
|
int replica_number; // replica where last event occured
|
||||||
|
int correlated_event; // 1 if last event was correlated, 0 otherwise
|
||||||
|
|
||||||
|
FixEvent(class LAMMPS *, int, char **);
|
||||||
|
~FixEvent();
|
||||||
|
int setmask();
|
||||||
|
|
||||||
|
double memory_usage();
|
||||||
|
void grow_arrays(int);
|
||||||
|
void copy_arrays(int, int);
|
||||||
|
int pack_exchange(int, double *);
|
||||||
|
int unpack_exchange(int, double *);
|
||||||
|
void write_restart(FILE *);
|
||||||
|
void restart(char *);
|
||||||
|
|
||||||
|
// methods specific to FixEvent, invoked by PRD
|
||||||
|
|
||||||
|
void store_event(int, int);
|
||||||
|
void store_state();
|
||||||
|
void restore_state();
|
||||||
|
|
||||||
|
private:
|
||||||
|
double **xevent; // atom coords at last event
|
||||||
|
double **xold; // atom coords for reset/restore
|
||||||
|
int *imageold; // image flags for reset/restore
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
772
src/PRD/prd.cpp
Normal file
772
src/PRD/prd.cpp
Normal file
@ -0,0 +1,772 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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: Mike Brown (SNL)
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "math.h"
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "prd.h"
|
||||||
|
#include "universe.h"
|
||||||
|
#include "update.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "domain.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "velocity.h"
|
||||||
|
#include "integrate.h"
|
||||||
|
#include "min.h"
|
||||||
|
#include "neighbor.h"
|
||||||
|
#include "modify.h"
|
||||||
|
#include "compute.h"
|
||||||
|
#include "fix.h"
|
||||||
|
#include "fix_event.h"
|
||||||
|
#include "force.h"
|
||||||
|
#include "pair.h"
|
||||||
|
#include "random_park.h"
|
||||||
|
#include "random_mars.h"
|
||||||
|
#include "output.h"
|
||||||
|
#include "dump.h"
|
||||||
|
#include "finish.h"
|
||||||
|
#include "timer.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
|
||||||
|
#define MAXINT 0x7FFFFFFF
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
perform PRD
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::command(int narg, char **arg)
|
||||||
|
{
|
||||||
|
int i,flag,allflag,ireplica;
|
||||||
|
|
||||||
|
// error checks
|
||||||
|
|
||||||
|
if (domain->box_exist == 0)
|
||||||
|
error->all("PRD command before simulation box is defined");
|
||||||
|
if (universe->nworlds != universe->nprocs &&
|
||||||
|
atom->map_style == 0)
|
||||||
|
error->all("Cannot use PRD command with multi-proc replicas "
|
||||||
|
"unless atom map exists");
|
||||||
|
if (universe->nworlds == 1 && comm->me == 0)
|
||||||
|
error->warning("Running PRD with only one replica");
|
||||||
|
|
||||||
|
if (narg < 7) error->universe_all("Illegal prd command");
|
||||||
|
|
||||||
|
nsteps = atoi(arg[0]);
|
||||||
|
t_event = atoi(arg[1]);
|
||||||
|
n_dephase = atoi(arg[2]);
|
||||||
|
t_dephase = atoi(arg[3]);
|
||||||
|
t_corr = atoi(arg[4]);
|
||||||
|
|
||||||
|
char id_compute[strlen(arg[5])+1];
|
||||||
|
strcpy(id_compute,arg[5]);
|
||||||
|
int seed = atoi(arg[6]);
|
||||||
|
|
||||||
|
options(narg-7,&arg[7]);
|
||||||
|
|
||||||
|
// total # of timesteps must be multiple of t_event
|
||||||
|
|
||||||
|
if (t_event == 0) error->universe_all("Invalid t_event in prd command");
|
||||||
|
if (nsteps % t_event)
|
||||||
|
error->universe_all("PRD nsteps must be multiple of t_event");
|
||||||
|
if (t_corr % t_event)
|
||||||
|
error->universe_all("PRD t_corr must be multiple of t_event");
|
||||||
|
|
||||||
|
// local storage
|
||||||
|
|
||||||
|
int me_universe = universe->me;
|
||||||
|
int nprocs_universe = universe->nprocs;
|
||||||
|
int nreplica = universe->nworlds;
|
||||||
|
int iworld = universe->iworld;
|
||||||
|
|
||||||
|
MPI_Comm_rank(world,&me);
|
||||||
|
MPI_Comm_size(world,&nprocs);
|
||||||
|
|
||||||
|
// comm_replica = communicator between same proc across replicas
|
||||||
|
// not used if replicas have unequal number of procs
|
||||||
|
// equal_size_replicas = 1 if all replicas have same # of procs
|
||||||
|
|
||||||
|
int color = me;
|
||||||
|
MPI_Comm_split(universe->uworld,color,0,&comm_replica);
|
||||||
|
|
||||||
|
flag = 0;
|
||||||
|
if (nreplica*nprocs == nprocs_universe) flag = 1;
|
||||||
|
MPI_Allreduce(&flag,&equal_size_replicas,1,MPI_INT,MPI_MIN,universe->uworld);
|
||||||
|
|
||||||
|
// workspace for inter-replica communication via gathers
|
||||||
|
|
||||||
|
natoms = static_cast<int> (atom->natoms);
|
||||||
|
|
||||||
|
displacements = NULL;
|
||||||
|
tagall = NULL;
|
||||||
|
xall = NULL;
|
||||||
|
imageall = NULL;
|
||||||
|
|
||||||
|
if (nreplica != nprocs_universe) {
|
||||||
|
displacements = new int[nprocs];
|
||||||
|
tagall = (int *) memory->smalloc(natoms*sizeof(int),"prd:tagall");
|
||||||
|
xall = memory->create_2d_double_array(natoms,3,"prd:xall");
|
||||||
|
imageall = (int *) memory->smalloc(natoms*sizeof(int),"prd:imageall");
|
||||||
|
}
|
||||||
|
|
||||||
|
// random_select = same RNG for each replica for multiple event selection
|
||||||
|
// random_dephase = unique RNG for each replica for dephasing
|
||||||
|
|
||||||
|
random_select = new RanPark(lmp,seed);
|
||||||
|
random_dephase = new RanMars(lmp,seed+iworld);
|
||||||
|
|
||||||
|
// create ComputeTemp class to monitor temperature
|
||||||
|
|
||||||
|
char **args = new char*[3];
|
||||||
|
args[0] = (char *) "prd_temp";
|
||||||
|
args[1] = (char *) "all";
|
||||||
|
args[2] = (char *) "temp";
|
||||||
|
modify->add_compute(3,args);
|
||||||
|
temperature = modify->compute[modify->ncompute-1];
|
||||||
|
|
||||||
|
// create Velocity class for velocity creation in dephasing
|
||||||
|
// pass it temperature compute, loop_setting, dist_setting settings
|
||||||
|
|
||||||
|
atom->check_mass();
|
||||||
|
velocity = new Velocity(lmp);
|
||||||
|
velocity->init_external("all");
|
||||||
|
|
||||||
|
args[0] = (char *) "temp";
|
||||||
|
args[1] = (char *) "prd_temp";
|
||||||
|
velocity->options(2,args);
|
||||||
|
args[0] = (char *) "loop";
|
||||||
|
args[1] = (char *) loop_setting;
|
||||||
|
if (loop_setting) velocity->options(2,args);
|
||||||
|
args[0] = (char *) "dist";
|
||||||
|
args[1] = (char *) dist_setting;
|
||||||
|
if (dist_setting) velocity->options(2,args);
|
||||||
|
|
||||||
|
// create FixEvent class to store event and pre-quench states
|
||||||
|
|
||||||
|
args[0] = (char *) "prd_event";
|
||||||
|
args[1] = (char *) "all";
|
||||||
|
args[2] = (char *) "EVENT";
|
||||||
|
modify->add_fix(3,args);
|
||||||
|
fix_event = (FixEvent *) modify->fix[modify->nfix-1];
|
||||||
|
|
||||||
|
// create Finish for timing output
|
||||||
|
|
||||||
|
finish = new Finish(lmp);
|
||||||
|
|
||||||
|
// string clean-up
|
||||||
|
|
||||||
|
delete [] args;
|
||||||
|
delete [] loop_setting;
|
||||||
|
delete [] dist_setting;
|
||||||
|
|
||||||
|
// assign FixEvent to event-detection compute
|
||||||
|
// necessary so it will know atom coords at last event
|
||||||
|
|
||||||
|
int icompute = modify->find_compute(id_compute);
|
||||||
|
if (icompute < 0) error->all("Could not find compute ID for PRD");
|
||||||
|
compute_event = modify->compute[icompute];
|
||||||
|
compute_event->reset_extra_compute_fix("prd_event");
|
||||||
|
|
||||||
|
// reset reneighboring criteria since will perform minimizations
|
||||||
|
|
||||||
|
neigh_every = neighbor->every;
|
||||||
|
neigh_delay = neighbor->delay;
|
||||||
|
neigh_dist_check = neighbor->dist_check;
|
||||||
|
|
||||||
|
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
|
||||||
|
if (me == 0)
|
||||||
|
error->warning("Resetting reneighboring criteria during PRD");
|
||||||
|
}
|
||||||
|
|
||||||
|
neighbor->every = 1;
|
||||||
|
neighbor->delay = 0;
|
||||||
|
neighbor->dist_check = 1;
|
||||||
|
|
||||||
|
// initialize PRD as if one long dynamics run
|
||||||
|
|
||||||
|
update->whichflag = 1;
|
||||||
|
update->nsteps = nsteps;
|
||||||
|
update->beginstep = update->firststep = update->ntimestep;
|
||||||
|
update->endstep = update->laststep = update->firststep + nsteps;
|
||||||
|
update->restrict_output = 1;
|
||||||
|
|
||||||
|
lmp->init();
|
||||||
|
|
||||||
|
// init minimizer settings and minimizer itself
|
||||||
|
|
||||||
|
update->etol = etol;
|
||||||
|
update->ftol = ftol;
|
||||||
|
update->max_eval = maxeval;
|
||||||
|
|
||||||
|
update->minimize->init();
|
||||||
|
|
||||||
|
// cannot use PRD with time-dependent fixes
|
||||||
|
|
||||||
|
for (int i = 0; i < modify->nfix; i++)
|
||||||
|
if (modify->fix[i]->time_depend)
|
||||||
|
error->all("Cannot use PRD with a time-dependent fix defined");
|
||||||
|
|
||||||
|
// perform PRD simulation
|
||||||
|
|
||||||
|
if (me_universe == 0 && universe->uscreen)
|
||||||
|
fprintf(universe->uscreen,"Setting up PRD ...\n");
|
||||||
|
|
||||||
|
if (me_universe == 0) {
|
||||||
|
if (universe->uscreen)
|
||||||
|
fprintf(universe->uscreen,"Step Clock Event Correlated Replica\n");
|
||||||
|
if (universe->ulogfile)
|
||||||
|
fprintf(universe->ulogfile,"Step Clock Event Correlated Replica\n");
|
||||||
|
}
|
||||||
|
|
||||||
|
// store hot state and quenched event for replica 0
|
||||||
|
// use share_event() to copy that info to all replicas
|
||||||
|
// this insures all start from same place
|
||||||
|
|
||||||
|
// need this line if quench() does only setup_minimal()
|
||||||
|
//update->minimize->setup();
|
||||||
|
|
||||||
|
fix_event->store_state();
|
||||||
|
quench();
|
||||||
|
share_event(0,0);
|
||||||
|
log_event();
|
||||||
|
|
||||||
|
// do full init/setup since are starting all replicas after event
|
||||||
|
// replica 0 bcasts temp to all replicas if user did not set temp_dephase
|
||||||
|
update->whichflag = 1;
|
||||||
|
lmp->init();
|
||||||
|
update->integrate->setup();
|
||||||
|
if (temp_flag == 0) {
|
||||||
|
if (universe->iworld == 0)
|
||||||
|
temp_dephase = temperature->compute_scalar();
|
||||||
|
MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[0],
|
||||||
|
universe->uworld);
|
||||||
|
}
|
||||||
|
|
||||||
|
// main loop: look for events until out of time
|
||||||
|
// (1) dephase independently on each proc after event
|
||||||
|
// (2) loop: dynamics, store state, quench, check event, restore state
|
||||||
|
// (3) share and record event
|
||||||
|
|
||||||
|
nbuild = ndanger = 0;
|
||||||
|
time_dephase = time_dynamics = time_quench = time_comm = time_output = 0.0;
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
double time_start = timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
while (update->ntimestep < update->endstep) {
|
||||||
|
dephase();
|
||||||
|
|
||||||
|
ireplica = -1;
|
||||||
|
while (update->ntimestep < update->endstep) {
|
||||||
|
dynamics();
|
||||||
|
fix_event->store_state();
|
||||||
|
quench();
|
||||||
|
ireplica = check_event();
|
||||||
|
if (ireplica >= 0) break;
|
||||||
|
fix_event->restore_state();
|
||||||
|
}
|
||||||
|
if (ireplica < 0) break;
|
||||||
|
|
||||||
|
// Can potentially be more efficient for correlated events by not
|
||||||
|
// sharing until correlated check has completed (this will complicate
|
||||||
|
// the dump (always on replica 0)).
|
||||||
|
share_event(ireplica,1);
|
||||||
|
log_event();
|
||||||
|
|
||||||
|
int restart_flag = 0;
|
||||||
|
if (output->restart_every && universe->iworld == 0)
|
||||||
|
if (fix_event->event_number % output->restart_every == 0)
|
||||||
|
restart_flag = 1;
|
||||||
|
|
||||||
|
// Correlated event loop
|
||||||
|
// -- We could have other procs doing dephasing during this time
|
||||||
|
int corr_endstep = update->ntimestep + t_corr;
|
||||||
|
while (update->ntimestep < corr_endstep) {
|
||||||
|
if (update->ntimestep == update->endstep) {
|
||||||
|
restart_flag = 0;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
dynamics();
|
||||||
|
fix_event->store_state();
|
||||||
|
quench();
|
||||||
|
int corr_event_check = check_event(ireplica);
|
||||||
|
if (corr_event_check >= 0) {
|
||||||
|
share_event(ireplica,2);
|
||||||
|
log_event();
|
||||||
|
corr_endstep = update->ntimestep + t_corr;
|
||||||
|
} else
|
||||||
|
fix_event->restore_state();
|
||||||
|
}
|
||||||
|
|
||||||
|
// do full init/setup since are starting all replicas after event
|
||||||
|
// -- also, need to get reneighbor before restart
|
||||||
|
update->whichflag = 1;
|
||||||
|
lmp->init();
|
||||||
|
update->integrate->setup();
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
if (t_corr > 0) replicate(ireplica);
|
||||||
|
|
||||||
|
// event replica bcasts temp to all replicas if user did not set temp_dephase
|
||||||
|
if (temp_flag == 0) {
|
||||||
|
if (ireplica == universe->iworld)
|
||||||
|
temp_dephase = temperature->compute_scalar();
|
||||||
|
MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[ireplica],
|
||||||
|
universe->uworld);
|
||||||
|
}
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_comm += timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
// write restart file of hot coords
|
||||||
|
if (restart_flag) {
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
output->write_restart(update->ntimestep);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_output += timer->array[TIME_LOOP];
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
// set total timers and counters so Finish() will process them
|
||||||
|
|
||||||
|
timer->array[TIME_LOOP] = time_start;
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
|
||||||
|
timer->array[TIME_PAIR] = time_dephase;
|
||||||
|
timer->array[TIME_BOND] = time_dynamics;
|
||||||
|
timer->array[TIME_KSPACE] = time_quench;
|
||||||
|
timer->array[TIME_COMM] = time_comm;
|
||||||
|
timer->array[TIME_OUTPUT] = time_output;
|
||||||
|
|
||||||
|
neighbor->ncalls = nbuild;
|
||||||
|
neighbor->ndanger = ndanger;
|
||||||
|
|
||||||
|
if (me_universe == 0) {
|
||||||
|
if (universe->uscreen)
|
||||||
|
fprintf(universe->uscreen,
|
||||||
|
"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
|
||||||
|
timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms);
|
||||||
|
if (universe->ulogfile)
|
||||||
|
fprintf(universe->ulogfile,
|
||||||
|
"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
|
||||||
|
timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms);
|
||||||
|
}
|
||||||
|
|
||||||
|
finish->end(2);
|
||||||
|
|
||||||
|
update->whichflag = 0;
|
||||||
|
update->firststep = update->laststep = 0;
|
||||||
|
update->beginstep = update->endstep = 0;
|
||||||
|
update->restrict_output = 0;
|
||||||
|
|
||||||
|
// reset reneighboring criteria
|
||||||
|
|
||||||
|
neighbor->every = neigh_every;
|
||||||
|
neighbor->delay = neigh_delay;
|
||||||
|
neighbor->dist_check = neigh_dist_check;
|
||||||
|
|
||||||
|
// clean up
|
||||||
|
|
||||||
|
delete [] displacements;
|
||||||
|
memory->sfree(tagall);
|
||||||
|
memory->destroy_2d_double_array(xall);
|
||||||
|
memory->sfree(imageall);
|
||||||
|
|
||||||
|
MPI_Comm_free(&comm_replica);
|
||||||
|
delete random_select;
|
||||||
|
delete random_dephase;
|
||||||
|
delete velocity;
|
||||||
|
delete finish;
|
||||||
|
modify->delete_compute("prd_temp");
|
||||||
|
modify->delete_fix("prd_event");
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
dephasing = one or more short runs with new random velocities
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::dephase()
|
||||||
|
{
|
||||||
|
int ntimestep_hold = update->ntimestep;
|
||||||
|
|
||||||
|
update->whichflag = 1;
|
||||||
|
update->nsteps = n_dephase*t_dephase;
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
|
||||||
|
for (int i = 0; i < n_dephase; i++) {
|
||||||
|
int seed = static_cast<int> (random_dephase->uniform() * MAXINT);
|
||||||
|
if (seed == 0) seed = 1;
|
||||||
|
velocity->create(temp_dephase,seed);
|
||||||
|
update->integrate->run(t_dephase);
|
||||||
|
if (temp_flag == 0) temp_dephase = temperature->compute_scalar();
|
||||||
|
}
|
||||||
|
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_dephase += timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
update->integrate->cleanup();
|
||||||
|
finish->end(0);
|
||||||
|
|
||||||
|
// reset timestep as if dephase did not occur
|
||||||
|
// clear timestep storage from computes, since now invalid
|
||||||
|
|
||||||
|
update->ntimestep = ntimestep_hold;
|
||||||
|
for (int i = 0; i < modify->ncompute; i++)
|
||||||
|
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
single short dynamics run
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::dynamics()
|
||||||
|
{
|
||||||
|
update->whichflag = 1;
|
||||||
|
update->nsteps = t_event;
|
||||||
|
|
||||||
|
lmp->init();
|
||||||
|
update->integrate->setup();
|
||||||
|
//modify->addstep_compute_all(update->ntimestep);
|
||||||
|
int ncalls = neighbor->ncalls;
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
update->integrate->run(t_event);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_dynamics += timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
nbuild += neighbor->ncalls - ncalls;
|
||||||
|
ndanger += neighbor->ndanger;
|
||||||
|
|
||||||
|
update->integrate->cleanup();
|
||||||
|
finish->end(0);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
quench minimization
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::quench()
|
||||||
|
{
|
||||||
|
int ntimestep_hold = update->ntimestep;
|
||||||
|
|
||||||
|
// need to change whichflag so that minimize->setup() calling
|
||||||
|
// modify->setup() will call fix->min_setup()
|
||||||
|
|
||||||
|
update->whichflag = 2;
|
||||||
|
update->nsteps = maxiter;
|
||||||
|
|
||||||
|
// these work
|
||||||
|
lmp->init();
|
||||||
|
update->minimize->setup();
|
||||||
|
|
||||||
|
// these do not work
|
||||||
|
//modify->addstep_compute_all(update->ntimestep);
|
||||||
|
//update->minimize->setup_minimal(1);
|
||||||
|
|
||||||
|
int ncalls = neighbor->ncalls;
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
update->minimize->run(maxiter);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_quench += timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
if (neighbor->ncalls == ncalls) quench_reneighbor = 0;
|
||||||
|
else quench_reneighbor = 1;
|
||||||
|
|
||||||
|
update->minimize->cleanup();
|
||||||
|
finish->end(0);
|
||||||
|
|
||||||
|
// reset timestep as if dephase did not occur
|
||||||
|
// clear timestep storage from computes, since now invalid
|
||||||
|
|
||||||
|
update->ntimestep = ntimestep_hold;
|
||||||
|
for (int i = 0; i < modify->ncompute; i++)
|
||||||
|
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
check for an event in any replica
|
||||||
|
if replica_num is non-negative only check for event on replica_num
|
||||||
|
if multiple events, choose one at random
|
||||||
|
return -1 if no event
|
||||||
|
else return ireplica = world in which event occured
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int PRD::check_event(int replica_num)
|
||||||
|
{
|
||||||
|
int worldflag,universeflag,scanflag,replicaflag,ireplica;
|
||||||
|
|
||||||
|
worldflag = 0;
|
||||||
|
if (compute_event->compute_scalar() > 0.0) worldflag = 1;
|
||||||
|
if (replica_num >= 0 && replica_num != universe->iworld) worldflag = 0;
|
||||||
|
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
if (me == 0) MPI_Allreduce(&worldflag,&universeflag,1,
|
||||||
|
MPI_INT,MPI_SUM,comm_replica);
|
||||||
|
MPI_Bcast(&universeflag,1,MPI_INT,0,world);
|
||||||
|
if (!universeflag)
|
||||||
|
ireplica = -1;
|
||||||
|
else {
|
||||||
|
if (universeflag > 1) {
|
||||||
|
int iwhich = static_cast<int> (universeflag*random_select->uniform()) + 1;
|
||||||
|
if (me == 0) MPI_Scan(&worldflag,&scanflag,1,
|
||||||
|
MPI_INT,MPI_SUM,comm_replica);
|
||||||
|
MPI_Bcast(&scanflag,1,MPI_INT,0,world);
|
||||||
|
if (scanflag != iwhich) worldflag = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (worldflag) replicaflag = universe->iworld;
|
||||||
|
else replicaflag = 0;
|
||||||
|
if (me == 0) MPI_Allreduce(&replicaflag,&ireplica,1,
|
||||||
|
MPI_INT,MPI_SUM,comm_replica);
|
||||||
|
MPI_Bcast(&ireplica,1,MPI_INT,0,world);
|
||||||
|
}
|
||||||
|
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_comm += timer->array[TIME_LOOP];
|
||||||
|
return ireplica;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
share quenched and hot coords owned by ireplica with all replicas
|
||||||
|
all replicas store event in fix_event
|
||||||
|
replica 0 dumps event snapshot
|
||||||
|
flag = 0 = called before PRD run
|
||||||
|
flag = 1 = called during PRD run = not correlated event
|
||||||
|
flag = 2 = called during PRD run = correlated event
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::share_event(int ireplica, int flag)
|
||||||
|
{
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
|
||||||
|
// communicate quenched coords to all replicas and store as event
|
||||||
|
// decrement event counter if flag = 0 since not really an event
|
||||||
|
|
||||||
|
replicate(ireplica);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_comm += timer->array[TIME_LOOP];
|
||||||
|
|
||||||
|
// Adjust time for last correlated event check (Not on first event)
|
||||||
|
int corr_adjust = t_corr;
|
||||||
|
if (fix_event->event_number<1 || flag == 2) corr_adjust = 0;
|
||||||
|
// Time since last correlated event check
|
||||||
|
int delta = update->ntimestep - fix_event->event_timestep - corr_adjust;
|
||||||
|
// If this is a correlated event, time was only on one partition
|
||||||
|
if (flag != 2) delta *= universe->nworlds;
|
||||||
|
delta += corr_adjust;
|
||||||
|
// Don't change the clock or timestep if this is a restart
|
||||||
|
if (flag == 0 && fix_event->event_number != 0)
|
||||||
|
fix_event->store_event(fix_event->event_timestep,0);
|
||||||
|
else {
|
||||||
|
fix_event->store_event(update->ntimestep,delta);
|
||||||
|
fix_event->replica_number = ireplica;
|
||||||
|
fix_event->correlated_event = 0;
|
||||||
|
if (flag == 2) fix_event->correlated_event = 1;
|
||||||
|
}
|
||||||
|
if (flag == 0) fix_event->event_number--;
|
||||||
|
|
||||||
|
// dump snapshot of quenched coords
|
||||||
|
// must reneighbor and compute forces before dumping
|
||||||
|
// since replica 0 possibly has new state from another replica
|
||||||
|
// addstep_compute_all insures eng/virial are calculated if needed
|
||||||
|
|
||||||
|
if (output->ndump && universe->iworld == 0) {
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
modify->addstep_compute_all(update->ntimestep);
|
||||||
|
update->integrate->setup_minimal(1);
|
||||||
|
output->write_dump(update->ntimestep);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_output += timer->array[TIME_LOOP];
|
||||||
|
}
|
||||||
|
|
||||||
|
// restore and communicate hot coords to all replicas
|
||||||
|
|
||||||
|
fix_event->restore_state();
|
||||||
|
timer->barrier_start(TIME_LOOP);
|
||||||
|
replicate(ireplica);
|
||||||
|
timer->barrier_stop(TIME_LOOP);
|
||||||
|
time_comm += timer->array[TIME_LOOP];
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
universe proc 0 prints event info
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::log_event()
|
||||||
|
{
|
||||||
|
if (universe->me == 0) {
|
||||||
|
if (universe->uscreen)
|
||||||
|
fprintf(universe->uscreen,"%d %d %d %d %d\n",
|
||||||
|
fix_event->event_timestep,fix_event->clock,
|
||||||
|
fix_event->event_number,fix_event->correlated_event,
|
||||||
|
fix_event->replica_number);
|
||||||
|
if (universe->ulogfile)
|
||||||
|
fprintf(universe->ulogfile,"%d %d %d %d %d\n",
|
||||||
|
fix_event->event_timestep,fix_event->clock,
|
||||||
|
fix_event->event_number,fix_event->correlated_event,
|
||||||
|
fix_event->replica_number);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
communicate atom coords and image flags in ireplica to all other replicas
|
||||||
|
one proc per replica:
|
||||||
|
direct overwrite via bcast
|
||||||
|
equal procs per replica and no replica reneighbored:
|
||||||
|
direct overwrite via bcast
|
||||||
|
unequal procs per replica or reneighboring occurred:
|
||||||
|
collect to root proc of event replica
|
||||||
|
bcast to roots of other replicas
|
||||||
|
bcast within each replica
|
||||||
|
each proc extracts info for atoms it owns using atom IDs
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::replicate(int ireplica)
|
||||||
|
{
|
||||||
|
int nreplica = universe->nworlds;
|
||||||
|
int nprocs_universe = universe->nprocs;
|
||||||
|
int i,m,flag,commflag;
|
||||||
|
int counts[nprocs];
|
||||||
|
|
||||||
|
if (nreplica == nprocs_universe) commflag = 0;
|
||||||
|
else if (equal_size_replicas) {
|
||||||
|
flag = 0;
|
||||||
|
if (quench_reneighbor) flag = 1;
|
||||||
|
MPI_Allreduce(&flag,&commflag,1,MPI_INT,MPI_MAX,universe->uworld);
|
||||||
|
} else commflag = 1;
|
||||||
|
|
||||||
|
if (commflag == 0) {
|
||||||
|
MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica);
|
||||||
|
MPI_Bcast(atom->x[0],3*atom->nlocal,MPI_DOUBLE,ireplica,comm_replica);
|
||||||
|
} else {
|
||||||
|
if (universe->iworld == ireplica) {
|
||||||
|
MPI_Gather(&atom->nlocal,1,MPI_INT,counts,1,MPI_INT,0,world);
|
||||||
|
displacements[0] = 0;
|
||||||
|
for (i = 0; i < nprocs-1; i++)
|
||||||
|
displacements[i+1] = displacements[i] + counts[i];
|
||||||
|
MPI_Gatherv(atom->tag,atom->nlocal,MPI_INT,
|
||||||
|
tagall,counts,displacements,MPI_INT,0,world);
|
||||||
|
MPI_Gatherv(atom->image,atom->nlocal,MPI_INT,
|
||||||
|
imageall,counts,displacements,MPI_INT,0,world);
|
||||||
|
for (i = 0; i < nprocs; i++) counts[i] *= 3;
|
||||||
|
for (i = 0; i < nprocs-1; i++)
|
||||||
|
displacements[i+1] = displacements[i] + counts[i];
|
||||||
|
MPI_Gatherv(atom->x[0],3*atom->nlocal,MPI_DOUBLE,
|
||||||
|
xall[0],counts,displacements,MPI_DOUBLE,0,world);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (me == 0) {
|
||||||
|
MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
|
||||||
|
MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
|
||||||
|
MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
|
||||||
|
}
|
||||||
|
|
||||||
|
MPI_Bcast(tagall,natoms,MPI_INT,0,world);
|
||||||
|
MPI_Bcast(imageall,natoms,MPI_INT,0,world);
|
||||||
|
MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world);
|
||||||
|
|
||||||
|
double **x = atom->x;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
for (i = 0; i < natoms; i++) {
|
||||||
|
m = atom->map(tagall[i]);
|
||||||
|
if (m >= 0 && m < nlocal) {
|
||||||
|
x[m][0] = xall[i][0];
|
||||||
|
x[m][1] = xall[i][1];
|
||||||
|
x[m][2] = xall[i][2];
|
||||||
|
atom->image[m] = imageall[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
parse optional parameters at end of PRD input line
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PRD::options(int narg, char **arg)
|
||||||
|
{
|
||||||
|
if (narg < 0) error->all("Illegal prd command");
|
||||||
|
|
||||||
|
// set defaults
|
||||||
|
|
||||||
|
etol = 0.1;
|
||||||
|
ftol = 0.1;
|
||||||
|
maxiter = 40;
|
||||||
|
maxeval = 50;
|
||||||
|
temp_flag = 0;
|
||||||
|
loop_setting = NULL;
|
||||||
|
dist_setting = NULL;
|
||||||
|
|
||||||
|
int iarg = 0;
|
||||||
|
while (iarg < narg) {
|
||||||
|
if (strcmp(arg[iarg],"min") == 0) {
|
||||||
|
if (iarg+5 > narg) error->all("Illegal prd command");
|
||||||
|
etol = atof(arg[iarg+1]);
|
||||||
|
ftol = atof(arg[iarg+2]);
|
||||||
|
maxiter = atoi(arg[iarg+3]);
|
||||||
|
maxeval = atoi(arg[iarg+4]);
|
||||||
|
if (maxiter < 0) error->all("Illegal prd command");
|
||||||
|
iarg += 5;
|
||||||
|
|
||||||
|
} else if (strcmp(arg[iarg],"temp") == 0) {
|
||||||
|
if (iarg+2 > narg) error->all("Illegal prd command");
|
||||||
|
temp_flag = 1;
|
||||||
|
temp_dephase = atof(arg[iarg+1]);
|
||||||
|
if (temp_dephase <= 0.0) error->all("Illegal prd command");
|
||||||
|
iarg += 2;
|
||||||
|
|
||||||
|
} else if (strcmp(arg[iarg],"vel") == 0) {
|
||||||
|
if (iarg+3 > narg) error->all("Illegal prd command");
|
||||||
|
|
||||||
|
if (strcmp(arg[iarg+1],"all") == 0) loop_setting = NULL;
|
||||||
|
else if (strcmp(arg[iarg+1],"local") == 0) loop_setting = NULL;
|
||||||
|
else if (strcmp(arg[iarg+1],"geom") == 0) loop_setting = NULL;
|
||||||
|
else error->all("Illegal prd command");
|
||||||
|
int n = strlen(arg[iarg+1]) + 1;
|
||||||
|
loop_setting = new char[n];
|
||||||
|
strcpy(loop_setting,arg[iarg+1]);
|
||||||
|
|
||||||
|
if (strcmp(arg[iarg+2],"uniform") == 0) dist_setting = NULL;
|
||||||
|
else if (strcmp(arg[iarg+2],"gaussian") == 0) dist_setting = NULL;
|
||||||
|
else error->all("Illegal prd command");
|
||||||
|
n = strlen(arg[iarg+2]) + 1;
|
||||||
|
dist_setting = new char[n];
|
||||||
|
strcpy(dist_setting,arg[iarg+2]);
|
||||||
|
|
||||||
|
iarg += 3;
|
||||||
|
} else error->all("Illegal prd command");
|
||||||
|
}
|
||||||
|
|
||||||
|
// Set defaults
|
||||||
|
if (loop_setting == NULL) {
|
||||||
|
loop_setting = new char[5];
|
||||||
|
strcpy(loop_setting,"geom");
|
||||||
|
}
|
||||||
|
if (dist_setting == NULL) {
|
||||||
|
dist_setting = new char[9];
|
||||||
|
strcpy(dist_setting,"gaussian");
|
||||||
|
}
|
||||||
|
}
|
||||||
65
src/PRD/prd.h
Normal file
65
src/PRD/prd.h
Normal file
@ -0,0 +1,65 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifndef PRD_H
|
||||||
|
#define PRD_H
|
||||||
|
|
||||||
|
#include "pointers.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class PRD : protected Pointers {
|
||||||
|
public:
|
||||||
|
PRD(class LAMMPS *);
|
||||||
|
~PRD() {}
|
||||||
|
void command(int, char **);
|
||||||
|
|
||||||
|
private:
|
||||||
|
int me,nprocs;
|
||||||
|
int nsteps,t_event,n_dephase,t_dephase,t_corr;
|
||||||
|
double etol,ftol,temp_dephase;
|
||||||
|
int maxiter,maxeval,temp_flag;
|
||||||
|
char *loop_setting,*dist_setting;
|
||||||
|
|
||||||
|
int equal_size_replicas,natoms;
|
||||||
|
int neigh_every,neigh_delay,neigh_dist_check;
|
||||||
|
int nbuild,ndanger;
|
||||||
|
int quench_reneighbor;
|
||||||
|
|
||||||
|
double time_dephase,time_dynamics,time_quench,time_comm,time_output;
|
||||||
|
|
||||||
|
MPI_Comm comm_replica;
|
||||||
|
int *tagall,*displacements,*imageall;
|
||||||
|
double **xall;
|
||||||
|
|
||||||
|
class RanPark *random_select;
|
||||||
|
class RanMars *random_dephase;
|
||||||
|
class Compute *compute_event;
|
||||||
|
class FixEvent *fix_event;
|
||||||
|
class Velocity *velocity;
|
||||||
|
class Compute *temperature;
|
||||||
|
class Finish *finish;
|
||||||
|
|
||||||
|
void dephase();
|
||||||
|
void dynamics();
|
||||||
|
void quench();
|
||||||
|
int check_event(int replica = -1);
|
||||||
|
void share_event(int, int);
|
||||||
|
void log_event();
|
||||||
|
void replicate(int);
|
||||||
|
void options(int, char **);
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
36
src/PRD/style_prd.h
Normal file
36
src/PRD/style_prd.h
Normal file
@ -0,0 +1,36 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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 CommandInclude
|
||||||
|
#include "prd.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef CommandClass
|
||||||
|
CommandStyle(prd,PRD)
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef ComputeInclude
|
||||||
|
#include "compute_event_displace.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef ComputeClass
|
||||||
|
ComputeStyle(event/displace,ComputeEventDisplace)
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef FixInclude
|
||||||
|
#include "fix_event.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef FixClass
|
||||||
|
FixStyle(EVENT,FixEvent)
|
||||||
|
#endif
|
||||||
@ -385,6 +385,7 @@ RegionStyle(union,RegUnion)
|
|||||||
#include "style_opt.h"
|
#include "style_opt.h"
|
||||||
#include "style_peri.h"
|
#include "style_peri.h"
|
||||||
#include "style_poems.h"
|
#include "style_poems.h"
|
||||||
|
#include "style_prd.h"
|
||||||
#include "style_reax.h"
|
#include "style_reax.h"
|
||||||
#include "style_xtc.h"
|
#include "style_xtc.h"
|
||||||
|
|
||||||
|
|||||||
@ -1,20 +0,0 @@
|
|||||||
/* ----------------------------------------------------------------------
|
|
||||||
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 ComputeInclude
|
|
||||||
#include "compute_ackland_atom.h"
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef ComputeClass
|
|
||||||
ComputeStyle(ackland/atom,ComputeAcklandAtom)
|
|
||||||
#endif
|
|
||||||
|
|||||||
@ -1,30 +0,0 @@
|
|||||||
/* ----------------------------------------------------------------------
|
|
||||||
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 KSpaceInclude
|
|
||||||
#include "ewald_n.h"
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef KSpaceClass
|
|
||||||
KSpaceStyle(ewald/n,EwaldN)
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef PairInclude
|
|
||||||
#include "pair_buck_coul.h"
|
|
||||||
#include "pair_lj_coul.h"
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef PairClass
|
|
||||||
PairStyle(buck/coul,PairBuckCoul)
|
|
||||||
PairStyle(lj/coul,PairLJCoul)
|
|
||||||
#endif
|
|
||||||
|
|||||||
Reference in New Issue
Block a user