Index: ase/__init__.py
===================================================================
--- ase/__init__.py	(revision 1)
+++ ase/__init__.py	(revision 2)
@@ -6,6 +6,7 @@
 
 from ase.atom import Atom
 from ase.atoms import Atoms
+from ase.embed import Embed
 
 _deprecate_things_from_ase_module = True
 
Index: ase/optimize/optimize.py
===================================================================
--- ase/optimize/optimize.py	(revision 1)
+++ ase/optimize/optimize.py	(revision 2)
@@ -9,7 +9,7 @@
 import numpy as np
 
 from ase.parallel import rank, barrier
-from ase.io.trajectory import PickleTrajectory
+#from ase.io.trajectory import PickleTrajectory
 
 
 class Dynamics:
Index: ase/svnversion.py
===================================================================
--- ase/svnversion.py	(revision 1)
+++ ase/svnversion.py	(revision 2)
@@ -1 +1 @@
-svnversion = "1765"
+svnversion = "4M"
Index: ase/embed.py
===================================================================
--- ase/embed.py	(revision 0)
+++ ase/embed.py	(revision 2)
@@ -0,0 +1,252 @@
+from ase import Atom, Atoms
+from ase.data import covalent_radii, atomic_numbers
+
+import numpy as np
+
+class Embed(Atoms):
+    #--- constructor of the Embed class ---
+    def __init__(self, system=None, cluster=None):
+        super(Embed, self).__init__()
+        # define the atom map
+        self.atom_map_sys_cl = []
+        self.linkatoms = []
+        # cluster dimensions
+        self.xyz_cl_min = None
+        self.xyz_cl_max = None
+        # set the search radius for link atoms
+        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 ---        
+    def set_cluster(self, atoms):
+        import copy
+        # set the min/max cluster dimensions
+        self.xyz_cl_min = atoms[0].get_position()
+        self.xyz_cl_max = atoms[0].get_position()
+        for atom in atoms:
+            # assign the label "Cluster (10)" in atom.TAG
+            atom.set_tag(10)
+            xyz=atom.get_position()
+            for i in xrange(3):
+                # set the min/max cluster dimensions
+                if xyz[i] < self.xyz_cl_min[i]:
+                    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]       
+        # 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        
+        # assign the label "Cluster (10)" in atom.TAG
+        for atom in atoms:
+            atom.set_tag(0)
+        # update search radius for link atoms
+        dx = 0
+        for atom in atoms:
+            r = covalent_radii[atom.get_atomic_number()]
+            if (r > dx):
+                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 ?
+        if self.atoms_cluster is None or self.atoms_system is None:
+            return
+        self.find_cluster()
+        self.set_linkatoms()
+        print "link atoms found: ", len(self.linkatoms)
+
+    def find_cluster(self):
+        # set tolerance
+        d = 0.001
+        #atoms
+        xyzs_cl=[]
+        for atom_cl in self.atoms_cluster:
+            xyzs_cl.append(atom_cl.get_position())
+        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                
+        for iat_sys in xrange(len(self.atoms_system)):
+            # get the coordinates of the system atom atom_sys
+            xyz_sys = xyzs_sys[iat_sys]
+            # bSysOnly: no identical atom has been found
+            bSysOnly = True
+            # loop over system atoms atom_cl
+            for iat_cl in xrange(len(self.atoms_cluster)):
+                # difference vector between both atoms
+                xyz_diff = np.abs(xyzs_sys[iat_sys]-xyzs_cl[iat_cl])
+                # identical atoms
+                if xyz_diff[0] < d and xyz_diff[1] < d and xyz_diff[2] < d:
+                    # 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
+                    # atom has been identified
+                    bSysOnly = False
+                    break
+            if bSysOnly:
+                self.atom_map_sys_cl[iat_sys] = -1
+
+    def set_linkatoms(self, tol=15., linkAtom=None, debug=False):
+        # local copies of xyz coordinates to avoid massive copying of xyz objects
+        xyzs_cl=[]
+        for atom_cl in self.atoms_cluster:
+            xyzs_cl.append(atom_cl.get_position())
+        xyzs_sys=[]
+        for atom_sys in self.atoms_system:
+            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:
+            linkAtom ='H'
+        # number of atoms in the cluster and the system
+        nat_cl=len(self.atoms_cluster)
+        nat_sys=len(self.atoms_system)
+        # system has 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
+        cells_L = [(0.0,  0.0,  0.0)]
+        # get the cell vectors
+        cell = self.get_cell()
+        # TODO: check 0D, 1D, 2D, 3D
+        if False:
+         for ix in xrange(-1, 2):
+            for iy in xrange(-1, 2):
+                for iz in xrange(-1, 2):
+                    if ix == 0 and iy == 0 and iz == 0:
+                        continue
+
+                    cells_L.append(np.dot([ix, iy, iz], cell))
+        # save the radius of system atoms
+        rs_sys = []
+        for atom in self.atoms_system:
+            rs_sys.append(covalent_radii[atom.get_atomic_number()])
+        # sum over cluster atoms (iat_cl)
+        for iat_cl in xrange(nat_cl):
+            # get the cluster atom
+            atom_cl=self.atoms_cluster[iat_cl]
+            # ignore link atoms
+            if atom_cl.get_tag() == 50:
+                continue
+            # 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
+                if self.atoms_system[iat_sys].get_tag()==10:
+                    continue
+                # sum over all cell_L
+                for cell_L in cells_L:
+                    # xyz coordinates and covalent radius of the system atom iat_sys
+                    xyz_sys = xyzs_sys[iat_sys]+cell_L
+                    # 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 
+                            xyz_sys[i] > self.xyz_cl_max[i]):
+                            lcont = False
+                            break
+                    if not lcont:
+                        continue
+                    # xyz coordinates and covalent radius of the system atom iat_sys
+                    r_sys = rs_sys[iat_sys]
+                    # diff vector
+                    xyz_diff = xyz_sys - xyz_cl
+                    # distance between the atoms
+                    r = np.sqrt(np.dot(xyz_diff, xyz_diff))
+                    # ratio of the distance to the sum of covalent radius
+                    f = r / (r_cl + r_sys)
+                    if debug:
+                        print "Covalent radii = ",r_cl, r_sys
+                        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
+                        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
+            # 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]
+            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]
+            # 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
+            self.linkatoms.append(s)
+        return
+        
+    def set_positions(self, positions_new):
+        # number of atoms
+        nat_sys=len(self.atoms_system)
+        # go over all pairs of atoms
+        for iat_sys in xrange(nat_sys):
+            xyz = positions_new[iat_sys]
+            iat_cl = self.atom_map_sys_cl[iat_sys]
+            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:
+            # determine position of the link atom
+            xyz_cl = self.atoms_cluster[iat_cl].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__()
Index: ase/atoms.py
===================================================================
--- ase/atoms.py	(revision 1)
+++ ase/atoms.py	(revision 2)
@@ -285,7 +285,6 @@
         """Add new array.
 
         If *shape* is not *None*, the shape of *a* will be checked."""
-        
         if dtype is not None:
             a = np.array(a, dtype)
         else:
@@ -311,6 +310,7 @@
 
         Returns a copy unless the optional argument copy is false.
         """
+
         if copy:
             return self.arrays[name].copy()
         else:
Index: ase/calculators/__init__.py
===================================================================
--- ase/calculators/__init__.py	(revision 1)
+++ ase/calculators/__init__.py	(revision 2)
@@ -14,12 +14,13 @@
 from ase.calculators.dftb import Dftb
 from ase.calculators.singlepoint import SinglePointCalculator
 from ase.calculators.test import numeric_force, numeric_forces, TestPotential
+from ase.calculators.qmx import Qmx
 
 if _deprecate_things_from_ase_module:
     from ase.utils.deprecate import Deprecate
     _locals = locals()
     for name in ['LennardJones', 'EMT', 'Siesta', 'Dacapo', 'Vasp',
-                 'Aims', 'AimsCube', 'Turbomole', 'Exciting', 'Dftb',
+                 'Aims', 'AimsCube', 'Turbomole', 'Exciting', 'Dftb', 'Qmx', 
                  'SinglePointCalculator', 'numeric_force', 'numeric_forces',
                  'TestPotential']:
         obj = _locals[name]
Index: ase/calculators/qmx.py
===================================================================
--- ase/calculators/qmx.py	(revision 0)
+++ ase/calculators/qmx.py	(revision 2)
@@ -0,0 +1,99 @@
+""" This is a QM:MM embedded system for ASE
+
+torsten.kerber@ens-lyon.fr
+"""
+
+import ase
+import ase.atoms
+import numpy as np
+from general import Calculator
+from ase.embed import Embed
+
+import sys, os
+
+class Qmx(Calculator):
+    def __init__(self, calculator_low, calculator_high):
+        self.string_params = {}
+        
+        self.forces=np.zeros((2,2))
+        self._constraints=None
+        
+        self.calculator_low = calculator_low
+        self.calculator_high = calculator_high
+                
+    def get_energy_subsystem(self, path, calculator, atoms, force_consistent):
+        os.chdir(path)
+        calculator.set_atoms(atoms)
+        energy = calculator.get_potential_energy(atoms)
+        os.chdir("..")
+        return energy
+        
+    def get_forces_subsystem(self, path, calculator, atoms):
+        os.chdir(path)
+        calculator.set_atoms(atoms)
+        forces = calculator.get_forces(atoms)
+        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)
+        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 "%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)
+        if force_consistent:
+            self.energy_free = energy
+            return self.energy_free
+        else:
+            self.energy_zero = energy
+            return self.energy_zero
+
+    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_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
+        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:
+            # calculate the bond distance (r_bond) at the border
+            xyz = embed[iat_sys].get_position() - embed.get_cluster()[iat_cl].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)
+            
+            for idir in xrange(3):
+                # correct the atom in the system
+                f_sys_lo[iat_sys][idir] += f*fH[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]
+        return f_sys_lo
+
+    def set_atoms(self, atoms):
+        return
+           
+            
+
Index: ase/calculators/turbomole.py
===================================================================
--- ase/calculators/turbomole.py	(revision 1)
+++ ase/calculators/turbomole.py	(revision 2)
@@ -9,10 +9,11 @@
 from ase.data import chemical_symbols
 from ase.units import Hartree, Bohr
 from ase.io.turbomole import write_turbomole
+from general import Calculator
 
 import numpy as np
 
-class Turbomole:
+class Turbomole(Calculator):
     """Class for doing Turbomole calculations.
     """
     def __init__(self, label='turbomole'):
Index: ase/calculators/vasp.py
===================================================================
--- ase/calculators/vasp.py	(revision 1)
+++ ase/calculators/vasp.py	(revision 2)
@@ -149,8 +149,9 @@
 ]
 
 class Vasp(Calculator):
-    def __init__(self, restart=None, output_template='vasp', track_output=False, 
+    def __init__(self, restart=None, output_template='vasp', track_output=False, write_input = True, 
                  **kwargs):
+        self.bool_write_input = write_input
         self.name = 'Vasp'
         self.float_params = {}
         self.exp_params = {}
@@ -358,9 +359,10 @@
         from ase.io.vasp import write_vasp
         self.initialize(atoms)
         write_vasp('POSCAR', self.atoms_sorted, symbol_count = self.symbol_count)
-        self.write_incar(atoms)
-        self.write_potcar()
-        self.write_kpoints()
+        if self.bool_write_input:
+            self.write_incar(atoms)
+            self.write_potcar()
+            self.write_kpoints()
         self.write_sort_file()
         
         # Execute VASP
