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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
#include "Molmodel.h"
#include <iostream>
#include <fstream>
using namespace SimTK;
// Propane is a three carbon linear alkane
// C(3)H(8), or CH3-CH2-CH3
class Propane : public Molecule
{
public:
// constructor
Propane()
{
setCompoundName("Propane");
// First atom
setBaseAtom( AliphaticCarbon("C1") );
setAtomBiotype("C1", "propane", "C1_or_C3");
convertInboardBondCenterToOutboard(); // this is the root of the parent compound
// Second atom
bondAtom( AliphaticCarbon("C2"), "C1/bond1" );
setAtomBiotype("C2", "propane", "C2");
// Third atom
bondAtom( AliphaticCarbon("C3"), "C2/bond2" );
setAtomBiotype("C3", "propane", "C1_or_C3");
// First methyl hydrogens
bondAtom( AliphaticHydrogen("H11"), "C1/bond2" );
bondAtom( AliphaticHydrogen("H12"), "C1/bond3" );
bondAtom( AliphaticHydrogen("H13"), "C1/bond4" );
setAtomBiotype("H11", "propane", "H1_or_H3");
setAtomBiotype("H12", "propane", "H1_or_H3");
setAtomBiotype("H13", "propane", "H1_or_H3");
// Second methylene hydrogens
bondAtom( AliphaticHydrogen("H21"), "C2/bond3" );
bondAtom( AliphaticHydrogen("H22"), "C2/bond4" );
setAtomBiotype("H21", "propane", "H2");
setAtomBiotype("H22", "propane", "H2");
// Third methyl hydrogens
bondAtom( AliphaticHydrogen("H31"), "C3/bond2" );
bondAtom( AliphaticHydrogen("H32"), "C3/bond3" );
bondAtom( AliphaticHydrogen("H33"), "C3/bond4" );
setAtomBiotype("H31", "propane", "H1_or_H3");
setAtomBiotype("H32", "propane", "H1_or_H3");
setAtomBiotype("H33", "propane", "H1_or_H3");
}
// create charged atom types
// ensure that charges sum to zero, unless molecule has a formal charge
// Must be called AFTER first propane is declared, so Biotypes and atom classes will be defined
static void setAmberLikeParameters(DuMMForceFieldSubsystem& forceField)
{
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"propane C1_or_C3",
forceField.getAtomClassIndex("CT"), // amber tetrahedral carbon
-0.060 // made up
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("propane C1_or_C3"),
Biotype::get("propane", "C1_or_C3").getIndex() );
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"propane C2",
forceField.getAtomClassIndex("CT"), // amber tetrahedral carbon
-0.040 // made up
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("propane C2"),
Biotype::get("propane", "C2").getIndex() );
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"propane H1_or_H3",
forceField.getAtomClassIndex("HC"), // amber tetrahedral carbon
0.020 // made up, use net neutral charge
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("propane H1_or_H3"),
Biotype::get("propane", "H1_or_H3").getIndex() );
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"propane H2",
forceField.getAtomClassIndex("HC"), // amber tetrahedral carbon
0.020 // made up, use net neutral charge
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("propane H2"),
Biotype::get("propane", "H2").getIndex() );
}
};
int main() {
try {
CompoundSystem system;
SimbodyMatterSubsystem matter(system);
DecorationSubsystem decorations(system);
DuMMForceFieldSubsystem forceField(system);
// Atom classes are available, but not charged atom types for propane
// in standard Amber force field
forceField.loadAmber99Parameters();
Propane propane1, propane2;
Propane::setAmberLikeParameters(forceField);
// place first propane, units are nanometers
// skew it a little to break strict symmetry
system.adoptCompound(propane1, Transform(Vec3(-0.5, 0, 0))
* Transform(Rotation(0.1, YAxis)) );
// place second propane, units are nanometers
system.adoptCompound(propane2, Vec3( 0.5, 0, 0));
// Show me a movie
Visualizer viz(system);
system.addEventReporter( new Visualizer::Reporter(viz, 0.100) );
system.modelCompounds(); // finalize multibody system
State state = system.realizeTopology();
VerletIntegrator integ(system);
TimeStepper ts(system, integ);
ts.initialize(state);
ts.stepTo(100.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;
}
}
|
