1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
#include "SimTKmolmodel.h"
#include <iostream>
using namespace std;
using namespace SimTK;
void verifyAngle(Real standardAngle, Real testAngle)
{
// Normalize to (-Pi:Pi] range
while (standardAngle <= -Pi) standardAngle += 2*Pi;
while (standardAngle > Pi) standardAngle -= 2*Pi;
// Normalize to (-Pi:Pi] range
while (testAngle <= -Pi) testAngle += 2*Pi;
while (testAngle > Pi) testAngle -= 2*Pi;
cout << "Desired angle = " << standardAngle*DuMM::Rad2Deg << "\t";
cout << "Test angle = " << testAngle*DuMM::Rad2Deg << endl;
// Make sure set and get values are "pretty much" the same
Real error = testAngle - standardAngle;
error *= error;
assert(error < 1e-6); // within one milliradian
}
int main()
{
// try {
int granularity = 11; // number of steps around the torsion
Real angleStep = 2 * Pi / granularity;
// Use integer as loop index to ensure that top value is hit
for (int i = -granularity; i <= granularity; ++i)
{
// Dec 10, 2007 - works with ethane, not with cytdine
RibonucleotideResidue::Cytidylate molecule;
String torsionName("chi");
// RiboseCore molecule("1", "1", '1');
// String torsionName("delta");
// molecule.assignBiotypes();
// Ethane molecule;
// String torsionName("torsion");
Real desiredAngle = i * angleStep;
// 1) Set default angle
molecule.setDefaultDihedralAngle(torsionName.c_str(), desiredAngle);
Real measuredAngle = molecule.calcDefaultDihedralAngle(torsionName.c_str());
// molecule.writeDefaultPdb(cout);
verifyAngle(desiredAngle, measuredAngle);
// 2) Set angle in simbody model
CompoundSystem system;
SimbodyMatterSubsystem matter(system);
DuMMForceFieldSubsystem dumm(system);
dumm.loadAmber99Parameters();
dumm.loadTestMoleculeParameters();
molecule.assignBiotypes();
system.adoptCompound( molecule );
system.modelCompounds();
State state = system.realizeTopology();
system.realizeModel(state);
Real initialAngle = molecule.calcDihedralAngle(state, torsionName.c_str());
verifyAngle(desiredAngle, initialAngle);
molecule.setDihedralAngle(state, torsionName.c_str(), desiredAngle);
Real angleAfterSet = molecule.calcDihedralAngle(state, torsionName.c_str());
verifyAngle(desiredAngle, angleAfterSet);
}
cout << "PASSED" << endl;
return 0;
}
//catch (const std::exception& e)
//{
//printf("EXCEPTION THROWN: %s\n", e.what());
//}
//catch (...)
//{
//printf("UNKNOWN EXCEPTION THROWN\n");
//} return 0;
//}
|
