LOOS  v2.3.2
solvate.py
1 #!/usr/bin/env python
2 
3 import sys
4 import os
5 import random
6 import stat
7 import shutil
8 import loos
9 import PSFGen
10 import WaterBox
11 
12 config = PSFGen.ReadConfig(sys.argv[1])
13 if len(sys.argv) > 2:
14  config.directory = sys.argv[2]
15 
16 
17 # make sure the working directory exists and is writeable
18 if not os.path.exists(config.directory):
19  sys.stderr.write("Creating directory %s\n" % config.directory)
20  os.makedirs(config.directory)
21 elif not os.access(config.directory, os.R_OK|os.W_OK|os.X_OK):
22  sys.stderr.write("Don't have access to directory %s\n"
23  % config.directory)
24  sys.stderr.write("Will attempt to change permissions\n")
25  try:
26  os.chmod(config.directory,
27  stat.S_IREAD|stat.S_IWRITE|stat.S_IEXEC)
28  except OSError:
29  sys.stderr.write(" Unable to change permissions on directory\n")
30  os.exit(-1)
31 
32 
33 
34 # We need a model containing just the protein and lipid, so we'll
35 # make a psfgen script, run it, use the psf to make the AtomicGroup,
36 temporary_psfname = os.path.join(config.directory, config.psfname)
37 psfgen_script = config.generate_psf(True, False, True, True, temporary_psfname)
38 psfgen = PSFGen.PSFGen(psfgen_script, config.psfgen_binary)
39 psfgen.run()
40 system = loos.createSystem(temporary_psfname)
41 
42 segments = []
43 for segment in config.segments:
44  s = loos.selectAtoms(system, 'segname == "' + segment.segname + '"')
45  if (len(s) == 0):
46  sys.stderr.write("Selection failed assembling system: segment %s doesn't exist\n" %
47  (segment.segname)
48  )
49  sys.stderr.write("Exiting...\n")
50  sys.exit(0)
51 
52  segments.append(s)
53 
54 
55 # copy the protein coordinates into the system
56 if config.protein is not None:
57  for s in config.protein.segments:
58  current_seg = s[0].segid()
59  # Don't need to trap failed selection here, because we
60  # already know this segment exists
61  seg = loos.selectAtoms(system, 'segname == "' + current_seg + '"')
62  seg.copyMappedCoordinatesFrom(s)
63 
64 
65 
66 sys.stderr.write("Beginning water box construction\n")
67 # now add water and salt
68 water_template = loos.GCoord(1.0, 1.0, 1.0)
69 water_template *= config.water.box_size
70 water_target = loos.GCoord(config.box.x(), config.box.y(),
71  config.box.z())
72 water = WaterBox.WaterBox(config.water.coords_filename,
73  water_template, water_target,
74  config.water.segname)
75 
76 
77 # Figure out how many ions we're going to add
78 total_salt = 0
79 for salt in config.salt:
80  total_salt += salt.numres
81 total_water_and_salt = total_salt + config.water.numres
82 
83 sys.stderr.write("Water box has %d waters before superposition\n" %
84  (len(water.full_system)/3))
85 sys.stderr.write("Final target: %d waters\n" % (config.water.numres))
86 
87 
88 # Verify we have enough water. We need enough to end up with
89 # the planned number of waters, even after we remove one water molecule
90 # for each ion we add.
91 if len(water.full_system)/3 < total_water_and_salt:
92  raise ValueError, "Too few waters before superposition: %d %d" % (
93  len(water.full_system)/3, total_water_and_salt)
94 
95 # translate so that the water box is centered at the origin
96 water.full_system.centerAtOrigin()
97 water_residues = water.full_system.splitByResidue()
98 
99 
100 for w in water_residues:
101  w.reimage()
102 
103 sys.stderr.write("Beginning bump-checking water against protein\n")
104 # bump-check the water against the protein
105 # First step is selecting water and protein heavy atoms
106 # These selections can't fail, so we won't trap for zero selection
107 # size.
108 water_oxygens = loos.selectAtoms(water.full_system, 'name =~ "^O"')
109 protein_heavy = loos.selectAtoms(system, '!(name =~ "^H")')
110 
111 # find water oxygens within 1.75 Ang of any lipid heavy atom
112 clashing_oxygens = water_oxygens.within(1.75, protein_heavy, config.box)
113 sys.stderr.write("Found %d clashing waters\n" % (len(clashing_oxygens)))
114 
115 # loop over the clashing oxygens, and find which residue each is in,
116 # and remove the corresponding residues from the full water box
117 for ox in clashing_oxygens:
118  i=0
119  #found = False
120  for w in water_residues:
121  if ox in w:
122  water.full_system.remove(w)
123  break
124  i+=1
125 
126 # verify we have enough water
127 if len(water.full_system)/3 < total_water_and_salt:
128  raise ValueError, "Too few waters after superposition: %d %d" % (
129  len(water.full_system)/3, total_water_and_salt)
130 
131 sys.stderr.write("Finished bump-checking water against protein\n")
132 sys.stderr.write("Current # water molecules: %d\n" %
133  (len(water.full_system)/3))
134 sys.stderr.write("Adding salt\n")
135 
136 # regenerate the list of oxygens
137 water_oxygens = loos.selectAtoms(water.full_system, 'name =~ "^O"')
138 water_residues = water.full_system.splitByResidue()
139 
140 # now replace waters with salt
141 salts = []
142 for salt in config.salt:
143  ions = loos.AtomicGroup()
144  for i in range(salt.numres):
145  a = loos.Atom()
146  a.resname(salt.resname)
147  a.segid(salt.segname)
148  a.resid(i+1)
149 
150  # pick a water oxygen at random, replace it with salt
151  ox = random.choice(water_oxygens)
152  a.coords(ox.coords())
153  ions.append(a)
154 
155  # remove this water
156  for w in water_residues:
157  if ox in w:
158  water.full_system.remove(w)
159  water_oxygens.remove(ox)
160  break
161 
162  salts.append(ions)
163 
164 # verify we have enough water
165 num_waters = len(water.full_system)/3
166 if num_waters < config.water.numres:
167  raise ValueError, "Too few waters after exchanging salt: %d %d" % (
168  num_waters, config.water.numres)
169 
170 # if we have too many waters, need to delete the difference
171 if num_waters > config.water.numres:
172  water_residues = water.full_system.splitByResidue()
173  diff = num_waters - config.water.numres
174  sys.stderr.write("Removing %d waters to get to target size\n" % diff)
175  for i in range(diff):
176  # remove this water
177  ox = random.choice(water_oxygens)
178  for w in water_residues:
179  if ox in w:
180  water.full_system.remove(w)
181  water_oxygens.remove(ox)
182  break
183 
184 # renumber the residues
185 for i in range(len(water.full_system)):
186  res = i/3 + 1
187  water.full_system[i].resid(res)
188 
189 # Replace some of the waters with the internal waters from the protein.
190 if config.protein is not None and config.protein.has_water:
191  if config.water.numres < len(config.protein.water_seg())/3:
192  raise ValueError, "Protein has more internal waters than the total target: %d %d" % (
193  config.water.numres, len(config.protein.water_seg())/3 )
194  water.full_system.copyCoordinatesFrom(config.protein.water_seg(), 0,
195  len(config.protein.water_seg()))
196 
197 
198 sys.stderr.write("Assembling final system and writing out coordinates\n")
199 # append water and salt to the full system
200 system.append(water.full_system)
201 for s in salts:
202  system.append(s)
203 
204 system.renumber()
205 system.periodicBox(config.box)
206 
207 # Note: all of the manipulation we do mangles the bond list. Since
208 # we've got a psf anyway, we can safely remove the bond
209 # information here.
210 system.clearBonds()
211 
212 # write out a final pdb file
213 final_pdbfile = open(os.path.join(config.directory,"final.pdb"), "w")
214 final_pdb = loos.PDB.fromAtomicGroup(system)
215 final_pdbfile.write(str(final_pdb))
216 final_pdbfile.close()
217 
218 # while we're at it, write out a script to generate a full psf as
219 # well, then run psfgen for them
220 psf_script = open(os.path.join(config.directory,
221  "generate_full_psf.inp"), "w")
222 
223 psfgen_script = config.generate_psf(True, True, True, False)
224 psf_script.write(psfgen_script)
225 psf_script.close()
226 
227 
228 # If there's a protein, the above psfgen script will depend on the protein's
229 # psf file, so we should copy that into the directory as well
230 if config.protein is not None:
231  shutil.copy(config.protein.psf_file, config.directory)
232 
233 # Now, run the psfgen script from the output directory
234 sys.stderr.write("Running psfgen to generate final psf\n")
235 os.chdir(config.directory)
236 psfgen = PSFGen.PSFGen(psfgen_script, config.psfgen_binary)
237 psfgen.run()
238 
239 
240 # now, read the psf back in and check the total charge
241 full_system = loos.createSystem(config.psfname)
242 total_charge = full_system.totalCharge()
243 if (abs(total_charge) > 1e-3):
244  sys.stderr.write("\nWARNING WARNING WARNING WARNING\n")
245  sys.stderr.write("System has a net charge of %f\n" % total_charge)
246  sys.stderr.write("This will likely cause pressure and area artifacts when you try to run with Ewald.\n")
247 
Basic Atom class for handling atom properties.
Definition: Atom.hpp:50
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
static PDB fromAtomicGroup(const AtomicGroup &)
Class method for creating a PDB from an AtomicGroup.
Definition: pdb.cpp:528
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)
Definition: sfactories.cpp:119