diff --git a/README.md b/README.md index 2dffa74..f351cea 100644 --- a/README.md +++ b/README.md @@ -111,5 +111,6 @@ and whatever other options are needed). Installing the python package is most e ## Done + - combine all move types into a single LAMMPS fix, for improved performance - restart from snapshots - refactor `NSConfig.n_quantities` to be a class rather than instance attribute (minor) diff --git a/lammps_patches/new/src/NS/fix_ns.cpp b/lammps_patches/new/src/NS/fix_ns.cpp new file mode 100644 index 0000000..2cb04d2 --- /dev/null +++ b/lammps_patches/new/src/NS/fix_ns.cpp @@ -0,0 +1,1021 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DEBUG +#include +#endif + +#include +#include +#include "fix_ns.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "pair.h" +#include "modify.h" +#include "compute.h" +#include "update.h" +#include "random_mars.h" +#include "domain.h" +#include "memory.h" +#include "respa.h" +#include "error.h" +#include "math_extra.h" +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define MOVE_UNDEF -1 +#define MOVE_NONE 0 +#define MOVE_POS 1 +#define MOVE_CELL 2 +#define MOVE_TYPE 3 +#define MOVE_VOL 20 +#define MOVE_STRETCH 21 +#define MOVE_SHEAR 22 + +/* ---------------------------------------------------------------------- */ + +void normalize_cumul_probs_3(double &p1, double &p2, double &p3) { + p1 = (p1 < 0.0) ? 0.0 : p1; + p2 = (p2 < 0.0) ? 0.0 : p2; + p3 = (p3 < 0.0) ? 0.0 : p3; + double pSum = p1 + p2 + p3; + p1 /= pSum; p2 /= pSum; p3 /= pSum; + + // make into cumulative probabilities + p2 += p1; + if (p3 == 0.0) { + p2 = 1.0; + } + p3 = 1.0; +} + +FixNS::FixNS(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), dx(NULL), prevx(NULL), prevtype(NULL), mu(NULL) +{ + // possible values: base 3 + global 2 + probabilities 3 + 1? + 9? + (1 + ntype)? + if (strcmp(style,"ns") != 0 && narg <= 3 + 2 + 3) + error->all(FLERR,"Illegal fix ns command"); + + // before ns-specific params + int iarg = 3; + + // global NS params + seed = utils::inumeric(FLERR,arg[iarg++],true,lmp); + Emax = utils::numeric(FLERR,arg[iarg++],true,lmp); + + // probabilities of move types + pPos = utils::numeric(FLERR,arg[iarg++],true,lmp); + pCell = utils::numeric(FLERR,arg[iarg++],true,lmp); + pType = utils::numeric(FLERR,arg[iarg++],true,lmp); + + pos_n_steps = cell_n_steps = type_n_steps = 0; + + // pos GMC params + if (pPos > 0.0) { + if (narg - iarg < 2) + error->all(FLERR,"Illegal fix ns command when parsing pos GMC arguments"); + pos_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); + gmc_step_size = utils::numeric(FLERR,arg[iarg++],true,lmp); + memory->create(dx,atom->nmax,3,"ns:dx"); + memory->create(prevx,atom->nmax,3,"ns:prevx"); + } + + // cell MC params + if (pCell > 0.0) { + if (narg - iarg < 3 + 6 + 1) + error->all(FLERR,"Illegal fix ns command when parsing cell MC arguments"); + + if (!domain->triclinic) + error->all(FLERR,"fix ns cell moves require triclinic box (for now)"); + + cell_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); + min_aspect_ratio = utils::numeric(FLERR,arg[iarg++],true,lmp); + pressure = utils::numeric(FLERR,arg[iarg++],true,lmp); + + pVol = utils::numeric(FLERR,arg[iarg++],true, lmp); + dVol = utils::numeric(FLERR,arg[iarg++],true, lmp); + pStretch = utils::numeric(FLERR,arg[iarg++],true, lmp); + dStretch = utils::numeric(FLERR,arg[iarg++],true, lmp); + pShear = utils::numeric(FLERR,arg[iarg++],true, lmp); + dShear = utils::numeric(FLERR,arg[iarg++],true, lmp); + flat_V_prior = 1; + if (strcmp(arg[iarg],"no") == 0) flat_V_prior = 0; + else if (strcmp(arg[iarg],"yes") != 0) + error->all(FLERR,"Illegal fix ns flat_V_prior value"); + iarg++; + +#ifdef DEBUG +std::cout << "DEBUG initial vol/stretch/shear probs " << pVol << " " << pStretch << " " << pShear << std::endl; +#endif + normalize_cumul_probs_3(pVol, pStretch, pShear); +#ifdef DEBUG +std::cout << "DEBUG normalized cumulative vol/stretch/shear probs " << pVol << " " << pStretch << " " << pShear << std::endl; +#endif + + if (!prevx) + memory->create(prevx,atom->nmax,3,"ns:prevx"); + + // from fix_nh.cpp + // who knows what other things may be needed? + box_change |= (BOX_CHANGE_X | BOX_CHANGE_Y | BOX_CHANGE_Z | + BOX_CHANGE_YZ | BOX_CHANGE_XZ | BOX_CHANGE_XY); + no_change_box = 1; + } + + unequal_cutoffs = false; + // type MC params + if (pType > 0.0) { + if (narg - iarg < 2) + error->all(FLERR,"Illegal fix ns command when parsing type MC arguments"); + + type_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); + semi_GC_flag = 1; + if (strcmp(arg[iarg],"no") == 0) semi_GC_flag = 0; + else if (strcmp(arg[iarg],"yes") != 0) + error->all(FLERR,"Illegal fix ns/type semi_GC_flag value"); + iarg++; + + if (semi_GC_flag) { + if (narg - iarg < atom->ntypes) + error->all(FLERR,"Illegal fix ns command when parsing type MC mu arguments"); + + mu = new double[atom->ntypes]; + for (int i=0; i < atom->ntypes; i++) { + mu[i] = utils::numeric(FLERR,arg[iarg++],true,lmp); + } + } else { + mu = 0; + } + + // based on fix_atom_swap.cpp + int ntypes = atom->ntypes; + double **cutsq = force->pair->cutsq; + for (int itype = 1; itype <= ntypes; itype++) + for (int jtype = 1; jtype <= ntypes; jtype++) + for (int ktype = 1; ktype <= ntypes; ktype++) + if (cutsq[itype][ktype] != cutsq[jtype][ktype]) + unequal_cutoffs = true; + if (unequal_cutoffs) + error->warning(FLERR,"unequal cutoffs may cause a problem for atom type moves"); + + memory->create(prevtype,atom->nmax,"ns:prevtype"); +#ifdef SWAP_POS + if (!prevx) + memory->create(prevx,atom->nmax,3,"ns:prevx"); +#endif + } + +#ifdef DEBUG +std::cout << "DEBUG initial pos/cell/type probs " << pPos << " " << pCell << " " << pType << std::endl; +#endif + normalize_cumul_probs_3(pPos, pCell, pType); +#ifdef DEBUG +std::cout << "DEBUG cumulative normalized pos/cell/type probs " << pPos << " " << pCell << " " << pType << std::endl; +#endif + + max_n_steps = (pos_n_steps > cell_n_steps) ? pos_n_steps : cell_n_steps; + max_n_steps = (max_n_steps > type_n_steps) ? max_n_steps : type_n_steps; + + if (narg - iarg < 0) + error->all(FLERR,"Illegal fix ns command - extra arguments left after parsing"); + + cell_cur_move = MOVE_UNDEF; + state_cur_move = MOVE_UNDEF; + state_traj_steps_remaining = 0; + + // zero accumulated energy shifts + cumulative_dPV = cumulative_dmuN = 0.0; + // attempt frequency counters + vector_flag = 1; + size_vector = 2 + 6 + 2; + global_freq = 1; + extvector = 0; + n_attempt_pos = n_success_pos = + n_attempt_vol = n_success_vol = + n_attempt_stretch = n_success_stretch = + n_attempt_shear = n_success_shear = + n_attempt_type = n_success_type = 0; + + // random_g will be used when all processes need to agree on the random + // number, so it uses the same seed on all. + // random_l will be used when each processes needs to do something different, + // and it may be called an unpredictable number of times, which will make it + // diverge between the different processes. + random_g = new RanMars(lmp,seed); + random_l = new RanMars(lmp,seed + 1 + comm->me); + + // copied from fix_gmc.cpp + dynamic_group_allow = 1; + time_integrate = 1; + // copied from atom_swap.cpp? + // set force_reneighbor so next_reneighbor is checked, but don't force it + // for the next step until we see what sort of step this is - only type + // steps require (?) reneighbors + force_reneighbor = 1; + // next_reneighbor = update->ntimestep + 1; + next_reneighbor = -1; + + // from MC/fix_atom_swap.cpp + // amount of data to be packed for inter-process comm (atom type) + if (atom->q_flag) + comm_forward = 2; + else + comm_forward = 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixNS::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= PRE_EXCHANGE; + mask |= FINAL_INTEGRATE; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + return mask; +} + +// from MC/fix_atom_swap.cpp +int FixNS::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + + int *type = atom->type; + double *q = atom->q; + + m = 0; + + if (atom->q_flag) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = type[j]; + buf[m++] = q[j]; + } + } else { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = type[j]; + } + } + + return m; +} + +// from MC/fix_atom_swap.cpp +void FixNS::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + int *type = atom->type; + double *q = atom->q; + + m = 0; + last = first + n; + + if (atom->q_flag) { + for (i = first; i < last; i++) { + type[i] = static_cast(buf[m++]); + q[i] = buf[m++]; + } + } else { + for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); + } +} + +void FixNS::init() +{ + + int id = modify->find_compute("thermo_pe"); + + modify->compute[id]->invoked_scalar = -1; + pe_compute = modify->compute[id]; + pe_compute->addstep(update->ntimestep+1); + + if (strstr(update->integrate_style,"respa")) + error->all(FLERR,"fix ns not compatible with RESPA"); +} + +void FixNS::pre_exchange() +{ + if (state_cur_move == MOVE_TYPE) { + // force reneighbor for swaps (maybe only needed on accept, + // or only with unequal cutoffs?) + next_reneighbor = update->ntimestep + 1; + } else { + next_reneighbor = -1; + } +} + +void FixNS::initial_integrate(int vflag) +{ + // state machine for move types + if (state_traj_steps_remaining <= 0) { + // for now, give up on moving if there aren't enough steps left + // to accommodate longest (across 3 step types) mini-traj. This + // should ensure correct ratio of step types + if (update->laststep + 1 - update->ntimestep < max_n_steps) { +#ifdef DEBUG +std::cout << "initial_integrate not enough steps left " << update->laststep << " + 1 - " << update->ntimestep << " < " << max_n_steps << std::endl; +#endif + // no time for another mini traj + state_cur_move = MOVE_NONE; + } else { + // pick new move type, all processes must agree + double rv = random_g->uniform(); +#ifdef DEBUG +std::cout << "initial_integrate pick rv " << rv << " probs " << pPos << " " << pCell << " " << pType << std::endl; +#endif + if (rv < pPos) { + state_cur_move = MOVE_POS; + pos_gmc_traj_prep(); + state_traj_steps_remaining = pos_n_steps; + } else if (rv < pCell) { + state_cur_move = MOVE_CELL; + state_traj_steps_remaining = cell_n_steps; + } else { + state_cur_move = MOVE_TYPE; + state_traj_steps_remaining = type_n_steps; + } + } + } +#ifdef DEBUG +else +std::cout << "initial_integrate continue with current move type" << std::endl; +#endif + + state_traj_steps_remaining--; + + if (state_cur_move == MOVE_POS) { + pos_gmc_initial_integrate(); + } else if (state_cur_move == MOVE_CELL) { + cell_initial_integrate(); + } else if (state_cur_move == MOVE_TYPE) { + type_initial_integrate(); + // move swapped or perturbed types, communicate to ghost atoms + comm->forward_comm(this); + } + + int id = modify->find_compute("thermo_pe"); + + modify->compute[id]->invoked_scalar = -1; + pe_compute = modify->compute[id]; + pe_compute->addstep(update->ntimestep+1); +} + +void FixNS::pos_gmc_traj_prep() +{ + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + double **x = atom->x; + + double dx2sum = 0; + + // if we ever want to maintain coherence of dx from traj to traj, we can do it + // using atomic velocities: + // double **v = atom->v; + + for(int i = 0; i < nlocal; i++) { + // random direction + dx[i][0] = random_l->gaussian(); + dx[i][1] = random_l->gaussian(); + dx[i][2] = random_l->gaussian(); + dx2sum += dx[i][0]*dx[i][0]+dx[i][1]*dx[i][1]+dx[i][2]*dx[i][2]; + + // save initial pos in case we have to reject move + prevx[i][0] = x[i][0]; + prevx[i][1] = x[i][1]; + prevx[i][2] = x[i][2]; + } + + // WARNING this really needs to be summed over all MPI processes + dx2sum = sqrt(dx2sum); + for(int i = 0; i < nlocal; i++) { + dx[i][0] /= dx2sum; + dx[i][1] /= dx2sum; + dx[i][2] /= dx2sum; + } + + n_attempt_pos += 1; +} + +void FixNS::pos_gmc_initial_integrate() +{ +#ifdef DEBUG +std::cout << "POS pos_gmc_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; +#endif + // update x of atoms in group + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + x[i][0] += gmc_step_size * dx[i][0]; + x[i][1] += gmc_step_size * dx[i][1]; + x[i][2] += gmc_step_size * dx[i][2]; + } + +} + +void FixNS::cell_initial_integrate() +{ +#ifdef DEBUG +std::cout << "CELL cell_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; +#endif + // do move (unles rejected) + + int natoms = atom->natoms; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + double **x = atom->x; + + // save previous cell and positions + for (int i=0; i < 3; i++) { + prev_boxlo[i] = domain->boxlo[i]; + prev_boxhi[i] = domain->boxhi[i]; + } + prev_xy = domain->xy; + prev_yz = domain->yz; + prev_xz = domain->xz; + + for (int iat=0; iat < nlocal; iat++) + for (int j=0; j < 3; j++) + prevx[iat][j] = x[iat][j]; + + double boxext[3]; + for (int i=0; i < 3; i++) + boxext[i] = domain->boxhi[i] - domain->boxlo[i]; + + cell_cur_move = MOVE_UNDEF; + cell_move_rejected_early = false; + + // keep track of change PV in this step + // default to 0, set otherwise when volume move is selected + dPV = 0.0; + + // default to accept, will be changed below to reject (-1.0) or probablistic (if flat_V_prior != 1) + double rv = random_g->uniform(); + double new_cell[3][3]; + // pick a step type + if (rv < pVol) { + n_attempt_vol++; + cell_cur_move = MOVE_VOL; + // volume step + double orig_V; + if (domain->dimension == 3) orig_V = domain->xprd * domain->yprd * domain->zprd; + else orig_V = domain->xprd * domain->yprd; + + double dV = random_g->gaussian(0.0, dVol); + double new_V = orig_V + dV; +#ifdef DEBUG + std::cout << "CELL VOLUME volumetric strain " << dV << " / " << orig_V << " = " << dV/orig_V << " "; +#endif + if (new_V/orig_V < 0.5) { + cell_move_rejected_early = true; +#ifdef DEBUG + std::cout << "CELL REJECT new_V/orig_V < 0.5" << std::endl; +#endif + } else { + if (flat_V_prior == 0 && new_V < orig_V && random_g->uniform() > pow(new_V / orig_V, natoms)) { + cell_move_rejected_early = true; +#ifdef DEBUG + std::cout << "CELL REJECT V probability " << pow(new_V / orig_V, natoms) << std::endl; +#endif + } else { + double transform_diag = pow(new_V / orig_V, 1.0/3.0); + new_cell[0][0] = boxext[0] * transform_diag; + new_cell[0][1] = 0.0; + new_cell[0][2] = 0.0; + new_cell[1][0] = domain->xy * transform_diag; + new_cell[1][1] = boxext[1] * transform_diag; + new_cell[1][2] = 0.0; + new_cell[2][0] = domain->xz * transform_diag; + new_cell[2][1] = domain->yz * transform_diag; + new_cell[2][2] = boxext[2] * transform_diag; + dPV = pressure * dV; + } + } + } else if (rv < pStretch) { + n_attempt_stretch++; + cell_cur_move = MOVE_STRETCH; + // stretch step + // pick directions to use v_ind and (v_ind+1)%3 +#ifdef DEBUG + std::cout << "CELL STRETCH "; +#endif + int v_ind = int(3 * random_g->uniform()) % 3; + rv = random_g->gaussian(0.0, dStretch); + double transform_diag[3]; + transform_diag[v_ind] = exp(rv); + transform_diag[(v_ind+1)%3] = exp(-rv); + transform_diag[(v_ind+2)%3] = 1.0; +#ifdef DEBUG + std::cout << transform_diag[0] << " " << transform_diag[1] << " " << transform_diag[2] << " "; +#endif + // create new cell for aspect ratio test and new domain + new_cell[0][0] = boxext[0] * transform_diag[0]; + new_cell[0][1] = 0.0; + new_cell[0][2] = 0.0; + new_cell[1][0] = domain->xy * transform_diag[0]; + new_cell[1][1] = boxext[1] * transform_diag[1]; + new_cell[1][2] = 0.0; + new_cell[2][0] = domain->xz * transform_diag[0]; + new_cell[2][1] = domain->yz * transform_diag[1]; + new_cell[2][2] = boxext[2] * transform_diag[2]; + if (cur_aspect_ratio(new_cell) < min_aspect_ratio) { + cell_move_rejected_early = true; +#ifdef DEBUG + std::cout << "CELL REJECT min_aspect_ratio " << cur_aspect_ratio(new_cell) << std::endl; +#endif + } + } else { + n_attempt_shear++; + cell_cur_move = MOVE_SHEAR; + // shear step + // save original cell and a temporary cell +#ifdef DEBUG + std::cout << "CELL SHEAR "; +#endif + double orig_cell[3][3], t_cell[3][3]; + t_cell[0][0] = orig_cell[0][0] = boxext[0]; + t_cell[0][1] = orig_cell[0][1] = 0.0; + t_cell[0][2] = orig_cell[0][2] = 0.0; + t_cell[1][0] = orig_cell[1][0] = domain->xy; + t_cell[1][1] = orig_cell[1][1] = boxext[1]; + t_cell[1][2] = orig_cell[1][2] = 0.0; + t_cell[2][0] = orig_cell[2][0] = domain->xz; + t_cell[2][1] = orig_cell[2][1] = domain->yz; + t_cell[2][2] = orig_cell[2][2] = boxext[2]; + + // pick vector to perturb + int vec_ind = int(3 * random_g->uniform()) % 3; + + // perturb t_cell[vec_ind] along other two vectors + double vhat[3]; + for (int di=1; di < 3; di++) { + for (int i=0; i < 3; i++) + vhat[i] = orig_cell[(vec_ind + di) % 3][i]; + MathExtra::norm3(vhat); + rv = random_g->gaussian(0.0, dShear); +#ifdef DEBUG + std::cout << (vec_ind + di) % 3 << " " << rv << " "; +#endif + for (int i=0; i < 3; i++) + t_cell[vec_ind][i] += rv * vhat[i]; + } + + if (cur_aspect_ratio(t_cell) < min_aspect_ratio) { + cell_move_rejected_early = true; +#ifdef DEBUG + std::cout << "CELL REJECT min_aspect_ratio " << cur_aspect_ratio(t_cell) << std::endl; +#endif + } else { + // rotate new_cell back to LAMMPS orientation + + double a0_norm = sqrt(MathExtra::dot3(t_cell[0], t_cell[0])); + new_cell[0][0] = a0_norm; + new_cell[0][1] = 0.0; + new_cell[0][2] = 0.0; + + double a0_hat[3] = {t_cell[0][0], t_cell[0][1], t_cell[0][2]}; MathExtra::norm3(a0_hat); + new_cell[1][0] = MathExtra::dot3(t_cell[1], a0_hat); + double a0_hat_cross_a1[3]; MathExtra::cross3(a0_hat, t_cell[1], a0_hat_cross_a1); + new_cell[1][1] = sqrt(MathExtra::dot3(a0_hat_cross_a1, a0_hat_cross_a1)); + new_cell[1][2] = 0.0; + + double a0_cross_a1_hat[3] = {a0_hat_cross_a1[0], a0_hat_cross_a1[1], a0_hat_cross_a1[2]}; + MathExtra::norm3(a0_cross_a1_hat); + + double a0_cross_a1_hat_cross_a0_hat[3]; + MathExtra::cross3(a0_cross_a1_hat, a0_hat, a0_cross_a1_hat_cross_a0_hat); + new_cell[2][0] = MathExtra::dot3(t_cell[2], a0_hat); + new_cell[2][1] = MathExtra::dot3(t_cell[2], a0_cross_a1_hat_cross_a0_hat); + new_cell[2][2] = MathExtra::dot3(t_cell[2], a0_cross_a1_hat); + } + } + + if (!cell_move_rejected_early) { + // apply move to box + double boxctr[3]; + for (int i=0; i < 3; i++) + boxctr[i] = 0.5*(prev_boxhi[0] + prev_boxlo[0]); + + domain->x2lamda(nlocal); + + boxext[0] = new_cell[0][0]; + domain->xy = new_cell[1][0]; + boxext[1] = new_cell[1][1]; + domain->xz = new_cell[2][0]; + domain->yz = new_cell[2][1]; + boxext[2] = new_cell[2][2]; + + for (int i=0; i < 3; i++) { + domain->boxlo[i] = boxctr[i] - boxext[i]/2.0; + domain->boxhi[i] = boxctr[i] + boxext[i]/2.0; + } + + domain->set_global_box(); + domain->set_local_box(); + + domain->lamda2x(nlocal); +#ifdef DEBUG + std::cout << std::endl; +#endif + } + +} + +void FixNS::type_initial_integrate() +{ +#ifdef DEBUG +std::cout << "TYPE type_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; +#endif + // do move (unless rejected) + + n_attempt_type++; + + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + int ntypes = atom->ntypes; + + int *type = atom->type; + for (int iat=0; iat < nlocal; iat++) + prevtype[iat] = type[iat]; + +#ifdef SWAP_POS + double **x = atom->x; + for (int iat=0; iat < nlocal; iat++) + for (int jat=0; jat < 3; jat++) + prevx[iat][jat] = x[iat][jat]; +#endif + + // WARNING move is restricted to only one proc - fine for semi-GC, but + // limiting for fixed composition swaps + double rv = random_g->uniform(); + int rnd_proc = static_cast(comm->nprocs * rv) % comm->nprocs; + if (comm->me != rnd_proc) + return; + + rv = random_l->uniform(); + int atom_0 = static_cast(nlocal * rv) % nlocal; + + if (semi_GC_flag) { + // perturb type randomly, store change to mu*N term + rv = random_l->uniform(); + int dtype = 1 + static_cast(rv * (ntypes-1)) % (ntypes-1); + int new_type = 1 + (((type[atom_0]-1) + dtype) % ntypes); + dmuN = mu[new_type-1] - mu[type[atom_0]-1]; + type[atom_0] = new_type; +#ifdef DEBUG + std::cout << "TYPE SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << + type[atom_0] << " mu " << mu[type[atom_0]-1] << " " << std::endl; +#endif + } else { + // not semi-GC, swap a pair of atoms + bool all_same = true; + for (int i=1; i < nlocal; i++) + if (type[0] != type[i]) { + all_same = false; + break; + } + if (all_same) + return; + + int atom_1 = atom_0; + while (type[atom_0] == type[atom_1]) { + rv = random_l->uniform(); + atom_1 = static_cast(nlocal * rv) % nlocal; + } + +#ifdef DEBUG +#ifdef SWAP_POS + std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << + " i " << atom_1 << " type " << type[atom_1] << " x " << x[atom_1][0] << " " << x[atom_1][1] << " " << x[atom_1][2] << " " << std::endl; +#else + std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " <-> i " << atom_1 << " type " << type[atom_1] << " " << std::endl; +#endif +#endif + +#ifdef SWAP_POS + double tx[3]; + tx[0] = x[atom_0][0]; + tx[1] = x[atom_0][1]; + tx[2] = x[atom_0][2]; + x[atom_0][0] = x[atom_1][0]; + x[atom_0][1] = x[atom_1][1]; + x[atom_0][2] = x[atom_1][2]; + x[atom_1][0] = tx[0]; + x[atom_1][1] = tx[1]; + x[atom_1][2] = tx[2]; +#else + int t_type = type[atom_0]; + type[atom_0] = type[atom_1]; + type[atom_1] = t_type; +#endif + + dmuN = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNS::final_integrate() +{ + if (state_cur_move == MOVE_POS) { + pos_gmc_final_integrate(); + } else if (state_cur_move == MOVE_CELL) { + cell_final_integrate(); + } else if (state_cur_move == MOVE_TYPE) { + type_final_integrate(); + } +} + +void FixNS::pos_gmc_final_integrate() +{ + // if potential energy is above Emax then want to modify dx with + // forces to change trajectory + + double **f = atom->f; + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "POS ecurrent " << ecurrent << " + " << cumulative_dPV << " - " << cumulative_dmuN; +#endif + ecurrent += cumulative_dPV - cumulative_dmuN; +#ifdef DEBUG +std::cout << " = " << ecurrent << std::endl; +#endif + + if (ecurrent < Emax) { +#ifdef DEBUG +std::cout << "POS ecurrent + dPV - dmuN = " << ecurrent << " < " << Emax << std::endl; +#endif + // accept + if (state_traj_steps_remaining == 0) { +#ifdef DEBUG +std::cout << "POS ACCEPT end of traj, increment n_success_pos" << std::endl; +#endif + // final step of mini traj, accept + n_success_pos += 1; + } +#ifdef DEBUG +std::cout << "POS mid traj, continuing" << std::endl; +#endif + // accept new position and continue straight + return; + } + +#ifdef DEBUG +std::cout << "POS ecurrent + dPV - dmuN = " << ecurrent << " >= " << Emax << std::endl; +#endif + + if (state_traj_steps_remaining == 0) { +#ifdef DEBUG +std::cout << "POS REJECT end of traj, revert positions to prevx" << std::endl; +#endif + // final step of mini traj, reject + for (int i = 0; i < nlocal; i++) { + x[i][0] = prevx[i][0]; + x[i][1] = prevx[i][1]; + x[i][2] = prevx[i][2]; + } + return; + } + +#ifdef DEBUG +std::cout << "POS reflect" << std::endl; +#endif + + // try to reflect from boundary + double fsum = 0; + double fhatdotdx = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + fsum += f[i][0]*f[i][0]+ + f[i][1]*f[i][1]+ + f[i][2]*f[i][2]; + fsum = sqrt(fsum); + + // should also detect nan + if (fsum != 0.0) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + fhatdotdx += f[i][0]/fsum*dx[i][0]; + fhatdotdx += f[i][1]/fsum*dx[i][1]; + fhatdotdx += f[i][2]/fsum*dx[i][2]; + } + } + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dx[i][0] -= 2*f[i][0]/fsum*fhatdotdx; + dx[i][1] -= 2*f[i][1]/fsum*fhatdotdx; + dx[i][2] -= 2*f[i][2]/fsum*fhatdotdx; + } + } + } + +} + + +void FixNS::cell_final_integrate() +{ + // if potential energy is above Emax then reject move + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "CELL ecurrent " << ecurrent << " - " << cumulative_dmuN; +#endif + ecurrent -= cumulative_dmuN; +#ifdef DEBUG +std::cout << " = " << ecurrent << std::endl; +#endif + + if (cell_move_rejected_early) { + return; + } + + // if potential energy + d(P V) is above Emax then reject move + // need to include previous steps' cumulative dPV contributions, as well as current ones + if (ecurrent + cumulative_dPV + dPV >= Emax) { +#ifdef DEBUG + std::cout << "CELL REJECT E == " << ecurrent << " + " << cumulative_dPV << " + " << dPV << " >= Emax == " << Emax << std::endl; +#endif + // reject move, so don't touch cumulative_dPV, since type change that led to current dPV was reverted + + for (int i=0; i < 3; i++) { + domain->boxlo[i] = prev_boxlo[i]; + domain->boxhi[i] = prev_boxhi[i]; + } + domain->xy = prev_xy; + domain->yz = prev_yz; + domain->xz = prev_xz; + + domain->set_global_box(); + domain->set_local_box(); + + double **x = atom->x; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + for (int iat=0; iat < nlocal; iat++) + for (int j=0; j < 3; j++) + x[iat][j] = prevx[iat][j]; + + } else { + // accept move, so accumulate dPV contribution from this step + cumulative_dPV += dPV; + switch (cell_cur_move) { + case MOVE_VOL: + n_success_vol++; + break; + case MOVE_STRETCH: + n_success_stretch++; + break; + case MOVE_SHEAR: + n_success_shear++; + break; + default: + error->all(FLERR,"Illegal value of cell_cur_move in increment n_success_*"); + } +#ifdef DEBUG +double new_cell[3][3]; +new_cell[0][0] = domain->boxhi[0] - domain->boxlo[0]; +new_cell[0][1] = 0.0; +new_cell[0][2] = 0.0; +new_cell[1][0] = domain->xy; +new_cell[1][1] = domain->boxhi[1] - domain->boxlo[1]; +new_cell[1][2] = 0.0; +new_cell[2][0] = domain->xz; +new_cell[2][1] = domain->yz; +new_cell[2][2] = domain->boxhi[2] - domain->boxlo[2]; + std::cout << "CELL ACCEPT E == " << ecurrent << " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << cur_aspect_ratio(new_cell) << std::endl; + // std::cout << "final cell " << new_cell[0][0] << " " << new_cell[0][1] << " " << new_cell[0][2] << std::endl; + // std::cout << " " << new_cell[1][0] << " " << new_cell[1][1] << " " << new_cell[1][2] << std::endl; + // std::cout << " " << new_cell[2][0] << " " << new_cell[2][1] << " " << new_cell[2][2] << std::endl; +#endif + } +} + + +void FixNS::type_final_integrate() +{ + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "TYPE ecurrent " << ecurrent << " + " << cumulative_dPV; +#endif + ecurrent += cumulative_dPV; +#ifdef DEBUG +std::cout << " = " << ecurrent << std::endl; +#endif + + // if potential energy - d(mu N) is above Emax then reject move + // need to include previous steps' cumulative dmuN contributions, as well as current ones + if (ecurrent - cumulative_dmuN - dmuN >= Emax) { +#ifdef DEBUG + std::cout << "TYPE REJECT E == " << ecurrent << " - " << cumulative_dmuN << " - " << dmuN << " >= Emax == " << Emax << " " << std::endl; +#endif + // reject move, so don't touch cumulative_dmuN, since type change that led to current dmuN was reverted + + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + int *type = atom->type; +// With SWAP_POS, type is only modified for semi-GC, so only has to be undone in that case, and x +// has to be reverted to undo swap move. +// By default (no SWAP_POS), type is modified for both swap and semi-GC, so undoing either step type +// uses prevtype +#ifdef SWAP_POS + if (semi_GC_flag) { +#endif + for (int iat=0; iat < nlocal; iat++) + type[iat] = prevtype[iat]; +#ifdef SWAP_POS + } else { + double **x = atom->x; + for (int iat=0; iat < nlocal; iat++) + for (int jat=0; jat < 3; jat++) + x[iat][jat] = prevx[iat][jat]; + } +#endif + + // modified type by restoring previous, communicate to ghost atoms + if (state_cur_move == MOVE_TYPE) { + comm->forward_comm(this); + } + + } else { + // accept move, so accumulate dmuN contribution from this step + cumulative_dmuN += dmuN; + n_success_type++; +#ifdef DEBUG + std::cout << "TYPE ACCEPT E == " << ecurrent << " - " << cumulative_dmuN << " < Emax == " << Emax << std::endl; +#endif + } +} + +FixNS::~FixNS() +{ + delete random_g; + delete random_l; + + if (dx) + memory->destroy(dx); + if (prevx) + memory->destroy(prevx); + if (prevtype) + memory->destroy(prevtype); + if (mu) + delete mu; +} + +double FixNS::cur_aspect_ratio(double cell[3][3]) +{ + double min_val = std::numeric_limits::max(); + for (int i=0; i < 3; i++) { + double vnorm_hat[3]; + MathExtra::cross3(cell[(i+1)%3], cell[(i+2)%3], vnorm_hat); + MathExtra::norm3(vnorm_hat); + double val = fabs(MathExtra::dot3(vnorm_hat, cell[i])); + min_val = (val < min_val) ? val : min_val; + } + + double V; + if (domain->dimension == 3) V = domain->xprd * domain->yprd * domain->zprd; + else V = domain->xprd * domain->yprd; + + return min_val / pow(V, 1.0/3.0); +} + +double FixNS::compute_vector(int n) +{ + switch (n) { + case 0: return n_attempt_pos; + case 1: return n_success_pos; + case 2: return n_attempt_vol; + case 3: return n_success_vol; + case 4: return n_attempt_stretch; + case 5: return n_success_stretch; + case 6: return n_attempt_shear; + case 7: return n_success_shear; + case 8: return n_attempt_type; + case 9: return n_success_type; + } + return -1.0; +} diff --git a/lammps_patches/new/src/NS/fix_ns.h b/lammps_patches/new/src/NS/fix_ns.h new file mode 100644 index 0000000..0f7ac9f --- /dev/null +++ b/lammps_patches/new/src/NS/fix_ns.h @@ -0,0 +1,106 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(ns,FixNS) + +#else + +#ifndef LMP_FIX_NS_H +#define LMP_FIX_NS_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNS : public Fix { + public: + FixNS(class LAMMPS *, int, char **); + ~FixNS() override; + int setmask() override; + void init() override; + void pre_exchange() override; + void initial_integrate(int) override; + void final_integrate() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + double compute_vector(int) override; + + protected: + class Compute *pe_compute; + double **dx, **prevx; + int *prevtype; + double prev_boxlo[3], prev_boxhi[3], prev_xy, prev_yz, prev_xz; + // global params + int seed; + double Emax; + double pPos, pCell, pType; + int state_cur_move; + int state_traj_steps_remaining; + // pos GMC params + int pos_n_steps; + double gmc_step_size; + // cell MC params + int cell_n_steps; + double min_aspect_ratio, pressure; + int flat_V_prior; + double pVol, dVol, pStretch, dStretch, pShear, dShear; + // type MC params + int type_n_steps; + int semi_GC_flag; + double *mu; + // other data + class RanMars *random_g, *random_l; + int peflag; + char *id_pe; + int max_n_steps; + + // pos GMC internal data + int n_attempt_pos, n_success_pos; + // cell MC internal data + double dPV, cumulative_dPV; + bool cell_move_rejected_early; + int cell_cur_move; + int n_attempt_vol, n_attempt_stretch, n_attempt_shear, n_success_vol, n_success_stretch, n_success_shear; + // type MC internal data + double dmuN, cumulative_dmuN; + bool unequal_cutoffs; + int n_attempt_type, n_success_type; + + private: + double cur_aspect_ratio(double cell[3][3]); + void pos_gmc_traj_prep(); + + void pos_gmc_initial_integrate(); + void cell_initial_integrate(); + void type_initial_integrate(); + void pos_gmc_final_integrate(); + void cell_final_integrate(); + void type_final_integrate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/lammps_patches/new/src/NS/fix_ns_cellmc.cpp b/lammps_patches/new/src/NS/fix_ns_cellmc.cpp index ce3a707..7f05b60 100644 --- a/lammps_patches/new/src/NS/fix_ns_cellmc.cpp +++ b/lammps_patches/new/src/NS/fix_ns_cellmc.cpp @@ -364,6 +364,9 @@ void FixNSCellMC::final_integrate() { // if potential energy is above Emax then reject move double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "CELL ecurrent " << ecurrent << std::endl; +#endif if (move_rejected_early) { return; @@ -373,7 +376,7 @@ void FixNSCellMC::final_integrate() // need to include previous steps' cumulative dPV contributions, as well as current ones if (ecurrent + cumulative_dPV + dPV >= Emax) { #ifdef DEBUG - std::cout << "REJECT E == " << ecurrent << " + " << cumulative_dPV + dPV << " >= Emax == " << Emax << std::endl; + std::cout << "REJECT E == " << ecurrent << " + " << cumulative_dPV << " + " << dPV << " >= Emax == " << Emax << std::endl; #endif // reject move, so don't touch cumulative_dPV, since type change that led to current dPV was reverted @@ -421,7 +424,7 @@ new_cell[1][2] = 0.0; new_cell[2][0] = domain->xz; new_cell[2][1] = domain->yz; new_cell[2][2] = domain->boxhi[2] - domain->boxlo[2]; - std::cout << "ACCEPT E == " << ecurrent " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << min_aspect_ratio_val(new_cell) << std::endl; + std::cout << "ACCEPT E == " << ecurrent << " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << min_aspect_ratio_val(new_cell) << std::endl; // std::cout << "final cell " << new_cell[0][0] << " " << new_cell[0][1] << " " << new_cell[0][2] << std::endl; // std::cout << " " << new_cell[1][0] << " " << new_cell[1][1] << " " << new_cell[1][2] << std::endl; // std::cout << " " << new_cell[2][0] << " " << new_cell[2][1] << " " << new_cell[2][2] << std::endl; diff --git a/lammps_patches/new/src/NS/fix_ns_cellmc.h b/lammps_patches/new/src/NS/fix_ns_cellmc.h index 778a268..b636603 100644 --- a/lammps_patches/new/src/NS/fix_ns_cellmc.h +++ b/lammps_patches/new/src/NS/fix_ns_cellmc.h @@ -27,12 +27,12 @@ namespace LAMMPS_NS { class FixNSCellMC : public Fix { public: FixNSCellMC(class LAMMPS *, int, char **); - ~FixNSCellMC(); - int setmask(); - virtual void init(); - virtual void initial_integrate(int); - virtual void final_integrate(); - double compute_vector(int); + ~FixNSCellMC() override; + int setmask() override; + void init() override; + void initial_integrate(int) override; + void final_integrate() override; + double compute_vector(int) override; protected: class Compute *pe_compute; @@ -46,7 +46,6 @@ class FixNSCellMC : public Fix { double pVol, dVol, pStretch, dStretch, pShear, dShear; int peflag; char *id_pe; - char str[64]; double dPV, cumulative_dPV; diff --git a/lammps_patches/new/src/NS/fix_ns_gmc.cpp b/lammps_patches/new/src/NS/fix_ns_gmc.cpp index 1669cb7..68eba4f 100644 --- a/lammps_patches/new/src/NS/fix_ns_gmc.cpp +++ b/lammps_patches/new/src/NS/fix_ns_gmc.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#ifdef DEBUG +#include +#endif + #include #include #include "fix_ns_gmc.h" @@ -150,11 +154,18 @@ void FixNSGMC::final_integrate() int nlocal = atom->nlocal; double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "POS ecurrent " << ecurrent << std::endl; +#endif if (ecurrent < Emax) // continue straight return; +#ifdef DEBUG +std::cout << "POS " << ecurrent << " >= " << Emax << " reflect" << std::endl; +#endif + // try to reflect from boundary double fsum = 0; double fhatdotdx = 0; diff --git a/lammps_patches/new/src/NS/fix_ns_gmc.h b/lammps_patches/new/src/NS/fix_ns_gmc.h index c085c0e..afd6571 100644 --- a/lammps_patches/new/src/NS/fix_ns_gmc.h +++ b/lammps_patches/new/src/NS/fix_ns_gmc.h @@ -27,11 +27,11 @@ namespace LAMMPS_NS { class FixNSGMC : public Fix { public: FixNSGMC(class LAMMPS *, int, char **); - ~FixNSGMC(); - int setmask(); - virtual void init(); - virtual void initial_integrate(int); - virtual void final_integrate(); + ~FixNSGMC() override; + int setmask() override; + void init() override; + void initial_integrate(int) override; + void final_integrate() override; protected: class Compute *pe_compute; @@ -42,7 +42,6 @@ class FixNSGMC : public Fix { double step_size; int peflag; char *id_pe; - char str[64]; }; } diff --git a/lammps_patches/new/src/NS/fix_ns_type.cpp b/lammps_patches/new/src/NS/fix_ns_type.cpp index d5c7fc2..09c80b8 100644 --- a/lammps_patches/new/src/NS/fix_ns_type.cpp +++ b/lammps_patches/new/src/NS/fix_ns_type.cpp @@ -98,6 +98,13 @@ FixNSType::FixNSType(LAMMPS *lmp, int narg, char **arg) : force_reneighbor = 1; next_reneighbor = update->ntimestep + 1; + // + // from MC/fix_atom_swap.cpp + // amount of data to be packed for inter-process comm (atom type) + if (atom->q_flag) + comm_forward = 2; + else + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -114,6 +121,53 @@ int FixNSType::setmask() return mask; } +// from MC/fix_atom_swap.cpp +int FixNSType::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + + int *type = atom->type; + double *q = atom->q; + + m = 0; + + if (atom->q_flag) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = type[j]; + buf[m++] = q[j]; + } + } else { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = type[j]; + } + } + + return m; +} + +// from MC/fix_atom_swap.cpp +void FixNSType::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + int *type = atom->type; + double *q = atom->q; + + m = 0; + last = first + n; + + if (atom->q_flag) { + for (i = first; i < last; i++) { + type[i] = static_cast(buf[m++]); + q[i] = buf[m++]; + } + } else { + for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); + } +} + /* ---------------------------------------------------------------------- */ void FixNSType::init() @@ -198,8 +252,8 @@ void FixNSType::initial_integrate(int vflag) dmuN = mu[new_type-1] - mu[type[atom_0]-1]; type[atom_0] = new_type; #ifdef DEBUG - std::cout << "SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << - type[atom_0] << " mu " << mu[type[atom_0]-1] << " "; + std::cout << "TYPE SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << + type[atom_0] << " mu " << mu[type[atom_0]-1] << std::endl; #endif } else { bool all_same = true; @@ -219,7 +273,7 @@ void FixNSType::initial_integrate(int vflag) #ifdef DEBUG #ifdef SWAP_POS - std::cout << "SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << + std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << " i " << atom_1 << " type " << type[atom_1] << " x " << x[atom_1][0] << " " << x[atom_1][1] << " " << x[atom_1][2] << " "; #else std::cout << "SWAP i " << atom_0 << " type " << type[atom_0] << " <-> i " << atom_1 << " type " << type[atom_1] << " "; @@ -246,6 +300,9 @@ void FixNSType::initial_integrate(int vflag) dmuN = 0.0; } + // communicate modified types to ghost atoms + comm->forward_comm(this); + // copied from fix_gmc.cpp // not sure if it's needed on rejections int id = modify->find_compute("thermo_pe"); @@ -260,6 +317,9 @@ void FixNSType::initial_integrate(int vflag) void FixNSType::final_integrate() { double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); +#ifdef DEBUG +std::cout << "TYPE ecurrent " << ecurrent << std::endl; +#endif // if potential energy + d(mu N) is above Emax then reject move // need to include previous steps' cumulative dmuN contributions, as well as current ones @@ -289,6 +349,9 @@ void FixNSType::final_integrate() } #endif + // communicate restored prev types to ghost atoms + comm->forward_comm(this); + // NOTE: no idea if this is necessary, or the right place to do this after rejecting a move // I suspect it's not necessary at all, since next step all these things will be done anyway if (unequal_cutoffs) { diff --git a/lammps_patches/new/src/NS/fix_ns_type.h b/lammps_patches/new/src/NS/fix_ns_type.h index 152031f..b328bf4 100644 --- a/lammps_patches/new/src/NS/fix_ns_type.h +++ b/lammps_patches/new/src/NS/fix_ns_type.h @@ -27,13 +27,15 @@ namespace LAMMPS_NS { class FixNSType : public Fix { public: FixNSType(class LAMMPS *, int, char **); - ~FixNSType(); - int setmask(); - virtual void init(); - virtual void pre_exchange(); - virtual void initial_integrate(int); - virtual void final_integrate(); - double compute_vector(int); + ~FixNSType() override; + int setmask() override; + void init() override; + void pre_exchange() override; + void initial_integrate(int) override; + void final_integrate() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + double compute_vector(int) override; protected: class Compute *pe_compute; @@ -46,7 +48,6 @@ class FixNSType : public Fix { double *mu; int peflag; char *id_pe; - char str[64]; double dmuN, cumulative_dmuN; bool unequal_cutoffs; diff --git a/pymatnext/__version__.py b/pymatnext/__version__.py index bbab024..1276d02 100644 --- a/pymatnext/__version__.py +++ b/pymatnext/__version__.py @@ -1 +1 @@ -__version__ = "0.1.4" +__version__ = "0.1.5" diff --git a/pymatnext/extras/lammps.patch b/pymatnext/extras/lammps.patch index 5a16288..8e0c472 100644 --- a/pymatnext/extras/lammps.patch +++ b/pymatnext/extras/lammps.patch @@ -1,7 +1,7 @@ diff -N -u -r orig/cmake/CMakeLists.txt new/cmake/CMakeLists.txt ---- orig/cmake/CMakeLists.txt 2023-07-20 17:31:29.420809479 -0400 -+++ new/cmake/CMakeLists.txt 2023-07-20 17:31:29.454809162 -0400 -@@ -305,7 +305,8 @@ +--- orig/cmake/CMakeLists.txt 2026-02-26 14:57:37.917018563 -0500 ++++ new/cmake/CMakeLists.txt 2026-02-26 14:57:37.959018179 -0500 +@@ -350,7 +350,8 @@ UEF VORONOI VTK @@ -13,8 +13,8 @@ diff -N -u -r orig/cmake/CMakeLists.txt new/cmake/CMakeLists.txt diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp --- orig/src/NS/fix_ns_cellmc.cpp 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_cellmc.cpp 2023-07-20 17:16:05.952402361 -0400 -@@ -0,0 +1,474 @@ ++++ new/src/NS/fix_ns_cellmc.cpp 2026-02-26 14:56:30.259636770 -0500 +@@ -0,0 +1,477 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories @@ -381,16 +381,1058 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp +{ + // if potential energy is above Emax then reject move + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "CELL ecurrent " << ecurrent << std::endl; ++#endif ++ ++ if (move_rejected_early) { ++ return; ++ } ++ ++ // if potential energy - d(P V) is above Emax then reject move ++ // need to include previous steps' cumulative dPV contributions, as well as current ones ++ if (ecurrent + cumulative_dPV + dPV >= Emax) { ++#ifdef DEBUG ++ std::cout << "REJECT E == " << ecurrent << " + " << cumulative_dPV << " + " << dPV << " >= Emax == " << Emax << std::endl; ++#endif ++ // reject move, so don't touch cumulative_dPV, since type change that led to current dPV was reverted ++ ++ for (int i=0; i < 3; i++) { ++ domain->boxlo[i] = prev_boxlo[i]; ++ domain->boxhi[i] = prev_boxhi[i]; ++ } ++ domain->xy = prev_xy; ++ domain->yz = prev_yz; ++ domain->xz = prev_xz; ++ ++ domain->set_global_box(); ++ domain->set_local_box(); ++ ++ double **x = atom->x; ++ int nlocal = atom->nlocal; ++ for (int iat=0; iat < nlocal; iat++) ++ for (int j=0; j < 3; j++) ++ x[iat][j] = prevx[iat][j]; ++ ++ } else { ++ // accept move, so accumulate dPV contribution from this step ++ cumulative_dPV += dPV; ++ switch (last_move_type) { ++ case MOVE_VOL: ++ n_success_vol++; ++ break; ++ case MOVE_STRETCH: ++ n_success_stretch++; ++ break; ++ case MOVE_SHEAR: ++ n_success_shear++; ++ break; ++ default: ++ error->all(FLERR,"Illegal value of last_move_type increment n_success_*"); ++ } ++#ifdef DEBUG ++double new_cell[3][3]; ++new_cell[0][0] = domain->boxhi[0] - domain->boxlo[0]; ++new_cell[0][1] = 0.0; ++new_cell[0][2] = 0.0; ++new_cell[1][0] = domain->xy; ++new_cell[1][1] = domain->boxhi[1] - domain->boxlo[1]; ++new_cell[1][2] = 0.0; ++new_cell[2][0] = domain->xz; ++new_cell[2][1] = domain->yz; ++new_cell[2][2] = domain->boxhi[2] - domain->boxlo[2]; ++ std::cout << "ACCEPT E == " << ecurrent << " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << min_aspect_ratio_val(new_cell) << std::endl; ++ // std::cout << "final cell " << new_cell[0][0] << " " << new_cell[0][1] << " " << new_cell[0][2] << std::endl; ++ // std::cout << " " << new_cell[1][0] << " " << new_cell[1][1] << " " << new_cell[1][2] << std::endl; ++ // std::cout << " " << new_cell[2][0] << " " << new_cell[2][1] << " " << new_cell[2][2] << std::endl; ++#endif ++ } ++} ++ ++FixNSCellMC::~FixNSCellMC() ++{ ++ ++ // delete temperature and pressure if fix created them ++ delete random; ++ ++ memory->destroy(prevx); ++ ++} ++ ++double FixNSCellMC::min_aspect_ratio_val(double cell[3][3]) { ++ double min_val = std::numeric_limits::max(); ++ for (int i=0; i < 3; i++) { ++ double vnorm_hat[3]; ++ MathExtra::cross3(cell[(i+1)%3], cell[(i+2)%3], vnorm_hat); ++ MathExtra::norm3(vnorm_hat); ++ double val = fabs(MathExtra::dot3(vnorm_hat, cell[i])); ++ min_val = (val < min_val) ? val : min_val; ++ } ++ ++ double V; ++ if (domain->dimension == 3) V = domain->xprd * domain->yprd * domain->zprd; ++ else V = domain->xprd * domain->yprd; ++ ++ return min_val / pow(V, 1.0/3.0); ++} ++ ++/* ---------------------------------------------------------------------- ++ return acceptance numbers ++------------------------------------------------------------------------- */ ++ ++double FixNSCellMC::compute_vector(int n) ++{ ++ switch (n) { ++ case 0: return n_attempt_vol; ++ case 1: return n_success_vol; ++ case 2: return n_attempt_stretch; ++ case 3: return n_success_stretch; ++ case 4: return n_attempt_shear; ++ case 5: return n_success_shear; ++ } ++ return -1.0; ++} +diff -N -u -r orig/src/NS/fix_ns_cellmc.h new/src/NS/fix_ns_cellmc.h +--- orig/src/NS/fix_ns_cellmc.h 1969-12-31 19:00:00.000000000 -0500 ++++ new/src/NS/fix_ns_cellmc.h 2026-02-26 14:56:30.259636770 -0500 +@@ -0,0 +1,73 @@ ++/* -*- c++ -*- ---------------------------------------------------------- ++ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator ++ http://lammps.sandia.gov, Sandia National Laboratories ++ Steve Plimpton, sjplimp@sandia.gov ++ ++ Copyright (2003) Sandia Corporation. Under the terms of Contract ++ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains ++ certain rights in this software. This software is distributed under ++ the GNU General Public License. ++ ++ See the README file in the top-level LAMMPS directory. ++------------------------------------------------------------------------- */ ++ ++#ifdef FIX_CLASS ++ ++FixStyle(ns/cellmc,FixNSCellMC) ++ ++#else ++ ++#ifndef LMP_FIX_NS_CELLMC_H ++#define LMP_FIX_NS_CELLMC_H ++ ++#include "fix.h" ++ ++namespace LAMMPS_NS { ++ ++class FixNSCellMC : public Fix { ++ public: ++ FixNSCellMC(class LAMMPS *, int, char **); ++ ~FixNSCellMC() override; ++ int setmask() override; ++ void init() override; ++ void initial_integrate(int) override; ++ void final_integrate() override; ++ double compute_vector(int) override; ++ ++ protected: ++ class Compute *pe_compute; ++ double **prevx; ++ double prev_boxlo[3], prev_boxhi[3], prev_xy, prev_yz, prev_xz; ++ int seed; ++ double Emax; ++ double min_aspect_ratio, pressure; ++ int flat_V_prior; ++ class RanMars *random; ++ double pVol, dVol, pStretch, dStretch, pShear, dShear; ++ int peflag; ++ char *id_pe; ++ ++ double dPV, cumulative_dPV; ++ ++ bool move_rejected_early; ++ int last_move_type; ++ int n_attempt_vol, n_attempt_stretch, n_attempt_shear, n_success_vol, n_success_stretch, n_success_shear; ++ ++ private: ++ double min_aspect_ratio_val(double cell[3][3]); ++}; ++ ++} ++ ++#endif ++#endif ++ ++/* ERROR/WARNING messages: ++ ++E: Illegal ... command ++ ++Self-explanatory. Check the input script syntax and compare to the ++documentation for the command. You can use -echo screen as a ++command-line option when running LAMMPS to see the offending line. ++ ++*/ +diff -N -u -r orig/src/NS/fix_ns.cpp new/src/NS/fix_ns.cpp +--- orig/src/NS/fix_ns.cpp 1969-12-31 19:00:00.000000000 -0500 ++++ new/src/NS/fix_ns.cpp 2026-02-26 14:56:30.260636761 -0500 +@@ -0,0 +1,1021 @@ ++/* ---------------------------------------------------------------------- ++ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator ++ http://lammps.sandia.gov, Sandia National Laboratories ++ Steve Plimpton, sjplimp@sandia.gov ++ ++ Copyright (2003) Sandia Corporation. Under the terms of Contract ++ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains ++ certain rights in this software. This software is distributed under ++ the GNU General Public License. ++ ++ See the README file in the top-level LAMMPS directory. ++------------------------------------------------------------------------- */ ++ ++#ifdef DEBUG ++#include ++#endif ++ ++#include ++#include ++#include "fix_ns.h" ++#include "atom.h" ++#include "comm.h" ++#include "force.h" ++#include "pair.h" ++#include "modify.h" ++#include "compute.h" ++#include "update.h" ++#include "random_mars.h" ++#include "domain.h" ++#include "memory.h" ++#include "respa.h" ++#include "error.h" ++#include "math_extra.h" ++#include ++ ++using namespace LAMMPS_NS; ++using namespace FixConst; ++ ++#define MOVE_UNDEF -1 ++#define MOVE_NONE 0 ++#define MOVE_POS 1 ++#define MOVE_CELL 2 ++#define MOVE_TYPE 3 ++#define MOVE_VOL 20 ++#define MOVE_STRETCH 21 ++#define MOVE_SHEAR 22 ++ ++/* ---------------------------------------------------------------------- */ ++ ++void normalize_cumul_probs_3(double &p1, double &p2, double &p3) { ++ p1 = (p1 < 0.0) ? 0.0 : p1; ++ p2 = (p2 < 0.0) ? 0.0 : p2; ++ p3 = (p3 < 0.0) ? 0.0 : p3; ++ double pSum = p1 + p2 + p3; ++ p1 /= pSum; p2 /= pSum; p3 /= pSum; ++ ++ // make into cumulative probabilities ++ p2 += p1; ++ if (p3 == 0.0) { ++ p2 = 1.0; ++ } ++ p3 = 1.0; ++} ++ ++FixNS::FixNS(LAMMPS *lmp, int narg, char **arg) : ++ Fix(lmp, narg, arg), dx(NULL), prevx(NULL), prevtype(NULL), mu(NULL) ++{ ++ // possible values: base 3 + global 2 + probabilities 3 + 1? + 9? + (1 + ntype)? ++ if (strcmp(style,"ns") != 0 && narg <= 3 + 2 + 3) ++ error->all(FLERR,"Illegal fix ns command"); ++ ++ // before ns-specific params ++ int iarg = 3; ++ ++ // global NS params ++ seed = utils::inumeric(FLERR,arg[iarg++],true,lmp); ++ Emax = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ ++ // probabilities of move types ++ pPos = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ pCell = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ pType = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ ++ pos_n_steps = cell_n_steps = type_n_steps = 0; ++ ++ // pos GMC params ++ if (pPos > 0.0) { ++ if (narg - iarg < 2) ++ error->all(FLERR,"Illegal fix ns command when parsing pos GMC arguments"); ++ pos_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); ++ gmc_step_size = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ memory->create(dx,atom->nmax,3,"ns:dx"); ++ memory->create(prevx,atom->nmax,3,"ns:prevx"); ++ } ++ ++ // cell MC params ++ if (pCell > 0.0) { ++ if (narg - iarg < 3 + 6 + 1) ++ error->all(FLERR,"Illegal fix ns command when parsing cell MC arguments"); ++ ++ if (!domain->triclinic) ++ error->all(FLERR,"fix ns cell moves require triclinic box (for now)"); ++ ++ cell_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); ++ min_aspect_ratio = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ pressure = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ ++ pVol = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ dVol = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ pStretch = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ dStretch = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ pShear = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ dShear = utils::numeric(FLERR,arg[iarg++],true, lmp); ++ flat_V_prior = 1; ++ if (strcmp(arg[iarg],"no") == 0) flat_V_prior = 0; ++ else if (strcmp(arg[iarg],"yes") != 0) ++ error->all(FLERR,"Illegal fix ns flat_V_prior value"); ++ iarg++; ++ ++#ifdef DEBUG ++std::cout << "DEBUG initial vol/stretch/shear probs " << pVol << " " << pStretch << " " << pShear << std::endl; ++#endif ++ normalize_cumul_probs_3(pVol, pStretch, pShear); ++#ifdef DEBUG ++std::cout << "DEBUG normalized cumulative vol/stretch/shear probs " << pVol << " " << pStretch << " " << pShear << std::endl; ++#endif ++ ++ if (!prevx) ++ memory->create(prevx,atom->nmax,3,"ns:prevx"); ++ ++ // from fix_nh.cpp ++ // who knows what other things may be needed? ++ box_change |= (BOX_CHANGE_X | BOX_CHANGE_Y | BOX_CHANGE_Z | ++ BOX_CHANGE_YZ | BOX_CHANGE_XZ | BOX_CHANGE_XY); ++ no_change_box = 1; ++ } ++ ++ unequal_cutoffs = false; ++ // type MC params ++ if (pType > 0.0) { ++ if (narg - iarg < 2) ++ error->all(FLERR,"Illegal fix ns command when parsing type MC arguments"); ++ ++ type_n_steps = utils::inumeric(FLERR,arg[iarg++],true,lmp); ++ semi_GC_flag = 1; ++ if (strcmp(arg[iarg],"no") == 0) semi_GC_flag = 0; ++ else if (strcmp(arg[iarg],"yes") != 0) ++ error->all(FLERR,"Illegal fix ns/type semi_GC_flag value"); ++ iarg++; ++ ++ if (semi_GC_flag) { ++ if (narg - iarg < atom->ntypes) ++ error->all(FLERR,"Illegal fix ns command when parsing type MC mu arguments"); ++ ++ mu = new double[atom->ntypes]; ++ for (int i=0; i < atom->ntypes; i++) { ++ mu[i] = utils::numeric(FLERR,arg[iarg++],true,lmp); ++ } ++ } else { ++ mu = 0; ++ } ++ ++ // based on fix_atom_swap.cpp ++ int ntypes = atom->ntypes; ++ double **cutsq = force->pair->cutsq; ++ for (int itype = 1; itype <= ntypes; itype++) ++ for (int jtype = 1; jtype <= ntypes; jtype++) ++ for (int ktype = 1; ktype <= ntypes; ktype++) ++ if (cutsq[itype][ktype] != cutsq[jtype][ktype]) ++ unequal_cutoffs = true; ++ if (unequal_cutoffs) ++ error->warning(FLERR,"unequal cutoffs may cause a problem for atom type moves"); ++ ++ memory->create(prevtype,atom->nmax,"ns:prevtype"); ++#ifdef SWAP_POS ++ if (!prevx) ++ memory->create(prevx,atom->nmax,3,"ns:prevx"); ++#endif ++ } ++ ++#ifdef DEBUG ++std::cout << "DEBUG initial pos/cell/type probs " << pPos << " " << pCell << " " << pType << std::endl; ++#endif ++ normalize_cumul_probs_3(pPos, pCell, pType); ++#ifdef DEBUG ++std::cout << "DEBUG cumulative normalized pos/cell/type probs " << pPos << " " << pCell << " " << pType << std::endl; ++#endif ++ ++ max_n_steps = (pos_n_steps > cell_n_steps) ? pos_n_steps : cell_n_steps; ++ max_n_steps = (max_n_steps > type_n_steps) ? max_n_steps : type_n_steps; ++ ++ if (narg - iarg < 0) ++ error->all(FLERR,"Illegal fix ns command - extra arguments left after parsing"); ++ ++ cell_cur_move = MOVE_UNDEF; ++ state_cur_move = MOVE_UNDEF; ++ state_traj_steps_remaining = 0; ++ ++ // zero accumulated energy shifts ++ cumulative_dPV = cumulative_dmuN = 0.0; ++ // attempt frequency counters ++ vector_flag = 1; ++ size_vector = 2 + 6 + 2; ++ global_freq = 1; ++ extvector = 0; ++ n_attempt_pos = n_success_pos = ++ n_attempt_vol = n_success_vol = ++ n_attempt_stretch = n_success_stretch = ++ n_attempt_shear = n_success_shear = ++ n_attempt_type = n_success_type = 0; ++ ++ // random_g will be used when all processes need to agree on the random ++ // number, so it uses the same seed on all. ++ // random_l will be used when each processes needs to do something different, ++ // and it may be called an unpredictable number of times, which will make it ++ // diverge between the different processes. ++ random_g = new RanMars(lmp,seed); ++ random_l = new RanMars(lmp,seed + 1 + comm->me); ++ ++ // copied from fix_gmc.cpp ++ dynamic_group_allow = 1; ++ time_integrate = 1; ++ // copied from atom_swap.cpp? ++ // set force_reneighbor so next_reneighbor is checked, but don't force it ++ // for the next step until we see what sort of step this is - only type ++ // steps require (?) reneighbors ++ force_reneighbor = 1; ++ // next_reneighbor = update->ntimestep + 1; ++ next_reneighbor = -1; ++ ++ // from MC/fix_atom_swap.cpp ++ // amount of data to be packed for inter-process comm (atom type) ++ if (atom->q_flag) ++ comm_forward = 2; ++ else ++ comm_forward = 1; ++} ++ ++/* ---------------------------------------------------------------------- */ ++ ++int FixNS::setmask() ++{ ++ int mask = 0; ++ mask |= INITIAL_INTEGRATE; ++ mask |= PRE_EXCHANGE; ++ mask |= FINAL_INTEGRATE; ++ mask |= INITIAL_INTEGRATE_RESPA; ++ mask |= FINAL_INTEGRATE_RESPA; ++ return mask; ++} ++ ++// from MC/fix_atom_swap.cpp ++int FixNS::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) ++{ ++ int i, j, m; ++ ++ int *type = atom->type; ++ double *q = atom->q; ++ ++ m = 0; ++ ++ if (atom->q_flag) { ++ for (i = 0; i < n; i++) { ++ j = list[i]; ++ buf[m++] = type[j]; ++ buf[m++] = q[j]; ++ } ++ } else { ++ for (i = 0; i < n; i++) { ++ j = list[i]; ++ buf[m++] = type[j]; ++ } ++ } ++ ++ return m; ++} ++ ++// from MC/fix_atom_swap.cpp ++void FixNS::unpack_forward_comm(int n, int first, double *buf) ++{ ++ int i, m, last; ++ ++ int *type = atom->type; ++ double *q = atom->q; ++ ++ m = 0; ++ last = first + n; ++ ++ if (atom->q_flag) { ++ for (i = first; i < last; i++) { ++ type[i] = static_cast(buf[m++]); ++ q[i] = buf[m++]; ++ } ++ } else { ++ for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); ++ } ++} ++ ++void FixNS::init() ++{ ++ ++ int id = modify->find_compute("thermo_pe"); ++ ++ modify->compute[id]->invoked_scalar = -1; ++ pe_compute = modify->compute[id]; ++ pe_compute->addstep(update->ntimestep+1); ++ ++ if (strstr(update->integrate_style,"respa")) ++ error->all(FLERR,"fix ns not compatible with RESPA"); ++} ++ ++void FixNS::pre_exchange() ++{ ++ if (state_cur_move == MOVE_TYPE) { ++ // force reneighbor for swaps (maybe only needed on accept, ++ // or only with unequal cutoffs?) ++ next_reneighbor = update->ntimestep + 1; ++ } else { ++ next_reneighbor = -1; ++ } ++} ++ ++void FixNS::initial_integrate(int vflag) ++{ ++ // state machine for move types ++ if (state_traj_steps_remaining <= 0) { ++ // for now, give up on moving if there aren't enough steps left ++ // to accommodate longest (across 3 step types) mini-traj. This ++ // should ensure correct ratio of step types ++ if (update->laststep + 1 - update->ntimestep < max_n_steps) { ++#ifdef DEBUG ++std::cout << "initial_integrate not enough steps left " << update->laststep << " + 1 - " << update->ntimestep << " < " << max_n_steps << std::endl; ++#endif ++ // no time for another mini traj ++ state_cur_move = MOVE_NONE; ++ } else { ++ // pick new move type, all processes must agree ++ double rv = random_g->uniform(); ++#ifdef DEBUG ++std::cout << "initial_integrate pick rv " << rv << " probs " << pPos << " " << pCell << " " << pType << std::endl; ++#endif ++ if (rv < pPos) { ++ state_cur_move = MOVE_POS; ++ pos_gmc_traj_prep(); ++ state_traj_steps_remaining = pos_n_steps; ++ } else if (rv < pCell) { ++ state_cur_move = MOVE_CELL; ++ state_traj_steps_remaining = cell_n_steps; ++ } else { ++ state_cur_move = MOVE_TYPE; ++ state_traj_steps_remaining = type_n_steps; ++ } ++ } ++ } ++#ifdef DEBUG ++else ++std::cout << "initial_integrate continue with current move type" << std::endl; ++#endif ++ ++ state_traj_steps_remaining--; ++ ++ if (state_cur_move == MOVE_POS) { ++ pos_gmc_initial_integrate(); ++ } else if (state_cur_move == MOVE_CELL) { ++ cell_initial_integrate(); ++ } else if (state_cur_move == MOVE_TYPE) { ++ type_initial_integrate(); ++ // move swapped or perturbed types, communicate to ghost atoms ++ comm->forward_comm(this); ++ } ++ ++ int id = modify->find_compute("thermo_pe"); ++ ++ modify->compute[id]->invoked_scalar = -1; ++ pe_compute = modify->compute[id]; ++ pe_compute->addstep(update->ntimestep+1); ++} ++ ++void FixNS::pos_gmc_traj_prep() ++{ ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ double **x = atom->x; ++ ++ double dx2sum = 0; ++ ++ // if we ever want to maintain coherence of dx from traj to traj, we can do it ++ // using atomic velocities: ++ // double **v = atom->v; ++ ++ for(int i = 0; i < nlocal; i++) { ++ // random direction ++ dx[i][0] = random_l->gaussian(); ++ dx[i][1] = random_l->gaussian(); ++ dx[i][2] = random_l->gaussian(); ++ dx2sum += dx[i][0]*dx[i][0]+dx[i][1]*dx[i][1]+dx[i][2]*dx[i][2]; ++ ++ // save initial pos in case we have to reject move ++ prevx[i][0] = x[i][0]; ++ prevx[i][1] = x[i][1]; ++ prevx[i][2] = x[i][2]; ++ } ++ ++ // WARNING this really needs to be summed over all MPI processes ++ dx2sum = sqrt(dx2sum); ++ for(int i = 0; i < nlocal; i++) { ++ dx[i][0] /= dx2sum; ++ dx[i][1] /= dx2sum; ++ dx[i][2] /= dx2sum; ++ } ++ ++ n_attempt_pos += 1; ++} ++ ++void FixNS::pos_gmc_initial_integrate() ++{ ++#ifdef DEBUG ++std::cout << "POS pos_gmc_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; ++#endif ++ // update x of atoms in group ++ double **x = atom->x; ++ int *mask = atom->mask; ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ ++ for (int i = 0; i < nlocal; i++) ++ if (mask[i] & groupbit) { ++ x[i][0] += gmc_step_size * dx[i][0]; ++ x[i][1] += gmc_step_size * dx[i][1]; ++ x[i][2] += gmc_step_size * dx[i][2]; ++ } ++ ++} ++ ++void FixNS::cell_initial_integrate() ++{ ++#ifdef DEBUG ++std::cout << "CELL cell_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; ++#endif ++ // do move (unles rejected) ++ ++ int natoms = atom->natoms; ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ double **x = atom->x; ++ ++ // save previous cell and positions ++ for (int i=0; i < 3; i++) { ++ prev_boxlo[i] = domain->boxlo[i]; ++ prev_boxhi[i] = domain->boxhi[i]; ++ } ++ prev_xy = domain->xy; ++ prev_yz = domain->yz; ++ prev_xz = domain->xz; ++ ++ for (int iat=0; iat < nlocal; iat++) ++ for (int j=0; j < 3; j++) ++ prevx[iat][j] = x[iat][j]; ++ ++ double boxext[3]; ++ for (int i=0; i < 3; i++) ++ boxext[i] = domain->boxhi[i] - domain->boxlo[i]; ++ ++ cell_cur_move = MOVE_UNDEF; ++ cell_move_rejected_early = false; ++ ++ // keep track of change PV in this step ++ // default to 0, set otherwise when volume move is selected ++ dPV = 0.0; ++ ++ // default to accept, will be changed below to reject (-1.0) or probablistic (if flat_V_prior != 1) ++ double rv = random_g->uniform(); ++ double new_cell[3][3]; ++ // pick a step type ++ if (rv < pVol) { ++ n_attempt_vol++; ++ cell_cur_move = MOVE_VOL; ++ // volume step ++ double orig_V; ++ if (domain->dimension == 3) orig_V = domain->xprd * domain->yprd * domain->zprd; ++ else orig_V = domain->xprd * domain->yprd; ++ ++ double dV = random_g->gaussian(0.0, dVol); ++ double new_V = orig_V + dV; ++#ifdef DEBUG ++ std::cout << "CELL VOLUME volumetric strain " << dV << " / " << orig_V << " = " << dV/orig_V << " "; ++#endif ++ if (new_V/orig_V < 0.5) { ++ cell_move_rejected_early = true; ++#ifdef DEBUG ++ std::cout << "CELL REJECT new_V/orig_V < 0.5" << std::endl; ++#endif ++ } else { ++ if (flat_V_prior == 0 && new_V < orig_V && random_g->uniform() > pow(new_V / orig_V, natoms)) { ++ cell_move_rejected_early = true; ++#ifdef DEBUG ++ std::cout << "CELL REJECT V probability " << pow(new_V / orig_V, natoms) << std::endl; ++#endif ++ } else { ++ double transform_diag = pow(new_V / orig_V, 1.0/3.0); ++ new_cell[0][0] = boxext[0] * transform_diag; ++ new_cell[0][1] = 0.0; ++ new_cell[0][2] = 0.0; ++ new_cell[1][0] = domain->xy * transform_diag; ++ new_cell[1][1] = boxext[1] * transform_diag; ++ new_cell[1][2] = 0.0; ++ new_cell[2][0] = domain->xz * transform_diag; ++ new_cell[2][1] = domain->yz * transform_diag; ++ new_cell[2][2] = boxext[2] * transform_diag; ++ dPV = pressure * dV; ++ } ++ } ++ } else if (rv < pStretch) { ++ n_attempt_stretch++; ++ cell_cur_move = MOVE_STRETCH; ++ // stretch step ++ // pick directions to use v_ind and (v_ind+1)%3 ++#ifdef DEBUG ++ std::cout << "CELL STRETCH "; ++#endif ++ int v_ind = int(3 * random_g->uniform()) % 3; ++ rv = random_g->gaussian(0.0, dStretch); ++ double transform_diag[3]; ++ transform_diag[v_ind] = exp(rv); ++ transform_diag[(v_ind+1)%3] = exp(-rv); ++ transform_diag[(v_ind+2)%3] = 1.0; ++#ifdef DEBUG ++ std::cout << transform_diag[0] << " " << transform_diag[1] << " " << transform_diag[2] << " "; ++#endif ++ // create new cell for aspect ratio test and new domain ++ new_cell[0][0] = boxext[0] * transform_diag[0]; ++ new_cell[0][1] = 0.0; ++ new_cell[0][2] = 0.0; ++ new_cell[1][0] = domain->xy * transform_diag[0]; ++ new_cell[1][1] = boxext[1] * transform_diag[1]; ++ new_cell[1][2] = 0.0; ++ new_cell[2][0] = domain->xz * transform_diag[0]; ++ new_cell[2][1] = domain->yz * transform_diag[1]; ++ new_cell[2][2] = boxext[2] * transform_diag[2]; ++ if (cur_aspect_ratio(new_cell) < min_aspect_ratio) { ++ cell_move_rejected_early = true; ++#ifdef DEBUG ++ std::cout << "CELL REJECT min_aspect_ratio " << cur_aspect_ratio(new_cell) << std::endl; ++#endif ++ } ++ } else { ++ n_attempt_shear++; ++ cell_cur_move = MOVE_SHEAR; ++ // shear step ++ // save original cell and a temporary cell ++#ifdef DEBUG ++ std::cout << "CELL SHEAR "; ++#endif ++ double orig_cell[3][3], t_cell[3][3]; ++ t_cell[0][0] = orig_cell[0][0] = boxext[0]; ++ t_cell[0][1] = orig_cell[0][1] = 0.0; ++ t_cell[0][2] = orig_cell[0][2] = 0.0; ++ t_cell[1][0] = orig_cell[1][0] = domain->xy; ++ t_cell[1][1] = orig_cell[1][1] = boxext[1]; ++ t_cell[1][2] = orig_cell[1][2] = 0.0; ++ t_cell[2][0] = orig_cell[2][0] = domain->xz; ++ t_cell[2][1] = orig_cell[2][1] = domain->yz; ++ t_cell[2][2] = orig_cell[2][2] = boxext[2]; ++ ++ // pick vector to perturb ++ int vec_ind = int(3 * random_g->uniform()) % 3; ++ ++ // perturb t_cell[vec_ind] along other two vectors ++ double vhat[3]; ++ for (int di=1; di < 3; di++) { ++ for (int i=0; i < 3; i++) ++ vhat[i] = orig_cell[(vec_ind + di) % 3][i]; ++ MathExtra::norm3(vhat); ++ rv = random_g->gaussian(0.0, dShear); ++#ifdef DEBUG ++ std::cout << (vec_ind + di) % 3 << " " << rv << " "; ++#endif ++ for (int i=0; i < 3; i++) ++ t_cell[vec_ind][i] += rv * vhat[i]; ++ } ++ ++ if (cur_aspect_ratio(t_cell) < min_aspect_ratio) { ++ cell_move_rejected_early = true; ++#ifdef DEBUG ++ std::cout << "CELL REJECT min_aspect_ratio " << cur_aspect_ratio(t_cell) << std::endl; ++#endif ++ } else { ++ // rotate new_cell back to LAMMPS orientation ++ ++ double a0_norm = sqrt(MathExtra::dot3(t_cell[0], t_cell[0])); ++ new_cell[0][0] = a0_norm; ++ new_cell[0][1] = 0.0; ++ new_cell[0][2] = 0.0; ++ ++ double a0_hat[3] = {t_cell[0][0], t_cell[0][1], t_cell[0][2]}; MathExtra::norm3(a0_hat); ++ new_cell[1][0] = MathExtra::dot3(t_cell[1], a0_hat); ++ double a0_hat_cross_a1[3]; MathExtra::cross3(a0_hat, t_cell[1], a0_hat_cross_a1); ++ new_cell[1][1] = sqrt(MathExtra::dot3(a0_hat_cross_a1, a0_hat_cross_a1)); ++ new_cell[1][2] = 0.0; ++ ++ double a0_cross_a1_hat[3] = {a0_hat_cross_a1[0], a0_hat_cross_a1[1], a0_hat_cross_a1[2]}; ++ MathExtra::norm3(a0_cross_a1_hat); ++ ++ double a0_cross_a1_hat_cross_a0_hat[3]; ++ MathExtra::cross3(a0_cross_a1_hat, a0_hat, a0_cross_a1_hat_cross_a0_hat); ++ new_cell[2][0] = MathExtra::dot3(t_cell[2], a0_hat); ++ new_cell[2][1] = MathExtra::dot3(t_cell[2], a0_cross_a1_hat_cross_a0_hat); ++ new_cell[2][2] = MathExtra::dot3(t_cell[2], a0_cross_a1_hat); ++ } ++ } ++ ++ if (!cell_move_rejected_early) { ++ // apply move to box ++ double boxctr[3]; ++ for (int i=0; i < 3; i++) ++ boxctr[i] = 0.5*(prev_boxhi[0] + prev_boxlo[0]); ++ ++ domain->x2lamda(nlocal); ++ ++ boxext[0] = new_cell[0][0]; ++ domain->xy = new_cell[1][0]; ++ boxext[1] = new_cell[1][1]; ++ domain->xz = new_cell[2][0]; ++ domain->yz = new_cell[2][1]; ++ boxext[2] = new_cell[2][2]; ++ ++ for (int i=0; i < 3; i++) { ++ domain->boxlo[i] = boxctr[i] - boxext[i]/2.0; ++ domain->boxhi[i] = boxctr[i] + boxext[i]/2.0; ++ } ++ ++ domain->set_global_box(); ++ domain->set_local_box(); ++ ++ domain->lamda2x(nlocal); ++#ifdef DEBUG ++ std::cout << std::endl; ++#endif ++ } ++ ++} ++ ++void FixNS::type_initial_integrate() ++{ ++#ifdef DEBUG ++std::cout << "TYPE type_initial_integrate ecurrent " << modify->compute[modify->find_compute("thermo_pe")]->compute_scalar() << std::endl; ++#endif ++ // do move (unless rejected) ++ ++ n_attempt_type++; ++ ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ int ntypes = atom->ntypes; ++ ++ int *type = atom->type; ++ for (int iat=0; iat < nlocal; iat++) ++ prevtype[iat] = type[iat]; ++ ++#ifdef SWAP_POS ++ double **x = atom->x; ++ for (int iat=0; iat < nlocal; iat++) ++ for (int jat=0; jat < 3; jat++) ++ prevx[iat][jat] = x[iat][jat]; ++#endif ++ ++ // WARNING move is restricted to only one proc - fine for semi-GC, but ++ // limiting for fixed composition swaps ++ double rv = random_g->uniform(); ++ int rnd_proc = static_cast(comm->nprocs * rv) % comm->nprocs; ++ if (comm->me != rnd_proc) ++ return; ++ ++ rv = random_l->uniform(); ++ int atom_0 = static_cast(nlocal * rv) % nlocal; ++ ++ if (semi_GC_flag) { ++ // perturb type randomly, store change to mu*N term ++ rv = random_l->uniform(); ++ int dtype = 1 + static_cast(rv * (ntypes-1)) % (ntypes-1); ++ int new_type = 1 + (((type[atom_0]-1) + dtype) % ntypes); ++ dmuN = mu[new_type-1] - mu[type[atom_0]-1]; ++ type[atom_0] = new_type; ++#ifdef DEBUG ++ std::cout << "TYPE SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << ++ type[atom_0] << " mu " << mu[type[atom_0]-1] << " " << std::endl; ++#endif ++ } else { ++ // not semi-GC, swap a pair of atoms ++ bool all_same = true; ++ for (int i=1; i < nlocal; i++) ++ if (type[0] != type[i]) { ++ all_same = false; ++ break; ++ } ++ if (all_same) ++ return; ++ ++ int atom_1 = atom_0; ++ while (type[atom_0] == type[atom_1]) { ++ rv = random_l->uniform(); ++ atom_1 = static_cast(nlocal * rv) % nlocal; ++ } ++ ++#ifdef DEBUG ++#ifdef SWAP_POS ++ std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << ++ " i " << atom_1 << " type " << type[atom_1] << " x " << x[atom_1][0] << " " << x[atom_1][1] << " " << x[atom_1][2] << " " << std::endl; ++#else ++ std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " <-> i " << atom_1 << " type " << type[atom_1] << " " << std::endl; ++#endif ++#endif ++ ++#ifdef SWAP_POS ++ double tx[3]; ++ tx[0] = x[atom_0][0]; ++ tx[1] = x[atom_0][1]; ++ tx[2] = x[atom_0][2]; ++ x[atom_0][0] = x[atom_1][0]; ++ x[atom_0][1] = x[atom_1][1]; ++ x[atom_0][2] = x[atom_1][2]; ++ x[atom_1][0] = tx[0]; ++ x[atom_1][1] = tx[1]; ++ x[atom_1][2] = tx[2]; ++#else ++ int t_type = type[atom_0]; ++ type[atom_0] = type[atom_1]; ++ type[atom_1] = t_type; ++#endif ++ ++ dmuN = 0.0; ++ } ++} ++ ++/* ---------------------------------------------------------------------- */ ++ ++void FixNS::final_integrate() ++{ ++ if (state_cur_move == MOVE_POS) { ++ pos_gmc_final_integrate(); ++ } else if (state_cur_move == MOVE_CELL) { ++ cell_final_integrate(); ++ } else if (state_cur_move == MOVE_TYPE) { ++ type_final_integrate(); ++ } ++} ++ ++void FixNS::pos_gmc_final_integrate() ++{ ++ // if potential energy is above Emax then want to modify dx with ++ // forces to change trajectory ++ ++ double **f = atom->f; ++ double **x = atom->x; ++ int *mask = atom->mask; ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ ++ double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "POS ecurrent " << ecurrent << " + " << cumulative_dPV << " - " << cumulative_dmuN; ++#endif ++ ecurrent += cumulative_dPV - cumulative_dmuN; ++#ifdef DEBUG ++std::cout << " = " << ecurrent << std::endl; ++#endif ++ ++ if (ecurrent < Emax) { ++#ifdef DEBUG ++std::cout << "POS ecurrent + dPV - dmuN = " << ecurrent << " < " << Emax << std::endl; ++#endif ++ // accept ++ if (state_traj_steps_remaining == 0) { ++#ifdef DEBUG ++std::cout << "POS ACCEPT end of traj, increment n_success_pos" << std::endl; ++#endif ++ // final step of mini traj, accept ++ n_success_pos += 1; ++ } ++#ifdef DEBUG ++std::cout << "POS mid traj, continuing" << std::endl; ++#endif ++ // accept new position and continue straight ++ return; ++ } ++ ++#ifdef DEBUG ++std::cout << "POS ecurrent + dPV - dmuN = " << ecurrent << " >= " << Emax << std::endl; ++#endif ++ ++ if (state_traj_steps_remaining == 0) { ++#ifdef DEBUG ++std::cout << "POS REJECT end of traj, revert positions to prevx" << std::endl; ++#endif ++ // final step of mini traj, reject ++ for (int i = 0; i < nlocal; i++) { ++ x[i][0] = prevx[i][0]; ++ x[i][1] = prevx[i][1]; ++ x[i][2] = prevx[i][2]; ++ } ++ return; ++ } ++ ++#ifdef DEBUG ++std::cout << "POS reflect" << std::endl; ++#endif ++ ++ // try to reflect from boundary ++ double fsum = 0; ++ double fhatdotdx = 0; ++ for (int i = 0; i < nlocal; i++) ++ if (mask[i] & groupbit) ++ fsum += f[i][0]*f[i][0]+ ++ f[i][1]*f[i][1]+ ++ f[i][2]*f[i][2]; ++ fsum = sqrt(fsum); ++ ++ // should also detect nan ++ if (fsum != 0.0) { ++ for (int i = 0; i < nlocal; i++) { ++ if (mask[i] & groupbit) { ++ fhatdotdx += f[i][0]/fsum*dx[i][0]; ++ fhatdotdx += f[i][1]/fsum*dx[i][1]; ++ fhatdotdx += f[i][2]/fsum*dx[i][2]; ++ } ++ } ++ ++ for (int i = 0; i < nlocal; i++) { ++ if (mask[i] & groupbit) { ++ dx[i][0] -= 2*f[i][0]/fsum*fhatdotdx; ++ dx[i][1] -= 2*f[i][1]/fsum*fhatdotdx; ++ dx[i][2] -= 2*f[i][2]/fsum*fhatdotdx; ++ } ++ } ++ } + -+ if (move_rejected_early) { ++} ++ ++ ++void FixNS::cell_final_integrate() ++{ ++ // if potential energy is above Emax then reject move ++ double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "CELL ecurrent " << ecurrent << " - " << cumulative_dmuN; ++#endif ++ ecurrent -= cumulative_dmuN; ++#ifdef DEBUG ++std::cout << " = " << ecurrent << std::endl; ++#endif ++ ++ if (cell_move_rejected_early) { + return; + } + -+ // if potential energy - d(P V) is above Emax then reject move ++ // if potential energy + d(P V) is above Emax then reject move + // need to include previous steps' cumulative dPV contributions, as well as current ones + if (ecurrent + cumulative_dPV + dPV >= Emax) { +#ifdef DEBUG -+ std::cout << "REJECT E == " << ecurrent << " + " << cumulative_dPV + dPV << " >= Emax == " << Emax << std::endl; ++ std::cout << "CELL REJECT E == " << ecurrent << " + " << cumulative_dPV << " + " << dPV << " >= Emax == " << Emax << std::endl; +#endif + // reject move, so don't touch cumulative_dPV, since type change that led to current dPV was reverted + @@ -407,6 +1449,7 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp + + double **x = atom->x; + int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; + for (int iat=0; iat < nlocal; iat++) + for (int j=0; j < 3; j++) + x[iat][j] = prevx[iat][j]; @@ -414,7 +1457,7 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp + } else { + // accept move, so accumulate dPV contribution from this step + cumulative_dPV += dPV; -+ switch (last_move_type) { ++ switch (cell_cur_move) { + case MOVE_VOL: + n_success_vol++; + break; @@ -425,7 +1468,7 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp + n_success_shear++; + break; + default: -+ error->all(FLERR,"Illegal value of last_move_type increment n_success_*"); ++ error->all(FLERR,"Illegal value of cell_cur_move in increment n_success_*"); + } +#ifdef DEBUG +double new_cell[3][3]; @@ -438,7 +1481,7 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp +new_cell[2][0] = domain->xz; +new_cell[2][1] = domain->yz; +new_cell[2][2] = domain->boxhi[2] - domain->boxlo[2]; -+ std::cout << "ACCEPT E == " << ecurrent " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << min_aspect_ratio_val(new_cell) << std::endl; ++ std::cout << "CELL ACCEPT E == " << ecurrent << " + " << cumulative_dPV << " < Emax == " << Emax << " min_aspect " << cur_aspect_ratio(new_cell) << std::endl; + // std::cout << "final cell " << new_cell[0][0] << " " << new_cell[0][1] << " " << new_cell[0][2] << std::endl; + // std::cout << " " << new_cell[1][0] << " " << new_cell[1][1] << " " << new_cell[1][2] << std::endl; + // std::cout << " " << new_cell[2][0] << " " << new_cell[2][1] << " " << new_cell[2][2] << std::endl; @@ -446,131 +1489,115 @@ diff -N -u -r orig/src/NS/fix_ns_cellmc.cpp new/src/NS/fix_ns_cellmc.cpp + } +} + -+FixNSCellMC::~FixNSCellMC() -+{ -+ -+ // delete temperature and pressure if fix created them -+ delete random; + -+ memory->destroy(prevx); ++void FixNS::type_final_integrate() ++{ ++ double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "TYPE ecurrent " << ecurrent << " + " << cumulative_dPV; ++#endif ++ ecurrent += cumulative_dPV; ++#ifdef DEBUG ++std::cout << " = " << ecurrent << std::endl; ++#endif + -+} ++ // if potential energy - d(mu N) is above Emax then reject move ++ // need to include previous steps' cumulative dmuN contributions, as well as current ones ++ if (ecurrent - cumulative_dmuN - dmuN >= Emax) { ++#ifdef DEBUG ++ std::cout << "TYPE REJECT E == " << ecurrent << " - " << cumulative_dmuN << " - " << dmuN << " >= Emax == " << Emax << " " << std::endl; ++#endif ++ // reject move, so don't touch cumulative_dmuN, since type change that led to current dmuN was reverted + -+double FixNSCellMC::min_aspect_ratio_val(double cell[3][3]) { -+ double min_val = std::numeric_limits::max(); -+ for (int i=0; i < 3; i++) { -+ double vnorm_hat[3]; -+ MathExtra::cross3(cell[(i+1)%3], cell[(i+2)%3], vnorm_hat); -+ MathExtra::norm3(vnorm_hat); -+ double val = fabs(MathExtra::dot3(vnorm_hat, cell[i])); -+ min_val = (val < min_val) ? val : min_val; ++ int nlocal = atom->nlocal; ++ if (igroup == atom->firstgroup) nlocal = atom->nfirst; ++ int *type = atom->type; ++// With SWAP_POS, type is only modified for semi-GC, so only has to be undone in that case, and x ++// has to be reverted to undo swap move. ++// By default (no SWAP_POS), type is modified for both swap and semi-GC, so undoing either step type ++// uses prevtype ++#ifdef SWAP_POS ++ if (semi_GC_flag) { ++#endif ++ for (int iat=0; iat < nlocal; iat++) ++ type[iat] = prevtype[iat]; ++#ifdef SWAP_POS ++ } else { ++ double **x = atom->x; ++ for (int iat=0; iat < nlocal; iat++) ++ for (int jat=0; jat < 3; jat++) ++ x[iat][jat] = prevx[iat][jat]; + } ++#endif + -+ double V; -+ if (domain->dimension == 3) V = domain->xprd * domain->yprd * domain->zprd; -+ else V = domain->xprd * domain->yprd; -+ -+ return min_val / pow(V, 1.0/3.0); -+} -+ -+/* ---------------------------------------------------------------------- -+ return acceptance numbers -+------------------------------------------------------------------------- */ ++ // modified type by restoring previous, communicate to ghost atoms ++ if (state_cur_move == MOVE_TYPE) { ++ comm->forward_comm(this); ++ } + -+double FixNSCellMC::compute_vector(int n) -+{ -+ switch (n) { -+ case 0: return n_attempt_vol; -+ case 1: return n_success_vol; -+ case 2: return n_attempt_stretch; -+ case 3: return n_success_stretch; -+ case 4: return n_attempt_shear; -+ case 5: return n_success_shear; ++ } else { ++ // accept move, so accumulate dmuN contribution from this step ++ cumulative_dmuN += dmuN; ++ n_success_type++; ++#ifdef DEBUG ++ std::cout << "TYPE ACCEPT E == " << ecurrent << " - " << cumulative_dmuN << " < Emax == " << Emax << std::endl; ++#endif + } -+ return -1.0; +} -diff -N -u -r orig/src/NS/fix_ns_cellmc.h new/src/NS/fix_ns_cellmc.h ---- orig/src/NS/fix_ns_cellmc.h 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_cellmc.h 2023-07-20 17:16:05.953402352 -0400 -@@ -0,0 +1,74 @@ -+/* -*- c++ -*- ---------------------------------------------------------- -+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator -+ http://lammps.sandia.gov, Sandia National Laboratories -+ Steve Plimpton, sjplimp@sandia.gov -+ -+ Copyright (2003) Sandia Corporation. Under the terms of Contract -+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains -+ certain rights in this software. This software is distributed under -+ the GNU General Public License. -+ -+ See the README file in the top-level LAMMPS directory. -+------------------------------------------------------------------------- */ -+ -+#ifdef FIX_CLASS -+ -+FixStyle(ns/cellmc,FixNSCellMC) -+ -+#else -+ -+#ifndef LMP_FIX_NS_CELLMC_H -+#define LMP_FIX_NS_CELLMC_H -+ -+#include "fix.h" -+ -+namespace LAMMPS_NS { -+ -+class FixNSCellMC : public Fix { -+ public: -+ FixNSCellMC(class LAMMPS *, int, char **); -+ ~FixNSCellMC(); -+ int setmask(); -+ virtual void init(); -+ virtual void initial_integrate(int); -+ virtual void final_integrate(); -+ double compute_vector(int); -+ -+ protected: -+ class Compute *pe_compute; -+ double **prevx; -+ double prev_boxlo[3], prev_boxhi[3], prev_xy, prev_yz, prev_xz; -+ int seed; -+ double Emax; -+ double min_aspect_ratio, pressure; -+ int flat_V_prior; -+ class RanMars *random; -+ double pVol, dVol, pStretch, dStretch, pShear, dShear; -+ int peflag; -+ char *id_pe; -+ char str[64]; -+ -+ double dPV, cumulative_dPV; -+ -+ bool move_rejected_early; -+ int last_move_type; -+ int n_attempt_vol, n_attempt_stretch, n_attempt_shear, n_success_vol, n_success_stretch, n_success_shear; + -+ private: -+ double min_aspect_ratio_val(double cell[3][3]); -+}; ++FixNS::~FixNS() ++{ ++ delete random_g; ++ delete random_l; + ++ if (dx) ++ memory->destroy(dx); ++ if (prevx) ++ memory->destroy(prevx); ++ if (prevtype) ++ memory->destroy(prevtype); ++ if (mu) ++ delete mu; +} + -+#endif -+#endif -+ -+/* ERROR/WARNING messages: ++double FixNS::cur_aspect_ratio(double cell[3][3]) ++{ ++ double min_val = std::numeric_limits::max(); ++ for (int i=0; i < 3; i++) { ++ double vnorm_hat[3]; ++ MathExtra::cross3(cell[(i+1)%3], cell[(i+2)%3], vnorm_hat); ++ MathExtra::norm3(vnorm_hat); ++ double val = fabs(MathExtra::dot3(vnorm_hat, cell[i])); ++ min_val = (val < min_val) ? val : min_val; ++ } + -+E: Illegal ... command ++ double V; ++ if (domain->dimension == 3) V = domain->xprd * domain->yprd * domain->zprd; ++ else V = domain->xprd * domain->yprd; + -+Self-explanatory. Check the input script syntax and compare to the -+documentation for the command. You can use -echo screen as a -+command-line option when running LAMMPS to see the offending line. ++ return min_val / pow(V, 1.0/3.0); ++} + -+*/ ++double FixNS::compute_vector(int n) ++{ ++ switch (n) { ++ case 0: return n_attempt_pos; ++ case 1: return n_success_pos; ++ case 2: return n_attempt_vol; ++ case 3: return n_success_vol; ++ case 4: return n_attempt_stretch; ++ case 5: return n_success_stretch; ++ case 6: return n_attempt_shear; ++ case 7: return n_success_shear; ++ case 8: return n_attempt_type; ++ case 9: return n_success_type; ++ } ++ return -1.0; ++} diff -N -u -r orig/src/NS/fix_ns_gmc.cpp new/src/NS/fix_ns_gmc.cpp --- orig/src/NS/fix_ns_gmc.cpp 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_gmc.cpp 2023-07-20 17:16:05.953402352 -0400 -@@ -0,0 +1,201 @@ ++++ new/src/NS/fix_ns_gmc.cpp 2026-02-26 14:56:30.260636761 -0500 +@@ -0,0 +1,212 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories @@ -584,6 +1611,10 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.cpp new/src/NS/fix_ns_gmc.cpp + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + ++#ifdef DEBUG ++#include ++#endif ++ +#include +#include +#include "fix_ns_gmc.h" @@ -723,11 +1754,18 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.cpp new/src/NS/fix_ns_gmc.cpp + int nlocal = atom->nlocal; + + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "POS ecurrent " << ecurrent << std::endl; ++#endif + + if (ecurrent < Emax) + // continue straight + return; + ++#ifdef DEBUG ++std::cout << "POS " << ecurrent << " >= " << Emax << " reflect" << std::endl; ++#endif ++ + // try to reflect from boundary + double fsum = 0; + double fhatdotdx = 0; @@ -774,8 +1812,8 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.cpp new/src/NS/fix_ns_gmc.cpp + diff -N -u -r orig/src/NS/fix_ns_gmc.h new/src/NS/fix_ns_gmc.h --- orig/src/NS/fix_ns_gmc.h 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_gmc.h 2023-07-20 17:16:05.954402343 -0400 -@@ -0,0 +1,61 @@ ++++ new/src/NS/fix_ns_gmc.h 2026-02-26 14:56:30.261636752 -0500 +@@ -0,0 +1,60 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories @@ -805,11 +1843,11 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.h new/src/NS/fix_ns_gmc.h +class FixNSGMC : public Fix { + public: + FixNSGMC(class LAMMPS *, int, char **); -+ ~FixNSGMC(); -+ int setmask(); -+ virtual void init(); -+ virtual void initial_integrate(int); -+ virtual void final_integrate(); ++ ~FixNSGMC() override; ++ int setmask() override; ++ void init() override; ++ void initial_integrate(int) override; ++ void final_integrate() override; + + protected: + class Compute *pe_compute; @@ -820,7 +1858,116 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.h new/src/NS/fix_ns_gmc.h + double step_size; + int peflag; + char *id_pe; -+ char str[64]; ++}; ++ ++} ++ ++#endif ++#endif ++ ++/* ERROR/WARNING messages: ++ ++E: Illegal ... command ++ ++Self-explanatory. Check the input script syntax and compare to the ++documentation for the command. You can use -echo screen as a ++command-line option when running LAMMPS to see the offending line. ++ ++*/ +diff -N -u -r orig/src/NS/fix_ns.h new/src/NS/fix_ns.h +--- orig/src/NS/fix_ns.h 1969-12-31 19:00:00.000000000 -0500 ++++ new/src/NS/fix_ns.h 2026-02-26 14:56:30.261636752 -0500 +@@ -0,0 +1,106 @@ ++/* -*- c++ -*- ---------------------------------------------------------- ++ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator ++ http://lammps.sandia.gov, Sandia National Laboratories ++ Steve Plimpton, sjplimp@sandia.gov ++ ++ Copyright (2003) Sandia Corporation. Under the terms of Contract ++ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains ++ certain rights in this software. This software is distributed under ++ the GNU General Public License. ++ ++ See the README file in the top-level LAMMPS directory. ++------------------------------------------------------------------------- */ ++ ++#ifdef FIX_CLASS ++ ++FixStyle(ns,FixNS) ++ ++#else ++ ++#ifndef LMP_FIX_NS_H ++#define LMP_FIX_NS_H ++ ++#include "fix.h" ++ ++namespace LAMMPS_NS { ++ ++class FixNS : public Fix { ++ public: ++ FixNS(class LAMMPS *, int, char **); ++ ~FixNS() override; ++ int setmask() override; ++ void init() override; ++ void pre_exchange() override; ++ void initial_integrate(int) override; ++ void final_integrate() override; ++ int pack_forward_comm(int, int *, double *, int, int *) override; ++ void unpack_forward_comm(int, int, double *) override; ++ double compute_vector(int) override; ++ ++ protected: ++ class Compute *pe_compute; ++ double **dx, **prevx; ++ int *prevtype; ++ double prev_boxlo[3], prev_boxhi[3], prev_xy, prev_yz, prev_xz; ++ // global params ++ int seed; ++ double Emax; ++ double pPos, pCell, pType; ++ int state_cur_move; ++ int state_traj_steps_remaining; ++ // pos GMC params ++ int pos_n_steps; ++ double gmc_step_size; ++ // cell MC params ++ int cell_n_steps; ++ double min_aspect_ratio, pressure; ++ int flat_V_prior; ++ double pVol, dVol, pStretch, dStretch, pShear, dShear; ++ // type MC params ++ int type_n_steps; ++ int semi_GC_flag; ++ double *mu; ++ // other data ++ class RanMars *random_g, *random_l; ++ int peflag; ++ char *id_pe; ++ int max_n_steps; ++ ++ // pos GMC internal data ++ int n_attempt_pos, n_success_pos; ++ // cell MC internal data ++ double dPV, cumulative_dPV; ++ bool cell_move_rejected_early; ++ int cell_cur_move; ++ int n_attempt_vol, n_attempt_stretch, n_attempt_shear, n_success_vol, n_success_stretch, n_success_shear; ++ // type MC internal data ++ double dmuN, cumulative_dmuN; ++ bool unequal_cutoffs; ++ int n_attempt_type, n_success_type; ++ ++ private: ++ double cur_aspect_ratio(double cell[3][3]); ++ void pos_gmc_traj_prep(); ++ ++ void pos_gmc_initial_integrate(); ++ void cell_initial_integrate(); ++ void type_initial_integrate(); ++ void pos_gmc_final_integrate(); ++ void cell_final_integrate(); ++ void type_final_integrate(); +}; + +} @@ -839,8 +1986,8 @@ diff -N -u -r orig/src/NS/fix_ns_gmc.h new/src/NS/fix_ns_gmc.h +*/ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp --- orig/src/NS/fix_ns_type.cpp 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_type.cpp 2023-07-20 17:16:05.954402343 -0400 -@@ -0,0 +1,341 @@ ++++ new/src/NS/fix_ns_type.cpp 2026-02-26 14:56:30.262636743 -0500 +@@ -0,0 +1,404 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories @@ -941,6 +2088,13 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + + force_reneighbor = 1; + next_reneighbor = update->ntimestep + 1; ++ // ++ // from MC/fix_atom_swap.cpp ++ // amount of data to be packed for inter-process comm (atom type) ++ if (atom->q_flag) ++ comm_forward = 2; ++ else ++ comm_forward = 1; +} + +/* ---------------------------------------------------------------------- */ @@ -957,6 +2111,53 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + return mask; +} + ++// from MC/fix_atom_swap.cpp ++int FixNSType::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) ++{ ++ int i, j, m; ++ ++ int *type = atom->type; ++ double *q = atom->q; ++ ++ m = 0; ++ ++ if (atom->q_flag) { ++ for (i = 0; i < n; i++) { ++ j = list[i]; ++ buf[m++] = type[j]; ++ buf[m++] = q[j]; ++ } ++ } else { ++ for (i = 0; i < n; i++) { ++ j = list[i]; ++ buf[m++] = type[j]; ++ } ++ } ++ ++ return m; ++} ++ ++// from MC/fix_atom_swap.cpp ++void FixNSType::unpack_forward_comm(int n, int first, double *buf) ++{ ++ int i, m, last; ++ ++ int *type = atom->type; ++ double *q = atom->q; ++ ++ m = 0; ++ last = first + n; ++ ++ if (atom->q_flag) { ++ for (i = first; i < last; i++) { ++ type[i] = static_cast(buf[m++]); ++ q[i] = buf[m++]; ++ } ++ } else { ++ for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); ++ } ++} ++ +/* ---------------------------------------------------------------------- */ + +void FixNSType::init() @@ -1041,8 +2242,8 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + dmuN = mu[new_type-1] - mu[type[atom_0]-1]; + type[atom_0] = new_type; +#ifdef DEBUG -+ std::cout << "SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << -+ type[atom_0] << " mu " << mu[type[atom_0]-1] << " "; ++ std::cout << "TYPE SEMI-GC " << atom_0 << " t " << prevtype[atom_0] << " mu " << mu[prevtype[atom_0]-1] << " -> " << ++ type[atom_0] << " mu " << mu[type[atom_0]-1] << std::endl; +#endif + } else { + bool all_same = true; @@ -1062,7 +2263,7 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + +#ifdef DEBUG +#ifdef SWAP_POS -+ std::cout << "SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << ++ std::cout << "TYPE SWAP i " << atom_0 << " type " << type[atom_0] << " x " << x[atom_0][0] << " " << x[atom_0][1] << " " << x[atom_0][2] << " <-> " << + " i " << atom_1 << " type " << type[atom_1] << " x " << x[atom_1][0] << " " << x[atom_1][1] << " " << x[atom_1][2] << " "; +#else + std::cout << "SWAP i " << atom_0 << " type " << type[atom_0] << " <-> i " << atom_1 << " type " << type[atom_1] << " "; @@ -1089,6 +2290,9 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + dmuN = 0.0; + } + ++ // communicate modified types to ghost atoms ++ comm->forward_comm(this); ++ + // copied from fix_gmc.cpp + // not sure if it's needed on rejections + int id = modify->find_compute("thermo_pe"); @@ -1103,6 +2307,9 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp +void FixNSType::final_integrate() +{ + double ecurrent = modify->compute[modify->find_compute("thermo_pe")]->compute_scalar(); ++#ifdef DEBUG ++std::cout << "TYPE ecurrent " << ecurrent << std::endl; ++#endif + + // if potential energy + d(mu N) is above Emax then reject move + // need to include previous steps' cumulative dmuN contributions, as well as current ones @@ -1132,6 +2339,9 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp + } +#endif + ++ // communicate restored prev types to ghost atoms ++ comm->forward_comm(this); ++ + // NOTE: no idea if this is necessary, or the right place to do this after rejecting a move + // I suspect it's not necessary at all, since next step all these things will be done anyway + if (unequal_cutoffs) { @@ -1184,8 +2394,8 @@ diff -N -u -r orig/src/NS/fix_ns_type.cpp new/src/NS/fix_ns_type.cpp +} diff -N -u -r orig/src/NS/fix_ns_type.h new/src/NS/fix_ns_type.h --- orig/src/NS/fix_ns_type.h 1969-12-31 19:00:00.000000000 -0500 -+++ new/src/NS/fix_ns_type.h 2023-07-20 17:16:05.955402333 -0400 -@@ -0,0 +1,70 @@ ++++ new/src/NS/fix_ns_type.h 2026-02-26 14:56:30.262636743 -0500 +@@ -0,0 +1,71 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories @@ -1215,13 +2425,15 @@ diff -N -u -r orig/src/NS/fix_ns_type.h new/src/NS/fix_ns_type.h +class FixNSType : public Fix { + public: + FixNSType(class LAMMPS *, int, char **); -+ ~FixNSType(); -+ int setmask(); -+ virtual void init(); -+ virtual void pre_exchange(); -+ virtual void initial_integrate(int); -+ virtual void final_integrate(); -+ double compute_vector(int); ++ ~FixNSType() override; ++ int setmask() override; ++ void init() override; ++ void pre_exchange() override; ++ void initial_integrate(int) override; ++ void final_integrate() override; ++ int pack_forward_comm(int, int *, double *, int, int *) override; ++ void unpack_forward_comm(int, int, double *) override; ++ double compute_vector(int) override; + + protected: + class Compute *pe_compute; @@ -1234,7 +2446,6 @@ diff -N -u -r orig/src/NS/fix_ns_type.h new/src/NS/fix_ns_type.h + double *mu; + int peflag; + char *id_pe; -+ char str[64]; + double dmuN, cumulative_dmuN; + + bool unequal_cutoffs; diff --git a/pymatnext/ns_configs/ase_atoms/__init__.py b/pymatnext/ns_configs/ase_atoms/__init__.py index 7944209..bcc6f54 100644 --- a/pymatnext/ns_configs/ase_atoms/__init__.py +++ b/pymatnext/ns_configs/ase_atoms/__init__.py @@ -455,16 +455,23 @@ def _prep_walk(self, params, vol_per_atom=None): self.walk_func = {} if self.calc_type == "ASE": from .walks_ase_calc import walk_pos_gmc, walk_cell, walk_type + walk_combined = None elif self.calc_type == "LAMMPS": - from .walks_lammps import walk_pos_gmc, walk_cell, walk_type + from .walks_lammps import walk_pos_gmc, walk_cell, walk_type, walk_combined else: raise NotImplementedError(f"Unknown calculator type {self.calc_type}") - self.walk_func["gmc"] = walk_pos_gmc - self.walk_func["cell"] = walk_cell - self.walk_func["type"] = walk_type + if params["combined"]: + if not walk_combined: + raise RuntimeError(f"Got walk_combined but no walk_combined function for calc_type {self.calc_type}") + self.walk_func["combined"] = walk_combined + else: + self.walk_func["gmc"] = walk_pos_gmc + self.walk_func["cell"] = walk_cell + self.walk_func["type"] = walk_type - assert set(self.walk_func.keys()) == set(NSConfig_ASE_Atoms._walk_moves) + # make sure all moves are valid + assert set(self.walk_func.keys()) == set(NSConfig_ASE_Atoms._walk_moves) if self.walk_prob[NSConfig_ASE_Atoms._walk_moves.index("gmc")] > 0.0: self.atoms.new_array("NS_velocities", np.zeros(self.atoms.positions.shape)) @@ -713,19 +720,27 @@ def walk(self, Emax, walk_len, rng): dict(str param_name: ndarray([int n_attempts, int n_success])) dict containing number of attempted moves and number of successful moves for each move type """ - # if we fixed the number of steps for every move type, and only varied proportions, - # we could do all the rng in a single call self.atoms.calc = self.calc - walk_len_so_far = 0 - while walk_len_so_far < walk_len: - move = rng.choice(NSConfig_ASE_Atoms._walk_moves, p=self.walk_prob) + if self.walk_func.get("combined"): # returns list of tuples with move param attempt/success statistics - n_att_acc_walk = self.walk_func[move](self, Emax, rng) + n_att_acc_walk = self.walk_func["combined"](self, Emax, rng, walk_len) for param, n_att, n_acc in n_att_acc_walk: self.n_att_acc[param] += (n_att, n_acc) - - walk_len_so_far += self.walk_traj_len[move] + else: + walk_len_so_far = 0 + while walk_len_so_far < walk_len: + # if we fixed the number of steps for every move type, and only varied proportions, + # we could do all the rng move selection in a single call, because the number of moves + # we need would be known ahead of time + move = rng.choice(NSConfig_ASE_Atoms._walk_moves, p=self.walk_prob) + + # returns list of tuples with move param attempt/success statistics + n_att_acc_walk = self.walk_func[move](self, Emax, rng) + for param, n_att, n_acc in n_att_acc_walk: + self.n_att_acc[param] += (n_att, n_acc) + + walk_len_so_far += self.walk_traj_len[move] # move walks keep track of NS_energy_shift by accumulating changes. It would definitely # be more stable to recalculate shift from scratch, although there's a chance it might lead to diff --git a/pymatnext/ns_configs/ase_atoms/ase_atoms_params.py b/pymatnext/ns_configs/ase_atoms/ase_atoms_params.py index 8dbc7fd..b5f1648 100644 --- a/pymatnext/ns_configs/ase_atoms/ase_atoms_params.py +++ b/pymatnext/ns_configs/ase_atoms/ase_atoms_params.py @@ -55,5 +55,7 @@ "type": { "sGC": False, "mu": { "_IGNORE_": True } - } + }, + + "combined": False } diff --git a/pymatnext/ns_configs/ase_atoms/walks_lammps.py b/pymatnext/ns_configs/ase_atoms/walks_lammps.py index bdde260..03b7a4b 100644 --- a/pymatnext/ns_configs/ase_atoms/walks_lammps.py +++ b/pymatnext/ns_configs/ase_atoms/walks_lammps.py @@ -135,6 +135,88 @@ def extract_E_F(lmp, recalc): return (lmp.extract_compute("pe", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_SCALAR), lmp.numpy.extract_atom("f")[:nlocal]) +def walk_combined(ns_atoms, Emax, rng, walk_len): + ns_atoms.calc.command("unfix NS") + atoms = ns_atoms.atoms + # for LAMMPS RanMars RNG + lammps_seed = rng.integers(1, 900000000) + + set_lammps_from_atoms(ns_atoms) + + Emax -= atoms.info["NS_energy_shift"] + pPos = ns_atoms.walk_prob[ns_atoms._walk_moves.index("gmc")] + pCell = ns_atoms.walk_prob[ns_atoms._walk_moves.index("cell")] + pType = ns_atoms.walk_prob[ns_atoms._walk_moves.index("type")] + fix_cmd = f"fix NS all ns {lammps_seed} {Emax} {pPos} {pCell} {pType}" + + if pPos > 0.0: + fix_cmd += f" {ns_atoms.walk_traj_len['gmc']} {ns_atoms.step_size['pos_gmc_each_atom']}" + if pCell > 0.0: + move_params_cell = ns_atoms.move_params["cell"] + submove_probs = move_params_cell["submove_probabilities"] + N_atoms = len(atoms) + step_size_volume = ns_atoms.step_size["cell_volume_per_atom"] * N_atoms + step_size_shear = ns_atoms.step_size["cell_shear_per_rt3_atom"] * (N_atoms ** (1.0 / 3.0)) + step_size_stretch = ns_atoms.step_size["cell_stretch"] + fix_cmd += (f" {ns_atoms.walk_traj_len['cell']} " + f"{move_params_cell['min_aspect_ratio']} {ns_atoms.pressure} " + f"{submove_probs['volume']} {step_size_volume} " + f"{submove_probs['stretch']} {step_size_stretch} " + f"{submove_probs['shear']} {step_size_shear} " + f"{'yes' if ns_atoms.move_params['cell']['flat_V_prior'] else 'no'}") + if pType > 0.0: + move_params_type = ns_atoms.move_params["type"] + fix_cmd += f" {ns_atoms.walk_traj_len['type']}" + if move_params_type["sGC"]: + # NOTE: if we don't mind being unable to change mu, we can construct this + # string once when we set up the calculator + fix_cmd += " yes " + " ".join([f"{mu:.10f}" for mu in ns_atoms.mu[ns_atoms.Z_of_type[1:]]]) + else: + fix_cmd += " no" + + ns_atoms.calc.command(fix_cmd) + try: + ns_atoms.calc.command(f"run {walk_len} post no") + failed = False + except Exception as exc: + exc_str = str(exc) + warnings.warn(f"LAMMPS ns/gmc run raised exception {exc_str}") + if "Lost atoms" in exc_str: + # arrays will now be wrong size, and will fail in next call to set_lammps_from_atoms + ns_atoms.end_calculator() + # don't store new calculator results, since they may be from wrong positions + ns_atoms.init_calculator(skip_initial_store=True) + failed = True + + if not failed: + # recompute, in case last step inside LAMMPS was a rejection + E, F = extract_E_F(ns_atoms.calc, True) + # set atoms from current lammps internal state + types, pos, vel = get_pointers_from_lammps(ns_atoms) + set_atoms_from_lammps(ns_atoms, types, pos, vel, E, F, True) + # wrap to avoid atoms moving far enough in periodic images for lammps to lose them + atoms.wrap() + + # gather number of attempted and accepted moves from fix + # positions 0-1 + pos_n_att = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 0, 0)) + pos_n_acc = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 1, 0)) + # cell 2-7 + cell_n_att = {} + cell_n_acc = {} + for submove_i, submove_type in enumerate(["volume", "stretch", "shear"]): + cell_n_att[submove_type] = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 2 + 2 * submove_i + 0, 0)) + cell_n_acc[submove_type] = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 2 + 2 * submove_i + 1, 0)) + ## no need to extract for type + ## type 8-9 + ## type_n_att = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 2 + 2 * 3 + 0, 0)) + ## type_n_acc = int(ns_atoms.calc.extract_fix("NS", lammps.LMP_STYLE_GLOBAL, lammps.LMP_TYPE_VECTOR, 2 + 2 * 3 + 1, 0)) + + return [("pos_gmc_each_atom", pos_n_att, pos_n_acc), + ("cell_volume_per_atom", cell_n_att["volume"], cell_n_acc["volume"]), + ("cell_shear_per_rt3_atom", cell_n_att["shear"], cell_n_acc["shear"]), + ("cell_stretch", cell_n_att["stretch"], cell_n_acc["stretch"])] + def walk_pos_gmc(ns_atoms, Emax, rng): """walk atomic positions with GMC @@ -186,7 +268,6 @@ def walk_pos_gmc(ns_atoms, Emax, rng): return [("pos_gmc_each_atom", 1, 0 if reject else 1)] - def walk_cell(ns_atoms, Emax, rng): """walk atomic cell with MC diff --git a/tests/test_sample.py b/tests/test_sample.py index cc40839..063dad4 100644 --- a/tests/test_sample.py +++ b/tests/test_sample.py @@ -14,8 +14,10 @@ try: import lammps + lammps_toml_files = ["LAMMPS"] except: lammps = None + lammps_toml_files = [] # tests without MPI @@ -32,12 +34,19 @@ def test_EAM_LAMMPS_no_mpi(tmp_path, monkeypatch): do_EAM_LAMMPS(tmp_path, monkeypatch, using_mpi=False, max_iter=300) @pytest.mark.mpi_skip -def test_pressure_no_mpi(tmp_path, monkeypatch): - do_pressure(tmp_path, monkeypatch, using_mpi=False) +@pytest.mark.parametrize("toml_file", ["ASE"] + lammps_toml_files) +def test_pressure_no_mpi(tmp_path, monkeypatch, toml_file): + do_pressure(tmp_path, monkeypatch, toml_file, using_mpi=False) @pytest.mark.mpi_skip -def test_sGC_no_mpi(tmp_path, monkeypatch): - do_sGC(tmp_path, monkeypatch, using_mpi=False) +@pytest.mark.parametrize("toml_file", ["matscipy"] + lammps_toml_files) +def test_sGC_no_mpi(tmp_path, monkeypatch, toml_file): + do_sGC(tmp_path, monkeypatch, toml_file, using_mpi=False) + +@pytest.mark.mpi_skip +@pytest.mark.parametrize("toml_file", lammps_toml_files) +def test_sGC_combined_no_mpi(tmp_path, monkeypatch, toml_file): + do_sGC(tmp_path, monkeypatch, toml_file, using_mpi=False, combined=True) # tests with MPI @@ -59,12 +68,14 @@ def test_EAM_LAMMPS_mpi(mpi_tmp_path, monkeypatch): do_EAM_LAMMPS(mpi_tmp_path, monkeypatch, using_mpi=True, max_iter=300) @pytest.mark.mpi -def test_pressure_mpi(mpi_tmp_path, monkeypatch): - do_pressure(mpi_tmp_path, monkeypatch, using_mpi=True) +@pytest.mark.parametrize("toml_file", ["ASE"] + lammps_toml_files) +def test_pressure_mpi(mpi_tmp_path, monkeypatch, toml_file): + do_pressure(mpi_tmp_path, monkeypatch, toml_file, using_mpi=True) @pytest.mark.mpi -def test_sGC_mpi(mpi_tmp_path, monkeypatch): - do_sGC(mpi_tmp_path, monkeypatch, using_mpi=True) +@pytest.mark.parametrize("toml_file", ["matscipy"] + lammps_toml_files) +def test_sGC_mpi(mpi_tmp_path, monkeypatch, toml_file): + do_sGC(mpi_tmp_path, monkeypatch, toml_file, using_mpi=True) def do_Morse_ASE_restart(tmp_path, monkeypatch, using_mpi): @@ -264,7 +275,7 @@ def do_EAM_LAMMPS(tmp_path, monkeypatch, using_mpi, max_iter=None): assert False, f"test {fields} ref {fields_ref}" -def do_pressure(tmp_path, monkeypatch, using_mpi): +def do_pressure(tmp_path, monkeypatch, toml_file, using_mpi): if using_mpi: from mpi4py import MPI if MPI.COMM_WORLD.size != 2: @@ -272,50 +283,47 @@ def do_pressure(tmp_path, monkeypatch, using_mpi): assets_dir = Path(__file__).parent / 'assets' / 'do_pressure' - toml_files = ['params_ASE.toml'] - if lammps is not None: - toml_files += ['params_LAMMPS.toml'] + toml_file = f'params_{toml_file}.toml' - for toml_file in toml_files: - with open(assets_dir / toml_file) as fin, open(tmp_path / toml_file, 'w') as fout: - for l in fin: - # fix output_filename_prefix so everything is written to tmp_path - if "output_filename_prefix" in l: - fout.write(f'output_filename_prefix = "{tmp_path}/{toml_file}"\n') - continue + with open(assets_dir / toml_file) as fin, open(tmp_path / toml_file, 'w') as fout: + for l in fin: + # fix output_filename_prefix so everything is written to tmp_path + if "output_filename_prefix" in l: + fout.write(f'output_filename_prefix = "{tmp_path}/{toml_file}"\n') + continue - fout.write(l) + fout.write(l) - # add assets dir for Morse.py module - sys.path.insert(0, str(assets_dir)) - # add assets dir for Morse.py module - if not using_mpi: - # need PYMATNEXT_NO_MPI since otherwise rngs will vary depending on number of MPI processes, and - # numbers in output will vary - monkeypatch.setenv("PYMATNEXT_NO_MPI", "1") + # add assets dir for Morse.py module + sys.path.insert(0, str(assets_dir)) + # add assets dir for Morse.py module + if not using_mpi: + # need PYMATNEXT_NO_MPI since otherwise rngs will vary depending on number of MPI processes, and + # numbers in output will vary + monkeypatch.setenv("PYMATNEXT_NO_MPI", "1") - main_args = ['--override_param', '/global/random_seed', '5', '--override_param', '/global/output_filename_prefix_extra', '.test', - str(tmp_path / toml_file)] + main_args = ['--override_param', '/global/random_seed', '5', '--override_param', '/global/output_filename_prefix_extra', '.test', + str(tmp_path / toml_file)] - sample.main(main_args, mpi_finalize=False) - del sys.path[0] + sample.main(main_args, mpi_finalize=False) + del sys.path[0] - with open(tmp_path / (toml_file + ".test.NS_samples")) as fin: - lfirst = None - for l in fin: - if l.startswith("#"): - continue + with open(tmp_path / (toml_file + ".test.NS_samples")) as fin: + lfirst = None + for l in fin: + if l.startswith("#"): + continue - if lfirst is None: - lfirst = l - llast = l + if lfirst is None: + lfirst = l + llast = l - Vfirst = float(lfirst.strip().split()[2]) - Vlast = float(llast.strip().split()[2]) - assert Vfirst / Vlast > 10 + Vfirst = float(lfirst.strip().split()[2]) + Vlast = float(llast.strip().split()[2]) + assert Vfirst / Vlast > 10 -def do_sGC(tmp_path, monkeypatch, using_mpi): +def do_sGC(tmp_path, monkeypatch, toml_file, using_mpi, combined=False): if using_mpi: from mpi4py import MPI if MPI.COMM_WORLD.size != 2: @@ -323,51 +331,51 @@ def do_sGC(tmp_path, monkeypatch, using_mpi): assets_dir = Path(__file__).parent / 'assets' / 'do_sGC' - toml_files = ['params_matscipy.toml'] - if lammps is not None: - toml_files += ['params_LAMMPS.toml'] + toml_file = f'params_{toml_file}.toml' - for toml_file in toml_files: - with open(assets_dir / toml_file) as fin, open(tmp_path / toml_file, 'w') as fout: - for l in fin: - # rewrite to find potentials in right place - if "_POTENTIAL_DIR_" in l: - fout.write(l.replace("_POTENTIAL_DIR_", str(assets_dir.resolve()))) - continue + with open(assets_dir / toml_file) as fin, open(tmp_path / toml_file, 'w') as fout: + for l in fin: + # rewrite to find potentials in right place + if "_POTENTIAL_DIR_" in l: + fout.write(l.replace("_POTENTIAL_DIR_", str(assets_dir.resolve()))) + continue - # fix output_filename_prefix so everything is written to tmp_path - if "output_filename_prefix" in l: - fout.write(f'output_filename_prefix = "{tmp_path}/{toml_file}"\n') - continue + # fix output_filename_prefix so everything is written to tmp_path + if "output_filename_prefix" in l: + fout.write(f'output_filename_prefix = "{tmp_path}/{toml_file}"\n') + continue - fout.write(l) + fout.write(l) - # add assets dir for Morse.py module - sys.path.insert(0, str(assets_dir)) - # add assets dir for Morse.py module - if not using_mpi: - # need PYMATNEXT_NO_MPI since otherwise rngs will vary depending on number of MPI processes, and - # numbers in output will vary - monkeypatch.setenv("PYMATNEXT_NO_MPI", "1") + if combined and "[configs.walk]" in l: + fout.write("combined = true") - main_args = ['--override_param', '/global/random_seed', '5', '--override_param', '/global/output_filename_prefix_extra', '.test', - str(tmp_path / toml_file)] + # add assets dir for Morse.py module + sys.path.insert(0, str(assets_dir)) + # add assets dir for Morse.py module + if not using_mpi: + # need PYMATNEXT_NO_MPI since otherwise rngs will vary depending on number of MPI processes, and + # numbers in output will vary + monkeypatch.setenv("PYMATNEXT_NO_MPI", "1") - sample.main(main_args, mpi_finalize=False) - del sys.path[0] + main_args = ['--override_param', '/global/random_seed', '5', '--override_param', '/global/output_filename_prefix_extra', '.test', + str(tmp_path / toml_file)] - with open(tmp_path / (toml_file + ".test.NS_samples")) as fin: - lfirst = None - for l in fin: - if l.startswith("#"): - continue + sample.main(main_args, mpi_finalize=False) + del sys.path[0] + + with open(tmp_path / (toml_file + ".test.NS_samples")) as fin: + lfirst = None + for l in fin: + if l.startswith("#"): + continue - if lfirst is None: - lfirst = l - llast = l + if lfirst is None: + lfirst = l + llast = l - f_13_first, f_29_first = [float(f) for f in lfirst.strip().split()[5:7]] - f_13_last, f_29_last = [float(f) for f in llast.strip().split()[5:7]] + f_13_first, f_29_first = [float(f) for f in lfirst.strip().split()[5:7]] + f_13_last, f_29_last = [float(f) for f in llast.strip().split()[5:7]] - assert f_29_first / f_13_first == 1.0 - assert f_29_last / f_13_last > 4 + assert f_29_first / f_13_first == 1.0 + assert f_29_last / f_13_last > 4