psi4-data 1:0.3-5 / usr / share / psi / python / qcdb / molpro.py

#
#@BEGIN LICENSE
#
# PSI4: an ab initio quantum chemistry software package
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
#@END LICENSE
#

import math

import qcdb.exceptions
import qcformat
import molpro_basissets


class MolproIn(qcformat.InputFormat):

    def __init__(self, mem, mtd, bas, mol, sys, cast):
        qcformat.InputFormat.__init__(self, mem, mtd, bas, mol, sys, cast)

        # memory in MB --> MW
        self.memory = int(math.ceil(mem / 8.0))
        # auxiliary basis sets
        [self.unaugbasis, self.augbasis, self.auxbasis] = self.corresponding_aux_basis()

    def format_global_parameters(self):
        text = ''

        if self.method in ['mp2c', 'dft-sapt-shift', 'dft-sapt', 'dft-sapt-pbe0ac', 'dft-sapt-pbe0acalda']:
            text += """GTHRESH,ZERO=1.e-14,ONEINT=1.e-14,TWOINT=1.e-14,ENERGY=1.e-8,ORBITAL=1.e-8,GRID=1.e-8\n\n"""
        elif self.method in ['b3lyp', 'b3lyp-d', 'df-b3lyp', 'df-b3lyp-d']:
            text += """GTHRESH,ZERO=1.e-14,ONEINT=1.e-14,TWOINT=1.e-14,ENERGY=1.e-8,ORBITAL=1.e-7,GRID=1.e-8\n\n"""
        else:
            text += """GTHRESH,ZERO=1.e-14,ONEINT=1.e-14,TWOINT=1.e-14,ENERGY=1.e-9\n\n"""

        return text

    def format_basis(self):
        text = ''
        text += """basis={\n"""

        try:
            # jaxz, maxz, etc.
            for line in molpro_basissets.altbasis[self.basis]:
                text += """%s\n""" % (line)
            text += '\n'
        except KeyError:
            # haxz
            if self.basis.startswith('heavy-aug-'):
                text += """set,orbital; default,%s,H=%s\n""" % (self.basis[6:], self.unaugbasis)
            # xz, axz, 6-31g*
            else:
                text += """set,orbital; default,%s\n""" % (self.basis)

        if ('df-' in self.method) or ('f12' in self.method) or (self.method in ['mp2c', 'dft-sapt', 'dft-sapt-pbe0acalda']):
            if self.unaugbasis and self.auxbasis:

                text += """set,jkfit;   default,%s/jkfit\n""" % (self.auxbasis)
                text += """set,jkfitb;  default,%s/jkfit\n""" % (self.unaugbasis)
                text += """set,mp2fit;  default,%s/mp2fit\n""" % (self.auxbasis)
                text += """set,dflhf;   default,%s/jkfit\n""" % (self.auxbasis)
            else:
                raise qcdb.exceptions.ValidationError("""Auxiliary basis not predictable from orbital basis '%s'""" % (self.basis))

        text += """}\n\n"""
        return text

    def format_infile_string(self):
        text = ''

        # format comment and memory
        text += """***, %s %s\n""" % (self.index, self.molecule.tagline)
        text += """memory,%d,m\n""" % (self.memory)

        # format molecule, incl. charges and dummy atoms
        text += self.molecule.format_molecule_for_molpro()

        # format global convergence directions
        text += self.format_global_parameters()

        # format castup directions
        if self.castup is True:
            text += """basis=sto-3g\n"""
            text += """rhf\n"""
            text += '\n'

        # format basis set
        text += self.format_basis()

        # format method
        for line in qcmtdIN[self.method]:
            text += """%s\n""" % (line)
        text += """show[1,20f20.12],ee*,ce*,te*\n"""
        text += """show[1,60f20.12],_E*\n"""
        text += '\n'

        return text


qcmtdIN = {
'ccsd(t)-f12': [
    'rhf',
    'eehf=energy',
    'ccsd(t)-f12,df_basis=mp2fit,df_basis_exch=jkfitb,ri_basis=jkfitb',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'eemp3=emp3',
    'cemp3=eemp3-eehf',
    'eeccsd=energc',
    'ceccsd=eeccsd-eehf',
    'eeccsdt=energy',
    'ceccsdt=eeccsdt-eehf',
    'temp2=emp2_trip',
    'teccsd=ectrip'],

'ccsd(t)': [
    'rhf',
    'eehf=energy',
    'ccsd(t)',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'eemp3=emp3',
    'cemp3=eemp3-eehf',
    'eeccsd=energc',
    'ceccsd=eeccsd-eehf',
    'eeccsdt=energy',
    'ceccsdt=eeccsdt-eehf',
    'temp2=emp2_trip',
    'teccsd=ectrip'],

'mp3': [
    'gdirect',
    'rhf',
    'eehf=energy',
    'mp3',
    'eemp2=emp2',
    'eemp3=emp3',
    'eemp25=0.5*(eemp2+eemp3)',
    'cemp2=eemp2-eehf',
    'cemp3=eemp3-eehf',
    'cemp25=eemp25-eehf',
    'temp2=emp2_trip',
    'temp3=ectrip'],

'mp2': [
    'gdirect',
    'rhf',
    'eehf=energy',
    'mp2',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'df-hf-mp2': [
    'gdirect',
    '{df-hf,basis=jkfit}',
    'eehf=energy',
    'mp2',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'hf-df-mp2': [
    'gdirect',
    'rhf',
    'eehf=energy',
    '{df-mp2,basis_mp2=mp2fit}',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'hf': [
    'rhf',
    'eehf=energy'],

'mp2-f12': [
    'gdirect',
    'rhf',
    'eehf=energy',
    'mp2-f12',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'df-mp2-f12': [
    'gdirect',
    #'rhf',
    '{df-hf,basis=jkfit}',
    'eehf=energy',
    #'{df-mp2-f12,df_basis=mp2fit,df_basis_exch=jkfit,ri_basis=optrib}',
    '{df-mp2-f12,df_basis=mp2fit,df_basis_exch=jkfitb,ri_basis=jkfitb}',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'df-mp2': [
    'gdirect',
    '{df-hf,basis=jkfit}',
    'eehf=energy',
    '{df-mp2,basis_mp2=mp2fit}',
    'eemp2=emp2',
    'cemp2=eemp2-eehf',
    'temp2=emp2_trip'],

'df-hf': [
    'gdirect',
    '{df-hf,basis=jkfit}',
    'eehf=energy'],

'b3lyp-d': [
    'gdirect',
    'rks,b3lyp3',
    'eehf=energy',
    'dispcorr',
    'eehfd=eehf+edisp'],

'df-b3lyp-d': [
    'gdirect',
    '{df-rks,b3lyp3,basis=jkfit}',
    'eehf=energy',
    'dispcorr',
    'eehfd=eehf+edisp'],

'b3lyp': [
    'gdirect',
    'rks,b3lyp3',
    'eehf=energy'],

'df-b3lyp': [
    'gdirect',
    '{df-rks,b3lyp3,basis=jkfit}',
    'eehf=energy'],

#'mp2c': [ # this job computes one part [E_disp(TDDFT)] of the three parts of a MP2C calculation
#        # check that nfrag = 2
#         'gdirect',
#         'ga=1101.2; gb=1102.2',
#         'ca=2101.2; cb=2102.2\n',
#
#         $spin = $cgmp{MLPmol1} - 1;
#         'SET,CHARGE=$cgmp{CHGmol1}',
#         'SET,SPIN=$spin',
#         'dummy',
#         foreach $at (@monoBreal) { print $handle ",$at"; }
#         ''
#         '{df-hf,basis=jkfit,locorb=0; start,atdens; save,$ga}',
#         '{df-ks,lhf,df_basis=dflhf,basis_coul=jkfitb,basis_exch=jkfitb; dftfac,1.0; start,$ga; save,$ca}',
#         'eehfa=energy; sapt; monomerA',
#         '',
#
#         $spin = $cgmp{MLPmol2} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol2}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoAreal) { print $handle ",$at"; }
#         print $handle "\n{df-hf,basis=jkfit,locorb=0; start,atdens; save,\$gb}\n";
#         print $handle "{df-ks,lhf,df_basis=dflhf,basis_coul=jkfitb,basis_exch=jkfitb; dftfac,1.0; start,\$gb; save,\$cb}\n";
#         print $handle "eehfb=energy; sapt; monomerB\n\n";
#
#         $spin = $cgmp{MLPsyst} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGsyst}\nSET,SPIN=$spin\n";
#         print $handle "{sapt,SAPT_LEVEL=3; intermol,ca=\$ca,cb=\$cb,icpks=0,fitlevel=3,nlexfac=0.0,cfac=0.0\n";
#         print $handle "dfit,basis_coul=jkfit,basis_exch=jkfit,cfit_scf=3}\n";
#         print $handle "eedisp=E2disp\n\n";
#
#    ],
}

#'dft-sapt-shift': [
#
#         # this is written in an inflexible way (fixed basis, functional) so that it is computed
#         #  only once, then used when writing DFT-SAPT inputs, which we'll be more flexible with
#
#         print $handle "basis={\n";
#         print $handle "set,orbital; default,aug-cc-pVQZ\n";
#         print $handle "set,jkfit;   default,avqz/jkfit\n";
#         print $handle "set,dflhf;   default,avqz/jkfit\n";
#         print $handle "}\n";
#
#         if    ($handle eq "M1OUT") { $charge = $cgmp{CHGmol1}; $spin = $cgmp{MLPmol1} - 1; }
#         elsif ($handle eq "M2OUT") { $charge = $cgmp{CHGmol2}; $spin = $cgmp{MLPmol2} - 1; }
#
#         print $handle "\ngdirect\n";
#         print $handle "{df-ks,pbex,pw91c,lhf; dftfac,0.75,1.0,0.25}\n";
#         print $handle "basis=tzvpp\n";
#         print $handle "{ks,pbe0; orbprint,0}\n";
#         print $handle "eeneut=energy\n";
#         $charge += 1;
#         $spin += 1;
#         print $handle "SET,CHARGE=$charge\nSET,SPIN=$spin\n";
#         print $handle "{ks,pbe0}\n";
#         print $handle "eecat=energy\n";
#         print $handle "eeie=eecat-eeneut\n";
#         print $handle "show[1,20f20.12],ee*,ce*,te*\n";
#         print $handle "show[1,60f20.12],_E*\n";
#    ]
#'dft-sapt': [
#
#         if ( ($asyA eq '') || ($asyB eq '') ) {
#            print "ERROR: asymptotic correction not defined for one or more monomers in index $system.\n";
#            close(DIOUT);
#            unlink("$pathDIOUT");
#         }
#
#         print $handle "gdirect\n";
#         print $handle "ca=2101.2; cb=2102.2\n\n";
#
#         $spin = $cgmp{MLPmol1} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol1}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoBreal) { print $handle ",$at"; }
#         print $handle "\n{df-ks,pbex,pw91c,lhf,df_basis=dflhf,basis_coul=jkfitb,basis_exch=jkfitb; dftfac,0.75,1.0,0.25; asymp,$asyA; save,\$ca}\n";
#         print $handle "eehfa=energy; sapt; monomerA\n\n";
#
#         $spin = $cgmp{MLPmol2} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol2}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoAreal) { print $handle ",$at"; }
#         print $handle "\n{df-ks,pbex,pw91c,lhf,df_basis=dflhf,basis_coul=jkfitb,basis_exch=jkfitb; dftfac,0.75,1.0,0.25; asymp,$asyB; save,\$cb}\n";
#         print $handle "eehfb=energy; sapt; monomerB\n\n";
#
#         $spin = $cgmp{MLPsyst} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGsyst}\nSET,SPIN=$spin\n";
#         print $handle "{sapt,sapt_level=3; intermol,ca=\$ca,cb=\$cb,icpks=0,fitlevel=3,nlexfac=0.0\n";
#         print $handle "dfit,basis_coul=jkfit,basis_exch=jkfit,basis_mp2=mp2fit,cfit_scf=3}\n";
#         print $handle "eeelst=E1pol\n";
#         print $handle "eeexch=E1ex\n";
#         print $handle "eeind=E2ind\n";
#         print $handle "eeexind=E2exind\n";
#         print $handle "eedisp=E2disp\n";
#         print $handle "eeexdisp=E2exdisp\n\n";
#
#    ]
#'dft-sapt-pbe0ac': [
#
#         if ( ($asyA eq '') || ($asyB eq '') ) {
#            print "ERROR: asymptotic correction not defined for one or more monomers in index $system.\n";
#            close(DIOUT);
#            unlink("$pathDIOUT");
#         }
#
#         print $handle "ca=2101.2; cb=2102.2\n\n";
#
#         $spin = $cgmp{MLPmol1} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol1}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoBreal) { print $handle ",$at"; }
#         print $handle "\n{ks,pbe0; asymp,$asyA; save,\$ca}\n";
#         print $handle "eehfa=energy; sapt; monomerA\n\n";
#
#         $spin = $cgmp{MLPmol2} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol2}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoAreal) { print $handle ",$at"; }
#         print $handle "\n{ks,pbe0; asymp,$asyB; save,\$cb}\n";
#         print $handle "eehfb=energy; sapt; monomerB\n\n";
#
#         $spin = $cgmp{MLPsyst} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGsyst}\nSET,SPIN=$spin\n";
#         print $handle "{sapt; intermol,ca=\$ca,cb=\$cb,icpks=0}\n";
#         print $handle "eeelst=E1pol\n";
#         print $handle "eeexch=E1ex\n";
#         print $handle "eeind=E2ind\n";
#         print $handle "eeexind=E2exind\n";
#         print $handle "eedisp=E2disp\n";
#         print $handle "eeexdisp=E2exdisp\n\n";
#    ]
#'dft-sapt-pbe0acalda': [
#
#         if ( ($asyA eq '') || ($asyB eq '') ) {
#            print "ERROR: asymptotic correction not defined for one or more monomers in index $system.\n";
#            close(DIOUT);
#            unlink("$pathDIOUT");
#         }
#
#         print $handle "ca=2101.2; cb=2102.2\n\n";
#
#         $spin = $cgmp{MLPmol1} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol1}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoBreal) { print $handle ",$at"; }
#         print $handle "\n{ks,pbe0; asymp,$asyA; save,\$ca}\n";
#         print $handle "eehfa=energy; sapt; monomerA\n\n";
#
#         $spin = $cgmp{MLPmol2} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGmol2}\nSET,SPIN=$spin\ndummy";
#         foreach $at (@monoAreal) { print $handle ",$at"; }
#         print $handle "\n{ks,pbe0; asymp,$asyB; save,\$cb}\n";
#         print $handle "eehfb=energy; sapt; monomerB\n\n";
#
#         $spin = $cgmp{MLPsyst} - 1;
#         print $handle "SET,CHARGE=$cgmp{CHGsyst}\nSET,SPIN=$spin\n";
#         print $handle "{sapt,sapt_level=3; intermol,ca=\$ca,cb=\$cb,icpks=0,fitlevel=3,nlexfac=0.0\n";
#         print $handle "dfit,basis_coul=jkfit,basis_exch=jkfit,basis_mp2=mp2fit,cfit_scf=3}\n";
#         print $handle "eeelst=E1pol\n";
#         print $handle "eeexch=E1ex\n";
#         print $handle "eeind=E2ind\n";
#         print $handle "eeexind=E2exind\n";
#         print $handle "eedisp=E2disp\n";
#         print $handle "eeexdisp=E2exdisp\n\n";
#
#         print $handle "show[1,20f20.12],ee*,ce*,te*\n";
#         print $handle "show[1,60f20.12],_E*\n";
#      }
#