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

import numpy as np

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

    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)

        self.atoms = atoms

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

        lines = outfile.readlines()

        outfile.close()

        chargelines = 0

        if charges:

            nats = len(self.atoms)

            block = ''

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

        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)

                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

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