unwrap_pdb.py

Reads PDB file with wrapped coordinates (from simulation with periodic boundary conditions), unwraps them and generates PDB with it.
USAGE:
python unwrap_pdb.py out_pbc.pdb 40

Keywords:

Categories:

  • core/data/basic/Vec3Cubic

Program source:

 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
import sys, math,copy

sys.path.append('../../../../../bin/')
from pybioshell.core.data.io import find_pdb
from pybioshell.core.data.basic import Vec3Cubic

from pybioshell.std import vector_core_data_basic_Vec3Cubic

if len(sys.argv) < 3 :
    print("""

    Reads PDB file with wrapped coordinates (from simulation with periodic boundary conditions), unwraps them and generates PDB with it.

USAGE:
    python unwrap_pdb.py out_pbc.pdb 40

    CATEGORIES: core/data/basic/Vec3Cubic
    KEYWORDS:   PDB input; PBC; SURPASS; Vec3Cubic
    IMG:        unwrapped.gif

  """)
    sys.exit()


pdb = find_pdb(sys.argv[1], "./")
n_atoms = pdb.count_atoms(0)
cutoff = float(sys.argv[2])

Vec3Cubic.set_box_len(cutoff)
xyz = vector_core_data_basic_Vec3Cubic()
for i in range(n_atoms): xyz.append(Vec3Cubic())

for i_model in range(0, pdb.count_models()) :
  xyz = vector_core_data_basic_Vec3Cubic()
  for i in range(n_atoms): xyz.append(Vec3Cubic())
  pdb.fill_structure(i_model, xyz)
  structure = pdb.create_structure(i_model)
  n_res = 0
  print("MODEL          ",i_model+1 )
  for ia in range(structure.count_chains()):
    chain = structure[ia]
    #wrapping first atom of every chain to the first box 
    if xyz[n_res ].x > cutoff: xyz[n_res ].x -= cutoff
    if xyz[n_res ].x < 0: xyz[n_res ].x += cutoff
    if xyz[n_res ].y > cutoff: xyz[n_res ].y -= cutoff
    if xyz[n_res ].y < 0: xyz[n_res ].y += cutoff
    if xyz[n_res ].z > cutoff: xyz[n_res ].z -= cutoff
    if xyz[n_res ].z < 0: xyz[n_res ].z += cutoff
    #calculating unwraped coordinates
    for ir in range(chain.count_residues()-1):
      ax = xyz[n_res+ir+1].closest_delta_x(xyz[n_res+ir])
      ay = xyz[n_res+ir+1].closest_delta_y(xyz[n_res+ir])
      az = xyz[n_res+ir+1].closest_delta_z(xyz[n_res+ir])
      xyz[n_res+ir+1].x = xyz[n_res+ir].x + ax
      xyz[n_res+ir+1].y = xyz[n_res+ir].y + ay
      xyz[n_res+ir+1].z = xyz[n_res+ir].z + az

      chain[ir+1][0].set(xyz[n_res+ir+1])
    n_res+=ir+2
    #writing to PDB file
    for ir in range(chain.count_residues()) :
      resid = chain[ir]
      print(resid[0].to_pdb_line())
  print("ENDMDL")

../_images/unwrapped.gif