Index: modlib/modeller/automodel/automodel.py
===================================================================
--- modlib/modeller/automodel/automodel.py	(revision 6112)
+++ modlib/modeller/automodel/automodel.py	(revision 6113)
@@ -7,7 +7,7 @@
 from modeller.energy_data import energy_data
 from modeller.selection import selection
 import modeller.modfile as modfile
-from modeller import log, physical
+from modeller import log, physical, io_data
 from modeller.scripts import align_strs_seq
 from modeller.optimizers import conjugate_gradients, molecular_dynamics, actions
 import refine
@@ -233,11 +233,19 @@
             job.queue_task(modeltask(self, num, atmsel))
         self.outputs.extend(job.run_all_tasks())
 
+    def read_initial_model(self):
+        """Read the initial model from a file.
+           Note that you are counting on some model arrays not being deleted by
+           read() (ie the charge array generated by generate_topology())."""
+        # Make sure we read every atom, since when we write out the model, it
+        # writes out every atom, not just non-HET/non-hydrogen/non-water
+        io = io_data(copy=self.env.io)
+        io.hetatm = io.water = io.hydrogen = True
+        self.read(file=self.inifile, io=io)
+
     def randomize_initial_structure(self, atmsel):
-        """Vary the initial structure"""
-        # Note that you are counting on some MODEL arrays not being deleted by
-        # read() (ie the charge(1:natm) array generated by generate_topology())
-        self.read(file=self.inifile)
+        """Get and randomize the initial structure"""
+        self.read_initial_model()
         if self.rand_method:
             self.rand_method(atmsel)
 
@@ -431,18 +439,16 @@
         self._check_model_hetatm()
 
     def _check_model_hetatm(self):
-        """Check to make sure env.io.hetatm is set if we're using HETATMs"""
+        """Check to see if env.io.hetatm is set if we're using HETATMs"""
         if not self.env.io.hetatm:
             het = selection(self).only_het_residues()
             if len(het) > 0:
-                self.env.io.hetatm = True
                 log.warning("_check_model_hetatm",
 """You have at least one HETATM residue in your model, but
-              io_data.hetatm is False. (This prevents Modeller from reading
-              HETATM template information, or from reading the .ini file back
-              in correctly.) This flag has now been automatically turned on.
-              To avoid this warning in future, set env.io.hetatm = True before
-              creating the automodel or loopmodel object.""")
+              io_data.hetatm is False. (This means that Modeller will not read
+              any HETATM data from your templates, which is usually not what
+              you want. To fix this, set env.io.hetatm = True before creating
+              the automodel or loopmodel object.""")
 
     def build_charmm_restraints(self, atmsel, rsr, aln):
         """Build restraints from CHARMM libraries"""
@@ -728,7 +734,7 @@
 
     def cluster(self, cluster_cut=1.5):
         """Cluster all output models, and output an optimized cluster average"""
-        self.read(file=self.inifile)
+        self.read_initial_model()
         aln = alignment(self.env)
         self.align_models(aln)
         if len(aln) == 0: