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
|
#============================================================================
# BALL - Example for a energy evaluations as it was used in Althaus
# et al. "A branch and cut algorithm for the optimal solution of the
# side-chain placement problem", 2000
#
# This example reads a PDB file and calculates a bonded energy using a force
# field and a non bonded energy (electrostatics only) by solving the Poisson-
# Boltzmann equation.
#============================================================================
#
# Lara Schneider 2010-04-22
#
import sys
from BALL import *
#issue a usage hint if called without parameters
if (len(sys.argv) != 2):
print sys.argv[0], " <PDB infile>"
exit()
#open a PDB file with the name of the first argument
file = MolFileFactory.open(sys.argv[1])
if (not file):
#if file does not exist: complain and abort
print "error opening ", argv[1], " for input."
exit()
s = System()
file.read(s)
file.close()
#normalize the names and build all bonds
print "normalizing names and building bonds..."
db = FragmentDB("")
s.apply(db.normalize_names)
s.apply(db.build_bonds)
#create an AMBER force field without non-bonded interactions
FF = AmberFF(s);
#calculate the total energy
total_energy = FF.updateEnergy();
print FF.getResults()
print " total energy using force field evaluation: ", total_energy, " kJ/mol"
print "removing non bonded energy terms ..."
FF.removeComponent("Amber NonBonded")
#calculate the internal energy (neglecting internal VdW interactions)
internal_energy = FF.updateEnergy()
print FF.getResults()
print " internal energy: ", internal_energy, " kJ/mol"
#assign atom radii
radius_processor = AssignRadiusProcessor("radii/PARSE.siz")
s.apply(radius_processor)
#calculate the electrostatic part of the solvation energy
print "calculating electrostatic energy terms with FD-Poisson-Boltzmann ..."
fdpb = FDPB(s)
fdpb.solve();
print "... using dielectric constant in medium: ", fdpb.options[FDPB.Option.SOLVENT_DC].toFloat()
solvent_energy = fdpb.getEnergy()
fdpb.options.setReal(FDPB.Option.SOLVENT_DC, 1.0)
print "... using dielectric constant in vacuum: ", fdpb.options[FDPB.Option.SOLVENT_DC].toFloat()
fdpb.setup(s)
fdpb.solve()
vacuum_energy = fdpb.getEnergy()
print "\n electrostatic solvation free energy: ", solvent_energy - vacuum_energy
print " total energy using a combination of force field and FDPB evaluation: ", internal_energy - vacuum_energy + solvent_energy, " kJ/mol"
#that's it
print "done."
|
