From bec71a22a932b901e9afd97362a54a0d8b1660a9 Mon Sep 17 00:00:00 2001 From: "Fattebert J.-L" Date: Thu, 10 Oct 2024 11:04:43 -0400 Subject: [PATCH] Clean up class Ions, add test for it --- src/Control.cc | 1 + src/Ions.cc | 324 ++++++++++++++++----------------------- src/Ions.h | 23 ++- src/KBprojectorSparse.cc | 4 +- tests/CMakeLists.txt | 8 + tests/testIons.cc | 104 +++++++++++++ 6 files changed, 265 insertions(+), 199 deletions(-) create mode 100644 tests/testIons.cc diff --git a/src/Control.cc b/src/Control.cc index af3d8ed9..dd786678 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -70,6 +70,7 @@ Control::Control() dm_approx_ndigits = 1; dm_approx_power_maxits = 100; wf_extrapolation_ = 1; + verbose = 0; // undefined values dm_algo_ = -1; diff --git a/src/Ions.cc b/src/Ions.cc index 940e12e8..852151c7 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -584,11 +584,9 @@ void Ions::printPositionsLocal(std::ostream& os, const int root) const os.setf(std::ios::right, std::ios::adjustfield); os.setf(std::ios::fixed, std::ios::floatfield); - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - (*ion)->printPosition(os); - ion++; + ion->printPosition(os); } os << std::endl; @@ -627,14 +625,12 @@ void Ions::writeAtomicNumbers(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - assert((*ion)->atomic_number() > 0); - assert((*ion)->atomic_number() < 200); + assert(ion->atomic_number() > 0); + assert(ion->atomic_number() < 200); - data.push_back((*ion)->atomic_number()); - ion++; + data.push_back(ion->atomic_number()); } } @@ -685,17 +681,15 @@ void Ions::writeAtomNames(HDFrestart& h5f_file) void Ions::lockAtom(const std::string& name) { - std::vector::iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - std::string name_ion((*ion)->name()); + std::string name_ion(ion->name()); if (name.compare(name_ion) == 0) { - (*ion)->lock(); + ion->lock(); if (onpe0) (*MPIdata::sout) << "Lock atom " << name << std::endl; break; } - ion++; } } @@ -975,12 +969,10 @@ void Ions::readRestartPositions(HDFrestart& h5_file) std::vector data; h5_file.readAtomicPositions(data); - int i = 0; - std::vector::iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + int i = 0; + for (auto& ion : local_ions_) { - (*ion)->setPosition(data[3 * i], data[3 * i + 1], data[3 * i + 2]); - ion++; + ion->setPosition(data[3 * i], data[3 * i + 1], data[3 * i + 2]); i++; } } @@ -1109,14 +1101,12 @@ void Ions::removeMassCenterMotion() #endif } - ion = local_ions_.begin(); - i = 0; - while (ion != local_ions_.end()) + i = 0; + for (auto& ion : local_ions_) { const int threei = 3 * i; - (*ion)->setVelocity(velocities[threei] - tmp[0], + ion->setVelocity(velocities[threei] - tmp[0], velocities[threei + 1] - tmp[1], velocities[threei + 2] - tmp[2]); - ion++; i++; } @@ -1124,16 +1114,15 @@ void Ions::removeMassCenterMotion() mv[0] = 0.; mv[1] = 0.; mv[2] = 0.; - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - velocities[3 * i] = (*ion)->velocity(0); - velocities[3 * i + 1] = (*ion)->velocity(1); - velocities[3 * i + 2] = (*ion)->velocity(2); + velocities[3 * i] = ion->velocity(0); + velocities[3 * i + 1] = ion->velocity(1); + velocities[3 * i + 2] = ion->velocity(2); mv[0] += mass[i] * velocities[3 * i]; mv[1] += mass[i] * velocities[3 * i + 1]; mv[2] += mass[i] * velocities[3 * i + 2]; - ion++; i++; } @@ -1155,12 +1144,10 @@ void Ions::readRestartVelocities(HDFrestart& h5_file) std::vector data; h5_file.readAtomicVelocities(data); - std::vector::iterator ion = local_ions_.begin(); - int i = 0; - while (ion != local_ions_.end()) + int i = 0; + for (auto& ion : local_ions_) { - (*ion)->setVelocity(data[3 * i], data[3 * i + 1], data[3 * i + 2]); - ion++; + ion->setVelocity(data[3 * i], data[3 * i + 1], data[3 * i + 2]); i++; } } @@ -1174,12 +1161,10 @@ void Ions::readRestartRandomStates(HDFrestart& h5f_file) std::vector data; h5f_file.readRestartRandomStates(data); - std::vector::iterator ion = local_ions_.begin(); - int i = 0; - while (ion != local_ions_.end()) + int i = 0; + for (auto& ion : local_ions_) { - (*ion)->setRandomState(data[3 * i], data[3 * i + 1], data[3 * i + 2]); - ion++; + ion->setRandomState(data[3 * i], data[3 * i + 1], data[3 * i + 2]); i++; } } @@ -1202,17 +1187,14 @@ void Ions::writeForces(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { // get position of local ion double force[3]; - (*ion)->getForce(&force[0]); + ion->getForce(&force[0]); data.push_back(force[0]); data.push_back(force[1]); data.push_back(force[2]); - - ++ion; } } @@ -1364,64 +1346,54 @@ void Ions::printForcesLocal(std::ostream& os, const int root) const << "FX" << std::setw(10) << "FY" << std::setw(10) << "FZ" << std::endl; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { + ion->printPositionAndForce(os); - (*ion)->printPositionAndForce(os); - - if (!(*ion)->locked()) + if (!ion->locked()) { + avg_forces[0] += fabs(ion->force(0)); + avg_forces[1] += fabs(ion->force(1)); + avg_forces[2] += fabs(ion->force(2)); - avg_forces[0] += fabs((*ion)->force(0)); - avg_forces[1] += fabs((*ion)->force(1)); - avg_forces[2] += fabs((*ion)->force(2)); - - double ff = (*ion)->norm2F(); + double ff = ion->norm2F(); maxf[0] = std::max(maxf[0], ff); - max_forces[0] = std::max(max_forces[0], fabs((*ion)->force(0))); - max_forces[1] = std::max(max_forces[1], fabs((*ion)->force(1))); - max_forces[2] = std::max(max_forces[2], fabs((*ion)->force(2))); + max_forces[0] = std::max(max_forces[0], fabs(ion->force(0))); + max_forces[1] = std::max(max_forces[1], fabs(ion->force(1))); + max_forces[2] = std::max(max_forces[2], fabs(ion->force(2))); num_movable++; } for (short ii = 0; ii < 3; ii++) - sum_forces[ii] += (*ion)->force(ii); + sum_forces[ii] += ion->force(ii); num_atoms++; - - ion++; } } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - - if (!(*ion)->locked()) + if (!ion->locked()) { + avg_forces[0] += fabs(ion->force(0)); + avg_forces[1] += fabs(ion->force(1)); + avg_forces[2] += fabs(ion->force(2)); - avg_forces[0] += fabs((*ion)->force(0)); - avg_forces[1] += fabs((*ion)->force(1)); - avg_forces[2] += fabs((*ion)->force(2)); - - double ff = (*ion)->norm2F(); + double ff = ion->norm2F(); maxf[0] = std::max(maxf[0], ff); - max_forces[0] = std::max(max_forces[0], fabs((*ion)->force(0))); - max_forces[1] = std::max(max_forces[1], fabs((*ion)->force(1))); - max_forces[2] = std::max(max_forces[2], fabs((*ion)->force(2))); + max_forces[0] = std::max(max_forces[0], fabs(ion->force(0))); + max_forces[1] = std::max(max_forces[1], fabs(ion->force(1))); + max_forces[2] = std::max(max_forces[2], fabs(ion->force(2))); num_movable++; } for (short ii = 0; ii < 3; ii++) - sum_forces[ii] += (*ion)->force(ii); + sum_forces[ii] += ion->force(ii); num_atoms++; - - ion++; } } // global statistics @@ -1479,12 +1451,10 @@ int Ions::countIonsHere() const { return (int)local_ions_.size(); } int Ions::countProjectorsHere() const { - int count = 0; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + int count = 0; + for (auto& ion : local_ions_) { - count += (*ion)->nProjectors(); - ion++; + count += ion->nProjectors(); } return count; } @@ -1497,12 +1467,10 @@ int Ions::countProjectors() const Mesh* mymesh = Mesh::instance(); const pb::PEenv& myPEenv = mymesh->peenv(); - int nproj = 0; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + int nproj = 0; + for (auto& ion : local_ions_) { - nproj += (*ion)->nProjectors(); - ion++; + nproj += ion->nProjectors(); } int tmp = nproj; MPI_Allreduce(&tmp, &nproj, 1, MPI_INT, MPI_SUM, myPEenv.comm()); @@ -1569,13 +1537,10 @@ void Ions::setLocalPositions(const std::vector& tau) { assert(tau.size() == 3 * local_ions_.size()); - std::vector::iterator ion = local_ions_.begin(); - int ia = 0; - while (ion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { - (*ion)->setPosition(tau[3 * ia + 0], tau[3 * ia + 1], tau[3 * ia + 2]); - - ion++; + ion->setPosition(tau[3 * ia + 0], tau[3 * ia + 1], tau[3 * ia + 2]); ia++; } @@ -1780,45 +1745,8 @@ int Ions::setAtoms( if (ia < 1000) aname.append("0"); if (ia < 10000) aname.append("0"); aname.append(ss.str()); - Ion* new_ion - = new Ion(species_[isp], aname, &crds[3 * ia], velocity, locked); - new_ion->bcast(mmpi.commGlobal()); - - // Populate list_ions_ list - // std::cout<<"crds: "< 2) - (*MPIdata::sout) - << "Ion " << aname << " at position " << crds[3 * ia + 0] - << "," << crds[3 * ia + 1] << "," << crds[3 * ia + 2] - << " added to the list... on PE" << mmpi.mypeGlobal() - << std::endl; - // populate local_ions_ list - if (inLocalIons( - crds[3 * ia + 0], crds[3 * ia + 1], crds[3 * ia + 2])) - { - (new_ion)->set_here(true); - local_ions_.push_back(new_ion); - if (onpe0 && ct.verbose > 2) - (*MPIdata::sout) - << "Ion " << aname << " at position " - << crds[3 * ia + 0] << "," << crds[3 * ia + 1] << "," - << crds[3 * ia + 2] - << " added to the list of local ions... on PE" - << mmpi.mypeGlobal() << std::endl; - } - else - (new_ion)->set_here(false); - } - else - { - //(*MPIdata::sout)<<"Ion "<set_here(true); + local_ions_.push_back(new_ion); + if (onpe0 && ct.verbose > 2) + (*MPIdata::sout) << "Ion " << name << " at position " << crds[0] + << "," << crds[1] << "," << crds[2] + << " added to the list of local ions... on PE" + << mmpi.mypeGlobal() << std::endl; + } + else + (new_ion)->set_here(false); + } + else + { + // delete Ion if not put in list + delete new_ion; + } +} + int Ions::readNatoms(const std::string& filename, const bool cell_relative) { Control& ct(*(Control::instance())); @@ -2117,9 +2085,8 @@ void Ions::setVelocities(const std::vector& tau0, assert(tau0.size() == 3 * local_ions_.size()); assert(taup.size() == 3 * local_ions_.size()); - int ia = 0; - std::vector::iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { double v[3]; for (short i = 0; i < 3; i++) @@ -2128,8 +2095,7 @@ void Ions::setVelocities(const std::vector& tau0, v[i] -= tau0[3 * ia + i]; v[i] /= dt; } - (*iion)->setVelocity(v[0], v[1], v[2]); - iion++; + ion->setVelocity(v[0], v[1], v[2]); ia++; } } @@ -2138,12 +2104,10 @@ void Ions::getLocalPositions(std::vector& tau) const { assert(tau.size() == 3 * local_ions_.size()); - int ia = 0; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { - (*iion)->getPosition(&tau[3 * ia]); - iion++; + ion->getPosition(&tau[3 * ia]); ia++; } } @@ -2188,12 +2152,10 @@ void Ions::setTau0() { assert(tau0_.size() == 3 * local_ions_.size()); - int ia = 0; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { - (*iion)->getPosition(&tau0_[3 * ia]); - iion++; + ion->getPosition(&tau0_[3 * ia]); ia++; } } @@ -2202,13 +2164,11 @@ void Ions::setPositionsToTau0() { assert(tau0_.size() == 3 * local_ions_.size()); - int ia = 0; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { - (*iion)->setPosition( + ion->setPosition( tau0_[3 * ia + 0], tau0_[3 * ia + 1], tau0_[3 * ia + 2]); - iion++; ia++; } } @@ -2244,13 +2204,11 @@ void Ions::getLocalForces(std::vector& tau) const { assert(tau.size() == 3 * local_ions_.size()); - int ia = 0; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { assert(3 * ia + 2 < (int)tau.size()); - (*iion)->getForce(&tau[3 * ia]); - iion++; + ion->getForce(&tau[3 * ia]); ia++; } } @@ -2351,12 +2309,10 @@ double Ions::computeMaxVlRadius() const { double radius = 0.; - std::vector::const_iterator iion = local_ions_.begin(); - while (iion != local_ions_.end()) + for (auto& ion : local_ions_) { - double r = (*iion)->computeRadiusVl(); + double r = ion->computeRadiusVl(); radius = r > radius ? r : radius; - iion++; } MGmol_MPI& mmpi(*(MGmol_MPI::instance())); @@ -2393,13 +2349,11 @@ void Ions::gatherNames(std::map& names, const int root, std::vector data; std::vector local_names; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { // get local name and index - local_names.push_back((*ion)->name()); - local_indexes.push_back((*ion)->index()); - ++ion; + local_names.push_back(ion->name()); + local_indexes.push_back(ion->index()); } // gather data to PE root @@ -2695,15 +2649,13 @@ bool Ions::hasLockedAtoms() const return (flag == 1); } -double Ions::getMaxNLradius() const +double Ions::getSpeciesMaxNLradius() const { - double radius = 0; - std::vector::const_iterator spi = species_.begin(); - while (spi != species_.end()) + double radius = 0.; + for (const auto& spi : species_) { - const double nlradius = spi->nlradius(); + const double nlradius = spi.nlradius(); radius = radius > nlradius ? radius : nlradius; - spi++; } MGmol_MPI& mmpi(*(MGmol_MPI::instance())); mmpi.allreduce(&radius, 1, MPI_MAX); @@ -2711,15 +2663,13 @@ double Ions::getMaxNLradius() const return radius; } -double Ions::getMaxLradius() const +double Ions::getSpeciesMaxLradius() const { - double radius = 0; - std::vector::const_iterator spi = species_.begin(); - while (spi != species_.end()) + double radius = 0.; + for (const auto& spi : species_) { - const double lradius = spi->lradius(); + const double lradius = spi.lradius(); radius = radius > lradius ? radius : lradius; - spi++; } MGmol_MPI& mmpi(*(MGmol_MPI::instance())); mmpi.allreduce(&radius, 1, MPI_MAX); @@ -2918,8 +2868,8 @@ void Ions::computeNumIons(void) double Ions::getMaxListRadius() const { // get radius of projectors - const double nlradius = getMaxNLradius(); - const double lradius = getMaxLradius(); + const double nlradius = getSpeciesMaxNLradius(); + const double lradius = getSpeciesMaxLradius(); double rmax = nlradius > lradius ? nlradius : lradius; @@ -3289,21 +3239,15 @@ void Ions::updateIons() Control& ct(*(Control::instance())); // update local_ions data - std::vector::iterator ion = local_ions_.begin(); - int ia = 0; - while (ion != local_ions_.end()) + int ia = 0; + for (auto& ion : local_ions_) { - (*ion)->setPosition( + ion->setPosition( tau0_[3 * ia + 0], tau0_[3 * ia + 1], tau0_[3 * ia + 2]); - // (*ion)->setForce(fion_[3*ia+0], - // fion_[3*ia+1], - // fion_[3*ia+2]); - - (*ion)->setRandomState(rand_states_[3 * ia + 0], - rand_states_[3 * ia + 1], rand_states_[3 * ia + 2]); + ion->setRandomState(rand_states_[3 * ia + 0], rand_states_[3 * ia + 1], + rand_states_[3 * ia + 2]); - ion++; ia++; } diff --git a/src/Ions.h b/src/Ions.h index 772c51d1..ae4f6adb 100644 --- a/src/Ions.h +++ b/src/Ions.h @@ -71,13 +71,14 @@ class Ions void setMapVL(); double computeMaxVlRadius() const; double computeMaxNLprojRadius() const; - double getMaxNLradius() const; - double getMaxLradius() const; - void updateListIons(); - // compute boundary for box containing all atoms to be known on local - // processor that is values for list_boundary_left_, list_boundary_right_ - void setupListIonsBoundaries(const double rmax); + /* + * Evaluate maximum pseudopotential radius among all species in class + */ + double getSpeciesMaxNLradius() const; + double getSpeciesMaxLradius() const; + + void updateListIons(); void augmentIonsData(const int nsteps, const int dir, const int disp, const int locSize, const int maxLocSize, std::vector& data, @@ -163,6 +164,8 @@ class Ions bool hasLockedAtoms() const; void clearLists(); + void rescaleVelocities(const double factor); + public: Ions(const double lat[3], const std::vector& sp); @@ -172,6 +175,10 @@ class Ions void setup(); + // compute boundary for box containing all atoms to be known on local + // processor that is values for list_boundary_left_, list_boundary_right_ + void setupListIonsBoundaries(const double rmax); + std::vector& getLocalNames() { return local_names_; } std::vector& getTau0() { return tau0_; } std::vector& getTaup() { return taup_; } @@ -334,13 +341,15 @@ class Ions void updateForcesInteractingIons(); void updateTaupInteractingIons(); - void rescaleVelocities(const double factor); /*! * Calculate minimum distance between local pairs */ double computeMinLocalSpacing() const; + void addIonToList(const Species& sp, const std::string& name, + const double crds[3], const double velocity[3], const bool lock); + // void checkUnicityLocalIons(); }; diff --git a/src/KBprojectorSparse.cc b/src/KBprojectorSparse.cc index afd8871c..937b95d2 100644 --- a/src/KBprojectorSparse.cc +++ b/src/KBprojectorSparse.cc @@ -103,7 +103,7 @@ void KBprojectorSparse::setNLindex( for (int i = 0; i < size_nl; i++) { - assert(i < lnumpt); + // assert(i < lnumpt); if ((pvec[i] < iloc * lnumpt) || (pvec[i] >= (iloc + 1) * lnumpt)) { (*MPIdata::sout) << " iloc=" << iloc << ", i=" << i @@ -960,7 +960,7 @@ bool KBprojectorSparse::setIndexesAndProjectors() // get "pvec" and "is_in_domain" int icount = get_index_array(pvec, iloc, index_low, index_high); assert(icount <= nl3); - assert(icount <= mymesh->npointsPatch()); + // assert(icount <= mymesh->npointsPatch()); assert(icount > 0); setNLindex(iloc, icount, pvec); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index cc205167..6de34971 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -209,6 +209,8 @@ add_executable(testMGkernels ${CMAKE_SOURCE_DIR}/src/pb/FDkernels.cc ${CMAKE_SOURCE_DIR}/src/tools/Timer.cc ${CMAKE_SOURCE_DIR}/tests/ut_main.cc) +add_executable(testIons + ${CMAKE_SOURCE_DIR}/tests/testIons.cc) add_executable(testGramMatrix ${CMAKE_SOURCE_DIR}/tests/testGramMatrix.cc ${CMAKE_SOURCE_DIR}/src/GramMatrix.cc @@ -334,6 +336,10 @@ add_test(NAME testBatchLaph4 add_test(NAME testtMGkernels COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/testMGkernels) +add_test(NAME testIons + COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testIons + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testGramMatrix COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/testGramMatrix) @@ -541,6 +547,7 @@ target_include_directories(testPowerDistMatrix PRIVATE ${Boost_INCLUDE_DIRS}) target_include_directories(testDensityMatrix PRIVATE ${Boost_INCLUDE_DIRS}) target_include_directories(testGramMatrix PRIVATE ${Boost_INCLUDE_DIRS}) target_include_directories(testAndersonMix PRIVATE ${Boost_INCLUDE_DIRS}) +target_include_directories(testIons PRIVATE ${Boost_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS}) target_link_libraries(testMPI PRIVATE MPI::MPI_CXX) target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES} @@ -550,6 +557,7 @@ target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX) target_link_libraries(testEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src) +target_link_libraries(testIons PRIVATE mgmol_src) if(${MAGMA_FOUND}) target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES} diff --git a/tests/testIons.cc b/tests/testIons.cc new file mode 100644 index 00000000..30b2ad0c --- /dev/null +++ b/tests/testIons.cc @@ -0,0 +1,104 @@ +#include "Control.h" +#include "Ions.h" +#include "MGmol_MPI.h" +#include "Mesh.h" +#include "Species.h" + +#include + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + + MPI_Comm comm = MPI_COMM_WORLD; + + int myrank; + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + + MGmol_MPI::setup(comm, std::cout); + Control::setup(comm, false, 0.); + + // create a domain [0.10.]^3 + const double origin[3] = { 0., 0., 0. }; + const double ll = 10.; + const double lattice[3] = { ll, ll, ll }; + const unsigned ngpts[3] = { 32, 24, 20 }; + short lap_type = 0; + + Mesh::setup(comm, ngpts, origin, lattice, lap_type); + + const double h[3] = { ll / (double(ngpts[0])), ll / (double(ngpts[1])), + ll / (double(ngpts[2])) }; + + // random number generator + static std::random_device rd; + static std::mt19937 gen(rd()); + static std::uniform_real_distribution<> dis(0.0, 1.0); + + // create one species + Species sp(MPI_COMM_WORLD); + + // read species info from pseudopotential file + std::string file_path = argv[1]; + std::string filename(file_path + "/pseudo.C_ONCV_PBE_SG15"); + std::cout << "Potential = " << filename << std::endl; + + sp.read_1species(filename); + sp.set_dim_nl(h[0]); + sp.set_dim_l(h[0]); + sp.initPotentials('f', h[0], true); + + // put species into a vector + std::vector vsp; + vsp.push_back(sp); + + Ions ions(lattice, vsp); + ions.setupListIonsBoundaries(10000.); + + double velocity[3] = { 0., 0., 0. }; + + // set "na" atoms coordinates and add them to "ions" + const int na = 10; + for (int i = 0; i < na; i++) + { + double x[3] = { origin[0] + lattice[0] * dis(gen), + origin[1] + lattice[1] * dis(gen), + origin[2] + lattice[2] * dis(gen) }; + if (myrank == 0) + std::cout << "x,y,z = " << x[0] << ", " << x[1] << ", " << x[2] + << std::endl; + + // set all x to the values of PE0 + MPI_Bcast(&x[0], 3, MPI_DOUBLE, 0, comm); + + // make a name for atom based on species and order of reading in + std::string stri = std::to_string(i); + std::string aname("C" + stri); + + ions.addIonToList(sp, aname, &x[0], velocity, false); + } + + ions.setup(); + + std::vector& new_local_ions(ions.local_ions()); + + int nlocal = new_local_ions.size(); + std::cout << "PE " << myrank << ", nlocal = " << nlocal << std::endl; + + int ntotal = 0; + MPI_Allreduce(&nlocal, &ntotal, 1, MPI_INT, MPI_SUM, comm); + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + return 1; + } + + if (ntotal != na) + { + std::cout << "ntotal = " << ntotal << std::endl; + return 1; + } + + return 0; +}