Chimie Théorique » QMX

class LJCalcDefinition(Definition):

    def __init__(self):

        self.name='Lennard Jones'

        self.system='LJ'

        self.keywords = {}

        self.keywords['import']='ase.calculators.lj'

        self.keywords['class']='LennardJones'

#-------------------------------------------------------------------------------

class EMTCalcDefinition(Definition):

    def __init__(self):

        self.name='EMT'

        self.system='EMT'

        self.keywords = {}

        self.keywords['import']='ase.calculators.emt'

        self.keywords['class']='EMT'

#-------------------------------------------------------------------------------

class SiestaCalcDefinition(Definition):

    def __init__(self):

        self.name='Siesta'

        self.system='Siesta'

        self.keywords = {}

        self.keywords['import']='ase.calculators.siesta'

        self.keywords['class']='Siesta'

#-------------------------------------------------------------------------------

class DacapoCalcDefinition(Definition):

    def __init__(self):

        self.name='Dacapo'

        self.system='Dacapo'

        self.keywords = {}

        self.keywords['import']='ase.calculators.dacapo'

        self.keywords['class']='Dacapo'

#-------------------------------------------------------------------------------

class AimsCubeCalcDefinition(Definition):

    def __init__(self):

        self.name='AimsCube'

        self.system='AimsCube'

        self.keywords = {}

        self.keywords['import']='ase.calculators.ase'

        self.keywords['class']='AimsCube'

#-------------------------------------------------------------------------------

class AimsCalcDefinition(Definition):

    def __init__(self):

        self.name='Aims'

        self.system='Aims'

        self.keywords = {}

        self.keywords['import']='ase.calculators.aims'

        self.keywords['class']='Aims'

#-------------------------------------------------------------------------------

class ExcitingCalcDefinition(Definition):

    def __init__(self):

        self.name='Exciting'

        self.system='Exciting'

        self.keywords = {}

        self.keywords['import']='ase.calculators.exciting'

        self.keywords['class']='Exciting'

#-------------------------------------------------------------------------------

class DftbCalcDefinition(Definition):

    def __init__(self):

        self.name='Dftb'

        self.system='Dftb'

        self.keywords = {}

        self.keywords['import']='ase.calculators.dftb'

        self.keywords['class']='Dftb'

#-------------------------------------------------------------------------------

class GaussianCalcDefinition(Definition):

    def __init__(self):

        self.name='Gaussian'

        self.system='Gaussian'

        self.keywords = {}

        self.keywords['import']='ase.calculators.gaussian'

        self.keywords['class']='Gaussian'

calcDefinitions = \

  [ VASPCalcDefinition(), TURBOMOLECalcDefinition(), MOPACCalcDefinition(), LJCalcDefinition(), \

    EMTCalcDefinition(), SiestaCalcDefinition(), DacapoCalcDefinition(), AimsCubeCalcDefinition(), \

    AimsCalcDefinition(), ExcitingCalcDefinition(), DftbCalcDefinition(), GaussianCalcDefinition() ]

class CLUSTERStrDefintion(Definition):

    def __init__(self):

        self.name='CLUSTER'

        self.keywords = {}

#-------------------------------------------------------------------------------

strSystemDefinitions = [VASPStrDefinition(), TURBOMOLEStrDefintion()]

strClusterDefinitions = [VASPStrDefinition(), TURBOMOLEStrDefintion(), CLUSTERStrDefintion()]

    for name in ['high-level', 'low-level', 'hybrid']:

    extraOptions = ''

    for name in ['cluster', 'system', 'embed']:

                if name is 'embed':

                    if strOptions is not '':

                        strOptions += ', '

                    strOptions += extraOptions

                stream.write(strSystem+'='+strClass+"("+strOptions+")\n")

            if name is 'cluster' and definition.name is 'CLUSTER':

                stream.write(strSystem+'=['+strOptions+']'+'\n')

                extraOptions='cluster_pos=False'

                strOptions=definition.getValue('set_calculator')

                stream.write(strSystem+'.set_calculator('+strOptions+')\n')

        self.keywords['set_calculator']='qmx'

from qmxSTR  import strSystemDefinitions, strClusterDefinitions

    if section == "cluster":

        definitions.append(analyzeSection(strClusterDefinitions, section, lines_section))

    #structures (atoms)

    if section == "system":

        definitions.append(analyzeSection(strSystemDefinitions, section, lines_section))

import tkFileDialog

from qmxSTR import strClusterDefinitions, strSystemDefinitions

    def __init__(self, mainWnd, master, title, name, def_list, fields):

            definition_new.system = name

        if fields is None:

            fields = ['']

            self.entries.append((key, text, Entry(self, textvariable=text, width=25)))

            Label(self, text=key, width=15, anchor=W).grid(row=irow,  column=2)

            textfield.grid(row=irow, column=3)

        self.listbox=Listbox(self, yscrollcommand=scrollbar.set, height=5, width=15, selectmode=SINGLE, exportselection=0)

        self.listbox.grid(row=0, column=0, rowspan=irow, sticky=NS)

        scrollbar.grid(row=0, column=1, rowspan=irow, sticky=NS)

        fields=['import', 'class', 'class.options']

        irow=0

        for calcLabel, strLabel, strDefinitions in (

                ("high-level",  "cluster", strClusterDefinitions), 

                ("low-level",  "system", strSystemDefinitions)):

            #column

            icolLocal=0;

            #build outFrame

            outFrame = LabelFrame(self, text=calcLabel+" method - "+strLabel)

            outFrame.grid(row=irow, columnspan=3,  sticky=W, pady=5)

            irow+=1

            for title, def_list in ((calcLabel, calcDefinitions),  (strLabel, strDefinitions)):

                # build frame

                myFrame=MyLabelFrame(self, outFrame, title, title, def_list, fields)

                myFrame.doLayout()

                myFrame.grid(row=0, column=icolLocal, padx=5, pady=5)

                # add frame for analyssis

                self.frames.append(myFrame)

                # move column

                icolLocal += 1

        outFrame = LabelFrame(self, text="General")

        outFrame.grid(row=irow, columnspan=3, pady=5); irow+=1

        for title, name, def_list, extra_list,  irowLocal, icolLocal,  irspan in (

            ("Optimization Method",  "job", jobDefinitions, fields + ['method', 'method.options'], 0, 0, 2), 

            ("Embedding Method", "embed", embedDefinitions, fields + ['set_calculator'], 0, 1, 1), 

            ("Hybrid Method", "hybrid", [QmxCalcDefinition()], fields, 1, 1, 1)):

            # build frame

            myFrame=MyLabelFrame(self, outFrame, title, name, def_list, extra_list)

            myFrame.doLayout()

            myFrame.grid(row=irowLocal, column=icolLocal,  rowspan=irspan, padx=5, pady=5)

            # add frame for analyssis

            self.frames.append(myFrame)

        outFrame=LabelFrame(self, text="PREVIEW")

        scrollbar=Scrollbar(outFrame, orient=VERTICAL)

        self.text=Text(outFrame, yscrollcommand=scrollbar.set, width=80, heigh=10, state=DISABLED)

        outFrame.grid(row=irow, column=0, sticky=W, pady=5);

        myFrame = Frame(self)

        buttonExit=Button(myFrame, text="Quit!", command=self.quit)

        buttonSave=Button(myFrame, text='save settings', state=DISABLED)

        buttonScript=Button(myFrame, text='Write SCRIPT file', command=self.writeScript)

        buttonPython=Button(myFrame, text='Write PYTHON file', command=self.writePython)

        myFrame.grid(row=irow, column=2, sticky=SE)

        fileName = tkFileDialog.asksaveasfile(mode='w')

        if fileName is None:

            return

        file = open(fileName,"w")

        fileName = tkFileDialog.asksaveasfile(mode='w')

        if fileName is None:

            return

        file = open(fileName,"w")

    def __init__(self, system, cluster, cell_cluster = "Auto", cluster_pos = True, linkType = 'H'):

        self.embed(linkType)

    def embed(self, linkType):

        self.set_linkatoms(linkType)

    def set_linkatoms(self, linkType, tol=15.):

        r_h = covalent_radii[atomic_numbers[linkType]]

            atom = Atom(symbol=linkType, position=xyz_diff, tag=50)

from ase.calculators.gaussian import Gaussian

                 'Aims', 'AimsCube', 'Turbomole', 'Exciting', 'Dftb', 'Qmx', 'Gaussian', 

        # writing output

        line = "running energy in: " + path + "\n"

        sys.stderr.write(line)

        # return path and result

        # writing output

        line = "running forces in: " + path + "\n"

        sys.stderr.write(line)

        # return path and result

        print "%20s = %15s - %15s + %15s" %("E(C:S)", "E(S,LL)", "E(C,LL)", "E(C,HL)")

    def get_stress(self, atoms):

    	return None

"""

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