0

I am running a Monte Carlo simulation of a polymer. The entire configuration of the current state of the system is given by the object called Grid. This is my definition of Grid:

class Grid{

public:
    std::vector <Polymer> PolymersInGrid;                    // all the polymers in the grid 
    int x;                                                   // length of x-edge of grid 
    int y;                                                   // length of y-edge of grid 
    int z;                                                   // length of z-edge of grid 
    double kT;                                               // energy factor 
    double Emm_n ;                                           // monomer-solvent when Not aligned 
    double Emm_a ;                                           // monomer-solvent when Aligned
    double Ems;                                              // monomer-solvent interaction
    double Energy;                                           // energy of grid 
    std::map <std::vector <int>, Particle> OccupancyMap;     // a map that gives the particle given the location
    

    Grid(int xlen, int ylen, int zlen, double kT_, double Emm_a_, double Emm_n_, double Ems_): x (xlen), y (ylen), z (zlen), kT (kT_), Emm_n(Emm_n_), Emm_a (Emm_a_), Ems (Ems_) {        // Constructor of class
        // this->instantiateOccupancyMap(); 
    };


    // Destructor of class 
    ~Grid(){                                    

    }; 

    // assignment operator that allows for a correct transfer of properties. Important to functioning of program. 
    Grid& operator=(Grid other){
        std::swap(PolymersInGrid, other.PolymersInGrid); 
        std::swap(Energy, other.Energy); 
        std::swap(OccupancyMap, other.OccupancyMap);
        return *this; 
    } 
.
.
.
}

I can go into the details of the object Polymer and Particle, if required.

In my driver code, this is what I am going: Define maximum number of iterations.

  1. Defining a complete Grid G.
  2. Creating a copy of G called G_.
  3. I am perturbing the configuration of G_.
  4. If the perturbance on G_ is accepted per the Metropolis criterion, I assign G_ to G (G=G_).
  5. Repeat steps 1-4 until maximum number of iterations is achieved.

This is my driver code:

auto start = std::chrono::high_resolution_clock::now(); 
Grid G_ (G); 
    int acceptance_count = 0; 
    for (int i{1}; i< (Nmov+1); i++){

        // choose a move 
        G_ = MoveChooser(G, v);  

        if ( MetropolisAcceptance (G.Energy, G_.Energy, G.kT) ) {
            // accepted
            // replace old config with new config

            acceptance_count++; 
            std::cout << "Number of acceptances is " << acceptance_count << std::endl;
            G = G_;
        }


        else {
            // continue;
        }

        if (i % dfreq == 0){
            G.dumpPositionsOfPolymers (i, dfile) ;
            G.dumpEnergyOfGrid(i, efile, call) ; 
        }
        // G.PolymersInGrid.at(0).printChainCoords();


    }
    
    
    
    auto stop = std::chrono::high_resolution_clock::now(); 
    
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds> (stop-start); 

    std::cout << "\n\nTime taken for simulation: " << duration.count() << " milliseconds" << std::endl;

This is the interesting part: if I run the simulation using condition that do not have lots of "acceptances" (low temperatures, bad solvent), the simulation runs pretty fast. However, if there are a large number of acceptances, the simulation gets incredibly slow. My hypothesis is that my assignment operator = is slowing down my simulation. I ran some tests:

number of acceptances = 25365, wall-clock time = 717770 milliseconds (!)

number of acceptances = 2165, wall-clock time = 64412 milliseconds

number of acceptances = 3000, wall-clock time = 75550 milliseconds

And the trend continues. Could anyone advise me on how to possibly make this more efficient? Is there a way to bypass the slowdown I am experiencing, I think, due to the = operator?

I would really appreciate any advice you have for me!

bad_chemist
  • 283
  • 1
  • 7
  • You're copying the value, try `G = std::move(G_)` – Alan Birtles Jan 11 '22 at 07:37
  • `Grid& operator=(Grid other)` seems indeed inefficient as it call constructor, `Grid& operator=(const Grid&) = default; Grid& operator=(Grid&&) = default;` seems better (not sure why you don't copy all members, probably better to use named function instead). – Jarod42 Jan 11 '22 at 09:45
  • @Jarod42, thank you for your advice. could you possibly show an example of what you mean? – bad_chemist Jan 11 '22 at 14:42

1 Answers1

1

One thing that you can certainly do to improve performance is to force moving _G rather than coping it to G:

 G = std::move(G_);

After all, at this stage you don't need G_ any more.

Side remark. The fact that you don't need to copy all member data in operator= indicates that your design of Grid is far from perfect, but, well, keep it if the program is small and you're sure you control everything. Anyway, rather than using operator=, you should define and use a member function with a meaningful name, like "fast_and_dirty_swap" etc :-) Then you can define operator= the way suggested by @Jarod42, that is, using = default.


An alternative approach that I used before C++11 is to operate on pointers. In this scenario one would have two Grids, one "real" and one treated as a buffer, or sandbox, and on acceptance on would simply swap the pointers, so that the "buffer" filled with MoveChooser would become the real, current Grid.

A pseudocode:

  • Create two buffers, previous and current, each capable of storing a simulation state
  • Initialize current
  • Create two pointers, p_prev = &previous, p_curr = &currenrt
  • For as many steps as you wish
    • compute the next state from *p_curr and store it in *p_prev (e.g. monte_carlo_step(p_curr, p_prev)
    • swap the pointers: now the current system state is at p_curr and the previous at p_prev.
  • analyze the results stored at *p_curr
zkoza
  • 2,644
  • 3
  • 16
  • 24
  • thank you for your response! This is interesting stuff. I got rid of my operator= and destructor and use the default copy constructor. Things are faster now. However, what do you mean by pointer swapping? – bad_chemist Jan 14 '22 at 01:28
  • I create G_ after performing a move on G. In my mind G_ is the buffer I filled up with MoveChooser. On acceptance, how do I swap the pointers of G_ and G? I have these two objects floating around right now, G and G_. Do you suggest doing something along the lines of std::(&G, &G_)? – bad_chemist Jan 14 '22 at 01:32
  • No, no. I just remarked that before C++11 there was no && and people used pointers instead. You had two buffers and two pointers pointing to them. Buffers would contain previous and current states. After updating the state in one of the buffers, one would swap the pointers so that one of them would always point at the current and the other at the previous state. But this is only a remark showing an alternative approach. – zkoza Jan 14 '22 at 01:56
  • 1
    I added a pseudocode for a pointer-based alternative method. – zkoza Jan 14 '22 at 02:10
  • I can also help with self-avoiding walks, if you still need it, but I'd need more detail there: https://stackoverflow.com/questions/69913536/executing-multiple-self-avoiding-walks-recursively – zkoza Jan 14 '22 at 17:51