# Python ZeroDivisionError

**URL:** <https://community.unix.com/t/python-zerodivisionerror/303403>\
**Category:** Programming\
**Created:** [January 25, 2012, 5:09pm UTC](https://community.unix.com/t/python-zerodivisionerror/303403 "2012-01-25T17:09:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![chrisjorg](https://community.unix.com/letter_avatar/chrisjorg/32/5_5575768a8748004e209b776fc1b2916d.png) [@chrisjorg](https://community.unix.com/u/chrisjorg)\
**Post date:** [January 25, 2012, 5:09pm UTC](https://community.unix.com/t/python-zerodivisionerror/303403/1 "2012-01-25T17:09:14Z")

</div>

On my two data files c7ax\_FINAL.pdb and c7eq\_FINAL.pdb I wish to perform a root-mean-square-deviation calculation:

My two files are:

```nohighlight
< ATOM 1 C ALA A 1 3.953 7.886 1.035 0.00 1.00 ALA   
< ATOM 2 O ALA A 1 3.516 8.953 0.733 0.00 0.00 ALA   
< ATOM 3 NT ALA A 1 5.207 7.558 0.961 0.00 1.00 ALA   
< ATOM 4 HNT ALA A 1 5.464 6.607 1.102 0.00 0.00 ALA   
< ATOM 5 CAT ALA A 1 6.224 8.463 0.408 0.00 0.00 ALA   
< ATOM 6 HT1 ALA A 1 6.081 8.684 -0.723 0.00 0.00 ALA   
< ATOM 7 HT2 ALA A 1 7.281 8.023 0.476 0.00 0.00 ALA   
< ATOM 8 HT3 ALA A 1 6.242 9.421 0.971 0.00 0.00 ALA   
< ATOM 9 CAY ALA A 1 3.829 3.665 -0.589 0.00 1.00 ALA   
< ATOM 10 HY1 ALA A 1 4.149 4.117 -1.530 0.00 0.00 ALA   
< ATOM 11 HY2 ALA A 1 4.389 2.789 -0.273 0.00 0.00 ALA   
< ATOM 12 HY3 ALA A 1 2.715 3.393 -0.645 0.00 0.00 ALA   
< ATOM 13 CY ALA A 1 3.917 4.701 0.520 0.00 0.00 ALA   
< ATOM 14 OY ALA A 1 4.828 4.638 1.301 0.00 1.00 ALA   
< ATOM 15 N ALA A 1 2.980 5.688 0.542 0.00 0.00 ALA   
< ATOM 16 HN ALA A 1 2.165 5.607 0.046 0.00 1.00 ALA   
< ATOM 17 CA ALA A 1 3.006 6.799 1.485 0.00 1.00 ALA   
< ATOM 18 HA ALA A 1 2.073 7.253 1.419 0.00 0.00 ALA   
< ATOM 19 CB ALA A 1 3.241 6.464 2.973 0.00 0.00 ALA   
< ATOM 20 HB1 ALA A 1 3.464 7.440 3.349 0.00 0.00 ALA   
< ATOM 21 HB2 ALA A 1 4.125 5.884 3.182 0.00 0.00 ALA   
< ATOM 22 HB3 ALA A 1 2.252 6.095 3.299 0.00 0.00 ALA   
---
> ATOM 1 C ALA A 1 3.976 7.842 1.236 0.00 1.00 ALA   
> ATOM 2 O ALA A 1 4.013 8.378 2.325 0.00 0.00 ALA   
> ATOM 3 NT ALA A 1 5.027 7.787 0.418 0.00 1.00 ALA   
> ATOM 4 HNT ALA A 1 4.956 7.219 -0.411 0.00 0.00 ALA   
> ATOM 5 CAT ALA A 1 6.221 8.583 0.554 0.00 0.00 ALA   
> ATOM 6 HT1 ALA A 1 6.061 9.618 0.811 0.00 0.00 ALA   
> ATOM 7 HT2 ALA A 1 6.891 8.420 -0.307 0.00 0.00 ALA   
> ATOM 8 HT3 ALA A 1 6.814 8.138 1.406 0.00 0.00 ALA   
> ATOM 9 CAY ALA A 1 3.787 3.479 -0.092 0.00 1.00 ALA   
> ATOM 10 HY1 ALA A 1 4.869 3.206 0.041 0.00 0.00 ALA   
> ATOM 11 HY2 ALA A 1 3.160 3.112 0.762 0.00 0.00 ALA   
> ATOM 12 HY3 ALA A 1 3.421 2.957 -1.005 0.00 0.00 ALA   
> ATOM 13 CY ALA A 1 3.693 4.933 -0.222 0.00 0.00 ALA   
> ATOM 14 OY ALA A 1 4.325 5.508 -1.108 0.00 1.00 ALA   
> ATOM 15 N ALA A 1 2.926 5.567 0.715 0.00 0.00 ALA   
> ATOM 16 HN ALA A 1 2.642 5.047 1.543 0.00 1.00 ALA   
> ATOM 17 CA ALA A 1 2.783 7.045 0.789 0.00 1.00 ALA   
> ATOM 18 HA ALA A 1 2.453 7.408 -0.216 0.00 0.00 ALA   
> ATOM 19 CB ALA A 1 1.593 7.521 1.595 0.00 0.00 ALA   
> ATOM 20 HB1 ALA A 1 1.687 7.108 2.610 0.00 0.00 ALA   
> ATOM 21 HB2 ALA A 1 0.653 7.196 1.126 0.00 0.00 ALA   
> ATOM 22 HB3 ALA A 1 1.490 8.595 1.523 0.00 0.00 ALA  

```

The run command is:  
match.py -nr c7ax\_FINAL.pdb c7eq\_FINAL.pdb "[('A:1', 'A:1')]"

I get the following error:

```nohighlight
Aligning CA atoms of residues: A:1
No rotations
Traceback (most recent call last):
  File "Documents/RMSD_two_PDB/match.py", line 221, in <module>
    rmsd = get_raw_rmsd(pdb1, pdb2, segments1, segments2, ['CA'])
  File "Documents/RMSD_two_PDB/match.py", line 118, in get_raw_rmsd
    return sum_rmsd(atoms1, atoms2)
  File "Documents/RMSD_two_PDB/match.py", line 110, in sum_rmsd
    return math.sqrt(sum_squared/float(len(atoms1)))
ZeroDivisionError: float division

```

The program is

```nohighlight
#!/usr/bin/env python

"""
Set of routines to calculate the RMSD between two molecular structures.
The module can be run from the command line using PDB files as input.
"""

import math
import numpy
import vector3d, util, molecule, polymer

def rmsd(crds1, crds2):
  """Returns RMSD between 2 sets of [nx3] numpy array"""

  assert(crds1.shape[1] == 3)
  assert(crds1.shape == crds2.shape)

  n_vec = numpy.shape(crds1)[0]
  correlation_matrix = numpy.dot(numpy.transpose(crds1), crds2)
  v, s, w = numpy.linalg.svd(correlation_matrix)
  is_reflection = (numpy.linalg.det(v) * numpy.linalg.det(w)) < 0.0
  if is_reflection:
    s[-1] = - s[-1]
  E0 = sum(sum(crds1 * crds1)) + \
       sum(sum(crds2 * crds2))
  rmsd_sq = (E0 - 2.0*sum(s)) / float(n_vec)
  rmsd_sq = max([rmsd_sq, 0.0])
  return numpy.sqrt(rmsd_sq)

def optimal_superposition(crds1, crds2):
  """Returns best-fit rotation matrix as [3x3] numpy matrix"""
  assert(crds1.shape[1] == 3)
  assert(crds1.shape == crds2.shape)
  correlation_matrix = numpy.dot(numpy.transpose(crds1), crds2)
  v, s, w = numpy.linalg.svd(correlation_matrix)
  is_reflection = (numpy.linalg.det(v) * numpy.linalg.det(w)) < 0.0
  if is_reflection:
    v[-1,:] = -v[-1,:]
  return numpy.dot(v, w)

def get_i_residue(residues, tag):

  def get_tag(residue):
    tag = ""
    if residue.chain_id != " " and residue.chain_id != "":
      tag += residue.chain_id + ":"
    tag += str(residue.num)
    if residue.insert:
      tag += residue.insert
    return tag
  # clean up tag
  tag = tag.strip()
  if tag[0] == ":":
    tag = tag[1:]
  if not tag[0].isdigit() and tag[1].isdigit():
    tag = tag[0] + ":" + tag[1:]

  for i, residue in enumerate(residues):
    if tag.lower() == get_tag(residue).lower():
      return i
  raise "Can't find residue", tag

def get_superposable_atoms(polymer, segments,
           atom_types=['CA', 'N', 'C', 'CB']):
  result = []
  allowed_i = []
  residues = polymer.residues()
  for res_num_i, res_num_j in segments:
    i = get_i_residue(residues, str(res_num_i))
    j = get_i_residue(residues, str(res_num_j))
    allowed_i.extend(range(i,j))
  for i, residue in enumerate(residues):
    if i in allowed_i:
      result.extend([a for a in residue.atoms()
                     if a.type in atom_types])
  return result

def get_crds(atoms):
  crds = numpy.zeros((len(atoms), 3), float)
  for i, a in enumerate(atoms):
    crds[i,0] = a.pos.x
    crds[i,1] = a.pos.y
    crds[i,2] = a.pos.z
  return crds

def calculate_superposition_matrix(atoms1, atoms2):

  def convert_to_matrix3d(numpy_matrix3d):
    result = vector3d.Matrix3d()
    for i in range(3):
      for j in range(3):
        result.setElem(i, j, numpy_rotation[j, i])
    return result

  numpy_rotation = optimal_superposition(get_crds(atoms1), get_crds(atoms2))
  return convert_to_matrix3d(numpy_rotation)
def sum_rmsd(atoms1, atoms2):
  sum_squared = 0.0
  for atom1, atom2 in zip(atoms1, atoms2):
    sum_squared += vector3d.pos_distance(atom1.pos, atom2.pos)**2
  return math.sqrt(sum_squared/float(len(atoms1)))

def get_raw_rmsd(pdb1, pdb2, segments1, segments2, atom_types):
  polymer1 = polymer.Polymer(pdb1)
  polymer2 = polymer.Polymer(pdb2)
  atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
  atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)
  return sum_rmsd(atoms1, atoms2)

def get_best_alignment(pdb1, pdb2, segments1, segments2, atom_types):
  """Returns rmsd and filename of transformed pdb2."""
  polymer1 = polymer.Polymer(pdb1)
  atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
  polymer2 = polymer.Polymer(pdb2)
  atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)

  center1 = molecule.get_center(atoms1)
  polymer1.transform(vector3d.Translation(-center1))
  polymer2.transform(vector3d.Translation(-molecule.get_center(atoms2)))
  polymer2.transform(calculate_superposition_matrix(atoms1, atoms2))

  rmsd = sum_rmsd(atoms1, atoms2)

  temp_pdb2 = util.fname_variant(pdb2)
  polymer2.transform(vector3d.Translation(center1))
  polymer2.write_pdb(temp_pdb2)

  return rmsd, temp_pdb2

def get_rmsd(pdb1, pdb2, segments1, segments2, atom_types):
  polymer1 = polymer.Polymer(pdb1)
  atoms1 = get_superposable_atoms(polymer1, segments1, atom_types)
  polymer2 = polymer.Polymer(pdb2)
  atoms2 = get_superposable_atoms(polymer2, segments2, atom_types)

  center1 = molecule.get_center(atoms1)
  polymer1.transform(vector3d.Translation(-center1))
  polymer2.transform(vector3d.Translation(-molecule.get_center(atoms2)))

  crds1 = get_crds(atoms1)
  crds2 = get_crds(atoms2)
  return rmsd(crds1, crds2)

def segments_str(segments):
  residues = []
  for i, j in segments:
    if i == j:
   residues.append(str(i))
    else:
      residues.append("%s-%s" % (i,j))
  return ', '.join(residues)

```

Any ideas? It seems that float(len(atoms1))=0, i.e. the length of the string?

---

<div class="post-metadata">

**Author:** ![dhzdh](https://community.unix.com/letter_avatar/dhzdh/32/5_5575768a8748004e209b776fc1b2916d.png) [@dhzdh](https://community.unix.com/u/dhzdh)\
**Post date:** [January 26, 2012, 3:50am UTC](https://community.unix.com/t/python-zerodivisionerror/303403/2 "2012-01-26T03:50:48Z")

</div>

Looking at the program, it appeared to me that atoms1 is a list. len(atoms1) is the number of items in the list. From skimming the program it appears to me that it is possible under certain circumstances (I cannot tell which circumstances since I don't have the necessary background on the algorithm) that atoms1 will be an empty list, as can be seen in the function get\_superposable\_atoms. I therefore suggest to add a check against such a contingency in the code, thus refraining from an exception. .