Bdtjtk

I am writing a molecular dynamics program that needs to take the atoms in a molecule and find the possible ways they can bond. To do this, I have a vector of Atom objects and I generate combination pairs using the following algorithm:

    void CombinationKN(std::vector<std::vector<int>> &indices, int K, int N) {

        std::string bitmask(K, 1);

        bitmask.resize(N, 0);



        do {

            /* This loop takes forever with larger N values (approx. 3000) */

            std::vector<int> indexRow;



            for (int i = 0; i < N; i++)

            {

                if (bitmask[i]) indexRow.push_back(i);

            }



            indices.push_back(indexRow);

        } while (std::prev_permutation(bitmask.begin(), bitmask.end()));

    }

It is a simple N choose K algorithm (i.e. the indices returned could contain (1, 2) but not (2, 1)) where in my case N is the number of atoms in the molecule and K is 2.

I then call the algorithm like this:

void CalculateBondGraph(const std::vector<Atom *> &atoms, std::map<int, 

    std::map<int, double>> &bondGraph, ForceField *forceField) {

    int natoms = atoms.size();



    std::vector<std::vector<int>> indices;



    utils::CombinationKN(indices, 2, natoms);



    for (auto &v : indices) {

        int i = v[0];

        int j = v[1];



        /*... Check if atoms i and j are bonded based on their coordinates.*/

    }

}

The issue with this algorithm is that it takes forever to complete for large molecules that have 3000+ atoms. I have thought about parallelizing it (specifically with OpenMP), but even then, the work would have to be split among a few threads and it would still take a lot of time to complete. I need a way to optimize this algorithm so it doesn't take so long to compute combinations. Any help is appreciated.

Thank you,
Vikas

Your CombinationKN function is way more expensive than it needs to be, if K is much smaller than N -- and if N is large then of course K is much smaller than N or you will run out of memory very quickly.

Notice that every valid index_row is a strictly monotonically increasing sequence of K integers less than N and vice-versa. It's easy enough to generate these directly:

void CombinationKN(std::vector<std::vector<int>> &indices, int K, int N) {

    std::vector<int> index_row;

    // lexographically first valid row

    for (int i=0; i<K; ++i) {

        index_row.push_back(i);

    }



    for(;;) {

        // output current row

        indeces.push_back(index_row);



        // increment index_row the the lexically next valid sequence

        // find the right-most index we can increment

        // This loop does O(1) amortized iterations if K is not large.  O(K) worst case

        int inc_index=K-1;

        int index_limit=N-1;

        while(inc_index >= 0 && index_row[inc_index] >= index_limit) {

          --inc_index;

          --index_limit;

        }

        if (inc_index < 0) {

            break; //all done

        }

        // generate the lexically first valid row with matching prefix and

        // larger value at inc_index

        int val = index_row[inc_index]+1;

        for (;inc_index<K; ++inc_index, ++val) {

            index_row[inc_index] = val;

        }

    }

}

Also, if the only thing you're doing with these combinations is iterating through them, then there's no reason to waste the (possible very large amount of) memory required to store the whole list of them. The above function contains a procedure for generating the next one from the previous one when you need it.