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
|
#include "Molmodel.h"
#include <iostream>
#include <exception>
#include <fstream>
using namespace SimTK;
int main() {
try {
CompoundSystem system;
SimbodyMatterSubsystem matter(system);
DecorationSubsystem decorations(system);
DuMMForceFieldSubsystem forceField(system);
// Atom classes are available, but not charged atom types for ethane
// in standard Amber force field
forceField.loadAmber99Parameters();
if (! Biotype::exists("ethane", "C"))
Biotype::defineBiotype(new Element::Carbon(), 4, "ethane", "C");
if (! Biotype::exists("ethane", "H"))
Biotype::defineBiotype(new Element::Hydrogen(), 1, "ethane", "H");
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"ethane C",
forceField.getAtomClassIndex("CT"),
-0.060 // made up
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("ethane C"),
Biotype::get("ethane", "C").getIndex() );
forceField.defineChargedAtomType(
forceField.getNextUnusedChargedAtomTypeIndex(),
"ethane H",
forceField.getAtomClassIndex("HC"),
0.020 // made up, use net neutral charge
);
forceField.setBiotypeChargedAtomType(
forceField.getChargedAtomTypeIndex("ethane H"),
Biotype::get("ethane", "H").getIndex() );
Ethane ethane1, ethane2;
ethane1.writeDefaultPdb(std::cout);
// place first ethane, units are nanometers
// skew it a little to break strict symmetry
system.adoptCompound(ethane1, Transform(Vec3(-0.5, 0, 0))
* Transform(Rotation(0.1, YAxis)) );
// place second ethane, units are nanometers
system.adoptCompound(ethane2, Vec3( 0.5, 0, 0));
// Show me a movie
Visualizer viz(system);
system.addEventReporter( new Visualizer::Reporter(viz, 0.050) );
system.modelCompounds(); // finalize multibody system
State state = system.realizeTopology();
// Simulate it.
VerletIntegrator integ(system);
TimeStepper ts(system, integ);
ts.initialize(state);
ts.stepTo(200.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;
}
}
|
