THREADS keyword issue

[vvishwesh]

I have been trying to run MOPAC PM6 calculations on Linux for small molecules with THREADS=4. CPU usage still shows up as 130-140% which suggests there is no paralleliization. Is the THREADS keyword now defunct.

[godotalgorithm]

The current behavior of the THREADS keyword is to make an OpenMP and/or MKL API call to request a specific number of threads (omp_set_num_threads and mkl_set_num_threads). MOPAC itself does not use any threading and never has, the parallelism of MOPAC is limited to BLAS and LAPACK calls that are the cost bottlenecks of large, conventional (non-MOZYME) calculations. MOZYME calculations have no parallelism at all.

Whether or not multiple cores are being used may depend on the size of a calculation, since those determinations are made internally by MKL and/or OpenBLAS at runtime. They might not use all available cores if matrices are too small, and calculations of very small molecules may simply spend the majority of their time outside of BLAS/LAPACK calls and in regions of the code that have no parallelism.

Keep in mind that the default behavior of MOPAC is to use all available cores on a machine. The THREADS keyword is mainly intended for restricting the number of threads when running multiple MOPAC jobs simultaneously.

[DanielWicz]

Can MOZYME be parallelized ?

[godotalgorithm]

Conceptually yes, actually maybe, pragmatically no.

Yes, the specific sparse pseudo-diagonalization algorithm used by MOZYME applies a sequence of Jacobi transformations that only interact with a few localized molecular orbitals (LMO) at a time, and many of them can be performed concurrently in principle while sparsity is (approximately) maintained. More broadly, there are other electronic structure solver algorithms that exploit locality and have similar natural concurrency [DOI:10.1103/RevModPhys.71.1085], although their practical development and deployment has been very limited.

Maybe, in that I did briefly attempt to refactor the MOZYME algorithm in isolation. I made some changes that improved parallelizability, but they also increased matrix fill-in for reasons that I could not immediately explain or identify, so I abandoned it. The basic algorithm does not guarantee matrix sparsity preservation - it depends on both the physical system being solved and the order of operations in the algorithm itself. There was an attempt from the Rocha group a few years back to switch to a different linear-scaling solver that was easier to implement on a GPU (in cuSPARSE, I think). Its baseline performance was worse, and they did not contribute their work to the open-source codebase. There is an ongoing attempt between the Rocha group and a graduate from the group who now works at AMD to add OpenMP-based threading to the MOZYME implementation in MOPAC. They've reported some progress, but they don't give me regular progress updates, so I don't know the current status. I will assist in integrating a threaded MOZYME algorithm into the open-source codebase if they are eventually successful and willing to contribute it.

No, in that the intellectual, research, and development base that built semiempirical quantum chemistry and MOPAC is functionally long dead. I'm an "outsider" (computational condensed matter physics background) that stepped in to make sure that MOPAC didn't die outright from neglect. I've tried to envision a technical future for this research area, but it's a moot point. All the development resources devoted to atomistic simulation got polarized into supporting either MM or ab-initio QM, and the middle ground represented by semiempirical models was effectively starved to death. New, performant parallel implementations are a quality-of-life feature that MOPAC development probably cannot afford at this point. MOPAC development has been part of my job in recent years, but soon it will only be an unpaid hobby for me focused strictly on maintenance and support, and MOPAC has no other focused professional developers anymore.

[DanielWicz]

Thanks for your answer @godotalgorithm !

I also tried to parallelize MOZYME. I get linear performance scaling with threads to some degree and the same energy as the single threaded version. But I'm not sure if it does not break the physics in the end. I did try to parallelize thread density/Fock, DIAGG1 and DIAGG2 with OpenMP. But I had to lock most of the variables to make it thread-safe, making it not very efficient. It gets a ~300% boost when using a six threads, but not much better.

Parallel version:
https://github.com/DanielWicz/mopac/commit/264aa64f4693c4385593f0ae0cdd6bfd99280a7e

[godotalgorithm]

Thank you @DanielWicz for your interest in MOPAC development. I'm happy to provide whatever information and insight I can to aid your efforts, and I'm happy to assist with a Pull Request if/when your contribution is ready for production. A contribution like this would still be very worthwhile even if performance saturates beyond a modest number of threads, provided that there is enough wall time savings over the serial version for typical MOZYME calculations (around 20%-30% is probably enough to be worthwhile). It's possible that this might get replaced in the future with a more scalable threading strategy, or extra switching might be needed to select between threaded versions. Also, serial operation of MOZYME would still likely be the preferred mode of operation for high-throughput workloads where many independent MOZYME calculations are required. Many simultaneous, serial jobs would likely complete the work faster than fewer, threaded calculations unless a high degree of thread scaling is achieved or memory is prematurely exhausted by the serial jobs.

Can you provide more information about the benchmark data you have posted? Is it a full MOZYME calculation, a single SCF cycle, or a specific subroutine (e.g. diagg)? What is the test molecule?

A contribution like this might need more test coverage to validate, which might delay its incorporation.

I can offer some preliminary thoughts for now. I probably can't examine this carefully until March as I have deadlines on other work.

diagg1 should conceptually be very amenable to parallelization. It is projecting the Fock matrix into the LMO basis, which is a product of 3 sparse matrices in a variant of compressed row/column storage format. Each of the final matrix elements can be computed independently. I would similarly expect density_for_MOZYME to be amenable to parallelization. It is diagg2 specifically that is very poorly designed for parallelism, because Jacobi rotations cannot be freely reordered. Two sparse Jacobi rotations can be reordered or performed simultaneously if they don't have overlapping rows or columns, but I don't think that MOZYME has any structure for detecting that right now because it isn't relevant to the serial implementation. If collisions are avoided, reordered Jacobi rotations can still effectively decouple the occupied and virtual spaces, but the specific LMOs will differ in the converged solution. In exact arithmetic, these differences do not matter because the LMOs will span the same occupied subspace, which is all that matters here. However, with floating-point arithmetic and aggressive pruning of small matrix elements, the answers will start to differ and the numerics will need closer examination. Collisions between Jacobi rotations will cause more serious mathematical problems, and so I'll need to examine how you are avoiding that more carefully.

[DanielWicz]

@godotalgorithm Thanks for the detailed answer! I won’t have much time until next week, but here is a quick response:

1. It is a measure of the entire SCF calculation (wall time to reach SCF convergence within 1e-5). I also have a breakdown into different parts and their scaling (DIAGG1, DIAGG2, density matrix, Fock).
2. The target system is a protein with 2,850 atoms.
3. The scaling stops at 6 threads, which may be because my CPU has 6 performance cores and 10 efficiency cores.