jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added tad example 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5475 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2011-01-05 00:23:34 +00:00
parent 4e99c5022e
commit c69f3aea57
16 changed files with 1597 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/REPLICA/Install.sh
View File

@ -4,32 +4,44 @@ if (test $1 = 1) then
cp compute_event_displace.cpp ..
cp fix_event.cpp ..
cp fix_event_prd.cpp ..
cp fix_event_tad.cpp ..
cp fix_neb.cpp ..
cp neb.cpp ..
cp prd.cpp ..
cp tad.cpp ..
cp temper.cpp ..
cp compute_event_displace.h ..
cp fix_event.h ..
cp fix_event_prd.h ..
cp fix_event_tad.h ..
cp fix_neb.h ..
cp neb.h ..
cp prd.h ..
cp tad.h ..
cp temper.h ..
elif (test $1 = 0) then
rm ../compute_event_displace.cpp
rm ../fix_event.cpp
rm ../fix_event_prd.cpp
rm ../fix_event_tad.cpp
rm ../fix_neb.cpp
rm ../neb.cpp
rm ../prd.cpp
rm ../tad.cpp
rm ../temper.cpp
rm ../compute_event_displace.h
rm ../fix_event.h
rm ../fix_event_prd.h
rm ../fix_event_tad.h
rm ../fix_neb.h
rm ../neb.h
rm ../prd.h
rm ../tad.h
rm ../temper.h
fi

11
src/REPLICA/compute_event_displace.cpp
View File

@ -23,7 +23,7 @@
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "fix_event.h"
#include "memory.h"
#include "error.h"
#include "update.h"
@ -70,9 +70,10 @@ void ComputeEventDisplace::init()
if (id_event != NULL) {
int ifix = modify->find_fix(id_event);
if (ifix < 0) error->all("Could not find compute event/displace fix ID");
fix = modify->fix[ifix];
fix_event = (FixEvent*) modify->fix[ifix];
if (strcmp(fix->style,"EVENT") != 0)
if (strcmp(fix_event->style,"EVENT/PRD") != 0 &&
strcmp(fix_event->style,"EVENT/TAD") != 0)
error->all("Compute event/displace has invalid fix event assigned");
}
@ -90,7 +91,7 @@ double ComputeEventDisplace::compute_scalar()
if (id_event == NULL) return 0.0;
double event = 0.0;
double **xevent = fix->array_atom;
double **xevent = fix_event->array_atom;
double **x = atom->x;
int *mask = atom->mask;
@ -119,7 +120,6 @@ double ComputeEventDisplace::compute_scalar()
break;
}
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@ -138,6 +138,7 @@ double ComputeEventDisplace::compute_scalar()
}
MPI_Allreduce(&event,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
return scalar;
}

2
src/REPLICA/compute_event_displace.h
View File

@ -36,7 +36,7 @@ class ComputeEventDisplace : public Compute {
int triclinic;
double displace_distsq;
char *id_event;
class Fix *fix;
class FixEvent *fix_event;
};
}

108
src/REPLICA/fix_event.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mike Brown (SNL)
Contributing author: Mike Brown (SNL), Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "stdlib.h"
@ -43,13 +43,10 @@ FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) :
xevent = NULL;
xold = NULL;
vold = NULL;
imageold = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
event_number = 0;
event_timestep = update->ntimestep;
clock = 0;
}
/* ---------------------------------------------------------------------- */
@ -64,6 +61,7 @@ FixEvent::~FixEvent()
memory->destroy_2d_double_array(xevent);
memory->destroy_2d_double_array(xold);
memory->destroy_2d_double_array(vold);
memory->sfree(imageold);
}
@ -76,12 +74,10 @@ int FixEvent::setmask()
/* ----------------------------------------------------------------------
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
called when an event occurs
------------------------------------------------------------------------- */
void FixEvent::store_event(int timestep, int delta_clock)
void FixEvent::store_event()
{
double **x = atom->x;
int *image = atom->image;
@ -90,9 +86,41 @@ void FixEvent::store_event(int timestep, int delta_clock)
for (int i = 0; i < nlocal; i++)
domain->unmap(x[i],image[i],xevent[i]);
event_timestep = timestep;
clock += delta_clock;
event_number++;
// printf("store_event %g %d %g %d \n",
// x[8][1],image[8],xevent[8][1],0);
}
/* ----------------------------------------------------------------------
restore atom coords to quenched initial state
called prior to NEB calculation
------------------------------------------------------------------------- */
void FixEvent::restore_event()
{
double **x = atom->x;
int *image = atom->image;
int nlocal = atom->nlocal;
// printf("restore_event1 %g %d %g %d \n",
// x[8][1],image[8],xevent[8][1],0);
for (int i = 0; i < nlocal; i++) {
x[i][0] = xevent[i][0];
x[i][1] = xevent[i][1];
x[i][2] = xevent[i][2];
// Since xevent is unwrapped coordinate, need to
// adjust image flags when remapping
image[i] = (512 << 20) | (512 << 10) | 512;
domain->remap(x[i],image[i]);
// domain->remap(x[i]);
}
// printf("restore_event2 %g %d %g %d \n",
// x[8][1],image[8],xevent[8][1],0);
}
/* ----------------------------------------------------------------------
@ -104,14 +132,20 @@ void FixEvent::store_event(int timestep, int delta_clock)
void FixEvent::store_state()
{
double **x = atom->x;
double **f = atom->f;
double **v = atom->v;
int *image = atom->image;
int nlocal = atom->nlocal;
// printf("store_state %g %d %g %d \n",
// xold[8][1],imageold[8],x[8][1],image[8]);
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];
vold[i][0] = v[i][0];
vold[i][1] = v[i][1];
vold[i][2] = v[i][2];
imageold[i] = image[i];
}
}
@ -124,13 +158,20 @@ void FixEvent::store_state()
void FixEvent::restore_state()
{
double **x = atom->x;
double **v = atom->v;
int *image = atom->image;
int nlocal = atom->nlocal;
// printf("restore_state %g %d %g %d \n",
// xold[8][1],imageold[8],x[8][1],image[8]);
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];
v[i][0] = vold[i][0];
v[i][1] = vold[i][1];
v[i][2] = vold[i][2];
image[i] = imageold[i];
}
}
@ -154,6 +195,7 @@ 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");
vold = memory->grow_2d_double_array(vold,nmax,3,"event:vold");
imageold = (int *)
memory->srealloc(imageold,nmax*sizeof(int),"event:imageold");
@ -174,6 +216,9 @@ void FixEvent::copy_arrays(int i, int j)
xold[j][0] = xold[i][0];
xold[j][1] = xold[i][1];
xold[j][2] = xold[i][2];
vold[j][0] = vold[i][0];
vold[j][1] = vold[i][1];
vold[j][2] = vold[i][2];
imageold[j] = imageold[i];
}
@ -189,9 +234,12 @@ int FixEvent::pack_exchange(int i, double *buf)
buf[3] = xold[i][0];
buf[4] = xold[i][1];
buf[5] = xold[i][2];
buf[6] = imageold[i];
buf[6] = vold[i][0];
buf[7] = vold[i][1];
buf[8] = vold[i][2];
buf[9] = imageold[i];
return 7;
return 10;
}
/* ----------------------------------------------------------------------
@ -206,9 +254,12 @@ int FixEvent::unpack_exchange(int nlocal, double *buf)
xold[nlocal][0] = buf[3];
xold[nlocal][1] = buf[4];
xold[nlocal][2] = buf[5];
imageold[nlocal] = static_cast<int>(buf[6]);
vold[nlocal][0] = buf[6];
vold[nlocal][1] = buf[7];
vold[nlocal][2] = buf[8];
imageold[nlocal] = static_cast<int>(buf[9]);
return 7;
return 10;
}
/* ----------------------------------------------------------------------
@ -217,20 +268,6 @@ int FixEvent::unpack_exchange(int nlocal, double *buf)
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;
list[n++] = ncoincident;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
@ -239,13 +276,4 @@ void FixEvent::write_restart(FILE *fp)
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++]);
ncoincident = static_cast<int> (list[n++]);
}

25
src/REPLICA/fix_event.h
View File

@ -11,12 +11,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(EVENT,FixEvent)
#else
#ifndef LMP_FIX_EVENT_H
#define LMP_FIX_EVENT_H
@ -26,15 +20,9 @@ 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
int ncoincident; // # of simultaneous events on different replicas
FixEvent(class LAMMPS *, int, char **);
~FixEvent();
virtual ~FixEvent()=0; // Use destructor to make base class virtual
int setmask();
double memory_usage();
@ -42,22 +30,23 @@ class FixEvent : public Fix {
void copy_arrays(int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
void write_restart(FILE *);
void restart(char *);
virtual void write_restart(FILE *);
virtual void restart(char *);
// methods specific to FixEvent, invoked by PRD
// methods specific to FixEvent
void store_event(int, int);
virtual void store_event(); // base class stores quenched atoms
void restore_event(); // restore quenched atoms
void store_state();
void restore_state();
private:
double **xevent; // atom coords at last event
double **xold; // atom coords for reset/restore
double **vold; // atom vels for reset/restore
int *imageold; // image flags for reset/restore
};
}
#endif
#endif

107
src/REPLICA/fix_event_prd.cpp Normal file
View File

@ -0,0 +1,107 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mike Brown (SNL)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "fix_event_prd.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;
/* ---------------------------------------------------------------------- */
FixEventPRD::FixEventPRD(LAMMPS *lmp, int narg, char **arg) :
FixEvent(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal fix event command");
restart_global = 1;
event_number = 0;
event_timestep = update->ntimestep;
clock = 0;
}
/* ---------------------------------------------------------------------- */
FixEventPRD::~FixEventPRD()
{
}
/* ----------------------------------------------------------------------
save current atom coords as an event (via call to base class)
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 FixEventPRD::store_event(int timestep, int delta_clock)
{
FixEvent::store_event();
event_timestep = timestep;
clock += delta_clock;
event_number++;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixEventPRD::write_restart(FILE *fp)
{
int n = 0;
double list[6];
list[n++] = event_number;
list[n++] = event_timestep;
list[n++] = clock;
list[n++] = replica_number;
list[n++] = correlated_event;
list[n++] = ncoincident;
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 FixEventPRD::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
event_number = static_cast<int> (list[n++]);
printf("Event number restart = %d\n",event_number);
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++]);
ncoincident = static_cast<int> (list[n++]);
}

53
src/REPLICA/fix_event_prd.h Normal file
View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(EVENT/PRD,FixEventPRD)
#else
#ifndef LMP_FIX_EVENT_PRD_H
#define LMP_FIX_EVENT_PRD_H
#include "fix_event.h"
namespace LAMMPS_NS {
class FixEventPRD : public FixEvent {
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
int ncoincident; // # of simultaneous events on different replicas
FixEventPRD(class LAMMPS *, int, char **);
~FixEventPRD();
void write_restart(FILE *);
void restart(char *);
// methods specific to FixEventPRD, invoked by PRD
void store_event(int, int);
private:
};
}
#endif
#endif

100
src/REPLICA/fix_event_tad.cpp Normal file
View File

@ -0,0 +1,100 @@
/* ----------------------------------------------------------------------
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_tad.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;
/* ---------------------------------------------------------------------- */
FixEventTAD::FixEventTAD(LAMMPS *lmp, int narg, char **arg) :
FixEvent(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal fix event command");
restart_global = 1;
event_number = 0;
event_timestep = update->ntimestep;
tlo = 0.0;
ebarrier = 0.0;
}
/* ---------------------------------------------------------------------- */
FixEventTAD::~FixEventTAD()
{
}
/* ----------------------------------------------------------------------
save current atom coords as an event (via call to base class)
called when an event occurs in some replica
set event_timestep = when event occurred
------------------------------------------------------------------------- */
void FixEventTAD::store_event(int timestep)
{
FixEvent::store_event();
event_timestep = timestep;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixEventTAD::write_restart(FILE *fp)
{
int n = 0;
double list[4];
list[n++] = event_number;
list[n++] = event_timestep;
list[n++] = tlo;
list[n++] = ebarrier;
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 FixEventTAD::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
event_number = static_cast<int> (list[n++]);
event_timestep = static_cast<int> (list[n++]);
tlo = static_cast<double> (list[n++]);
ebarrier = list[n++];
}

51
src/REPLICA/fix_event_tad.h Normal file
View File

@ -0,0 +1,51 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(EVENT/TAD,FixEventTAD)
#else
#ifndef LMP_FIX_EVENT_TAD_H
#define LMP_FIX_EVENT_TAD_H
#include "fix_event.h"
namespace LAMMPS_NS {
class FixEventTAD : public FixEvent {
public:
int event_number; // event counter
int event_timestep; // timestep of last event
double tlo; // event time at low temperature
double ebarrier; // energy barrier for this event
FixEventTAD(class LAMMPS *, int, char **);
~FixEventTAD();
void write_restart(FILE *);
void restart(char *);
// methods specific to FixEventTAD, invoked by TAD
void store_event(int);
private:
};
}
#endif
#endif

31
src/REPLICA/fix_neb.cpp
View File

@ -132,7 +132,6 @@ void FixNEB::min_setup(int vflag)
void FixNEB::min_post_force(int vflag)
{
MPI_Status status;
MPI_Request request;
double vprev,vnext,vmax,vmin;
double delx,dely,delz;
double delta1[3],delta2[3];
@ -142,21 +141,10 @@ void FixNEB::min_post_force(int vflag)
veng = pe->compute_scalar();
if (ireplica == 0)
MPI_Send(&veng,1,MPI_DOUBLE,procnext,0,uworld);
else if (ireplica == nreplica-1) {
MPI_Irecv(&vprev,1,MPI_DOUBLE,procprev,0,uworld,&request);
MPI_Wait(&request,&status);
} else
if (ireplica < nreplica-1)
MPI_Sendrecv(&veng,1,MPI_DOUBLE,procnext,0,
&vprev,1,MPI_DOUBLE,procprev,0,uworld,&status);
if (ireplica == 0) {
MPI_Irecv(&vnext,1,MPI_DOUBLE,procnext,0,uworld,&request);
MPI_Wait(&request,&status);
} else if (ireplica == nreplica-1)
MPI_Send(&veng,1,MPI_DOUBLE,procprev,0,uworld);
else
if (ireplica > 0)
MPI_Sendrecv(&veng,1,MPI_DOUBLE,procprev,0,
&vnext,1,MPI_DOUBLE,procnext,0,uworld,&status);
@ -169,21 +157,10 @@ void FixNEB::min_post_force(int vflag)
int nlocal = atom->nlocal;
if (nlocal != nebatoms) error->one("Atom count changed in fix neb");
if (ireplica == 0)
MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);
else if (ireplica == nreplica-1) {
MPI_Irecv(xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request);
MPI_Wait(&request,&status);
} else
if (ireplica < nreplica-1)
MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procnext,0,
xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&status);
if (ireplica == 0) {
MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
MPI_Wait(&request,&status);
} else if (ireplica == nreplica-1)
MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
else
if (ireplica > 0)
MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procprev,0,
xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&status);

124
src/REPLICA/neb.cpp
View File

@ -40,6 +40,50 @@ using namespace LAMMPS_NS;
NEB::NEB(LAMMPS *lmp) : Pointers(lmp) {}
/* ----------------------------------------------------------------------
internal NEB constructor
------------------------------------------------------------------------- */
NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
: Pointers(lmp)
{
double delx,dely,delz;
etol = etol_in;
ftol = ftol_in;
n1steps = n1steps_in;
n2steps = n2steps_in;
nevery = nevery_in;
// replica info
nreplica = universe->nworlds;
ireplica = universe->iworld;
me_universe = universe->me;
uworld = universe->uworld;
MPI_Comm_rank(world,&me);
// generate linear interpolate replica
double fraction = ireplica/(nreplica-1.0);
double **x = atom->x;
int nlocal = atom->nlocal;
int ii = 0;
for (int i = 0; i < nlocal; i++) {
delx = buf_final[ii] - buf_init[ii];
dely = buf_final[ii+1] - buf_init[ii+1];
delz = buf_final[ii+2] - buf_init[ii+2];
domain->minimum_image(delx,dely,delz);
x[i][0] = buf_init[ii] + fraction*delx;
x[i][1] = buf_init[ii+1] + fraction*dely;
x[i][2] = buf_init[ii+2] + fraction*delz;
ii += 3;
}
}
/* ---------------------------------------------------------------------- */
NEB::~NEB()
@ -58,16 +102,18 @@ void NEB::command(int narg, char **arg)
if (domain->box_exist == 0)
error->all("NEB command before simulation box is defined");
if (narg != 5) error->universe_all("Illegal NEB command");
if (narg != 6) error->universe_all("Illegal NEB command");
double ftol = atof(arg[0]);
int n1steps = atoi(arg[1]);
int n2steps = atoi(arg[2]);
int nevery = atoi(arg[3]);
char *infile = arg[4];
etol = atof(arg[0]);
ftol = atof(arg[1]);
n1steps = atoi(arg[2]);
n2steps = atoi(arg[3]);
nevery = atoi(arg[4]);
infile = arg[5];
// error checks
if (etol < 0.0) error->all("Illegal NEB command");
if (ftol < 0.0) error->all("Illegal NEB command");
if (nevery == 0) error->universe_all("Illegal NEB command");
if (n1steps % nevery || n2steps % nevery)
@ -81,6 +127,22 @@ void NEB::command(int narg, char **arg)
uworld = universe->uworld;
MPI_Comm_rank(world,&me);
// read in file of final state atom coords and reset my coords
readfile(infile);
// run the NEB calculation
run();
}
/* ----------------------------------------------------------------------
run NEB on multiple replicas
------------------------------------------------------------------------- */
void NEB::run()
{
// create MPI communicator for root proc from each world
int color;
@ -104,13 +166,14 @@ void NEB::command(int narg, char **arg)
if (ineb == modify->nfix) error->all("NEB requires use of fix neb");
fneb = (FixNEB *) modify->fix[ineb];
all = memory->create_2d_double_array(nreplica,3,"neb:all");
nall = 4;
all = memory->create_2d_double_array(nreplica,nall,"neb:all");
rdist = new double[nreplica];
// initialize LAMMPS
update->whichflag = 2;
update->etol = 0.0;
update->etol = etol;
update->ftol = ftol;
update->multireplica = 1;
@ -119,11 +182,7 @@ void NEB::command(int narg, char **arg)
if (update->minimize->searchflag)
error->all("NEB requires damped dynamics minimizer");
// read in file of final state atom coords and reset my coords
readfile(infile);
// setup regular NEB minimizaiton
// setup regular NEB minimization
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up regular NEB ...\n");
@ -139,10 +198,12 @@ void NEB::command(int narg, char **arg)
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
@ -206,10 +267,12 @@ void NEB::command(int narg, char **arg)
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
@ -220,7 +283,6 @@ void NEB::command(int narg, char **arg)
// damped dynamic min styles insure all replicas converge together
timer->barrier_start(TIME_LOOP);
while (update->minimize->niter < n2steps) {
update->minimize->run(nevery);
print_status();
@ -372,38 +434,39 @@ void NEB::print_status()
double fmaxatom;
MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots);
double one[3];
double one[4];
one[0] = fneb->veng;
one[1] = fneb->plen;
one[2] = fneb->nlen;
one[nall-1] = fneb->gradvnorm;
if (output->thermo->normflag) one[0] /= atom->natoms;
if (me == 0) MPI_Allgather(one,3,MPI_DOUBLE,&all[0][0],3,MPI_DOUBLE,roots);
if (me == 0)
MPI_Allgather(one,nall,MPI_DOUBLE,&all[0][0],nall,MPI_DOUBLE,roots);
rdist[0] = 0.0;
for (int i = 1; i < nreplica; i++)
rdist[i] = rdist[i-1] + all[i][1];
double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2];
if (endpt > 0.0)
for (int i = 1; i < nreplica; i++)
rdist[i] /= endpt;
for (int i = 1; i < nreplica; i++)
rdist[i] /= endpt;
// look up GradV for the initial, final, and climbing replicas
// these should be identical to fnorm2
// but to be safe take them straight from fix neb
// these are identical to fnorm2, but to be safe we
// take them straight from fix_neb
double gradvnorm0, gradvnorm1, gradvnormc;
int irep;
irep = 0;
if (me_universe == irep) gradvnorm0 = fneb->gradvnorm;
MPI_Bcast(&gradvnorm0,1,MPI_DOUBLE,irep,uworld);
gradvnorm0 = all[irep][3];
irep = nreplica-1;
if (me_universe == irep) gradvnorm1 = fneb->gradvnorm;
MPI_Bcast(&gradvnorm1,1,MPI_DOUBLE,irep,uworld);
gradvnorm1 = all[irep][3];
irep = fneb->rclimber;
if (irep > -1) {
if (me_universe == irep) gradvnormc = fneb->gradvnorm;
MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld);
gradvnormc = all[irep][3];
ebf = all[irep][0]-all[0][0];
ebr = all[irep][0]-all[nreplica-1][0];
} else {
double vmax = all[0][0];
int top = 0;
@ -413,8 +476,9 @@ void NEB::print_status()
top = m;
}
irep = top;
if (me_universe == irep) gradvnormc = fneb->gradvnorm;
MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld);
gradvnormc = all[irep][3];
ebf = all[irep][0]-all[0][0];
ebr = all[irep][0]-all[nreplica-1][0];
}
if (me_universe == 0) {
@ -423,6 +487,7 @@ void NEB::print_status()
fmaxreplica,fmaxatom);
fprintf(universe->uscreen,"%g %g %g ",
gradvnorm0,gradvnorm1,gradvnormc);
fprintf(universe->uscreen,"%g %g %g ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++)
fprintf(universe->uscreen,"%g %g ",rdist[i],all[i][0]);
fprintf(universe->uscreen,"\n");
@ -432,6 +497,7 @@ void NEB::print_status()
fmaxreplica,fmaxatom);
fprintf(universe->ulogfile,"%g %g %g ",
gradvnorm0,gradvnorm1,gradvnormc);
fprintf(universe->ulogfile,"%g %g %g ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++)
fprintf(universe->ulogfile,"%g %g ",rdist[i],all[i][0]);
fprintf(universe->ulogfile,"\n");

14
src/REPLICA/neb.h
View File

@ -28,8 +28,12 @@ namespace LAMMPS_NS {
class NEB : protected Pointers {
public:
NEB(class LAMMPS *);
NEB(class LAMMPS *, double, double, int, int, int, double *, double *);
~NEB();
void command(int, char **);
void command(int, char **); // process neb command
void run(); // run NEB
double ebf,ebr; // forward and reverse energy barriers
private:
int me,me_universe; // my proc ID in world and universe
@ -38,9 +42,15 @@ class NEB : protected Pointers {
MPI_Comm roots; // MPI comm with 1 root proc from each world
FILE *fp;
int compressed;
double etol; // energy tolerance convergence criterion
double ftol; // force tolerance convergence criterion
int n1steps, n2steps; // number of steps in stage 1 and 2
int nevery; // output interval
char *infile; // name of file containing final state
class FixNEB *fneb;
double **all; // PE,plen,nlen from each replica
int nall; // per-replica dimension of array all
double **all; // PE,plen,nlen,gradvnorm from each replica
double *rdist; // normalize reaction distance, 0 to 1
void readfile(char *);

10
src/REPLICA/prd.cpp
View File

@ -34,7 +34,7 @@
#include "modify.h"
#include "compute.h"
#include "fix.h"
#include "fix_event.h"
#include "fix_event_prd.h"
#include "force.h"
#include "pair.h"
#include "random_park.h"
@ -162,13 +162,13 @@ void PRD::command(int narg, char **arg)
args[1] = (char *) dist_setting;
if (dist_setting) velocity->options(2,args);
// create FixEvent class to store event and pre-quench states
// create FixEventPRD class to store event and pre-quench states
args[0] = (char *) "prd_event";
args[1] = (char *) "all";
args[2] = (char *) "EVENT";
args[2] = (char *) "EVENT/PRD";
modify->add_fix(3,args);
fix_event = (FixEvent *) modify->fix[modify->nfix-1];
fix_event = (FixEventPRD *) modify->fix[modify->nfix-1];
// create Finish for timing output
@ -180,7 +180,7 @@ void PRD::command(int narg, char **arg)
delete [] loop_setting;
delete [] dist_setting;
// assign FixEvent to event-detection compute
// assign FixEventPRD to event-detection compute
// necessary so it will know atom coords at last event
int icompute = modify->find_compute(id_compute);

2
src/REPLICA/prd.h
View File

@ -54,7 +54,7 @@ class PRD : protected Pointers {
class RanPark *random_select;
class RanMars *random_dephase;
class Compute *compute_event;
class FixEvent *fix_event;
class FixEventPRD *fix_event;
class Velocity *velocity;
class Compute *temperature;
class Finish *finish;

987
src/REPLICA/tad.cpp Normal file
View File

@ -0,0 +1,987 @@
// To do:
// Mysterious problem with // if (universe->iworld == 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "tad.h"
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "velocity.h"
#include "integrate.h"
#include "min.h"
#include "neighbor.h"
#include "modify.h"
#include "neb.h"
#include "compute.h"
#include "fix.h"
#include "fix_event_tad.h"
#include "fix_store_state.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;
/* ---------------------------------------------------------------------- */
TAD::TAD(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
TAD::~TAD() {
memory->sfree(fix_event_list);
if (neb_logfilename != NULL) delete [] neb_logfilename;
delete [] min_style;
delete [] min_style_neb;
}
/* ----------------------------------------------------------------------
perform TAD simulation on root proc
other procs only used for NEB calcs
------------------------------------------------------------------------- */
void TAD::command(int narg, char **arg)
{
fix_event_list = NULL;
n_event_list = 0;
nmax_event_list = 0;
nmin_event_list = 10;
// error checks
if (domain->box_exist == 0)
error->all("tad command before simulation box is defined");
if (universe->nworlds == 1) error->all("Cannot use TAD with a single replica for NEB");
if (universe->nworlds != universe->nprocs)
error->all("Can only use TAD with 1-processor replicas for NEB");
if (atom->sortfreq > 0)
error->all("Cannot use TAD with atom_modify sort enabled for NEB");
if (atom->map_style == 0)
error->all("Cannot use TAD unless atom map exists for NEB");
if (narg < 7) error->universe_all("Illegal tad command");
nsteps = atoi(arg[0]);
t_event = atoi(arg[1]);
templo = atof(arg[2]);
temphi = atof(arg[3]);
delta_conf = atof(arg[4]);
tmax = atof(arg[5]);
char id_compute[strlen(arg[6])+1];
strcpy(id_compute,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 tad command");
if (nsteps % t_event)
error->universe_all("TAD nsteps must be multiple of t_event");
if (delta_conf <= 0.0 || delta_conf >= 1.0)
error->universe_all("Invalid delta_conf in tad command");
if (tmax <= 0.0)
error->universe_all("Invalid tmax in tad command");
// deltconf = (ln(1/delta))/freq_min (timestep units)
deltconf = -log(delta_conf)*tmax/update->dt;
// local storage
int me_universe = universe->me;
int nprocs_universe = universe->nprocs;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
delta_beta = (1.0/templo - 1.0/temphi) / force->boltz;
ratio_beta = templo/temphi;
// create FixEventTAD object to store last event
int narg2 = 3;
char **args = new char*[narg2];
args[0] = (char *) "tad_event";
args[1] = (char *) "all";
args[2] = (char *) "EVENT/TAD";
modify->add_fix(narg2,args);
fix_event = (FixEventTAD *) modify->fix[modify->nfix-1];
delete [] args;
// create FixStoreState object to store revert state
narg2 = 13;
args = new char*[narg2];
args[0] = (char *) "tad_revert";
args[1] = (char *) "all";
args[2] = (char *) "store/state";
args[3] = (char *) "0";
args[4] = (char *) "x";
args[5] = (char *) "y";
args[6] = (char *) "z";
args[7] = (char *) "ix";
args[8] = (char *) "iy";
args[9] = (char *) "iz";
args[10] = (char *) "vx";
args[11] = (char *) "vy";
args[12] = (char *) "vz";
modify->add_fix(narg2,args);
fix_revert = (FixStoreState *) modify->fix[modify->nfix-1];
delete [] args;
// create Finish for timing output
finish = new Finish(lmp);
// assign FixEventTAD 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 TAD");
compute_event = modify->compute[icompute];
compute_event->reset_extra_compute_fix("tad_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_universe == 0)
error->warning("Resetting reneighboring criteria during TAD");
}
neighbor->every = 1;
neighbor->delay = 0;
neighbor->dist_check = 1;
// initialize TAD 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();
// set minimize style for quench
narg2 = 1;
args = new char*[narg2];
args[0] = min_style;
update->create_minimize(narg2,args);
delete [] args;
// init minimizer settings and minimizer itself
update->etol = etol;
update->ftol = ftol;
update->max_eval = maxeval;
update->minimize->init();
// perform TAD simulation
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up TAD ...\n");
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step CPU Clock Event "
"\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step CPU Clock Event "
"\n");
}
ulogfile_lammps = universe->ulogfile;
uscreen_lammps = universe->uscreen;
ulogfile_neb = NULL;
uscreen_neb = NULL;
if (me_universe == 0 && neb_logfilename)
ulogfile_neb = fopen(neb_logfilename,"w");
// store hot state and quenched event, only on replica 0
// need this line if quench() does only setup_minimal()
// update->minimize->setup();
// This should work with if uncommented, but does not
// if (universe->iworld == 0) {
fix_event->store_state();
quench();
timer->barrier_start(TIME_LOOP);
time_start = timer->array[TIME_LOOP];
fix_event->store_event(update->ntimestep);
log_event();
fix_event->restore_state();
// do full init/setup
update->whichflag = 1;
lmp->init();
update->integrate->setup();
// }
// main loop: look for events until out of time
// (1) dynamics, store state, quench, check event, restore state
// (2) if event, perform NEB, record in fix_event_list
// (3) if confident, pick earliest event
nbuild = ndanger = 0;
time_neb = time_dynamics = time_quench = time_comm = time_output = 0.0;
timer->barrier_start(TIME_LOOP);
time_start = timer->array[TIME_LOOP];
int confident_flag, event_flag;
if (universe->iworld == 0) {
while (update->ntimestep < update->endstep) {
// initialize list of possible events
initialize_event_list();
confident_flag = 0;
while (update->ntimestep < update->endstep) {
event_flag = 0;
while (update->ntimestep < update->endstep) {
dynamics();
fix_event->store_state();
quench();
event_flag = check_event();
MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);
if (event_flag) break;
// restore hot state
fix_event->restore_state();
// store hot state in revert
fix_revert->end_of_step();
}
if (!event_flag) break;
add_event();
perform_neb(n_event_list-1);
compute_tlo(n_event_list-1);
confident_flag = check_confidence();
MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
if (confident_flag) break;
if (universe->iworld == 0) revert();
}
if (!confident_flag) break;
perform_event(event_first);
// need to sync timestep with TAD
MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
int restart_flag = 0;
if (output->restart_every && universe->iworld == 0)
if (fix_event->event_number % output->restart_every == 0)
restart_flag = 1;
// full init/setup since are starting after event
update->whichflag = 1;
lmp->init();
update->integrate->setup();
// 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];
}
}
} else {
while (update->ntimestep < update->endstep) {
confident_flag = 0;
while (update->ntimestep < update->endstep) {
event_flag = 0;
while (update->ntimestep < update->endstep) {
update->ntimestep += t_event;
MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);
if (event_flag) break;
}
if (!event_flag) break;
perform_neb(-1);
MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
if (confident_flag) break;
}
if (!confident_flag) break;
// need to sync timestep with TAD
MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
}
}
// 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_neb;
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);
}
if (me_universe == 0) fclose(ulogfile_neb);
finish->end(3);
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;
delete finish;
modify->delete_fix("tad_event");
modify->delete_fix("tad_revert");
delete_event_list();
compute_event->reset_extra_compute_fix(NULL);
}
/* ----------------------------------------------------------------------
single short dynamics run
------------------------------------------------------------------------- */
void TAD::dynamics()
{
update->whichflag = 1;
update->nsteps = t_event;
lmp->init();
update->integrate->setup();
// this may be needed if don't do full init
//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 TAD::quench()
{
int ntimestep_hold = update->ntimestep;
int endstep_hold = update->endstep;
// need to change whichflag so that minimize->setup() calling
// modify->setup() will call fix->min_setup()
update->whichflag = 2;
update->nsteps = maxiter;
update->endstep = update->laststep = update->firststep + maxiter;
// full init works
lmp->init();
update->minimize->setup();
// partial init does 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(1);
// reset timestep as if quench did not occur
// clear timestep storage from computes, since now invalid
update->ntimestep = ntimestep_hold;
update->endstep = update->laststep = endstep_hold;
for (int i = 0; i < modify->ncompute; i++)
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
}
/* ----------------------------------------------------------------------
check for an event
return 0 if no event
return 1 if event
------------------------------------------------------------------------- */
int TAD::check_event()
{
int flag;
flag = 0;
if (compute_event->compute_scalar() > 0.0) flag = 1;
return flag;
}
/* ----------------------------------------------------------------------
universe proc 0 prints event info
------------------------------------------------------------------------- */
void TAD::log_event()
{
timer->array[TIME_LOOP] = time_start;
if (universe->me == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"%d %.3f %.3f %d\n",
fix_event->event_timestep,
timer->elapsed(TIME_LOOP),
fix_event->tlo,
fix_event->event_number);
if (universe->ulogfile)
fprintf(universe->ulogfile,"%d %.3f %.3f %d\n",
fix_event->event_timestep,
timer->elapsed(TIME_LOOP),
fix_event->tlo,
fix_event->event_number);
}
// dump snapshot of quenched coords
// must reneighbor and compute forces before dumping
// 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];
}
}
/* ----------------------------------------------------------------------
parse optional parameters at end of TAD input line
------------------------------------------------------------------------- */
void TAD::options(int narg, char **arg)
{
if (narg < 0) error->all("Illegal tad command");
// set defaults
etol = 0.1;
ftol = 0.1;
maxiter = 40;
maxeval = 50;
etol_neb = 0.01;
ftol_neb = 0.01;
n1steps_neb = 100;
n2steps_neb = 100;
nevery_neb = 10;
min_style = new char[3];
strcpy(min_style,"cg");
min_style_neb = new char[9];
strcpy(min_style_neb,"quickmin");
neb_logfilename = NULL;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"min") == 0) {
if (iarg+5 > narg) error->all("Illegal tad command");
etol = atof(arg[iarg+1]);
ftol = atof(arg[iarg+2]);
maxiter = atoi(arg[iarg+3]);
maxeval = atoi(arg[iarg+4]);
if (maxiter < 0 || maxeval < 0 ||
etol < 0.0 || ftol < 0.0 )
error->all("Illegal tad command");
iarg += 5;
} else if (strcmp(arg[iarg],"neb") == 0) {
if (iarg+6 > narg) error->all("Illegal tad command");
etol_neb = atof(arg[iarg+1]);
ftol_neb = atof(arg[iarg+2]);
n1steps_neb = atoi(arg[iarg+3]);
n2steps_neb = atoi(arg[iarg+4]);
nevery_neb = atoi(arg[iarg+5]);
if (etol_neb < 0.0 || ftol_neb < 0.0 ||
n1steps_neb < 0 || n2steps_neb < 0 ||
nevery_neb < 0) error->all("Illegal tad command");
iarg += 6;
} else if (strcmp(arg[iarg],"min_style") == 0) {
if (iarg+2 > narg) error->all("Illegal tad command");
int n = strlen(arg[iarg+1]) + 1;
delete [] min_style;
min_style = new char[n];
strcpy(min_style,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"neb_style") == 0) {
if (iarg+2 > narg) error->all("Illegal tad command");
int n = strlen(arg[iarg+1]) + 1;
delete [] min_style_neb;
min_style_neb = new char[n];
strcpy(min_style_neb,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"neb_log") == 0) {
delete [] neb_logfilename;
if (iarg+2 > narg) error->all("Illegal tad command");
if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = NULL;
else {
int n = strlen(arg[iarg+1]) + 1;
neb_logfilename = new char[n];
strcpy(neb_logfilename,arg[iarg+1]);
}
iarg += 2;
} else error->all("Illegal tad command");
}
}
/* ----------------------------------------------------------------------
perform NEB calculation
------------------------------------------------------------------------- */
void TAD::perform_neb(int ievent)
{
double **x = atom->x;
int nlocal = atom->nlocal;
double *buf_final = (double *)
memory->smalloc(3*nlocal*sizeof(double),"tad:buffinal");
// set system to quenched state of event ievent
if (universe->iworld == 0) {
fix_event_list[ievent]->restore_event();
int ii = 0;
for (int i = 0; i < nlocal; i++) {
buf_final[ii++] = x[i][0];
buf_final[ii++] = x[i][1];
buf_final[ii++] = x[i][2];
}
}
MPI_Bcast(buf_final,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
universe->uworld);
double *buf_init = (double *)
memory->smalloc(3*nlocal*sizeof(double),"tad:bufinit");
// set system to quenched state of fix_event
if (universe->iworld == 0) {
fix_event->restore_event();
int ii = 0;
for (int i = 0; i < nlocal; i++) {
buf_init[ii++] = x[i][0];
buf_init[ii++] = x[i][1];
buf_init[ii++] = x[i][2];
}
}
MPI_Bcast(buf_init,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
universe->uworld);
// create FixNEB object to support NEB
int narg2 = 4;
char **args = new char*[narg2];
args[0] = (char *) "neb";
args[1] = (char *) "all";
args[2] = (char *) "neb";
char str[128];
args[3] = str;
double kspring = 1.0;
sprintf(args[3],"%f",kspring);
modify->add_fix(narg2,args);
fix_neb = (Fix *) modify->fix[modify->nfix-1];
delete [] args;
// switch minimize style to quickmin for NEB
narg2 = 1;
args = new char*[narg2];
args[0] = min_style_neb;
update->create_minimize(narg2,args);
delete [] args;
// create NEB object
neb = new NEB(lmp,etol_neb,ftol_neb,n1steps_neb,
n2steps_neb,nevery_neb,buf_init,buf_final);
// free up temporary arrays
memory->sfree(buf_init);
memory->sfree(buf_final);
// Run NEB
double beginstep_hold = update->beginstep;
double endstep_hold = update->endstep;
double ntimestep_hold = update->ntimestep;
double nsteps_hold = update->nsteps;
if (universe->me == 0) {
universe->ulogfile = ulogfile_neb;
universe->uscreen = uscreen_neb;
}
// Had to bypass timer interface
// because timer->array is reset
// inside neb->run()
// timer->barrier_start(TIME_LOOP);
// neb->run();
// timer->barrier_stop(TIME_LOOP);
// time_neb += timer->array[TIME_LOOP];
MPI_Barrier(world);
double time_tmp = MPI_Wtime();
neb->run();
MPI_Barrier(world);
time_neb += MPI_Wtime() - time_tmp;
if (universe->me == 0) {
universe->ulogfile = ulogfile_lammps;
universe->uscreen = uscreen_lammps;
}
// Extract barrier energy from NEB
if (universe->iworld == 0)
fix_event_list[ievent]->ebarrier = neb->ebf;
update->beginstep = update->firststep = beginstep_hold;
update->endstep = update->laststep = endstep_hold;
update->ntimestep = ntimestep_hold;
update->nsteps = nsteps_hold;
// switch minimize style back for quench
narg2 = 1;
args = new char*[narg2];
args[0] = min_style;
update->create_minimize(narg2,args);
update->etol = etol;
update->ftol = ftol;
delete [] args;
// Clean up
modify->delete_fix("neb");
delete neb;
}
/* ----------------------------------------------------------------------
check if confidence criterion for tstop is satisfied
return 0 if not satisfied
return 1 if satisfied
------------------------------------------------------------------------- */
int TAD::check_confidence()
{
int flag;
// update stopping time
deltstop = deltconf*pow(deltfirst/deltconf, ratio_beta);
flag = 0;
if (deltstop < update->ntimestep - fix_event->event_timestep) flag = 1;
return flag;
}
/* ----------------------------------------------------------------------
reflect back in to starting state
------------------------------------------------------------------------- */
void TAD::revert()
{
double **x = atom->x;
double **v = atom->v;
int *image = atom->image;
int nlocal = atom->nlocal;
double **array_atom = fix_revert->array_atom;
for (int i = 0; i < nlocal; i++) {
x[i][0] = array_atom[i][0];
x[i][1] = array_atom[i][1];
x[i][2] = array_atom[i][2];
image[i] = ((int(array_atom[i][5]) + 512 & 1023) << 20) |
((int(array_atom[i][4]) + 512 & 1023) << 10) |
(int(array_atom[i][3]) + 512 & 1023);
v[i][0] = -array_atom[i][6];
v[i][1] = -array_atom[i][7];
v[i][2] = -array_atom[i][8];
}
}
/* ----------------------------------------------------------------------
Initialize list of possible events
------------------------------------------------------------------------- */
void TAD::initialize_event_list() {
// First delete old events, if any
delete_event_list();
// Create new list
n_event_list = 0;
grow_event_list(nmin_event_list);
}
/* ----------------------------------------------------------------------
Delete list of possible events
------------------------------------------------------------------------- */
void TAD::delete_event_list() {
for (int i = 0; i < n_event_list; i++) {
char str[128];
sprintf(str,"tad_event_%d",i);
modify->delete_fix(str);
}
memory->sfree(fix_event_list);
fix_event_list = NULL;
n_event_list = 0;
nmax_event_list = 0;
}
/* ----------------------------------------------------------------------
add event
------------------------------------------------------------------------- */
void TAD::add_event()
{
// create FixEventTAD object to store possible event
int narg = 3;
char **args = new char*[narg];
char str[128];
sprintf(str,"tad_event_%d",n_event_list);
args[0] = str;
args[1] = (char *) "all";
args[2] = (char *) "EVENT/TAD";
modify->add_fix(narg,args);
if (n_event_list == nmax_event_list)
grow_event_list(nmax_event_list+nmin_event_list);
n_event_list += 1;
int ievent = n_event_list-1;
fix_event_list[ievent] = (FixEventTAD *) modify->fix[modify->nfix-1];
// store quenched state for new event
fix_event_list[ievent]->store_event(update->ntimestep);
// store hot state for new event
fix_event->restore_state();
fix_event_list[ievent]->store_state();
// string clean-up
delete [] args;
}
/* ----------------------------------------------------------------------
compute cold time for event ievent
------------------------------------------------------------------------- */
void TAD::compute_tlo(int ievent)
{
double deltlo,delthi,ebarrier;
ebarrier = fix_event_list[ievent]->ebarrier;
delthi = fix_event_list[ievent]->event_timestep
- fix_event->event_timestep;
deltlo = delthi*exp(ebarrier*delta_beta);
fix_event_list[ievent]->tlo = fix_event->tlo + deltlo;
// first-replica output about each event
if (universe->me == 0) {
char str[128];
double tfrac = 0.0;
if (ievent > 0) tfrac = delthi/deltstop;
// sprintf(str,
// "New event: ievent = %d eb = %g t_lo = %g t_hi = %g t_hi/t_stop = %g \n",
// ievent,ebarrier,deltlo,delthi,tfrac);
// error->warning(str);
if (screen)
fprintf(screen,
"New event: t_hi = %d ievent = %d eb = %g dt_lo = %g dt_hi/t_stop = %g \n",
fix_event_list[ievent]->event_timestep,
ievent,ebarrier,deltlo,tfrac);
if (logfile)
fprintf(logfile,
"New event: t_hi = %d ievent = %d eb = %g dt_lo = %g dt_hi/t_stop = %g \n",
fix_event_list[ievent]->event_timestep,
ievent,ebarrier,deltlo,tfrac);
}
// update first event
if (ievent == 0) {
deltfirst = deltlo;
event_first = ievent;
} else if (deltlo < deltfirst) {
deltfirst = deltlo;
event_first = ievent;
}
}
/* ----------------------------------------------------------------------
perform event
------------------------------------------------------------------------- */
void TAD::perform_event(int ievent)
{
// reset timestep to that of event
update->ntimestep = fix_event_list[ievent]->event_timestep;
// Copy event to current event
fix_event->tlo = fix_event_list[ievent]->tlo;
fix_event->event_number++;
fix_event->event_timestep = update->ntimestep;
fix_event_list[ievent]->restore_event();
fix_event->store_event(fix_event_list[ievent]->event_timestep);
// output stats and dump for quenched state
log_event();
// load and store hot state
fix_event_list[ievent]->restore_state();
fix_event->store_state();
}
/* ----------------------------------------------------------------------
Allocate list of pointers to events
------------------------------------------------------------------------- */
void TAD::grow_event_list(int nmax) {
if (nmax_event_list > nmax) return;
fix_event_list = (FixEventTAD **) memory->srealloc(fix_event_list,nmax*sizeof(FixEventTAD *),"tad:eventlist");
nmax_event_list = nmax;
}

88
src/REPLICA/tad.h Normal file
View File

@ -0,0 +1,88 @@
/* ----------------------------------------------------------------------
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 COMMAND_CLASS
CommandStyle(tad,TAD)
#else
#ifndef LMP_TAD_H
#define LMP_TAD_H
#include "pointers.h"
namespace LAMMPS_NS {
class TAD : protected Pointers {
public:
TAD(class LAMMPS *);
~TAD();
void command(int, char **);
private:
int me,nprocs;
int nsteps,t_event;
double templo,temphi,delta_conf,tmax;
double etol,ftol,etol_neb,ftol_neb;
int maxiter,maxeval,n1steps_neb,n2steps_neb,nevery_neb;
char *min_style, *min_style_neb;
double delta_beta,ratio_beta;
double deltconf,deltstop,deltfirst; // Times since last event
int event_first;
int neigh_every,neigh_delay,neigh_dist_check;
int nbuild,ndanger;
int quench_reneighbor;
double time_dynamics,time_quench,time_neb,time_comm,time_output;
double time_start;
class NEB *neb; // NEB object
class Fix *fix_neb; // FixNEB object
class Compute *compute_event; // compute to detect event
class FixEventTAD *fix_event; // current event/state
class FixStoreState *fix_revert; // revert state
FixEventTAD **fix_event_list; // list of possible events
int n_event_list; // number of events
int nmax_event_list; // allocated events
int nmin_event_list; // minimum allocation
char *neb_logfilename; // filename for ulogfile_neb
FILE *uscreen_neb; // neb universe screen output
FILE *ulogfile_neb; // neb universe logfile
FILE *uscreen_lammps; // lammps universe screen output
FILE *ulogfile_lammps; // lammps universe logfile
class Finish *finish;
void dynamics();
void quench();
int check_event();
int check_confidence();
void perform_neb(int);
void log_event();
void options(int, char **);
void revert();
void add_event();
void perform_event(int);
void compute_tlo(int);
void grow_event_list(int);
void initialize_event_list();
void delete_event_list();
};
}
#endif
#endif