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
|
#include "SimTKmolmodel.h"
#include <iostream>
using namespace std;
using namespace SimTK;
int main()
{
// Start with a simple protein
// NOTE that the first residue, number 0, will actually be an acetyl end cap.
// The first real amino acid, an alanine, will be residue number 1
// Protein protein("ACDEFGHIKLMNPQRSTVWY");
// In order to extract body ids, we need to model the protein in a system
CompoundSystem system;
SimbodyMatterSubsystem matter(system);
// DecorationSubsystem decorations(system);
DuMMForceFieldSubsystem dumm(system);
// Use hard-coded force field parameters
dumm.loadAmber99Parameters();
std::vector<Protein> proteins;
{
Protein protein1("QISQAIKYLQNNIKGFIIRQRVNDEMKVNCATLLQAAYRGHSIRANVF");
// Protein protein1("ACDEFGHIKLMNPQRSTVWY");
protein1.assignBiotypes();
proteins.push_back(protein1);
Protein protein2("AAAA");
protein2.assignBiotypes();
proteins.push_back(protein2);
// adoptCompound must be done after proteins vector has stabilized
system.adoptCompound(proteins[0]);
system.adoptCompound( proteins[1], Vec3(-2, 0, 0) );
}
// Now we can model
system.modelCompounds();
const Protein& protein = proteins[0];
// Next extract the body ids and locations
int numResidues = protein.getNumResidues() - 1;
for (ResidueInfo::Index r(0); r < numResidues; ++r)
{
const ResidueInfo& residue = protein.getResidue(ResidueInfo::Index(r+1));
for (ResidueInfo::AtomIndex ra(0); ra < residue.getNumAtoms(); ++ra)
{
Compound::AtomIndex a = residue.getAtomIndex(ra);
cout << residue.getPdbResidueName();
cout << "\t";
cout << residue.getPdbResidueNumber();
cout << "\t";
cout << residue.getAtomName(ra);
cout << "\t";
cout << "Body Index = " << protein.getAtomMobilizedBodyIndex(a);
cout << "\t";
cout << protein.getAtomLocationInMobilizedBodyFrame(a);
const MobilizedBody& body = matter.getMobilizedBody(protein.getAtomMobilizedBodyIndex(a));
const MobilizedBody& parentBody = body.getParentMobilizedBody();
cout << "\t";
cout << "parent body index = " << parentBody.getMobilizedBodyIndex();
cout << endl;
}
}
}
|
