simtk.org

Dear,
I am running a code for "AAA" by modifying the file:
ExampleAdenylateMobilitiesVTK.cpp. Each residue is treated as a rigid body and the torsion is allowed at the connecting joint.
I need to get the inertia properties of the each residue. When I used "getBodyMass", I get different masses for different residues which doesnt make sense (266.234, 328.199, 329.207).
In the following I attach my code.
It would be grateful to have your opinion on this problem.

#include "SimTKmolmodel.h"
#include "SimTKsimbody_aux.h"
#include <iostream>
#include <fstream>

using namespace std;
using namespace SimTK;


class PositionReporter : public PeriodicEventReporter {
public:
PositionReporter(const MultibodySystem& system, const MobilizedBody& body, Real interval) :
PeriodicEventReporter(interval), system(system), body(body) {
}
void handleEvent(const State& state) const {
system.realize(state, Stage::Position);

SpatialVec HCol[3];
SpatialMat Inert[3];
Inertia In[3];
Real mass[3];

for (int i = 0; i < 3; ++i) {
MobilizedBody body;
body = system.getMatterSubsystem().getMobilizedBody(MobilizedBodyIndex(i+1));

mass = body.getBodyMass(state);
std::cout<<mass<<std::endl;
}

char a ;
std::cin>>a;
}
private:
const MultibodySystem& system;
const MobilizedBody& body;
};

int main()
{
try {
// Initialize simbody objects
CompoundSystem system;
SimbodyMatterSubsystem matter(system);
DecorationSubsystem decoration(system);
DuMMForceFieldSubsystem forces(system);
forces.loadAmber99Parameters();


// Initialize molmodel objects
RNA rna("AAA");
rna.setCompoundBondMobility(BondMobility::Torsion);
rna.setRNABondMobility(BondMobility::Rigid,0,0);
rna.setRNABondMobility(BondMobility::Rigid,1,1);
rna.setRNABondMobility(BondMobility::Rigid,2,2);

system.adoptCompound(rna);
system.modelCompounds();

// Maintain temperature
system.updDefaultSubsystem().addEventHandler(new VelocityRescalingThermostat(
system, 293.15, 0.1));

MobilizedBody body1;
body1 = system.getMatterSubsystem().getMobilizedBody(MobilizedBodyIndex(1));
system.updDefaultSubsystem().addEventReporter(new PositionReporter(system, body1, 0.1));

system.realizeTopology();
State& state = system.updDefaultState();

// Relax the structure before dynamics run
LocalEnergyMinimizer::minimizeEnergy(system, state, 15.0);

// Prepare for molecular dynamics
VerletIntegrator integrator(system);
integrator.setAccuracy(0.001);
TimeStepper timeStepper(system, integrator);
timeStepper.initialize(state);

// Start simulation
timeStepper.stepTo(2.0);

return 0;

}
catch(const std::exception& e) {
std::cerr << "ERROR: " << e.what() << std::endl;
return 1;
}
catch(...) {
std::cerr << "ERROR: An unknown exception was raised" << std::endl;
return 1;
}

}

The mass differences make sense to me.

The first residue is lighter because it is capped at the 5' end with a hydroxyl group instead of a heavier phosphate group. The final residue is heavier because it has an additional hydrogen atom capping the 3' oxygen, where subsequent residues would otherwise be.

Thanks for your help. I didnt know that.

Regards