Index: qmx.py
===================================================================
--- qmx.py	(revision 2)
+++ qmx.py	(revision 3)
@@ -8,43 +8,48 @@
 import numpy as np
 from general import Calculator
 from ase.embed import Embed
+from ase.units import Hartree
 
+from copy import deepcopy
+
 import sys, os
 
 class Qmx(Calculator):
-    def __init__(self, calculator_low, calculator_high):
-        self.string_params = {}
-        
-        self.forces=np.zeros((2,2))
+    def __init__(self, calculator_high,  calculator_low):
         self._constraints=None
-        
-        self.calculator_low = calculator_low
-        self.calculator_high = calculator_high
-                
+        self.calculator_low_cluster = calculator_low
+        self.calculator_low_system = deepcopy(calculator_low)
+        self.calculator_high = calculator_high            
+
     def get_energy_subsystem(self, path, calculator, atoms, force_consistent):
+        # go to directory and calculate energies
+        print "running energy in: ", path
         os.chdir(path)
-        calculator.set_atoms(atoms)
-        energy = calculator.get_potential_energy(atoms)
+        atoms.set_calculator(calculator)
+        energy = atoms.get_potential_energy()
         os.chdir("..")
         return energy
-        
+
     def get_forces_subsystem(self, path, calculator, atoms):
+        # go to directory and calculate forces
+        print "running forces in: ", path
         os.chdir(path)
-        calculator.set_atoms(atoms)
-        forces = calculator.get_forces(atoms)
+        atoms.set_calculator(calculator)
+        forces = atoms.get_forces()
         os.chdir("..")
         return forces
-        
-        
+
     def get_potential_energy(self, embed, force_consistent=False):
-        # perform energy calculations 
-        e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low, embed.get_system(), force_consistent)
-        e_cl_lo  = self.get_energy_subsystem("cluster.low-level", self.calculator_low, embed.get_cluster(), force_consistent)
+        # perform energy calculations
+        e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low_system, embed.get_system(), force_consistent)
+        e_cl_lo  = self.get_energy_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster(), force_consistent)
         e_cl_hi  = self.get_energy_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster(), force_consistent)
         # calculate energies
         energy = e_sys_lo - e_cl_lo + e_cl_hi
+        # print energies
         print "%20s=%15s - %15s + %15s" %("E(C:S)", "E(S-MM)", "E(C-MM)", "E(C-QM)")
         print "%20f=%15f - %15f + %15f" %(energy, e_sys_lo, e_cl_lo, e_cl_hi)
+        # set energies and return
         if force_consistent:
             self.energy_free = energy
             return self.energy_free
@@ -55,45 +60,38 @@
     def get_forces(self, embed):
         atom_map_sys_cl = embed.atom_map_sys_cl
         # get forces for the three systems
-        f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low, embed.get_system())
-        f_cl_lo  = self.get_forces_subsystem("cluster.low-level", self.calculator_low, embed.get_cluster())
+        f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low_system, embed.get_system())
+        f_cl_lo  = self.get_forces_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster())
         f_cl_hi  = self.get_forces_subsystem("cluster.high-level", self.calculator_high, embed.get_cluster())
         # forces correction for the atoms
         f_cl = f_cl_hi - f_cl_lo
-        #number of atoms
+        # number of atoms
         nat_sys = len(embed)
         # lo-sys + (hi-lo)
         for iat_sys in xrange(nat_sys):
             iat_cl = atom_map_sys_cl[iat_sys]
             if iat_cl > -1:
                 f_sys_lo[iat_sys] += f_cl[iat_cl]
-        
+
         # correct gradients
         # Reference: Eichler, Koelmel, Sauer, J. of Comput. Chem., 18(4). 1997, 463-477.
-        for cell_L, iat_cl, iat_sys, r, iat_link in embed.linkatoms:
+        for cell_L, iat_cl_sys, iat_sys, r, iat_link in embed.linkatoms:
             # calculate the bond distance (r_bond) at the border
-            xyz = embed[iat_sys].get_position() - embed.get_cluster()[iat_cl].get_position() + cell_L
+            xyz = embed[iat_sys].get_position() - embed[iat_cl_sys].get_position() + cell_L
             
             # calculate the bond lenght and the factor f
             rbond = np.sqrt(np.dot(xyz, xyz))
             f = r / rbond
             #normalize xyz
             xyz /= rbond
-            
             # receive the gradients for the link atom
-            fH = f_cl[iat_link]
-            # Skalarprodukt fH, xyz
-            fs = np.dot(xyz, fH)
-            
+            fL = f_cl[iat_link]
+            # dot product fL, xyz
+            fs = np.dot(xyz, fL)
+            # apply corrections for each direction
             for idir in xrange(3):
                 # correct the atom in the system
-                f_sys_lo[iat_sys][idir] += f*fH[idir] - f*fs*xyz[idir]
+                f_sys_lo[iat_sys][idir] += f*fL[idir] - f*fs*xyz[idir]
                 # correct the atom in the cluster
-                f_sys_lo[iat_cl][idir] += (1-f)*fH[idir] + f*fs*xyz[idir]
+                f_sys_lo[iat_cl_sys][idir] += (1-f)*fL[idir] + f*fs*xyz[idir]
         return f_sys_lo
-
-    def set_atoms(self, atoms):
-        return
-           
-            
-
