diff --git a/propka/input.py b/propka/input.py
index 675a770..3971f12 100644
--- a/propka/input.py
+++ b/propka/input.py
@@ -68,11 +68,6 @@ def read_molecule_file(input_file, mol_container):
             mol_container.conformations[name].sort_atoms()
         # find coupled groups
         mol_container.find_covalently_coupled_groups()
-        # write out the input file
-        # TODO - figure out why this I/O has to happen here
-        output_path = Path(input_path.name.replace(
-            input_file_extension, '.propka_input'))
-        write_propka(mol_container, output_path)
     elif input_file_extension.lower() == '.propka_input':
         # input is a propka_input file
         conformations, conformation_names = read_propka(
diff --git a/propka/lib.py b/propka/lib.py
index 518f7e6..f7f5774 100644
--- a/propka/lib.py
+++ b/propka/lib.py
@@ -279,6 +279,9 @@ def build_parser(parser=None):
     group.add_argument(
         "-q", "--quiet", action="store_const", const="WARNING",
         dest="log_level", help="suppress non-warning messages")
+    group.add_argument(
+        "--generate-propka-input", action="store_true",
+        help="Generate a PROPKA input file")
     group.add_argument(
         "--protonate-all", dest="protonate_all", action="store_true",
         help="Protonate all atoms (will not influence pKa calculation)",
diff --git a/propka/molecular_container.py b/propka/molecular_container.py
index 43ad372..d1492ae 100644
--- a/propka/molecular_container.py
+++ b/propka/molecular_container.py
@@ -1,6 +1,7 @@
 """Molecular container for storing all contents of PDB files."""
 import os
 import propka.version
+from propka.output import write_propka, write_pka, print_header, print_result
 from propka.conformation_container import ConformationContainer
 from propka.lib import info, warning, make_grid
 
@@ -26,7 +27,7 @@ class MolecularContainer:
             options:  options object
         """
         # printing out header before parsing input
-        propka.output.print_header()
+        print_header()
         self.conformation_names = []
         self.conformations = {}
         self.options = options
@@ -82,7 +83,7 @@ class MolecularContainer:
         # find the average of the conformations
         self.average_of_conformations()
         # print out the conformation-average results
-        propka.output.print_result(self, 'AVR', self.version.parameters)
+        print_result(self, 'AVR', self.version.parameters)
 
     def average_of_conformations(self):
         """Generate an average of conformations."""
@@ -117,6 +118,16 @@ class MolecularContainer:
             self.conformation_names[0]].chains
         self.conformations['AVR'] = avr_conformation
 
+    def write_propka(self, filename=None):
+        """Write PROPKA input file.
+
+        Args:
+            filename:  file to write to
+        """
+        if filename is None:
+            filename = os.path.join('{0:s}.propka_input'.format(self.name))
+        write_propka(self, filename)
+
     def write_pka(self, filename=None, reference="neutral",
                   direction="folding", options=None):
         """Write pKa information to a file.
@@ -127,8 +138,8 @@ class MolecularContainer:
             direction:  folding vs. unfolding
             options:  options object
         """
-        # write out the average conformation
-        filename = os.path.join('{0:s}.pka'.format(self.name))
+        if filename is None:
+            filename = os.path.join('{0:s}.pka'.format(self.name))
         # if the display_coupled_residues option is true, write the results out
         # to an alternative pka file
         if self.options.display_coupled_residues:
@@ -138,7 +149,7 @@ class MolecularContainer:
             filename = os.path.join(
                 '{0:s}_{1:s}.pka'.format(
                     self.name, self.version.parameters.output_file_tag))
-        propka.output.write_pka(
+        write_pka(
             self, self.version.parameters, filename=filename,
             conformation='AVR', reference=reference)
 
diff --git a/propka/run.py b/propka/run.py
index 633e2c5..0769715 100644
--- a/propka/run.py
+++ b/propka/run.py
@@ -21,6 +21,8 @@ def main(optargs=None):
         my_molecule = read_molecule_file(pdbfile, my_molecule)
         my_molecule.calculate_pka()
         my_molecule.write_pka()
+        if options.generate_propka_input:
+            my_molecule.write_propka()
 
 
 def single(pdbfile, optargs=None):
@@ -44,4 +46,6 @@ def single(pdbfile, optargs=None):
     my_molecule = read_molecule_file(pdbfile, my_molecule)
     my_molecule.calculate_pka()
     my_molecule.write_pka()
+    if options.generate_propka_input:
+        my_molecule.write_propka()
     return my_molecule
diff --git a/scripts/propka31.py b/scripts/propka31.py
index 8a90114..faa5b46 100755
--- a/scripts/propka31.py
+++ b/scripts/propka31.py
@@ -28,6 +28,8 @@ def main():
         my_molecule = read_molecule_file(pdbfile, my_molecule)
         my_molecule.calculate_pka()
         my_molecule.write_pka()
+        if options.generate_propka_input:
+            my_molecule.write_propka()
 
 
 if __name__ == '__main__':
diff --git a/tests/test_basic_regression.py b/tests/test_basic_regression.py
index 9a63854..fc78a21 100644
--- a/tests/test_basic_regression.py
+++ b/tests/test_basic_regression.py
@@ -78,6 +78,8 @@ def run_propka(options, pdb_path, tmp_path):
         molecule = read_molecule_file(str(pdb_path), molecule)
         molecule.calculate_pka()
         molecule.write_pka()
+        if args.generate_propka_input:
+            molecule.write_propka()
     finally:
         os.chdir(cwd)