From b76e645182464a53f287f00c8b0e04567a41d85b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 Dec 2022 16:14:29 -0500 Subject: [PATCH] remove optional code --- src/MC/fix_sgcmc.cpp | 945 +++---------------------------------------- src/MC/fix_sgcmc.h | 73 +--- 2 files changed, 58 insertions(+), 960 deletions(-) diff --git a/src/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index 07b9b7310e..280c0b6f53 100644 --- a/src/MC/fix_sgcmc.cpp +++ b/src/MC/fix_sgcmc.cpp @@ -56,15 +56,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -// This is for debugging purposes. The ASSERT() macro is used in the code to check -// if everything runs as expected. -#if SGCMC_DEBUG - inline void my_noop() {} - #define ASSERT(cond) ((!(cond)) ? error->one(FLERR, "Assertion failure.") : my_noop()) -#else - #define ASSERT(cond) -#endif - /********************************************************************* * Constructs the fix object and parses the input parameters * that control the Monte Carlo routine. @@ -90,17 +81,6 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *lmp, int narg, char **a this->comm_forward = 4; this->comm_reverse = 3; -#if SGCMC_DEBUG - this->peratom_flag = 1; - this->size_peratom_cols = 0; - this->peratom_freq = 1; - this->create_attribute = 1; - trialCounters = nullptr; - grow_arrays(atom->nmax); - ASSERT(trialCounters != nullptr || atom->nlocal == 0); - memset(trialCounters, 0, sizeof(trialCounters[0]) * atom->nlocal); -#endif - if(domain->triclinic) error->all(FLERR, "Fix sgcmc does not support non-orthogonal simulation boxes."); @@ -245,58 +225,8 @@ FixSemiGrandCanonicalMC::~FixSemiGrandCanonicalMC() delete random; delete localRandom; -#if SGCMC_DEBUG - memory->sfree(trialCounters); -#endif } -#if SGCMC_DEBUG - -/********************************************************************* - * Allocate atom-based array. - *********************************************************************/ -void FixSemiGrandCanonicalMC::grow_arrays(int nmax) -{ - trialCounters = (double*)memory->srealloc(trialCounters, atom->nmax * sizeof(trialCounters[0]), "sgcmc:trialcounters"); - vector_atom = trialCounters; -} - -/********************************************************************* - * Copy values within local atom-based array. - *********************************************************************/ -void FixSemiGrandCanonicalMC::copy_arrays(int i, int j) -{ - trialCounters[j] = trialCounters[i]; -} - -/********************************************************************* - * Initialize one atom's array values, called when atom is created. - *********************************************************************/ -void FixSemiGrandCanonicalMC::set_arrays(int i) -{ - trialCounters[i] = 0; -} - -/********************************************************************* - * Pack values in local atom-based array for exchange with another proc. - *********************************************************************/ -int FixSemiGrandCanonicalMC::pack_exchange(int i, double *buf) -{ - buf[0] = trialCounters[i]; - return 1; -} - -/********************************************************************* - * Unpack values in local atom-based array from exchange with another proc. - *********************************************************************/ -int FixSemiGrandCanonicalMC::unpack_exchange(int nlocal, double *buf) -{ - trialCounters[nlocal] = buf[0]; - return 1; -} - -#endif - /********************************************************************* * The return value of this method specifies at which points the * fix is invoked during the simulation. @@ -325,18 +255,8 @@ void FixSemiGrandCanonicalMC::init() if(count > 1) error->all(FLERR, "More than one fix sgcmc defined."); // Save a pointer to the EAM potential. - if((pairEAM = dynamic_cast(force->pair))) { -#if CDEAM_MC_SUPPORT - if((pairCDEAM = dynamic_cast(pairEAM))) { - if(pairCDEAM->cdeamVersion != 1) - error->all(FLERR, "The sgcmc fix works only with the one-site concentration version of the CD-EAM potential."); - } -#endif - } -#if TERSOFF_MC_SUPPORT - else if((pairTersoff = dynamic_cast(force->pair))) {} -#endif - else { + pairEAM = dynamic_cast(force->pair); + if (!pairEAM) { if (comm->me == 0) utils::logmesg(lmp, " SGC - Using naive total energy calculation for MC -> SLOW!\n"); @@ -357,14 +277,14 @@ void FixSemiGrandCanonicalMC::init() const int *mask = atom->mask; std::vector localSpeciesCounts(atom->ntypes+1, 0); for(int i = 0; i < atom->nlocal; i++, ++type) { - ASSERT(*type >= 1 && *type <= atom->ntypes); if(mask[i] & groupbit) localSpeciesCounts[*type]++; } // MPI sum to get global concentrations. speciesCounts.resize(atom->ntypes+1); - MPI_Allreduce(&localSpeciesCounts.front(), &speciesCounts.front(), localSpeciesCounts.size(), MPI_INT, MPI_SUM, world); + MPI_Allreduce(&localSpeciesCounts.front(), &speciesCounts.front(), localSpeciesCounts.size(), + MPI_INT, MPI_SUM, world); } /********************************************************************* @@ -406,23 +326,6 @@ void FixSemiGrandCanonicalMC::doMC() fetchGhostAtomElectronDensities(); const int *mask = atom->mask; -#if SGCMC_DEBUG - // This check is for debugging only! Can be safely removed. - // Check the global concentration counter if it is still in sync on all nodes. - const int *type = atom->type; - std::vector localSpeciesCounts(atom->ntypes+1, 0); - std::vector globalSpeciesCounts(atom->ntypes+1, 0); - for(int i = 0; i < atom->nlocal; i++, ++type) { - if(mask[i] & groupbit) - localSpeciesCounts[*type]++; - } - MPI_Allreduce(&localSpeciesCounts.front(), &globalSpeciesCounts.front(), localSpeciesCounts.size(), MPI_INT, MPI_SUM, world); - for(int i = 1; i <= atom->ntypes; i++) { - // Note: If this test fails it might be an error in the code or - // because LAMMPS has lost atoms. This can happen for bad atom systems. - ASSERT(globalSpeciesCounts[i] == speciesCounts[i]); - } -#endif // End of debugging code // Reset counters. int nAcceptedSwapsLocal = 0; @@ -474,20 +377,12 @@ void FixSemiGrandCanonicalMC::doMC() // Choose a random atom from the pool of atoms that are inside the sampling window. int index = (int)(localRandom->uniform() * (double)samplingWindowAtoms.size()); - ASSERT(index < samplingWindowAtoms.size()); selectedAtomNL = samplingWindowAtoms[index]; // Get the real atom index. - ASSERT(selectedAtomNL < neighborList->inum); selectedAtom = neighborList->ilist[selectedAtomNL]; - ASSERT(selectedAtom < atom->nlocal); - ASSERT(selectedAtom == selectedAtomNL); // This assumption may be wrong. oldSpecies = atom->type[selectedAtom]; -#if SGCMC_DEBUG - trialCounters[selectedAtom]++; -#endif - // Choose the new type for the swapping atom by random. if(atom->ntypes > 2) { // Use a random number to choose the new species if there are three or more atom types. @@ -498,25 +393,13 @@ void FixSemiGrandCanonicalMC::doMC() // If there are only two atom types, then the decision is clear. newSpecies = (oldSpecies == 1) ? 2 : 1; } - ASSERT(newSpecies >= 1 && newSpecies <= atom->ntypes && newSpecies != oldSpecies); deltaN[oldSpecies] = -1; deltaN[newSpecies] = +1; // Compute the energy difference that swapping this atom would cost or gain. if(pairEAM) { -#if CDEAM_MC_SUPPORT - if(pairCDEAM) // Concentration dependent EAM case: - deltaE = computeEnergyChangeCDEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - else // Standard EAM case: -#endif - deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - } -#if TERSOFF_MC_SUPPORT - else if(pairTersoff) { - deltaE = computeEnergyChangeTersoff(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - } -#endif - else { + deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); + } else { // Generic case: deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies); } @@ -537,46 +420,15 @@ void FixSemiGrandCanonicalMC::doMC() dm += deltamu[i] * deltaN[i]; } double deltaB = -(deltaE + dm) * beta; -#if SGCMC_DEBUG - std::cout << "Old species: " << oldSpecies << " new species: " << newSpecies << " deltaE: " << deltaE << " dm: " << dm << " deltaB: " << deltaB; -#endif if(deltaB < 0.0) { if(deltaB < log(localRandom->uniform())) { std::fill(deltaN.begin(), deltaN.end(), 0); selectedAtom = -1; deltaE = 0; -#if SGCMC_DEBUG - std::cout << " REJECTED"; -#endif } } - -#if SGCMC_DEBUG - cout << endl; -#endif } -#if 0 - // This is for debugging purposes only to check the energy change calculation routine. - // Check the return value by calculating the energy change the slow way: difference between total energy before and after the swap. - // The following code should be deactivated for full performance. - // WARNING: Calling this function screws up the stored forces and will therefore change the time evolution of the system in - // the MD simulation. - double deltaECheck = computeEnergyChangeGeneric(selectedAtomDebug, oldSpecies, newSpecies); - double deltaETotal, deltaECheckTotal; - MPI_Allreduce(&deltaEDebug, &deltaETotal, 1, MPI_DOUBLE, MPI_SUM, world); - MPI_Allreduce(&deltaECheck, &deltaECheckTotal, 1, MPI_DOUBLE, MPI_SUM, world); - double posx = selectedAtomDebug >= 0 ? atom->x[selectedAtomDebug][0] : 0.0; - double posy = selectedAtomDebug >= 0 ? atom->x[selectedAtomDebug][1] : 0.0; - double posz = selectedAtomDebug >= 0 ? atom->x[selectedAtomDebug][2] : 0.0; - printLog("Checking MC routine. DeltaE=%f DeltaE(check)=%f Atom: %i [%f %f %f] Old species: %i New species: %i\n", - deltaETotal, deltaECheckTotal, selectedAtomDebug, posx, posy, posz, oldSpecies, newSpecies); - if(fabs(deltaECheckTotal - deltaETotal) > 1e-6) { // Delta E must be equal to the total energy difference. - error->one(FLERR, "Error in MC energy routine detected. Computed energy change deviates from correct value."); - } -#endif // End of debugging code - - if(kappa != 0.0 && serialMode == false) { // What follows is the second rejection test for the variance-constrained @@ -613,14 +465,8 @@ void FixSemiGrandCanonicalMC::doMC() // Make accepted atom swap permanent. if(selectedAtom >= 0) { - if(pairEAM) { -#if CDEAM_MC_SUPPORT - if(pairCDEAM) // Concentration dependent EAM case: - flipAtomCDEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - else // Standard EAM case: -#endif - flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - } + if(pairEAM) + flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); else flipAtomGeneric(selectedAtom, oldSpecies, newSpecies); nAcceptedSwapsLocal++; @@ -656,7 +502,6 @@ void FixSemiGrandCanonicalMC::doMC() const int *type = atom->type; std::vector localSpeciesCounts(atom->ntypes+1, 0); for(int i = 0; i < atom->nlocal; i++, ++type) { - ASSERT(*type >= 1 && *type <= atom->ntypes); if(mask[i] & groupbit) localSpeciesCounts[*type]++; } @@ -707,56 +552,20 @@ int FixSemiGrandCanonicalMC::pack_forward_comm(int n, int* list, double* buf, in int m = 0; if(communicationStage == 1) { // Send electron densities of local atoms to neighbors. - ASSERT(pairEAM->rho != nullptr); -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = 0; i < n; i++) - buf[m++] = pairEAM->rho[list[i]]; -#if CDEAM_MC_SUPPORT + for(int i = 0; i < n; i++) buf[m++] = pairEAM->rho[list[i]]; + } else if(communicationStage == 3) { + if(pairEAM) { + // Send types and rhos of real atoms to the ghost atoms of the neighbor proc. + for(int i = 0; i < n; i++) { + buf[m++] = atom->type[list[i]]; + buf[m++] = pairEAM->rho[list[i]]; } - else { - // In case of the CD-EAM model we have to send the RhoB values and D values as well. - for(int i = 0; i < n; i++) { - buf[m++] = pairCDEAM->rho[list[i]]; - buf[m++] = pairCDEAM->rhoB[list[i]]; - buf[m++] = pairCDEAM->D_values[list[i]]; - } + } else { + // Generic potential case: + for(int i = 0; i < n; i++) { + buf[m++] = atom->type[list[i]]; } -#endif - } - else if(communicationStage == 3) { - if(pairEAM) { - // Send types and rhos of real atoms to the ghost atoms of the neighbor proc. -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = 0; i < n; i++) { - buf[m++] = atom->type[list[i]]; - buf[m++] = pairEAM->rho[list[i]]; - } -#if CDEAM_MC_SUPPORT - } - else { - // In case of the CD-EAM model we have to send the RhoB values and D values as well. - for(int i = 0; i < n; i++) { - buf[m++] = atom->type[list[i]]; - buf[m++] = pairCDEAM->rho[list[i]]; - buf[m++] = pairCDEAM->rhoB[list[i]]; - buf[m++] = pairCDEAM->D_values[list[i]]; - } - } -#endif - } - else { - // Generic potential case: - for(int i = 0; i < n; i++) { - buf[m++] = atom->type[list[i]]; - } - } - } - else { - ASSERT(false); + } } return m; } @@ -766,70 +575,29 @@ int FixSemiGrandCanonicalMC::pack_forward_comm(int n, int* list, double* buf, in *********************************************************************/ void FixSemiGrandCanonicalMC::unpack_forward_comm(int n, int first, double* buf) { - if(communicationStage == 1) { - // Receive electron densities of ghost atoms from neighbors. - int last = first + n; -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = first; i < last; i++) - pairEAM->rho[i] = *buf++; -#if CDEAM_MC_SUPPORT - } - else { - // Also receive the partial densities and D values when using the CD-EAM model. - for(int i = first; i < last; i++) { - pairCDEAM->rho[i] = *buf++; - pairCDEAM->rhoB[i] = *buf++; - pairCDEAM->D_values[i] = *buf++; - } - } -#endif - } - else if(communicationStage == 3) { - int last = first + n; - if(pairEAM) { - // Receive types and rhos of real atoms of the neighbor proc and assign them - // to the local ghost atoms. -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = first; i < last; i++, buf += 2) { - atom->type[i] = (int)buf[0]; - // We have to make sure that rhos changed locally do not get overridden by the rhos - // sent by the neighbor procs. - ASSERT(i < (int)changedAtoms.size()); - if(!changedAtoms[i]) - pairEAM->rho[i] = buf[1]; - } -#if CDEAM_MC_SUPPORT - } - else { - // Also receive the partial densities and D values when using the CD-EAM model. - for(int i = first; i < last; i++, buf += 4) { - atom->type[i] = (int)buf[0]; - // We have to make sure that values changed locally do not get overridden by the values - // sent by the neighbor procs. - ASSERT(i < (int)changedAtoms.size()); - if(!changedAtoms[i]) { - pairCDEAM->rho[i] = buf[1]; - pairCDEAM->rhoB[i] = buf[2]; - pairCDEAM->D_values[i] = buf[3]; - } - } - } -#endif - } - else { - // Generic potential case: - for(int i = first; i < last; i++, buf += 1) { - atom->type[i] = (int)buf[0]; - } - } - } - else { - ASSERT(false); + if(communicationStage == 1) { + // Receive electron densities of ghost atoms from neighbors. + int last = first + n; + for(int i = first; i < last; i++) pairEAM->rho[i] = *buf++; + } else if(communicationStage == 3) { + int last = first + n; + if(pairEAM) { + // Receive types and rhos of real atoms of the neighbor proc and assign them + // to the local ghost atoms. + for(int i = first; i < last; i++, buf += 2) { + atom->type[i] = (int)buf[0]; + // We have to make sure that rhos changed locally do not get overridden by the rhos + // sent by the neighbor procs. + if(!changedAtoms[i]) + pairEAM->rho[i] = buf[1]; + } + } else { + // Generic potential case: + for(int i = first; i < last; i++, buf += 1) { + atom->type[i] = (int)buf[0]; + } } + } } /********************************************************************* @@ -837,29 +605,12 @@ void FixSemiGrandCanonicalMC::unpack_forward_comm(int n, int first, double* buf) *********************************************************************/ int FixSemiGrandCanonicalMC::pack_reverse_comm(int n, int first, double* buf) { - ASSERT(communicationStage == 2); - int m = 0; + int m = 0; - // Send changed electron densities of ghost atoms to the real atoms of neighbor procs. - ASSERT(pairEAM->rho != nullptr); - int last = first + n; -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = first; i < last; i++) - buf[m++] = pairEAM->rho[i]; -#if CDEAM_MC_SUPPORT - } - else { - // In case of the CD-EAM model we have to send the RhoB values and D values as well. - for(int i = first; i < last; i++) { - buf[m++] = pairCDEAM->rho[i]; - buf[m++] = pairCDEAM->rhoB[i]; - buf[m++] = pairCDEAM->D_values[i]; - } - } -#endif - return m; + // Send changed electron densities of ghost atoms to the real atoms of neighbor procs. + int last = first + n; + for(int i = first; i < last; i++) buf[m++] = pairEAM->rho[i]; + return m; } /********************************************************************* @@ -867,36 +618,14 @@ int FixSemiGrandCanonicalMC::pack_reverse_comm(int n, int first, double* buf) *********************************************************************/ void FixSemiGrandCanonicalMC::unpack_reverse_comm(int n, int *list, double* buf) { - ASSERT(communicationStage == 2); - - // Received changed electron densities of ghost atoms of neighbor procs and assign them to our - // real atoms. - ASSERT(pairEAM->rho != nullptr); -#if CDEAM_MC_SUPPORT - if(pairCDEAM == nullptr) { -#endif - for(int i = 0; i < n; i++, buf++) { - // We have to make sure that rhos changed locally do not get overridden by the rhos - // sent by the neighbor procs. - ASSERT(list[i] < (int)changedAtoms.size()); - if(!changedAtoms[list[i]]) - pairEAM->rho[list[i]] = *buf; - } -#if CDEAM_MC_SUPPORT - } - else { - for(int i = 0; i < n; i++, buf += 3) { - // We have to make sure that rhos changed locally do not get overridden by the rhos - // sent by the neighbor procs. - ASSERT(list[i] < (int)changedAtoms.size()); - if(!changedAtoms[list[i]]) { - pairCDEAM->rho[list[i]] = buf[0]; - pairCDEAM->rhoB[list[i]] = buf[1]; - pairCDEAM->D_values[list[i]] = buf[2]; - } - } - } -#endif + // Received changed electron densities of ghost atoms of neighbor procs and assign them to our + // real atoms. + for(int i = 0; i < n; i++, buf++) { + // We have to make sure that rhos changed locally do not get overridden by the rhos + // sent by the neighbor procs. + if(!changedAtoms[list[i]]) + pairEAM->rho[list[i]] = *buf; + } } /********************************************************************* @@ -904,9 +633,6 @@ void FixSemiGrandCanonicalMC::unpack_reverse_comm(int n, int *list, double* buf) *********************************************************************/ bool FixSemiGrandCanonicalMC::placeSamplingWindow() { - ASSERT(neighborList != nullptr); - ASSERT(neighborList->inum == atom->nlocal); - // By default the size of the sampling window is the size of the processor bounds minus two cutoff radii. // This ensures that changing atoms in the sampling windows of two adjacent processors cannot affect // the same atoms in the region between the two sampling windows. @@ -949,7 +675,6 @@ bool FixSemiGrandCanonicalMC::placeSamplingWindow() numSamplingWindowAtoms = 0; numFixAtomsLocal = 0; - ASSERT(atom->nlocal == neighborList->inum); const int *mask = atom->mask; for(int ii = 0; ii < neighborList->inum; ii++) { int i = neighborList->ilist[ii]; @@ -1023,7 +748,6 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEAM(int flipAtom, int flipAto double zi = atom->x[flipAtom][2]; // Loop over all neighbors of the selected atom. - ASSERT(flipAtomNL < neighborList->inum); int* jlist = neighborList->firstneigh[flipAtomNL]; int jnum = neighborList->numneigh[flipAtomNL]; for(int jj = 0; jj < jnum; jj++) { @@ -1105,414 +829,6 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEAM(int flipAtom, int flipAto return deltaE; } -#if CDEAM_MC_SUPPORT - -/********************************************************************* - * Calculates the change in energy that swapping the given - * atom would produce. - * This routine is for the concentration-dependent EAM potential. - * It has to do some extra work in comparison to the standard EAM routine above. - * - * Parameters: - * - * flipAtom [in] - * This specifies the atom to be swapped. It's an index into the local list of atoms. - * - * flipAtomNL [in] - * This specifies the atom to be swapped. It's an index into the neighbor list. - * - * oldSpecies [in] - * The current species of the atom before the routine is called. - * - * newSpecies [in] - * The new species of the atom. The atom's type is not changed by this routine. It only computes the induced energy change. - * - * Return value: - * The expected change in total potential energy. - *********************************************************************/ -double FixSemiGrandCanonicalMC::computeEnergyChangeCDEAM(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies) -{ - ASSERT(pairCDEAM != nullptr); // Make sure we have a CD-EAM potential in use. - - double p; - int m; - double* rho = pairEAM->rho; - double* coeff; - double new_total_rho_i = 0.0; - double new_total_rhoB_i = 0.0; - - // The energy change to calculate. - double deltaE = 0.0; - - // Calculate each change of electron density (and partial density) at the - // surrounding sites induced by the swapped atom. Also calculate the change of pair interaction energy. - // Then calculate the change of embedding energy for each neighbor atom. - - double xi = atom->x[flipAtom][0]; - double yi = atom->x[flipAtom][1]; - double zi = atom->x[flipAtom][2]; - - /// The change in atom i's D value. - double deltaD_i = 0.0; - - // Loop over all neighbors of the selected atom. - int* jlist = neighborList->firstneigh[flipAtomNL]; - int jnum = neighborList->numneigh[flipAtomNL]; - for(int jj = 0; jj < jnum; jj++) { - int j = jlist[jj]; - - double delx = xi - atom->x[j][0]; - double dely = yi - atom->x[j][1]; - double delz = zi - atom->x[j][2]; - double rsq = delx*delx + dely*dely + delz*delz; - if(rsq >= pairEAM->cutforcesq) continue; - - int jtype = atom->type[j]; - double r = sqrt(rsq); - - p = r * pairEAM->rdr + 1.0; - m = static_cast(p); - m = MIN(m, pairEAM->nr - 1); - p -= m; - p = MIN(p, 1.0); - - // Calculate change of electron density at site j. - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[oldSpecies][jtype]][m]; - double oldrho_contr = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[newSpecies][jtype]][m]; - double newrho_contr = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - double delta_rho = newrho_contr - oldrho_contr; - - // Sum total rho at site of swapped atom. - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[jtype][newSpecies]][m]; - double rho_ji = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - new_total_rho_i += rho_ji; - if(jtype == pairCDEAM->speciesB) - new_total_rhoB_i += rho_ji; - - // Determine change of partial electron density at site j. - // It increases if the atom i becomes a B atom and it decreases if it was an B atom. - double delta_rhoB_j = 0; - if(newSpecies == pairCDEAM->speciesB) delta_rhoB_j = newrho_contr; - else if(oldSpecies == pairCDEAM->speciesB) delta_rhoB_j = -oldrho_contr; - - // Now we can calculate the new concentration x_j at site j. - double new_x_j = (pairCDEAM->rhoB[j] + delta_rhoB_j) / (rho[j] + delta_rho); - // Calculate the old concentration x_j at site j as well. - double old_x_j = pairCDEAM->rhoB[j] / rho[j]; - - // Calculate change of pair energy ij. - coeff = pairEAM->z2r_spline[pairEAM->type2z2r[oldSpecies][jtype]][m]; - double oldz2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - coeff = pairEAM->z2r_spline[pairEAM->type2z2r[newSpecies][jtype]][m]; - double newz2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - /// Was the old i-j interaction concentration dependent? - if((jtype == pairCDEAM->speciesA && oldSpecies == pairCDEAM->speciesB) - || (jtype == pairCDEAM->speciesB && oldSpecies == pairCDEAM->speciesA)) { - - // The old pair interaction was concentration dependent. - // Please note that it must now become concentration independent that - // either the A or the B atom has changed to another species. - // We here require that only the A-B interactions are concentration dependent. - - // Add the new static pair interaction term. - deltaE += newz2 / r; - - // The D values of i and j decrease due to the removed concentration dependent cross interaction. - double deltaD_j = -oldz2 / r; - deltaD_i += deltaD_j; - - // Since the concentration x_j at site j changes, the mixed pair interaction between atom j and - // other atoms k is also affected. - // The energy change of site j due to the cross term is: (new_h*new_D - old_h*old_D)/2 - double old_h_j = pairCDEAM->evalH(old_x_j); - double new_h_j = pairCDEAM->evalH(new_x_j); - deltaE += 0.5 * (new_h_j * (pairCDEAM->D_values[j] + deltaD_j) - old_h_j * pairCDEAM->D_values[j]); - } - else { - - // The old pair interaction was not concentration dependent. Now it might have - // become dependent. Do check: - if((jtype == pairCDEAM->speciesA && newSpecies == pairCDEAM->speciesB) - || (jtype == pairCDEAM->speciesB && newSpecies == pairCDEAM->speciesA)) { - - // The new pair interaction is concentration dependent. It's an AB interaction. - - // Subtract old static pair interaction. - deltaE -= oldz2 / r; - - // The D values of i and j increase due to the created AB cross interaction. - double deltaD_j = newz2 / r; - deltaD_i += deltaD_j; - - // Since the concentration x_j at site j changes, the mixed pair interaction between atom j and - // other atoms k is also affected. - // The energy change of site j due to the cross term is: (new_h*new_D - old_h*old_D)/2 - double old_h_j = pairCDEAM->evalH(old_x_j); - double new_h_j = pairCDEAM->evalH(new_x_j); - deltaE += 0.5 * (new_h_j * (pairCDEAM->D_values[j] + deltaD_j) - old_h_j * pairCDEAM->D_values[j]); - } - else { - - // The pair interaction stays concentration independent. - // This is like standard EAM. - deltaE += (newz2 - oldz2) / r; - - double old_h_j = pairCDEAM->evalH(old_x_j); - double new_h_j = pairCDEAM->evalH(new_x_j); - deltaE += 0.5 * (new_h_j * (pairCDEAM->D_values[j]) - old_h_j * pairCDEAM->D_values[j]); - } - } - - // Calculate old embedding energy of atom j. - p = rho[j] * pairEAM->rdrho + 1.0; - m = static_cast(p); - m = MAX(1, MIN(m, pairEAM->nrho - 1)); - p -= m; - p = MIN(p, 1.0); - coeff = pairEAM->frho_spline[pairEAM->type2frho[jtype]][m]; - double oldF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - // Calculate new embedding energy of atom j. - p = (rho[j] + delta_rho) * pairEAM->rdrho + 1.0; - m = static_cast(p); - m = MAX(1, MIN(m, pairEAM->nrho - 1)); - p -= m; - p = MIN(p, 1.0); - coeff = pairEAM->frho_spline[pairEAM->type2frho[jtype]][m]; - double newF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - deltaE += newF - oldF; - } - - ASSERT(rho[flipAtom] > 0.0); - ASSERT(new_total_rho_i > 0.0); - - // Things are easier if the rho(r) functional does not depend on the type of both atoms I and J (as for Finnis/Sinclair type potentials). - if(rho[flipAtom] == new_total_rho_i && pairCDEAM->rhoB[flipAtom] == new_total_rhoB_i) { - // Calculate local concentration at site i. This did not change. - double x_i = pairCDEAM->rhoB[flipAtom] / rho[flipAtom]; - - // Calculate h(x_i) polynomial function. - double h_i = pairCDEAM->evalH(x_i); - - // This is the energy change at site i due to the cross term: - deltaE += 0.5 * h_i * deltaD_i; - - // Compute the change in embedding energy of the swapping atom. - p = rho[flipAtom] * pairEAM->rdrho + 1.0; - m = static_cast(p); - m = MAX(1, MIN(m, pairEAM->nrho - 1)); - p -= m; - p = MIN(p, 1.0); - coeff = pairEAM->frho_spline[pairEAM->type2frho[oldSpecies]][m]; - double oldF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - coeff = pairEAM->frho_spline[pairEAM->type2frho[newSpecies]][m]; - double newF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - deltaE += newF - oldF; - } - else { - // Calculate the new and old concentration at site i. - double x_i_old = pairCDEAM->rhoB[flipAtom] / rho[flipAtom]; - double x_i_new = new_total_rhoB_i / new_total_rho_i; - - // Calculate h(x_i) polynomial function. - double old_h_i = pairCDEAM->evalH(x_i_old); - double new_h_i = pairCDEAM->evalH(x_i_new); - - // This is the energy change at site i due to the cross term: - deltaE += 0.5 * (new_h_i * (pairCDEAM->D_values[flipAtom] + deltaD_i) - old_h_i * pairCDEAM->D_values[flipAtom]); - - // Compute the change in embedding energy of the swapping atom. - p = rho[flipAtom] * pairEAM->rdrho + 1.0; - m = static_cast(p); - m = MAX(1, MIN(m, pairEAM->nrho - 1)); - p -= m; - p = MIN(p, 1.0); - coeff = pairEAM->frho_spline[pairEAM->type2frho[oldSpecies]][m]; - double oldF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - p = new_total_rho_i * pairEAM->rdrho + 1.0; - m = static_cast(p); - m = MAX(1, MIN(m, pairEAM->nrho - 1)); - p -= m; - p = MIN(p, 1.0); - coeff = pairEAM->frho_spline[pairEAM->type2frho[newSpecies]][m]; - double newF = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - deltaE += newF - oldF; - } - - return deltaE; -} - -#endif - -#if TERSOFF_MC_SUPPORT - -/********************************************************************* - * Calculates the change in energy that swapping the given - * atom would produce. This routine is for the Tersoff potential. - * - * Parameters: - * - * flipAtom [in] - * This specifies the atom to be swapped. It's an index into the local list of atoms. - * - * flipAtomNL [in] - * This specifies the atom to be swapped. It's an index into the neighbor list. - * - * oldSpecies [in] - * The current species of the atom before the routine is called. - * - * newSpecies [in] - * The new species of the atom. The atom's type is not changed by this routine. It only computes the induced energy change. - * - * Return value: - * The expected change in total potential energy. - *********************************************************************/ -double FixSemiGrandCanonicalMC::computeEnergyChangeTersoff(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies) -{ - // This routine is called even when no trial move is being performed during the - // the current iteration to keep the parallel processors in sync. If no trial - // move is performed then the energy is calculated twice for the same state of the system. - if(flipAtom >= 0) { - // Change system. Perform trial move. - atom->type[flipAtom] = newSpecies; - } - // Transfer changed atom types of the real atoms to the ghost atoms. - communicationStage = 3; - comm->forward_comm(this); - - // Calculate new total energy. - double newEnergy = 0; - if(flipAtom >= 0) { - newEnergy = computeEnergyTersoff(flipAtom); - } - - // Undo trial move. Restore old system state. - if(flipAtom >= 0) { - atom->type[flipAtom] = oldSpecies; - } - // Transfer changed atom types of the real atoms to the ghost atoms. - communicationStage = 3; - comm->forward_comm(this); - - // Calculate old total energy. - double oldEnergy = 0; - if(flipAtom >= 0 || totalPotentialEnergy == 0) { - totalPotentialEnergy = oldEnergy = computeEnergyTersoff(flipAtom); - } - else oldEnergy = totalPotentialEnergy; - - return newEnergy - oldEnergy; -} - -/// Computes the energy of the atom group around the flipped atom using the Tersoff potential. -double FixSemiGrandCanonicalMC::computeEnergyTersoff(int flipAtom) -{ - double totalEnergy = 0; - for(int ii = 0; ii < atom->nlocal; ii++) - totalEnergy += computeAtomicEnergyTersoff(ii); - return totalEnergy; -} - -/// Computes the energy of an atom using the Tersoff potential. -double FixSemiGrandCanonicalMC::computeAtomicEnergyTersoff(int i) -{ - double **x = atom->x; - int *tag = atom->tag; - int *type = atom->type; - - int* numneigh = neighborList->numneigh; - int** firstneigh = neighborList->firstneigh; - - double atomicEnergy = 0; - - int itype = pairTersoff->map[type[i]]; - int itag = tag[i]; - double xtmp = x[i][0]; - double ytmp = x[i][1]; - double ztmp = x[i][2]; - - // two-body interactions, skip half of them - int* jlist = firstneigh[i]; - int jnum = numneigh[i]; - for(int jj = 0; jj < jnum; jj++) { - int j = jlist[jj]; - j &= NEIGHMASK; - int jtype = pairTersoff->map[type[j]]; - int jtag = tag[j]; - - if (itag > jtag) { - if ((itag+jtag) % 2 == 0) continue; - } else if (itag < jtag) { - if ((itag+jtag) % 2 == 1) continue; - } else { - if (x[j][2] < x[i][2]) continue; - if (x[j][2] == ztmp && x[j][1] < ytmp) continue; - if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; - } - - double delx = xtmp - x[j][0]; - double dely = ytmp - x[j][1]; - double delz = ztmp - x[j][2]; - double rsq = delx*delx + dely*dely + delz*delz; - - int iparam_ij = pairTersoff->elem2param[itype][jtype][jtype]; - if (rsq > pairTersoff->params[iparam_ij].cutsq) continue; - - double r = sqrt(rsq); - double tmp_fc = pairTersoff->ters_fc(r, &pairTersoff->params[iparam_ij]); - double tmp_exp = exp(-pairTersoff->params[iparam_ij].lam1 * r); - atomicEnergy += tmp_fc * pairTersoff->params[iparam_ij].biga * tmp_exp; - } - - // three-body interactions - // skip immediately if I-J is not within cutoff - for(int jj = 0; jj < jnum; jj++) { - int j = jlist[jj]; - j &= NEIGHMASK; - int jtype = pairTersoff->map[type[j]]; - int iparam_ij = pairTersoff->elem2param[itype][jtype][jtype]; - - double delr1[3]; - delr1[0] = x[j][0] - xtmp; - delr1[1] = x[j][1] - ytmp; - delr1[2] = x[j][2] - ztmp; - double rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - if(rsq1 > pairTersoff->params[iparam_ij].cutsq) continue; - - // accumulate bondorder zeta for each i-j interaction via loop over k - double zeta_ij = 0.0; - for (int kk = 0; kk < jnum; kk++) { - if (jj == kk) continue; - int k = jlist[kk]; - k &= NEIGHMASK; - int ktype = pairTersoff->map[type[k]]; - int iparam_ijk = pairTersoff->elem2param[itype][jtype][ktype]; - double delr2[3]; - delr2[0] = x[k][0] - xtmp; - delr2[1] = x[k][1] - ytmp; - delr2[2] = x[k][2] - ztmp; - double rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - if(rsq2 > pairTersoff->params[iparam_ijk].cutsq) continue; - zeta_ij += pairTersoff->zeta(&pairTersoff->params[iparam_ijk],rsq1,rsq2,delr1,delr2); - } - - double r = sqrt(rsq1); - double fa = pairTersoff->ters_fa(r,&pairTersoff->params[iparam_ij]); - double bij = pairTersoff->ters_bij(zeta_ij,&pairTersoff->params[iparam_ij]); - atomicEnergy += 0.5*bij*fa; - } - - return atomicEnergy; -} - -#endif - /********************************************************************* * Calculates the change in energy that swapping the given atom would produce. * This routine is for the general case of an arbitrary potential and @@ -1573,8 +889,6 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeGeneric(int flipAtom, int old *********************************************************************/ double FixSemiGrandCanonicalMC::computeTotalEnergy() { - ASSERT(compute_pe != nullptr); - int eflag = 1; int vflag = 0; @@ -1639,7 +953,6 @@ void FixSemiGrandCanonicalMC::flipAtomEAM(int flipAtom, int flipAtomNL, int oldS int jnum = neighborList->numneigh[flipAtomNL]; for(int jj = 0; jj < jnum; jj++) { int j = jlist[jj]; - ASSERT(j < (int)changedAtoms.size()); double delx = xi - atom->x[j][0]; double dely = yi - atom->x[j][1]; @@ -1677,151 +990,6 @@ void FixSemiGrandCanonicalMC::flipAtomEAM(int flipAtom, int flipAtomNL, int oldS rho[flipAtom] = new_total_rho_i; } -#if CDEAM_MC_SUPPORT - -/********************************************************************* - * Flips the type of one atom and changes the electron densities and - * D values of nearby atoms accordingly. - * This routine is for the case of the concentration dependent CD-EAM potential. - * - * Parameters: - * - * flipAtom [in] - * This specifies the atom to be swapped. It's an index into the local list of atoms. - * - * flipAtomNL [in] - * This specifies the atom to be swapped. It's an index into the neighbor list. - * - * oldSpecies [in] - * The current species of the atom before the routine is called. - * - * newSpecies [in] - * The new type to be assigned to the atom. - *********************************************************************/ -void FixSemiGrandCanonicalMC::flipAtomCDEAM(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies) -{ - double p; - int m; - double* rho = pairEAM->rho; - double* coeff; - double new_total_rho_i = 0.0; - double new_total_rhoB_i = 0.0; - - atom->type[flipAtom] = newSpecies; - - // Rescale particle velocity vector to conserve kinetic energy. - double vScaleFactor = sqrt(atom->mass[oldSpecies] / atom->mass[newSpecies]); - atom->v[flipAtom][0] *= vScaleFactor; - atom->v[flipAtom][1] *= vScaleFactor; - atom->v[flipAtom][2] *= vScaleFactor; - - double xi = atom->x[flipAtom][0]; - double yi = atom->x[flipAtom][1]; - double zi = atom->x[flipAtom][2]; - - /// The change in atom i's D value. - double deltaD_i = 0.0; - - // Loop over all neighbors of the selected atom. - int* jlist = neighborList->firstneigh[flipAtomNL]; - int jnum = neighborList->numneigh[flipAtomNL]; - for(int jj = 0; jj < jnum; jj++) { - int j = jlist[jj]; - ASSERT(j < (int)changedAtoms.size()); - - double delx = xi - atom->x[j][0]; - double dely = yi - atom->x[j][1]; - double delz = zi - atom->x[j][2]; - double rsq = delx*delx + dely*dely + delz*delz; - if(rsq >= pairEAM->cutforcesq) continue; - - int jtype = atom->type[j]; - double r = sqrt(rsq); - p = r * pairEAM->rdr + 1.0; - m = static_cast(p); - m = MIN(m, pairEAM->nr - 1); - p -= m; - p = MIN(p, 1.0); - - // Calculate change of electron density at site j. - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[oldSpecies][jtype]][m]; - double oldrho_contr = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[newSpecies][jtype]][m]; - double newrho_contr = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - double delta_rho = newrho_contr - oldrho_contr; - - rho[j] += delta_rho; - // Sum total rho at site of swapped atom. - coeff = pairEAM->rhor_spline[pairEAM->type2rhor[jtype][newSpecies]][m]; - double rho_ji = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - new_total_rho_i += rho_ji; - if(jtype == pairCDEAM->speciesB) - new_total_rhoB_i += rho_ji; - - // Determine change of partial electron density at site j. - // It increases if the atom i becomes a B atom and it decreases if it was a B atom. - if(newSpecies == pairCDEAM->speciesB) - pairCDEAM->rhoB[j] += newrho_contr; - else if(oldSpecies == pairCDEAM->speciesB) - pairCDEAM->rhoB[j] -= oldrho_contr; - - /// The change in atom j's D value. - double deltaD_j; - - /// Was the old i-j interaction concentration dependent? - if((jtype == pairCDEAM->speciesA && oldSpecies == pairCDEAM->speciesB) - || (jtype == pairCDEAM->speciesB && oldSpecies == pairCDEAM->speciesA)) { - - // The old pair interaction was concentration dependent. - // Please note that it must now become concentration independent that - // either the A or the B atom has changed to another species. - // We here require that only the A-B interactions are concentration dependent. - - // Calculate change of pair energy ij. - coeff = pairEAM->z2r_spline[pairEAM->type2z2r[oldSpecies][jtype]][m]; - double oldz2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - // The D value of j decreases due to the removed cross interaction. - deltaD_j = -oldz2 / r; - } - else { - - // The old pair interaction was not concentration dependent. Now it might have - // become dependent. Do check: - if((jtype == pairCDEAM->speciesA && newSpecies == pairCDEAM->speciesB) - || (jtype == pairCDEAM->speciesB && newSpecies == pairCDEAM->speciesA)) { - - // The new pair interaction is concentration dependent. It's an AB interaction. - - // Calculate change of pair energy ij. - coeff = pairEAM->z2r_spline[pairEAM->type2z2r[newSpecies][jtype]][m]; - double newz2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - - // The D value of j increases due to the created cross interaction. - deltaD_j = newz2 / r; - } - else deltaD_j = 0.0; - } - pairCDEAM->D_values[j] += deltaD_j; - - // D_i changes in the same way as D_j changes. - deltaD_i += deltaD_j; - - // Set the flag for this atom to indicate that its rho has changed and needs - // to be transfered at end of MC step. - changedAtoms[j] = true; - } - - // Store newly calculated electron densities and D value at swapped atom site. - rho[flipAtom] = new_total_rho_i; - pairCDEAM->rhoB[flipAtom] = new_total_rhoB_i; - pairCDEAM->D_values[flipAtom] += deltaD_i; - - changedAtoms[flipAtom] = true; -} - -#endif - /********************************************************************* * Flips the type of one atom. * This routine is for the generic case. @@ -1858,12 +1026,11 @@ double FixSemiGrandCanonicalMC::compute_vector(int index) if(index == 0) return nAcceptedSwaps; if(index == 1) return nRejectedSwaps; index -= 1; - ASSERT(index >= 1 && index < (int)speciesCounts.size()); int totalAtoms = 0; for(int i = 0; i < (int)speciesCounts.size(); i++) totalAtoms += speciesCounts[i]; if(index <= atom->ntypes) return (double)speciesCounts[index] / (totalAtoms > 0 ? totalAtoms : 1); - ASSERT(false); return 0.0; + return 0.0; } /********************************************************************* diff --git a/src/MC/fix_sgcmc.h b/src/MC/fix_sgcmc.h index f74bc105e3..fd1f67baee 100644 --- a/src/MC/fix_sgcmc.h +++ b/src/MC/fix_sgcmc.h @@ -32,24 +32,6 @@ FixStyle(sgcmc,FixSemiGrandCanonicalMC); #include "fix.h" -// Setting this to 1 enables support for the concentration-dependent EAM potential (pair_style eam/cd) -// in the Monte Carlo routine. Setting to 0 limits support to standard EAM only and removes all dependencies -// on the CD-EAM potential code. -#ifndef CDEAM_MC_SUPPORT -#define CDEAM_MC_SUPPORT 0 -#endif - -// Setting this to 1 enables support for Tersoff-like potentials (pair_style tersoff) -// in the Monte Carlo routine. -#ifndef TERSOFF_MC_SUPPORT -#define TERSOFF_MC_SUPPORT 0 -#endif -// Setting this to 1 enables additional debugging/sanity checks (with a small performance penalty). -#ifndef SGCMC_DEBUG -#define SGCMC_DEBUG 0 -#endif - -#include #include namespace LAMMPS_NS { @@ -104,24 +86,6 @@ class FixSemiGrandCanonicalMC : public Fix { /// This routine is for the case of a standard EAM potential. double computeEnergyChangeEAM(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies); -#if CDEAM_MC_SUPPORT - /// Calculates the change in energy that swapping the given atom would produce. - /// This routine is for the case of the concentration dependent CD-EAM potential. - double computeEnergyChangeCDEAM(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies); -#endif - -#if TERSOFF_MC_SUPPORT - /// Calculates the change in energy that swapping the given atom would produce. - /// This routine is for the Tersoff potential. - double computeEnergyChangeTersoff(int flipAtom, int flipAtomNL, int oldSpecies, int newSpecies); - - /// Computes the energy of the atom group around the flipped atom using the Tersoff potential. - double computeEnergyTersoff(int flipAtom); - - /// Computes the energy of an atom using the Tersoff potential. - double computeAtomicEnergyTersoff(int i); -#endif - /// Calculates the change in energy that swapping the given atom would produce. /// This routine is for the general case of an arbitrary potential and /// IS VERY SLOW! It computes the total energies of the system for the unmodified state @@ -147,19 +111,6 @@ class FixSemiGrandCanonicalMC : public Fix { /// Transfers the locally changed electron densities and atom types to the neighbors. void communicateRhoAndTypes(); -#if SGCMC_DEBUG - /// Allocate atom-based array. - void grow_arrays(int) override; - /// Copy values within local atom-based array. - void copy_arrays(int, int) override; - /// Initialize one atom's array values, called when atom is created. - void set_arrays(int) override; - /// Pack values in local atom-based array for exchange with another proc. - int pack_exchange(int, double *) override; - /// Unpack values in local atom-based array from exchange with another proc. - int unpack_exchange(int, double *) override; -#endif - private: /// The number of MD steps between each MC cycle. int nevery_mdsteps; @@ -237,20 +188,6 @@ class FixSemiGrandCanonicalMC : public Fix { /// This is required to access the Rho arrays calculated by the potential class and its potential tables. class PairEAM *pairEAM; -#if CDEAM_MC_SUPPORT - /// Pointer to the CD-EAM potential class. - /// This is required to access the RhoB arrays calculated by the potential class. - /// The pointer is NULL if only the standard EAM model is used in the simulation. - class PairEAMCD *pairCDEAM; -#endif - -#if TERSOFF_MC_SUPPORT - /// Pointer to the Tersoff potential class. - /// This is required to access the parameters of the potential when computing the - /// change in energy. - class PairTersoff *pairTersoff; -#endif - /// This array contains a boolean value per atom (real and ghosts) that indicates whether /// the electron density or another property at that site has been affected by one of the accepted MC swaps. std::vector changedAtoms; @@ -272,14 +209,8 @@ class FixSemiGrandCanonicalMC : public Fix { /// A compute used to compute the total potential energy of the system. class Compute *compute_pe; - -#if SGCMC_DEBUG - /// This per-atom array counts how often each atom is picked for a trial move. - /// This is only used for debugging purposes. - double *trialCounters; -#endif }; } // namespace LAMMPS_NS -#endif // FIX_SEMIGRANDCANONICAL_MC_H -#endif // FIX_CLASS +#endif +#endif