OpenMP parallelization of energy

[directysj]

Hello @mlund

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]

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).

[directysj]

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]

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]

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]

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]

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).