Programming new simulation methods

Introduction

 

This tutorial will guide you through the different stages required to update the state of a set of atoms using a state updater. State updater are parts of the simulation engine in SAMSON. They can be used to implement integrators (like Verlet) or optimisers (like the steepest descent algorithm).

Relations between the elements inside a simulator

 

Inside a simulator the state updater and the interaction model share the same dynamical model. Each time the simulator requires an update through the updateState method, the state updater moves the atoms in the dynamical model according to the forces given by the interaction model.

 

How to add a new state updater?

 

Use the Helper to create a new element in your SE folder (4 steps):

Step 1: select your SE folder

Step 2: set the name/description fields for your state updater

 

Step 3: choose state updater for the class type

 

Step 4: generate this new element

Once this element has been generated in your SE folder update your project using the cmake tools. Then build this new element using your favourite IDE.

Start SAMSON and open a file containing atoms (samx/pdb/xyz) then create a new simulator (ctrl/cmd+shift+m) : your new element should be in the list of the state updater available.

Create a new simulator

 

How to move atoms using this state updater?

1 101 to make things move

According to the graphic presented at the beginning of this tutorial we need to modify the method called updateState in the file SE*yourname*StateUpdater.cpp. Normally, this method contains some commented code that allows you to implement a steepest descent if you uncomment it. Leave this code commented for the moment.

In order to move the atoms we need to modify the position of these atoms. To do so, we first need to get the numbers of atoms that are inside the dynamical model.

    SBPointerIndexer const* particleIndexer = (*dynamicalModel)->getStructuralParticleIndexer();
    unsigned int nParticles = particleIndexer->size(); // number of particles in the particle system

Now we can iterate over the atoms and retrieve the positions of the particles inside the dynamical model. If you work on a integrator for MD you can also retrieve / modify the momentum of the particle using the get/setMomentum method contained in the dynamical model.

        SBPosition3 currentPosition = (*dynamicalModel)->getPosition(i); // get the position of particle i

Then we can create a displacement using SAMSON quantities:

        SBPosition3 displacement = SBPosition3(SBQuantity::length(1),
                                               SBQuantity::length(0),
                                               SBQuantity::length(0));

Finally we add this shift to the current position and we update the position of this particle:

        SBPosition3 newPosition = currentPosition + displacement; // change the current position
 
        (*dynamicalModel)->setPosition(i, newPosition); // set the new position of particle i

We can now compile this code and run SAMSON. If we create a new simulator with this state updater the atoms will be shifted on the x axis.

2 Steepest Descent

For this example just copy paste the piece of code below (some sdk mistakes are fixed in the sample below as SBPointerIndex became SBPointerIndexer). You can see that as we use the forces that come from the interaction model we need to multiply the SBForce3 by the time step two times in order to be coherent in terms of units and get a SBPosition3.

SBQuantity::time stepSize = getStepSize();
 
SBPointerIndexer const* particleIndexer = (*dynamicalModel)->getStructuralParticleIndex();
unsigned int nParticles = particleIndexer->;size();
 
unsigned int nSteps = getNumberOfSteps();
for (unsigned int k = 0; k < nSteps; k++) { (*interactionModel)->updateInteractions();
    //(*dynamicalModel)->flushPositionBuffer();
 
    for (unsigned int i = 0; i < nParticles; i++) { SBPosition3 oldPosition = (*dynamicalModel)->;getPosition(i);
        SBForce3 force = (*interactionModel)->getForce(i);
        SBPosition3 deltaPosition = 1.0 / SBQuantity::mass(1)*stepSize*stepSize*force;
        SBPosition3 newPosition = oldPosition + deltaPosition;
        (*dynamicalModel)->setPosition(i, newPosition);
 
    }
 
    //(*interactionModel)->flushForceBuffer();
 
}

 

Case study: a simple monte carlo simulator

In this last example we will implement a simple exploration method based on a Metropolis Monte Carlo approach. To do so we will need to declare some variable in the file SE*yourname*StateUpdater.hpp.

private:
    // Monte Carlo variables
    unsigned int numberOfRejects;
    unsigned int numberOfTrials;
    SBRandom randomNumberGenerator;
 
    // Metropolis variables
    unsigned int numberOfMovingParticles;
    SBQuantity::angstrom maximumDisplacement;
    SBQuantity::temperature temperature;
    SBQuantityInverse::Type beta;
 
    // Dynamical model variables
    SBPointerIndexer const* particleIndexer;
    unsigned int nParticles;

We will init some of these variables in the constructor of the class in the file SE*yourname*StateUpdater.cpp. Pay attention to the fact that for debug we set the seed of the PRNG SBRandom to 0 in order to get reproductive results.

    // get the number of particles in the particle system
    particleIndexer = (*dynamicalModel)->getStructuralParticleIndexer();
    nParticles = particleIndexer->;size();
 
    // init
    randomNumberGenerator.seed(0);//SAMSON::getTime());
    numberOfRejects = 0;
    numberOfTrials = 0;
 
    maximumDisplacement = SBQuantity::angstrom(0.1);
    temperature = SBQuantity::kelvin(50);
    numberOfMovingParticles = 5;

Finally you can copy paste the code below inside the method updateState in order to implement the behaviour of the state updater.

    // update
    numberOfMovingParticles = (numberOfMovingParticles > nParticles)?nParticles:numberOfMovingParticles;
    beta = 1.0 / ((*SBConstant::boltzmannConstant) * temperature);
 
    // compute the current energy
    (*interactionModel)->updateInteractions(); // compute the new energy and forces based on the updated positions
    SBQuantity::energy currentEnergy = (*interactionModel)->getEnergy();
 
    SBIndexer movingParticleIndexer;
    std::vector currentPositionVector;
 
    // generate the list of moving particles
    while (movingParticleIndexer.size() != numberOfMovingParticles)
        movingParticleIndexer.push_back((unsigned int)(randomNumberGenerator.randUnsignedLong() % nParticles));
 
    // move particles incrementally
    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 displacement
        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);
    }
 
    // compute the new energy
    (*interactionModel)->updateInteractions();
    SBQuantity::energy newEnergy = (*interactionModel)->getEnergy();
 
    // accept or reject
    if (newEnergy > currentEnergy) {
        double p = randomNumberGenerator.randDouble1();
        double b = exp((beta*(currentEnergy - newEnergy)).getValue());
 
        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++;
 
    if (numberOfTrials == 50) {
        if (((double)numberOfRejects / (double)numberOfTrials) > 0.5) {
            maximumDisplacement *= 0.95;
            maximumDisplacement = (maximumDisplacement < SBQuantity::angstrom(0.01))?SBQuantity::angstrom(0.01):maximumDisplacement;
        }
        else
            maximumDisplacement *= 1.05;
 
        numberOfRejects = 0;
        numberOfTrials = 0;
    }

Finally if you modify a little bit the SE*yourname*StateUpdaterProperties.hpp and SE*yourname*StateUpdaterProperties.ui you should obtain easily something like that:

If you want to store some interesting conformations along the path you can add this piece of code in order to create a conformation and populate it with the position corresponding to a new energy record. This code is the else side of the accept / reject test:

    else {
        SBNodeIndexer atomIndexer;
        SAMSON::getActiveLayer()->getNodes(atomIndexer, SBNode::IsType(SBNode::Atom));
 
        SBConformation* conformation = new SBConformation("Pafff", atomIndexer);
        for (unsigned int i = 0; i < nParticles; ++i) // conformation->setPosition(i, (*dynamicalModel)->getPosition(i));
        conformation->create();
        SAMSON::getActiveDocument()->addChild(conformation);
    }

 

Comments are closed