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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4589 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2010-08-25 17:03:55 +00:00
parent 02adb67715
commit a5f54fa1e2
16 changed files with 160 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -16,7 +16,7 @@ OBJ = $(SRC:.cpp=.o)
PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
kspace manybody meam molecule opt peri poems prd reax shock xtc
PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm \
PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-eff \
user-ewaldn user-imd user-smd
PACKALL = $(PACKAGE) $(PACKUSER)

8
src/atom.cpp
View File

@ -74,7 +74,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
x0 = NULL;
spin = NULL;
eradius = evel = eforce = NULL;
eradius = ervel = erforce = NULL;
maxspecial = 1;
nspecial = NULL;
@ -101,7 +101,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
q_flag = mu_flag = 0;
quat_flag = omega_flag = angmom_flag = torque_flag = 0;
radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
spin_flag = eradius_flag = evel_flag = eforce_flag = 0;
spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
// ntype-length arrays
@ -182,8 +182,8 @@ Atom::~Atom()
memory->sfree(spin);
memory->sfree(eradius);
memory->sfree(evel);
memory->sfree(eforce);
memory->sfree(ervel);
memory->sfree(erforce);
memory->sfree(molecule);

4
src/atom.h
View File

@ -53,7 +53,7 @@ class Atom : protected Pointers {
double **x0;
int *spin;
double *eradius,*evel,*eforce;
double *eradius,*ervel,*erforce;
int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
int **special; // IDs of 1-2,1-3,1-4 neighs of each atom
@ -83,7 +83,7 @@ class Atom : protected Pointers {
int q_flag,mu_flag;
int quat_flag,omega_flag,angmom_flag,torque_flag;
int radius_flag,density_flag,rmass_flag,vfrac_flag;
int spin_flag,eradius_flag,evel_flag,eforce_flag;
int spin_flag,eradius_flag,ervel_flag,erforce_flag;
// extra peratom info in restart file destined for fix & diag

24
src/compute_property_atom.cpp
View File

@ -207,16 +207,16 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
error->all("Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_eradius;
} else if (strcmp(arg[iarg],"evel") == 0) {
if (!atom->evel_flag)
} else if (strcmp(arg[iarg],"ervel") == 0) {
if (!atom->ervel_flag)
error->all("Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_evel;
} else if (strcmp(arg[iarg],"eforce") == 0) {
if (!atom->eforce_flag)
pack_choice[i] = &ComputePropertyAtom::pack_ervel;
} else if (strcmp(arg[iarg],"erforce") == 0) {
if (!atom->erforce_flag)
error->all("Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_eforce;
pack_choice[i] = &ComputePropertyAtom::pack_erforce;
} else error->all("Invalid keyword in compute property/atom command");
}
@ -1073,14 +1073,14 @@ void ComputePropertyAtom::pack_eradius(int n)
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_evel(int n)
void ComputePropertyAtom::pack_ervel(int n)
{
double *evel = atom->evel;
double *ervel = atom->ervel;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = evel[i];
if (mask[i] & groupbit) buf[n] = ervel[i];
else buf[n] = 0.0;
n += nvalues;
}
@ -1088,14 +1088,14 @@ void ComputePropertyAtom::pack_evel(int n)
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_eforce(int n)
void ComputePropertyAtom::pack_erforce(int n)
{
double *eforce = atom->eforce;
double *erforce = atom->erforce;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = eforce[i];
if (mask[i] & groupbit) buf[n] = erforce[i];
else buf[n] = 0.0;
n += nvalues;
}

4
src/compute_property_atom.h
View File

@ -92,8 +92,8 @@ class ComputePropertyAtom : public Compute {
void pack_tqz(int);
void pack_spin(int);
void pack_eradius(int);
void pack_evel(int);
void pack_eforce(int);
void pack_ervel(int);
void pack_erforce(int);
};
}

107
src/dump_custom.cpp
View File

@ -38,7 +38,7 @@ enum{ID,MOL,TYPE,MASS,
X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,EVEL,EFORCE,
QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
COMPUTE,FIX,VARIABLE};
enum{LT,LE,GT,GE,EQ,NEQ};
enum{INT,DOUBLE};
@ -693,15 +693,15 @@ int DumpCustom::count()
error->all("Threshhold for an atom property that isn't allocated");
ptr = atom->eradius;
nstride = 1;
} else if (thresh_array[ithresh] == EVEL) {
if (!atom->evel_flag)
} else if (thresh_array[ithresh] == ERVEL) {
if (!atom->ervel_flag)
error->all("Threshhold for an atom property that isn't allocated");
ptr = atom->evel;
ptr = atom->ervel;
nstride = 1;
} else if (thresh_array[ithresh] == EFORCE) {
if (!atom->eforce_flag)
} else if (thresh_array[ithresh] == ERFORCE) {
if (!atom->erforce_flag)
error->all("Threshhold for an atom property that isn't allocated");
ptr = atom->eforce;
ptr = atom->erforce;
nstride = 1;
} else if (thresh_array[ithresh] == COMPUTE) {
@ -1008,15 +1008,15 @@ void DumpCustom::parse_fields(int narg, char **arg)
error->all("Dumping an atom quantity that isn't allocated");
pack_choice[i] = &DumpCustom::pack_eradius;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"evel") == 0) {
if (!atom->evel_flag)
} else if (strcmp(arg[iarg],"ervel") == 0) {
if (!atom->ervel_flag)
error->all("Dumping an atom quantity that isn't allocated");
pack_choice[i] = &DumpCustom::pack_evel;
pack_choice[i] = &DumpCustom::pack_ervel;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"eforce") == 0) {
if (!atom->evel_flag)
} else if (strcmp(arg[iarg],"erforce") == 0) {
if (!atom->erforce_flag)
error->all("Dumping an atom quantity that isn't allocated");
pack_choice[i] = &DumpCustom::pack_eforce;
pack_choice[i] = &DumpCustom::pack_erforce;
vtype[i] = DOUBLE;
// compute value = c_ID
@ -1306,8 +1306,8 @@ int DumpCustom::modify_param(int narg, char **arg)
else if (strcmp(arg[1],"spin") == 0) thresh_array[nthresh] = SPIN;
else if (strcmp(arg[1],"eradius") == 0) thresh_array[nthresh] = ERADIUS;
else if (strcmp(arg[1],"evel") == 0) thresh_array[nthresh] = EVEL;
else if (strcmp(arg[1],"eforce") == 0) thresh_array[nthresh] = EFORCE;
else if (strcmp(arg[1],"ervel") == 0) thresh_array[nthresh] = ERVEL;
else if (strcmp(arg[1],"erforce") == 0) thresh_array[nthresh] = ERFORCE;
// compute value = c_ID
// if no trailing [], then arg is set to 0, else arg is between []
@ -2246,71 +2246,44 @@ void DumpCustom::pack_spin(int n)
}
}
/* ----------------------------------------------------------------------
different interpretation of electron radius
depending on dynamics vs minimization
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_eradius(int n)
{
double *eradius = atom->eradius;
int *spin = atom->spin;
int nlocal = atom->nlocal;
if (update->whichflag == 1) {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = eradius[i];
n += size_one;
}
} else {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
if (spin[i]) buf[n] = exp(eradius[i]);
else buf[n] = 0.0;
n += size_one;
}
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_evel(int n)
{
double *evel = atom->evel;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = evel[i];
buf[n] = eradius[i];
n += size_one;
}
}
/* ----------------------------------------------------------------------
different interpretation of electron radial force
depending on dynamics vs minimization
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_eforce(int n)
void DumpCustom::pack_ervel(int n)
{
double *eradius = atom->eradius;
double *eforce = atom->eforce;
int *spin = atom->spin;
double *ervel = atom->ervel;
int nlocal = atom->nlocal;
if (update->whichflag == 1) {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = eforce[i];
n += size_one;
}
} else {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
if (spin[i]) buf[n] = eforce[i]/exp(eradius[i]);
else buf[n] = 0.0;
n += size_one;
}
}
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = ervel[i];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_erforce(int n)
{
double *erforce = atom->erforce;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = erforce[i];
n += size_one;
}
}

4
src/dump_custom.h
View File

@ -150,8 +150,8 @@ class DumpCustom : public Dump {
void pack_tqz(int);
void pack_spin(int);
void pack_eradius(int);
void pack_evel(int);
void pack_eforce(int);
void pack_ervel(int);
void pack_erforce(int);
void pack_compute(int);
void pack_fix(int);

39
src/min.cpp
View File

@ -135,12 +135,12 @@ void Min::init()
ev_setup();
// set flags for what arrays to clear in force_clear()
// need to clear torques,eforce if arrays exists
// need to clear torques,erforce if arrays exists
torqueflag = 0;
if (atom->torque_flag) torqueflag = 1;
eforceflag = 0;
if (atom->eforce_flag) eforceflag = 1;
erforceflag = 0;
if (atom->erforce_flag) erforceflag = 1;
// orthogonal vs triclinic simulation box
@ -246,6 +246,12 @@ void Min::setup()
if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles
if (nextra_atom)
for (int m = 0; m < nextra_atom; m++)
requestor[m]->min_xf_get(m);
modify->setup(vflag);
output->setup(1);
@ -310,6 +316,12 @@ void Min::setup_minimal(int flag)
if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles
if (nextra_atom)
for (int m = 0; m < nextra_atom; m++)
requestor[m]->min_xf_get(m);
modify->setup(vflag);
// stats for Finish to print
@ -468,6 +480,12 @@ double Min::energy_force(int resetflag)
timer->stamp(TIME_COMM);
}
// update per-atom minimization variables stored by pair styles
if (nextra_atom)
for (int m = 0; m < nextra_atom; m++)
requestor[m]->min_xf_get(m);
// fixes that affect minimization
if (modify->n_min_post_force) modify->min_post_force(vflag);
@ -523,19 +541,21 @@ void Min::force_clear()
}
}
if (eforceflag) {
double *eforce = atom->eforce;
if (erforceflag) {
double *erforce = atom->erforce;
for (i = 0; i < nall; i++)
eforce[i] = 0.0;
erforce[i] = 0.0;
}
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
setup and clear other arrays as needed
pair style makes request to add a per-atom variables to minimization
requestor stores callback to pair class to invoke during min
to get current variable and forces on it and to update the variable
return flag that pair can use if it registers multiple variables
------------------------------------------------------------------------- */
void Min::request(Pair *pair, int peratom, double maxvalue)
int Min::request(Pair *pair, int peratom, double maxvalue)
{
int n = nextra_atom + 1;
xextra_atom = (double **) memory->srealloc(xextra_atom,n*sizeof(double *),
@ -555,6 +575,7 @@ void Min::request(Pair *pair, int peratom, double maxvalue)
extra_peratom[nextra_atom] = peratom;
extra_max[nextra_atom] = maxvalue;
nextra_atom++;
return nextra_atom-1;
}
/* ---------------------------------------------------------------------- */

18
src/min.h
View File

@ -34,7 +34,7 @@ class Min : protected Pointers {
void setup_minimal(int);
void run(int,int);
void cleanup();
void request(class Pair *, int, double);
int request(class Pair *, int, double);
double memory_usage() {return 0.0;}
void modify_params(int, char **);
@ -58,7 +58,7 @@ class Min : protected Pointers {
class Compute **vlist_atom;
int pairflag;
int torqueflag,eforceflag;
int torqueflag,erforceflag;
int triclinic; // 0 if domain is orthog, 1 if triclinic
int narray; // # of arrays stored by fix_minimize
@ -77,13 +77,13 @@ class Min : protected Pointers {
double *fextra; // force vector for extra global dof
// xextra is stored by fix
int nextra_atom; // # of sets of extra per-atom dof
double **xextra_atom; // variables for extra per-atom dof sets
double **fextra_atom; // force vectors for extra per-atom dof sets
int *extra_peratom; // # of per-atom values in each set
int *extra_nlen; // total local length of each set, e.g 3*nlocal
double *extra_max; // max allowed change in one iter for each set
class Pair **requestor; // Pair that requested each set
int nextra_atom; // # of extra per-atom variables
double **xextra_atom; // ptr to the variable
double **fextra_atom; // ptr to the force on the variable
int *extra_peratom; // # of values in variable, e.g. 3 in x
int *extra_nlen; // total local length of variable, e.g 3*nlocal
double *extra_max; // max allowed change per iter for atom's var
class Pair **requestor; // Pair that stores/manipulates the variable
int neigh_every,neigh_delay,neigh_dist_check; // neighboring params

28
src/min_hftn.cpp
View File

@ -180,7 +180,7 @@ void MinHFTN::reset_vectors()
int n = NUM_HFTN_ATOM_BASED_VECTORS;
for (int m = 0; m < nextra_atom; m++) {
extra_nlen[m] = extra_peratom[m] * atom->nlocal;
requestor[m]->min_pointers(&xextra_atom[m], &fextra_atom[m]);
requestor[m]->min_xf_pointers(m,&xextra_atom[m],&fextra_atom[m]);
for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++)
_daExtraAtom[i][m] = fix_minimize->request_vector (n++);
}
@ -323,6 +323,7 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = xkAtom[i];
requestor[m]->min_x_set(m);
}
}
dFinalEnergy = energy_force (0);
@ -366,6 +367,7 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = xkAtom[i] + pAtom[i];
requestor[m]->min_x_set(m);
}
}
if (nextra_global)
@ -485,6 +487,7 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = xkAtom[i];
requestor[m]->min_x_set(m);
}
}
bHaveEvaluatedAtX = false;
@ -651,12 +654,13 @@ bool MinHFTN::compute_inner_cg_step_(const double dTrustRadius,
}
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
double * fatom = fextra_atom[m];
double * rAtom = _daExtraAtom[VEC_CG_R][m];
double * dAtom = _daExtraAtom[VEC_CG_D][m];
int n = extra_nlen[m];
for (int i = 0; i < n; i++) {
rAtom[i] = fextra_atom[m][i];
dAtom[i] = fextra_atom[m][i];
rAtom[i] = fatom[i];
dAtom[i] = fatom[i];
}
}
}
@ -716,6 +720,7 @@ bool MinHFTN::compute_inner_cg_step_(const double dTrustRadius,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = xkAtom[i];
requestor[m]->min_x_set(m);
}
}
dEnergyAtX = energy_force (0);
@ -947,10 +952,11 @@ double MinHFTN::calc_grad_dot_v_using_mpi_(const int nIx) const
dGradDotVLocal += - _daAVectors[nIx][i] * fvec[i];
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
double * fatom = fextra_atom[m];
double * iAtom = _daExtraAtom[nIx][m];
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
dGradDotVLocal += - iAtom[i] * fextra_atom[m][i];
dGradDotVLocal += - iAtom[i] * fatom[i];
}
}
@ -1396,6 +1402,7 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] += dEps * iAtom[i];
requestor[m]->min_x_set(m);
}
}
energy_force (0);
@ -1410,10 +1417,11 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
_daAVectors[VEC_DIF2][i] = fvec[i];
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
double * fatom = fextra_atom[m];
double * d2Atom = _daExtraAtom[VEC_DIF2][m];
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
d2Atom[i] = fextra_atom[m][i];
d2Atom[i] = fatom[i];
}
}
@ -1431,6 +1439,7 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] += d1Atom[i];
requestor[m]->min_x_set(m);
}
}
dNewEnergy = energy_force (0);
@ -1490,6 +1499,7 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] += dEps * iAtom[i];
requestor[m]->min_x_set(m);
}
}
energy_force (0);
@ -1504,10 +1514,11 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
_daAVectors[VEC_DIF2][i] = fvec[i];
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
double * fatom = fextra_atom[m];
double * d2Atom = _daExtraAtom[VEC_DIF2][m];
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
d2Atom[i] = fextra_atom[m][i];
d2Atom[i] = fatom[i];
}
}
@ -1525,6 +1536,7 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = d1Atom[i] - dEps * iAtom[i];
requestor[m]->min_x_set(m);
}
}
energy_force (0);
@ -1544,11 +1556,12 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
(fvec[i] - _daAVectors[VEC_DIF2][i]) / (2.0 * dEps);
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
double * fatom = fextra_atom[m];
double * iAtom = _daExtraAtom[nIxResult][m];
double * d2Atom = _daExtraAtom[VEC_DIF2][m];
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
iAtom[i] = (fextra_atom[m][i] - d2Atom[i]) / (2.0 + dEps);
iAtom[i] = (fatom[i] - d2Atom[i]) / (2.0 + dEps);
}
}
@ -1567,6 +1580,7 @@ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs,
int n = extra_nlen[m];
for (int i = 0; i < n; i++)
xatom[i] = d1Atom[i];
requestor[m]->min_x_set(m);
}
}
dNewEnergy = energy_force (0);

16
src/min_linesearch.cpp
View File

@ -146,7 +146,7 @@ void MinLineSearch::reset_vectors()
int n = 3;
for (int m = 0; m < nextra_atom; m++) {
extra_nlen[m] = extra_peratom[m] * atom->nlocal;
requestor[m]->min_pointers(&xextra_atom[m],&fextra_atom[m]);
requestor[m]->min_xf_pointers(m,&xextra_atom[m],&fextra_atom[m]);
x0extra_atom[m] = fix_minimize->request_vector(n++);
gextra_atom[m] = fix_minimize->request_vector(n++);
hextra_atom[m] = fix_minimize->request_vector(n++);
@ -215,9 +215,9 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha)
alpha = MIN(ALPHA_MAX,dmax/hmaxall);
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
hme = 0.0;
fatom = fextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
hme = 0.0;
for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(alpha,extra_max[m]/hmax);
@ -364,9 +364,9 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
alphamax = MIN(ALPHA_MAX,dmax/hmaxall);
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
hme = 0.0;
fatom = fextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
hme = 0.0;
for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alphamax = MIN(alphamax,extra_max[m]/hmax);
@ -402,7 +402,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
engprev = eoriginal;
alphaprev = 0.0;
// Important diagnostic: test the gradient against energy
// important diagnostic: test the gradient against energy
// double etmp;
// double alphatmp = alphamax*1.0e-4;
// etmp = alpha_step(alphatmp,1);
@ -423,7 +423,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
}
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
xatom = xextra_atom[m];
fatom = fextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) {
@ -521,6 +521,7 @@ double MinLineSearch::alpha_step(double alpha, int resetflag)
x0atom = x0extra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] = x0atom[i];
requestor[m]->min_x_set(m);
}
// step forward along h
@ -534,6 +535,7 @@ double MinLineSearch::alpha_step(double alpha, int resetflag)
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i];
requestor[m]->min_x_set(m);
}
}

10
src/pair.h
View File

@ -71,7 +71,11 @@ class Pair : protected Pointers {
virtual void compute_outer(int, int) {}
virtual double single(int, int, int, int,
double, double, double, double &) {return 0.0;}
double, double
, double, double &) {return 0.0;}
virtual void settings(int, char **) = 0;
virtual void coeff(int, char **) = 0;
@ -96,7 +100,9 @@ class Pair : protected Pointers {
virtual void *extract(char *) {return NULL;}
virtual void swap_eam(double *, double **) {}
virtual void reset_dt() {}
virtual void min_pointers(double **, double **) {}
virtual void min_xf_pointers(int, double **, double **) {}
virtual void min_xf_get(int) {}
virtual void min_x_set(int) {}
virtual int pre_adapt(char *, int, int, int, int) {return -1;}
virtual void adapt(int, int, int, int, int, double) {}

12
src/respa.cpp
View File

@ -278,12 +278,12 @@ void Respa::init()
ev_setup();
// set flags for what arrays to clear in force_clear()
// need to clear torques,eforce if arrays exists
// need to clear torques,erforce if arrays exists
torqueflag = 0;
if (atom->torque_flag) torqueflag = 1;
eforceflag = 0;
if (atom->eforce_flag) eforceflag = 1;
erforceflag = 0;
if (atom->erforce_flag) erforceflag = 1;
// step[] = timestep for each level
@ -614,10 +614,10 @@ void Respa::force_clear(int newtonflag)
}
}
if (eforceflag) {
double *eforce = atom->eforce;
if (erforceflag) {
double *erforce = atom->erforce;
for (i = 0; i < nall; i++)
eforce[i] = 0.0;
erforce[i] = 0.0;
}
}

2
src/respa.h
View File

@ -54,7 +54,7 @@ class Respa : public Integrate {
private:
int triclinic; // 0 if domain is orthog, 1 if triclinic
int torqueflag; // zero out arrays every step
int eforceflag;
int erforceflag;
int *newton; // newton flag at each level
class FixRespa *fix_respa; // Fix to store the force level array

24
src/verlet.cpp
View File

@ -63,12 +63,12 @@ void Verlet::init()
ev_setup();
// set flags for what arrays to clear in force_clear()
// need to clear torques,eforce if arrays exists
// need to clear torques,erforce if arrays exists
torqueflag = 0;
if (atom->torque_flag) torqueflag = 1;
eforceflag = 0;
if (atom->eforce_flag) eforceflag = 1;
erforceflag = 0;
if (atom->erforce_flag) erforceflag = 1;
// orthogonal vs triclinic simulation box
@ -319,10 +319,10 @@ void Verlet::force_clear()
}
}
if (eforceflag) {
double *eforce = atom->eforce;
if (erforceflag) {
double *erforce = atom->erforce;
for (i = 0; i < nall; i++)
eforce[i] = 0.0;
erforce[i] = 0.0;
}
// neighbor includegroup flag is set
@ -348,10 +348,10 @@ void Verlet::force_clear()
}
}
if (eforceflag) {
double *eforce = atom->eforce;
if (erforceflag) {
double *erforce = atom->erforce;
for (i = 0; i < nall; i++)
eforce[i] = 0.0;
erforce[i] = 0.0;
}
if (force->newton) {
@ -372,10 +372,10 @@ void Verlet::force_clear()
}
}
if (eforceflag) {
double *eforce = atom->eforce;
if (erforceflag) {
double *erforce = atom->erforce;
for (i = atom->nlocal; i < nall; i++)
eforce[i] = 0.0;
erforce[i] = 0.0;
}
}
}

2
src/verlet.h
View File

@ -36,7 +36,7 @@ class Verlet : public Integrate {
private:
int triclinic; // 0 if domain is orthog, 1 if triclinic
int torqueflag; // zero out arrays every step
int eforceflag;
int erforceflag;
void force_clear();
};