diff --git a/cmake/gmxManageMimic.cmake b/cmake/gmxManageMimic.cmake index c0c98980bf..b5d28965b3 100644 --- a/cmake/gmxManageMimic.cmake +++ b/cmake/gmxManageMimic.cmake @@ -33,7 +33,7 @@ # the research papers on the package. Check out http://www.gromacs.org. if(GMX_MIMIC) - find_package(MimicCommLib REQUIRED) - include_directories(SYSTEM ${MimicCommLib_INCLUDE_DIRS}) - list(APPEND GMX_COMMON_LIBRARIES Upstream::mimiccomm) + find_package(MiMiCCommLib REQUIRED CXX) + include_directories(SYSTEM ${MiMiCCommLib_INCLUDE_DIRS}) + list(APPEND GMX_COMMON_LIBRARIES ${MiMiCCommLib_LIBRARIES}) endif() diff --git a/src/gromacs/gmxlib/network.cpp b/src/gromacs/gmxlib/network.cpp index 66023edac8..6b83ce7c03 100644 --- a/src/gromacs/gmxlib/network.cpp +++ b/src/gromacs/gmxlib/network.cpp @@ -114,12 +114,12 @@ void done_commrec(t_commrec* cr) // structure of the commrec and domdec initialization code makes // it hard to avoid both leaks and double frees. bool mySimIsMyGroup = (cr->mpi_comm_mysim == cr->mpi_comm_mygroup); - if (cr->mpi_comm_mysim != MPI_COMM_NULL && cr->mpi_comm_mysim != MPI_COMM_WORLD) + if (cr->mpi_comm_mysim != MPI_COMM_NULL && cr->mpi_comm_mysim != MPI_WORLD_STUB) { // TODO see above // MPI_Comm_free(&cr->mpi_comm_mysim); } - if (!mySimIsMyGroup && cr->mpi_comm_mygroup != MPI_COMM_NULL && cr->mpi_comm_mygroup != MPI_COMM_WORLD) + if (!mySimIsMyGroup && cr->mpi_comm_mygroup != MPI_COMM_NULL && cr->mpi_comm_mygroup != MPI_WORLD_STUB) { // TODO see above // MPI_Comm_free(&cr->mpi_comm_mygroup); @@ -485,7 +485,7 @@ void gmx_fatal_collective(int f_errno, #if GMX_MPI int result; /* Check if we are calling on all processes in MPI_COMM_WORLD */ - MPI_Comm_compare(comm, MPI_COMM_WORLD, &result); + MPI_Comm_compare(comm, MPI_WORLD_STUB, &result); /* Any result except MPI_UNEQUAL allows us to call MPI_Finalize */ bFinalize = (result != MPI_UNEQUAL); #else diff --git a/src/gromacs/hardware/detecthardware.cpp b/src/gromacs/hardware/detecthardware.cpp index d808249678..8b2aa4cd3b 100644 --- a/src/gromacs/hardware/detecthardware.cpp +++ b/src/gromacs/hardware/detecthardware.cpp @@ -252,7 +252,7 @@ static void gmx_collect_hardware_mpi(const gmx::CpuInfo& cpuInfo, } MPI_Allreduce(countsLocal.data(), countsReduced.data(), countsLocal.size(), MPI_INT, - MPI_SUM, MPI_COMM_WORLD); + MPI_SUM, physicalNodeComm.comm_); } constexpr int numElementsMax = 11; @@ -275,7 +275,7 @@ static void gmx_collect_hardware_mpi(const gmx::CpuInfo& cpuInfo, maxMinLocal[10] = (cpuIsAmdZen1 ? 1 : 0); MPI_Allreduce(maxMinLocal.data(), maxMinReduced.data(), maxMinLocal.size(), MPI_INT, - MPI_MAX, MPI_COMM_WORLD); + MPI_MAX, physicalNodeComm.comm_); } hardwareInfo->nphysicalnode = countsReduced[0]; diff --git a/src/gromacs/hardware/printhardware.cpp b/src/gromacs/hardware/printhardware.cpp index 0720a12bbe..398a03b8a2 100644 --- a/src/gromacs/hardware/printhardware.cpp +++ b/src/gromacs/hardware/printhardware.cpp @@ -213,7 +213,7 @@ static std::string detected_hardware_string(const gmx_hw_info_t* hwinfo, bool bF gmx_gethostname(host, STRLEN); - MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_rank(MPI_WORLD_STUB, &rank); // TODO Use a wrapper around MPI_Get_processor_name instead. s += gmx::formatString("Hardware detected on host %s (the node of MPI rank %d):\n", host, rank); diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index e37d921c5a..e49649281a 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -50,6 +50,8 @@ #include #include +#include "Requests.h" + #include "gromacs/applied_forces/awh/awh.h" #include "gromacs/commandline/filenm.h" #include "gromacs/domdec/collect.h" @@ -209,12 +211,11 @@ void gmx::LegacySimulator::do_mimic() int nstglobalcomm = 1; const bool bNS = true; + MimicCommunicator mimicCommunicator{}; if (MASTER(cr)) { - MimicCommunicator::init(); - auto nonConstGlobalTopology = const_cast(top_global); - MimicCommunicator::sendInitData(nonConstGlobalTopology, state_global->x); - ir->nsteps = MimicCommunicator::getStepNumber(); + mimicCommunicator.init(); + ir->nsteps=2; } ir->nstxout_compressed = 0; @@ -266,7 +267,10 @@ void gmx::LegacySimulator::do_mimic() imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr, nrnb, nullptr, FALSE); shouldCheckNumberOfBondedInteractions = true; - gmx_bcast(sizeof(ir->nsteps), &ir->nsteps, cr->mpi_comm_mygroup); + if (PAR(cr)) + { + gmx_bcast(sizeof(ir->nsteps), &ir->nsteps, cr->mpi_comm_mygroup); + } } else { @@ -361,122 +365,124 @@ void gmx::LegacySimulator::do_mimic() isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps)); while (!isLastStep) { - isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel == ir->nsteps)); - wallcycle_start(wcycle, ewcSTEP); - - t = step; - + int command = -1; if (MASTER(cr)) { - MimicCommunicator::getCoords(&state_global->x, state_global->natoms); + command = mimicCommunicator.getCommand(); } - if (ir->efep != efepNO) + if (PAR(cr)) { - state->lambda = currentLambdas(step, *(ir->fepvals), state_global->fep_state); + gmx_barrier(cr->mpi_comm_mygroup); + gmx_bcast(sizeof(int), &command, cr->mpi_comm_mygroup); } - if (MASTER(cr)) + t = step; + + auto* nonConstGlobalTopology = const_cast(top_global); + if (command == MCL_EXIT) { - const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites); - if (constructVsites && DOMAINDECOMP(cr)) + isLastStep = true; + } + else if (command == MCL_SEND_ATOM_SPECIES) + { + if (MASTER(cr)) { - gmx_fatal(FARGS, - "Vsite recalculation with -rerun is not implemented with domain " - "decomposition, " - "use a single rank"); + mimicCommunicator.sendAtomSpecies(nonConstGlobalTopology); } - if (constructVsites) + } + else if (command == MCL_SEND_NUM_FRAGMENTS) + { + if (MASTER(cr)) { - wallcycle_start(wcycle, ewcVSITECONSTR); - vsite->construct(state->x, ir->delta_t, state->v, state->box); - wallcycle_stop(wcycle, ewcVSITECONSTR); + mimicCommunicator.sendFragmentCount(nonConstGlobalTopology); } } - - if (DOMAINDECOMP(cr)) + else if (command == MCL_SEND_BOND_ATOMS) { - /* Repartition the domain decomposition */ - const bool bMasterState = true; - dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global, - *top_global, ir, imdSession, pull_work, state, &f, mdAtoms, &top, - fr, vsite, constr, nrnb, wcycle, mdrunOptions.verbose); - shouldCheckNumberOfBondedInteractions = true; + if (MASTER(cr)) + { + mimicCommunicator.sendBondAtoms(nonConstGlobalTopology); + } } - - if (MASTER(cr)) + else if (command == MCL_SEND_BOND_LENGTHS) { - EnergyOutput::printHeader(fplog, step, t); /* can we improve the information printed here? */ + if (MASTER(cr)) + { + mimicCommunicator.sendBondLengths(nonConstGlobalTopology); + } } - - if (ir->efep != efepNO) + else if (command == MCL_SEND_ATOM_ELEMENTS) { - update_mdatoms(mdatoms, state->lambda[efptMASS]); + if (MASTER(cr)) + { + mimicCommunicator.sendAtomElements(nonConstGlobalTopology); + } } - - force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES - | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved - GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); - - if (shellfc) + else if (command == MCL_SEND_NUM_ATOMS_IN_FRAGMENTS) { - /* Now is the time to relax the shells */ - relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, - imdSession, pull_work, bNS, force_flags, &top, constr, enerd, - state->natoms, state->x.arrayRefWithPadding(), - state->v.arrayRefWithPadding(), state->box, state->lambda, - &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc, - fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); + if (MASTER(cr)) + { + mimicCommunicator.sendFragmentAtomNum(nonConstGlobalTopology); + } } - else + else if (command == MCL_SEND_ATOMS_IN_FRAGMENTS) { - /* The coordinates (x) are shifted (to get whole molecules) - * in do_force. - * This is parallellized as well, and does communication too. - * Check comments in sim_util.c - */ - Awh* awh = nullptr; - gmx_edsam* ed = nullptr; - do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, - wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, - vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); + if (MASTER(cr)) + { + mimicCommunicator.sendFragmentAtomIds(nonConstGlobalTopology); + } } - - /* Now we have the energies and forces corresponding to the - * coordinates at time t. - */ + else if (command == MCL_SEND_ATOM_MASSES) { - const bool isCheckpointingStep = false; - const bool doRerun = false; - const bool bSumEkinhOld = false; - do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, - state_global, observablesHistory, top_global, fr, outf, - energyOutput, ekind, f.view().force(), isCheckpointingStep, - doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); + if (MASTER(cr)) + { + mimicCommunicator.sendAtomMasses(nonConstGlobalTopology); + } } - - stopHandler->setSignal(); - + else if (command == MCL_SEND_ATOM_MULTIPOLES) + { + if (MASTER(cr)) + { + mimicCommunicator.sendCharges(nonConstGlobalTopology); + } + } + else if (command == MCL_SEND_ATOM_COORDS) { - const bool doInterSimSignal = false; - const bool doIntraSimSignal = true; - bool bSumEkinhOld = false; - t_vcm* vcm = nullptr; - SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal); - - compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x), - makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle, - enerd, nullptr, nullptr, nullptr, nullptr, constr, &signaller, - state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld, - CGLO_GSTAT | CGLO_ENERGY - | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS - : 0)); - checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, - &top, makeConstArrayRef(state->x), state->box, - &shouldCheckNumberOfBondedInteractions); + if (MASTER(cr)) + { + mimicCommunicator.sendCoords(state_global->x); + } } - + else if (command == MCL_SEND_NUM_ATOMS) + { + if (MASTER(cr)) + { + mimicCommunicator.sendAtomNum(nonConstGlobalTopology); + } + } + else if (command == MCL_SEND_NUM_ATOM_TYPES) + { + if (MASTER(cr)) + { + mimicCommunicator.sendSpeciesNum(nonConstGlobalTopology); + } + } + else if (command == MCL_SEND_NUM_BONDS) + { + if (MASTER(cr)) + { + mimicCommunicator.sendBondsNumber(nonConstGlobalTopology); + } + } + else if (command == MCL_SEND_NUM_ANGLES) + { + if (MASTER(cr)) + { + mimicCommunicator.sendAnglesNumber(nonConstGlobalTopology); + } + } + else if (command == MCL_SEND_ATOM_FORCES) { gmx::HostVector fglobal(top_global->natoms); gmx::ArrayRef ftemp; @@ -494,73 +500,206 @@ void gmx::LegacySimulator::do_mimic() if (MASTER(cr)) { - MimicCommunicator::sendEnergies(enerd->term[F_EPOT]); - MimicCommunicator::sendForces(ftemp, state_global->natoms); + mimicCommunicator.sendForces(ftemp, state_global->natoms); } } - - - /* Note: this is OK, but there are some numerical precision issues with using the convergence of - the virial that should probably be addressed eventually. state->veta has better properies, - but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could - generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ - - if (ir->efep != efepNO) + else if (command == MCL_SEND_ENERGY) { - /* Sum up the foreign energy and dhdl terms for md and sd. - Currently done every step so that dhdl is correct in the .edr */ - accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals); + if (MASTER(cr)) + { + mimicCommunicator.sendEnergies(enerd->term[F_EPOT]); + } } - - /* Output stuff */ - if (MASTER(cr)) + else if (command == MCL_SEND_MULTIPOLE_ORDER) { - const bool bCalcEnerStep = true; - energyOutput.addDataAtEnergyStep( - doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals, - ir->expandedvals, state->box, - PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi, - state->nhpres_xi, state->nhpres_vxi }), - state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr); - - const bool do_ene = true; - const bool do_log = true; - Awh* awh = nullptr; - const bool do_dr = ir->nstdisreout != 0; - const bool do_or = ir->nstorireout != 0; - - EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts)); - energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, - do_log ? fplog : nullptr, step, t, fr->fcdata.get(), awh); - - if (do_per_step(step, ir->nstlog)) - { - if (fflush(fplog) != 0) + if (MASTER(cr)) + { + mimicCommunicator.sendMultipoleOrder(nonConstGlobalTopology); + } + } + else if (command == MCL_RECV_ATOM_COORDS) + { + wallcycle_start(wcycle, ewcSTEP); + if (MASTER(cr)) + { + mimicCommunicator.getCoords(&state_global->x, state_global->natoms); + } + + if (ir->efep != efepNO) + { + state->lambda = currentLambdas(step, *(ir->fepvals), state_global->fep_state); + } + + if (MASTER(cr)) + { + const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites); + if (constructVsites && DOMAINDECOMP(cr)) { - gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + gmx_fatal(FARGS, + "Vsite recalculation with -rerun is not implemented with domain " + "decomposition, " + "use a single rank"); + } + if (constructVsites) + { + wallcycle_start(wcycle, ewcVSITECONSTR); + vsite->construct(state->x, ir->delta_t, state->v, state->box); + wallcycle_stop(wcycle, ewcVSITECONSTR); } } - } - /* Print the remaining wall clock time for the run */ - if (isMasterSimMasterRank(ms, MASTER(cr)) && (mdrunOptions.verbose || gmx_got_usr_signal())) + if (DOMAINDECOMP(cr)) + { + /* Repartition the domain decomposition */ + const bool bMasterState = true; + dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global, + *top_global, ir, imdSession, pull_work, state, &f, mdAtoms, + &top, fr, vsite, constr, nrnb, wcycle, mdrunOptions.verbose); + shouldCheckNumberOfBondedInteractions = true; + } + + if (MASTER(cr)) + { + EnergyOutput::printHeader(fplog, step, t); /* can we improve the information printed here? */ + } + + if (ir->efep != efepNO) + { + update_mdatoms(mdatoms, state->lambda[efptMASS]); + } + } + else if (command == MCL_COMPUTE_FORCE_ENERGY) { + force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES + | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved + GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0)); + if (shellfc) { - fprintf(stderr, "\n"); + /* Now is the time to relax the shells */ + relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir, + imdSession, pull_work, bNS, force_flags, &top, constr, enerd, + state->natoms, state->x.arrayRefWithPadding(), + state->v.arrayRefWithPadding(), state->box, state->lambda, + &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc, + fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); + } + else + { + /* The coordinates (x) are shifted (to get whole molecules) + * in do_force. + * This is parallellized as well, and does communication too. + * Check comments in sim_util.c + */ + Awh* awh = nullptr; + gmx_edsam* ed = nullptr; + do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, + wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, + &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork, + vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); + } + + /* Now we have the energies and forces corresponding to the + * coordinates at time t. + */ + { + const bool isCheckpointingStep = false; + const bool doRerun = false; + const bool bSumEkinhOld = false; + do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, + state_global, observablesHistory, top_global, fr, outf, + energyOutput, ekind, f.view().force(), isCheckpointingStep, + doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld); + } + + stopHandler->setSignal(); + + { + const bool doInterSimSignal = false; + const bool doIntraSimSignal = true; + bool bSumEkinhOld = false; + t_vcm* vcm = nullptr; + SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal); + + compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x), + makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle, + enerd, nullptr, nullptr, nullptr, nullptr, constr, &signaller, + state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld, + CGLO_GSTAT | CGLO_ENERGY + | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS + : 0)); + checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, + top_global, &top, makeConstArrayRef(state->x), + state->box, &shouldCheckNumberOfBondedInteractions); + } + + /* Note: this is OK, but there are some numerical precision issues with using the convergence of + the virial that should probably be addressed eventually. state->veta has better properies, + but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could + generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ + + if (ir->efep != efepNO) + { + /* Sum up the foreign energy and dhdl terms for md and sd. + Currently done every step so that dhdl is correct in the .edr */ + accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals); } - print_time(stderr, walltime_accounting, step, ir, cr); - } - cycles = wallcycle_stop(wcycle, ewcSTEP); - if (DOMAINDECOMP(cr) && wcycle) + /* Output stuff */ + if (MASTER(cr)) + { + const bool bCalcEnerStep = true; + energyOutput.addDataAtEnergyStep( + doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd, + ir->fepvals, ir->expandedvals, state->box, + PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi, + state->nhpres_xi, state->nhpres_vxi }), + state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr); + + const bool do_ene = true; + const bool do_log = true; + Awh* awh = nullptr; + const bool do_dr = ir->nstdisreout != 0; + const bool do_or = ir->nstorireout != 0; + + EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts)); + energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, + do_log ? fplog : nullptr, step, t, + fr->fcdata.get(), awh); + + if (do_per_step(step, ir->nstlog)) + { + if (fflush(fplog) != 0) + { + gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?"); + } + } + } + + /* Print the remaining wall clock time for the run */ + if (isMasterSimMasterRank(ms, MASTER(cr)) && (mdrunOptions.verbose || gmx_got_usr_signal())) + { + if (shellfc) + { + fprintf(stderr, "\n"); + } + print_time(stderr, walltime_accounting, step, ir, cr); + } + + cycles = wallcycle_stop(wcycle, ewcSTEP); + if (DOMAINDECOMP(cr) && wcycle) + { + dd_cycles_add(cr->dd, cycles, ddCyclStep); + } + + /* increase the MD step number */ + step++; + step_rel++; + } + else { - dd_cycles_add(cr->dd, cycles, ddCyclStep); + gmx_fatal(FARGS, "Unrecognized MiMiC command: %d", command); } - - /* increase the MD step number */ - step++; - step_rel++; } /* End of main MD loop */ @@ -573,7 +712,7 @@ void gmx::LegacySimulator::do_mimic() if (MASTER(cr)) { - MimicCommunicator::finalize(); + mimicCommunicator.finalize(); } if (!thisRankHasDuty(cr, DUTY_PME)) diff --git a/src/gromacs/mdrunutility/threadaffinity.cpp b/src/gromacs/mdrunutility/threadaffinity.cpp index 8cf377acbf..f299888cd8 100644 --- a/src/gromacs/mdrunutility/threadaffinity.cpp +++ b/src/gromacs/mdrunutility/threadaffinity.cpp @@ -542,7 +542,7 @@ static bool detectDefaultAffinityMask(const int nthreads_hw_avail) if (mpiIsInitialized) { bool maskToReduce = detectedDefaultAffinityMask; - MPI_Allreduce(&maskToReduce, &detectedDefaultAffinityMask, 1, MPI_C_BOOL, MPI_LAND, MPI_COMM_WORLD); + MPI_Allreduce(&maskToReduce, &detectedDefaultAffinityMask, 1, MPI_C_BOOL, MPI_LAND, MPI_WORLD_STUB); } #endif diff --git a/src/gromacs/mimic/communicator.cpp b/src/gromacs/mimic/communicator.cpp index 12dcf874d3..28a6087fcf 100644 --- a/src/gromacs/mimic/communicator.cpp +++ b/src/gromacs/mimic/communicator.cpp @@ -46,6 +46,7 @@ #if GMX_MIMIC # include # include +# include #endif // When not built in a configuration with QMMM support, much of this @@ -102,26 +103,61 @@ void gmx::MimicCommunicator::init() { char path[GMX_PATH_MAX]; gmx_getcwd(path, GMX_PATH_MAX); - MCL_init_client(path); + char del = ';'; + MCL_handshake(path, &del, 0); } -void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector coords) +void gmx::MimicCommunicator::sendAtomSpecies(gmx_mtop_t* mtop) { + std::vector atomTypes; + atomTypes.reserve(static_cast(mtop->natoms)); + + int offset = 0; + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + for (int at = 0; at < type->atoms.nr; ++at) + { + int atomtype = type->atoms.atom[at].type; + atomTypes.push_back(atomtype + 1); + } + } + } + MCL_send(&*atomTypes.begin(), mtop->natoms, TYPE_INT, MCL_DATA, 0); +} +void gmx::MimicCommunicator::sendFragmentCount(gmx_mtop_t* mtop) { - MCL_send(&mtop->natoms, 1, TYPE_INT, 0); - MCL_send(&mtop->atomtypes.nr, 1, TYPE_INT, 0); + int nMolecules = 0; + for (const gmx_molblock_t& molblock : mtop->molblock) + { + nMolecules += molblock.nmol; + } + MCL_send(&nMolecules, 1, TYPE_INT, MCL_DATA, 0); +} - std::vector atomTypes; - std::vector nAtomsMol; - std::vector idOrder; - std::vector charges; - std::vector masses(mtop->atomtypes.nr, -1); - std::vector elements(mtop->atomtypes.nr, -1); - std::vector bonds; - std::vector bondLengths; - std::unordered_set existingTypes; +void gmx::MimicCommunicator::sendBondsNumber(gmx_mtop_t *mtop) +{ + int bondsNumber = 0; + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + int nconstr = type->ilist[F_CONSTR].size() / 3; + int nconstrc = type->ilist[F_CONSTRNC].size() / 3; + int nsettle = type->ilist[F_SETTLE].size() / 4; - atomTypes.reserve(static_cast(mtop->natoms)); - charges.reserve(static_cast(mtop->natoms)); + bondsNumber += nconstr + nconstrc + nsettle * 3; + } + } + + MCL_send(&bondsNumber, 1, TYPE_INT, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendBondAtoms(gmx_mtop_t *mtop) +{ + std::vector bonds; int offset = 0; for (const gmx_molblock_t& molblock : mtop->molblock) @@ -135,13 +171,10 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVectorilist[F_CONSTR].iatoms[0]; int at1 = type->ilist[F_CONSTR].iatoms[1]; int at2 = type->ilist[F_CONSTR].iatoms[2]; bonds.push_back(offset + at1 + 1); bonds.push_back(offset + at2 + 1); - bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dA) - / BOHR2NM); } for (int ncon = 0; ncon < nsettle; ++ncon) @@ -150,8 +183,6 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVectorilist[F_SETTLE].iatoms[0]; - ox = type->ilist[F_SETTLE].iatoms[1]; h1 = type->ilist[F_SETTLE].iatoms[2]; h2 = type->ilist[F_SETTLE].iatoms[3]; @@ -164,6 +195,38 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVectoratoms.nr; + } + } + + MCL_send(&*bonds.begin(), bonds.size(), TYPE_INT, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendBondLengths(gmx_mtop_t* mtop) +{ + std::vector bondLengths; + + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + int nconstr = type->ilist[F_CONSTR].size() / 3; + int nconstrc = type->ilist[F_CONSTRNC].size() / 3; + int nsettle = type->ilist[F_SETTLE].size() / 4; + + for (int ncon = 0; ncon < nconstr + nconstrc; ++ncon) + { + int contype = type->ilist[F_CONSTR].iatoms[0]; + bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dA) + / BOHR2NM); + } + + for (int ncon = 0; ncon < nsettle; ++ncon) + { + int contype = type->ilist[F_SETTLE].iatoms[0]; bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dA) / BOHR2NM); bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dA) @@ -172,65 +235,128 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVectoratoms.nr); + } + } + MCL_send(&*bondLengths.begin(), bondLengths.size(), TYPE_DOUBLE, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendAnglesNumber(gmx_mtop_t* mtop) +{ + int dummyAngles = 0; + MCL_send(&dummyAngles, 1, TYPE_INT, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendAtomElements(gmx_mtop_t* mtop) +{ + std::vector elements(mtop->atomtypes.nr, -1); + std::unordered_set existingTypes; + + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { for (int at = 0; at < type->atoms.nr; ++at) { int atomtype = type->atoms.atom[at].type; - auto charge = static_cast(type->atoms.atom[at].q); - idOrder.push_back(offset + 1); - offset++; - atomTypes.push_back(atomtype + 1); - charges.push_back(charge); if (existingTypes.insert(atomtype).second) { - masses[atomtype] = type->atoms.atom[at].m; elements[atomtype] = type->atoms.atom[at].atomnumber; } } } } - // sending atom types - MCL_send(&*atomTypes.begin(), mtop->natoms, TYPE_INT, 0); - int max_multipole_order = 0; - // sending multipole orders - MCL_send(&max_multipole_order, 1, TYPE_INT, 0); - - int nMolecules = nAtomsMol.size(); - // sending molecule number - MCL_send(&nMolecules, 1, TYPE_INT, 0); + MCL_send(&*elements.begin(), mtop->atomtypes.nr, TYPE_INT, MCL_DATA, 0); +} - // sending number of atoms per molecules - MCL_send(&*nAtomsMol.begin(), nAtomsMol.size(), TYPE_INT, 0); +void gmx::MimicCommunicator::sendAtomMasses(gmx_mtop_t* mtop) +{ + std::vector masses(mtop->atomtypes.nr, -1); + std::unordered_set existingTypes; - int nBonds = bonds.size() / 2; - // sending number of bond constraints - MCL_send(&nBonds, 1, TYPE_INT, 0); + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + for (int at = 0; at < type->atoms.nr; ++at) + { + int atomtype = type->atoms.atom[at].type; + if (existingTypes.insert(atomtype).second) + { + masses[atomtype] = type->atoms.atom[at].m; + } + } + } + } - // sending number of angle constraints - MCL_send(&max_multipole_order, 1, TYPE_INT, 0); + MCL_send(&*masses.begin(), masses.size(), TYPE_DOUBLE, MCL_DATA, 0); +} - if (nBonds > 0) +void gmx::MimicCommunicator::sendFragmentAtomNum(gmx_mtop_t* mtop) +{ + std::vector nAtomsMol; + for (const gmx_molblock_t& molblock : mtop->molblock) { - // sending bonded atoms indices - MCL_send(&*bonds.begin(), bonds.size(), TYPE_INT, 0); - - // sending bond lengths - MCL_send(&*bondLengths.begin(), bondLengths.size(), TYPE_DOUBLE, 0); + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + nAtomsMol.push_back(type->atoms.nr); + } } + // sending number of atoms per molecules + MCL_send(&*nAtomsMol.begin(), nAtomsMol.size(), TYPE_INT, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendFragmentAtomIds(gmx_mtop_t* mtop) +{ + std::vector idOrder; - // sending array of atomic charges - MCL_send(&*charges.begin(), mtop->natoms, TYPE_DOUBLE, 0); + int offset = 0; + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + for (int at = 0; at < type->atoms.nr; ++at) + { + idOrder.push_back(offset + 1); + offset++; + } + } + } + MCL_send(&*idOrder.begin(), idOrder.size(), TYPE_INT, MCL_DATA, 0); +} - // sending array of atomic masses - MCL_send(&*masses.begin(), mtop->atomtypes.nr, TYPE_DOUBLE, 0); +void gmx::MimicCommunicator::sendMultipoleOrder(gmx_mtop_t* mtop) +{ + int multipleOrder = 0; + MCL_send(&multipleOrder, 1, TYPE_INT, MCL_DATA, 0); +} - // sending ids of atoms per molecule - MCL_send(&*idOrder.begin(), idOrder.size(), TYPE_INT, 0); +void gmx::MimicCommunicator::sendCharges(gmx_mtop_t* mtop) +{ + std::vector charges; + charges.reserve(static_cast(mtop->natoms)); - // sending list of elements - MCL_send(&*elements.begin(), mtop->atomtypes.nr, TYPE_INT, 0); + for (const gmx_molblock_t& molblock : mtop->molblock) + { + gmx_moltype_t* type = &mtop->moltype[molblock.type]; + for (int mol = 0; mol < molblock.nmol; ++mol) + { + for (int at = 0; at < type->atoms.nr; ++at) + { + auto charge = static_cast(type->atoms.atom[at].q); + charges.push_back(charge); + } + } + } + MCL_send(&*charges.begin(), mtop->natoms, TYPE_DOUBLE, MCL_DATA, 0); +} +void gmx::MimicCommunicator::sendCoords(PaddedHostVector coords) +{ std::vector convertedCoords; for (auto& coord : coords) { @@ -240,20 +366,27 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVectornatoms, TYPE_DOUBLE, 0); + MCL_send(&*convertedCoords.begin(), convertedCoords.size(), TYPE_DOUBLE, MCL_DATA, 0); } -int64_t gmx::MimicCommunicator::getStepNumber() -{ - int steps; - MCL_receive(&steps, 1, TYPE_INT, 0); - return steps; +void gmx::MimicCommunicator::sendAtomNum(gmx_mtop_t* mtop) { + MCL_send(&mtop->natoms, 1, TYPE_INT, MCL_DATA, 0); +} + +void gmx::MimicCommunicator::sendSpeciesNum(gmx_mtop_t* mtop) { + MCL_send(&mtop->atomtypes.nr, 1, TYPE_INT, MCL_DATA, 0); +} + +int gmx::MimicCommunicator::getCommand() { + int command{}; + MCL_receive(&command, 1, TYPE_INT, MCL_COMMAND, 0); + return command; } void gmx::MimicCommunicator::getCoords(PaddedHostVector* x, const int natoms) { std::vector coords(natoms * 3); - MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, 0); + MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, MCL_DATA, 0); for (int j = 0; j < natoms; ++j) { (*x)[j][0] = static_cast(coords[j * 3] * BOHR2NM); @@ -265,7 +398,7 @@ void gmx::MimicCommunicator::getCoords(PaddedHostVector* x, const int nato void gmx::MimicCommunicator::sendEnergies(real energy) { double convertedEnergy = energy / (HARTREE2KJ * AVOGADRO); - MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, 0); + MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, MCL_DATA, 0); } void gmx::MimicCommunicator::sendForces(gmx::ArrayRef forces, int natoms) @@ -277,7 +410,7 @@ void gmx::MimicCommunicator::sendForces(gmx::ArrayRef forces, int nat convertedForce.push_back(static_cast(forces[j][1]) / HARTREE_BOHR2MD); convertedForce.push_back(static_cast(forces[j][2]) / HARTREE_BOHR2MD); } - MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, 0); + MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, MCL_DATA, 0); } void gmx::MimicCommunicator::finalize() diff --git a/src/gromacs/mimic/communicator.h b/src/gromacs/mimic/communicator.h index 6d2b9c2761..436a13c784 100644 --- a/src/gromacs/mimic/communicator.h +++ b/src/gromacs/mimic/communicator.h @@ -60,27 +60,7 @@ public: /*! \brief * Initializes the communicator */ - static void init(); - - /*! \brief - * Sends the data needed for MiMiC initialization - * - * That includes number of atoms, element numbers, charges, masses, - * maximal order of multipoles (0 for point-charge forcefields), - * number of molecules, number of atoms per each molecule, - * bond constraints data - * - * @param mtop global topology data - * @param coords coordinates of all atoms - */ - static void sendInitData(gmx_mtop_t* mtop, PaddedHostVector coords); - - /*! \brief - * Gets the number of MD steps to perform from MiMiC - * - * @return nsteps the number of MD steps to perform - */ - static int64_t getStepNumber(); + void init(); /*! \brief * Receive and array of updated atomic coordinates from MiMiC @@ -88,14 +68,14 @@ public: * @param x array of coordinates to fill * @param natoms number of atoms in the system */ - static void getCoords(PaddedHostVector* x, int natoms); + void getCoords(PaddedHostVector* x, int natoms); /*! \brief * Send the potential energy value to MiMiC * * @param energy energy value to send */ - static void sendEnergies(real energy); + void sendEnergies(real energy); /*! \brief * Send classical forces acting on all atoms in the system @@ -104,12 +84,44 @@ public: * @param forces array of forces to send * @param natoms number of atoms in the system */ - static void sendForces(ArrayRef forces, int natoms); + void sendForces(ArrayRef forces, int natoms); /*! \brief * Finish communications and disconnect from the server */ - static void finalize(); + void finalize(); + + void sendAtomSpecies(gmx_mtop_t *mtop); + + void sendFragmentCount(gmx_mtop_t *mtop); + + void sendBondsNumber(gmx_mtop_t *mtop); + + void sendBondAtoms(gmx_mtop_t *mtop); + + void sendBondLengths(gmx_mtop_t *mtop); + + void sendAnglesNumber(gmx_mtop_t *mtop); + + void sendAtomElements(gmx_mtop_t *mtop); + + void sendAtomMasses(gmx_mtop_t *mtop); + + void sendFragmentAtomNum(gmx_mtop_t *mtop); + + void sendFragmentAtomIds(gmx_mtop_t *mtop); + + void sendMultipoleOrder(gmx_mtop_t *mtop); + + void sendCharges(gmx_mtop_t *mtop); + + void sendCoords(PaddedHostVector coords); + + int getCommand(); + + void sendAtomNum(gmx_mtop_t *mtop); + + void sendSpeciesNum(gmx_mtop_t *mtop); }; } // namespace gmx diff --git a/src/gromacs/tools/pme_error.cpp b/src/gromacs/tools/pme_error.cpp index ba1fb3b142..7f2c31fecc 100644 --- a/src/gromacs/tools/pme_error.cpp +++ b/src/gromacs/tools/pme_error.cpp @@ -1121,7 +1121,7 @@ int gmx_pme_error(int argc, char* argv[]) #define NFILE asize(fnm) - CommrecHandle commrecHandle = init_commrec(MPI_COMM_WORLD); + CommrecHandle commrecHandle = init_commrec(MPI_WORLD_STUB); t_commrec* cr = commrecHandle.get(); PCA_Flags = PCA_NOEXIT_ON_ARGS; diff --git a/src/gromacs/utility/basenetwork.cpp b/src/gromacs/utility/basenetwork.cpp index 6fade03e57..bae56d5259 100644 --- a/src/gromacs/utility/basenetwork.cpp +++ b/src/gromacs/utility/basenetwork.cpp @@ -74,7 +74,7 @@ int gmx_node_num() } # endif int i; - (void)MPI_Comm_size(MPI_COMM_WORLD, &i); + (void)MPI_Comm_size(MPI_WORLD_STUB, &i); return i; #endif } @@ -91,7 +91,7 @@ int gmx_node_rank() } # endif int i; - (void)MPI_Comm_rank(MPI_COMM_WORLD, &i); + (void)MPI_Comm_rank(MPI_WORLD_STUB, &i); return i; #endif } @@ -153,7 +153,7 @@ int gmx_physicalnode_id_hash() void gmx_broadcast_world(int size, void* buffer) { #if GMX_MPI - MPI_Bcast(buffer, size, MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(buffer, size, MPI_BYTE, 0, MPI_WORLD_STUB); #else GMX_UNUSED_VALUE(size); GMX_UNUSED_VALUE(buffer); @@ -163,7 +163,7 @@ void gmx_broadcast_world(int size, void* buffer) #if GMX_LIB_MPI void gmx_abort(int errorno) { - MPI_Abort(MPI_COMM_WORLD, errorno); + MPI_Abort(MPI_WORLD_STUB, errorno); std::abort(); } #endif diff --git a/src/gromacs/utility/fatalerror.cpp b/src/gromacs/utility/fatalerror.cpp index f29b400a80..d09225c000 100644 --- a/src/gromacs/utility/fatalerror.cpp +++ b/src/gromacs/utility/fatalerror.cpp @@ -186,7 +186,7 @@ void gmx_exit_on_fatal_error(ExitType exitType, int returnValue) case ExitType_NonMasterAbort: // Let all other processes wait till the master has printed // the error message and issued MPI_Abort. - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(MPI_WORLD_STUB); break; } } diff --git a/src/gromacs/utility/gmxmpi.h b/src/gromacs/utility/gmxmpi.h index 5c2a247526..f226e3644c 100644 --- a/src/gromacs/utility/gmxmpi.h +++ b/src/gromacs/utility/gmxmpi.h @@ -60,6 +60,7 @@ /* disable bindings for SGI MPT also */ # define MPI_NO_CPPBIND 1 # include +extern MPI_Comm MPI_WORLD_STUB; /* Starting with 2.2 MPI_INT64_T is required. Earlier version still might have it. In theory MPI_Datatype doesn't have to be a #define, but current available MPI implementations (OpenMPI + MPICH (+derivates)) use #define and future versions diff --git a/src/gromacs/utility/init.cpp b/src/gromacs/utility/init.cpp index 1f6e07414b..f2ca2a3986 100644 --- a/src/gromacs/utility/init.cpp +++ b/src/gromacs/utility/init.cpp @@ -53,6 +53,10 @@ # include "gromacs/utility/gmxmpi.h" #endif +#include "MessageApi.h" + +MPI_Comm MPI_WORLD_STUB = MPI_COMM_WORLD; + namespace gmx { @@ -114,6 +118,7 @@ void init(int* argc, char*** argv) // NOLINT(readability-non-const-parameter) # else MPI_Init(argc, argv); # endif + MCL_init(&MPI_WORLD_STUB); } // Bump the counter to record this initialization event g_initializationCounter++; @@ -137,7 +142,7 @@ void finalize() * with buggy MPI implementations that could cause * unfinished processes to terminate. */ - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(MPI_WORLD_STUB); /* Apparently certain mpich implementations cause problems * with MPI_Finalize. In that case comment out MPI_Finalize. diff --git a/src/programs/mdrun/mdrun.cpp b/src/programs/mdrun/mdrun.cpp index 77a25f1fdf..73ded62fd4 100644 --- a/src/programs/mdrun/mdrun.cpp +++ b/src/programs/mdrun/mdrun.cpp @@ -70,6 +70,10 @@ #include "gromacs/utility/basenetwork.h" #include "gromacs/utility/physicalnodecommunicator.h" +#if GMX_MIMIC +# include +#endif + #include "mdrun_main.h" namespace gmx @@ -79,7 +83,7 @@ int gmx_mdrun(int argc, char* argv[]) { // Set up the communicator, where possible (see docs for // SimulationContext). - MPI_Comm communicator = GMX_LIB_MPI ? MPI_COMM_WORLD : MPI_COMM_NULL; + MPI_Comm communicator = GMX_LIB_MPI ? MPI_WORLD_STUB : MPI_COMM_NULL; PhysicalNodeCommunicator physicalNodeCommunicator(communicator, gmx_physicalnode_id_hash()); std::unique_ptr hwinfo = gmx_detect_hardware(physicalNodeCommunicator); return gmx_mdrun(communicator, *hwinfo, argc, argv);