diff --git a/.gitignore b/.gitignore
index a987b49..58f5f60 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,3 +14,5 @@ simfiles/*
 .vscode/*
 .idea/*
 *.pkl
+
+*.xyz
diff --git a/control.example.yml b/control.example.yml
index fddc56c..0830d4d 100644
--- a/control.example.yml
+++ b/control.example.yml
@@ -11,13 +11,13 @@ diceplayer:
 
     dice:
         nmol: [1, 100]
-        dens: 1.0
+        dens: 1.5
         nstep: [2000, 3000]
         isave: 1000
         outname: 'phb'
         progname: '~/.local/bin/dice'
         ljname: 'phb.ljc'
-        randominit: 'first'
+        randominit: 'always'
 
     gaussian:
         qmprog: 'g16'
diff --git a/diceplayer/__main__.py b/diceplayer/__main__.py
index 2be65fc..1e4fa0e 100644
--- a/diceplayer/__main__.py
+++ b/diceplayer/__main__.py
@@ -51,6 +51,7 @@ def main():
         player.read_potentials()
 
         player.create_simulation_dir()
+        player.create_geoms_file()
 
     player.print_keywords()
 
diff --git a/diceplayer/player.py b/diceplayer/player.py
index 11814d5..677ec3d 100644
--- a/diceplayer/player.py
+++ b/diceplayer/player.py
@@ -10,6 +10,7 @@ from diceplayer.shared.utils.misc import weekday_date_time
 from diceplayer.shared.environment.molecule import Molecule
 from diceplayer.shared.environment.system import System
 from diceplayer.shared.environment.atom import Atom
+from diceplayer.shared.utils.ptable import atomsymb
 from diceplayer import logger
 
 from dataclasses import fields
@@ -19,6 +20,7 @@ import yaml
 import sys
 import os
 
+
 ENV = ["OMP_STACKSIZE"]
 
 
@@ -28,10 +30,8 @@ class Player:
             raise ValueError("Must specify either infile or optimization")
 
         elif infile is not None:
-            config_data = self.read_keywords(infile)
-
             self.config = self.set_config(
-                config_data.get("diceplayer")
+                self.read_keywords(infile)
             )
 
             self.system = System()
@@ -45,7 +45,7 @@ class Player:
 
             self.system = save[1]
 
-            self.initial_cycle = save[2]
+            self.initial_cycle = save[2] + 1
 
         else:
             raise ValueError("Must specify either infile or config")
@@ -60,7 +60,7 @@ class Player:
             "==========================================================================================\n"
         )
 
-        for cycle in range(self.initial_cycle, self.config.maxcyc + 1):
+        for cycle in range(self.initial_cycle, self.initial_cycle + self.config.maxcyc + 1):
 
             logger.info(
                 f"------------------------------------------------------------------------------------------\n"
@@ -72,8 +72,7 @@ class Player:
 
             try:
                 self.gaussian_start(cycle)
-                pass
-            except StopIteration as e:
+            except StopIteration:
                 break
 
             self.save_run_in_pickle(cycle)
@@ -105,6 +104,15 @@ class Player:
             )
         simulation_dir_path.mkdir()
 
+    def create_geoms_file(self):
+        geoms_file_path = Path(self.config.geoms_file)
+        if geoms_file_path.exists():
+            raise FileExistsError(
+                f"Error: a file or a directory {self.config.geoms_file} already exists,"
+                f" move or delete the simfiles directory to continue."
+            )
+        geoms_file_path.touch()
+
     def print_keywords(self) -> None:
 
         def log_keywords(config: Dataclass, dto: Type[Dataclass]):
@@ -295,6 +303,7 @@ class Player:
                 .update_charges(result['charges'])
 
             self.system.print_charges_and_dipole(cycle)
+            self.print_geoms(cycle)
 
             if diff < self.config.gaussian.chg_tol:
                 logger.info(
@@ -302,6 +311,18 @@ class Player:
                 )
                 raise StopIteration()
 
+    def print_geoms(self, cycle: int):
+        with open(self.config.geoms_file, 'a') as file:
+            file.write(f'Cycle # {cycle}\n')
+
+            for atom in self.system.molecule[0].atom:
+                symbol = atomsymb[atom.na]
+                file.write(
+                    f'{symbol:<2s}    {atom.rx:>10.6f}  {atom.ry:>10.6f}  {atom.rz:>10.6f}\n'
+                )
+
+            file.write('\n')
+
     @staticmethod
     def validate_atom_dict(molecule_type, molecule_site, atom_dict: dict) -> dict:
         molecule_type += 1
@@ -442,7 +463,14 @@ class Player:
     @staticmethod
     def read_keywords(infile) -> dict:
         with open(infile, 'r') as yml_file:
-            return yaml.load(yml_file, Loader=yaml.SafeLoader)
+            config = yaml.load(yml_file, Loader=yaml.SafeLoader)
+
+        if "diceplayer" in config:
+            return config.get("diceplayer")
+
+        raise RuntimeError(
+            f'Could not find diceplayer section in {infile}.'
+        )
 
     @classmethod
     def from_file(cls, infile: str) -> 'Player':
diff --git a/diceplayer/shared/config/player_config.py b/diceplayer/shared/config/player_config.py
index 75e2b12..ce1f7d9 100644
--- a/diceplayer/shared/config/player_config.py
+++ b/diceplayer/shared/config/player_config.py
@@ -23,6 +23,7 @@ class PlayerConfig(Dataclass):
     switchcyc: int = 3
     qmprog: str = 'g16'
     altsteps: int = 20000
+    geoms_file = 'geoms.xyz'
     simulation_dir = 'simfiles'
 
     def __post_init__(self):
diff --git a/diceplayer/shared/interface/dice_interface.py b/diceplayer/shared/interface/dice_interface.py
index 5b754a6..f79bc0a 100644
--- a/diceplayer/shared/interface/dice_interface.py
+++ b/diceplayer/shared/interface/dice_interface.py
@@ -166,11 +166,14 @@ class DiceInterface(Interface):
     def _make_init_file(self, proc_dir: Path, last_xyz_file: TextIO):
         xyz_lines = last_xyz_file.readlines()
 
-        nsites_mm = 0
-        for i in range(1, len(self.step.dice.nmol)):
-            nsites_mm += self.step.dice.nmol[i] * len(self.system.molecule[i].atom)
+        HEADER_LENGTH = 2
+        MAIN_MOLECULE_LENGTH = len(self.system.molecule[0].atom)
 
-        xyz_lines = xyz_lines[-nsites_mm:]
+        SECONDARY_MOLECULES_LENGTH = HEADER_LENGTH + MAIN_MOLECULE_LENGTH
+        for i in range(1, len(self.step.dice.nmol)):
+            SECONDARY_MOLECULES_LENGTH += self.step.dice.nmol[i] * len(self.system.molecule[i].atom)
+
+        xyz_lines = xyz_lines[HEADER_LENGTH + MAIN_MOLECULE_LENGTH: SECONDARY_MOLECULES_LENGTH]
 
         input_file = Path(proc_dir, self.step.dice.outname + ".xy")
         with open(input_file, 'w') as f: