LOOS  v2.3.2
OptimalMembraneGenerator.py
1 #!/usr/bin/env python
2 
3 import random
4 import math
5 
6 import loos
7 import LipidLibrary
8 import PSFGen
9 import NAMD
10 import WaterBox
11 
12 if __name__ == '__main__':
13  import sys
14  import shutil
15  import os
16  import stat
17 
18 
19  config = PSFGen.ReadConfig(sys.argv[1])
20  if len(sys.argv) > 2:
21  config.directory = sys.argv[2]
22 
23 
24  box_scaling = 3.0
25  box = box_scaling*config.box
26 
27  huge_z = 1000.0 # ridiculously large z dimension used initially,
28  # chosen to be effectively infinite
29 
30  box.z(huge_z)
31  scaled_box = box.x()
32  half_box = 0.5*scaled_box
33 
34 
35  overlap_dist = 3.0 # the distance that is considered a "contact"
36  overlap_threshold = 3 # the number of contacts required for a molecule
37  # to be rejected
38 
39  # set how big the "water" region of the box will be
40  water_width = config.box.z()
41 
42  # Don't need a blank spot in the middle if there are no lipids
43  if len(config.segments) > 0:
44  water_width -= 29.0 # TODO: need a smarter algorithm
45  # for this, probably using
46  # the phos values
47 
48  # set the number of minimization and dynamics NAMD will do on each cycle
49  number_of_steps = 100
50 
51 
52  # make sure the working directory exists and is writeable
53  if not os.path.exists(config.directory):
54  sys.stderr.write("Creating directory %s\n" % config.directory)
55  os.makedirs(config.directory)
56  elif not os.access(config.directory, os.R_OK|os.W_OK|os.X_OK):
57  sys.stderr.write("Don't have access to directory %s\n"
58  % config.directory)
59  sys.stderr.write("Will attempt to change permissions\n")
60  try:
61  os.chmod(config.directory,
62  stat.S_IREAD|stat.S_IWRITE|stat.S_IEXEC)
63  except OSError:
64  sys.stderr.write(" Unable to change permissions on directory\n")
65  os.exit(-1)
66 
67  # We need a model containing just the protein and lipid, so we'll
68  # make a psfgen script, run it, use the psf to make the AtomicGroup,
69  temporary_psfname = os.path.join(config.directory, config.psfname)
70  psfgen_script = config.generate_psf(True, False, True, True, temporary_psfname)
71  psfgen = PSFGen.PSFGen(psfgen_script, config.psfgen_binary)
72  psfgen.run()
73  system = loos.createSystem(temporary_psfname)
74 
75  # If the "protein" is actually a bunch of independent molecules (e.g. a bunch of peptides),
76  # we'll want to scale them in x & y to match the expanded box.
77  if config.protein is not None and config.protein.scale:
78  protein_molecules = config.protein.model.splitByMolecule()
79  for m in protein_molecules:
80  centroid = m.centroid()
81  scaled = box_scaling*centroid
82  diff = scaled - centroid
83  diff.z(0.0) # we only want to translate in the xy plane
84  m.translate(diff)
85 
86  molecules = system.splitByMolecule()
87 
88  segments = []
89  for segment in config.segments:
90  s = loos.selectAtoms(system, 'segname == "' + segment.segname + '"')
91  if (len(s) == 0):
92  sys.stderr.write("Selection failed assembling system: segment %s doesn't exist\n" %
93  (segment.segname)
94  )
95  sys.stderr.write("Exiting...\n")
96  sys.exit(0)
97 
98  segments.append(s)
99 
100  x_axis = loos.GCoord(1,0,0)
101  z_axis = loos.GCoord(0,0,1)
102 
103  for j in range(len(segments)):
104  seg_ag = segments[j]
105  seg_ag_arr = seg_ag.splitByMolecule()
106  seg_conf = config.segments[j]
107  i = 0
108 
109  while i < seg_conf.numres:
110  lipid = seg_conf.library.pick_structure()
111  phos = loos.selectAtoms(lipid, 'name == "' + seg_conf.phos_atom +
112  '"')
113  if (len(phos) == 0):
114  sys.stderr.write("Selection failed: \"phos\" atom %s doesn't exist in lipid %s\n"
115  % (seg_conf.phos_atom, seg_conf.resname))
116  sys.stderr.write("Exiting...\n")
117  sys.exit(0)
118 
119 
120  # put the molecule at the origin
121  lipid.centerAtOrigin()
122 
123  # perform random rotation about z axis
124  lipid.rotate(z_axis, random.uniform(0.,360.))
125 
126  if seg_conf.placement < 0:
127  lipid.rotate(x_axis, 180)
128 
129  # put the phosphate back at the origin
130  lipid.translate(-phos[0].coords())
131 
132  # generate the new coordinates
133  x = (random.random() * scaled_box) - half_box
134  y = (random.random() * scaled_box) - half_box
135 
136  # set the z location, with +/- 0.5 random noise
137  z = seg_conf.phos_height + (random.random() - 0.5)
138 
139  if seg_conf.placement < 0:
140  z = -z
141 
142  vec = loos.GCoord(x,y,z)
143  lipid.translate(vec)
144 
145 
146  # do a bump-check against previously placed lipids
147  # If at any time we realize we've got too many clashes,
148  # we'll break out of the loop early, and try to place
149  # this lipid again by bouncing to the top of the while
150  # loop without incrementing the index i
151 
152  # first, bump check against the protein, if any
153  num_overlap = 0
154  if config.protein is not None:
155  for seg in config.protein.segments:
156  overlap = lipid.within(overlap_dist, seg, box)
157  num_overlap += overlap.size()
158  if num_overlap > overlap_threshold: break
159  if num_overlap > overlap_threshold: continue
160 
161  # second, check this lipid against all previous segments
162  for k in range(j):
163  overlap = lipid.within(overlap_dist, segments[k], box)
164  num_overlap += overlap.size()
165  if num_overlap > overlap_threshold: break
166  if num_overlap > overlap_threshold: continue
167 
168  # now, check earlier lipids in this segment
169  for k in range(i):
170  overlap = lipid.within(overlap_dist, seg_ag_arr[k], box)
171  num_overlap += overlap.size()
172  if num_overlap > overlap_threshold: break
173  if num_overlap > overlap_threshold: continue
174 
175  # copy the coordinates back into the real object
176  #seg_ag_arr[i].copyCoordinatesFrom(lipid, 0, lipid.size())
177  seg_ag_arr[i].copyMappedCoordinatesFrom(lipid)
178 
179  i += 1
180 
181  # copy the protein coordinates into the system
182  if config.protein is not None:
183  for s in config.protein.segments:
184  current_seg = s[0].segid()
185  # Don't need to trap failed selection here, because we
186  # already know this segment exists
187  seg = loos.selectAtoms(system, 'segname == "' + current_seg + '"')
188  seg.copyMappedCoordinatesFrom(s)
189 
190  # TODO: decide if we want to put this back. Assume
191  # the "protein" is already centered"
192  #system.centerAtOrigin()
193  pdb = loos.PDB.fromAtomicGroup(system)
194  pdb_out = open(os.path.join(config.directory,"lipid_only.pdb"), "w")
195  pdb_out.write(str(pdb))
196  pdb_out.close()
197 
198  # Now start scaling the box down with successive minimizations
199  num_iter = 100
200  delta = (scaled_box - config.box.x())/num_iter
201 
202  # only need to do the scaling and minimization if we've got lipids
203  if len(config.segments) > 0:
204 
205  next_pdb = "lipid_only.pdb"
206  dim = scaled_box
207  for i in range(num_iter+1):
208  old_dim = dim
209  dim = scaled_box - (i*delta)
210 
211  # scale the centers of mass of individual molecules
212  # TODO: this will have to change to handle "fixed" molecules
213  scale_factor = dim/old_dim
214  for j in range(len(molecules)):
215  centroid = molecules[j].centroid()
216  centroid.z(0.0) # only scale in plane of membrane
217  new_centroid = centroid * scale_factor
218  diff = new_centroid - centroid
219  molecules[j].translate(diff)
220 
221  # output the new coordinates so we can minimize with NAMD
222  pdb_out = open(os.path.join(config.directory,next_pdb), "w")
223  pdb_out.write(str(pdb))
224  pdb_out.close()
225 
226 
227 
228 
229  # TODO: long term, I'd like to change this so that instead of
230  # running a new instance of NAMD, we're reusing the same
231  # process, and adding new commands
232  current_box = loos.GCoord(dim, dim, huge_z)
233  core_name = "lipid_shrink_" + str(i)
234  full_core = os.path.join(config.directory, core_name)
235 
236  # copy the psf file to the working directory
237  local_psfname = os.path.basename(config.psfname)
238  #shutil.copyfile(config.psfname,
239  # os.path.join(config.directory, local_psfname))
240 
241  namd = NAMD.NAMD(local_psfname, next_pdb,
242  core_name, config.parameters,
243  current_box, config.namd_binary)
244  namd_inputfilename = full_core + ".inp"
245  namd_outputfilename = full_core + ".out"
246  namd.write_inputfile(namd_inputfilename, number_of_steps)
247 
248  namd.write_restraintfile(config.directory, system)
249  namd.run_namd(namd_inputfilename, namd_outputfilename)
250 
251  # copy the final coordinates to a pdb file
252  shutil.copyfile(full_core + ".coor",
253  os.path.join(config.directory,next_pdb))
254 
255  # read in the new structure
256  new_system = loos.createSystem(os.path.join(config.directory,next_pdb))
257 
258  # copy coordinates back into old AtomicGroup
259  system.copyCoordinatesFrom(new_system, 0, system.size())
260 
261  # TODO: this will have to change to handle "fixed" molecules
262  system.centerAtOrigin()
263  sys.stderr.write("lipid only minimization cycle: %d\n" % i)
264 
265  pdb_out = open(os.path.join(config.directory,"lipid_min.pdb"), "w")
266  pdb_out.write(str(pdb))
267  pdb_out.close()
268 
269  sys.stderr.write("Beginning water box construction\n")
270  # now add water and salt
271  water_template = loos.GCoord(1.0, 1.0, 1.0)
272  water_template *= config.water.box_size
273  water_target = loos.GCoord(config.box.x(), config.box.y(),
274  water_width)
275  water = WaterBox.WaterBox(config.water.coords_filename,
276  water_template, water_target,
277  config.water.segname)
278 
279 
280  # Figure out how many ions we're going to add
281  total_salt = 0
282  for salt in config.salt:
283  total_salt += salt.numres
284  total_water_and_salt = total_salt + config.water.numres
285 
286  sys.stderr.write("Water box has %d waters before superposition\n" %
287  (len(water.full_system)/3))
288  sys.stderr.write("Final target: %d waters\n" % (config.water.numres))
289 
290 
291  # Verify we have enough water. We need enough to end up with
292  # the planned number of waters, even after we remove one water molecule
293  # for each ion we add.
294  if len(water.full_system)/3 < total_water_and_salt:
295  raise ValueError, "Too few waters before superposition: %d %d" % (
296  len(water.full_system)/3, total_water_and_salt)
297 
298  # translate so that the water box is centered on the periodic boundary,
299  # then set its periodic box to the full system box size and reimage
300  water.full_system.periodicBox(config.box)
301  trans = loos.GCoord(0.0, 0.0, 0.5*config.box.z())
302  water.full_system.translate(trans)
303  water_residues = water.full_system.splitByResidue()
304 
305 
306  for w in water_residues:
307  w.reimage()
308 
309  sys.stderr.write("Beginning bump-checking water against lipid\n")
310  # bump-check the water against the lipids
311  # First step is selecting water and lipid heavy atoms
312  # These selections can't fail, so we won't trap for zero selection
313  # size.
314  water_oxygens = loos.selectAtoms(water.full_system, 'name =~ "^O"')
315  lipid_heavy = loos.selectAtoms(system, '!(name =~ "^H")')
316 
317  # find water oxygens within 1.75 Ang of any lipid heavy atom
318  clashing_oxygens = water_oxygens.within(1.75, lipid_heavy, config.box)
319  sys.stderr.write("Found %d clashing waters\n" % (len(clashing_oxygens)))
320 
321  # loop over the clashing oxygens, and find which residue each is in,
322  # and remove the corresponding residues from the full water box
323  for ox in clashing_oxygens:
324  i=0
325  #found = False
326  for w in water_residues:
327  if ox in w:
328  water.full_system.remove(w)
329  break
330  i+=1
331 
332  # verify we have enough water
333  if len(water.full_system)/3 < total_water_and_salt:
334  raise ValueError, "Too few waters after superposition: %d %d" % (
335  len(water.full_system)/3, total_water_and_salt)
336 
337  sys.stderr.write("Finished bump-checking water against lipid\n")
338  sys.stderr.write("Current # water molecules: %d\n" %
339  (len(water.full_system)/3))
340  sys.stderr.write("Adding salt\n")
341 
342  # regenerate the list of oxygens
343  water_oxygens = loos.selectAtoms(water.full_system, 'name =~ "^O"')
344  water_residues = water.full_system.splitByResidue()
345 
346  # now replace waters with salt
347  salts = []
348  for salt in config.salt:
349  ions = loos.AtomicGroup()
350  for i in range(salt.numres):
351  a = loos.Atom()
352  a.resname(salt.resname)
353  a.segid(salt.segname)
354  a.resid(i+1)
355 
356  # pick a water oxygen at random, replace it with salt
357  ox = random.choice(water_oxygens)
358  a.coords(ox.coords())
359  ions.append(a)
360 
361  # remove this water
362  for w in water_residues:
363  if ox in w:
364  water.full_system.remove(w)
365  water_oxygens.remove(ox)
366  break
367 
368  salts.append(ions)
369 
370  # verify we have enough water
371  num_waters = len(water.full_system)/3
372  if num_waters < config.water.numres:
373  raise ValueError, "Too few waters after exchanging salt: %d %d" % (
374  num_waters, config.water.numres)
375 
376  # if we have too many waters, need to delete the difference
377  if num_waters > config.water.numres:
378  water_residues = water.full_system.splitByResidue()
379  diff = num_waters - config.water.numres
380  sys.stderr.write("Removing %d waters to get to target size\n" % diff)
381  for i in range(diff):
382  # remove this water
383  ox = random.choice(water_oxygens)
384  for w in water_residues:
385  if ox in w:
386  water.full_system.remove(w)
387  water_oxygens.remove(ox)
388  break
389 
390  # renumber the residues
391  for i in range(len(water.full_system)):
392  res = i/3 + 1
393  water.full_system[i].resid(res)
394 
395  # Replace some of the waters with the internal waters from the protein.
396  if config.protein is not None and config.protein.has_water:
397  if config.water.numres < len(config.protein.water_seg())/3:
398  raise ValueError, "Protein has more internal waters than the total target: %d %d" % (
399  config.water.numres, len(config.protein.water_seg())/3 )
400  water.full_system.copyCoordinatesFrom(config.protein.water_seg(), 0,
401  len(config.protein.water_seg()))
402 
403 
404  sys.stderr.write("Assembling final system and writing out coordinates\n")
405  # append water and salt to the full system
406  system.append(water.full_system)
407  for s in salts:
408  system.append(s)
409 
410  system.renumber()
411  system.periodicBox(config.box)
412 
413  # Note: all of the manipulation we do mangles the bond list. Since
414  # we've got a psf anyway, we can safely remove the bond
415  # information here.
416  system.clearBonds()
417 
418  # write out a final pdb file
419  final_pdbfile = open(os.path.join(config.directory,"final.pdb"), "w")
420  final_pdb = loos.PDB.fromAtomicGroup(system)
421  final_pdbfile.write(str(final_pdb))
422  final_pdbfile.close()
423 
424  # while we're at it, write out a script to generate a full psf as
425  # well, then run psfgen for them
426  psf_script = open(os.path.join(config.directory,
427  "generate_full_psf.inp"), "w")
428 
429  psfgen_script = config.generate_psf(True, True, True, False)
430  psf_script.write(psfgen_script)
431  psf_script.close()
432 
433 
434  # If there's a protein, the above psfgen script will depend on the protein's
435  # psf file, so we should copy that into the directory as well
436  if config.protein is not None:
437  shutil.copy(config.protein.psf_file, config.directory)
438 
439  # Now, run the psfgen script from the output directory
440  sys.stderr.write("Running psfgen to generate final psf\n")
441  os.chdir(config.directory)
442  psfgen = PSFGen.PSFGen(psfgen_script, config.psfgen_binary)
443  psfgen.run()
444 
445 
446  # now, read the psf back in and check the total charge
447  full_system = loos.createSystem(config.psfname)
448  total_charge = full_system.totalCharge()
449  if (abs(total_charge) > 1e-3):
450  sys.stderr.write("\nWARNING WARNING WARNING WARNING\n")
451  sys.stderr.write("System has a net charge of %f\n" % total_charge)
452  sys.stderr.write("This will likely cause pressure and area artifacts when you try to run with Ewald.\n")
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