diff --git a/src/Makefile b/src/Makefile index aa39dff1c6..7b32442c39 100755 --- a/src/Makefile +++ b/src/Makefile @@ -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) diff --git a/src/atom.cpp b/src/atom.cpp index 2ea85dd33b..cb65eabbc2 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -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); diff --git a/src/atom.h b/src/atom.h index 87d75e3b45..b12b337062 100644 --- a/src/atom.h +++ b/src/atom.h @@ -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 diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index 232264e4d9..4447d9f06f 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -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; } diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h index 46124cc30a..1d8cd07f23 100644 --- a/src/compute_property_atom.h +++ b/src/compute_property_atom.h @@ -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); }; } diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 20bbe1f9dc..b385da6577 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -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; + } } diff --git a/src/dump_custom.h b/src/dump_custom.h index 03828adb63..158b52699b 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -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); diff --git a/src/min.cpp b/src/min.cpp index bc4bfc3c79..3539c2e7dc 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -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; } /* ---------------------------------------------------------------------- */ diff --git a/src/min.h b/src/min.h index 53457acc85..64e5453761 100644 --- a/src/min.h +++ b/src/min.h @@ -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 diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index d4322ef06a..f31160e288 100644 --- a/src/min_hftn.cpp +++ b/src/min_hftn.cpp @@ -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); diff --git a/src/min_linesearch.cpp b/src/min_linesearch.cpp index e1de4e350f..92c61841a1 100644 --- a/src/min_linesearch.cpp +++ b/src/min_linesearch.cpp @@ -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); } } diff --git a/src/pair.h b/src/pair.h index 9f1a10cfea..7f7da4927c 100644 --- a/src/pair.h +++ b/src/pair.h @@ -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) {} diff --git a/src/respa.cpp b/src/respa.cpp index f3f0d2b36a..13da9d9e4a 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -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; } } diff --git a/src/respa.h b/src/respa.h index dbb981068a..e0981c764e 100644 --- a/src/respa.h +++ b/src/respa.h @@ -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 diff --git a/src/verlet.cpp b/src/verlet.cpp index 523a03839a..844234171f 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -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; } } } diff --git a/src/verlet.h b/src/verlet.h index 9257691640..02a5c1704e 100644 --- a/src/verlet.h +++ b/src/verlet.h @@ -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(); };