forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5719 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
2bd4074d8c
commit
2204858474
|
|
@ -1,255 +1,255 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "lmptype.h"
|
||||
#include "mpi.h"
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_qeq_comb.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "respa.h"
|
||||
#include "pair_comb.h"
|
||||
#include "update.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
|
||||
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 5) error->all("Illegal fix qeq/comb command");
|
||||
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = 0;
|
||||
peratom_freq = 1;
|
||||
|
||||
nevery = force->inumeric(arg[3]);
|
||||
precision = force->numeric(arg[4]);
|
||||
|
||||
if (nevery <= 0 || precision <= 0.0)
|
||||
error->all("Illegal fix qeq/comb command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
// optional args
|
||||
|
||||
fp = NULL;
|
||||
|
||||
int iarg = 5;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"file") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix qeq/comb command");
|
||||
if (me == 0) {
|
||||
fp = fopen(arg[iarg+1],"w");
|
||||
if (fp == NULL) {
|
||||
char str[128];
|
||||
sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
|
||||
error->one(str);
|
||||
}
|
||||
}
|
||||
iarg += 2;
|
||||
} else error->all("Illegal fix qeq/comb command");
|
||||
}
|
||||
|
||||
nmax = atom->nmax;
|
||||
qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
|
||||
q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
|
||||
q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
|
||||
vector_atom = qf;
|
||||
|
||||
// zero the vector since dump may access it on timestep 0
|
||||
// zero the vector since a variable may access it before first run
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEQComb::~FixQEQComb()
|
||||
{
|
||||
if (me == 0 && fp) fclose(fp);
|
||||
memory->sfree(qf);
|
||||
memory->sfree(q1);
|
||||
memory->sfree(q2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixQEQComb::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
mask |= POST_FORCE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::init()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all("Fix qeq/comb requires atom attribute q");
|
||||
|
||||
comb = (PairComb *) force->pair_match("comb",1);
|
||||
if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");
|
||||
|
||||
if (strcmp(update->integrate_style,"respa") == 0)
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
ngroup = group->count(igroup);
|
||||
if (ngroup == 0) error->all("Fix qeq/comb group has no atoms");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::setup(int vflag)
|
||||
{
|
||||
firstflag = 1;
|
||||
if (strcmp(update->integrate_style,"verlet") == 0)
|
||||
post_force(vflag);
|
||||
else {
|
||||
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
|
||||
post_force_respa(vflag,nlevels_respa-1,0);
|
||||
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
|
||||
}
|
||||
firstflag = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::post_force(int vflag)
|
||||
{
|
||||
int i,iloop,loopmax;
|
||||
double heatpq,qmass,dtq,dtq2;
|
||||
double enegchkall,enegmaxall;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// reallocate work arrays if necessary
|
||||
// qf = charge force
|
||||
// q1 = charge displacement
|
||||
// q2 = tmp storage of charge force for next iteration
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->sfree(qf);
|
||||
memory->sfree(q1);
|
||||
memory->sfree(q2);
|
||||
nmax = atom->nmax;
|
||||
qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
|
||||
q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
|
||||
q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
|
||||
vector_atom = qf;
|
||||
}
|
||||
|
||||
// more loops for first-time charge equilibrium
|
||||
|
||||
iloop = 0;
|
||||
if (firstflag) loopmax = 5000;
|
||||
else loopmax = 2000;
|
||||
|
||||
// charge-equilibration loop
|
||||
|
||||
if (me == 0 && fp)
|
||||
fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
|
||||
update->ntimestep);
|
||||
|
||||
heatpq = 0.05;
|
||||
qmass = 0.000548580;
|
||||
dtq = 0.0006;
|
||||
dtq2 = 0.5*dtq*dtq/qmass;
|
||||
|
||||
double enegchk = 0.0;
|
||||
double enegtot = 0.0;
|
||||
double enegmax = 0.0;
|
||||
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
q1[i] = q2[i] = qf[i] = 0.0;
|
||||
|
||||
for (iloop = 0; iloop < loopmax; iloop ++ ) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q1[i] += qf[i]*dtq2 - heatpq*q1[i];
|
||||
q[i] += q1[i];
|
||||
}
|
||||
|
||||
enegtot = comb->yasu_char(qf,igroup);
|
||||
enegtot /= ngroup;
|
||||
enegchk = enegmax = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal ; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q2[i] = enegtot-qf[i];
|
||||
enegmax = MAX(enegmax,fabs(q2[i]));
|
||||
enegchk += fabs(q2[i]);
|
||||
qf[i] = q2[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
enegchk = enegchkall/ngroup;
|
||||
MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
|
||||
enegmax = enegmaxall;
|
||||
|
||||
if (enegchk <= precision && enegmax <= 100.0*precision) break;
|
||||
|
||||
if (me == 0 && fp)
|
||||
fprintf(fp," iteration: %d, enegtot %.6g, "
|
||||
"enegmax %.6g, fq deviation: %.6g\n",
|
||||
iloop,enegtot,enegmax,enegchk);
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
q1[i] += qf[i]*dtq2 - heatpq*q1[i];
|
||||
}
|
||||
|
||||
if (me == 0 && fp) {
|
||||
if (iloop == loopmax)
|
||||
fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
|
||||
else
|
||||
fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
|
||||
iloop,enegchk);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_force(vflag);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixQEQComb::memory_usage()
|
||||
{
|
||||
double bytes = atom->nmax*3 * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
/* ----------------------------------------------------------------------
|
||||
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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "lmptype.h"
|
||||
#include "mpi.h"
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_qeq_comb.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "respa.h"
|
||||
#include "pair_comb.h"
|
||||
#include "update.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
|
||||
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 5) error->all("Illegal fix qeq/comb command");
|
||||
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = 0;
|
||||
peratom_freq = 1;
|
||||
|
||||
nevery = force->inumeric(arg[3]);
|
||||
precision = force->numeric(arg[4]);
|
||||
|
||||
if (nevery <= 0 || precision <= 0.0)
|
||||
error->all("Illegal fix qeq/comb command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
// optional args
|
||||
|
||||
fp = NULL;
|
||||
|
||||
int iarg = 5;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"file") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix qeq/comb command");
|
||||
if (me == 0) {
|
||||
fp = fopen(arg[iarg+1],"w");
|
||||
if (fp == NULL) {
|
||||
char str[128];
|
||||
sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
|
||||
error->one(str);
|
||||
}
|
||||
}
|
||||
iarg += 2;
|
||||
} else error->all("Illegal fix qeq/comb command");
|
||||
}
|
||||
|
||||
nmax = atom->nmax;
|
||||
qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
|
||||
q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
|
||||
q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
|
||||
vector_atom = qf;
|
||||
|
||||
// zero the vector since dump may access it on timestep 0
|
||||
// zero the vector since a variable may access it before first run
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEQComb::~FixQEQComb()
|
||||
{
|
||||
if (me == 0 && fp) fclose(fp);
|
||||
memory->sfree(qf);
|
||||
memory->sfree(q1);
|
||||
memory->sfree(q2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixQEQComb::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
mask |= POST_FORCE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::init()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all("Fix qeq/comb requires atom attribute q");
|
||||
|
||||
comb = (PairComb *) force->pair_match("comb",1);
|
||||
if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");
|
||||
|
||||
if (strcmp(update->integrate_style,"respa") == 0)
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
ngroup = group->count(igroup);
|
||||
if (ngroup == 0) error->all("Fix qeq/comb group has no atoms");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::setup(int vflag)
|
||||
{
|
||||
firstflag = 1;
|
||||
if (strcmp(update->integrate_style,"verlet") == 0)
|
||||
post_force(vflag);
|
||||
else {
|
||||
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
|
||||
post_force_respa(vflag,nlevels_respa-1,0);
|
||||
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
|
||||
}
|
||||
firstflag = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::post_force(int vflag)
|
||||
{
|
||||
int i,iloop,loopmax;
|
||||
double heatpq,qmass,dtq,dtq2;
|
||||
double enegchkall,enegmaxall;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// reallocate work arrays if necessary
|
||||
// qf = charge force
|
||||
// q1 = charge displacement
|
||||
// q2 = tmp storage of charge force for next iteration
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->sfree(qf);
|
||||
memory->sfree(q1);
|
||||
memory->sfree(q2);
|
||||
nmax = atom->nmax;
|
||||
qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
|
||||
q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
|
||||
q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
|
||||
vector_atom = qf;
|
||||
}
|
||||
|
||||
// more loops for first-time charge equilibrium
|
||||
|
||||
iloop = 0;
|
||||
if (firstflag) loopmax = 5000;
|
||||
else loopmax = 2000;
|
||||
|
||||
// charge-equilibration loop
|
||||
|
||||
if (me == 0 && fp)
|
||||
fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
|
||||
update->ntimestep);
|
||||
|
||||
heatpq = 0.05;
|
||||
qmass = 0.000548580;
|
||||
dtq = 0.0006;
|
||||
dtq2 = 0.5*dtq*dtq/qmass;
|
||||
|
||||
double enegchk = 0.0;
|
||||
double enegtot = 0.0;
|
||||
double enegmax = 0.0;
|
||||
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
q1[i] = q2[i] = qf[i] = 0.0;
|
||||
|
||||
for (iloop = 0; iloop < loopmax; iloop ++ ) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q1[i] += qf[i]*dtq2 - heatpq*q1[i];
|
||||
q[i] += q1[i];
|
||||
}
|
||||
|
||||
enegtot = comb->yasu_char(qf,igroup);
|
||||
enegtot /= ngroup;
|
||||
enegchk = enegmax = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal ; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q2[i] = enegtot-qf[i];
|
||||
enegmax = MAX(enegmax,fabs(q2[i]));
|
||||
enegchk += fabs(q2[i]);
|
||||
qf[i] = q2[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
enegchk = enegchkall/ngroup;
|
||||
MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
|
||||
enegmax = enegmaxall;
|
||||
|
||||
if (enegchk <= precision && enegmax <= 100.0*precision) break;
|
||||
|
||||
if (me == 0 && fp)
|
||||
fprintf(fp," iteration: %d, enegtot %.6g, "
|
||||
"enegmax %.6g, fq deviation: %.6g\n",
|
||||
iloop,enegtot,enegmax,enegchk);
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
q1[i] += qf[i]*dtq2 - heatpq*q1[i];
|
||||
}
|
||||
|
||||
if (me == 0 && fp) {
|
||||
if (iloop == loopmax)
|
||||
fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
|
||||
else
|
||||
fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
|
||||
iloop,enegchk);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_force(vflag);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixQEQComb::memory_usage()
|
||||
{
|
||||
double bytes = atom->nmax*3 * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
592
src/update.cpp
592
src/update.cpp
|
|
@ -1,296 +1,296 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "update.h"
|
||||
#include "style_integrate.h"
|
||||
#include "style_minimize.h"
|
||||
#include "neighbor.h"
|
||||
#include "force.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "domain.h"
|
||||
#include "region.h"
|
||||
#include "compute.h"
|
||||
#include "output.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
Update::Update(LAMMPS *lmp) : Pointers(lmp)
|
||||
{
|
||||
int n;
|
||||
char *str;
|
||||
|
||||
ntimestep = 0;
|
||||
first_update = 0;
|
||||
|
||||
whichflag = 0;
|
||||
firststep = laststep = 0;
|
||||
beginstep = endstep = 0;
|
||||
setupflag = 0;
|
||||
multireplica = 0;
|
||||
|
||||
restrict_output = 0;
|
||||
|
||||
eflag_global = vflag_global = -1;
|
||||
|
||||
unit_style = NULL;
|
||||
set_units("lj");
|
||||
|
||||
str = (char *) "verlet";
|
||||
n = strlen(str) + 1;
|
||||
integrate_style = new char[n];
|
||||
strcpy(integrate_style,str);
|
||||
integrate = new Verlet(lmp,0,NULL);
|
||||
|
||||
str = (char *) "cg";
|
||||
n = strlen(str) + 1;
|
||||
minimize_style = new char[n];
|
||||
strcpy(minimize_style,str);
|
||||
minimize = new MinCG(lmp);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
Update::~Update()
|
||||
{
|
||||
delete [] unit_style;
|
||||
|
||||
delete [] integrate_style;
|
||||
delete integrate;
|
||||
|
||||
delete [] minimize_style;
|
||||
delete minimize;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::init()
|
||||
{
|
||||
// init the appropriate integrate and/or minimize class
|
||||
// if neither (e.g. from write_restart) then just return
|
||||
|
||||
if (whichflag == 0) return;
|
||||
if (whichflag == 1) integrate->init();
|
||||
else if (whichflag == 2) minimize->init();
|
||||
|
||||
// only set first_update if a run or minimize is being performed
|
||||
|
||||
first_update = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::set_units(const char *style)
|
||||
{
|
||||
// physical constants from:
|
||||
// http://physics.nist.gov/cuu/Constants/Table/allascii.txt
|
||||
// using thermochemical calorie = 4.184 J
|
||||
|
||||
if (strcmp(style,"lj") == 0) {
|
||||
force->boltz = 1.0;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 0.005;
|
||||
neighbor->skin = 0.3;
|
||||
|
||||
} else if (strcmp(style,"real") == 0) {
|
||||
force->boltz = 0.0019872067;
|
||||
force->mvv2e = 48.88821291 * 48.88821291;
|
||||
force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
|
||||
force->mv2d = 1.0 / 0.602214179;
|
||||
force->nktv2p = 68568.415;
|
||||
force->qqr2e = 332.06371;
|
||||
force->qe2f = 23.060549;
|
||||
force->vxmu2f = 1.4393264316e4;
|
||||
force->xxt2kmu = 0.1;
|
||||
dt = 1.0;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else if (strcmp(style,"metal") == 0) {
|
||||
force->boltz = 8.617343e-5;
|
||||
force->mvv2e = 1.0364269e-4;
|
||||
force->ftm2v = 1.0 / 1.0364269e-4;
|
||||
force->mv2d = 1.0 / 0.602214179;
|
||||
force->nktv2p = 1.6021765e6;
|
||||
force->qqr2e = 14.399645;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 0.6241509647;
|
||||
force->xxt2kmu = 1.0e-4;
|
||||
dt = 0.001;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else if (strcmp(style,"si") == 0) {
|
||||
force->boltz = 1.3806504e-23;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 8.9876e9;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 1.0e-8;
|
||||
neighbor->skin = 0.001;
|
||||
|
||||
} else if (strcmp(style,"cgs") == 0) {
|
||||
force->boltz = 1.3806504e-16;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 1.0e-8;
|
||||
neighbor->skin = 0.1;
|
||||
|
||||
} else if (strcmp(style,"electron") == 0) {
|
||||
force->boltz = 3.16681534e-6;
|
||||
force->mvv2e = 1.06657236;
|
||||
force->ftm2v = 0.937582899;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 2.94210108e13;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.94469051e-10;
|
||||
force->vxmu2f = 3.39893149e1;
|
||||
force->xxt2kmu = 3.13796367e-2;
|
||||
dt = 0.001;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else error->all("Illegal units command");
|
||||
|
||||
delete [] unit_style;
|
||||
int n = strlen(style) + 1;
|
||||
unit_style = new char[n];
|
||||
strcpy(unit_style,style);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::create_integrate(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1) error->all("Illegal run_style command");
|
||||
|
||||
delete [] integrate_style;
|
||||
delete integrate;
|
||||
|
||||
if (0) return; // dummy line to enable else-if macro expansion
|
||||
|
||||
#define INTEGRATE_CLASS
|
||||
#define IntegrateStyle(key,Class) \
|
||||
else if (strcmp(arg[0],#key) == 0) integrate = new Class(lmp,narg-1,&arg[1]);
|
||||
#include "style_integrate.h"
|
||||
#undef INTEGRATE_CLASS
|
||||
|
||||
else error->all("Illegal run_style command");
|
||||
|
||||
int n = strlen(arg[0]) + 1;
|
||||
integrate_style = new char[n];
|
||||
strcpy(integrate_style,arg[0]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::create_minimize(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1) error->all("Illegal min_style command");
|
||||
|
||||
delete [] minimize_style;
|
||||
delete minimize;
|
||||
|
||||
if (0) return; // dummy line to enable else-if macro expansion
|
||||
|
||||
#define MINIMIZE_CLASS
|
||||
#define MinimizeStyle(key,Class) \
|
||||
else if (strcmp(arg[0],#key) == 0) minimize = new Class(lmp);
|
||||
#include "style_minimize.h"
|
||||
#undef MINIMIZE_CLASS
|
||||
|
||||
else error->all("Illegal min_style command");
|
||||
|
||||
int n = strlen(arg[0]) + 1;
|
||||
minimize_style = new char[n];
|
||||
strcpy(minimize_style,arg[0]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
reset timestep from input script
|
||||
do not allow dump files or a restart to be defined
|
||||
do not allow any timestep-dependent fixes to be defined
|
||||
do not allow any dynamic regions to be defined
|
||||
reset eflag/vflag global so nothing will think eng/virial are current
|
||||
reset invoked flags of computes,
|
||||
so nothing will think they are current between runs
|
||||
clear timestep list of computes that store future invocation times
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Update::reset_timestep(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1) error->all("Illegal reset_timestep command");
|
||||
|
||||
for (int i = 0; i < output->ndump; i++)
|
||||
if (output->last_dump[i] >= 0)
|
||||
error->all("Cannot reset timestep with dump file already written to");
|
||||
if (output->restart && output->last_restart >= 0)
|
||||
error->all("Cannot reset timestep with restart file already written");
|
||||
|
||||
for (int i = 0; i < modify->nfix; i++)
|
||||
if (modify->fix[i]->time_depend)
|
||||
error->all("Cannot reset timestep with a time-dependent fix defined");
|
||||
|
||||
for (int i = 0; i < domain->nregion; i++)
|
||||
if (domain->regions[i]->dynamic_check())
|
||||
error->all("Cannot reset timestep with a dynamic region defined");
|
||||
|
||||
eflag_global = vflag_global = -1;
|
||||
|
||||
for (int i = 0; i < modify->ncompute; i++) {
|
||||
modify->compute[i]->invoked_scalar = -1;
|
||||
modify->compute[i]->invoked_vector = -1;
|
||||
modify->compute[i]->invoked_array = -1;
|
||||
modify->compute[i]->invoked_peratom = -1;
|
||||
modify->compute[i]->invoked_local = -1;
|
||||
}
|
||||
|
||||
for (int i = 0; i < modify->ncompute; i++)
|
||||
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
|
||||
|
||||
ntimestep = ATOBIGINT(arg[0]);
|
||||
if (ntimestep < 0) error->all("Timestep must be >= 0");
|
||||
if (ntimestep > MAXBIGINT) error->all("Too big a timestep");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of update and integrate/minimize
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double Update::memory_usage()
|
||||
{
|
||||
double bytes = 0.0;
|
||||
if (whichflag == 1) bytes += integrate->memory_usage();
|
||||
else if (whichflag == 2) bytes += minimize->memory_usage();
|
||||
return bytes;
|
||||
}
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "update.h"
|
||||
#include "style_integrate.h"
|
||||
#include "style_minimize.h"
|
||||
#include "neighbor.h"
|
||||
#include "force.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "domain.h"
|
||||
#include "region.h"
|
||||
#include "compute.h"
|
||||
#include "output.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
Update::Update(LAMMPS *lmp) : Pointers(lmp)
|
||||
{
|
||||
int n;
|
||||
char *str;
|
||||
|
||||
ntimestep = 0;
|
||||
first_update = 0;
|
||||
|
||||
whichflag = 0;
|
||||
firststep = laststep = 0;
|
||||
beginstep = endstep = 0;
|
||||
setupflag = 0;
|
||||
multireplica = 0;
|
||||
|
||||
restrict_output = 0;
|
||||
|
||||
eflag_global = vflag_global = -1;
|
||||
|
||||
unit_style = NULL;
|
||||
set_units("lj");
|
||||
|
||||
str = (char *) "verlet";
|
||||
n = strlen(str) + 1;
|
||||
integrate_style = new char[n];
|
||||
strcpy(integrate_style,str);
|
||||
integrate = new Verlet(lmp,0,NULL);
|
||||
|
||||
str = (char *) "cg";
|
||||
n = strlen(str) + 1;
|
||||
minimize_style = new char[n];
|
||||
strcpy(minimize_style,str);
|
||||
minimize = new MinCG(lmp);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
Update::~Update()
|
||||
{
|
||||
delete [] unit_style;
|
||||
|
||||
delete [] integrate_style;
|
||||
delete integrate;
|
||||
|
||||
delete [] minimize_style;
|
||||
delete minimize;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::init()
|
||||
{
|
||||
// init the appropriate integrate and/or minimize class
|
||||
// if neither (e.g. from write_restart) then just return
|
||||
|
||||
if (whichflag == 0) return;
|
||||
if (whichflag == 1) integrate->init();
|
||||
else if (whichflag == 2) minimize->init();
|
||||
|
||||
// only set first_update if a run or minimize is being performed
|
||||
|
||||
first_update = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::set_units(const char *style)
|
||||
{
|
||||
// physical constants from:
|
||||
// http://physics.nist.gov/cuu/Constants/Table/allascii.txt
|
||||
// using thermochemical calorie = 4.184 J
|
||||
|
||||
if (strcmp(style,"lj") == 0) {
|
||||
force->boltz = 1.0;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 0.005;
|
||||
neighbor->skin = 0.3;
|
||||
|
||||
} else if (strcmp(style,"real") == 0) {
|
||||
force->boltz = 0.0019872067;
|
||||
force->mvv2e = 48.88821291 * 48.88821291;
|
||||
force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
|
||||
force->mv2d = 1.0 / 0.602214179;
|
||||
force->nktv2p = 68568.415;
|
||||
force->qqr2e = 332.06371;
|
||||
force->qe2f = 23.060549;
|
||||
force->vxmu2f = 1.4393264316e4;
|
||||
force->xxt2kmu = 0.1;
|
||||
dt = 1.0;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else if (strcmp(style,"metal") == 0) {
|
||||
force->boltz = 8.617343e-5;
|
||||
force->mvv2e = 1.0364269e-4;
|
||||
force->ftm2v = 1.0 / 1.0364269e-4;
|
||||
force->mv2d = 1.0 / 0.602214179;
|
||||
force->nktv2p = 1.6021765e6;
|
||||
force->qqr2e = 14.399645;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 0.6241509647;
|
||||
force->xxt2kmu = 1.0e-4;
|
||||
dt = 0.001;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else if (strcmp(style,"si") == 0) {
|
||||
force->boltz = 1.3806504e-23;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 8.9876e9;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 1.0e-8;
|
||||
neighbor->skin = 0.001;
|
||||
|
||||
} else if (strcmp(style,"cgs") == 0) {
|
||||
force->boltz = 1.3806504e-16;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
dt = 1.0e-8;
|
||||
neighbor->skin = 0.1;
|
||||
|
||||
} else if (strcmp(style,"electron") == 0) {
|
||||
force->boltz = 3.16681534e-6;
|
||||
force->mvv2e = 1.06657236;
|
||||
force->ftm2v = 0.937582899;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 2.94210108e13;
|
||||
force->qqr2e = 1.0;
|
||||
force->qe2f = 1.94469051e-10;
|
||||
force->vxmu2f = 3.39893149e1;
|
||||
force->xxt2kmu = 3.13796367e-2;
|
||||
dt = 0.001;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else error->all("Illegal units command");
|
||||
|
||||
delete [] unit_style;
|
||||
int n = strlen(style) + 1;
|
||||
unit_style = new char[n];
|
||||
strcpy(unit_style,style);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::create_integrate(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1) error->all("Illegal run_style command");
|
||||
|
||||
delete [] integrate_style;
|
||||
delete integrate;
|
||||
|
||||
if (0) return; // dummy line to enable else-if macro expansion
|
||||
|
||||
#define INTEGRATE_CLASS
|
||||
#define IntegrateStyle(key,Class) \
|
||||
else if (strcmp(arg[0],#key) == 0) integrate = new Class(lmp,narg-1,&arg[1]);
|
||||
#include "style_integrate.h"
|
||||
#undef INTEGRATE_CLASS
|
||||
|
||||
else error->all("Illegal run_style command");
|
||||
|
||||
int n = strlen(arg[0]) + 1;
|
||||
integrate_style = new char[n];
|
||||
strcpy(integrate_style,arg[0]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void Update::create_minimize(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1) error->all("Illegal min_style command");
|
||||
|
||||
delete [] minimize_style;
|
||||
delete minimize;
|
||||
|
||||
if (0) return; // dummy line to enable else-if macro expansion
|
||||
|
||||
#define MINIMIZE_CLASS
|
||||
#define MinimizeStyle(key,Class) \
|
||||
else if (strcmp(arg[0],#key) == 0) minimize = new Class(lmp);
|
||||
#include "style_minimize.h"
|
||||
#undef MINIMIZE_CLASS
|
||||
|
||||
else error->all("Illegal min_style command");
|
||||
|
||||
int n = strlen(arg[0]) + 1;
|
||||
minimize_style = new char[n];
|
||||
strcpy(minimize_style,arg[0]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
reset timestep from input script
|
||||
do not allow dump files or a restart to be defined
|
||||
do not allow any timestep-dependent fixes to be defined
|
||||
do not allow any dynamic regions to be defined
|
||||
reset eflag/vflag global so nothing will think eng/virial are current
|
||||
reset invoked flags of computes,
|
||||
so nothing will think they are current between runs
|
||||
clear timestep list of computes that store future invocation times
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Update::reset_timestep(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1) error->all("Illegal reset_timestep command");
|
||||
|
||||
for (int i = 0; i < output->ndump; i++)
|
||||
if (output->last_dump[i] >= 0)
|
||||
error->all("Cannot reset timestep with dump file already written to");
|
||||
if (output->restart && output->last_restart >= 0)
|
||||
error->all("Cannot reset timestep with restart file already written");
|
||||
|
||||
for (int i = 0; i < modify->nfix; i++)
|
||||
if (modify->fix[i]->time_depend)
|
||||
error->all("Cannot reset timestep with a time-dependent fix defined");
|
||||
|
||||
for (int i = 0; i < domain->nregion; i++)
|
||||
if (domain->regions[i]->dynamic_check())
|
||||
error->all("Cannot reset timestep with a dynamic region defined");
|
||||
|
||||
eflag_global = vflag_global = -1;
|
||||
|
||||
for (int i = 0; i < modify->ncompute; i++) {
|
||||
modify->compute[i]->invoked_scalar = -1;
|
||||
modify->compute[i]->invoked_vector = -1;
|
||||
modify->compute[i]->invoked_array = -1;
|
||||
modify->compute[i]->invoked_peratom = -1;
|
||||
modify->compute[i]->invoked_local = -1;
|
||||
}
|
||||
|
||||
for (int i = 0; i < modify->ncompute; i++)
|
||||
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
|
||||
|
||||
ntimestep = ATOBIGINT(arg[0]);
|
||||
if (ntimestep < 0) error->all("Timestep must be >= 0");
|
||||
if (ntimestep > MAXBIGINT) error->all("Too big a timestep");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of update and integrate/minimize
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double Update::memory_usage()
|
||||
{
|
||||
double bytes = 0.0;
|
||||
if (whichflag == 1) bytes += integrate->memory_usage();
|
||||
else if (whichflag == 2) bytes += minimize->memory_usage();
|
||||
return bytes;
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue