diff --git a/DPpack/MolHandling.py b/DPpack/MolHandling.py
index d47d817..6f5d78e 100644
--- a/DPpack/MolHandling.py
+++ b/DPpack/MolHandling.py
@@ -1,4 +1,3 @@
-from DPpack.MolHandling import total_mass
 import os, sys
 import math
 import shutil
@@ -11,7 +10,11 @@ from numpy import linalg
 
 from DPpack.Misc import *
 from DPpack.PTable import *
-from DPpack.SetGlobals import *
+
+env = ["OMP_STACKSIZE"]
+
+bohr2ang = 0.52917721092
+ang2bohr = 1/bohr2ang
 
 # Usaremos uma nova classe que ira conter toda interação entre moleculas
 
@@ -173,9 +176,9 @@ class System:
 
 	def print_geom(self, cycle, fh):
 	
-		fh.write("{}\n".format(len(self.molecule[0])))
+		fh.write("{}\n".format(len(self.molecule[0].atom)))
 		fh.write("Cycle # {}\n".format(cycle))
-		for atom in self.molecule[0].atoms:
+		for atom in self.molecule[0].atom:
 			symbol = atomsymb[atom.na]
 			fh.write("{:<2s}    {:>10.6f}  {:>10.6f}  {:>10.6f}\n".format(symbol, 
 														atom.rx, atom.ry, atom.rz))
@@ -201,6 +204,8 @@ class Molecule:
 		self.atom.append(a)			# Inserção de um novo atomo
 		self.total_mass += a.mass
 
+		self.center_of_mass()
+
 	def center_of_mass(self):
 	
 		com = np.zeros(3)
@@ -221,9 +226,9 @@ class Molecule:
 
 		for atom in self.atom:
 
-			atom.rx -= com[0]
-			atom.ry -= com[1]
-			atom.rz -= com[2]
+			atom.rx -= self.com[0]
+			atom.ry -= self.com[1]
+			atom.rz -= self.com[2]
 	
 	def charges_and_dipole(self):
 	
@@ -285,15 +290,14 @@ class Molecule:
 
 	def inertia_tensor(self):
 	
-		com = self.center_of_mass()
 		Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0
 
 		for atom in self.atom:
 
 			####  Obtain the displacement from the center of mass
-			dx = atom.rx - com[0]
-			dy = atom.ry - com[1]
-			dz = atom.rz - com[2]
+			dx = atom.rx - self.com[0]
+			dy = atom.ry - self.com[1]
+			dz = atom.rz - self.com[2]
 			####  Update the diagonal components of the tensor
 			Ixx += atom.mass * (dy**2 + dz**2)
 			Iyy += atom.mass * (dz**2 + dx**2)
@@ -332,7 +336,7 @@ class Molecule:
 		
 		mat1 = 1/np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient)
 		mat2 = 1/np.dot(step, np.matmul(self.hessian, step.T).T)
-		mat2 *= np.matmul( np.matmul(self.hessian, step.T), np.matmul(step, hessian) )
+		mat2 *= np.matmul( np.matmul(self.hessian, step.T), np.matmul(step, self.hessian) )
 		
 		self.hessian += mat1 - mat2
 
@@ -393,10 +397,9 @@ class Molecule:
 		return new_molecule
 
 	def print_mol_info(self, fh):
-	
-		com = self.center_of_mass()
-		fh.write("    Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(com[0], 
-																			com[1], com[2]))
+
+		fh.write("    Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(self.com[0], 
+																			self.com[1], self.com[2]))
 		inertia = self.inertia_tensor()
 		evals, evecs = self.principal_axes()
 		
@@ -413,60 +416,13 @@ class Molecule:
 		sizes = self.sizes_of_molecule()
 		fh.write("    Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
 																sizes[0], sizes[1], sizes[2]))
-		mol_mass = self.total_mass()
-		fh.write("    Total mass = {:>8.2f} au\n".format(mol_mass))
+		fh.write("    Total mass = {:>8.2f} au\n".format(self.total_mass))
 		
 		chg_dip = self.charges_and_dipole()
 		fh.write("    Total charge = {:>8.4f} e\n".format(chg_dip[0]))
 		fh.write("    Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} )     Total = {:>9.4f} Debye\n\n".format(
 											chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]))
 
-	def calculate_step(self, fh):
-		
-		invhessian = linalg.inv(self.hessian)
-		pre_step = -1 * np.matmul(invhessian, self.gradient.T).T
-		maxstep = np.amax(np.absolute(pre_step))
-		factor = min(1, player['maxstep']/maxstep)
-		step = factor * pre_step
-		
-		fh.write("\nCalculated step:\n")
-		pre_step_list = pre_step.tolist()
-		
-		fh.write("-----------------------------------------------------------------------\n"
-				"Center     Atomic                          Step (Bohr)\n"
-				"Number     Number              X                Y                Z\n"
-				"-----------------------------------------------------------------------\n")
-		for i in range(len(molecules[0])):
-			fh.write("  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
-																i + 1, molecules[0][i]['na'], 
-							pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0)))
-		
-		fh.write("-----------------------------------------------------------------------\n")
-		
-		fh.write("Maximum step is {:>11.6}\n".format(maxstep))
-		fh.write("Scaling factor = {:>6.4f}\n".format(factor))
-		fh.write("\nFinal step (Bohr):\n")
-		step_list = step.tolist()
-		
-		fh.write("-----------------------------------------------------------------------\n"
-				"Center     Atomic                          Step (Bohr)\n"
-				"Number     Number              X                Y                Z\n"
-				"-----------------------------------------------------------------------\n")
-		for i in range(len(molecules[0])):
-			fh.write("  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
-																i + 1, molecules[0][i]['na'], 
-										step_list.pop(0), step_list.pop(0), step_list.pop(0)))
-		
-		fh.write("-----------------------------------------------------------------------\n")
-		
-		step_max = np.amax(np.absolute(step))
-		step_rms = np.sqrt(np.mean(np.square(step)))
-		
-		fh.write("  Max Step = {:>14.9f}      RMS Step = {:>14.9f}\n\n".format(
-																		step_max, step_rms))
-		
-		return step
-
 class Atom:
 
 	def __init__(self, lbl,na,rx,ry,rz,chg,eps,sig):
diff --git a/DPpack/SetGlobals.py b/DPpack/SetGlobals.py
index e582316..6fd3c29 100644
--- a/DPpack/SetGlobals.py
+++ b/DPpack/SetGlobals.py
@@ -6,16 +6,13 @@ from DPpack.MolHandling import *
 from DPpack.PTable import *
 from DPpack.Misc import *
 
-env = ["OMP_STACKSIZE"]
-
-bohr2ang = 0.52917721092
-ang2bohr = 1/bohr2ang
-
 class Internal:
 
-	def __init__(self, infile):
+	def __init__(self, infile, outfile):
 
 		self.infile = infile
+		self.outfile = outfile
+
 		self.system = System()
 
 		self.player = self.Player()
@@ -45,10 +42,9 @@ class Internal:
 	def read_keywords(self):
 
 		try:
-			with open(self.infile) as fh:
-				controlfile = fh.readlines()
+			controlfile = self.infile.readlines()
 		except EnvironmentError:
-			sys.exit("Error: cannot open file {}".format(self.infile))
+			sys.exit("Error: cannot read file {}".format(self.infile))
 		
 		for line in controlfile:
 			
@@ -77,7 +73,7 @@ class Internal:
 				elif key in ('readhessian', 'vdwforces') and value[0].lower() in ("yes", "no"):
 					setattr(self.player, key, value[0].lower())
 				
-				elif key in ('maxcyc', 'initcyc', 'nprocs', 'altsteps', 'switchcyc'):
+				elif key in ('maxcyc', 'nprocs', 'altsteps', 'switchcyc'):
 					err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile)
 					try:
 						new_value = int(value[0])
@@ -95,7 +91,7 @@ class Internal:
 						if new_value < 0.01:
 							sys.exit(err)
 						else:
-							setattr(self.player, key).append(new_value)
+							setattr(self.player, key, new_value)
 					except ValueError:
 						sys.exit(err)
 						
@@ -136,7 +132,7 @@ class Internal:
 							if new_value < 1:
 								sys.exit(err)
 							else:
-								setattr(self.dice, key, new_value)
+								getattr(self.dice, key).append(new_value)
 						elif i == 0:
 							sys.exit(err)
 						else:
@@ -153,7 +149,7 @@ class Internal:
 							if new_value < 1:
 								sys.exit(err)
 							else:
-								setattr(self.dice, key, new_value)
+								getattr(self.dice, key).append(new_value)
 						elif i < 2:
 							sys.exit(err)
 						else:
@@ -179,7 +175,7 @@ class Internal:
 						sys.exit(err)
 					for i in range (2):
 						try:
-							setattr(self.gaussian, key)[i] = int(value[i])
+							getattr(self.gaussian, key)[i] = int(value[i])
 						except ValueError:
 							sys.exit(err)
 				
@@ -192,24 +188,24 @@ class Internal:
 				elif key == 'pop' and value[0].lower() in ("chelpg", "mk", "nbo"):
 					setattr(self.gaussian, key, value[0].lower())
 			
-			####  Read the Molcas related keywords
-			elif key in self.molcas_keywords and len(value) != 0:  ##  'value' is not empty!
+			# ####  Read the Molcas related keywords
+			# elif key in self.molcas_keywords and len(value) != 0:  ##  'value' is not empty!
 				
-				if key == 'root': # If defined, must be well defined (only positive integer values)
-					err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile)
-					if not value[0].isdigit():
-						sys.exit(err)
-					new_value = int(value[0])
-					if new_value >= 1:
-						setattr(self.molcas, key, new_value)
+			# 	if key == 'root': # If defined, must be well defined (only positive integer values)
+			# 		err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile)
+			# 		if not value[0].isdigit():
+			# 			sys.exit(err)
+			# 		new_value = int(value[0])
+			# 		if new_value >= 1:
+			# 			setattr(self.molcas, key, new_value)
 				
-				elif key in ('mbottom', 'orbfile'):
-					setattr(self.molcas, key, value[0])
+			# 	elif key in ('mbottom', 'orbfile'):
+			# 		setattr(self.molcas, key, value[0])
 				
-				elif key == 'basis':
-					setattr(self.molcas ,key, value[0])
+			# 	elif key == 'basis':
+			# 		setattr(self.molcas ,key, value[0])
 			
-			#### End
+			# #### End
 
 		
 	def check_keywords(self):
@@ -225,10 +221,10 @@ class Internal:
 		if self.dice.dens == None:
 			sys.exit("Error: 'dens' keyword not specified in file {}".format(self.infile))
 		
-		if len(self.dice.nmol) == 0:
+		if self.dice.nmol == 0:
 			sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(self.infile))
 		
-		if len(self.dice.nstep) == 0:
+		if self.dice.nstep == 0:
 			sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(self.infile))
 		
 		## Check only if QM program is Gaussian:
@@ -290,22 +286,22 @@ class Internal:
 		# isave value is always the nearest multiple of 100
 		self.dice.isave = round(self.dice.isave / 100) * 100
 
-	def print_keywords(self, fh):
+	def print_keywords(self):
 		
-		fh.write("##########################################################################################\n"
+		self.outfile.write("##########################################################################################\n"
 				"#############               Welcome to DICEPLAYER version 1.0                #############\n"
 				"##########################################################################################\n"
 				"\n")
-		fh.write("Your python version is {}\n".format(sys.version))
-		fh.write("\n")
-		fh.write("Program started on {}\n".format(weekday_date_time()))
-		fh.write("\n")
-		fh.write("Environment variables:\n")
+		self.outfile.write("Your python version is {}\n".format(sys.version))
+		self.outfile.write("\n")
+		self.outfile.write("Program started on {}\n".format(weekday_date_time()))
+		self.outfile.write("\n")
+		self.outfile.write("Environment variables:\n")
 		for var in env:
-			fh.write("{} = {}\n".format(var, 
+			self.outfile.write("{} = {}\n".format(var, 
 									(os.environ[var] if var in os.environ else "Not set")))
 		
-		fh.write("\n==========================================================================================\n"
+		self.outfile.write("\n==========================================================================================\n"
 				"                         CONTROL variables being used in this run:\n"
 				"------------------------------------------------------------------------------------------\n"
 				"\n")
@@ -314,13 +310,13 @@ class Internal:
 			if getattr(self.player,key) != None:
 				if isinstance(getattr(self.player,key), list):
 					string = " ".join(str(x) for x in getattr(self.player,key))
-					fh.write("{} = {}\n".format(key, string))
+					self.outfile.write("{} = {}\n".format(key, string))
 				else:	
-					fh.write("{} = {}\n".format(key, getattr(self.player,key)))
+					self.outfile.write("{} = {}\n".format(key, getattr(self.player,key)))
 		
-		fh.write("\n")
+		self.outfile.write("\n")
 
-		fh.write("------------------------------------------------------------------------------------------\n"
+		self.outfile.write("------------------------------------------------------------------------------------------\n"
 				"                         DICE variables being used in this run:\n"
 				"------------------------------------------------------------------------------------------\n"
 				"\n")
@@ -329,15 +325,15 @@ class Internal:
 			if getattr(self.dice,key) != None:
 				if isinstance(getattr(self.dice,key), list):
 					string = " ".join(str(x) for x in getattr(self.dice,key))
-					fh.write("{} = {}\n".format(key, string))
+					self.outfile.write("{} = {}\n".format(key, string))
 				else:	
-					fh.write("{} = {}\n".format(key, getattr(self.dice,key)))
+					self.outfile.write("{} = {}\n".format(key, getattr(self.dice,key)))
 		
-		fh.write("\n")
+		self.outfile.write("\n")
 		
 		if self.player.qmprog in ("g03", "g09", "g16"):
 
-			fh.write("------------------------------------------------------------------------------------------\n"
+			self.outfile.write("------------------------------------------------------------------------------------------\n"
 					"                         GAUSSIAN variables being used in this run:\n"
 					"------------------------------------------------------------------------------------------\n"
 					"\n")
@@ -346,15 +342,15 @@ class Internal:
 				if getattr(self.gaussian,key) != None:
 					if isinstance(getattr(self.gaussian,key), list):
 						string = " ".join(str(x) for x in getattr(self.gaussian,key))
-						fh.write("{} = {}\n".format(key, string))
+						self.outfile.write("{} = {}\n".format(key, string))
 					else:	
-						fh.write("{} = {}\n".format(key, getattr(self.gaussian,key)))
+						self.outfile.write("{} = {}\n".format(key, getattr(self.gaussian,key)))
 			
-			fh.write("\n")
+			self.outfile.write("\n")
 		
 		# elif self.player.qmprog == "molcas":
 
-		# 	fh.write("------------------------------------------------------------------------------------------\n"
+		# 	self.outfile.write("------------------------------------------------------------------------------------------\n"
 		# 			"                         MOLCAS variables being used in this run:\n"
 		# 			"------------------------------------------------------------------------------------------\n"
 		# 			"\n")
@@ -363,11 +359,11 @@ class Internal:
 		# 		if molcas[key] != None:
 		# 			if isinstance(molcas[key], list):
 		# 				string = " ".join(str(x) for x in molcas[key])
-		# 				fh.write("{} = {}\n".format(key, string))
+		# 				self.outfile.write("{} = {}\n".format(key, string))
 		# 			else:	
-		# 				fh.write("{} = {}\n".format(key, molcas[key]))
+		# 				self.outfile.write("{} = {}\n".format(key, molcas[key]))
 			
-		# 	fh.write("\n")
+		# 	self.outfile.write("\n")
 
 	def read_potential(self): # Deve ser atualizado para o uso de 
 		
@@ -390,7 +386,6 @@ class Internal:
 		if ntypes != len(self.dice.nmol):
 			sys.exit("Error: number of molecule types in file {} must match that of 'nmol' keyword in file {}".format(
 																	self.dice.ljname, self.infile))
-		
 		line = 2
 		for i in range(ntypes):
 
@@ -400,7 +395,7 @@ class Internal:
 				sys.exit("Error: expected an integer in line {} of file {}".format(line, self.dice.ljname))
 			
 			nsites = int(nsites)
-			self.system.add_type(Molecule())
+			self.system.add_type(nsites,Molecule())
 			
 			for j in range(nsites):
 
@@ -410,8 +405,6 @@ class Internal:
 				if len(new_atom) < 8:
 					sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, self.dice.ljname))
 				
-				self.system.molecule[i].add_atom()
-				
 				if not new_atom[0].isdigit():
 					sys.exit("Error: expected an integer in field 1, line {} of file {}".format(line, self.dice.ljname))
 				lbl = int(new_atom[0])
@@ -468,7 +461,7 @@ class Internal:
 							"Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format(
 																							line, self.dice.ljname))
 
-				self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig,mass))
+				self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig))
 				
 		to_delete = ['lbl','na','rx','ry','rz','chg','eps','sig','mass']
 		for _var in to_delete:
@@ -476,44 +469,44 @@ class Internal:
 				exec(f'del {_var}')
 
 	
-	def print_potential(self, fh):
+	def print_potential(self):
 		
 		formatstr = "{:<3d} {:>3d}  {:>10.5f} {:>10.5f} {:>10.5f}  {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n"
-		fh.write("\n"
+		self.outfile.write("\n"
 				"==========================================================================================\n")
-		fh.write("                    Potential parameters from file {}:\n".format(self.dice.ljname))
-		fh.write("------------------------------------------------------------------------------------------\n"
+		self.outfile.write("                    Potential parameters from file {}:\n".format(self.dice.ljname))
+		self.outfile.write("------------------------------------------------------------------------------------------\n"
 				"\n")
 		
-		fh.write("Combination rule: {}\n".format(self.dice.combrule))
-		fh.write("Types of molecules: {}\n\n".format(len(self.system.molecule)))
+		self.outfile.write("Combination rule: {}\n".format(self.dice.combrule))
+		self.outfile.write("Types of molecules: {}\n\n".format(len(self.system.molecule)))
 		
 		i = 0
 		for mol in self.system.molecule:
 			i += 1
-			fh.write("{} atoms in molecule type {}:\n".format(len(mol), i))
-			fh.write("---------------------------------------------------------------------------------\n"
+			self.outfile.write("{} atoms in molecule type {}:\n".format(len(mol.atom), i))
+			self.outfile.write("---------------------------------------------------------------------------------\n"
 					"Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass\n")
-			fh.write("---------------------------------------------------------------------------------\n")
+			self.outfile.write("---------------------------------------------------------------------------------\n")
 
 			for atom in mol.atom:
 
-				fh.write(formatstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz,
+				self.outfile.write(formatstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz,
 										atom.chg, atom.eps, atom.sig, atom.mass))
 			
-			fh.write("\n")
+			self.outfile.write("\n")
 			
 		if self.player.ghosts == "yes" or self.player.lps == "yes":
-			fh.write("\n"
+			self.outfile.write("\n"
 				"------------------------------------------------------------------------------------------\n"
 				"                    Aditional potential parameters:\n"
 				"------------------------------------------------------------------------------------------\n")
 		
 		# if player['ghosts'] == "yes":
 			
-		# 	fh.write("\n")
-		# 	fh.write("{} ghost atoms appended to molecule type 1 at:\n".format(len(ghost_types)))
-		# 	fh.write("---------------------------------------------------------------------------------\n")
+		# 	self.outfile.write("\n")
+		# 	self.outfile.write("{} ghost atoms appended to molecule type 1 at:\n".format(len(ghost_types)))
+		# 	self.outfile.write("---------------------------------------------------------------------------------\n")
 			
 		# 	atoms_string = ""
 		# 	for ghost in ghost_types:
@@ -522,19 +515,19 @@ class Internal:
 		# 			atoms_string += "{}{} ".format(atom_sym,atom)
 				
 		# 		if ghost['type'] == "g":
-		# 			fh.write(textwrap.fill("* Geometric center of atoms {}".format(atoms_string), 80))
+		# 			self.outfile.write(textwrap.fill("* Geometric center of atoms {}".format(atoms_string), 80))
 		# 		elif ghost['type'] == "m":
-		# 			fh.write(textwrap.fill("* Center of mass of atoms {}".format(atoms_string), 80))
+		# 			self.outfile.write(textwrap.fill("* Center of mass of atoms {}".format(atoms_string), 80))
 		# 		elif ghost['type'] == "z":
-		# 			fh.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80))
+		# 			self.outfile.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80))
 				
-		# 		fh.write("\n")
+		# 		self.outfile.write("\n")
 		
 		# if player['lps'] == 'yes':
 			
-		# 	fh.write("\n")
-		# 	fh.write("{} lone pairs appended to molecule type 1:\n".format(len(lp_types)))
-		# 	fh.write("---------------------------------------------------------------------------------\n")
+		# 	self.outfile.write("\n")
+		# 	self.outfile.write("{} lone pairs appended to molecule type 1:\n".format(len(lp_types)))
+		# 	self.outfile.write("---------------------------------------------------------------------------------\n")
 			
 		# 	for lp in lp_types:
 		# 		# LP type 1 or 2
@@ -546,36 +539,33 @@ class Internal:
 		# 			atom3_num = lp['numbers'][2]
 		# 			atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip()
 					
-		# 			fh.write(textwrap.fill(
+		# 			self.outfile.write(textwrap.fill(
 		# 			"* Type {} on atom {}{} with {}{} {}{}. Alpha = {:<5.1f} Deg and D = {:<4.2f} Angs".format(
 		# 				lp['type'], atom1_sym, atom1_num, atom2_sym, atom2_num, atom3_sym, atom3_num, lp['alpha'],
 		# 				lp['dist']), 86))
-		# 			fh.write("\n")
+		# 			self.outfile.write("\n")
 					
 		# 		# Other LP types
 				
-		fh.write("\n"
+		self.outfile.write("\n"
 				"==========================================================================================\n")
-		
-		return
 
-
-	def check_executables(self, fh):
+	def check_executables(self):
 	
-		fh.write("\n")
-		fh.write(90 * "=")
-		fh.write("\n\n")
+		self.outfile.write("\n")
+		self.outfile.write(90 * "=")
+		self.outfile.write("\n\n")
 		
 		dice_path = shutil.which(self.dice.progname)
 		if dice_path != None:
-			fh.write("Program {} found at {}\n".format(self.dice.progname, dice_path))
+			self.outfile.write("Program {} found at {}\n".format(self.dice.progname, dice_path))
 			self.dice.path = dice_path
 		else:
 			sys.exit("Error: cannot find dice executable")
 
 		qmprog_path = shutil.which(self.gaussian.qmprog)
 		if qmprog_path != None:
-			fh.write("Program {} found at {}\n".format(self.gaussian.qmprog, qmprog_path))
+			self.outfile.write("Program {} found at {}\n".format(self.gaussian.qmprog, qmprog_path))
 			self.gaussian.path = qmprog_path
 		else:
 			sys.exit("Error: cannot find {} executable".format(self.gaussian.qmprog))
@@ -583,16 +573,211 @@ class Internal:
 		if self.gaussian.qmprog in ("g03", "g09", "g16"):
 			formchk_path = shutil.which("formchk")
 			if formchk_path != None:
-				fh.write("Program formchk found at {}\n".format(formchk_path))
+				self.outfile.write("Program formchk found at {}\n".format(formchk_path))
 			else:
 				sys.exit("Error: cannot find formchk executable")
 
+
+	def calculate_step(self):
+		
+		invhessian = linalg.inv(self.system.molecule[0].hessian)
+		pre_step = -1 * np.matmul(invhessian, self.system.molecule[0].gradient.T).T
+		maxstep = np.amax(np.absolute(pre_step))
+		factor = min(1, self.player.maxstep/maxstep)
+		step = factor * pre_step
+		
+		self.player.outfile.write("\nCalculated step:\n")
+		pre_step_list = pre_step.tolist()
+		
+		self.player.outfile.write("-----------------------------------------------------------------------\n"
+				"Center     Atomic                          Step (Bohr)\n"
+				"Number     Number              X                Y                Z\n"
+				"-----------------------------------------------------------------------\n")
+		for i in range(len(self.system.molecule[0].atom)):
+			self.player.outfile.write("  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
+																i + 1, self.system.molecule[0].atom[i].na, 
+							pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0)))
+		
+		self.player.outfile.write("-----------------------------------------------------------------------\n")
+		
+		self.player.outfile.write("Maximum step is {:>11.6}\n".format(maxstep))
+		self.player.outfile.write("Scaling factor = {:>6.4f}\n".format(factor))
+		self.player.outfile.write("\nFinal step (Bohr):\n")
+		step_list = step.tolist()
+		
+		self.player.outfile.write("-----------------------------------------------------------------------\n"
+				"Center     Atomic                          Step (Bohr)\n"
+				"Number     Number              X                Y                Z\n"
+				"-----------------------------------------------------------------------\n")
+		for i in range(len(self.system.molecule[0].atom)):
+			self.player.outfile.write("  {:>5d}     {:>3d}        {:>14.9f}   {:>14.9f}   {:>14.9f}\n".format(
+																i + 1, self.system.molecule[0].atom[i].na, 
+										step_list.pop(0), step_list.pop(0), step_list.pop(0)))
+		
+		self.player.outfile.write("-----------------------------------------------------------------------\n")
+		
+		step_max = np.amax(np.absolute(step))
+		step_rms = np.sqrt(np.mean(np.square(step)))
+		
+		self.player.outfile.write("  Max Step = {:>14.9f}      RMS Step = {:>14.9f}\n\n".format(
+																		step_max, step_rms))
+		
+		return step
+
+	def read_initial_cicle(self):
+	
+		try:
+			with open(self.infile) as self.outfile:
+				controlfile = self.outfile.readlines()
+		except EnvironmentError:
+			sys.exit("Error: cannot open file {}".format(self.infile))
+		
+		for line in controlfile:
+
+			pass
+
+		
+	def populate_asec_vdw(self, cycle):
+			
+		asec_charges = []  	# (rx, ry, rz, chg)
+		vdw_meanfield = [] 	# (rx, ry, rz, eps, sig)
+		
+		if self.dice.nstep[-1] % self.dice.isave == 0:
+			nconfigs = round(self.dice.nstep[-1] / self.dice.isave)
+		else:
+			nconfigs = int(self.dice.nstep[-1] / self.dice.isave)
+		
+		norm_factor = nconfigs * self.player.nprocs
+		
+		nsitesref = len(self.system.molecule[0]) + len(ghost_atoms) + len(lp_atoms)
+		
+		nsites_total = dice['nmol'][0] * nsitesref
+		for i in range(1, len(dice['nmol'])):
+			nsites_total += dice['nmol'][i] * len(molecules[i])
+		
+		thickness = []
+		picked_mols = []
+			
+		for proc in range(1, player['nprocs'] + 1):  ## Run over folders
+			
+			path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
+			file = path + os.sep + dice['outname'] + ".xyz" 
+			if not os.path.isfile(file):
+				sys.exit("Error: cannot find file {}".format(file))
+			try:
+				with open(file) as xyzfh:
+					xyzfile = xyzfh.readlines()
+			except:
+				sys.exit("Error: cannot open file {}".format(file))
+			
+			for config in range(nconfigs):  ## Run over configs in a folder
+				
+				if int( xyzfile.pop(0).split()[0] ) != nsites_total:
+					sys.exit("Error: wrong number of sites in file {}".format(file))
+				
+				box = xyzfile.pop(0).split()[-3:]
+				box = [ float(box[0]), float(box[1]), float(box[2]) ]
+				sizes = sizes_of_molecule(molecules[0])
+				thickness.append( min([ (box[0] - sizes[0])/2, (box[1] - sizes[1])/2, 
+															(box[2] - sizes[2])/2 ]) )
+				
+				xyzfile = xyzfile[nsitesref:]  ## Skip the first (reference) molecule
+				mol_count = 0
+				for type in range(len(dice['nmol'])):  ## Run over types of molecules
+					
+					if type == 0:
+						nmols = dice['nmol'][0] - 1
+					else:
+						nmols = dice['nmol'][type]
+					
+					for mol in range(nmols):  ## Run over molecules of each type
+					
+						new_molecule = []
+						for site in range(len(molecules[type])):  ## Run over sites of each molecule
+						
+							new_molecule.append({})
+							line = xyzfile.pop(0).split()
+							
+							if line[0].title() != atomsymb[molecules[type][site]['na']].strip():
+								sys.exit("Error reading file {}".format(file))
+							
+							new_molecule[site]['na'] = molecules[type][site]['na']
+							new_molecule[site]['rx'] = float(line[1])
+							new_molecule[site]['ry'] = float(line[2])
+							new_molecule[site]['rz'] = float(line[3])
+							new_molecule[site]['chg'] = molecules[type][site]['chg']
+							new_molecule[site]['eps'] = molecules[type][site]['eps']
+							new_molecule[site]['sig'] = molecules[type][site]['sig']
+							
+						dist = minimum_distance(molecules[0], new_molecule)
+						if dist < thickness[-1]:
+							mol_count += 1
+							for atom in new_molecule:
+								asec_charges.append({})
+								vdw_meanfield.append({})
+								
+								asec_charges[-1]['rx'] = atom['rx']
+								asec_charges[-1]['ry'] = atom['ry']
+								asec_charges[-1]['rz'] = atom['rz']
+								asec_charges[-1]['chg'] = atom['chg'] / norm_factor
+								
+								if player['vdwforces'] == "yes":
+									vdw_meanfield[-1]['rx'] = atom['rx']
+									vdw_meanfield[-1]['ry'] = atom['ry']
+									vdw_meanfield[-1]['rz'] = atom['rz']
+									vdw_meanfield[-1]['eps'] = atom['eps']
+									vdw_meanfield[-1]['sig'] = atom['sig']
+						
+						####  Read lines with ghosts or lps in molecules of type 0 (reference)
+						####  and, if dist < thickness, appends to asec
+						if type == 0:
+							for ghost in ghost_atoms:
+								line = xyzfile.pop(0).split()
+								if line[0] != dice_ghost_label:
+									sys.exit("Error reading file {}".format(file))
+								if dist < thickness[-1]:
+									asec_charges.append({})
+									asec_charges[-1]['rx'] = float(line[1])
+									asec_charges[-1]['ry'] = float(line[2])
+									asec_charges[-1]['rz'] = float(line[3])
+									asec_charges[-1]['chg'] = ghost['chg'] / norm_factor
+							
+							for lp in lp_atoms:
+								line = xyzfile.pop(0).split()
+								if line[0] != dice_ghost_label:
+									sys.exit("Error reading file {}".format(file))
+								if dist < thickness[-1]:
+									asec_charges.append({})
+									asec_charges[-1]['rx'] = float(line[1])
+									asec_charges[-1]['ry'] = float(line[2])
+									asec_charges[-1]['rz'] = float(line[3])
+									asec_charges[-1]['chg'] = lp['chg'] / norm_factor
+				
+				picked_mols.append(mol_count)
+		
+		self.player.outfile.write("Done\n")
+		
+		string =  "In average, {:^7.2f} molecules ".format(sum(picked_mols)/norm_factor)
+		string += "were selected from each of the {} configurations ".format(len(picked_mols))
+		string += "of the production simulations to form the ASEC, comprising a shell with "
+		string += "minimum thickness of {:>6.2f} Angstrom\n".format(sum(thickness)/norm_factor)
+		
+		self.player.outfile.write(textwrap.fill(string, 86))
+		self.player.outfile.write("\n")
+		
+		otherfh = open("ASEC.dat", "w")
+		for charge in asec_charges:
+			otherfh.write("{:>10.5f}   {:>10.5f}   {:>10.5f}     {:>11.8f}\n".format(
+									charge['rx'], charge['ry'], charge['rz'], charge['chg']))
+		otherfh.close()
+		
+		return asec_charges
+
 	class Player:
 
 		def __init__(self):
 
 			self.maxcyc = None
-			# self.initcyc = 1 # Eliminated
 			self.nprocs = 1
 			self.switchcyc = 3
 			self.altsteps = 20000
@@ -604,9 +789,10 @@ class Internal:
 			self.ghosts = "no"
 			self.vdwforces = "no"
 			self.tol_factor = 1.2
-		
-
+			self.qmprog = "g16"
 
+			self.cyc = 1
+	
 	class Dice:
 
 		def __init__(self):
diff --git a/DPpack/__pycache__/Dice.cpython-38.pyc b/DPpack/__pycache__/Dice.cpython-38.pyc
new file mode 100644
index 0000000..b171a96
Binary files /dev/null and b/DPpack/__pycache__/Dice.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/Gaussian.cpython-38.pyc b/DPpack/__pycache__/Gaussian.cpython-38.pyc
new file mode 100644
index 0000000..e482efe
Binary files /dev/null and b/DPpack/__pycache__/Gaussian.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/Misc.cpython-38.pyc b/DPpack/__pycache__/Misc.cpython-38.pyc
new file mode 100644
index 0000000..b5da372
Binary files /dev/null and b/DPpack/__pycache__/Misc.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/MolHandling.cpython-38.pyc b/DPpack/__pycache__/MolHandling.cpython-38.pyc
new file mode 100644
index 0000000..95c9d61
Binary files /dev/null and b/DPpack/__pycache__/MolHandling.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/PTable.cpython-38.pyc b/DPpack/__pycache__/PTable.cpython-38.pyc
new file mode 100644
index 0000000..f0daee8
Binary files /dev/null and b/DPpack/__pycache__/PTable.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/SetGlobals.cpython-38.pyc b/DPpack/__pycache__/SetGlobals.cpython-38.pyc
new file mode 100644
index 0000000..00317ea
Binary files /dev/null and b/DPpack/__pycache__/SetGlobals.cpython-38.pyc differ
diff --git a/DPpack/__pycache__/__init__.cpython-38.pyc b/DPpack/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000..60f10af
Binary files /dev/null and b/DPpack/__pycache__/__init__.cpython-38.pyc differ
diff --git a/control.in b/control.in
index 98a38a1..36ed7f6 100644
--- a/control.in
+++ b/control.in
@@ -1,5 +1,5 @@
 # diceplayer
-initcyc = 1 
+initcyc = 1
 maxcyc = 3
 opt = NO
 nprocs = 2
@@ -10,12 +10,12 @@ altsteps = 20000
 
 # dice
 ncores = 1
-nmol = 1 100 
-dens = 0.75 
+nmol = 1 100
+dens = 0.75
 nstep = 40000 60000 50000
 isave = 1000
 ljname = phb.pot
 outname = phb
 
 # Gaussian
-level = MP2/aug-cc-pVTZ
+level = MP2/aug-cc-pVTZ
\ No newline at end of file
diff --git a/diceplayer.py b/diceplayer.py
index fb79f77..821fdc9 100644
--- a/diceplayer.py
+++ b/diceplayer.py
@@ -18,6 +18,7 @@ if __name__ == '__main__':
 ####  and set the usage and help messages                 ####
 
 	parser = argparse.ArgumentParser(prog='Diceplayer')
+	parser.add_argument('--continue', dest='opt_continue' , default=False, action='store_true')
 	parser.add_argument('--version', action='version', version='%(prog)s 1.0')
 	parser.add_argument('-i', dest='infile', default='control.in', metavar='INFILE', 
     			                   help='input file of diceplayer [default = control.in]')
@@ -27,270 +28,300 @@ if __name__ == '__main__':
 
 	args = parser.parse_args()
 
-####  Read and check the keywords in INFILE
-
-	read_keywords(args.infile)
-	check_keywords(args.infile)
-
 ####  Open OUTFILE for writing and print keywords and initial info
 
 	try:
-		if player['initcyc'] > 1 and os.path.exists(args.outfile):
-			oldname = args.outfile + ".old"
-			os.replace(args.outfile, oldname)
-		logfh = open(args.outfile, 'w', 1)
+
+		if args.opt_continue and os.path.exists(args.outfile):
+			
+			outfile = open(args.outfile,'r')
+			run_file = outfile.readlines()
+			control_sequence = '                                         Step # '
+			
+			for line in run_file:
+				if control_sequence in line:
+					cyc = int(line[-2]) + 1
+
+			outfile.close()
+			os.rename(os.path.abspath(args.outfile),os.path.abspath(args.outfile)+".backup")
+			outfile = open(args.outfile,'w')
+
+
+		if os.path.exists(args.outfile):		
+			os.rename(os.path.abspath(args.outfile),os.path.abspath(args.outfile)+".backup")
+			outfile = open(args.outfile,'w')
+		else:
+			outfile = open(args.outfile,"w")
+
 	except EnvironmentError as err:
 		sys.exit(err)
 
-	print_keywords(logfh)
+	try:
 
-####  Check whether the executables are in the path
+		if os.path.exists(args.infile):
+			infile = open(args.infile,"r")
 
-	check_executables(logfh)
+	except EnvironmentError as err:
+		sys.exit(err)
 
-####  Read the potential, store the info in 'molecules' and prints the info in OUTFILE
+####  Read and check the keywords in INFILE
 
-	read_potential(args.infile)
+	internal = Internal(infile, outfile)
 
-	if player['lps'] == "yes":
-		read_lps()
+	internal.read_keywords()
 
-	if player['ghosts'] == "yes":
-		read_ghosts()
+	if args.opt_continue:
+		internal.player.cyc = cyc
 
-	print_potential(logfh)
+	internal.check_keywords()
+	internal.print_keywords()
+
+# ####  Check whether the executables are in the path
+
+	# internal.check_executables()
+
+# ####  Read the potential, store the info in 'molecules' and prints the info in OUTFILE
+
+	internal.read_potential()
+
+	# if internal.player.lps == "yes":
+	# 	read_lps()
+
+	# if internal.player.ghosts == "yes":
+	# 	read_ghosts()
+
+	internal.print_potential()
 
 ####  Bring the molecules to standard orientation and prints info about them
 
-	for i in range(len(molecules)):
-		logfh.write("\nMolecule type {}:\n\n".format(i + 1))
-		print_mol_info(molecules[i], logfh)
-		logfh.write("    Translating and rotating molecule to standard orientation...")
-		standard_orientation(molecules[i])
-		logfh.write(" Done\n\n    New values:\n")
-		print_mol_info(molecules[i], logfh)
+	for i in range(len(internal.system.molecule)):
+		internal.outfile.write("\nMolecule type {}:\n\n".format(i + 1))
+		internal.system.molecule[i].print_mol_info(internal.outfile)
+		internal.outfile.write("    Translating and rotating molecule to standard orientation...")
+		internal.system.molecule[i].standard_orientation()
+		internal.outfile.write(" Done\n\n    New values:\n")
+		internal.system.molecule[i].print_mol_info(internal.outfile)
 	
-	logfh.write(90 * "=")
-	logfh.write("\n")
+	internal.outfile.write(90 * "=")
+	internal.outfile.write("\n")
 
 ####  Open the geoms.xyz file and prints the initial geometry if starting from zero
 
-	if player['initcyc'] == 1:
+	if internal.player.cyc == 1:
 		try:
 			geomsfh = open("geoms.xyz", "w", 1)
 		except EnvironmentError as err:
 			sys.exit(err)
-		print_geom(0, geomsfh)
+		internal.system.print_geom(0, geomsfh)
 	else:
 		try:
 			geomsfh = open("geoms.xyz", "A", 1)
 		except EnvironmentError as err:
 			sys.exit(err)
 
-
-	logfh.write("\nStarting the iterative process.\n")
+	# internal.outfile.write("\nStarting the iterative process.\n")
 	
-	## Initial position (in Bohr)
-	position = read_position(molecules[0])
+	# ## Initial position (in Bohr)
+	# position = internal.system.molecule[0].read_position()
 	
-	## If restarting, read the last gradient and hessian
-	if player['initcyc'] > 1:
-		if player['qmprog'] in ("g03", "g09", "g16"):
-			Gaussian.read_forces("grad_hessian.dat")
-			Gaussian.read_hessian_fchk("grad_hessian.dat")
+	# ## If restarting, read the last gradient and hessian
+	# if internal.player.cyc > 1:
+	# 	if internal.player.qmprog in ("g03", "g09", "g16"):
+	# 		Gaussian.read_forces("grad_hessian.dat")
+	# 		Gaussian.read_hessian_fchk("grad_hessian.dat")
 	
-		#if player['qmprog'] == "molcas":
-			#Molcas.read_forces("grad_hessian.dat")
-			#Molcas.read_hessian("grad_hessian.dat")
+	# 	#if player['qmprog'] == "molcas":
+	# 		#Molcas.read_forces("grad_hessian.dat")
+	# 		#Molcas.read_hessian("grad_hessian.dat")
 	
-	####
-	####  Start the iterative process
-	####
+	# ####
+	# ####  Start the iterative process
+	# ####
 	
-	for cycle in range(player['initcyc'], player['initcyc'] + player['maxcyc']):
+# 	for cycle in range(internal.player.cyc, internal.player.cyc + internal.player.maxcyc):
 	
-		logfh.write("\n" + 90 * "-" + "\n")
-		logfh.write("{} Step # {}\n".format(40 * " ", cycle))
-		logfh.write(90 * "-" + "\n\n")
+# 		internal.outfile.write("\n" + 90 * "-" + "\n")
+# 		internal.outfile.write("{} Step # {}\n".format(40 * " ", cycle))
+# 		internal.outfile.write(90 * "-" + "\n\n")
 	
-		make_step_dir(cycle)
+# 		make_step_dir(cycle)
 	
-		if player['altsteps'] == 0 or cycle == 1:
-			dice['randominit'] = True
-		else:
-			dice['randominit'] = False
+# 		if internal.player.altsteps == 0 or cycle == 1:
+# 			internal.dice.randominit = True
+# 		else:
+# 			internal.dice.randominit = False
 	
-		####
-		####  Start block of parallel simulations
-		####
+# 		####
+# 		####  Start block of parallel simulations
+# 		####
 		
-		procs = []
-		sentinels = []
-		for proc in range(1, player['nprocs'] + 1):
+# 		procs = []
+# 		sentinels = []
+# 		for proc in range(1, internal.player.nprocs + 1):
 		
-			p = Process(target=Dice.simulation_process, args=(cycle, proc, logfh))
-			p.start()
-			procs.append(p)
-			sentinels.append(p.sentinel)
+# 			p = Process(target=Dice.simulation_process, args=(cycle, proc, internal.outfile))
+# 			p.start()
+# 			procs.append(p)
+# 			sentinels.append(p.sentinel)
 			
-		while procs:
-			finished = connection.wait(sentinels)
-			for proc_sentinel in finished:
-				i = sentinels.index(proc_sentinel)
-				status = procs[i].exitcode
-				procs.pop(i)
-				sentinels.pop(i)
-				if status != 0:
-					for p in procs:
-						p.terminate()
-					sys.exit(status)
+# 		while procs:
+# 			finished = connection.wait(sentinels)
+# 			for proc_sentinel in finished:
+# 				i = sentinels.index(proc_sentinel)
+# 				status = procs[i].exitcode
+# 				procs.pop(i)
+# 				sentinels.pop(i)
+# 				if status != 0:
+# 					for p in procs:
+# 						p.terminate()
+# 					sys.exit(status)
 		
-		for proc in range(1, player['nprocs'] + 1):
-			Dice.print_last_config(cycle, proc)
+# 		for proc in range(1, internal.player.nprocs + 1):
+# 			Dice.print_last_config(cycle, proc)
 		
-		####
-		####  End of parallel simulations block
-		####	
+# 		####
+# 		####  End of parallel simulations block
+# 		####	
 		
-		## Make ASEC
-		logfh.write("\nBuilding the ASEC and vdW meanfields... ")
-		asec_charges = populate_asec_vdw(cycle, logfh)
+# 		## Make ASEC
+# 		internal.outfile.write("\nBuilding the ASEC and vdW meanfields... ")
+# 		asec_charges = internal.populate_asec_vdw(cycle)
 		
-		## After ASEC is built, compress files bigger than 1MB
-		for proc in range(1, player['nprocs'] + 1):
-			path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
-			compress_files_1mb(path)
+# 		## After ASEC is built, compress files bigger than 1MB
+# 		for proc in range(1, internal.player.nprocs + 1):
+# 			path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
+# 			compress_files_1mb(path)
 		
-		####
-		####  Start QM calculation
-		####
+# 		####
+# 		####  Start QM calculation
+# 		####
 		
-		make_qm_dir(cycle)
+# 		make_qm_dir(cycle)
 		
-		if player['opt'] == "yes":
+# 		if internal.player.opt == "yes":
 			
-			##
-			##  Gaussian block
-			##
-			if player['qmprog'] in ("g03", "g09", "g16"):
+# 			##
+# 			##  Gaussian block
+# 			##
+# 			if internal.player.qmprog in ("g03", "g09", "g16"):
 				
-				if cycle > 1:
-					src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
-					dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
-					shutil.copyfile(src, dst)
+# 				if cycle > 1:
+# 					src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
+# 					dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
+# 					shutil.copyfile(src, dst)
 				
-				Gaussian.make_force_input(cycle, asec_charges)
-				Gaussian.run_gaussian(cycle, "force", logfh)
-				Gaussian.run_formchk(cycle, logfh)
+# 				Gaussian.make_force_input(cycle, asec_charges)
+# 				Gaussian.run_gaussian(cycle, "force", internal.outfile)
+# 				Gaussian.run_formchk(cycle, internal.outfile)
 				
-				## Read the gradient
-				file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk"
-				gradient = Gaussian.read_forces(file, logfh)
-				if len(cur_gradient) > 0:
-					old_gradient = cur_gradient
-				cur_gradient = gradient
+# 				## Read the gradient
+# 				file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk"
+# 				gradient = Gaussian.read_forces(file, internal.outfile)
+# 				if len(cur_gradient) > 0:
+# 					old_gradient = cur_gradient
+# 				cur_gradient = gradient
 				
-				## If 1st step, read the hessian
-				if cycle == 1:
-					if player['readhessian'] == "yes":
-						file = "grad_hessian.dat"
-						logfh.write("\nReading the hessian matrix from file {}\n".format(file))
-						hessian = Gaussian.read_hessian_fchk(file)
-					else:
-						file = "step01" + os.sep + "qm" + os.sep + "asec.fchk"
-						logfh.write("\nReading the hessian matrix from file {}\n".format(file))
-						hessian = Gaussian.read_hessian(file)
+# 				## If 1st step, read the hessian
+# 				if cycle == 1:
+# 					if internal.player.readhessian == "yes":
+# 						file = "grad_hessian.dat"
+# 						internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file))
+# 						hessian = Gaussian.read_hessian_fchk(file)
+# 					else:
+# 						file = "step01" + os.sep + "qm" + os.sep + "asec.fchk"
+# 						internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file))
+# 						hessian = internal.gaussian.read_hessian(file)
 				
-				## From 2nd step on, update the hessian
-				else:
-					logfh.write("\nUpdating the hessian matrix using the BFGS method... ")
-					hessian = update_hessian(step, cur_gradient, old_gradient, hessian)
-					logfh.write("Done\n")
+# 				## From 2nd step on, update the hessian
+# 				else:
+# 					internal.outfile.write("\nUpdating the hessian matrix using the BFGS method... ")
+# 					hessian = internal.system.molecule[0].update_hessian(step, cur_gradient, old_gradient, hessian)
+# 					internal.outfile.write("Done\n")
 				
-				## Save gradient and hessian
-				Gaussian.print_grad_hessian(cycle, cur_gradient, hessian)
+# 				## Save gradient and hessian
+# 				internal.gaussian.print_grad_hessian(cycle, cur_gradient, hessian)
 				
-				## Calculate the step and update the position
-				step = calculate_step(cur_gradient, hessian, logfh)
-				position += step
+# 				## Calculate the step and update the position
+# 				step = internal.calculate_step(cur_gradient, hessian, internal.outfile)
+# 				position += step
 				
-				## Update the geometry of the reference molecule
-				update_molecule(position, logfh)
+# 				## Update the geometry of the reference molecule
+# 				internal.system.molecule[0].update_molecule(position, internal.outfile)
 				
-				## If needed, calculate the charges
-				if cycle < player['switchcyc']:
+# 				## If needed, calculate the charges
+# 				if cycle < internal.player.switchcyc:
 					
-					Gaussian.make_charge_input(cycle, asec_charges)
-					Gaussian.run_gaussian(cycle, "charge", logfh)
+# 					internal.gaussian.make_charge_input(cycle, asec_charges)
+# 					internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
 				
-				## Read the new charges and update molecules[0]
-				if cycle < player['switchcyc']:
-					file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
-					Gaussian.read_charges(file, logfh)
-				else:
-					file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log"
-					Gaussian.read_charges(file, logfh)
+# 				## Read the new charges and update molecules[0]
+# 				if cycle < internal.player.switchcyc:
+# 					file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
+# 					internal.gaussian.read_charges(file, internal.outfile)
+# 				else:
+# 					file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log"
+# 					internal.gaussian.read_charges(file, internal.outfile)
 				
-				## Print new info for molecule[0]
-				logfh.write("\nNew values for molecule type 1:\n\n")
-				print_mol_info(molecules[0], logfh)
+# 				## Print new info for molecule[0]
+# 				internal.outfile.write("\nNew values for molecule type 1:\n\n")
+# 				internal.system.molecule[0].print_mol_info()
 				
-				## Print new geometry in geoms.xyz
-				print_geom(cycle, geomsfh)
+# 				## Print new geometry in geoms.xyz
+# 				internal.system.molecule[0].print_geom(cycle, geomsfh)
 				
-			##
-			##  Molcas block
-			##
-			#if player['qmprog'] == "molcas":
+# 			##
+# 			##  Molcas block
+# 			##
+# 			#if player['qmprog'] == "molcas":
 		
 		
-		#elif player['opt'] == "ts":
+# 		#elif player['opt'] == "ts":
 			
-			##
-			##  Gaussian block
-			##
-			#if player['qmprog'] in ("g03", "g09", "g16"):
+# 			##
+# 			##  Gaussian block
+# 			##
+# 			#if player['qmprog'] in ("g03", "g09", "g16"):
 				
 				
 				
-			##
-			##  Molcas block
-			##
-			#if player['qmprog'] == "molcas":
+# 			##
+# 			##  Molcas block
+# 			##
+# 			#if player['qmprog'] == "molcas":
 				
 		
-		else:    ## Only relax the charge distribution
+# 		else:    ## Only relax the charge distribution
 			
-			if player['qmprog'] in ("g03", "g09", "g16"):
+# 			if internal.player.qmprog in ("g03", "g09", "g16"):
 				
-				if cycle > 1:
-					src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
-					dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
-					shutil.copyfile(src, dst)
+# 				if cycle > 1:
+# 					src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
+# 					dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
+# 					shutil.copyfile(src, dst)
 				
-				Gaussian.make_charge_input(cycle, asec_charges)
-				Gaussian.run_gaussian(cycle, "charge", logfh)
+# 				Gaussian.make_charge_input(cycle, asec_charges)
+# 				Gaussian.run_gaussian(cycle, "charge", internal.outfile)
 				
-				file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
-				Gaussian.read_charges(file)
+# 				file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
+# 				Gaussian.read_charges(file)
 				
-				## Print new info for molecule[0]
-				logfh.write("\nNew values for molecule type 1:\n\n")
-				print_mol_info(molecules[0], logfh)
+# 				## Print new info for molecule[0]
+# 				internal.outfile.write("\nNew values for molecule type 1:\n\n")
+# 				internal.system.molecule[0].print_mol_info()
 				
-			#if player['qmprog'] == "molcas":
+# 			#if player['qmprog'] == "molcas":
 				
 		
 	
-	####
-	####  End of the iterative process
-	####
+# 	####
+# 	####  End of the iterative process
+# 	####
 
-## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
-## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
+# ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
+# ## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
 
-	logfh.write("\nDiceplayer finished normally!\n")
-	logfh.close()
-####
-####  End of the program
-####
\ No newline at end of file
+# 	internal.outfile.write("\nDiceplayer finished normally!\n")
+# 	internal.outfile.close()
+# ####
+# ####  End of the program
+# ####
\ No newline at end of file
diff --git a/run.log b/run.log
index 0be58a8..45c1227 100644
--- a/run.log
+++ b/run.log
@@ -2,34 +2,38 @@
 #############               Welcome to DICEPLAYER version 1.0                #############
 ##########################################################################################
 
-Your python version is 3.6.1 (default, Apr 24 2017, 06:18:27) 
-[GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)]
+Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) 
+[GCC 7.3.0]
 
-Program started on Wednesday, 14 Jun 2017 at 12:30:02
+Program started on Saturday, 25 Sep 2021 at 15:24:31
 
 Environment variables:
-OMP_STACKSIZE = 32M
+OMP_STACKSIZE = Not set
 
 ==========================================================================================
                          CONTROL variables being used in this run:
 ------------------------------------------------------------------------------------------
 
 altsteps = 20000
+cyc = 1
+freq = no
 ghosts = no
-initcyc = 1
 lps = no
 maxcyc = 3
 maxstep = 0.3
 nprocs = 2
 opt = no
 qmprog = g09
-switch_step = 3
-zipprog = gzip
+readhessian = no
+switchcyc = 3
+tol_factor = 1.2
+vdwforces = no
 
 ------------------------------------------------------------------------------------------
                          DICE variables being used in this run:
 ------------------------------------------------------------------------------------------
 
+combrule = *
 dens = 0.75
 isave = 1000
 ljname = phb.pot
@@ -38,6 +42,7 @@ nmol = 1 100
 nstep = 40000 60000 50000
 outname = phb
 press = 1.0
+progname = dice
 temp = 300.0
 title = Diceplayer run
 
@@ -45,15 +50,13 @@ title = Diceplayer run
                          GAUSSIAN variables being used in this run:
 ------------------------------------------------------------------------------------------
 
+chglevel = MP2/aug-cc-pVTZ
 chgmult = 0 1
 level = MP2/aug-cc-pVTZ
 pop = chelpg
+qmprog = g09
 
 
-==========================================================================================
-
-Program dice found at /usr/local/bin/dice
-
 ==========================================================================================
                     Potential parameters from file phb.pot:
 ------------------------------------------------------------------------------------------
@@ -136,11 +139,11 @@ Molecule type 1:
     Translating and rotating molecule to standard orientation... Done
 
     New values:
-    Center of mass = (    -0.0000 ,     0.0000 ,    -0.0000 )
+    Center of mass = (    -0.0000 ,     0.0000 ,     0.0000 )
     Moments of inertia =  3.144054E+02  2.801666E+03  3.027366E+03
     Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
     Inter principal axis = (  -0.000000 ,   1.000000 ,  -0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
+    Minor principal axis = (  -0.000000 ,   0.000000 ,   1.000000 )
     Characteristic lengths = (  11.97 ,   5.27 ,   2.99 )
     Total mass =   226.28 au
     Total charge =  -0.0000 e
@@ -162,68 +165,14 @@ Molecule type 2:
     Translating and rotating molecule to standard orientation... Done
 
     New values:
-    Center of mass = (     0.0000 ,     0.0000 ,     0.0000 )
-    Moments of inertia =  1.205279E+02  2.859254E+02  4.033761E+02
-    Major principal axis = (  -1.000000 ,   0.000000 ,   0.000000 )
-    Inter principal axis = (   0.000000 ,   1.000000 ,   0.000000 )
-    Minor principal axis = (   0.000000 ,   0.000000 ,   1.000000 )
-    Characteristic lengths = (   5.67 ,   5.13 ,   1.51 )
+    Center of mass = (     1.6067 ,    -1.0929 ,     0.0026 )
+    Moments of inertia =  1.561638E+02  6.739391E+02  8.270247E+02
+    Major principal axis = (   0.899747 ,  -0.436411 ,   0.000541 )
+    Inter principal axis = (   0.436411 ,   0.899747 ,  -0.001219 )
+    Minor principal axis = (   0.000045 ,   0.001333 ,   0.999999 )
+    Characteristic lengths = (   5.67 ,   5.05 ,   1.51 )
     Total mass =   112.20 au
     Total charge =  -0.0000 e
-    Dipole moment = (    1.7575 ,    1.5369 ,   -0.0000 )     Total =    2.3347 Debye
+    Dipole moment = (    1.8148 ,    1.4689 ,   -0.0016 )     Total =    2.3347 Debye
 
 ==========================================================================================
-
-Starting the iterative process.
-
-------------------------------------------------------------------------------------------
-                                         Step # 1
-------------------------------------------------------------------------------------------
-
-Simulation process p01 initiated with pid 6822
-Simulation process p02 initiated with pid 6823
-p01> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02
-p02> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02
-p02> NPT thermalization initiated on 14 Jun 2017 at 12:41:16
-p01> NPT thermalization initiated on 14 Jun 2017 at 12:41:17
-p02> NPT production initiated on 14 Jun 2017 at 13:01:25
-p01> NPT production initiated on 14 Jun 2017 at 13:01:27
-p02> ----- NPT production finished on 14 Jun 2017 at 13:57:41
-p01> ----- NPT production finished on 14 Jun 2017 at 13:57:53
-
-Building the ASEC and vdW meanfields... In average, 99.97 molecules were selected from the production simulations to form the
-ASEC comprising a shell with minimum thickness of 15.352263400006278 AngstromDone.
-
-------------------------------------------------------------------------------------------
-                                         Step # 2
-------------------------------------------------------------------------------------------
-
-Simulation process p01 initiated with pid 6903
-Simulation process p02 initiated with pid 6904
-p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52
-p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52
-p02> NPT production initiated on 14 Jun 2017 at 14:05:14
-p01> NPT production initiated on 14 Jun 2017 at 14:05:20
-p02> ----- NPT production finished on 14 Jun 2017 at 14:23:58
-p01> ----- NPT production finished on 14 Jun 2017 at 14:24:05
-
-Building the ASEC and vdW meanfields... In average, 99.66 molecules were selected from the production simulations to form the
-ASEC comprising a shell with minimum thickness of 14.942583400006276 AngstromDone.
-
-------------------------------------------------------------------------------------------
-                                         Step # 3
-------------------------------------------------------------------------------------------
-
-Simulation process p01 initiated with pid 6945
-p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45
-Simulation process p02 initiated with pid 6946
-p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45
-p02> NPT production initiated on 14 Jun 2017 at 14:31:05
-p01> NPT production initiated on 14 Jun 2017 at 14:31:06
-p02> ----- NPT production finished on 14 Jun 2017 at 14:48:15
-p01> ----- NPT production finished on 14 Jun 2017 at 14:48:15
-
-Building the ASEC and vdW meanfields... In average, 99.0 molecules were selected from the production simulations to form the
-ASEC comprising a shell with minimum thickness of 14.79446440000628 AngstromDone.
-
-Diceplayer finished normally!
\ No newline at end of file
diff --git a/run.log.backup b/run.log.backup
new file mode 100644
index 0000000..1f9ac41
--- /dev/null
+++ b/run.log.backup
@@ -0,0 +1,180 @@
+##########################################################################################
+#############               Welcome to DICEPLAYER version 1.0                #############
+##########################################################################################
+
+Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) 
+[GCC 7.3.0]
+
+Program started on Saturday, 25 Sep 2021 at 15:23:44
+
+Environment variables:
+OMP_STACKSIZE = Not set
+
+==========================================================================================
+                         CONTROL variables being used in this run:
+------------------------------------------------------------------------------------------
+
+altsteps = 20000
+cyc = 1
+freq = no
+ghosts = no
+lps = no
+maxcyc = 3
+maxstep = 0.3
+nprocs = 2
+opt = no
+qmprog = g09
+readhessian = no
+switchcyc = 3
+tol_factor = 1.2
+vdwforces = no
+
+------------------------------------------------------------------------------------------
+                         DICE variables being used in this run:
+------------------------------------------------------------------------------------------
+
+combrule = *
+dens = 0.75
+isave = 1000
+ljname = phb.pot
+ncores = 1
+nmol = 1 100
+nstep = 40000 60000 50000
+outname = phb
+press = 1.0
+progname = dice
+temp = 300.0
+title = Diceplayer run
+
+------------------------------------------------------------------------------------------
+                         GAUSSIAN variables being used in this run:
+------------------------------------------------------------------------------------------
+
+chglevel = MP2/aug-cc-pVTZ
+chgmult = 0 1
+level = MP2/aug-cc-pVTZ
+pop = chelpg
+qmprog = g09
+
+
+==========================================================================================
+                    Potential parameters from file phb.pot:
+------------------------------------------------------------------------------------------
+
+Combination rule: *
+Types of molecules: 2
+
+31 atoms in molecule type 1:
+---------------------------------------------------------------------------------
+Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
+---------------------------------------------------------------------------------
+1     6    -4.40034    0.66551    0.41431   -0.387593   0.07000  3.5500   12.0110
+1     6    -4.49957   -0.67279   -0.15327    0.804001   0.07000  3.5500   12.0110
+1     6    -3.25630   -1.28922   -0.61351   -0.338491   0.07000  3.5500   12.0110
+1     6    -2.05849   -0.66807   -0.48674   -0.220783   0.07000  3.5500   12.0110
+1     6    -1.97114    0.65148    0.10511    0.563527   0.07000  3.5500   12.0110
+1     6    -3.20601    1.29684    0.49712   -0.145372   0.07000  3.5500   12.0110
+2     8    -5.59453   -1.25309   -0.28765   -0.714461   0.21000  2.9600   15.9990
+3     7    -0.87404    1.35870    0.27636   -0.721792   0.17000  3.2500   14.0070
+4     6     0.38785    0.81512    0.14410    0.607123   0.07000  3.5500   12.0110
+4     6     1.40776    1.65800   -0.33261   -0.295476   0.07000  3.5500   12.0110
+4     6     0.75950   -0.46815    0.58791   -0.346661   0.07000  3.5500   12.0110
+4     6     2.71021    1.21519   -0.46119   -0.298546   0.07000  3.5500   12.0110
+4     6     2.07354   -0.89884    0.50831   -0.234284   0.07000  3.5500   12.0110
+4     6     3.08076   -0.08936   -0.06112    0.294651   0.07000  3.5500   12.0110
+5     7     4.37238   -0.53979   -0.20183   -0.190497   0.17000  3.2500   14.0070
+6     6     5.41980    0.43824   -0.43836   -0.213183   0.06600  3.5000   12.0110
+6     6     4.76361   -1.73522    0.52225   -0.225420   0.06600  3.5000   12.0110
+7     1     6.36700   -0.08536   -0.51962    0.109840   0.03000  2.5000    1.0079
+7     1     5.25196    0.96612   -1.37421    0.077140   0.03000  2.5000    1.0079
+7     1     5.49517    1.17412    0.36749    0.125758   0.03000  2.5000    1.0079
+7     1     4.15838   -2.58442    0.21446    0.083229   0.03000  2.5000    1.0079
+7     1     5.79695   -1.96254    0.28041    0.115361   0.03000  2.5000    1.0079
+7     1     4.67416   -1.61707    1.60661    0.121440   0.03000  2.5000    1.0079
+7     1     0.02813   -1.09554    1.08007    0.179434   0.03000  2.4200    1.0079
+7     1     2.31876   -1.86576    0.91913    0.177250   0.03000  2.4200    1.0079
+7     1     3.45102    1.89880   -0.84669    0.198844   0.03000  2.4200    1.0079
+7     1     1.14759    2.66909   -0.61498    0.173388   0.03000  2.4200    1.0079
+7     1    -3.11930    2.29679    0.90060    0.151534   0.03000  2.4200    1.0079
+7     1    -5.31582    1.13611    0.74479    0.205289   0.03000  2.4200    1.0079
+7     1    -3.33813   -2.25037   -1.10231    0.184367   0.03000  2.4200    1.0079
+7     1    -1.16374   -1.11791   -0.89380    0.160382   0.03000  2.4200    1.0079
+
+16 atoms in molecule type 2:
+---------------------------------------------------------------------------------
+Lbl  AN       X          Y          Z         Charge    Epsilon   Sigma     Mass
+---------------------------------------------------------------------------------
+1     6     0.67203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
+1     6     2.07203   -2.82345    0.00263   -0.115000   0.07000  3.5500   12.0110
+1     6     2.76823   -1.61764    0.00263   -0.115000   0.07000  3.5500   12.0110
+1     6     2.06824   -0.40521    0.00264   -0.115000   0.07000  3.5500   12.0110
+1    14     0.67589   -0.40522    0.00264    0.150000   0.07000  3.5500   28.0860
+1   104    -0.02420   -1.61760    0.00263   -0.115000   0.07000  3.5500    0.0000
+2     1     0.13203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
+2     1     2.61203   -3.75875    0.00263    0.115000   0.03000  2.4200    1.0079
+2     1     2.60824    0.53010    0.00264    0.115000   0.03000  2.4200    1.0079
+2     1    -1.10420   -1.61760    0.00263    0.115000   0.03000  2.4200    1.0079
+3     8    -0.00411    0.77257    0.00264   -0.585000   0.17000  3.0700   15.9990
+3     1     0.61978    1.50220    0.00264    0.435000   0.00000  0.0000    1.0079
+4     6     4.27823   -1.61764    0.00263    0.115000   0.17000  3.8000   12.0110
+4     1     4.63490   -0.74399    0.50704    0.000000   0.00000  0.0000    1.0079
+4     1     4.63490   -2.49130    0.50704    0.000000   0.00000  0.0000    1.0079
+4     1     4.63490   -1.61764   -1.00617    0.000000   0.00000  0.0000    1.0079
+
+
+==========================================================================================
+
+Molecule type 1:
+
+    Center of mass = (    -0.0000 ,     0.0000 ,     0.0000 )
+    Moments of inertia =  3.144054E+02  2.801666E+03  3.027366E+03
+    Major principal axis = (   0.999972 ,   0.007210 ,   0.002184 )
+    Inter principal axis = (  -0.007218 ,   0.999967 ,   0.003660 )
+    Minor principal axis = (  -0.002157 ,  -0.003676 ,   0.999991 )
+    Characteristic lengths = (  11.96 ,   5.25 ,   2.98 )
+    Total mass =   226.28 au
+    Total charge =  -0.0000 e
+    Dipole moment = (    9.8367 ,    0.6848 ,    0.8358 )     Total =    9.8959 Debye
+
+    Translating and rotating molecule to standard orientation... Done
+
+    New values:
+    Center of mass = (    -0.0000 ,     0.0000 ,     0.0000 )
+    Moments of inertia =  3.144054E+02  2.801666E+03  3.027366E+03
+    Major principal axis = (   1.000000 ,   0.000000 ,   0.000000 )
+    Inter principal axis = (  -0.000000 ,   1.000000 ,  -0.000000 )
+    Minor principal axis = (  -0.000000 ,   0.000000 ,   1.000000 )
+    Characteristic lengths = (  11.97 ,   5.27 ,   2.99 )
+    Total mass =   226.28 au
+    Total charge =  -0.0000 e
+    Dipole moment = (    9.8432 ,    0.6168 ,    0.8120 )     Total =    9.8959 Debye
+
+
+Molecule type 2:
+
+    Center of mass = (     1.6067 ,    -1.0929 ,     0.0026 )
+    Moments of inertia =  1.205279E+02  2.859254E+02  4.033761E+02
+    Major principal axis = (   0.795913 ,  -0.605411 ,  -0.000001 )
+    Inter principal axis = (   0.605411 ,   0.795913 ,   0.000001 )
+    Minor principal axis = (   0.000000 ,  -0.000002 ,   1.000000 )
+    Characteristic lengths = (   5.74 ,   5.26 ,   1.51 )
+    Total mass =   112.20 au
+    Total charge =  -0.0000 e
+    Dipole moment = (    2.3293 ,    0.1593 ,   -0.0000 )     Total =    2.3347 Debye
+
+    Translating and rotating molecule to standard orientation... Done
+
+    New values:
+    Center of mass = (     1.6067 ,    -1.0929 ,     0.0026 )
+    Moments of inertia =  1.561638E+02  6.739391E+02  8.270247E+02
+    Major principal axis = (   0.899747 ,  -0.436411 ,   0.000541 )
+    Inter principal axis = (   0.436411 ,   0.899747 ,  -0.001219 )
+    Minor principal axis = (   0.000045 ,   0.001333 ,   0.999999 )
+    Characteristic lengths = (   5.67 ,   5.05 ,   1.51 )
+    Total mass =   112.20 au
+    Total charge =  -0.0000 e
+    Dipole moment = (    1.8148 ,    1.4689 ,   -0.0016 )     Total =    2.3347 Debye
+
+==========================================================================================
+
+Starting the iterative process.