Programming new force fields: Springs interaction model

In this tutorial, we will create a simple spring model force field such that a virtual spring links two pairs of atoms connected by a covalent bond.

In SAMSON, a force field is presented as an interaction model.

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 Springs 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 Springs and add some description;
• add an Interaction Model class (called SESpringsInteractionModel);
• generate the SAMSON Extension.

The SAMSON Extension generator generates an interaction model class (derived from the SBMInteractionModelParticleSystem class) and a property widget class for the interaction model.

Setting the description of the Extension

Open the header file of your interaction model descriptor (SESpringsInteractionModelDescriptor.hpp) and modify the class description:

SB_CLASS_DESCRIPTION("DevTutorial: Springs");

Do the same for the property window descriptor (SESpringsInteractionModelPropertiesDescriptor.hpp). And modify the getName function of the property window (SESpringsInteractionModelProperties.cpp) to have a user-friendly description of your interaction model inside SAMSON:

QString SESpringsInteractionModelProperties::getName() const { 
 
    return "DevTutorial: Springs properties";
 
}

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 interaction models.

Initializing the interaction model

Open the header file of your interaction model (SESpringsInteractionModel.hpp) and include SBBond.hpp.

#include "SBBond.hpp"

Add the following functions and variables in the SESpringsInteractionModel class that will be used in the force field.

class SESpringsInteractionModel : public SBMInteractionModelParticleSystem {
 
public:
 
    // ...
    
    void                                            initializeBondSpringsModel();       ///< set up the bond spring model
 
private:
 
    void                                            updateBondSpringsInteractions();    ///< update the bond spring interactions
 
    SBPointerIndexer<SBAtom> const*                 particleIndexer;
 
    double                                          bondSpringStiffness;    ///< the stiffness of bond springs
 
    SBVector<SBQuantity::length>                    springLengthVector;     ///< the equilibrium length
    SBVector<SBBond*>                               springBondVector;       ///< the vector of all bonds in the selected system
 
};

Here:

• initializeBondSpringsModel is a function that initializes the bond spring interaction model;
• updateBondSpringsInteractions is a function that updates the bond spring interactions;
• particleIndexer is an indexer that will point to the atoms (structural particles) of your system;
• bondSpringStiffness regulates the stiffness of bond springs;
• springLengthVector stores the equilibrium length for the springs of your model;
• springBondVector stores the bonds of the springs of your model.

In both constructors of the interaction model (SESpringsInteractionModel.cpp) initialize the bondSpringStiffness variable to 1.0:

bondSpringStiffness = 1.0;

Modify the initializeInteractions (SESpringsInteractionModel.cpp). This function is called to set your interaction model up and initialize the energy and forces.

void SESpringsInteractionModel::initializeInteractions() {
 
    // Initialize the bond springs model
    
    initializeBondSpringsModel();
    
    // Signal that the energy and forces of the interaction model have been changed
    
    changed();
 
}

Here, we invoke the initialization function for the bond spring interaction model and emit the changed signal to let SAMSON know that the energy and forces have been changed.

Implement the initializeBondSpringsModel function as follows:

void SESpringsInteractionModel::initializeBondSpringsModel() {
 
    // Clear vectors
    springBondVector.clear();
    springLengthVector.clear();
 
    // Initialize the set of atoms: the selected atoms are stored in particleIndexer
    particleIndexer = (*particleSystem)->getAtomIndexer();
 
    // Get the all the bonds in the active document
    SBNodeIndexer nodeIndexer;
    SAMSON::getActiveDocument()->getNodes(nodeIndexer, SBNode::IsType(SBNode::Bond));
 
    // Initialize bond springs
    SB_FOR(SBNode* node, nodeIndexer) {
 
        SBPointer<SBBond> bond = static_cast<SBBond*>(node);
        if (!bond.isValid()) continue;
 
        SBPointer<SBAtom> atomI = bond->getLeftAtom();
        SBPointer<SBAtom> atomJ = bond->getRightAtom();
        if (!atomI.isValid() || !atomJ.isValid()) continue;
 
        // The spring is considered only if the two bonded atoms are selected
        // Check if the selected atoms are connected through a bond
        if (!particleIndexer->hasIndex(atomI()) || !particleIndexer->hasIndex(atomJ())) continue;
 
        // Add the bond to the bond vector
        springBondVector.push_back(bond());
 
        // Add the equilibrium length to the springLength vector
        springLengthVector.push_back(bond->getLength());
 
    }
 
    // Set energy and forces to zero
     *energy = SBQuantity::energy(0.0);
    for (unsigned int i = 0; i < particleIndexer->size(); ++i)
        setForce(i, SBForce3(SBQuantity::force(0)));
 
}

First, we clear the vectors used to store the bonds and their equilibrium lengths. Then we initialize the particle indexer (particleIndexer) based on the particle system (particleSystem) the interaction model is attached to. Then we get an indexer of all the nodes in the active document and loop through this indexer checking whether the atoms connected by the bond are present in the particle system and adding bonds and their equilibrium length to the according vectors. The equilibrium spring length is taken as the covalent bond distances in the initial conformation of the model. Then we initialize the energy and forces to zero.

Updating the interaction model

Now, we implement the update of the interaction model.

Modify the updateInteractions function (SESpringsInteractionModel.cpp). This function is called to update the energy and forces.

void SESpringsInteractionModel::updateInteractions() {
 
    // Reset energy and forces of the system
     *energy = SBQuantity::energy(0.0);
    for (unsigned int i = 0; i < particleIndexer->size(); ++i)
        setForce(i, SBForce3(SBQuantity::force(0.0)));
 
    // Apply the bond spring interaction model
    updateBondSpringsInteractions();
 
    // Signal that the energy and forces of the interaction model has been changed
    changed();
 
}

Here, we reset the energy and forces, then apply the bond spring interaction model, and emit the changed signal to let SAMSON know that the energy and forces have been changed.

Implement the updateBondSpringsInteractions function as follows:

void SESpringsInteractionModel::updateBondSpringsInteractions() {
 
    SBQuantity::kJPerMol energyBondSprings = SBQuantity::energy(0.0); // Store the bond stretch energy
 
    // Conversion for the bond spring stiffness coefficient
    SBQuantity::forcePerLength forceFactor = bondSpringStiffness * SBQuantity::zeptojoule(690.0) / SBQuantity::squareAngstrom(0.5);
 
    for (unsigned int i = 0; i < springBondVector.size(); ++i) {
 
        SBPointer<SBBond> bond = static_cast<SBBond*>(springBondVector[i]);
        if (!bond.isValid()) continue;
 
        // Get atoms connected through the bond
        SBPointer<SBAtom> atomI = springBondVector[i]->getLeftAtom();
        SBPointer<SBAtom> atomJ = springBondVector[i]->getRightAtom();
        if (!atomI.isValid() || !atomJ.isValid()) continue;
 
        unsigned int indexI = particleIndexer->getIndex(atomI());
        unsigned int indexJ = particleIndexer->getIndex(atomJ());
 
        const SBPosition3& positionI = (*particleSystem)->getPosition(indexI);
        const SBPosition3& positionJ = (*particleSystem)->getPosition(indexJ);
 
        // The force intensity depends on the shift with respect to the equilibrium length
        SBQuantity::length forceIntensity(0);
        forceIntensity = (positionJ - positionI).norm() - springLengthVector[i];
 
        // Compute the forces acting on both atoms connected through the bond
        SBForce3 force = forceFactor * forceIntensity * (positionJ - positionI).normalizedVersion();
        SBForce3 forceI = getForce(indexI);
        SBForce3 forceJ = getForce(indexJ);
 
        setForce(indexI, forceI + force);
        setForce(indexJ, forceJ - force);
 
        // Compute the energy
        energyBondSprings += 0.5 * forceFactor * forceIntensity * forceIntensity;
 
    }
 
    // Update the energy of the system
     *energy += energyBondSprings;
 
}

In the loop, for each spring:

• get positions of atoms connected via a bond,
• compute the force intencity as the difference of distance between the current and reference bond string (the equilibrium bond length),
• compute the force acting on atoms connected via the bond,
• sum up the energy.

You have implemented a simple interaction model. Build the Extension and launch SAMSON. Create or load any system, and apply a simulator to it (Simulation > Add or Ctrl / Cmd⌘ + Shift + M ) with the newly created interaction model and the Interactive modeling state updater:

Adding a simulator

You will see property windows of the interaction model and the state updater:

Property windows of the interaction model and the state updater

Run a simulator: press Space or the play icon in the toolbar:

Press play to run a simulator

The simulator should now work properly, updating the energy and forces along the simulation. Try to slightly pull an atom connected via a bond to another atom. You should see how the bond spring tries to bring another atom along with the one being pulled.

Try pulling an atom

Customizing the interaction model interface

Changing the spring stiffness parameter

We want to allow the user to control the bond spring stiffness (bondSpringStiffness) through the property window of the interaction model.

Open the SESpringsInteractionModelProperties.ui, add a label and a QDoubleSpinBox, name the spin box as doubleSpinBoxSpringStiffness, set its value to 1.0, and minimum to zero, as shown on an image below:

Add a spin box for the spring stiffness parameter

Add a slot onBondSpringStiffnessChanged(double):

Add a slot

Connect the valueChanged(double) signal of the doubleSpinBoxSpringStiffness to the onBondSpringStiffnessChanged(double) slot:

Connect a signal to a slot

Add the onBondSpringStiffnessChanged(double) slot in the SESpringsInteractionModelProperties class (SESpringsInteractionModelProperties.hpp):

public slots:
 
    void    onBondSpringStiffnessChanged(double);   ///< Change the bond spring stiffness

and implement it as follows (SESpringsInteractionModelProperties.cpp):

void SESpringsInteractionModelProperties::onBondSpringStiffnessChanged(double stiffness) {
 
    interactionModel->setBondSpringStiffness(stiffness);
 
}

Add the setBondSpringStiffness function to the SESpringsInteractionModel class (SESpringsInteractionModel.hpp):

void setBondSpringStiffness(double stiffness) { bondSpringStiffness = stiffness; }  ///< Set the bond spring stiffness parameter

Now, when you launch a simulator with this interaction model in SAMSON you will be able to modify the bond spring stiffness parameter through the property window. Try to pull an atom with different bond spring stiffness parameter:

Try pulling an atom with different bond spring stiffness parameter

If we want to save the spring stiffness parameter value for the next launch, we need to implement the saveSettings and loadSettings functions for the SESpringsInteractionModelProperties class:

void SESpringsInteractionModelProperties::saveSettings( SBGSettings *settings ) {
 
    if ( settings == 0 ) return;
 
    settings->saveValue("bondSpringStiffness", ui.doubleSpinBoxSpringStiffness->value());
 
}

 void SESpringsInteractionModelProperties::loadSettings( SBGSettings *settings ) {
 
    if ( settings == 0 ) return;
    
    double bondSpringStiffness = settings->loadDoubleValue("bondSpringStiffness", 1.0);
    ui.doubleSpinBoxSpringStiffness->setValue(bondSpringStiffness);
    onBondSpringStiffnessChanged(bondSpringStiffness);
 
}

Reinitializing the equilibrium bond distance

Let's now provide the user the possibility to reinitialize the springs equillibrium length.

For that, add a QPushButton to the SESpringsInteractionModelProperties.ui:

Add a button

Add a slot onBondSpringStiffnessChanged(double):

Add a slot

Connect the valueChanged(double) signal of the doubleSpinBoxSpringStiffness to the onBondSpringStiffnessChanged(double) slot:

Connect a signal to a slot

Add the onResetSpringsClicked slot in the SESpringsInteractionModelProperties class (SESpringsInteractionModelProperties.hpp):

public slots:
 
    void    onResetSpringsClicked();    ///< Reset the spring model to the current conformation

and implement it as follows (SESpringsInteractionModelProperties.cpp):

void SESpringsInteractionModelProperties::onResetSpringsClicked() {
 
    // Reinitialize the spring model with the current bond lengths
    
    interactionModel->initializeBondSpringsModel();
 
}

Now, when you launch the simulator with this interaction model in SAMSON you will be able to reset the spring equillibrium distance through the property window. Try to pull an atom and reset the springs equillibrium distance:

Try to pull an atom and reset the springs equillibrium distance