Loading...
Searching...
No Matches
Programming new state updaters: Monte Carlo

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 Extension created in this tutorial can be found at https://github.com/1A-OneAngstrom/SAMSON-Developer-Tutorials.

Generating an Extension

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

The SAMSON Extension 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 Extension

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

SB_CLASS_DESCRIPTION("DevTutorial: Monte Carlo");
#define SB_CLASS_DESCRIPTION(CLASS_DESCRIPTION)
Declares the description of the class.
Definition: SBCClassProxy.hpp:202

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

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 Extension 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;
};
This template class defines physical quantity types.
Definition: SBDQuantityType.hpp:43
This class implements a random number generator.
Definition: SBDTypeRandom.hpp:17
This class is the base class for state updaters of particle systems simulators.
Definition: SBSStateUpdaterParticleSystem.hpp:21

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);
static SBCTime getTime()
Returns SAMSON's internal time.
Definition: SAMSON.cpp:713
SBDQuantityType< SBDQuantityUnitType< SBUnitSystemSI, -10, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 > > angstrom
The angstrom type.
Definition: SBDQuantity.hpp:81
SBDQuantityType< SBDQuantityUnitType< SBUnitSystemSI, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 > > kelvin
The kelvin type.
Definition: SBDQuantity.hpp:352

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<SBAtom> const* particleIndexer = (*dynamicalModel)->getAtomIndexer();
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();
SBDQuantityType< SBDQuantityUnitType< typename Unit::SystemType, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 >, Value > exp(const SBDQuantityType< Unit, Value > &q)
Returns the exponential of physical quantity q (for dimensionless physical quantities only)
Definition: SBDQuantityType.hpp:1136
static const SBDQuantityType< SBDQuantityUnitType< SBUnitSystemSI, -12, 2, -27, 1, -15, -2, 0, 0, 0, -1, 0, 0, 0, 0 > > *const boltzmannConstant
The Boltzmann constant.
Definition: SBDQuantityConstant.hpp:32

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 Extension, 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;
}
kelvin temperature
The temperature type.
Definition: SBDQuantity.hpp:709

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 Extension 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);
};
This class describes a widget that allows users to view and potentially edit properties of data graph...
Definition: SBGDataGraphNodeProperties.hpp:12

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());
}
This class is used to save and load settings.
Definition: SBGSettings.hpp:17
int loadIntValue(const QString &name, int defaultValue=0)
Loads an integer.
Definition: SBGSettings.cpp:82
void saveValue(const QString &name, QVariant value)
Saves a value.
Definition: SBGSettings.cpp:58
qreal loadDoubleValue(const QString &name, double defaultValue=0)
Loads a double.
Definition: SBGSettings.cpp:90

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 Extension, 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