Index: turbomole.py
===================================================================
--- turbomole.py	(revision 2)
+++ turbomole.py	(revision 3)
@@ -2,21 +2,18 @@
 
 http://www.turbomole.com/
 """
-
-
 import os, sys, string
 
-from ase.data import chemical_symbols
 from ase.units import Hartree, Bohr
 from ase.io.turbomole import write_turbomole
-from general import Calculator
+from ase.calculators.general import Calculator
 
 import numpy as np
 
 class Turbomole(Calculator):
     """Class for doing Turbomole calculations.
     """
-    def __init__(self, label='turbomole'):
+    def __init__(self, label='turbomole', calculate_energy="dscf", calculate_forces="grad", post_HF = False):
         """Construct TURBOMOLE-calculator object.
 
         Parameters
@@ -59,56 +56,28 @@
         e = a.get_potential_energy()
         
         """
-        
-        # get names of executables for turbomole energy and forces
-        # get the path
-
-        os.system('rm -f sysname.file; sysname > sysname.file')
-        f = open('sysname.file')
-        architechture = f.readline()[:-1]
-        f.close()
-        tmpath = os.environ['TURBODIR']
-        pre = tmpath + '/bin/' + architechture + '/'
-
-
-        if os.path.isfile('control'):
-            f = open('control')
-        else:
-            print 'File control is missing'
-            raise RuntimeError, \
-                'Please run Turbomole define and come thereafter back'
-        lines = f.readlines()
-        f.close()
-        self.tm_program_energy=pre+'dscf'
-        self.tm_program_forces=pre+'grad'        
-        for line in lines:
-            if line.startswith('$ricore'):
-                self.tm_program_energy=pre+'ridft'
-                self.tm_program_forces=pre+'rdgrad'
-
+       
         self.label = label
         self.converged = False
-        #clean up turbomole energy file
-        os.system('rm -f energy; touch energy')
+        
+        # set calculators for energy and forces
+        self.calculate_energy = calculate_energy
+        self.calculate_forces = calculate_forces
 
-        #turbomole has no stress
+        # turbomole has no stress
         self.stress = np.empty((3, 3))
         
+        # storage for energy and forces
+        self.e_total = None
+        self.forces = None
+        self.updated = False
+        
+        # atoms must be set
+        self.atoms = None
+        
+        # POST-HF method 
+        self.post_HF  = post_HF
 
-    def update(self, atoms):
-        """Energy and forces are calculated when atoms have moved
-        by calling self.calculate
-        """
-        if (not self.converged or
-            len(self.numbers) != len(atoms) or
-            (self.numbers != atoms.get_atomic_numbers()).any()):
-            self.initialize(atoms)
-            self.calculate(atoms)
-        elif ((self.positions != atoms.get_positions()).any() or
-              (self.pbc != atoms.get_pbc()).any() or
-              (self.cell != atoms.get_cell()).any()):
-            self.calculate(atoms)
-
     def initialize(self, atoms):
         self.numbers = atoms.get_atomic_numbers().copy()
         self.species = []
@@ -116,55 +85,66 @@
             self.species.append(Z)
         self.converged = False
         
+    def execute(self, command):
+        from subprocess import Popen, PIPE
+        try:
+            proc = Popen([command], shell=True, stderr=PIPE)
+            error = proc.communicate()[1]
+            if "abnormally" in error:
+                raise OSError(error)
+            print "TM command: ", command, "successfully executed"
+        except OSError, e:
+            print >>sys.stderr, "Execution failed:", e
+            sys.exit(1)
+
     def get_potential_energy(self, atoms):
-        self.update(atoms)
-        return self.etotal
+        # update atoms
+        self.set_atoms(atoms)
+        # if update of energy is neccessary
+        if self.update_energy:
+            # calculate energy
+            self.execute(self.calculate_energy + " > ASE.TM.energy.out")
+            # check for convergence of dscf cycle
+            if os.path.isfile('dscf_problem'):
+                print 'Turbomole scf energy calculation did not converge'
+                raise RuntimeError, \
+                    'Please run Turbomole define and come thereafter back'
+            # read energy
+            self.read_energy()
+        else :
+            print "taking old values (E)"
+        self.update_energy = False
+        return self.e_total
 
     def get_forces(self, atoms):
-        self.update(atoms)
+        # update atoms
+        self.set_atoms(atoms)
+        # complete energy calculations
+        if self.update_energy:
+            self.get_potential_energy(atoms)
+        # if update of forces is neccessary
+        if self.update_forces:
+            # calculate forces
+            self.execute(self.calculate_forces + " > ASE.TM.forces.out")
+            # read forces
+            self.read_forces()
+        else :
+            print "taking old values (F)"
+        self.update_forces = False
         return self.forces.copy()
     
     def get_stress(self, atoms):
-        self.update(atoms)
         return self.stress.copy()
-
-    def calculate(self, atoms):
-        """Total Turbomole energy is calculated (to file 'energy'
-        also forces are calculated (to file 'gradient')
-        """
-        self.positions = atoms.get_positions().copy()
-        self.cell = atoms.get_cell().copy()
-        self.pbc = atoms.get_pbc().copy()
-
-        #write current coordinates to file 'coord' for Turbomole
+        
+    def set_atoms(self, atoms):
+        if self.atoms == atoms:
+            return
+        # performs an update of the atoms 
+        super(Turbomole, self).set_atoms(atoms)
         write_turbomole('coord', atoms)
-
-
-        #Turbomole energy calculation
-        os.system('rm -f output.energy.dummy; \
-                      '+ self.tm_program_energy +'> output.energy.dummy')
-
-        #check that the energy run converged
-        if os.path.isfile('dscf_problem'):
-            print 'Turbomole scf energy calculation did not converge'
-            print 'issue command t2x -c > last.xyz'
-            print 'and check geometry last.xyz and job.xxx or statistics'
-            raise RuntimeError, \
-                'Please run Turbomole define and come thereafter back'
-
-
-        self.read_energy()
-
-        #Turbomole atomic forces calculation
-        #killing the previous gradient file because 
-        #turbomole gradients are affected by the previous values
-        os.system('rm -f gradient; rm -f output.forces.dummy; \
-                      '+ self.tm_program_forces +'> output.forces.dummy')
-
-        self.read_forces(atoms)
-
-        self.converged = True
-
+        # energy and forces must be re-calculated
+        self.update_energy = True
+        self.update_forces = True
         
     def read_energy(self):
         """Read Energy from Turbomole energy file."""
@@ -178,40 +158,39 @@
             elif line.startswith('$'):
                 pass
             else:
-                #print 'LINE',line
+                # print 'LINE',line
                 energy_tmp = float(line.split()[1])
-        #print 'energy_tmp',energy_tmp
-        self.etotal = energy_tmp * Hartree
+                if self.post_HF:
+                    energy_tmp += float(line.split()[4])
+        self.e_total = energy_tmp * Hartree
+        
 
-
-    def read_forces(self,atoms):
+    def read_forces(self):
         """Read Forces from Turbomole gradient file."""
-
         file = open('gradient','r')
-        line=file.readline()
-        line=file.readline()
-        tmpforces = np.array([[0, 0, 0]])
-        while line:
-            if 'cycle' in line:
-                for i, dummy in enumerate(atoms):
-                            line=file.readline()
-                forces_tmp=[]
-                for i, dummy in enumerate(atoms):
-                            line=file.readline()
-                            line2=string.replace(line,'D','E')
-                            #tmp=np.append(forces_tmp,np.array\
-                            #      ([[float(f) for f in line2.split()[0:3]]]))
-                            tmp=np.array\
-                                ([[float(f) for f in line2.split()[0:3]]])
-                            tmpforces=np.concatenate((tmpforces,tmp))  
-            line=file.readline()
-            
+        lines=file.readlines()
+        file.close()
 
-        #note the '-' sign for turbomole, to get forces
-        self.forces = (-np.delete(tmpforces, np.s_[0:1], axis=0))*Hartree/Bohr
+        forces = np.array([[0, 0, 0]])
+        
+        nline = len(lines)
+        iline = -1
+        
+        for i in xrange(nline):
+            if 'cycle' in lines[i]:
+                iline = i
+        
+        if iline < 0:
+            print 'Gradients not found'
+            raise RuntimeError("Please check TURBOMOLE gradients")
 
-        #print 'forces', self.forces
-
-    def read(self):
-        """Dummy stress for turbomole"""
-        self.stress = np.empty((3, 3))
+        iline += len(self.atoms) + 1
+        nline -= 1
+        for i in xrange(iline, nline):
+            line = string.replace(lines[i],'D','E')
+            #tmp=np.append(forces_tmp,np.array\
+            #      ([[float(f) for f in line.split()[0:3]]]))
+            tmp=np.array([[float(f) for f in line.split()[0:3]]])
+            forces=np.concatenate((forces,tmp))  
+        #note the '-' sign for turbomole, to get forces
+        self.forces = (-np.delete(forces, np.s_[0:1], axis=0))*Hartree/Bohr
