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 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):
- launch SAMSON and run the SAMSON Extension Generator;
- specify the path to a folder where you want to develop your SAMSON Extensions;
- name the SAMSON Extension as MonteCarlo and add some description;
- add a StateUpdater class (called
SEMonteCarloStateUpdater); - generate the SAMSON Extension.
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:
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:
And modify the constructor for the state updater class (class: SEMonteCarloStateUpdater, file: SEMonteCarloStateUpdater.cpp) by modifying an argument for the setName function as follows:
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 (Edit > 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:
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<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();
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 (Edit > Start simulation) 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 next 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 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
Modify the spin box: change its name tp spinBoxNumberOfMovingParticles modify the default and maximum value, and the single step size as follows
Add two slots: onTemperatureChanged(double) and onNumberOfMovingParticlesChanged(int):
And connect signal from two spin boxes to these 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 == nullptr) return;
ui.spinBoxNumberOfMovingParticles->setValue(
settings->loadIntValue("spinBoxNumberOfMovingParticles", 1));
ui.doubleSpinBoxTemperature->setValue(
settings->loadDoubleValue("doubleSpinBoxTemperature", 1.0));
}
void SEMonteCarloStateUpdaterProperties::saveSettings(SBGSettings *settings) {
if (settings == nullptr) 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 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 (Edit > Start simulation) and try slighlty pulling an atom.
