Index: ase/calculators/gaussian.py
===================================================================
--- ase/calculators/gaussian.py	(revision 18)
+++ ase/calculators/gaussian.py	(revision 19)
@@ -1,251 +1,167 @@
 """
-This module is the Gaussian module for ASE,  based on the PyMD Gaussian Module by R. Bulo,  
-implemented by T. Kerber
+@author:       Torsten Kerber and Paul Fleurat-Lessard
+@organization: Ecole Normale Superieure de Lyon
+@contact:      Paul.Fleurat-Lessard@ens-lyon.fr
 
- @author:       Rosa Bulo
- @organization: Vrije Universiteit Amsterdam
- @contact:      bulo@few.vu.nl
+@data: June, 6, 2012
 
-Torsten Kerber,  ENS LYON: 2011,  07,  11
-
-This work is supported by Award No. UK-C0017,  made by King Abdullah
-University of Science and Technology(KAUST)
+This work is supported by Award No. UK-C0017 made by King Abdullah University of Science and Technology(KAUST)
 """
-import os,  sys,  string
-
+import os, sys
 import numpy as np
-from ase.units import Hartree,  Bohr
+
+from ase.units import Hartree, Bohr
 from ase.calculators.general import Calculator
-import os, sys, re, math, commands, subprocess, time
 
-class Gaussian(Calculator):
-    """
-    Class for representing the settings of a ReaxffForceJob
-        
-    It contains the following variables:
-                        
-    basis:  
-        A string representing the basisset to be used
+str_keys = [
+    'functional', # functional
+    'basis',      # basis set
+    'flags'       # more options
+    'memory',     # RAM memory used
+    'disk',       # disk space available
+    'filename',   # name of the gaussian files
+    'command'     # command to run Gaussian
+    ]
 
-    functional:
-        A string representing the DFT functional or wavefunction
-        method to be used.
+int_keys = [
+    'multiplicity', # multiplicity of the system
+    'nproc',        # number of processors
+    ]
 
-    """
-
-    def __init__(self, functional='BP86', basis='TZVP', atoms=None, inputfile=None):
-        self.functional = functional
-        self.basis = basis
-        self.restartdata = None
-        self.forces = None
-        self.atoms_old = None
-        if inputfile != None:
-            self = self.get_settings(inputfile)                
+class Gaussian(Calculator):
+    def __init__(self, atoms=None, **kwargs):
+        self.set_standard()
+        
+        self.set(**kwargs)
         self.atoms = atoms
+        
+        self.potential_energy = 0.0
+        self.forces = None
 
-    def set_functional(self, functional):
-        self.functional = functional
+        self.new_calculation = True
+        if atoms != None:
+            self.forces = np.zeros((len(atoms), 3))
+        self.stress = None
+        
+    def set_standard(self):
+        self.str_params = {}
+        self.int_params = {}
+        
+        for key in str_keys:
+            self.str_params['key'] = None
+        
+        for key in int_keys:
+            self.int_params['key'] = None
+        
+        self.str_params['functional'] = "B3LYP"
+        self.str_params['basis'] = '6-31G*'
+        self.str_params['flags'] = ' Force '
+        self.str_params['memory'] = '256MB'
+        self.str_params['disk'] = '16GB'
+        
+        self.str_params['command'] = 'g09'
+        self.str_params['filename'] = 'ASE-gaussian'
 
-    def set_basis(self, basis):
-        self.basis = basis
+        self.int_params['multiplicity'] = 1
+        self.int_params['nproc'] = 8
+        
+        
+    def set(self, **kwargs):
+        for key in kwargs:
+            if self.str_params.has_key(key):
+                self.str_params[key] = kwargs[key]
+            if self.int_params.has_key(key):
+                self.int_params[key] = kwargs[key]
+                
+                
+    def set_atoms(self, atoms):
+        if self.atoms != atoms:
+            self.atoms = atoms.copy()
+            self.new_calculation = True
+    
+    def write_file_line(self, file, line):
+        file.write(line)
+        file.write('\n')
+    
+    def write_input(self):
+        name = self.str_params['filename']
+        
+        inputfile = open(name + '.com', 'w')
+        self.write_file_line(inputfile, '%chk=' + name + '.chk')
+        self.write_file_line(inputfile, '%mem=' + self.str_params['memory'])
+        self.write_file_line(inputfile, '%nproc=' + ('%d' % self.int_params['nproc']))
+        self.write_file_line(inputfile, 'MaxDisk=' + self.str_params['disk'])
 
-    def setup_job(self,  type='force'):
-        ForceJob.setup_job(self, type=type)
- 
-        os.chdir(self.dirname)
-        self.myfiles = GaussianFileManager()
-        # Add the restart restuls to the filemanager. The checksum will be empty('')
-        if self.restartdata != None:
-            self.r_m = GaussianResults(self.myfiles)
-            self.myfiles.add_results(self.r_m, dir=self.restartdata+'/')
+        line = '#' + self.str_params['functional'] + '/' + self.str_params['basis'] + ' ' + self.str_params['flags']
+        self.write_file_line(inputfile, line)
+        
+        self.write_file_line(inputfile, '')
+        self.write_file_line(inputfile, 'Gaussian job created by ASE (supported by King Abdullah University of Science and Technology ,KAUST)')
+        self.write_file_line(inputfile, '')
 
-        os.chdir(self.dir)
+        line = '%d %d' % (sum(self.atoms.get_charges()), self.int_params['multiplicity'])
+        self.write_file_line(inputfile, line)
 
-    def write_parfile(self, filename):
-        parfile = open(filename, 'w')
-        parfile.write(self.prm)
-        parfile.close()
+        coords = self.atoms.get_positions()
+        for i, at in enumerate(self.atoms.get_chemical_symbols()):
+            line = at
+            coord = coords[i]
+            line += '  %16.5f %16.5f %16.5f'%(coord[0], coord[1], coord[2])
+            self.write_file_line(inputfile, line)
 
-    def write_moldata(self,  type='force'):
-        out = commands.getoutput('mkdir SETUP')
-        out = commands.getoutput('ls SETUP/coord.*').split()
-        pdb = False
-        psf = False
-        prm = False
+        if self.atoms.get_pbc().any(): 
+            cell = self.atoms.get_cell()
+            line = ''
+            for v in cell: 
+                line += 'TV %8.3f %8.3f %8.3f\n'%(v[0], v[1], v[2])
+            self.write_file_line(inputfile, line)
+            
+        self.write_file_line(inputfile, '')
+        
+    def update(self, atoms):
+        if (self.atoms is None or self.atoms.positions.shape != atoms.positions.shape):
+            self.atoms = atoms.copy()
+            self.calculate(atoms)
 
-        # Change labels variable if necessary
-        norderats = 0
-#        for el in self.labels:
-#            norderats += 1
-#            if norderats != self.atoms.pdb.get_number_of_atoms():
-#                self.labels = self.get_labels()
+    def calculate(self, atoms):
+        self.write_input()
+        line = self.str_params['command'] + ' ' + self.str_params['filename'] + '.com'
+        exitcode = os.system(line)
+        if exitcode != 0:
+            raise RuntimeError('Gaussian exited with exit code: %d.  ' % exitcode)
+        self.read_output()
+        self.new_calculation = False
+            
+    def get_potential_energy(self, atoms=None, force_consistent=False):
+        self.set_atoms(atoms)
+        if self.new_calculation:
+            self.calculate(atoms)
+        return self.potential_energy
 
-        self.write_input(type)
+    def get_forces(self, atoms):
+        self.set_atoms(atoms)
+        if self.new_calculation:
+            self.calculate(atoms)
+        return self.forces
 
-    def write_input(self, type):
-        input = self.get_input_block(type)
-        pyfile = open('ASE-gaussian.com', 'w')
-        pyfile.write(input)
-        pyfile.close()
-
-    def get_inp(self,  filename):
-        input = None
-        lis = commands.getoutput('ls %s'%(filename)).split('\n')
-        if not re.search('ls:', lis[0]):
-            infile = open(filename)
-            input = infile.readlines()
-            infile.close()
-        return input
-
-    def run(self, type='force'):
-        startrun = time.time()
-        # I have to write the input every timestep,  because that is where the coordinates are
-        # self.write_input(type)
-        self.write_moldata(type)
-
-        # Make sure that restart data is copied to the current directory
-#        if self.r_m != None:
-#            self.r_m.files.copy_job_result_files(self.r_m.fileid)
-
-        out = commands.getoutput('g09 ASE-gaussian.com')
-        endrun = time.time()
-        #print 'timings: ', endrun-startrun
-        # Write the output
-#        outputfile = open('ASE-gaussian.log')
-#        outfile = open('ASE-gaussian.out',  'w')
-#        for line in outputfile:
-#            outfile.write(line)
-#        outputfile.close()
-#        outfile.close()
-        # End write output
-        self.energy = self.read_energy()
-        if type == 'force':
-            self.forces = self.read_forces()
-        self.energy_zero = self.energy
-
-        # Change the results object to the new results
-        input = self.get_input_block(type)
-
-    def read_energy(self,  charges=True):
-        hartree2kcal = 627.5095
-
-        outfile = open('ASE-gaussian.log')
+    def read_output(self):
+        outfile = open(self.str_params['filename'] + '.log')
         lines = outfile.readlines()
         outfile.close()
 
-        chargelines = 0
-        if charges:
-            nats = len(self.atoms)
-            block = ''
+        factor = Hartree
         for line in lines:
-            if re.search('SCF Done', line):
-                words = line.split()
-                energy = float(words[4]) * hartree2kcal
-            if charges:
-                if re.search(' Mulliken atomic charges:', line):
-                    chargelines += 1
-            if chargelines > 0 and chargelines <= nats+2:
-                chargelines += 1
-                block += line
-   
-        if charges:
-            chargefile = open('charge.out', 'a')
-            chargefile.write(block)
-            chargefile.close()
+            if 'SCF Done' in line:
+                line = line.split()
+                self.potential_energy = float(line[4]) * factor
 
-        return energy
-
-    def read_forces(self):
         factor = Hartree / Bohr
-        factor = 1.0
-
-        outfile = open('ASE-gaussian.log')
-        lines = outfile.readlines()
-        outfile.close()
-
-        outputforces = None
-        forces = None
         nats = len(self.atoms)
-
         for iline, line in enumerate(lines):
-            if re.search('Forces \(Ha', line):
-                forces = np.zeros((nats,  3),  float)
+            if 'Forces (Hartrees/Bohr)' in line:
+                self.forces = np.zeros((nats, 3), float)
                 for iat in range(nats):
-                    forceline = lines[iline+iat+3]
-                    words = forceline.split()
-                    for idir,  word in enumerate(words[2:5]):
-                        forces[iat,  idir] = float(word) * factor
+                    line = lines[iline+iat+3].split()
+                    for idir, val in enumerate(line[2:5]):
+                        self.forces[iat, idir] = float(val) * factor
                 break
-        return forces
-
-    def get_input_block(self, type):
-        block = ''
-        block += '%chk=ASE-gaussian.chk\n'
-        block += '%Mem=256MB\n'
-        block += '%nproc=32\n'
-        block += 'MaxDisk=16gb\n'
-
-        block += '#'
-        block += '%s'%(self.functional)
-        block += '/'
-        block += '%s'%(self.basis)
-
-        if type == 'force':
-            block += 'Force '
-#        if self.r_m != None:
-#            block += 'Guess=Read'
-
-        block += '\n\nGaussian job\n\n'
-
-        charge = sum(self.atoms.get_charges())
-        block += '%i '%(charge)
-        block += '%i '%(1)
-        block += '\n'
-
-        coords = self.atoms.get_positions()
-        for i, at in enumerate(self.atoms.get_chemical_symbols()):
-            block += ' %2s'%(at)
-            coord = coords[i]
-            block += '  %16.5f %16.5f %16.5f'%(coord[0], coord[1], coord[2])
-            block += '\n'
-
-        if self.atoms.get_pbc().any(): 
-            cell = self.atoms.get_cell()
-           
-            for v in cell: 
-                block += 'TV %8.3f %8.3f %8.3f\n'%(v[0], v[1], v[2])
-        block += '\n'
-
-        return block
-
-
-    def get_settings(self,  filename=None):
-        settings = GaussianSettings()
-
-        if filename != None or filename == '':
-            input = self.get_inp(filename)
-        else:
-            settings.set_functional("")
-            return settings
-    
-        if input == None:
-            settings.set_functional("")
-            return settings
-
-        for i, line in enumerate(input):
-            if len(line.strip()) == 0:
-                continue
-            if line[0] == '#':
-                keyline = i
-                words = line[1:].split()
-            else:
-                settings.functional = words[0].split('/')[0]
-                settings.basis = words[0].split('/')[1]
-            break
-            return settings
-                
-    def update(self,  atoms):
-        if self.atoms_old != atoms:
-            self.atoms = atoms.copy()
-            self.atoms_old = atoms.copy()
-            self.run()
