diff --git a/modules/mol/alg/src/pdbize.cc b/modules/mol/alg/src/pdbize.cc
index d6371e639c5519f21e2a4bab4de9a254584e19aa..e99366831c01bf2f4c8e4b422a4d617c1722d654 100644
--- a/modules/mol/alg/src/pdbize.cc
+++ b/modules/mol/alg/src/pdbize.cc
@@ -106,13 +106,20 @@ void PDBize::Add(EntityView asu, const geom::Mat4List& transforms,
       }
       if (chain.GetType() == CHAINTYPE_WATER) {
         if (!water_chain_) {
+          last_water_rnum_ = ResNum(0, 'Z');
           water_chain_ = edi.InsertChain(WATER_CHAIN_NAME);
           edi.SetChainDescription(water_chain_, chain.GetDescription());
           edi.SetChainType(water_chain_, chain.GetType());
         }
         for (ResidueViewList::const_iterator k = chain.GetResidueList().begin(),
              e3 = chain.GetResidueList().end(); k != e3; ++k) {
-          ResidueHandle new_res = edi.AppendResidue(water_chain_, k->GetName());
+          if (last_water_rnum_.GetInsCode() == 'Z') {
+            last_water_rnum_ = ResNum(last_water_rnum_.GetNum()+1, 'A');
+          } else {
+            last_water_rnum_.SetInsCode(last_water_rnum_.GetInsCode()+1);
+          }
+          ResidueHandle new_res = edi.AppendResidue(water_chain_, k->GetName(), 
+                                                    last_water_rnum_);
           transfer_residue_properties(*k, new_res);
           new_res.SetStringProp("type", StringFromChainType(chain.GetType()));
           new_res.SetStringProp("description", chain.GetDescription());
diff --git a/modules/mol/alg/src/pdbize.hh b/modules/mol/alg/src/pdbize.hh
index f3afec693524f12e3c691db8fe7b65ba318dc050..b344cfca415f809d39a5164396f7ae844a444945 100644
--- a/modules/mol/alg/src/pdbize.hh
+++ b/modules/mol/alg/src/pdbize.hh
@@ -52,6 +52,7 @@ private:
   const char*                           curr_chain_name_;
   bool                                  needs_adjustment_;
   int                                   last_rnum_;
+  ResNum                                last_water_rnum_;
   std::map<ResidueHandle,ResidueHandle> dst_to_src_map_;
 };
 
diff --git a/modules/mol/alg/tests/CMakeLists.txt b/modules/mol/alg/tests/CMakeLists.txt
index 758bbbe10df71687391cac11a0425a253b94dc27..8258cb3d1398734c19b1e1f73051444dc87711ab 100644
--- a/modules/mol/alg/tests/CMakeLists.txt
+++ b/modules/mol/alg/tests/CMakeLists.txt
@@ -2,6 +2,7 @@ set(OST_MOL_ALG_UNIT_TESTS
   test_superposition.cc
   tests.cc
   test_consistency_checks.cc
+  test_pdbize.py
   test_convenient_superpose.py
 )
 
diff --git a/modules/mol/alg/tests/test_pdbize.py b/modules/mol/alg/tests/test_pdbize.py
new file mode 100644
index 0000000000000000000000000000000000000000..abf2bc892f2618aa85f37b1ffd79c93d106b4d5d
--- /dev/null
+++ b/modules/mol/alg/tests/test_pdbize.py
@@ -0,0 +1,51 @@
+from ost import io, mol, geom, seq
+import unittest
+import os
+import random
+
+class TestPDBize(unittest.TestCase):
+  
+  def test_numbers_water_molecules_with_ins_codes(self):
+    m = mol.CreateEntity()
+    e = m.EditXCS(mol.BUFFERED_EDIT)
+    c = e.InsertChain("A");
+    e.SetChainType(c, mol.CHAINTYPE_WATER)
+    for i in range(27):
+      e.AppendResidue(c, "HOH")
+    pdbizer = mol.alg.PDBize()
+    transformations = geom.Mat4List()
+    transformations.append(geom.Mat4())
+    seqs = seq.CreateSequenceList()
+    pdbizer.Add(m.Select(''), transformations, seqs)
+    pdbized = pdbizer.Finish()
+    self.assertEquals([c.name for c in pdbized.chains], ["-"])
+    residues = pdbized.residues
+    for i in range(26):
+      self.assertEquals(residues[i].number.num, 1)
+      self.assertEquals(residues[i].number.ins_code, chr(ord('A')+i))
+    self.assertEquals(residues[26].number.num, 2)
+    self.assertEquals(residues[26].number.ins_code, 'A')
+  def test_starts_from_last_water_rnum(self):
+    m = mol.CreateEntity()
+    e = m.EditXCS(mol.BUFFERED_EDIT)
+    c = e.InsertChain("A");
+    e.SetChainType(c, mol.CHAINTYPE_WATER)
+    e.AppendResidue(c, "HOH")
+    pdbizer = mol.alg.PDBize()
+    transformations = geom.Mat4List()
+    transformations.append(geom.Mat4())
+    seqs = seq.CreateSequenceList()
+    pdbizer.Add(m.Select(''), transformations,seqs)
+    pdbizer.Add(m.Select(''), transformations,seqs)
+    pdbized = pdbizer.Finish()
+    self.assertEquals([c.name for c in pdbized.chains], ["-"])
+    residues = pdbized.residues
+    self.assertEquals([r.number for r in residues],
+                      [mol.ResNum(1, 'A'), mol.ResNum(1, 'B')])
+
+
+
+
+if __name__ == "__main__":
+  from ost import testutils
+  testutils.RunTests()