Index: embed.py
===================================================================
--- embed.py	(revision 2)
+++ embed.py	(revision 3)
@@ -9,20 +9,21 @@
         super(Embed, self).__init__()
         # define the atom map
         self.atom_map_sys_cl = []
+        self.atom_map_cl_sys = []
         self.linkatoms = []
         # cluster dimensions
         self.xyz_cl_min = None
         self.xyz_cl_max = None
         # set the search radius for link atoms
-        self.d = 10        
+        self.d = 10
         # define the systems for calculations
         self.set_system(system)
         self.set_cluster(cluster)
         #
         self.set_cell([10, 10, 10])
         return
-                
-    #--- set the cluster ---        
+
+    #--- set the cluster ---
     def set_cluster(self, atoms):
         import copy
         # set the min/max cluster dimensions
@@ -38,17 +39,17 @@
                     self.xyz_cl_min[i] = xyz[i]
                 if xyz[i] > self.xyz_cl_max[i]:
                     self.xyz_cl_max[i] = xyz[i]
-                
+
         # add self.d around min/max cluster dimensions
         self.xyz_cl_min -= [self.d, self.d, self.d]
-        self.xyz_cl_max += [self.d, self.d, self.d]       
+        self.xyz_cl_max += [self.d, self.d, self.d]
         # set the cluster for low and high level calculation
         self.atoms_cluster = atoms
         return
 
     #--- set the system ---
     def set_system(self, atoms):
-        self.atoms_system = atoms        
+        self.atoms_system = atoms
         # assign the label "Cluster (10)" in atom.TAG
         for atom in atoms:
             atom.set_tag(0)
@@ -60,14 +61,14 @@
                 dx = r
         self.d = dx * 2.1
         return
-        
+
     #--- return cluster ---
     def get_cluster(self):
         return self.atoms_cluster
-        
+
     def get_system(self):
         return self.atoms_system
-        
+
     #--- Embedding ---
     def embed(self):
         # is the cluster and the host system definied ?
@@ -87,9 +88,10 @@
         xyzs_sys=[]
         for atom_sys in self.atoms_system:
             xyzs_sys.append(atom_sys.get_position())
-        
-        self.atom_map_sys_cl=np.zeros(len(self.atoms_system), int)
-        # loop over cluster atoms atom_sys                
+
+        self.atom_map_sys_cl = np.zeros(len(self.atoms_system), int)
+        self.atom_map_cl_sys = np.zeros(len(self.atoms_cluster), int)
+        # loop over cluster atoms atom_sys
         for iat_sys in xrange(len(self.atoms_system)):
             # get the coordinates of the system atom atom_sys
             xyz_sys = xyzs_sys[iat_sys]
@@ -104,7 +106,8 @@
                     # set tag (CLUSTER+HOST: 10) to atom_sys
                     self.atoms_system[iat_sys].set_tag(10)
                     # map the atom in the atom list
-                    self.atom_map_sys_cl[iat_sys]=iat_cl
+                    self.atom_map_sys_cl[iat_sys] = iat_cl
+                    self.atom_map_cl_sys[iat_cl] = iat_sys
                     # atom has been identified
                     bSysOnly = False
                     break
@@ -118,7 +121,7 @@
             xyzs_cl.append(atom_cl.get_position())
         xyzs_sys=[]
         for atom_sys in self.atoms_system:
-            xyzs_sys.append(atom_sys.get_position())        
+            xyzs_sys.append(atom_sys.get_position())
         # FIXME: mapping does not work when cluster atoms are outside of the box
         # set the standard link atom
         if linkAtom is None:
@@ -127,7 +130,7 @@
         nat_cl=len(self.atoms_cluster)
         nat_sys=len(self.atoms_system)
         # system has pbc?
-        pbc = self.get_pbc()        
+        pbc = self.get_pbc()
         # set the bond table
         bonds = []
         # set the 27 cell_vec, starting with the (0,0,0) vector for the unit cell
@@ -157,7 +160,7 @@
             # xyz coordinates and covalent radius of the cluster atom iat_cl
             xyz_cl = xyzs_cl[iat_cl]
             r_cl = covalent_radii[atom_cl.get_atomic_number()]
-            
+
             # sum over system atoms (iat_sys)
             for iat_sys in xrange(nat_sys):
                 # avoid cluster atoms
@@ -170,7 +173,7 @@
                     # go only in distance self.d around the cluster
                     lcont = True
                     for i in xrange(3):
-                        if (xyz_sys[i] < self.xyz_cl_min[i] or 
+                        if (xyz_sys[i] < self.xyz_cl_min[i] or
                             xyz_sys[i] > self.xyz_cl_max[i]):
                             lcont = False
                             break
@@ -189,33 +192,33 @@
                         print "Distance ", f
                         print "tol = ",(1+tol/100.),(1-tol/100.),(1-2*tol/100.)
                     if f <= (1+tol/100.) and f >= (1-2*tol/100.):
-                        s = cell_L, iat_cl, iat_sys, r_cl
+                        s = cell_L, self.atom_map_cl_sys[iat_cl], iat_sys, r_cl
                         bonds.append(s)
                         break
                     if f <= (1-2*tol/100.):
                         raise RuntimeError("QMX: The cluster atom", iat_cl, " and the system atom", iat_sys, "came too close")
-        
+
         r_h = covalent_radii[atomic_numbers[linkAtom]]
         for bond in bonds:
-            cell_L, iat_cl, iat_sys, r_cl = bond
+            cell_L, iat_cl_sys, iat_sys, r_cl = bond
             # assign the tags for the border atoms
             atom_sys.set_tag(1)
             atom_cl.set_tag(11)
             #difference vector for the link atom, scaling
-            xyz_diff = xyzs_sys[iat_sys] + cell_L - xyzs_cl[iat_cl]
+            xyz_diff = xyzs_sys[iat_sys] + cell_L - xyzs_sys[iat_cl_sys]
             r = (r_cl + r_h)
             xyz_diff *=  r / np.sqrt(np.dot(xyz_diff, xyz_diff))
             # determine position of the link atom
-            xyz_diff += xyzs_cl[iat_cl]
+            xyz_diff += xyzs_sys[iat_cl_sys]
             # create link atom
             atom = Atom(symbol=linkAtom, position=xyz_diff, tag=50)
             # add atom to cluster
             self.atoms_cluster.append(atom)
             # add atom to the linkatoms
-            s = cell_L, iat_cl, iat_sys, r, len(self.atoms_cluster)-1
+            s = cell_L, iat_cl_sys, iat_sys, r, len(self.atoms_cluster)-1
             self.linkatoms.append(s)
         return
-        
+
     def set_positions(self, positions_new):
         # number of atoms
         nat_sys=len(self.atoms_system)
@@ -226,27 +229,30 @@
             self.atoms_system[iat_sys].set_position(xyz)
             if iat_cl > -1:
                 self.atoms_cluster[iat_cl].set_position(xyz)
-                
-        for cell_L, iat_cl, iat_sys, r, iat in self.linkatoms:
+
+        for cell_L, iat_cl_sys, iat_sys, r, iat in self.linkatoms:
             # determine position of the link atom
-            xyz_cl = self.atoms_cluster[iat_cl].get_position()
+            xyz_cl = self.atoms_system[iat_cl_sys].get_position()
             xyz = self.atoms_system[iat_sys].get_position() - xyz_cl + cell_L
             xyz *= r / np.sqrt(np.dot(xyz, xyz))
             xyz += xyz_cl
             # update xyz coordinates of the cluster
             self.atoms_cluster[iat].set_position(xyz)
-            
+
     def __getitem__(self, i):
         return self.atoms_system.__getitem__(i)
 
     def get_positions(self):
         return self.atoms_system.get_positions()
-        
+
     def __add__(self, other):
         return self.atoms_system.__add__(other)
-        
+
     def __delitem__(self, i):
         return self.atoms_system.__delitem__(i)
-    
+
     def __len__(self):
         return self.atoms_system.__len__()
+        
+    def get_chemical_symbols(self, reduce=False):
+        return self.atoms_system.get_chemical_symbols(reduce)
