Programming new force fields

This tutorial will guide you to create a simple spring model force field such that a virtual spring links two pairs of atoms connected by a covalent bond.

Using the Element Generator

The first step is to use the ElementGenerator  to generate a module where you will implement your force field.

  • Lauch the ElementGenerator.
  • Select the path for your new module, press next.
  • Give a name to your module (e.g. MySpringModel), provide a description, press next.
  • Choose Interaction Model in the combo box describing the type of module to be created and press next.
  • Press Generate and press finish

Setting the description of the module

  • Compile your solution. You should now see your new module appear
  • In the header file of your  interaction model descriptor e.g. SEMySpringModelInteractoinModelDescriptor.hpp), modify the class description. Put for example My Spring Interaction Model.
 SB_CLASS_DESCRIPTION("My Spring Interaction Model");
  • In the header file  of your property model descriptor  (e.g. SEMySpringModelInteractionModelPropertiesDescriptor.hpp), modify the class description. Put for example  My Spring Interaction Model Properties.
 SB_CLASS_DESCRIPTION("My Spring InteractionModel Properties");
  • Compile and install your solution.
  • Launch Samson and load a system
  • Run a simulator: in the menu, select Simulation->Add simulator select the name of your interaction model in the list of the interaction and press OK.
  • You can see the property window if your interaction model that shows the current energy (0 eV) and the property window of the state updater you chose.
  • Close SAMSON

Initializing the interaction model

  • In the header of the file of your interaction model (e.g. SEMySpringModelInteractionModel.hpp), include SBAtom.hpp.
 #include "SBAtom.hpp";
  • In the header  file of your interaction model (e.g. SEMySpringModelInteractionModel.hpp), add as private the variables that will be used by your force field.
private:
 
	SBPointerIndexer<SBStructuralParticle> const*		        particleIndexer;
 
	SBVector<SBQuantity::length>					springLengthVector;
 
	SBVector<SBAtom*>						springAtomIVector;
	SBVector<SBAtom*>						springAtomJVector;
  • The particleIndexer is an indexer that will point to the atoms (structural particles) of your systems.
  • The springLengthVector will store the equilibrium lengths for the springs of our model.
  • The springAtomIVector and springAtomJVector will store the left and right atoms of the springs of our model.

 

  • In the .cpp file of your interaction model (e.g. SEMySpringModelInteractionModel.cpp), modify the initializeInteractions function:
  • Initialize the particle indexer
//initialize the particle indexer
	particleIndexer = (*particleSystem)->getStructuralParticleIndexer();
  • Next, get the indexer of the model bonds
//get the indexer of the model bonds 
	SBNodePredicate* nodePredicate = SAMSON::makeNodePredicate("node.type bond");
	SBNodeIndexer nodeIndexer;
	SAMSON::getActiveDocument()->getNodes(nodeIndexer, *nodePredicate);
	SAMSON::setStatusMessage(QString::number(nodeIndexer.size()), 0);
  • Next, loop on all the bonds to set springAtomIVector, springAtomJVector as well as springLengthVector. The equilibrium distance is taken as the covalent bond distances in the initial conformation of the model.
SB_FOR(SBNode* node, nodeIndexer) {
 
		SBBond* bond = static_cast<SBBond*>(node);
 
		SBAtom* atomI = bond->getLeftAtom();
		SBAtom* atomJ = bond->getRightAtom();
 
		unsigned int atomIIndex = particleIndexer->getIndex(atomI);
		unsigned int atomJIndex = particleIndexer->getIndex(atomJ);
 
		//add the atoms to the atoms Vectors
		springAtomIVector.push_back(atomI);
		springAtomJVector.push_back(atomJ);
 
		//add the equilibrium lengths to the springLength vector
		SBQuantity::length distance = ((*particleSystem)->getPosition(atomIIndex) -
			(*particleSystem)->getPosition(atomJIndex)).norm();
		springLengthVector.push_back(distance);
	}
  • Next, initialize the energy and forces to 0.
//initialize energy and forces
	*energy = SBQuantity::energy(0.0);	
	unsigned int nParticles = particleIndexer->size();
	for (unsigned int i = 0; i < nParticles; ++i)
		setForce(i, SBForce3(SBQuantity::force(0)));
  • Finally, add changed() to signal that the energy and forces of the interaction model has been changed.
 changed();

Updating the interaction model

  • In the .cpp file of your interaction model (e.g. SEMySpringModelInteractionModel.cpp), modify the updateInteractions function:
  • Initialize the energy and forces
unsigned int nParticles = particleIndexer->size();
 
		*energy = SBQuantity::energy(0.0);
		for (unsigned int i = 0; i < nParticles; ++i)
			setForce(i, SBForce3(SBQuantity::force(0)));
  • Next, for each spring, compute the current positions, and compute the force intencity as the difference of distance between the current and reference bond string.
unsigned int nSprings = springLengthVector.size();
for (unsigned int i = 0; i < nSprings; ++i) {
	SBAtom* leftAtom = springAtomIVector[i];
	SBAtom* rightAtom = springAtomJVector[i];
	unsigned int leftAtomIndex = particleIndexer->getIndex(leftAtom);
	unsigned int rightAtomIndex = particleIndexer->getIndex(rightAtom);
 
	const SBPosition3& leftAtomPosition = (*particleSystem)->getPosition(leftAtomIndex);
	const SBPosition3& rightAtomPosition = (*particleSystem)->getPosition(rightAtomIndex);
 
	//the force intensity depends on the shift respect to the equilibrium
	SBQuantity::length lengthDif = (rightAtomPosition - leftAtomPosition).norm() - springLengthVector[i];
  • In the same loop, compute the forces on the I and the J atoms and set these forces in the interaction model
	SBQuantity::forcePerLength forceFactor(100);
 
	SBForce3  force = forceFactor * lengthDif* (rightAtomPosition - leftAtomPosition).normalizedVersion();
 
	SBForce3 forceI = getForce(leftAtomIndex) + force;
	SBForce3 forceJ = getForce(rightAtomIndex) - force;
 
	setForce(leftAtomIndex, forceI);
	setForce(rightAtomIndex, forceJ);
  • Still in this loop, sum up the corresponding energy. Finally, close the loop
	*energy += 0.5 * forceFactor * lengthDif* lengthDif;
}
  • Finally, add changed() to signal that the energy and forces of the interaction model has been changed.
 changed();
  • Launch Samson and load a system.
  • Run a simulator: in the menu, select Simulation->Add simulator select the name of your interaction model in the list of the interaction and press OK.
  • Launch your simulator: press play or use the corresponding  icon   . The simulator should now work properly, updating the force and energies along the simulation.

Customizing the interaction model interface

  • We want to allow to control the force factor through the property window of the interaction model.
  • Open the interface of the interaction model property window with Qt Designer (e.g. SEMySpringModelInteractionModelProperties.ui).
  • Add  a label with text Force factor and a spin Box named spinBoxForceFactor. Set the current spin box value to 100 and its maximum value to 100000.

  • Create a slot named onSpinBoxForceFactorChanged(int) and connect it to the signal valueChanged(int) of  spinBoxForceFactor.
  • In the header of the interaction model property, declare onSpinBoxForceFactorChanged(int) as public slot.
public slots: 
 
void						onSpinBoxForceFactorValueChanged(int);
  • In the .cpp of the Interaction model property window, define the onSpinBoxForceFactorChanged(int) function.
  • Modify the .cpp and .hpp of the interaction model, to create a function that sets the forceFactor variable.
  • .hpp:
 void	setForceFactor(int _forceFactor);
 private:
 
	SBQuantity::forcePerLength									forceFactor;
  •  in .cpp :
void SEMyTestSpringModelInteractionModel::setForceFactor(int _forceFactor) {
 
	forceFactor = SBQuantity::forcePerLength(_forceFactor);
}
  • in initializeInteractions of .cpp
forceFactor = SBQuantity::forcePerLength(100);
  • And remove the local variable forceFactor in updateInteractions().
  • Finally, call the setForceFactor function in the onSpinBoxForceFactorValueChanged slot of the uinteraction model property.
void SEMyTestSpringModelInteractionModelProperties::onSpinBoxForceFactorValueChanged(int forceFactor) {
 
	interactionModel->setForceFactor(forceFactor);
}
  • Now, you can launch the simulator in SAMSON and modify the Force factor through the property window.

Comments are closed.