diff --git a/modules/base/doc/generic.rst b/modules/base/doc/generic.rst
index 919497f390775a68deee2fcc973a22f29beace75..2d71160dae7dd0640fa506bbea95976e2bea0d2d 100644
--- a/modules/base/doc/generic.rst
+++ b/modules/base/doc/generic.rst
@@ -54,7 +54,7 @@ already exists, it will be overwritten. To check if it exists, use:
 .. code-block:: python
   
   exists=atom.HasProp("myfloatprop")
-  print exists
+  print(exists)
     
 To access the value of a generic property, we first check if the property exists
 and then access it, using the method suitable for the data type of the property. 
@@ -65,7 +65,7 @@ level:
   
   for atom in entity.GetAtomList(): 
     if atom.HasProp("myfloatprop"):
-      print atom.GetFloatProp("myfloatprop")
+      print(atom.GetFloatProp("myfloatprop"))
         
 When trying to access a property that has not been set, or one that has been 
 set, but at a different level, an error is thrown. The same is true when trying 
@@ -75,13 +75,13 @@ to access a property of a different data type, e.g.:
 
   # all of the following lines will throw errors
   # error because the property does not exist 
-  print atom.GetFloatProp("unknownprop")
+  print(atom.GetFloatProp("unknownprop"))
   
   # error because the property was set at another level
-  print entity.GetFloatProp("myfloatprop")
+  print(entity.GetFloatProp("myfloatprop"))
   
   # error because the data type of the property is different
-  print atom.GetStringProp("myfloatprop")
+  print(atom.GetStringProp("myfloatprop"))
       
 
 Use of Generic Properties in Queries
diff --git a/modules/base/doc/logging.rst b/modules/base/doc/logging.rst
index adb5a57682981fec8c588b7bb9aa04d9839c3b15..e9cff1360835ef23e9880f3c5949a5d810e3d41b 100644
--- a/modules/base/doc/logging.rst
+++ b/modules/base/doc/logging.rst
@@ -106,7 +106,7 @@ You can change the verbosity level  with the following two methods:
           ost.LogLevel.Info
           # Outputs: ost._ost_base.LogLevel.Info
 
-          print ost.LogLevel.Info
+          print(ost.LogLevel.Info)
           # Outputs: Info
           
           int(ost.LogLevel.Info)
@@ -267,7 +267,7 @@ terminal (or the python shell in DNG). The logger also prints the current time.
       levels=['ERROR', 'WARNING', 'SCRIPT', 'INFO', 
               'VERBOSE', 'DEBUG', 'TRACE']
       level=levels[severity]
-      print '%s[%s]: %s' % (level, str(datetime.datetime.now()), message),
+      print('%s[%s]: %s' % (level, str(datetime.datetime.now()), message))
 
   py_logger=PyLogger()
   ost.PushLogSink(py_logger)
diff --git a/modules/bindings/doc/dssp.rst b/modules/bindings/doc/dssp.rst
index bb9d2b57b0f081bd0f6ec0d7b82e0ef9bcf59089..3a51a870c9c740b495baa5caa53c0d1475818d77 100644
--- a/modules/bindings/doc/dssp.rst
+++ b/modules/bindings/doc/dssp.rst
@@ -39,7 +39,7 @@ using the mmCIF interface.
   for chain in ent.chains:
     if chain.is_polypeptide:
       for res in chain.residues:
-        print res.GetFloatProp('relative_solvent_accessibility')
+        print(res.GetFloatProp('relative_solvent_accessibility'))
 
 
 DSSP bindings Usage
diff --git a/modules/bindings/doc/hhblits.rst b/modules/bindings/doc/hhblits.rst
index 9ffb1c0c6f5d81b2a35b4299e54a30165bf656e1..62b3262a0a79bfbc30f9b0ce537cf161d52cc4cb 100644
--- a/modules/bindings/doc/hhblits.rst
+++ b/modules/bindings/doc/hhblits.rst
@@ -54,7 +54,7 @@ First query by sequence:
   with open(hit_file) as hit_fh:
       header, hits = hhblits.ParseHHblitsOutput(hit_fh)
   for hit in hits:
-      print hit.aln
+      print(hit.aln)
 
   # cleanup
   hh.Cleanup()
@@ -82,7 +82,7 @@ Very similar going by file:
   with open(hit_file) as hit_fh:
       header, hits = hhblits.ParseHHblitsOutput(hit_fh)
   for hit in hits:
-      print hit.aln
+      print(hit.aln)
 
   # cleanup
   hh.Cleanup()
@@ -106,7 +106,7 @@ so one may want to extract them:
   # note that ParseA3M is not a class method but a module function
   output = hhblits.ParseA3M(open(a3m_file))
 
-  print output['msa']
+  print(output['msa'])
 
   # cleanup
   hh.Cleanup()
diff --git a/modules/bindings/doc/tmtools.rst b/modules/bindings/doc/tmtools.rst
index 7a39ba987ece4c6ba688ae781802758dad6bd73b..51b786efced05180a0a0aacc8309a42f448005a1 100644
--- a/modules/bindings/doc/tmtools.rst
+++ b/modules/bindings/doc/tmtools.rst
@@ -61,9 +61,9 @@ structures and print the RMSD as well as the GDT_TS and GDT_HA similarity measur
   pdb1=io.LoadPDB('1ake.pdb', restrict_chains='A')
   pdb2=io.LoadPDB('4ake.pdb', restrict_chains='A')
   result=tmtools.TMScore(pdb1, pdb2)
-  print result.rmsd_below_five # 1.9
-  print result.gdt_ha # 0.41
-  print result.gdt_ts # 0.56
+  print(result.rmsd_below_five) # 1.9
+  print(result.gdt_ha) # 0.41
+  print(result.gdt_ts) # 0.56
 
 Usage of TMalign
 --------------------------------------------------------------------------------
@@ -96,8 +96,8 @@ generated in order to call the executable.
   pdb2=io.LoadPDB('4ake.pdb').Select("peptide=true")
   result = bindings.WrappedTMAlign(pdb1.chains[0], pdb2.chains[0], 
                                    fast=True)
-  print result.tm_score
-  print result.alignment.ToString(80)
+  print(result.tm_score)
+  print(result.alignment.ToString(80))
 
 
 .. class:: TMAlignResult(rmsd, tm_score, aligned_length, transform, alignment)
diff --git a/modules/conop/doc/compoundlib.rst b/modules/conop/doc/compoundlib.rst
index 195614fbd6d962dda1f2dff29f98515debc8b897..194979cfd388b4c1c129e1904cf83e1d26f01925 100644
--- a/modules/conop/doc/compoundlib.rst
+++ b/modules/conop/doc/compoundlib.rst
@@ -194,7 +194,7 @@ In this example we will translate the three-letter-codes given in the SEQRES rec
     compound=compound_lib.FindCompound(tlc)
     if compound:
        sequence+=compound.one_letter_code
-  print sequence # prints 'AGMVF'
+  print(sequence) # prints 'AGMVF'
 
 .. _mmcif-convert:
 
diff --git a/modules/doc/external.rst b/modules/doc/external.rst
index 3d75fa80b3fa1cc2217892159e0e4e9ca8f03d93..0d7f7d2047bc7319e0d168fb3b3b65a11afcf7b8 100644
--- a/modules/doc/external.rst
+++ b/modules/doc/external.rst
@@ -26,7 +26,7 @@ As an example, we would like to obtain the full path of the msms executable (a p
   from ost import settings
   exe_path = settings.Locate('msms', search_paths=['/opt/app','/home/app'],
               env_name='MSMS', search_system_paths=True)
-  print exe_path
+  print(exe_path)
   
 The :func:`~ost.settings.Locate` command looks for the program with the name 
 `msms`. If env_name is set, it first looks if an environment variable with the 
@@ -52,12 +52,12 @@ An example how to generate a temporary directory, open a file in this directory
   
   # generate a temporary directory
   tmp_dir_name = tempfile.mkdtemp()
-  print 'temporary directory:',tmp_dir_name
+  print('temporary directory:',tmp_dir_name)
   
   # generate and open a file in the temp directory
   tmp_file_name = os.path.join(tmp_dir_name,"entity")
   tmp_file_handle = open(tmp_file_name, 'w')
-  print 'temporary file:',tmp_file_handle
+  print('temporary file:',tmp_file_handle)
   
   # write position and radius of all atoms to file
   for a in entity.GetAtomList():
@@ -82,7 +82,7 @@ To run the external program msms from the above example, with the temporary file
   # set the command to execute
   command = "%s -if %s -of %s" % (exe_path,
             tmp_file_name, tmp_file_name)
-  print 'command:',command
+  print('command:',command)
 
   # run the executable with the command
   proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE)
@@ -90,11 +90,11 @@ To run the external program msms from the above example, with the temporary file
 
   # check for successful completion of msms
   if proc.returncode != 0:
-    print "WARNING: msms error\n", stdout_value
+    print("WARNING: msms error\n", stdout_value)
     raise subprocess.CalledProcessError(proc.returncode, command)
 
   # print everything written to the command line (stdout)
-  print stdout_value
+  print(stdout_value)
     
 Read Generated Output
 --------------------------------------------------------------------------------
@@ -106,8 +106,8 @@ Here we first print the command line output and then load the generated msms sur
 .. code-block:: python
 
   # print everything written to the command line (stdout)
-  print stdout_value
+  print(stdout_value)
   
   # read msms surface from file
   surface = io.LoadSurface(tmp_file_name, "msms")
-  print 'number of vertices:',len(surface.GetVertexIDList())
+  print('number of vertices:',len(surface.GetVertexIDList()))
diff --git a/modules/doc/intro-01.rst b/modules/doc/intro-01.rst
index 25c092f8ae5801f258b206d2b409498ea0e3ae8d..98c3427e7a445adab8f4e5ad020996ad0e18d50b 100644
--- a/modules/doc/intro-01.rst
+++ b/modules/doc/intro-01.rst
@@ -47,9 +47,9 @@ The loaded structure is an instance of :class:`~ost.mol.EntityHandle` which offe
 
 .. code-block:: python
 
-   print len(fragment.chains), fragment.chains
-   print len(fragment.residues), fragment.residues
-   print len(fragment.atoms), fragment.atoms
+   print(len(fragment.chains), fragment.chains)
+   print(len(fragment.residues), fragment.residues)
+   print(len(fragment.atoms), fragment.atoms)
 
 As you can see, our fragment consists of one peptide chain of 12 amino acids and 
 has 81 atoms in total. Now let's examine our fragment in more detail. Enter the 
@@ -58,9 +58,9 @@ command
 .. code-block:: python
 
   for residue in fragment.residues:
-    print residue, 'has', len(residue.atoms), 'atom(s).'
+    print(residue, 'has', len(residue.atoms), 'atom(s).')
     for atom in residue.atoms:
-      print ' ', atom.name, atom.pos
+      print(' ', atom.name, atom.pos)
 
 
 This will group the atoms by residue. And, for completeness, we will first group them by chain, then by residues.
@@ -68,11 +68,11 @@ This will group the atoms by residue. And, for completeness, we will first group
 .. code-block:: python
 
   for chain in fragment.chains:
-    print 'chain', chain.name, 'has', len(chain.residues), 'residue(s)'
+    print('chain', chain.name, 'has', len(chain.residues), 'residue(s)')
     for residue in chain.residues:
-      print ' ', residue, 'has', len(residue.atoms), 'atom(s).'
+      print(' ', residue, 'has', len(residue.atoms), 'atom(s).')
       for atom in residue.atoms:
-        print '    ', atom.name, atom.pos
+        print('    ', atom.name, atom.pos)
 
 A protein fragment would not be complete without bonds. Let's see 
 what bonds we have in there:
@@ -80,7 +80,7 @@ what bonds we have in there:
 .. code-block:: python
   
   for bond in fragment.bonds:
-    print bond
+    print(bond)
     
 Let There Be Shiny Graphics
 --------------------------------------------------------------------------------
diff --git a/modules/doc/intro-02.rst b/modules/doc/intro-02.rst
index 4470a9d8c21746abe5c2a3f1b9bcc8deacebe33f..67bfc4ddbed9cefbf05f52e871c1d835578cc6d2 100644
--- a/modules/doc/intro-02.rst
+++ b/modules/doc/intro-02.rst
@@ -30,7 +30,7 @@ Now let's inspect what we just loaded:
 
 .. code-block:: python
 
-  print map.GetPixelSampling(), map.GetSize()
+  print(map.GetPixelSampling(), map.GetSize())
     
 We can see that the sampling is set to 1.0 Angstrom in all three dimensions. The loaded map is an instance of :class:`~ost.img.ImageHandle`, a class to represent images in one, two and three dimensions.
 
diff --git a/modules/doc/table.rst b/modules/doc/table.rst
index d0dd45691b89094bba4815fab8ded65e7badd206..4a8748e50bd9a72c7fa1ab01d02c2f926bc33d3b 100644
--- a/modules/doc/table.rst
+++ b/modules/doc/table.rst
@@ -42,14 +42,14 @@ Iterating over table items:
   # iterate over all rows
   for row in tab.rows:
     # print complete row
-    print row
+    print(row)
 
     # print value for column 'foo'
-    print row[idx]
+    print(row[idx])
 
   # iterate over all rows of selected columns
   for foo, bar in tab.Zip('foo','bar'):
-    print foo, bar
+    print(foo, bar)
 
 Doing element wise mathematical operations on entire colums:
 
@@ -63,7 +63,7 @@ Doing element wise mathematical operations on entire colums:
   # addition of column foo and column bar
   tab.AddCol('qux', 'f', tab['foo']+tab['bar'])
 
-  print tab
+  print(tab)
 
 
 Select part of the table based on a query:
@@ -76,11 +76,11 @@ Select part of the table based on a query:
 
   # select all rows where foo>=2 and bar<10
   subtab = tab.Select('foo>=2 and bar<10')
-  print subtab
+  print(subtab)
 
   # select all rows where foo>3 or bar=1
   subtab = tab.Select('foo>3 or bar=1')
-  print subtab
+  print(subtab)
 
 
 Functions You Might be Interested In
diff --git a/modules/geom/doc/mat.rst b/modules/geom/doc/mat.rst
index c8736733a7a18533cf18e6fcf40ddecf6b3a52c7..a05dab46569994c22ba46d1064e1ba0b52e7b444 100644
--- a/modules/geom/doc/mat.rst
+++ b/modules/geom/doc/mat.rst
@@ -11,8 +11,8 @@ the following code examples:
 .. code-block:: python
 
   m=geom.Mat2(1, 2, 3, 4)
-  print m # will print {{1,2},{3,4}}
-  print m[(0,0)], m[(0,1)], m[(1,0)], m[(1,1)] # will print 1, 2, 3, 4
+  print(m) # will print {{1,2},{3,4}}
+  print(m[(0,0)], m[(0,1)], m[(1,0)], m[(1,1)]) # will print 1, 2, 3, 4
 
 Matrices support arithmetic via overloaded operators. The following operations are 
 supported:
@@ -164,4 +164,4 @@ Functions Operating on Matrices
   (0, 0, 0).
   
   :param vec: A vector of arbitrary length
-  :type vec: :class:`Vec3`
\ No newline at end of file
+  :type vec: :class:`Vec3`
diff --git a/modules/geom/doc/vec.rst b/modules/geom/doc/vec.rst
index bc88b3f2b0e020bdfa0b96480c561414e8f2fa7f..c55dcbb658b82859514addc548884e7fcd6edc6d 100644
--- a/modules/geom/doc/vec.rst
+++ b/modules/geom/doc/vec.rst
@@ -16,9 +16,9 @@ This is shown in the following example:
    
    vec_a=geom.Vec2(1, 0)
    vec_b=geom.Vec2(0, 1)
-   print vec_a, vec_b
-   print vec_a+vec_b
-   print vec_a*3-vec_b
+   print(vec_a, vec_b)
+   print(vec_a+vec_b)
+   print(vec_a*3-vec_b)
 
 The standard vector operations are implemented as :ref:`free standing functions 
 <vector-functions>`:
@@ -29,12 +29,12 @@ The standard vector operations are implemented as :ref:`free standing functions
   vec_a=geom.Vec3(1, 0, 0)
   vec_b=geom.Vec3(0, 1, 0)
   
-  print geom.Dot(vec_a, vec_b)
-  print geom.Cross(vec_a, vec_b)
+  print(geom.Dot(vec_a, vec_b))
+  print(geom.Cross(vec_a, vec_b))
   
-  print geom.Normalize(geom.Vec3(1, 1, 0))
+  print(geom.Normalize(geom.Vec3(1, 1, 0)))
   
-  print geom.Length(geom.Vec3(1, 1, 1))
+  print(geom.Length(geom.Vec3(1, 1, 1)))
 
 
 Vector Classes
diff --git a/modules/gui/doc/python_cpp.rst b/modules/gui/doc/python_cpp.rst
index 0f4be50a63797040a7965daa68b143f9b9b8204c..f7f726251dbae3740036eaa3a06d6cae0b397ecb 100644
--- a/modules/gui/doc/python_cpp.rst
+++ b/modules/gui/doc/python_cpp.rst
@@ -12,7 +12,7 @@ which wraps the Object into a Python SIP Object.
    
   seq_viewer = gui.SequenceViewer() # Create SequenceViewer Object
   qobj = seq_viewer.qobject #Get Python SIP Object
-  print qobj.size()   # Call function on QWidget
+  print(qobj.size())   # Call function on QWidget
 
 The other way around, each boost::python Qt Object accepts Python objects as 
 input for Qt Objects. It handles the cast to a C++ Qt Object internally.
diff --git a/modules/img/base/doc/point-size-extent.rst b/modules/img/base/doc/point-size-extent.rst
index af2731adb21d7ffae75bb859d78f16835b392965..344f0ea288ea09c7be22544a230a2735b8d4a1a1 100644
--- a/modules/img/base/doc/point-size-extent.rst
+++ b/modules/img/base/doc/point-size-extent.rst
@@ -39,7 +39,7 @@ Point
     ..code-block:: python
     
       p=img.Point(1,2,3)
-      print p
+      print(p)
       p[1]=5
 
   .. method:: ToVec3()
@@ -140,4 +140,4 @@ Extent
     # visits all the pixels in the extent and
     # prints out their values
     for pixel in ei:
-      print i.GetReal(pixel)
\ No newline at end of file
+      print(i.GetReal(pixel))
diff --git a/modules/io/doc/io.rst b/modules/io/doc/io.rst
index 995950ed29fcf6a4c0d22a8606f64463f841ca01..10be9a05532493ccaf23afa06cd49be89a8887a1 100644
--- a/modules/io/doc/io.rst
+++ b/modules/io/doc/io.rst
@@ -478,28 +478,6 @@ Saving Density Maps
     io.SaveImage(image, 'new_map.map', CCP4())
 
 
-
-
-
-.. testsetup:: io
-  from ost import io
-
-.. testcode:: io
-  :hide:
-
-  from ost import io,seq
-  ent = io.LoadPDB('./examples/entity/fragment.pdb')
-  print ent.atom_count 
-  myseq = seq.SequenceFromChain('t',ent.chains[0])
-  print myseq.GetLength()
-
-.. testoutput:: io
-  :hide:
-
-  81
-  12
-
-
 Stereochemical Parameters
 --------------------------------------------------------------------------------
 
diff --git a/modules/mol/alg/doc/lddt.rst b/modules/mol/alg/doc/lddt.rst
index 3201243bdaa59551cc73e543e16100c161bfb5f4..711c06f9f89a39ff1d8f9a394270756d6dcd1eb2 100644
--- a/modules/mol/alg/doc/lddt.rst
+++ b/modules/mol/alg/doc/lddt.rst
@@ -240,7 +240,7 @@ One can replicate the binary using simple python script:
     #
     # Check consistency
     is_cons = ResidueNamesMatch(model_view, references[0], True)
-    print "Consistency check: ", "OK" if is_cons else "ERROR"
+    print("Consistency check: ", "OK" if is_cons else "ERROR")
     #
     # Calculate lDDT
     LocalDistDiffTest(model_view,
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 300852d954f6869c02ed7d2410ba9d3f8f958bf7..5262c236a20247a08e1bd5c9db0114447ffe1ac3 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -576,7 +576,7 @@ Local Distance Test scores (lDDT, DRMSD)
     #
     # Calculate lDDT
     scorer = lDDTScorer(references=references, model=model_view, settings=settings)
-    print "Global score:", scorer.global_score
+    print("Global score:", scorer.global_score)
     scorer.PrintPerResidueStats()
   
   :param references: Sets :attr:`references`
diff --git a/modules/mol/base/doc/entity.rst b/modules/mol/base/doc/entity.rst
index fd99d80cc0286876d056095826fe646e28e93762..1e7b008861df783c0970fbc2af4a14975f939f38 100644
--- a/modules/mol/base/doc/entity.rst
+++ b/modules/mol/base/doc/entity.rst
@@ -58,7 +58,7 @@ The Handle Classes
     .. code-block:: python
     
       for res in ent.residues:
-        print res.name, res.atom_count
+        print(res.name, res.atom_count)
     
     Also available as :meth:`GetResidueList`. To access a single residue, use 
     :meth:`FindResidue`. 
@@ -380,7 +380,7 @@ The Handle Classes
      
        chain=ent.FindChain("A")
        for res in chain.residues:
-         print res.name, res.atom_count
+         print(res.name, res.atom_count)
    
      Also available as :meth:`GetResidueList`. To access a single residue, use 
      :meth:`FindResidue`. 
@@ -394,12 +394,12 @@ The Handle Classes
      .. code-block:: python
      
        chain=ent.FindChain("A")
-       print chain.residues # [A.GLY1, A.GLY2, A.GLY4A, A.GLY4B]
-       print chain.in_sequence # prints true
+       print(chain.residues) # [A.GLY1, A.GLY2, A.GLY4A, A.GLY4B]
+       print(chain.in_sequence) # prints true
        
        chain=ent.FindChain("B")
-       print chain.residues # [B.GLY1, B.GLY4, B.GLY3]
-       print chain.in_sequence # prints false
+       print(chain.residues) # [B.GLY1, B.GLY4, B.GLY3]
+       print(chain.in_sequence) # prints false
 
   .. attribute:: residue_count
 
@@ -537,7 +537,7 @@ The Handle Classes
     
     .. code-block:: python
     
-      print ''.join([r.one_letter_code for r in chain.residues])
+      print(''.join([r.one_letter_code for r in chain.residues]))
     
     :type: str
     
@@ -1063,7 +1063,7 @@ The View Classes
     .. code-block:: python
     
       for res in ent.residues:
-        print res.name, res.atom_count
+        print(res.name, res.atom_count)
     
     Also available as :meth:`GetResidueList`. To access a single residue, use 
     :meth:`FindResidue`. 
@@ -1410,7 +1410,7 @@ The View Classes
      
        chain=ent.FindChain("A")
        for res in chain.residues:
-         print res.name, res.atom_count
+         print(res.name, res.atom_count)
    
      Also available as :meth:`GetResidueList`. To access a single residue, use
      :meth:`FindResidue`.
@@ -1424,12 +1424,12 @@ The View Classes
      .. code-block:: python
      
        chain=ent.FindChain("A")
-       print chain.residues # [A.GLY1, A.GLY2, A.GLY4A, A.GLY4B]
-       print chain.in_sequence # prints true
+       print(chain.residues) # [A.GLY1, A.GLY2, A.GLY4A, A.GLY4B]
+       print(chain.in_sequence) # prints true
        
        chain=ent.FindChain("B")
-       print chain.residues # [B.GLY1, B.GLY4, B.GLY3]
-       print chain.in_sequence # prints false
+       print(chain.residues) # [B.GLY1, B.GLY4, B.GLY3]
+       print(chain.in_sequence) # prints false
 
   .. attribute:: atoms
 
@@ -1653,7 +1653,7 @@ The View Classes
   
     .. code-block:: python
   
-      print ''.join([r.one_letter_code for r in chain.residues])
+      print(''.join([r.one_letter_code for r in chain.residues]))
   
     :type: str
   
diff --git a/modules/seq/alg/doc/seqalg.rst b/modules/seq/alg/doc/seqalg.rst
index 265892530eb97227216348523af3a1830c656c7e..5e5b5a7c864906d37c202768746d4c38f41b6c32 100644
--- a/modules/seq/alg/doc/seqalg.rst
+++ b/modules/seq/alg/doc/seqalg.rst
@@ -60,14 +60,14 @@ Algorithms for Alignments
        aln_a = seq.CreateAlignment()
        aln_a.AddSequence(seq_a1)
        aln_a.AddSequence(seq_a2)
-       print aln_a
+       print(aln_a)
        # >>> A1  acdefghikl-mn
        # >>> A2  atd-fghikllmn
 
        aln_b = seq.CreateAlignment()
        aln_b.AddSequence(seq_b1)
        aln_b.AddSequence(seq_b2)
-       print aln_b
+       print(aln_b)
        # >>> B1  acdefg-hiklmn
        # >>> B2  acd---qhirlmn
 
@@ -76,7 +76,7 @@ Algorithms for Alignments
        aln_list.append(aln_b)
 
        merged_aln = ost.seq.alg.MergePairwiseAlignments(aln_list, ref_seq)
-       print merged_aln
+       print(merged_aln)
        # >>> ref  acdefg-hikl-mn
        # >>> A2   atd-fg-hikllmn
        # >>> B2   acd---qhirl-mn
@@ -121,7 +121,7 @@ Algorithms for Alignments
     seq_a = seq.CreateSequence('A', 'acdefghiklmn')
     seq_b = seq.CreateSequence('B', 'acdhiklmn')
     alns = seq.alg.LocalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
-    print alns[0].ToString(80)
+    print(alns[0].ToString(80))
     # >>> A acdefghiklmn
     # >>> B acd---hiklmn
 
@@ -150,7 +150,7 @@ Algorithms for Alignments
     seq_a = seq.CreateSequence('A', 'acdefghiklmn')
     seq_b = seq.CreateSequence('B', 'acdhiklmn')
     alns = seq.alg.GlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
-    print alns[0].ToString(80)
+    print(alns[0].ToString(80))
     # >>> A acdefghiklmn
     # >>> B acd---hiklmn
 
@@ -193,11 +193,11 @@ Algorithms for Alignments
     seq_a = seq.CreateSequence('A', 'abcdefghijklmnok')
     seq_b = seq.CreateSequence('B', 'cdehijk')
     alns = seq.alg.GlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
-    print alns[0].ToString(80)
+    print(alns[0].ToString(80))
     # >>> A abcdefghijklmnok
     # >>> B --cde--hi-----jk
     alns = seq.alg.SemiGlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
-    print alns[0].ToString(80)
+    print(alns[0].ToString(80))
     # >>> A abcdefghijklmnok
     # >>> B --cde--hijk-----
 
@@ -370,7 +370,7 @@ differences between the structures.
   dist_to_mean = seq.alg.CreateDist2Mean(d_map)
 
   # report min. and max. variances
-  print "MIN-MAX:", var_map.Min(), "-", var_map.Max()
+  print("MIN-MAX:", var_map.Min(), "-", var_map.Max())
   # get data and json-strings for further processing
   var_map_data = var_map.GetData()
   var_map_json = var_map.GetJsonString()
diff --git a/modules/seq/base/doc/seq.rst b/modules/seq/base/doc/seq.rst
index 8e7e91f57722c1d28ab3a37a444a4a4c4c7a7d9a..510241c8b95dddf82466800e33bc2c67c6e1cc59 100644
--- a/modules/seq/base/doc/seq.rst
+++ b/modules/seq/base/doc/seq.rst
@@ -56,7 +56,7 @@ input  methods, sequences can also be loaded from a string:
   seq_string = '''>sequence
   abcdefghiklmnop'''
   s = io.SequenceFromString(seq_string, 'fasta')
-  print s.name, s # will print "sequence abcdefghiklmnop"
+  print(s.name, s) # will print "sequence abcdefghiklmnop"
   
 Note that, in that case specifying the format is mandatory.
 
@@ -92,8 +92,8 @@ The SequenceHandle
     .. code-block:: python
       
       s=seq.CreateSequence("A", "abc---def")
-      print s.GetPos(1) # prints 1
-      print s.GetPos(3) # prints 6
+      print(s.GetPos(1)) # prints 1
+      print(s.GetPos(3)) # prints 6
     
     The reverse mapping, that is from position in the sequence to residue index 
     can be achieved with :meth:`GetResidueIndex`.
@@ -108,12 +108,12 @@ The SequenceHandle
     .. code-block:: python
       
       s=seq.CreateSequence("A", "abc--def")
-      print s.GetResidueIndex(1) # prints 1
-      print s.GetResidueIndex(6) # prints 4
+      print(s.GetResidueIndex(1)) # prints 1
+      print(s.GetResidueIndex(6)) # prints 4
       # the following line raises an exception of type
       # Error with the message "requested position contains 
       # a gap"
-      print s.GetResidueIndex(3)
+      print(s.GetResidueIndex(3))
 
   .. method:: GetResidue(pos)
      
@@ -319,11 +319,11 @@ an alignment:
   aln=io.LoadAlignment('aln.fasta')
   # iterate over the columns
   for col in aln:
-    print col
+    print(col)
 
   # iterate over the sequences
   for s in aln.sequences:
-    print s
+    print(s)
 
 .. function:: CreateAlignment()
 
@@ -376,7 +376,7 @@ an alignment:
       aln.AddSequence(seq.CreateSequence("A", "abcdefghik"))
       aln.AddSequence(seq.CreateSequence("B", "1234567890"))
       # The following command will print the output given below
-      print aln.ToString(7)
+      print(aln.ToString(7))
       # A abcde
       # B 12345
       #
@@ -431,7 +431,7 @@ an alignment:
       aln.AddSequence(seq.CreateSequence("B", "1234567890"))
       aln.Cut(4, 7)
       
-      print aln.ToString(80)
+      print(aln.ToString(80))
       # will print
       # A abcdhik
       # B 1234890