In this tutorial, we show how to implement a simple exploration method based on a Metropolis Monte Carlo approach.

Note: This tutorial shows a simplified implementation.

The source code for an Element created in this tutorial can be found at https://github.com/1A-OneAngstrom/SAMSON-Developer-Tutorials.

Generating an Element

We will start by creating a new SAMSON Element called MonteCarlo thanks to the SAMSON Element generator (please, refer to the SAMSON Element generator tutorial for a reminder on how to use it):

launch SAMSON and run the SAMSON Element generator;
specify the path to a folder where you want to develop your SAMSON Elements;
name the SAMSON Element as MonteCarlo and add some description;
add a StateUpdater class (called SEMonteCarloStateUpdater);
generate the SAMSON Element.

The SAMSON Element generator generates a state updater class (derived from the SBSStateUpdaterParticleSystem class) and a property widget class for the state updater.

Setting the description of the Element

Open the header file of your state updater descriptor (file: SEMonteCarloStateUpdaterDescriptor.hpp) and modify the class description:

SB_CLASS_DESCRIPTION("DevTutorial: Monte Carlo");

Do the same for the property window descriptor (file: SEMonteCarloStateUpdaterPropertiesDescriptor.hpp). Modify the getName function of the property window (file: SEMonteCarloStateUpdaterProperties.cpp) to have a user-friendly description of your interaction model inside SAMSON:

QString SEMonteCarloStateUpdaterProperties::getName() const {
    return "DevTutorial: Monte Carlo";
}

And modify the constructor for the state updater class (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdater.cpp) by modifying an argument for the setName function as follows:

setName("DevTutorial: Monte Carlo");

Now, try building the SAMSON Element.

Start SAMSON and open any file containing atoms or add atoms in the data graph, then create a new simulator (Simulation > Add simulator or ctrl/ Cmd⌘ + shift + M): your new Element should be in the list of the available state updaters.

Implementation

Open the state updater header file (SEMonteCarloStateUpdaterProperties.hpp) and include the SBRandom.hpp header:

#include "SBRandom.hpp"

#include "SAMSON.hpp"

Add the following private variables to the state updater class definition (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdaterProperties.hpp):

class SEMonteCarloStateUpdater : public SBSStateUpdaterParticleSystem {
    // ...
private:
    // Monte Carlo variables
    unsigned int                numberOfRejects;
    unsigned int                numberOfTrials;
    SBRandom                    randomNumberGenerator;
    // Metropolis variables
    unsigned int                numberOfMovingParticles;
    SBQuantity::angstrom        maximumDisplacement;
    SBQuantity::temperature     temperature;
};

We initialize some of these variables in the constructor of the state updater class (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdaterProperties.cpp):

numberOfRejects = numberOfTrials = 0;
randomNumberGenerator.seed(SAMSON::getTime());
numberOfMovingParticles = 5;
maximumDisplacement = SBQuantity::angstrom(0.1);
temperature = SBQuantity::kelvin(50);

Note: We set the seed for the random number generator based on the current SAMSON time. If you want to get reproductive results (e.g. for Debugging purposes) you can set the seed to 0 or any other constant.

Now, let's implement the updateState function of the state updater class (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdaterProperties.cpp).

First, we get the number of particles in the particle system:

SBPointerIndexer<SBStructuralParticle> const* particleIndexer = (*dynamicalModel)->getStructuralParticleIndexer();
unsigned int nParticles = particleIndexer->size(); // number of particles in the particle system
// check the maximum number of particles
numberOfMovingParticles = (numberOfMovingParticles > nParticles) ? nParticles : numberOfMovingParticles;

We request the interaction model to update the interactions and get the energy of the particle system, which then will be used in the algorithm. Since we do not require forces we ask the interaction model to flush the force buffer.

// compute the current energy
(*interactionModel)->updateInteractions();
SBQuantity::energy currentEnergy = (*interactionModel)->getEnergy();
// notify the interaction model that we have used the list of updated forces since we do not need it
(*interactionModel)->flushForceBuffer();

We randomly determine the list of particles that will be moved by the algorithm:

// randomly generate the list of moving particles
std::vector<unsigned int> movingParticleIndexer;
while (movingParticleIndexer.size() != numberOfMovingParticles) {
    // generate an index of a particle and add it in the indexer of moving particles
    unsigned int idx = (unsigned int)(randomNumberGenerator.randUnsignedLong() % nParticles);
    movingParticleIndexer.push_back(idx);
}

We remember the current positions of these particles and compute the new ones based on a random displacement normalized by the maximum possible displacement:

// move particles incrementally
std::vector<SBPosition3> currentPositionVector;
for (unsigned int i = 0; i < numberOfMovingParticles; i++) {
    unsigned int currentMovingParticleIndex = movingParticleIndexer[i];
    // store the current position
    SBPosition3 currentPosition = (*dynamicalModel)->getPosition(currentMovingParticleIndex);
    currentPositionVector.push_back(currentPosition);
    // create the random displacement normalized by the maximumDisplacement
    SBPosition3 displacement((randomNumberGenerator.randDouble1() - 0.5) * maximumDisplacement,
                             (randomNumberGenerator.randDouble1() - 0.5) * maximumDisplacement,
                             (randomNumberGenerator.randDouble1() - 0.5) * maximumDisplacement);
    SBPosition3 newPosition = currentPosition + displacement;
    // set the new position of particle i
    (*dynamicalModel)->setPosition(currentMovingParticleIndex, newPosition);
}

We request the interaction model to update the interactions and get the new energy of the particle system. If the new energy is less then the current one, we accept all the displacements, else we reject the move for each displaced particle based on the probability.

// compute the new energy
(*interactionModel)->updateInteractions();
SBQuantity::energy newEnergy = (*interactionModel)->getEnergy();
// accept or reject
if (newEnergy > currentEnergy) {
    double p = randomNumberGenerator.randDouble1();
    double b = exp( ( (currentEnergy - newEnergy) /
                      ((*SBConstant::boltzmannConstant) * temperature) ).getValue() );
    // check the probability
    if (p > b) {
        // reject the move
        for (unsigned int i = 0; i < numberOfMovingParticles; i++) {
            unsigned int currentMovingParticleIndex = movingParticleIndexer[i];
            (*dynamicalModel)->setPosition(currentMovingParticleIndex, currentPositionVector[i]);
        }
        numberOfRejects++;
    }
}
numberOfTrials++;
// update the structural model based on the new positions in the dynamical model
(*dynamicalModel)->updateStructuralNodes();

After the certain number of trials (e.g. 100), we decide whether to increase or to decrease the maximum possible displacement for a particle based on the ratio of the number of rejects to the number of trials; and we reset the number of trials and rejects.

// check for the rejects ratio after a certain number of trials
// and modify the maximum displacement
if (numberOfTrials == 100) {
    // check the rejects ratio
    if ( (double)numberOfRejects / (double)numberOfTrials > 0.5 ) {
        // rejects constitute to more than a half of trials
        // decrease the maximum displacement
        maximumDisplacement *= 0.95;
        // but not less than 0.01 Å
        maximumDisplacement = (maximumDisplacement > SBQuantity::angstrom(0.01)) ?
                    SBQuantity::angstrom(0.01) : maximumDisplacement;
    }
    else {
        // rejects constitute to less (or equal) than a half of trials
        // increase the maximum displacement by 5%
        maximumDisplacement *= 1.05;
    }
    // reset the number of trials and rejects
    numberOfRejects = 0;
    numberOfTrials = 0;
}

Now, you should have the basic implementation of a simple exploration method based on a Metropolis Monte Carlo approach. You can build the Element, launch SAMSON, load or create any molecule, and add a new simulator with this state updater and with, for example, the springs interaction model. Run the simulation (Simulation > Start simulation or Space or click on the Play icon) and try slighlty pulling an atom.

But for now, we have not exposed the possibility for the user to modify parameters of the method (e.g. temperature and number of moving particles). Let's do this in the Exposing parameters section.

Exposing parameters

We want to make it possible for the user to modify temperature and the number of moving particles. For that we will need to modify both state updater and state updater property window classes (SEMonteCarloStateUpdater and SEMonteCarloStateUpdaterProperties classes).

First, let's create setter functions in the state updater class (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdater.hpp)

class SEMonteCarloStateUpdater : public SBSStateUpdaterParticleSystem {
    // ...
    
    /// \name Setters
    //@{
    void    setNumberOfMovingParticles(int n);
    void    setTemperature(SBQuantity::temperature T);
    //@}
    
};

We implement these functions as follows (with checks for negative values):

void SEMonteCarloStateUpdater::setNumberOfMovingParticles(int n) {
    numberOfMovingParticles = n < 0 ? numberOfMovingParticles : n;
}
void SEMonteCarloStateUpdater::setTemperature(SBQuantity::temperature T){
    temperature = T < SBQuantity::temperature(0.0) ? temperature : T;
}

We check for the maximum number of particles in the updateState function of the state updater.

Let's now modify the property window class.

In the current GUI of the Element we have the step size and the number of steps which we do not actually use in our state updater, we will modify these to allow the user to change the temperature and the number of moving particles. You might need to modify the maximum size of the widget and lables.

Modify the double spin box: change its name to doubleSpinBoxTemperature, modify the default and maximum value, and the single step size as follows

A double spin box for temperature

Modify the spin box: change its name tp spinBoxNumberOfMovingParticles modify the default and maximum value, and the single step size as follows

A spin box for the number of moving particles

Add two slots: onTemperatureChanged(double) and onNumberOfMovingParticlesChanged(int):

Add slots

And connect signal from two spin boxes to these slots:

Connect signals to slots

Add these public slots to the state updater property window class declaration (class: SEMonteCarloStateUpdaterProperties, file: SEMonteCarloStateUpdaterProperties.hpp):

class SEMonteCarloStateUpdaterProperties : public SBGDataGraphNodeProperties {
    // ...
    
public slots:
    void    onNumberOfMovingParticlesChanged(int);
    void    onTemperatureChanged(double);
    
};

And implement these functions as follows (file: SEMonteCarloStateUpdaterProperties.cpp):

void SEMonteCarloStateUpdaterProperties::onNumberOfMovingParticlesChanged(int i) {
    stateUpdater->setNumberOfMovingParticles(i);
}
void SEMonteCarloStateUpdaterProperties::onTemperatureChanged(double d) {
    stateUpdater->setTemperature(SBQuantity::kelvin(d));
}

Modify the loadSettings and loadSettings functions to save the GUI state from one session to another:

void SEMonteCarloStateUpdaterProperties::loadSettings(SBGSettings *settings) {
    if (settings == 0) return;
    ui.spinBoxNumberOfMovingParticles->setValue(settings->loadIntValue("spinBoxNumberOfMovingParticles", 1));
    double v = settings->loadDoubleValue("doubleSpinBoxTemperature", 1.0);
    ui.doubleSpinBoxTemperature->setValue(v);
}
void SEMonteCarloStateUpdaterProperties::saveSettings(SBGSettings *settings) {
    if (settings == 0) return;
    settings->saveValue("doubleSpinBoxTemperature", ui.doubleSpinBoxTemperature->value());
    settings->saveValue("spinBoxNumberOfMovingParticles", ui.spinBoxNumberOfMovingParticles->value());
}

In both setup functions of the state updater property window class (class: SEMonteCarloStateUpdaterProperties) modify the setting of steps size and number of steps to set temperature and number of moving particles as follows:

stateUpdater->setTemperature(SBQuantity::kelvin(ui.doubleSpinBoxTemperature->value()));

stateUpdater->setNumberOfMovingParticles(ui.spinBoxNumberOfMovingParticles->value());

Now, when the user changes temperature or the number of moving particles in the property window it will automaticaly use the new parameters in the updateState function of the state updater.

Build the Element, launch SAMSON, load or create any molecule and add a new simulator with this state updater and with, for example, the springs interaction model. Run the simulation (Simulation > Start simulation or Space or click on the Play icon) and try slighlty pulling an atom.

Applying a simulator with the developed state updater