From 5b6e8884b498988694d72447452063026f07114b Mon Sep 17 00:00:00 2001 From: Taylor Barnes Date: Thu, 6 Apr 2023 18:11:38 +0000 Subject: [PATCH 1/8] Add support for non-periodic systems in the MDI fix --- src/MDI/fix_mdi_qm.cpp | 74 ++++++++++++++++++++++++++++++++++-------- src/MDI/fix_mdi_qm.h | 2 +- 2 files changed, 61 insertions(+), 15 deletions(-) diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index ee98429ccb..e9f608c53e 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/src/MDI/fix_mdi_qm.cpp @@ -325,6 +325,14 @@ void FixMDIQM::init() ierr = MDI_Check_command_exists("@DEFAULT", "all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: CELL_DISPL", mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); - ierr = MDI_Send(qm_cell_displ, 3, MDI_DOUBLE, mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data"); - } + // Only send the cell dimensions if this is a periodic simulation + if ( domain->xperiodic == 1 && + domain->yperiodic == 1 && + domain->zperiodic == 1 ) { - ierr = MDI_Send_command(">CELL", mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL command"); - ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL data"); + if (celldispl_exists) { + ierr = MDI_Send_command(">CELL_DISPL", mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); + ierr = MDI_Send(qm_cell_displ, 3, MDI_DOUBLE, mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data"); + } + + ierr = MDI_Send_command(">CELL", mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL command"); + ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL data"); + + } + else { + + if ( domain->xperiodic == 1 || + domain->yperiodic == 1 || + domain->zperiodic == 1 ) { + + error->all(FLERR, "MDI: The QM driver does not work for simulations with both periodic and non-periodic dimensions."); + + } + + } } /* ---------------------------------------------------------------------- diff --git a/src/MDI/fix_mdi_qm.h b/src/MDI/fix_mdi_qm.h index b86fe90d93..fd31012cac 100644 --- a/src/MDI/fix_mdi_qm.h +++ b/src/MDI/fix_mdi_qm.h @@ -62,7 +62,7 @@ class FixMDIQM : public Fix { MDI_Comm mdicomm; int natoms_exists,celldispl_exists,elements_exists,types_exists; - int stress_exists; + int stress_exists, pe_exists, keelec_exists; int nmax; From 276e1dd12b2e1724d1daa614b292d9a2cc580b2b Mon Sep 17 00:00:00 2001 From: Taylor Barnes Date: Fri, 7 Apr 2023 14:24:23 +0000 Subject: [PATCH 2/8] Add support for the FORCES"); @@ -639,6 +644,7 @@ void MDIEngine::mdi_commands() MDI_Register_command("@ENDSTEP", "all(FLERR, "MDI: all(FLERR, "MDI: Date: Fri, 7 Apr 2023 18:27:08 +0000 Subject: [PATCH 3/8] Add support for non-periodic calculations to MDI QMMM --- src/MDI/fix_mdi_qmmm.cpp | 96 +++++++++++++++++++++++++++++++--------- src/MDI/fix_mdi_qmmm.h | 4 +- 2 files changed, 78 insertions(+), 22 deletions(-) diff --git a/src/MDI/fix_mdi_qmmm.cpp b/src/MDI/fix_mdi_qmmm.cpp index 196dbde655..a7c6defdc1 100644 --- a/src/MDI/fix_mdi_qmmm.cpp +++ b/src/MDI/fix_mdi_qmmm.cpp @@ -360,6 +360,14 @@ void FixMDIQMMM::init() if (ierr) error->all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: CELL_DISPL", mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); - ierr = MDI_Send(qm_cell_displ, 3, MDI_DOUBLE, mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data"); + // Only send the cell dimensions if this is a periodic simulation + if ( domain->xperiodic == 1 && + domain->yperiodic == 1 && + domain->zperiodic == 1 ) { + + if (celldispl_exists) { + ierr = MDI_Send_command(">CELL_DISPL", mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); + ierr = MDI_Send(qm_cell_displ, 3, MDI_DOUBLE, mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data"); + } + + ierr = MDI_Send_command(">CELL", mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL command"); + ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL data"); + + } + else { + + if ( domain->xperiodic == 1 || + domain->yperiodic == 1 || + domain->zperiodic == 1 ) { + + error->all(FLERR, "MDI: The QM driver does not work for simulations with both periodic and non-periodic dimensions."); + + } + } - ierr = MDI_Send_command(">CELL", mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL command"); - ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); - if (ierr) error->all(FLERR, "MDI: >CELL data"); +} + +/* ---------------------------------------------------------------------- + get the energy of the QM system +------------------------------------------------------------------------- */ + +void FixMDIQMMM::get_qm_energy() +{ + int ierr; + + if (pe_exists && keelec_exists) { + int pe_energy, keelec_energy; + + // get the total potential energy + ierr = MDI_Send_command("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: Date: Fri, 7 Apr 2023 18:41:12 +0000 Subject: [PATCH 4/8] Add periodicity warning to MDI documentation --- doc/src/fix_mdi_qm.rst | 2 ++ doc/src/fix_mdi_qmmm.rst | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/src/fix_mdi_qm.rst b/doc/src/fix_mdi_qm.rst index 46b62fe8db..007c2c275a 100644 --- a/doc/src/fix_mdi_qm.rst +++ b/doc/src/fix_mdi_qm.rst @@ -294,6 +294,8 @@ LAMMPS can also be used as an MDI driver in other unit choices it supports, e.g. *lj*, but then no unit conversion to MDI units is performed. +If this fix is used in conjuction with a QM code that does not support periodic boundary conditions (more specifically, a QM code that does not support the ``>CELL`` MDI command), the LAMMPS input file must specify non-periodic boundary conditions. + Related commands """""""""""""""" diff --git a/doc/src/fix_mdi_qmmm.rst b/doc/src/fix_mdi_qmmm.rst index 3485e5152b..4f81fbb999 100644 --- a/doc/src/fix_mdi_qmmm.rst +++ b/doc/src/fix_mdi_qmmm.rst @@ -261,6 +261,8 @@ used to specify *real* or *metal* units. This will ensure the correct unit conversions between LAMMPS and MDI units. The other code will also perform similar unit conversions into its preferred units. +If this fix is used in conjuction with a QM code that does not support periodic boundary conditions (more specifically, a QM code that does not support the ``>CELL`` MDI command), the LAMMPS input file must specify non-periodic boundary conditions. + Related commands """""""""""""""" From 0f07c5e809745b8d2084dbf19cd318b3c11ecb04 Mon Sep 17 00:00:00 2001 From: Taylor Barnes Date: Fri, 7 Apr 2023 18:56:40 +0000 Subject: [PATCH 5/8] Run clang on MDI fixes --- src/MDI/fix_mdi_qm.cpp | 183 +++++++++++------------ src/MDI/fix_mdi_qm.h | 16 +- src/MDI/fix_mdi_qmmm.cpp | 311 +++++++++++++++++++-------------------- src/MDI/fix_mdi_qmmm.h | 51 ++++--- 4 files changed, 279 insertions(+), 282 deletions(-) diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index e9f608c53e..79d85ebced 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/src/MDI/fix_mdi_qm.cpp @@ -43,7 +43,7 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (atom->tag_enable == 0) error->all(FLERR, "Cannot use MDI engine without atom IDs"); if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Fix mdi/qm requires an atom map be defined"); + error->all(FLERR, "Fix mdi/qm requires an atom map be defined"); // confirm LAMMPS is being run as a driver @@ -99,18 +99,14 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg], "elements") == 0) { const char *symbols[] = { - "H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne", - "Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca", - "Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", - "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y" , "Zr", - "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", - "Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", - "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", - "Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg", - "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", - "Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", - "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", - "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og", + "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", + "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", + "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", + "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", + "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", + "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", + "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", + "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og", }; int ntypes = atom->ntypes; @@ -120,18 +116,17 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) for (int i = 1; i <= ntypes; i++) { int anum; for (anum = 0; anum < MAXELEMENT; anum++) - if (strcmp(arg[iarg + i],symbols[anum]) == 0) break; - if (anum == MAXELEMENT) - error->all(FLERR,"Invalid chemical element in fix mdi/qm command"); + if (strcmp(arg[iarg + i], symbols[anum]) == 0) break; + if (anum == MAXELEMENT) error->all(FLERR, "Invalid chemical element in fix mdi/qm command"); elements[i] = anum + 1; } iarg += ntypes + 1; - } else if (strcmp(arg[iarg],"mc") == 0) { - if (iarg+2 > narg) error->all(FLERR, "Illegal fix mdi/qm command"); + } else if (strcmp(arg[iarg], "mc") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix mdi/qm command"); mcflag = 1; delete[] id_mcfix; - id_mcfix = utils::strdup(arg[iarg+1]); + id_mcfix = utils::strdup(arg[iarg + 1]); iarg += 2; } else @@ -198,7 +193,7 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // per-atom data nmax = atom->nmax; - memory->create(array_atom,nmax,3,"mdi/qm:array_atom"); + memory->create(array_atom, nmax, 3, "mdi/qm:array_atom"); // initialize outputs @@ -210,8 +205,7 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) sumflag = 0; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - array_atom[i][0] = array_atom[i][1] = array_atom[i][2] = 0.0; + for (int i = 0; i < nlocal; i++) array_atom[i][0] = array_atom[i][1] = array_atom[i][2] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -341,11 +335,12 @@ void FixMDIQM::init() if (mcflag) { Fix *f_mc = modify->get_fix_by_id(id_mcfix); - if (!f_mc) error->all(FLERR,"Fix mdi/qm could not find Monte Carlo fix ID {}", id_mcfix); + if (!f_mc) error->all(FLERR, "Fix mdi/qm could not find Monte Carlo fix ID {}", id_mcfix); int dim; mc_active_ptr = (int *) f_mc->extract("mc_active", dim); if (!mc_active_ptr || dim != 0) - error->all(FLERR,"Fix mdi/qm could not query mc_active from Monte Carlo fix ID {}", id_mcfix); + error->all(FLERR, "Fix mdi/qm could not query mc_active from Monte Carlo fix ID {}", + id_mcfix); exclusion_group_ptr = (int *) f_mc->extract("exclusion_group", dim); } @@ -371,11 +366,12 @@ void FixMDIQM::init() // check if box has changed - if (new_system) set_box(); + if (new_system) + set_box(); else { - double old_cell[9],old_cell_displ[3]; - memcpy(old_cell,qm_cell,9*sizeof(double)); - memcpy(old_cell_displ,qm_cell_displ,3*sizeof(double)); + double old_cell[9], old_cell_displ[3]; + memcpy(old_cell, qm_cell, 9 * sizeof(double)); + memcpy(old_cell_displ, qm_cell_displ, 3 * sizeof(double)); set_box(); for (int i = 0; i < 9; i++) if (qm_cell[i] != old_cell[i]) new_system = 1; @@ -386,11 +382,12 @@ void FixMDIQM::init() // check if atom elements or types have changed if (elements && elements_exists) { - if (new_system) set_eqm(); + if (new_system) + set_eqm(); else { int *eqm_old; - memory->create(eqm_old,nqm,"mdi/qm:eqm_old"); - memcpy(eqm_old,eqm,nqm*sizeof(int)); + memory->create(eqm_old, nqm, "mdi/qm:eqm_old"); + memcpy(eqm_old, eqm, nqm * sizeof(int)); set_eqm(); for (int i = 0; i < nqm; i++) if (eqm[i] != eqm_old[i]) new_system = 1; @@ -398,11 +395,12 @@ void FixMDIQM::init() } } else if (types_exists) { - if (new_system) set_tqm(); + if (new_system) + set_tqm(); else { int *tqm_old; - memory->create(tqm_old,nqm,"mdi/qm:tqm_old"); - memcpy(tqm_old,tqm,nqm*sizeof(int)); + memory->create(tqm_old, nqm, "mdi/qm:tqm_old"); + memcpy(tqm_old, tqm, nqm * sizeof(int)); set_tqm(); for (int i = 0; i < nqm; i++) if (tqm[i] != tqm_old[i]) new_system = 1; @@ -417,8 +415,10 @@ void FixMDIQM::init() if (new_system) { send_natoms(); send_box(); - if (elements && elements_exists) send_elements(); - else if (types_exists) send_types(); + if (elements && elements_exists) + send_elements(); + else if (types_exists) + send_types(); nqm_last = nqm; } } @@ -459,7 +459,7 @@ void FixMDIQM::post_force(int vflag) if (atom->nmax > nmax) { nmax = atom->nmax; memory->destroy(array_atom); - memory->create(array_atom,nmax,3,"mdi/qm:array_atom"); + memory->create(array_atom, nmax, 3, "mdi/qm:array_atom"); } // determine whether a new vs incremental QM calc is needed @@ -468,8 +468,10 @@ void FixMDIQM::post_force(int vflag) // incremental when a system is slowly evolving (AIMD) int new_system = 0; - if (nqm != nqm_last) new_system = 1; - else if (mcflag && *mc_active_ptr) new_system = 1; + if (nqm != nqm_last) + new_system = 1; + else if (mcflag && *mc_active_ptr) + new_system = 1; // send new system info to MDI engine: atom count and elements/types // reset QM data structures if atom count has changed @@ -493,8 +495,8 @@ void FixMDIQM::post_force(int vflag) nqm_last = nqm; - // incremental system - // if simulation box dynamically changes, send current box to MDI engine + // incremental system + // if simulation box dynamically changes, send current box to MDI engine } else if (domain->box_change_size || domain->box_change_shape) { set_box(); @@ -530,8 +532,7 @@ void FixMDIQM::post_force(int vflag) if (ierr) error->all(FLERR, "MDI: all(FLERR, "MDI: dimension == 2) volume = domain->xprd * domain->yprd; else if (domain->dimension == 3) - volume = domain->xprd * domain->yprd * domain->zprd; + volume = domain->xprd * domain->yprd * domain->zprd; for (int i = 0; i < 6; i++) virial[i] = qm_virial_symmetric[i] * volume / nprocs; } } @@ -670,10 +671,10 @@ double FixMDIQM::compute_vector(int n) double FixMDIQM::memory_usage() { double bytes = 0.0; - bytes += nqm * sizeof(tagint); // qmIDs - bytes += nqm * sizeof(int); // qm2owned - bytes += 4 * nqm * sizeof(int); // int QM arrays/vecs - bytes += 3*3 * nqm * sizeof(double); // double QM arrays/vecs + bytes += nqm * sizeof(tagint); // qmIDs + bytes += nqm * sizeof(int); // qm2owned + bytes += 4 * nqm * sizeof(int); // int QM arrays/vecs + bytes += 3 * 3 * nqm * sizeof(double); // double QM arrays/vecs return bytes; } @@ -703,17 +704,17 @@ void FixMDIQM::reallocate() memory->destroy(tqm_mine); memory->destroy(xqm_mine); - memory->create(qmIDs,max_nqm,"mdi/qm:qmIDs"); - memory->create(qm2owned,max_nqm,"mdi/qm:qm2owned"); + memory->create(qmIDs, max_nqm, "mdi/qm:qmIDs"); + memory->create(qm2owned, max_nqm, "mdi/qm:qm2owned"); - memory->create(eqm,max_nqm,"mdi/qm:eqm"); - memory->create(tqm,max_nqm,"mdi/qm:tqm"); - memory->create(xqm,max_nqm,3,"mdi/qm:xqm"); - memory->create(fqm,max_nqm,3,"mdi/qm:fqm"); + memory->create(eqm, max_nqm, "mdi/qm:eqm"); + memory->create(tqm, max_nqm, "mdi/qm:tqm"); + memory->create(xqm, max_nqm, 3, "mdi/qm:xqm"); + memory->create(fqm, max_nqm, 3, "mdi/qm:fqm"); - memory->create(eqm_mine,max_nqm,"mdi/qm:eqm_mine"); - memory->create(tqm_mine,max_nqm,"mdi/qm:tqm_mine"); - memory->create(xqm_mine,max_nqm,3,"mdi/qm:xqm_mine"); + memory->create(eqm_mine, max_nqm, "mdi/qm:eqm_mine"); + memory->create(tqm_mine, max_nqm, "mdi/qm:tqm_mine"); + memory->create(xqm_mine, max_nqm, 3, "mdi/qm:xqm_mine"); } /* ---------------------------------------------------------------------- @@ -726,8 +727,7 @@ int FixMDIQM::set_nqm() { // require 3*nqm be a small INT, so can MPI_Allreduce xqm - if (3*atom->natoms > MAXSMALLINT) - error->all(FLERR,"Fix mdi/qm has too many atoms"); + if (3 * atom->natoms > MAXSMALLINT) error->all(FLERR, "Fix mdi/qm has too many atoms"); int ncount = atom->natoms; nexclude = 0; @@ -746,7 +746,7 @@ int FixMDIQM::set_nqm() if (mask[i] & excludebit) nexclude_mine = 1; } - MPI_Allreduce(&nexclude_mine,&nexclude,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&nexclude_mine, &nexclude, 1, MPI_INT, MPI_SUM, world); } ncount -= nexclude; @@ -776,31 +776,35 @@ void FixMDIQM::create_qm_list() int nqm_mine = 0; for (int i = 0; i < nlocal; i++) { - if (!nexclude) nqm_mine++; - else if (!(mask[i] & excludebit)) nqm_mine++; + if (!nexclude) + nqm_mine++; + else if (!(mask[i] & excludebit)) + nqm_mine++; } tagint *qmIDs_mine; - memory->create(qmIDs_mine,nqm_mine,"mdi/qm:qmIDs_mine"); + memory->create(qmIDs_mine, nqm_mine, "mdi/qm:qmIDs_mine"); nqm_mine = 0; for (int i = 0; i < nlocal; i++) { - if (!nexclude) qmIDs_mine[nqm_mine++] = tag[i]; - else if (!(mask[i] & excludebit)) qmIDs_mine[nqm_mine++] = tag[i]; + if (!nexclude) + qmIDs_mine[nqm_mine++] = tag[i]; + else if (!(mask[i] & excludebit)) + qmIDs_mine[nqm_mine++] = tag[i]; } int *recvcounts, *displs; - memory->create(recvcounts,nprocs,"mdi/qm:recvcounts"); - memory->create(displs,nprocs,"mdi/qm:displs"); + memory->create(recvcounts, nprocs, "mdi/qm:recvcounts"); + memory->create(displs, nprocs, "mdi/qm:displs"); - MPI_Allgather(&nqm_mine,1,MPI_INT,recvcounts,1,MPI_INT,world); + MPI_Allgather(&nqm_mine, 1, MPI_INT, recvcounts, 1, MPI_INT, world); displs[0] = 0; for (int iproc = 1; iproc < nprocs; iproc++) - displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; + displs[iproc] = displs[iproc - 1] + recvcounts[iproc - 1]; - MPI_Allgatherv(qmIDs_mine,nqm_mine,MPI_LMP_TAGINT,qmIDs,recvcounts,displs, - MPI_LMP_TAGINT,world); + MPI_Allgatherv(qmIDs_mine, nqm_mine, MPI_LMP_TAGINT, qmIDs, recvcounts, displs, MPI_LMP_TAGINT, + world); memory->destroy(qmIDs_mine); memory->destroy(recvcounts); @@ -811,15 +815,15 @@ void FixMDIQM::create_qm_list() int *order; tagint *qmIDs_sort; - memory->create(order,nqm,"mdi/qm:order"); - memory->create(qmIDs_sort,nqm,"mdi/qm:qmIDs_sort"); + memory->create(order, nqm, "mdi/qm:order"); + memory->create(qmIDs_sort, nqm, "mdi/qm:qmIDs_sort"); for (int i = 0; i < nqm; i++) { qmIDs_sort[i] = qmIDs[i]; order[i] = i; } - utils::merge_sort(order,nqm,(void *) qmIDs_sort,compare_IDs); + utils::merge_sort(order, nqm, (void *) qmIDs_sort, compare_IDs); int j; for (int i = 0; i < nqm; i++) { @@ -827,7 +831,7 @@ void FixMDIQM::create_qm_list() qmIDs_sort[i] = qmIDs[j]; } - memcpy(qmIDs,qmIDs_sort,nqm*sizeof(tagint)); + memcpy(qmIDs, qmIDs_sort, nqm * sizeof(tagint)); memory->destroy(order); memory->destroy(qmIDs_sort); @@ -846,8 +850,10 @@ void FixMDIQM::set_qm2owned() for (int i = 0; i < nqm; i++) { index = atom->map(qmIDs[i]); - if (index >= nlocal) qm2owned[i] = -1; - else qm2owned[i] = index; + if (index >= nlocal) + qm2owned[i] = -1; + else + qm2owned[i] = index; } } @@ -902,7 +908,7 @@ void FixMDIQM::set_xqm() } } - MPI_Allreduce(&xqm_mine[0][0],&xqm[0][0],3*nqm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&xqm_mine[0][0], &xqm[0][0], 3 * nqm, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -919,7 +925,7 @@ void FixMDIQM::set_eqm() if (ilocal >= 0) eqm_mine[i] = elements[type[ilocal]]; } - MPI_Allreduce(eqm_mine,eqm,nqm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(eqm_mine, eqm, nqm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -936,7 +942,7 @@ void FixMDIQM::set_tqm() if (ilocal >= 0) tqm_mine[i] = type[ilocal]; } - MPI_Allreduce(tqm_mine,tqm,nqm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(tqm_mine, tqm, nqm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- @@ -1003,9 +1009,7 @@ void FixMDIQM::send_box() int ierr; // Only send the cell dimensions if this is a periodic simulation - if ( domain->xperiodic == 1 && - domain->yperiodic == 1 && - domain->zperiodic == 1 ) { + if (domain->xperiodic == 1 && domain->yperiodic == 1 && domain->zperiodic == 1) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); @@ -1019,17 +1023,14 @@ void FixMDIQM::send_box() ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL data"); - } - else { + } else { - if ( domain->xperiodic == 1 || - domain->yperiodic == 1 || - domain->zperiodic == 1 ) { - - error->all(FLERR, "MDI: The QM driver does not work for simulations with both periodic and non-periodic dimensions."); + if (domain->xperiodic == 1 || domain->yperiodic == 1 || domain->zperiodic == 1) { + error->all(FLERR, + "MDI: The QM driver does not work for simulations with both periodic and " + "non-periodic dimensions."); } - } } diff --git a/src/MDI/fix_mdi_qm.h b/src/MDI/fix_mdi_qm.h index fd31012cac..fff23560bb 100644 --- a/src/MDI/fix_mdi_qm.h +++ b/src/MDI/fix_mdi_qm.h @@ -50,18 +50,18 @@ class FixMDIQM : public Fix { int nprocs; int every, virialflag, addflag, connectflag; int plugin; - int sumflag,mcflag; + int sumflag, mcflag; char *id_mcfix; - int *mc_active_ptr,*exclusion_group_ptr; + int *mc_active_ptr, *exclusion_group_ptr; int *elements; - double qm_cell[9],qm_cell_displ[3]; + double qm_cell[9], qm_cell_displ[3]; double qm_energy; double qm_virial[9], qm_virial_symmetric[6]; MDI_Comm mdicomm; - int natoms_exists,celldispl_exists,elements_exists,types_exists; + int natoms_exists, celldispl_exists, elements_exists, types_exists; int stress_exists, pe_exists, keelec_exists; int nmax; @@ -77,15 +77,15 @@ class FixMDIQM : public Fix { // QM atom data structs - int nqm,nqm_last,max_nqm; + int nqm, nqm_last, max_nqm; int nexclude; tagint *qmIDs; int *qm2owned; - int *eqm,*eqm_mine; - int *tqm,*tqm_mine; - double **xqm,**xqm_mine; + int *eqm, *eqm_mine; + int *tqm, *tqm_mine; + double **xqm, **xqm_mine; double **fqm; // methods diff --git a/src/MDI/fix_mdi_qmmm.cpp b/src/MDI/fix_mdi_qmmm.cpp index a7c6defdc1..81d97dfbf5 100644 --- a/src/MDI/fix_mdi_qmmm.cpp +++ b/src/MDI/fix_mdi_qmmm.cpp @@ -44,11 +44,10 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // check requirements for LAMMPS to work with MDI for a QMMM engine // atom IDs do not need to be consecutive - if (!atom->tag_enable) - error->all(FLERR,"Fix mdi/qmmm requires atom IDs be defined"); + if (!atom->tag_enable) error->all(FLERR, "Fix mdi/qmmm requires atom IDs be defined"); if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Fix mdi/qmmm requires an atom map be defined"); + error->all(FLERR, "Fix mdi/qmmm requires an atom map be defined"); // confirm LAMMPS is being run as a driver @@ -59,9 +58,12 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // mode arg - if (strcmp(arg[3],"direct") == 0) mode = DIRECT; - else if (strcmp(arg[3],"potential") == 0) mode = POTENTIAL; - else error->all(FLERR,"Illegal fix mdi/qmmm command"); + if (strcmp(arg[3], "direct") == 0) + mode = DIRECT; + else if (strcmp(arg[3], "potential") == 0) + mode = POTENTIAL; + else + error->all(FLERR, "Illegal fix mdi/qmmm command"); // optional args @@ -92,30 +94,26 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg], "elements") == 0) { const char *symbols[] = { - "H" , "He", "Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne", - "Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", "K" , "Ca", - "Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", - "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y" , "Zr", - "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", - "Sb", "Te", "I" , "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", - "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", - "Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg", - "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", - "Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", - "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", - "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og", + "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", + "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", + "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", + "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", + "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", + "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", + "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", + "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og", }; int ntypes = atom->ntypes; if (iarg + ntypes + 1 > narg) error->all(FLERR, "Illegal fix mdi/qmmm command"); delete[] elements; elements = new int[ntypes + 1]; - for (int i = 1; i <= ntypes; i++) { + for (int i = 1; i <= ntypes; i++) { int anum; for (anum = 0; anum < MAXELEMENT; anum++) - if (strcmp(arg[iarg + i],symbols[anum]) == 0) break; + if (strcmp(arg[iarg + i], symbols[anum]) == 0) break; if (anum == MAXELEMENT) - error->all(FLERR,"Invalid chemical element in fix mdi/qmmm command"); + error->all(FLERR, "Invalid chemical element in fix mdi/qmmm command"); elements[i] = anum + 1; } iarg += ntypes + 1; @@ -210,7 +208,7 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // per-atom data nmax = atom->nmax; - memory->create(array_atom,nmax,3,"mdi/qmmm:array_atom"); + memory->create(array_atom, nmax, 3, "mdi/qmmm:array_atom"); // initialize outputs @@ -222,8 +220,7 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) sumflag = 0; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - array_atom[i][0] = array_atom[i][1] = array_atom[i][2] = 0.0; + for (int i = 0; i < nlocal; i++) array_atom[i][0] = array_atom[i][1] = array_atom[i][2] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -371,7 +368,7 @@ void FixMDIQMMM::init() // check if MDI engine supports commands that match QMMM mode if (mode == POTENTIAL) { - int check1,check2; + int check1, check2; ierr = MDI_Check_command_exists("@DEFAULT", ">POTENTIAL_AT_NUCLEI", mdicomm, &check1); if (ierr) error->all(FLERR, "MDI: >POTENTIAL_AT_NUCLEI command check"); @@ -382,11 +379,11 @@ void FixMDIQMMM::init() MPI_Bcast(&check2, 1, MPI_INT, 0, world); if (!check1 || !check2) - error->all(FLERR,"Fix mdi/qmmm potential mode not supported by MDI engine"); + error->all(FLERR, "Fix mdi/qmmm potential mode not supported by MDI engine"); } if (mode == DIRECT) { - int check1,check2,check3,check4; + int check1, check2, check3, check4; ierr = MDI_Check_command_exists("@DEFAULT", ">NLATTICE", mdicomm, &check1); if (ierr) error->all(FLERR, "MDI: >NLATTICE command check"); @@ -405,7 +402,7 @@ void FixMDIQMMM::init() MPI_Bcast(&check4, 1, MPI_INT, 0, world); if (!check1 || !check2 || !check3 || !check4) - error->all(FLERR,"Fix mdi/qmmm direct mode not supported by MDI engine"); + error->all(FLERR, "Fix mdi/qmmm direct mode not supported by MDI engine"); ierr = MDI_Check_command_exists("@DEFAULT", ">LATTICE_ELEMENTS", mdicomm, &check1); if (ierr) error->all(FLERR, "MDI: >LATTICE_ELEMENTS command check"); @@ -416,9 +413,9 @@ void FixMDIQMMM::init() MPI_Bcast(&check2, 1, MPI_INT, 0, world); if (elements_exists && !check1) - error->all(FLERR,"Fix mdi/qmmm direct mode elements not supported by MDI engine"); + error->all(FLERR, "Fix mdi/qmmm direct mode elements not supported by MDI engine"); if (types_exists && !check2) - error->all(FLERR,"Fix mdi/qmmm direct mode types not supported by MDI engine"); + error->all(FLERR, "Fix mdi/qmmm direct mode types not supported by MDI engine"); } } @@ -426,14 +423,15 @@ void FixMDIQMMM::init() // POTENTIAL mode requires pair style to calculate only Coulombic interactions // can be in conjunction with KSpace solver - if (!atom->q_flag) error->all(FLERR,"Fix mdi/qmmm requires per-atom charge"); + if (!atom->q_flag) error->all(FLERR, "Fix mdi/qmmm requires per-atom charge"); if (mode == POTENTIAL) { - if (!force->pair) error->all(FLERR,"Fix mdi/qmmm potential requires a pair style"); - pair_coul = force->pair_match("coul/cut",1,0); - if (!pair_coul) pair_coul = force->pair_match("coul/long",1,0); - if (!pair_coul) pair_coul = force->pair_match("coul/msm",1,0); - if (!pair_coul) error->all(FLERR,"Fix mdi/qmmm potential requires Coulomb-only pair sub-style"); + if (!force->pair) error->all(FLERR, "Fix mdi/qmmm potential requires a pair style"); + pair_coul = force->pair_match("coul/cut", 1, 0); + if (!pair_coul) pair_coul = force->pair_match("coul/long", 1, 0); + if (!pair_coul) pair_coul = force->pair_match("coul/msm", 1, 0); + if (!pair_coul) + error->all(FLERR, "Fix mdi/qmmm potential requires Coulomb-only pair sub-style"); } // determine whether a new vs incremental QMMM calc is needed @@ -468,11 +466,12 @@ void FixMDIQMMM::init() // check if box has changed - if (new_system) set_box(); + if (new_system) + set_box(); else { - double old_cell[9],old_cell_displ[3]; - memcpy(old_cell,qm_cell,9*sizeof(double)); - memcpy(old_cell_displ,qm_cell_displ,3*sizeof(double)); + double old_cell[9], old_cell_displ[3]; + memcpy(old_cell, qm_cell, 9 * sizeof(double)); + memcpy(old_cell_displ, qm_cell_displ, 3 * sizeof(double)); set_box(); for (int i = 0; i < 9; i++) if (qm_cell[i] != old_cell[i]) new_system = 1; @@ -488,8 +487,8 @@ void FixMDIQMMM::init() if (mode == DIRECT) set_emm(); } else { int *eqm_old; - memory->create(eqm_old,nqm,"mdi/qmmm:eqm_old"); - memcpy(eqm_old,eqm,nqm*sizeof(int)); + memory->create(eqm_old, nqm, "mdi/qmmm:eqm_old"); + memcpy(eqm_old, eqm, nqm * sizeof(int)); set_eqm(); for (int i = 0; i < nqm; i++) if (eqm[i] != eqm_old[i]) new_system = 1; @@ -497,8 +496,8 @@ void FixMDIQMMM::init() if (mode == DIRECT) { int *emm_old; - memory->create(emm_old,nmm,"mdi/qmmm:emm_old"); - memcpy(emm_old,emm,nmm*sizeof(int)); + memory->create(emm_old, nmm, "mdi/qmmm:emm_old"); + memcpy(emm_old, emm, nmm * sizeof(int)); set_emm(); for (int i = 0; i < nmm; i++) if (emm[i] != emm_old[i]) new_system = 1; @@ -512,8 +511,8 @@ void FixMDIQMMM::init() if (mode == DIRECT) set_tmm(); } else { int *tqm_old; - memory->create(tqm_old,nqm,"mdi/qmmm:tqm_old"); - memcpy(tqm_old,tqm,nqm*sizeof(int)); + memory->create(tqm_old, nqm, "mdi/qmmm:tqm_old"); + memcpy(tqm_old, tqm, nqm * sizeof(int)); set_tqm(); for (int i = 0; i < nqm; i++) if (tqm[i] != tqm_old[i]) new_system = 1; @@ -521,8 +520,8 @@ void FixMDIQMMM::init() if (mode == DIRECT) { int *tmm_old; - memory->create(tmm_old,nmm,"mdi/qmmm:tmm_old"); - memcpy(tmm_old,tmm,nmm*sizeof(int)); + memory->create(tmm_old, nmm, "mdi/qmmm:tmm_old"); + memcpy(tmm_old, tmm, nmm * sizeof(int)); set_tmm(); for (int i = 0; i < nmm; i++) if (tmm[i] != tmm_old[i]) new_system = 1; @@ -538,14 +537,18 @@ void FixMDIQMMM::init() if (new_system) { send_natoms_qm(); send_box(); - if (elements && elements_exists) send_elements_qm(); - else if (types_exists) send_types_qm(); + if (elements && elements_exists) + send_elements_qm(); + else if (types_exists) + send_types_qm(); nqm_last = nqm; if (mode == DIRECT) { send_natoms_mm(); - if (elements && elements_exists) send_elements_mm(); - else if (types_exists) send_types_mm(); + if (elements && elements_exists) + send_elements_mm(); + else if (types_exists) + send_types_mm(); set_qmm(); send_charges_mm(); } @@ -597,12 +600,12 @@ void FixMDIQMMM::pre_force(int vflag) // invoke pair hybrid sub-style pair coulomb and Kspace directly // set eflag = 2 so they calculate per-atom energy - pair_coul->compute(2,0); + pair_coul->compute(2, 0); double *eatom_pair = pair_coul->eatom; double *eatom_kspace = nullptr; if (force->kspace) { - force->kspace->compute(2,0); + force->kspace->compute(2, 0); eatom_kspace = force->kspace->eatom; } @@ -612,7 +615,7 @@ void FixMDIQMMM::pre_force(int vflag) if (atom->nmax > ncoulmax) { memory->destroy(ecoul); ncoulmax = atom->nmax; - memory->create(ecoul,ncoulmax,"mdi/qmmm:ecoul"); + memory->create(ecoul, ncoulmax, "mdi/qmmm:ecoul"); } // ecoul = per-atom energy for my owned atoms @@ -622,14 +625,12 @@ void FixMDIQMMM::pre_force(int vflag) int ntotal = nlocal; if (force->newton_pair) ntotal += atom->nghost; - for (int i = 0; i < ntotal; i++) - ecoul[i] = eatom_pair[i]; + for (int i = 0; i < ntotal; i++) ecoul[i] = eatom_pair[i]; if (force->newton_pair) comm->reverse_comm(this); if (force->kspace) { - for (int i = 0; i < nlocal; i++) - ecoul[i] += eatom_kspace[i]; + for (int i = 0; i < nlocal; i++) ecoul[i] += eatom_kspace[i]; } // setup QM inputs: xqm and qpotential @@ -652,8 +653,10 @@ void FixMDIQMMM::pre_force(int vflag) for (int i = 0; i < nqm; i++) { ilocal = qm2owned[i]; if (ilocal >= 0) { - if (q[ilocal] == 0.0) qpotential_mine[i] = 0.0; - else qpotential_mine[i] = 2.0 * ecoul[ilocal] / q[ilocal]; + if (q[ilocal] == 0.0) + qpotential_mine[i] = 0.0; + else + qpotential_mine[i] = 2.0 * ecoul[ilocal] / q[ilocal]; } } @@ -666,13 +669,13 @@ void FixMDIQMMM::pre_force(int vflag) double dely = xqm[i][1] - xqm[j][1]; double delz = xqm[i][2] - xqm[j][2]; domain->minimum_image(delx, dely, delz); - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; qpotential_mine[i] -= qqrd2e * qqm[j] / sqrt(rsq); } } } - MPI_Allreduce(qpotential_mine,qpotential,nqm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(qpotential_mine, qpotential, nqm, MPI_DOUBLE, MPI_SUM, world); // unit conversion from LAMMPS to MDI // must be done here, rather than in set_xqm() @@ -829,8 +832,10 @@ void FixMDIQMMM::pre_force(int vflag) void FixMDIQMMM::post_force(int vflag) { - if (mode == DIRECT) post_force_direct(vflag); - else if (mode == POTENTIAL) post_force_potential(vflag); + if (mode == DIRECT) + post_force_direct(vflag); + else if (mode == POTENTIAL) + post_force_potential(vflag); } /* ---------------------------------------------------------------------- @@ -1121,10 +1126,9 @@ void FixMDIQMMM::min_post_force(int vflag) /* ---------------------------------------------------------------------- */ -int FixMDIQMMM::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixMDIQMMM::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - int i,j,m; + int i, j, m; double *q = atom->q; @@ -1140,7 +1144,7 @@ int FixMDIQMMM::pack_forward_comm(int n, int *list, double *buf, void FixMDIQMMM::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -1154,7 +1158,7 @@ void FixMDIQMMM::unpack_forward_comm(int n, int first, double *buf) int FixMDIQMMM::pack_reverse_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -1166,7 +1170,7 @@ int FixMDIQMMM::pack_reverse_comm(int n, int first, double *buf) void FixMDIQMMM::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i, j, m; m = 0; for (i = 0; i < n; i++) { @@ -1200,10 +1204,10 @@ double FixMDIQMMM::compute_vector(int n) double FixMDIQMMM::memory_usage() { double bytes = 0.0; - bytes += (double)ncoulmax * sizeof(double); - bytes += (double)(3*3 + 4) * nqm * sizeof(double); // fpoint QM arrays/vecs - bytes += nqm * sizeof(tagint); // qmIDs - bytes += nqm * sizeof(int); // qm2owned + bytes += (double) ncoulmax * sizeof(double); + bytes += (double) (3 * 3 + 4) * nqm * sizeof(double); // fpoint QM arrays/vecs + bytes += nqm * sizeof(tagint); // qmIDs + bytes += nqm * sizeof(int); // qm2owned return bytes; } @@ -1237,21 +1241,21 @@ void FixMDIQMMM::reallocate_qm() memory->destroy(qqm_mine); memory->destroy(qpotential_mine); - memory->create(qmIDs,max_nqm,"mdi/qmmm:qmIDs"); - memory->create(qm2owned,max_nqm,"mdi/qmmm:qm2owned"); + memory->create(qmIDs, max_nqm, "mdi/qmmm:qmIDs"); + memory->create(qm2owned, max_nqm, "mdi/qmmm:qm2owned"); - memory->create(eqm,max_nqm,"mdi/qmmm:eqm"); - memory->create(tqm,max_nqm,"mdi/qmmm:tqm"); - memory->create(xqm,max_nqm,3,"mdi/qmmm:xqm"); - memory->create(fqm,max_nqm,3,"mdi/qmmm:fqm"); - memory->create(qqm,max_nqm,"mdi/qmmm:qqm"); - memory->create(qpotential,max_nqm,"mdi/qmmm:qpotential"); + memory->create(eqm, max_nqm, "mdi/qmmm:eqm"); + memory->create(tqm, max_nqm, "mdi/qmmm:tqm"); + memory->create(xqm, max_nqm, 3, "mdi/qmmm:xqm"); + memory->create(fqm, max_nqm, 3, "mdi/qmmm:fqm"); + memory->create(qqm, max_nqm, "mdi/qmmm:qqm"); + memory->create(qpotential, max_nqm, "mdi/qmmm:qpotential"); - memory->create(eqm_mine,max_nqm,"mdi/qmmm:eqm_mine"); - memory->create(tqm_mine,max_nqm,"mdi/qmmm:tqm_mine"); - memory->create(xqm_mine,max_nqm,3,"mdi/qmmm:xqm_mine"); - memory->create(qqm_mine,max_nqm,"mdi/qmmm:qqm_mine"); - memory->create(qpotential_mine,max_nqm,"mdi/qmmm:qpotential_mine"); + memory->create(eqm_mine, max_nqm, "mdi/qmmm:eqm_mine"); + memory->create(tqm_mine, max_nqm, "mdi/qmmm:tqm_mine"); + memory->create(xqm_mine, max_nqm, 3, "mdi/qmmm:xqm_mine"); + memory->create(qqm_mine, max_nqm, "mdi/qmmm:qqm_mine"); + memory->create(qpotential_mine, max_nqm, "mdi/qmmm:qpotential_mine"); } /* ---------------------------------------------------------------------- @@ -1276,19 +1280,19 @@ void FixMDIQMMM::reallocate_mm() memory->destroy(xmm_mine); memory->destroy(qmm_mine); - memory->create(mmIDs,max_nmm,"mdi/qmmm:mmIDs"); - memory->create(mm2owned,max_nmm,"mdi/qmmm:mm2owned"); + memory->create(mmIDs, max_nmm, "mdi/qmmm:mmIDs"); + memory->create(mm2owned, max_nmm, "mdi/qmmm:mm2owned"); - memory->create(emm,max_nmm,"mdi/qmmm:emm"); - memory->create(tmm,max_nmm,"mdi/qmmm:tmm"); - memory->create(xmm,max_nmm,3,"mdi/qmmm:xmm"); - memory->create(qmm,max_nmm,"mdi/qmmm:qmm"); - memory->create(fmm,max_nmm,3,"mdi/qmmm:fmm"); + memory->create(emm, max_nmm, "mdi/qmmm:emm"); + memory->create(tmm, max_nmm, "mdi/qmmm:tmm"); + memory->create(xmm, max_nmm, 3, "mdi/qmmm:xmm"); + memory->create(qmm, max_nmm, "mdi/qmmm:qmm"); + memory->create(fmm, max_nmm, 3, "mdi/qmmm:fmm"); - memory->create(emm_mine,max_nmm,"mdi/qmmm:emm_mine"); - memory->create(tmm_mine,max_nmm,"mdi/qmmm:tmm_mine"); - memory->create(xmm_mine,max_nmm,3,"mdi/qmmm:xmm_mine"); - memory->create(qmm_mine,max_nmm,"mdi/qmmm:qmm_mine"); + memory->create(emm_mine, max_nmm, "mdi/qmmm:emm_mine"); + memory->create(tmm_mine, max_nmm, "mdi/qmmm:tmm_mine"); + memory->create(xmm_mine, max_nmm, 3, "mdi/qmmm:xmm_mine"); + memory->create(qmm_mine, max_nmm, "mdi/qmmm:qmm_mine"); } /* ---------------------------------------------------------------------- @@ -1302,14 +1306,13 @@ int FixMDIQMMM::set_nqm() // require 3*nqm be a small INT, so can MPI_Allreduce QM values - if (3*ngroup > MAXSMALLINT) - error->all(FLERR,"Fix mdi/qmmm has too many quantum atoms"); + if (3 * ngroup > MAXSMALLINT) error->all(FLERR, "Fix mdi/qmmm has too many quantum atoms"); // error if nqm = 0 // error if nqm = natoms, should use fix mdi/qm instead - if (ngroup == 0) error->all(FLERR,"Fix mdi/qmmm has no atoms in quantum group"); - if (ngroup == atom->natoms) error->all(FLERR,"Fix mdi/qmmm has all atoms in quantum group"); + if (ngroup == 0) error->all(FLERR, "Fix mdi/qmmm has no atoms in quantum group"); + if (ngroup == atom->natoms) error->all(FLERR, "Fix mdi/qmmm has all atoms in quantum group"); int ncount = ngroup; return ncount; @@ -1324,8 +1327,8 @@ int FixMDIQMMM::set_nmm() { // require 3*nmm be a small INT, so can MPI_Allreduce xmm - if (3*(atom->natoms-nqm) > MAXSMALLINT) - error->all(FLERR,"Fix mdi/qmmm has too many classical atoms"); + if (3 * (atom->natoms - nqm) > MAXSMALLINT) + error->all(FLERR, "Fix mdi/qmmm has too many classical atoms"); int ncount = atom->natoms - nqm; return ncount; @@ -1351,24 +1354,24 @@ void FixMDIQMMM::create_qm_list() if (mask[i] & groupbit) nqm_mine++; tagint *qmIDs_mine; - memory->create(qmIDs_mine,nqm_mine,"mdi/qmmm:qmIDs_mine"); + memory->create(qmIDs_mine, nqm_mine, "mdi/qmmm:qmIDs_mine"); nqm_mine = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) qmIDs_mine[nqm_mine++] = tag[i]; int *recvcounts, *displs; - memory->create(recvcounts,nprocs,"mdi/qmmm:recvcounts"); - memory->create(displs,nprocs,"mdi/qmmm:displs"); + memory->create(recvcounts, nprocs, "mdi/qmmm:recvcounts"); + memory->create(displs, nprocs, "mdi/qmmm:displs"); - MPI_Allgather(&nqm_mine,1,MPI_INT,recvcounts,1,MPI_INT,world); + MPI_Allgather(&nqm_mine, 1, MPI_INT, recvcounts, 1, MPI_INT, world); displs[0] = 0; for (int iproc = 1; iproc < nprocs; iproc++) - displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; + displs[iproc] = displs[iproc - 1] + recvcounts[iproc - 1]; - MPI_Allgatherv(qmIDs_mine,nqm_mine,MPI_LMP_TAGINT,qmIDs,recvcounts,displs, - MPI_LMP_TAGINT,world); + MPI_Allgatherv(qmIDs_mine, nqm_mine, MPI_LMP_TAGINT, qmIDs, recvcounts, displs, MPI_LMP_TAGINT, + world); memory->destroy(qmIDs_mine); memory->destroy(recvcounts); @@ -1379,15 +1382,15 @@ void FixMDIQMMM::create_qm_list() int *order; tagint *qmIDs_sort; - memory->create(order,nqm,"mdi/qmmm:order"); - memory->create(qmIDs_sort,nqm,"mdi/qmmm:qmIDs_sort"); + memory->create(order, nqm, "mdi/qmmm:order"); + memory->create(qmIDs_sort, nqm, "mdi/qmmm:qmIDs_sort"); for (int i = 0; i < nqm; i++) { qmIDs_sort[i] = qmIDs[i]; order[i] = i; } - utils::merge_sort(order,nqm,(void *) qmIDs_sort,compare_IDs); + utils::merge_sort(order, nqm, (void *) qmIDs_sort, compare_IDs); int j; for (int i = 0; i < nqm; i++) { @@ -1395,7 +1398,7 @@ void FixMDIQMMM::create_qm_list() qmIDs_sort[i] = qmIDs[j]; } - memcpy(qmIDs,qmIDs_sort,nqm*sizeof(tagint)); + memcpy(qmIDs, qmIDs_sort, nqm * sizeof(tagint)); memory->destroy(order); memory->destroy(qmIDs_sort); @@ -1415,27 +1418,27 @@ void FixMDIQMMM::create_mm_list() int nmm_mine = 0; for (int i = 0; i < nlocal; i++) - if (!(mask[i] & groupbit)) nmm_mine++; + if (!(mask[i] & groupbit)) nmm_mine++; tagint *mmIDs_mine; - memory->create(mmIDs_mine,nmm_mine,"mdi/qmmm:mmIDs_mine"); + memory->create(mmIDs_mine, nmm_mine, "mdi/qmmm:mmIDs_mine"); nmm_mine = 0; for (int i = 0; i < nlocal; i++) if (!(mask[i] & groupbit)) mmIDs_mine[nmm_mine++] = tag[i]; int *recvcounts, *displs; - memory->create(recvcounts,nprocs,"mdi/qmmm:recvcounts"); - memory->create(displs,nprocs,"mdi/qmmm:displs"); + memory->create(recvcounts, nprocs, "mdi/qmmm:recvcounts"); + memory->create(displs, nprocs, "mdi/qmmm:displs"); - MPI_Allgather(&nmm_mine,1,MPI_INT,recvcounts,1,MPI_INT,world); + MPI_Allgather(&nmm_mine, 1, MPI_INT, recvcounts, 1, MPI_INT, world); displs[0] = 0; for (int iproc = 1; iproc < nprocs; iproc++) - displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; + displs[iproc] = displs[iproc - 1] + recvcounts[iproc - 1]; - MPI_Allgatherv(mmIDs_mine,nmm_mine,MPI_LMP_TAGINT,mmIDs,recvcounts,displs, - MPI_LMP_TAGINT,world); + MPI_Allgatherv(mmIDs_mine, nmm_mine, MPI_LMP_TAGINT, mmIDs, recvcounts, displs, MPI_LMP_TAGINT, + world); memory->destroy(mmIDs_mine); memory->destroy(recvcounts); @@ -1446,15 +1449,15 @@ void FixMDIQMMM::create_mm_list() int *order; tagint *mmIDs_sort; - memory->create(order,nmm,"mdi/qmmm:order"); - memory->create(mmIDs_sort,nmm,"mdi/qmmm:mmIDs_sort"); + memory->create(order, nmm, "mdi/qmmm:order"); + memory->create(mmIDs_sort, nmm, "mdi/qmmm:mmIDs_sort"); for (int i = 0; i < nmm; i++) { mmIDs_sort[i] = mmIDs[i]; order[i] = i; } - utils::merge_sort(order,nmm,(void *) mmIDs_sort,compare_IDs); + utils::merge_sort(order, nmm, (void *) mmIDs_sort, compare_IDs); int j; for (int i = 0; i < nmm; i++) { @@ -1462,7 +1465,7 @@ void FixMDIQMMM::create_mm_list() mmIDs_sort[i] = mmIDs[j]; } - memcpy(mmIDs,mmIDs_sort,nmm*sizeof(tagint)); + memcpy(mmIDs, mmIDs_sort, nmm * sizeof(tagint)); memory->destroy(order); memory->destroy(mmIDs_sort); @@ -1481,8 +1484,10 @@ void FixMDIQMMM::set_qm2owned() for (int i = 0; i < nqm; i++) { index = atom->map(qmIDs[i]); - if (index >= nlocal) qm2owned[i] = -1; - else qm2owned[i] = index; + if (index >= nlocal) + qm2owned[i] = -1; + else + qm2owned[i] = index; } } @@ -1499,8 +1504,10 @@ void FixMDIQMMM::set_mm2owned() for (int i = 0; i < nmm; i++) { index = atom->map(mmIDs[i]); - if (index >= nlocal) mm2owned[i] = -1; - else mm2owned[i] = index; + if (index >= nlocal) + mm2owned[i] = -1; + else + mm2owned[i] = index; } } @@ -1527,7 +1534,6 @@ void FixMDIQMMM::set_box() for (int icell = 0; icell < 9; icell++) qm_cell[icell] *= lmp2mdi_length; } - /* ---------------------------------------------------------------------- fill xqm with QM atom coords if convert, perform LAMMPS to MDI unit conversion @@ -1562,7 +1568,7 @@ void FixMDIQMMM::set_xqm(int convert) } } - MPI_Allreduce(&xqm_mine[0][0],&xqm[0][0],3*nqm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&xqm_mine[0][0], &xqm[0][0], 3 * nqm, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1579,7 +1585,7 @@ void FixMDIQMMM::set_eqm() if (ilocal >= 0) eqm_mine[i] = elements[type[ilocal]]; } - MPI_Allreduce(eqm_mine,eqm,nqm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(eqm_mine, eqm, nqm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1596,7 +1602,7 @@ void FixMDIQMMM::set_tqm() if (ilocal >= 0) tqm_mine[i] = type[ilocal]; } - MPI_Allreduce(tqm_mine,tqm,nqm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(tqm_mine, tqm, nqm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1613,7 +1619,7 @@ void FixMDIQMMM::set_qqm() if (ilocal >= 0) qqm_mine[i] = q[ilocal]; } - MPI_Allreduce(qqm_mine,qqm,nqm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(qqm_mine, qqm, nqm, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1644,7 +1650,7 @@ void FixMDIQMMM::set_xmm() } } - MPI_Allreduce(&xmm_mine[0][0],&xmm[0][0],3*nmm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&xmm_mine[0][0], &xmm[0][0], 3 * nmm, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1661,7 +1667,7 @@ void FixMDIQMMM::set_emm() if (ilocal >= 0) emm_mine[i] = elements[type[ilocal]]; } - MPI_Allreduce(emm_mine,emm,nmm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(emm_mine, emm, nmm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1678,7 +1684,7 @@ void FixMDIQMMM::set_tmm() if (ilocal >= 0) tmm_mine[i] = type[ilocal]; } - MPI_Allreduce(tmm_mine,tmm,nmm,MPI_INT,MPI_SUM,world); + MPI_Allreduce(tmm_mine, tmm, nmm, MPI_INT, MPI_SUM, world); } /* ---------------------------------------------------------------------- */ @@ -1695,7 +1701,7 @@ void FixMDIQMMM::set_qmm() if (ilocal >= 0) qmm_mine[i] = q[ilocal]; } - MPI_Allreduce(qmm_mine,qmm,nmm,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(qmm_mine, qmm, nmm, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- @@ -1715,7 +1721,6 @@ void FixMDIQMMM::send_natoms_qm() ierr = MDI_Send(&nqm, 1, MDI_INT, mdicomm); if (ierr) error->all(FLERR, "MDI: >NATOMS data"); - } else { ierr = MDI_Send_command("all(FLERR, "MDI: xperiodic == 1 && - domain->yperiodic == 1 && - domain->zperiodic == 1 ) { + if (domain->xperiodic == 1 && domain->yperiodic == 1 && domain->zperiodic == 1) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); @@ -1779,19 +1782,15 @@ void FixMDIQMMM::send_box() ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL data"); - } - else { + } else { - if ( domain->xperiodic == 1 || - domain->yperiodic == 1 || - domain->zperiodic == 1 ) { - - error->all(FLERR, "MDI: The QM driver does not work for simulations with both periodic and non-periodic dimensions."); + if (domain->xperiodic == 1 || domain->yperiodic == 1 || domain->zperiodic == 1) { + error->all(FLERR, + "MDI: The QM driver does not work for simulations with both periodic and " + "non-periodic dimensions."); } - } - } /* ---------------------------------------------------------------------- @@ -1818,8 +1817,7 @@ void FixMDIQMMM::get_qm_energy() if (ierr) error->all(FLERR, "MDI: all(FLERR, "MDI: Date: Wed, 19 Apr 2023 10:22:10 -0600 Subject: [PATCH 6/8] cosmetic changes to comments and code structure --- doc/src/fix_mdi_qm.rst | 5 +- doc/src/fix_mdi_qmmm.rst | 5 +- src/MDI/fix_mdi_qm.cpp | 84 +++++++++++++++++--------------- src/MDI/fix_mdi_qm.h | 2 + src/MDI/fix_mdi_qmmm.cpp | 101 +++++++++++++++++++-------------------- src/MDI/fix_mdi_qmmm.h | 2 +- src/MDI/mdi_engine.cpp | 7 ++- 7 files changed, 109 insertions(+), 97 deletions(-) diff --git a/doc/src/fix_mdi_qm.rst b/doc/src/fix_mdi_qm.rst index 007c2c275a..6dcf2cd1e8 100644 --- a/doc/src/fix_mdi_qm.rst +++ b/doc/src/fix_mdi_qm.rst @@ -294,7 +294,10 @@ LAMMPS can also be used as an MDI driver in other unit choices it supports, e.g. *lj*, but then no unit conversion to MDI units is performed. -If this fix is used in conjuction with a QM code that does not support periodic boundary conditions (more specifically, a QM code that does not support the ``>CELL`` MDI command), the LAMMPS input file must specify non-periodic boundary conditions. +If this fix is used in conjuction with a QM code that does not support +periodic boundary conditions (more specifically, a QM code that does +not support the ``>CELL`` MDI command), the LAMMPS system must be +fully non-periodic. I.e. no dimension of the system can be periodic. Related commands """""""""""""""" diff --git a/doc/src/fix_mdi_qmmm.rst b/doc/src/fix_mdi_qmmm.rst index 4f81fbb999..16e1c1eaae 100644 --- a/doc/src/fix_mdi_qmmm.rst +++ b/doc/src/fix_mdi_qmmm.rst @@ -261,7 +261,10 @@ used to specify *real* or *metal* units. This will ensure the correct unit conversions between LAMMPS and MDI units. The other code will also perform similar unit conversions into its preferred units. -If this fix is used in conjuction with a QM code that does not support periodic boundary conditions (more specifically, a QM code that does not support the ``>CELL`` MDI command), the LAMMPS input file must specify non-periodic boundary conditions. +If this fix is used in conjuction with a QM code that does not support +periodic boundary conditions (more specifically, a QM code that does +not support the ``>CELL`` MDI command), the LAMMPS system must be +fully non-periodic. I.e. no dimension of the system can be periodic. Related commands """""""""""""""" diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index 79d85ebced..bdcd5daacf 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/src/MDI/fix_mdi_qm.cpp @@ -512,35 +512,11 @@ void FixMDIQM::post_force(int vflag) ierr = MDI_Send(&xqm[0][0], 3 * nqm, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >COORDS data"); - // request potential energy from MDI engine + // request QM energy from MDI engine // this triggers engine to perform QM calculation - // qm_energy = fix output for global QM energy + // sets qm_energy = fix output for global QM energy - if (pe_exists && keelec_exists) { - int pe_energy, keelec_energy; - - // get the total potential energy - ierr = MDI_Send_command("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: xperiodic == 1 && domain->yperiodic == 1 && domain->zperiodic == 1) { - + // only send cell dimensions if fully periodic simulation + + if (domain->nonperiodic == 0) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); @@ -1023,17 +999,49 @@ void FixMDIQM::send_box() ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL data"); - } else { - - if (domain->xperiodic == 1 || domain->yperiodic == 1 || domain->zperiodic == 1) { - - error->all(FLERR, - "MDI: The QM driver does not work for simulations with both periodic and " - "non-periodic dimensions."); - } + } else if (domain->xperiodic == 1 || domain->yperiodic == 1 || + domain->zperiodic == 1) { + error->all(FLERR,"MDI requires fully periodic or fully non-periodic system"); } } +/* ---------------------------------------------------------------------- + request QM energy from MDI engine + set qm_energy = fix output for global QM energy +------------------------------------------------------------------------- */ + +void FixMDIQM::request_qm_energy() +{ + int ierr; + + // QM energy = all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: >POTENTIAL_AT_NUCLEI data"); - // request QM potential energy from MDI engine + // request QM energy from MDI engine // this triggers engine to perform QM calculation - // qm_energy = fix output for global QM energy + // sets qm_energy = fix output for global QM energy - get_qm_energy(); + request_qm_energy(); // request forces on QM atoms from MDI engine @@ -886,11 +886,11 @@ void FixMDIQMMM::post_force_direct(int vflag) ierr = MDI_Send(&xmm[0][0], 3 * nmm, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >CLATTICE data"); - // request QM potential energy from MDI engine + // request QM energy from MDI engine // this triggers engine to perform QM calculation - // qm_energy = fix output for global QM energy + // sets qm_energy = fix output for global QM energy - get_qm_energy(); + request_qm_energy(); // request forces on QM atoms from MDI engine @@ -1767,9 +1767,9 @@ void FixMDIQMMM::send_box() { int ierr; - // Only send the cell dimensions if this is a periodic simulation - if (domain->xperiodic == 1 && domain->yperiodic == 1 && domain->zperiodic == 1) { - + // only send cell dimensions if fully periodic simulation + + if (domain->nonperiodic == 0) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command"); @@ -1782,52 +1782,12 @@ void FixMDIQMMM::send_box() ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm); if (ierr) error->all(FLERR, "MDI: >CELL data"); - } else { - - if (domain->xperiodic == 1 || domain->yperiodic == 1 || domain->zperiodic == 1) { - - error->all(FLERR, - "MDI: The QM driver does not work for simulations with both periodic and " - "non-periodic dimensions."); - } + } else if (domain->xperiodic == 1 || domain->yperiodic == 1 || + domain->zperiodic == 1) { + error->all(FLERR,"MDI requires fully periodic or fully non-periodic system"); } } -/* ---------------------------------------------------------------------- - get the energy of the QM system -------------------------------------------------------------------------- */ - -void FixMDIQMMM::get_qm_energy() -{ - int ierr; - - if (pe_exists && keelec_exists) { - int pe_energy, keelec_energy; - - // get the total potential energy - ierr = MDI_Send_command("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: >LATTICE data"); } +/* ---------------------------------------------------------------------- + request QM energy from MDI engine + set qm_energy = fix output for global QM energy +------------------------------------------------------------------------- */ + +void FixMDIQMMM::request_qm_energy() +{ + int ierr; + + // QM energy = all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: Date: Wed, 19 Apr 2023 13:27:29 -0400 Subject: [PATCH 7/8] whitespace --- src/MDI/fix_mdi_qm.cpp | 8 ++++---- src/MDI/fix_mdi_qm.h | 2 +- src/MDI/fix_mdi_qmmm.cpp | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index bdcd5daacf..9c47d9c05a 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/src/MDI/fix_mdi_qm.cpp @@ -985,7 +985,7 @@ void FixMDIQM::send_box() int ierr; // only send cell dimensions if fully periodic simulation - + if (domain->nonperiodic == 0) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); @@ -1018,7 +1018,7 @@ void FixMDIQM::request_qm_energy() if (pe_exists && keelec_exists) { int pe_energy, keelec_energy; - + ierr = MDI_Send_command("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: nonperiodic == 0) { if (celldispl_exists) { ierr = MDI_Send_command(">CELL_DISPL", mdicomm); @@ -1849,7 +1849,7 @@ void FixMDIQMMM::request_qm_energy() if (pe_exists && keelec_exists) { int pe_energy, keelec_energy; - + ierr = MDI_Send_command("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: Date: Wed, 19 Apr 2023 13:28:02 -0400 Subject: [PATCH 8/8] fix typo --- examples/COUPLE/library/README | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/COUPLE/library/README b/examples/COUPLE/library/README index 03ec2379d0..aa94eed8aa 100644 --- a/examples/COUPLE/library/README +++ b/examples/COUPLE/library/README @@ -6,7 +6,7 @@ SPPARKS applications in 2 sister directories. The library dir has a Makefile (which you may need to edit for your box). If you type -g++ -f Makefile.mpi +make -f Makefile.mpi you should create libcouple.a, which the other coupled applications link to.