Chimie Théorique » QMX

This module is the Gaussian module for ASE,  based on the PyMD Gaussian Module by R. Bulo,  

implemented by T. Kerber

 @author:       Rosa Bulo

 @organization: Vrije Universiteit Amsterdam

 @contact:      bulo@few.vu.nl

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)

import os,  sys,  string

from ase.units import Hartree,  Bohr

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

    functional:

        A string representing the DFT functional or wavefunction

        method to be used.

    """

    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)                

    def set_functional(self, functional):

        self.functional = functional

    def set_basis(self, basis):

        self.basis = basis

    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+'/')

        os.chdir(self.dir)

    def write_parfile(self, filename):

        parfile = open(filename, 'w')

        parfile.write(self.prm)

        parfile.close()

    def write_moldata(self,  type='force'):

        out = commands.getoutput('mkdir SETUP')

        out = commands.getoutput('ls SETUP/coord.*').split()

        pdb = False

        psf = False

        prm = False

        # 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()

        self.write_input(type)

    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')

        chargelines = 0

        if charges:

            nats = len(self.atoms)

            block = ''

            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()

        return energy

    def read_forces(self):

        factor = 1.0

        outfile = open('ASE-gaussian.log')

        lines = outfile.readlines()

        outfile.close()

        outputforces = None

        forces = None

            if re.search('Forces \(Ha', line):

                forces = np.zeros((nats,  3),  float)

                    forceline = lines[iline+iat+3]

                    words = forceline.split()

                    for idir,  word in enumerate(words[2:5]):

                        forces[iat,  idir] = float(word) * factor

        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()