OpenMP parallelization of energy #364

directysj · 2021-05-12T14:25:42Z

directysj
May 12, 2021

I would like to simulate really large system (approximately 20000 particles). Therefore, I made parallel energy.cpp of Faunus, but when running the code. It is using a single thread. It would be really helpful for your suggestions in the code, providing below section wise.

void PolicyIonIon::updateBox(EwaldData &d, const Point &box) const {
assert(d.policy == EwaldData::PBC or d.policy == EwaldData::PBCEigen);
d.box_length = box;
int n_cutoff_ceil = ceil(d.n_cutoff);
d.check_k2_zero = 0.1 * std::pow(2 * pc::pi / d.box_length.maxCoeff(), 2);
int k_vector_size = (2 * n_cutoff_ceil + 1) * (2 * n_cutoff_ceil + 1) * (2 * n_cutoff_ceil + 1) - 1;
if (k_vector_size == 0) {
d.k_vectors.resize(3, 1);
d.Aks.resize(1);
d.k_vectors.col(0) = Point(1, 0, 0); // Just so it is not the zero-vector
d.Aks[0] = 0;
d.num_kvectors = 1;
d.Q_ion.resize(1);
d.Q_dipole.resize(1);
} else {
double nc2 = d.n_cutoff * d.n_cutoff;
d.k_vectors.resize(3, k_vector_size);
d.Aks.resize(k_vector_size);
d.num_kvectors = 0;
d.k_vectors.setZero();
d.Aks.setZero();
int start_value = 1;

			double factor,dnx2,dny2,k2;
			int ny,nz;

    #pragma omp parallel for firstprivate(nc2,start_value,n_cutoff_ceil)\
			private(ny,nz,factor,dnx2,dny2,k2) schedule(dynamic)
    
			for (int nx = 0; nx <= n_cutoff_ceil; nx++) {
        dnx2 = double(nx * nx);
        factor = (nx > 0) ? 2.0 : 1.0; // optimization of PBC Ewald (and
                                              // always the case for IPBC Ewald)
        //for (int ny = -n_cutoff_ceil * start_value; ny <= n_cutoff_ceil; ny++) {
        for ( ny = -n_cutoff_ceil * start_value; ny <= n_cutoff_ceil; ny++) {
            dny2 = double(ny * ny);
            for ( nz = -n_cutoff_ceil * start_value; nz <= n_cutoff_ceil; nz++) {
                Point kv = 2 * pc::pi * Point(nx, ny, nz).cwiseQuotient(d.box_length);
                k2 = kv.squaredNorm() + d.kappa_squared; // last term is only for Yukawa-Ewald
                if (k2 < d.check_k2_zero)                       // Check if k2 != 0
                    continue;
                if (d.use_spherical_sum) {
                    double dnz2 = double(nz * nz);
                    if ((dnx2 + dny2 + dnz2) / nc2 > 1)
                        continue;
                }
                d.k_vectors.col(d.num_kvectors) = kv;
                d.Aks[d.num_kvectors] = factor * exp(-k2 / (4 * d.alpha * d.alpha)) / k2;
                d.num_kvectors++;
            }
        }
    }
    d.Q_ion.resize(d.num_kvectors);
    d.Q_dipole.resize(d.num_kvectors);
    d.Aks.conservativeResize(d.num_kvectors);
    d.k_vectors.conservativeResize(3, d.num_kvectors);
}

}

void PolicyIonIon::updateComplex(EwaldData &data, Space::Tgvec &groups) const {

double qr;

#pragma omp parallel for private(qr) schedule(dynamic)

for (int k = 0; k < data.k_vectors.cols(); k++) {
    const Point &q = data.k_vectors.col(k);
    EwaldData::Tcomplex Q(0, 0);
    for (auto &g : groups) {       // loop over molecules
        for (auto &particle : g) { // loop over active particles
            //double qr = q.dot(particle.pos);
            qr = q.dot(particle.pos);
            Q += particle.charge * EwaldData::Tcomplex(std::cos(qr),
                                                       std::sin(qr)); // 'Q^q', see eq. 25 in ref.
        }
        data.Q_ion[k] = Q;
    }
}

}

void PolicyIonIon::updateComplex(EwaldData &d, Change &change, Space::Tgvec &groups, Space::Tgvec &oldgroups) const {
assert(groups.size() == oldgroups.size());

#pragma omp parallel for schedule(dynamic)

for (int k = 0; k < d.k_vectors.cols(); k++) {
    auto &Q = d.Q_ion[k];
    const Point &q = d.k_vectors.col(k);

    for (auto &changed_group : change.groups) {
        auto &g_new = groups.at(changed_group.index);
        auto &g_old = oldgroups.at(changed_group.index);
        for (auto i : changed_group.atoms) {
            if (i < g_new.size()) {
                double qr = q.dot(g_new[i].pos);
                Q += g_new[i].charge * EwaldData::Tcomplex(std::cos(qr), std::sin(qr));
            }
            if (i < g_old.size()) {
                double qr = q.dot(g_old[i].pos);
                Q -= g_old[i].charge * EwaldData::Tcomplex(std::cos(qr), std::sin(qr));
            }
        }
    }
}

}

void PolicyIonIonIPBC::updateBox(EwaldData &data, const Point &box) const {
assert(data.policy == EwaldData::IPBC or data.policy == EwaldData::IPBCEigen);
data.box_length = box;
int ncc = std::ceil(data.n_cutoff);
data.check_k2_zero = 0.1 * std::pow(2 * pc::pi / data.box_length.maxCoeff(), 2);
int k_vector_size = (2 * ncc + 1) * (2 * ncc + 1) * (2 * ncc + 1) - 1;
if (k_vector_size == 0) {
data.k_vectors.resize(3, 1);
data.Aks.resize(1);
data.k_vectors.col(0) = Point(1, 0, 0); // Just so it is not the zero-vector
data.Aks[0] = 0;
data.num_kvectors = 1;
data.Q_ion.resize(1);
data.Q_dipole.resize(1);
} else {
double nc2 = data.n_cutoff * data.n_cutoff;
data.k_vectors.resize(3, k_vector_size);
data.Aks.resize(k_vector_size);
data.num_kvectors = 0;
data.k_vectors.setZero();
data.Aks.setZero();
int start_value = 0;

			double xfactor,yfactor,factor;
			double dnx2,dny2,dnz2,k2;

    #pragma omp parallel for firstprivate(nc2,start_value)\
			private(xfactor,yfactor,factor,dnx2,dny2,dnz2,k2) schedule(dynamic)

    for (int nx = 0; nx <= ncc; nx++) {
        dnx2 = double(nx * nx);
        xfactor = (nx > 0) ? 2.0 : 1.0; // optimization of PBC Ewald
        for (int ny = -ncc * start_value; ny <= ncc; ny++) {
            dny2 = double(ny * ny);
            yfactor = (ny > 0) ? 2.0 : 1.0; // optimization of PBC Ewald
            for (int nz = -ncc * start_value; nz <= ncc; nz++) {
                factor = xfactor * yfactor;
                if (nz > 0)
                    factor *= 2;
                Point kv = 2 * pc::pi * Point(nx, ny, nz).cwiseQuotient(data.box_length);
                k2 = kv.squaredNorm() + data.kappa_squared; // last term is only for Yukawa-Ewald
                if (k2 < data.check_k2_zero)                       // Check if k2 != 0
                    continue;
                if (data.use_spherical_sum) {
                    dnz2 = double(nz * nz);
                    if ((dnx2 + dny2 + dnz2) / nc2 > 1)
                        continue;
                }
                data.k_vectors.col(data.num_kvectors) = kv;
                data.Aks[data.num_kvectors] = factor * exp(-k2 / (4 * data.alpha * data.alpha)) / k2;
                data.num_kvectors++;
            }
        }
    }
    data.Q_ion.resize(data.num_kvectors);
    data.Q_dipole.resize(data.num_kvectors);
    data.Aks.conservativeResize(data.num_kvectors);
    data.k_vectors.conservativeResize(3, data.num_kvectors);
}

}

void PolicyIonIonIPBC::updateComplex(EwaldData &d, Space::Tgvec &groups) const {
assert(d.policy == EwaldData::IPBC or d.policy == EwaldData::IPBCEigen);

#pragma omp parallel for schedule(dynamic)

for (int k = 0; k < d.k_vectors.cols(); k++) {
    const Point &q = d.k_vectors.col(k);
    EwaldData::Tcomplex Q(0, 0);
    for (auto &g : groups) {
        for (auto &particle : g) {
            Q += q.cwiseProduct(particle.pos).array().cos().prod() * particle.charge; // see eq. 2 in doi:10/css8
        }
    }
    d.Q_ion[k] = Q;
}

}

void PolicyIonIonIPBC::updateComplex(EwaldData &d, Change &change, Space::Tgvec &groups,
Space::Tgvec &oldgroups) const {
assert(d.policy == EwaldData::IPBC or d.policy == EwaldData::IPBCEigen);
assert(groups.size() == oldgroups.size());

#pragma omp parallel for schedule(dynamic)

for (int k = 0; k < d.k_vectors.cols(); k++) {
    auto &Q = d.Q_ion[k];
    const Point &q = d.k_vectors.col(k);
    for (auto &changed_group : change.groups) {
        auto &g_new = groups.at(changed_group.index);
        auto &g_old = oldgroups.at(changed_group.index);
        for (auto i : changed_group.atoms) {
            if (i < g_new.size())
                Q += q.cwiseProduct(g_new[i].pos).array().cos().prod() * g_new[i].charge;
            if (i < g_old.size())
                Q -= q.cwiseProduct(g_old[i].pos).array().cos().prod() * g_old[i].charge;
        }
    }
}

}
`

mlund · 2021-05-14T07:26:08Z

mlund
May 14, 2021
Maintainer

Thanks for the suggested OpenMP additions to the Ewald code. I wonder if you could provide a bit more detail about the system you want to simulate? My immediate thought is that the bottle-neck is the real-space part, i.e. Energy::Nonbonded. Could you try to profile this, for example by running your system with and without Ewald summation? (coulomb type plain vs. ewald) Another major speed increase would likely come from a cell-list which is currently in the works -- see pull-request (#346).

48 replies

directysj Aug 25, 2021
Author

Interestingly, If only OPENMP is ON (with TBB OFF), then also I am getting same warnings as

[faunus 19057] W: 'openmp' summation unavailable; falling back to 'serial'
[faunus 19057] D: setting summation policy to "serial"

...........

mlund Aug 25, 2021
Maintainer

Yes, I have seen this before. I think it's because some versions of GCC require TBB for it's openmp support. The warning is triggered in src/energy.h:1335 when _OPENMP (set by the compiler) is not defined.

directysj Aug 30, 2021
Author

Yes, I have seen this before. I think it's because some versions of GCC require TBB for it's openmp support. The warning is triggered in src/energy.h:1335 when _OPENMP (set by the compiler) is not defined.

Hello @mlund ,

I have installed LLVM 12 (which gives CLANG/CLANG++) and Intel's TBB, locally. I have modified the CMAKE command from my previous build as:

$CC=$HOME/llvm/bin/clang CXX=$HOME/llvm/bin/clang++ cmake -B build -DCMAKE_CXX_FLAGS=-stdlib=libc++ -DCMAKE_EXE_LINKER_FLAGS=-stdlib=libc++ -DENABLE_OPENMP=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX:PATH="..." -DZLIB_LIBRARY=/usr/lib/x86_64-linux-gnu/libz.so

I am using CMAKE 3.20 and Anaconda3 for python. Before CMAKE, I am exporting compiler flags as:
$export LDFLAGS="-L$HOME/llvm/lib $LDFLAGS"; export CFLAGS="-I$HOME/llvm/include $CFLAGS"; export CC=$HOME/llvm/bin/clang; export CXX=$HOME/llvm/bin/clang++

Now, I can configure and make like shown below:
........

[ 98%] Building CXX object src/CMakeFiles/faunus.dir/faunus.cpp.o
[100%] Linking CXX executable ../faunus
/usr/bin/ld: CMakeFiles/faunus.dir/faunus.cpp.o: in function `createProgressTracker(bool, unsigned int)':
faunus.cpp:(.text+0xcc6f): undefined reference to `ProgressIndicator::ProgressBar::ProgressBar(unsigned int, std::__1::basic_ostream<char, std::__1::char_traits<char> >&, unsigned int, char, char)'

...............................
And, then exit with clang-12: error: linker command failed with exit code 1 (use -v to see invocation).

I don't have any idea, how to go from here. Can I ask you to suggest something?
Thank you in advance!!

~Jalim Singh

directysj Sep 9, 2021
Author

Yes, I have seen this before. I think it's because some versions of GCC require TBB for it's openmp support. The warning is triggered in src/energy.h:1335 when _OPENMP (set by the compiler) is not defined.

Hello @mlund ,

I have installed LLVM 12 (which gives CLANG/CLANG++) and Intel's TBB, locally. I have modified the CMAKE command from my previous build as:

$CC=$HOME/llvm/bin/clang CXX=$HOME/llvm/bin/clang++ cmake -B build -DCMAKE_CXX_FLAGS=-stdlib=libc++ -DCMAKE_EXE_LINKER_FLAGS=-stdlib=libc++ -DENABLE_OPENMP=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX:PATH="..." -DZLIB_LIBRARY=/usr/lib/x86_64-linux-gnu/libz.so

I am using CMAKE 3.20 and Anaconda3 for python. Before CMAKE, I am exporting compiler flags as:
$export LDFLAGS="-L$HOME/llvm/lib $LDFLAGS"; export CFLAGS="-I$HOME/llvm/include $CFLAGS"; export CC=$HOME/llvm/bin/clang; export CXX=$HOME/llvm/bin/clang++

Now, I can configure and make like shown below:
........
[ 98%] Building CXX object src/CMakeFiles/faunus.dir/faunus.cpp.o
[100%] Linking CXX executable ../faunus
/usr/bin/ld: CMakeFiles/faunus.dir/faunus.cpp.o: in function `createProgressTracker(bool, unsigned int)':
faunus.cpp:(.text+0xcc6f): undefined reference to `ProgressIndicator::ProgressBar::ProgressBar(unsigned int, std::__1::basic_ostream<char, std::__1::char_traits<char> >&, unsigned int, char, char)'
...............................
And, then exit with clang-12: error: linker command failed with exit code 1 (use -v to see invocation).

I don't have any idea, how to go from here. Can I ask you to suggest something?
Thank you in advance!!

~Jalim Singh

Dear @mlund ,

I have looked the problem deeper and it seems that cmake is not passing some of library flags to the linker. This I found when Cmake BOOL option is kept on and in the detail it shows many times -std=gnu++17 instead of -std=c++17. Is this the problem of cmake OR the files related to cmake in this project? Can I ask you to have a look at it?

Thanking you in anticipation!!

mlund Feb 27, 2023
Maintainer

The development environment can unfortunately be tricky to setup. I have not been able to reproduce the errors above. If nothing else works, I recommend to use the scripts/Dockerfile (see Installation about Docker in the manual). Alternatively, use Remote Containers with Visual Code Studio which will setup a working environment in a Docker image on all platforms that can run Docker (one can now click in a badge on the Faunus github site).

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

OpenMP parallelization of energy #364

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 48 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

OpenMP parallelization of energy #364

directysj May 12, 2021

Replies: 1 comment · 48 replies

mlund May 14, 2021 Maintainer

directysj Aug 25, 2021 Author

mlund Aug 25, 2021 Maintainer

directysj Aug 30, 2021 Author

directysj Sep 9, 2021 Author

mlund Feb 27, 2023 Maintainer

directysj
May 12, 2021

Replies: 1 comment 48 replies

mlund
May 14, 2021
Maintainer

directysj Aug 25, 2021
Author

mlund Aug 25, 2021
Maintainer

directysj Aug 30, 2021
Author

directysj Sep 9, 2021
Author

mlund Feb 27, 2023
Maintainer