<commit_before>#!/usr/bin/env python
# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
#

'''
Internal-contracted MPS perturbation method.  You can contact <NAME> for the "icpt" program required by this module.  If this method
is used in your work, please cite
<NAME> and <NAME>,  J. Chem. Phys., 136 (2012), 124121
<NAME>, <NAME>, and <NAME>,  J. Chem. Phys., 144 (2016), 034103
'''

import pyscf
import os
import time
import tempfile
from functools import reduce
import numpy
import h5py
from pyscf import lib
from pyscf.lib import logger
from pyscf import mcscf
from pyscf import ao2mo
from pyscf import scf
from pyscf.ao2mo import _ao2mo
from pyscf.dmrgscf import dmrgci
from pyscf.dmrgscf import dmrg_sym
from pyscf import tools
import sys

float_precision = numpy.dtype('Float64')
mpiprefix=""
executable="/home/mussard/softwares/icpt/icpt"
#executable="/home/mussard/softwares/icpt_before_me/icpt"

if not os.path.isfile(executable):
    msg = ('MPSLCC executable %s not found.  Please specify "executable" in %s'
           % (executable, __file__))
    raise ImportError(msg)

#in state average calculationg E1eff will be different than E1
#this means that the h1eff in the fock operator which is stored in eris_sp['h1eff'] will be
#calculated using the E1eff and will in general not result in diagonal matrices
def writeNEVPTIntegrals(mc, E1, E2, E1eff, aaavsplit, nfro, fully_ic=False, third_order=False):
    # Initializations
    ncor = mc.ncore
    nact = mc.ncas
    norb = mc.mo_coeff.shape[1]
    nvir = norb-ncor-nact
    nocc = ncor+nact
    mo   = mc.mo_coeff
    intfolder=mc.fcisolver.scratchDirectory+'/int/'



    # (Note: Integrals are in chemistry notation)
    start = time.time()
    print('Producing the integrals')
    eris = _ERIS(mc, mo)
    eris_sp={}
    if (third_order):
      eris['pcpc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nfro:ncor], mo, mo[:,nfro:ncor]), compact=False)
      eris['ppcc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nfro:ncor], mo[:,nfro:ncor]), compact=False)
      eris['ppee'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nocc:], mo[:,nocc:]), compact=False)
      eris['pcpc'].shape=(norb, ncor-nfro, norb, ncor-nfro)
      eris['ppcc'].shape=(norb, norb, ncor-nfro, ncor-nfro)
      eris['ppee'].shape=(norb, norb, nvir, nvir)

    # h1eff
    eris_sp['h1eff']= eris['h1eff']
    eris_sp['h1eff'][:ncor,:ncor] += numpy.einsum('abcd,cd', eris['ppaa'][:ncor,:ncor,:,:], E1eff)
    eris_sp['h1eff'][:ncor,:ncor] -= numpy.einsum('abcd,bd', eris['papa'][:ncor,:,:ncor,:], E1eff)*0.5
    eris_sp['h1eff'][nocc:,nocc:] += numpy.einsum('abcd,cd', eris['ppaa'][nocc:,nocc:,:,:], E1eff)
    eris_sp['h1eff'][nocc:,nocc:] -= numpy.einsum('abcd,bd', eris['papa'][nocc:,:,nocc:,:], E1eff)*0.5
    numpy.save(intfolder+"int1eff",numpy.asfortranarray(eris_sp['h1eff'][nfro:,nfro:]))

    # CVCV
    eriscvcv = eris['cvcv']
    if (not isinstance(eris_sp['h1eff'], type(eris['cvcv']))):
      eriscvcv = lib.chkfile.load(eris['cvcv'].name, "eri_mo")#h5py.File(eris['cvcv'].name,'r')["eri_mo"]
    eris_sp['cvcv'] = eriscvcv.reshape(ncor, nvir, ncor, nvir)
    end = time.time()
    print('......production of INT took %10.2f sec' %(end-start))
    print('')

    # energy_core
    hcore  = mc.get_hcore()
    dmcore = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    energy_core = numpy.einsum('ij,ji', dmcore, hcore) \
                + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5

    # energyE0
    energyE0 = 1.0*numpy.einsum('ij,ij',     E1, eris_sp['h1eff'][ncor:nocc,ncor:nocc])\
             + 0.5*numpy.einsum('ijkl,ijkl', E2, eris['ppaa'][ncor:nocc,ncor:nocc,:,:].transpose(0,2,1,3))
    energyE0 += energy_core
    energyE0 += mc.mol.energy_nuc()

    print("Energy_nuc  = %13.8f"%(mc.mol.energy_nuc()))
    print("Energy_core = %13.8f"%(energy_core))
    print("Energy      = %13.8f"%(energyE0))
    print("")

    # offdiagonal warning
    offdiagonal = 0.0
    for k in range(ncor):
      for l in range(ncor):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    for k in range(nocc, norb):
      for l in range(nocc,norb):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    if (abs(offdiagonal) > 1e-6):
      print("WARNING: Have to use natural orbitals from CAASCF")
      print("         offdiagonal elements: {:13.6f}".format(offdiagonal))
      print("")

    # Write out ingredients to intfolder
    # 2 "C"
    start = time.time()
    print("Basic ingredients written to "+intfolder,nfro,ncor,nocc,norb)
    numpy.save(intfolder+"W:ccae", numpy.asfortranarray(eris['pacv'][nfro:ncor,     :    , nfro:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eecc", numpy.asfortranarray(eris_sp['cvcv'][nfro: ,     :    , nfro:    ,     :    ].transpose(1,3,0,2)))
    if (third_order):
      numpy.save(intfolder+"W:ceec", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:cccc", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cece", numpy.asfortranarray(eris['ppcc'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (third_order):
      numpy.save(intfolder+"W:ccca", numpy.asfortranarray(eris['ppcc'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:ccce", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cace", numpy.asfortranarray(eris['ppcc'][ncor:nocc, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    # 2 "A"
    numpy.save(intfolder+"W:caac", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,3,1,2)))
    numpy.save(intfolder+"W:ccaa", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aeca", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:eeaa", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aaaa", numpy.asfortranarray(eris['ppaa'][ncor:nocc, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eaca", numpy.asfortranarray(eris['ppaa'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:caca", numpy.asfortranarray(eris['ppaa'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    if (third_order):
      numpy.save(intfolder+"W:aeea", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:aeae", numpy.asfortranarray(eris['ppaa'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (fully_ic):
      numpy.save(intfolder+"W:eaaa", numpy.asfortranarray(eris['ppaa'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:caaa", numpy.asfortranarray(eris['ppaa'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    # 2 "E"
    numpy.save(intfolder+"W:eeca", numpy.asfortranarray(eris['pacv'][nocc:    ,     :    , nfro:    ,     :    ].transpose(3,0,2,1)))
    if (third_order):
      numpy.save(intfolder+"W:aece", numpy.asfortranarray(eris['ppee'][ncor:nocc, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:eeee", numpy.asfortranarray(eris['ppee'][nocc:    ,nocc :    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeec", numpy.asfortranarray(eris['ppee'][nocc:    ,     :    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeea", numpy.asfortranarray(eris['ppee'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))

    # OUTPUT EVERYTHING (for debug of PT3)
    #int2=ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo, mo), compact=False)
    #int2.shape=(norb,norb,norb,norb)
    ##dom=['c','a','e']
    #inout={}
    #inout['c']=[nfro,ncor]
    #inout['a']=[ncor,nocc]
    #inout['e']=[nocc,norb]
    ##for p in range(3):
    ##  for q in range(3):
    ##    for r in range(3):
    ##      for s in range(3):
    ##        name="W:"+dom[p]+dom[q]+dom[r]+dom[s]
    ##        test=int2[inout[dom[p]][0]:inout[dom[p]][1],\
    ##                  inout[dom[r]][0]:inout[dom[r]][1],\
    ##                  inout[dom[q]][0]:inout[dom[q]][1],\
    ##                  inout[dom[s]][0]:inout[dom[s]][1]].transpose(0,2,1,3)
    ##        print("Output: {:} Shape: {:}".format(name,test.shape))
    ##        numpy.save(intfolder+name, numpy.asfortranarray(test))
    #list=['aaaa','aeae','aeca','aece','aeea',\
    #      'caaa','caac','caca','cace','ccaa','ccae','ccca','cccc','ccce','cece','ceec',\
    #      'eaaa','eaca','eeaa','eeca','eecc','eeea','eeec']
    #for i in range(len(list)):
    #  p=list[i][0]
    #  q=list[i][1]
    #  r=list[i][2]
    #  s=list[i][3]
    #  name="W:"+p+q+r+s
    #  test=int2[inout[p][0]:inout[p][1],\
    #            inout[r][0]:inout[r][1],\
    #            inout[q][0]:inout[q][1],\
    #            inout[s][0]:inout[s][1]].transpose(0,2,1,3)
    #  print("Output: {:} Shape: {:}".format(name,test.shape))
    #  numpy.save(intfolder+name, numpy.asfortranarray(test))

    end = time.time()
    print('......savings of INGREDIENTS took %10.2f sec' %(end-start))
    print("")

    # Write "FCIDUMP_aaav0" and "FCIDUMP_aaac"
    if (not fully_ic):
      writeFCIDUMPs_NEVPT(mc,eris,eris_sp,aaavsplit,energy_core,energyE0,nfro)

    return norb, energyE0


def writeMRLCCIntegrals(mc, E1, E2, nfro, fully_ic=False, third_order=False):
    # Initializations
    ncor = mc.ncore
    nact = mc.ncas
    norb = mc.mo_coeff.shape[1]
    nvir = norb-ncor-nact
    nocc = ncor+nact
    mo   = mc.mo_coeff
    intfolder=mc.fcisolver.scratchDirectory+'/int/'



    # (Note: Integrals are in chemistry notation)
    start = time.time()
    print('Producing the integrals')
    eris={}
    eris_sp={}
    eris['pcpc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nfro:ncor], mo, mo[:,nfro:ncor]), compact=False)
    eris['ppcc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nfro:ncor], mo[:,nfro:ncor]), compact=False)
    eris['ppee'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nocc:], mo[:,nocc:]), compact=False)
    eris['pcpc'].shape=(norb, ncor-nfro, norb, ncor-nfro)
    eris['ppcc'].shape=(norb, norb, ncor-nfro, ncor-nfro)
    eris['ppee'].shape=(norb, norb, nvir, nvir)
    eris['papa'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,ncor:nocc], mo, mo[:,ncor:nocc]), compact=False)
    eris['ppaa'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,ncor:nocc], mo[:,ncor:nocc]), compact=False)
    eris['pepe'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nocc:], mo, mo[:,nocc:]), compact=False)
    eris['papa'].shape=(norb, nact, norb, nact)
    eris['ppaa'].shape=(norb, norb, nact, nact)
    eris['pepe'].shape=(norb, nvir, norb, nvir)
    end = time.time()
    print('......production of INT took %10.2f sec' %(end-start))
    print('')

    # h1eff
    hcore    = mc.get_hcore()
    dmcore   = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk   = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore  = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    eris_sp['h1eff'] = reduce(numpy.dot, (mo.T, hcore    , mo))+vhfcore
    numpy.save(intfolder+"int1eff",numpy.asfortranarray(eris_sp['h1eff'][nfro:,nfro:]))

    # int1
    hcore    = mc.get_hcore()
    dmcore   = numpy.dot(mo[:,:nfro], mo[:,:nfro].T)*2
    vj, vk   = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore  = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    int1     = reduce(numpy.dot, (mo.T, hcore    , mo))+vhfcore
    numpy.save(intfolder+"int1",   numpy.asfortranarray(int1[nfro:,nfro:]))

    # energy_core
    hcore  = mc.get_hcore()
    dmcore = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    energy_core = numpy.einsum('ij,ji', dmcore, hcore) \
                + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5

    # energyE0
    energyE0 = 1.0*numpy.einsum('ij,ij',     E1, eris_sp['h1eff'][ncor:nocc,ncor:nocc])\
             + 0.5*numpy.einsum('ijkl,ijkl', E2, eris['ppaa'][ncor:nocc,ncor:nocc,:,:].transpose(0,2,1,3))
    energyE0 += energy_core
    energyE0 += mc.mol.energy_nuc()

    print("Energy_nuc  = %13.8f"%(mc.mol.energy_nuc()))
    print("Energy_core = %13.8f"%(energy_core))
    print("Energy      = %13.8f"%(energyE0))
    print("")

    # offdiagonal warning
    offdiagonal = 0.0
    for k in range(ncor):
      for l in range(ncor):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    for k in range(nocc, norb):
      for l in range(nocc,norb):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    if (abs(offdiagonal) > 1e-6):
      print("WARNING: Have to use natural orbitals from CAASCF")
      print("         offdiagonal elements: %13.8f"%(offdiagonal))
      print("")

    # Write out ingredients to intfolder
    # 2 "C"
    start = time.time()
    print("Basic ingredients written to "+intfolder)
    numpy.save(intfolder+"W:ccae", numpy.asfortranarray(eris['pcpc'][ncor:nocc,     :    , nocc:    ,     :    ].transpose(1,3,0,2)))
    numpy.save(intfolder+"W:eecc", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    if (third_order):
      numpy.save(intfolder+"W:ceec", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:cccc", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cece", numpy.asfortranarray(eris['ppcc'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (third_order):
      numpy.save(intfolder+"W:ccca", numpy.asfortranarray(eris['ppcc'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:ccce", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cace", numpy.asfortranarray(eris['ppcc'][ncor:nocc, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    # 2 "A"
    numpy.save(intfolder+"W:caac", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,3,1,2)))
    numpy.save(intfolder+"W:ccaa", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aeca", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:eeaa", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aaaa", numpy.asfortranarray(eris['ppaa'][ncor:nocc, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eaca", numpy.asfortranarray(eris['ppaa'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:caca", numpy.asfortranarray(eris['ppaa'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    if (third_order):
      numpy.save(intfolder+"W:aeea", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:aeae", numpy.asfortranarray(eris['ppaa'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (fully_ic):
      numpy.save(intfolder+"W:eaaa", numpy.asfortranarray(eris['ppaa'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:caaa", numpy.asfortranarray(eris['ppaa'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    # 2 "E"
    numpy.save(intfolder+"W:eeca", numpy.asfortranarray(eris['pepe'][nfro:ncor,     :    , ncor:nocc,     :    ].transpose(1,3,0,2)))
    if (third_order):
      numpy.save(intfolder+"W:aece", numpy.asfortranarray(eris['ppee'][ncor:nocc, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:eeec", numpy.asfortranarray(eris['ppee'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeea", numpy.asfortranarray(eris['ppee'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))

    feri = h5py.File(intfolder+"int2eeee.hdf5", 'w')
    ao2mo.full(mc.mol, mo[:,nocc:], feri, compact=False)
    for o in range(nvir):
      int2eee = feri['eri_mo'][o*(norb-nocc):(o+1)*(norb-nocc),:]
      numpy.asfortranarray(int2eee).tofile(intfolder+"W:eeee%04d"%(o))
<commit_msg>Fix integrals writing for MRLCC2<commit_after>#!/usr/bin/env python
# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
#

'''
Internal-contracted MPS perturbation method.  You can contact Sandeep
Sharma for the "icpt" program required by this module.  If this method
is used in your work, please cite
<NAME> and <NAME>,  J. Chem. Phys., 136 (2012), 124121
<NAME>, <NAME>, and <NAME>,  J. Chem. Phys., 144 (2016), 034103
'''

import pyscf
import os
import time
import tempfile
from functools import reduce
import numpy
import h5py
from pyscf import lib
from pyscf.lib import logger
from pyscf import mcscf
from pyscf import ao2mo
from pyscf import scf
from pyscf.ao2mo import _ao2mo
from pyscf.dmrgscf import dmrgci
from pyscf.dmrgscf import dmrg_sym
from pyscf import tools
import sys

float_precision = numpy.dtype('Float64')
mpiprefix=""
executable="/home/mussard/softwares/icpt/icpt"
#executable="/home/mussard/softwares/icpt_before_me/icpt"

if not os.path.isfile(executable):
    msg = ('MPSLCC executable %s not found.  Please specify "executable" in %s'
           % (executable, __file__))
    raise ImportError(msg)

#in state average calculationg E1eff will be different than E1
#this means that the h1eff in the fock operator which is stored in eris_sp['h1eff'] will be
#calculated using the E1eff and will in general not result in diagonal matrices
def writeNEVPTIntegrals(mc, E1, E2, E1eff, aaavsplit, nfro, fully_ic=False, third_order=False):
    # Initializations
    ncor = mc.ncore
    nact = mc.ncas
    norb = mc.mo_coeff.shape[1]
    nvir = norb-ncor-nact
    nocc = ncor+nact
    mo   = mc.mo_coeff
    intfolder=mc.fcisolver.scratchDirectory+'/int/'



    # (Note: Integrals are in chemistry notation)
    start = time.time()
    print('Producing the integrals')
    eris = _ERIS(mc, mo)
    eris_sp={}
    if (third_order):
      eris['pcpc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nfro:ncor], mo, mo[:,nfro:ncor]), compact=False)
      eris['ppcc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nfro:ncor], mo[:,nfro:ncor]), compact=False)
      eris['ppee'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nocc:], mo[:,nocc:]), compact=False)
      eris['pcpc'].shape=(norb, ncor-nfro, norb, ncor-nfro)
      eris['ppcc'].shape=(norb, norb, ncor-nfro, ncor-nfro)
      eris['ppee'].shape=(norb, norb, nvir, nvir)

    # h1eff
    eris_sp['h1eff']= eris['h1eff']
    eris_sp['h1eff'][:ncor,:ncor] += numpy.einsum('abcd,cd', eris['ppaa'][:ncor,:ncor,:,:], E1eff)
    eris_sp['h1eff'][:ncor,:ncor] -= numpy.einsum('abcd,bd', eris['papa'][:ncor,:,:ncor,:], E1eff)*0.5
    eris_sp['h1eff'][nocc:,nocc:] += numpy.einsum('abcd,cd', eris['ppaa'][nocc:,nocc:,:,:], E1eff)
    eris_sp['h1eff'][nocc:,nocc:] -= numpy.einsum('abcd,bd', eris['papa'][nocc:,:,nocc:,:], E1eff)*0.5
    numpy.save(intfolder+"int1eff",numpy.asfortranarray(eris_sp['h1eff'][nfro:,nfro:]))

    # CVCV
    eriscvcv = eris['cvcv']
    if (not isinstance(eris_sp['h1eff'], type(eris['cvcv']))):
      eriscvcv = lib.chkfile.load(eris['cvcv'].name, "eri_mo")#h5py.File(eris['cvcv'].name,'r')["eri_mo"]
    eris_sp['cvcv'] = eriscvcv.reshape(ncor, nvir, ncor, nvir)
    end = time.time()
    print('......production of INT took %10.2f sec' %(end-start))
    print('')

    # energy_core
    hcore  = mc.get_hcore()
    dmcore = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    energy_core = numpy.einsum('ij,ji', dmcore, hcore) \
                + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5

    # energyE0
    energyE0 = 1.0*numpy.einsum('ij,ij',     E1, eris_sp['h1eff'][ncor:nocc,ncor:nocc])\
             + 0.5*numpy.einsum('ijkl,ijkl', E2, eris['ppaa'][ncor:nocc,ncor:nocc,:,:].transpose(0,2,1,3))
    energyE0 += energy_core
    energyE0 += mc.mol.energy_nuc()

    print("Energy_nuc  = %13.8f"%(mc.mol.energy_nuc()))
    print("Energy_core = %13.8f"%(energy_core))
    print("Energy      = %13.8f"%(energyE0))
    print("")

    # offdiagonal warning
    offdiagonal = 0.0
    for k in range(ncor):
      for l in range(ncor):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    for k in range(nocc, norb):
      for l in range(nocc,norb):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    if (abs(offdiagonal) > 1e-6):
      print("WARNING: Have to use natural orbitals from CAASCF")
      print("         offdiagonal elements: {:13.6f}".format(offdiagonal))
      print("")

    # Write out ingredients to intfolder
    # 2 "C"
    start = time.time()
    print("Basic ingredients written to "+intfolder,nfro,ncor,nocc,norb)
    numpy.save(intfolder+"W:ccae", numpy.asfortranarray(eris['pacv'][nfro:ncor,     :    , nfro:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eecc", numpy.asfortranarray(eris_sp['cvcv'][nfro: ,     :    , nfro:    ,     :    ].transpose(1,3,0,2)))
    if (third_order):
      numpy.save(intfolder+"W:ceec", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:cccc", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cece", numpy.asfortranarray(eris['ppcc'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (third_order):
      numpy.save(intfolder+"W:ccca", numpy.asfortranarray(eris['ppcc'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:ccce", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cace", numpy.asfortranarray(eris['ppcc'][ncor:nocc, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    # 2 "A"
    numpy.save(intfolder+"W:caac", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,3,1,2)))
    numpy.save(intfolder+"W:ccaa", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aeca", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:eeaa", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aaaa", numpy.asfortranarray(eris['ppaa'][ncor:nocc, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eaca", numpy.asfortranarray(eris['ppaa'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:caca", numpy.asfortranarray(eris['ppaa'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    if (third_order):
      numpy.save(intfolder+"W:aeea", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
      numpy.save(intfolder+"W:aeae", numpy.asfortranarray(eris['ppaa'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (fully_ic):
      numpy.save(intfolder+"W:eaaa", numpy.asfortranarray(eris['ppaa'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:caaa", numpy.asfortranarray(eris['ppaa'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    # 2 "E"
    numpy.save(intfolder+"W:eeca", numpy.asfortranarray(eris['pacv'][nocc:    ,     :    , nfro:    ,     :    ].transpose(3,0,2,1)))
    if (third_order):
      numpy.save(intfolder+"W:aece", numpy.asfortranarray(eris['ppee'][ncor:nocc, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:eeee", numpy.asfortranarray(eris['ppee'][nocc:    ,nocc :    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeec", numpy.asfortranarray(eris['ppee'][nocc:    ,     :    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeea", numpy.asfortranarray(eris['ppee'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))

    # OUTPUT EVERYTHING (for debug of PT3)
    #int2=ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo, mo), compact=False)
    #int2.shape=(norb,norb,norb,norb)
    ##dom=['c','a','e']
    #inout={}
    #inout['c']=[nfro,ncor]
    #inout['a']=[ncor,nocc]
    #inout['e']=[nocc,norb]
    ##for p in range(3):
    ##  for q in range(3):
    ##    for r in range(3):
    ##      for s in range(3):
    ##        name="W:"+dom[p]+dom[q]+dom[r]+dom[s]
    ##        test=int2[inout[dom[p]][0]:inout[dom[p]][1],\
    ##                  inout[dom[r]][0]:inout[dom[r]][1],\
    ##                  inout[dom[q]][0]:inout[dom[q]][1],\
    ##                  inout[dom[s]][0]:inout[dom[s]][1]].transpose(0,2,1,3)
    ##        print("Output: {:} Shape: {:}".format(name,test.shape))
    ##        numpy.save(intfolder+name, numpy.asfortranarray(test))
    #list=['aaaa','aeae','aeca','aece','aeea',\
    #      'caaa','caac','caca','cace','ccaa','ccae','ccca','cccc','ccce','cece','ceec',\
    #      'eaaa','eaca','eeaa','eeca','eecc','eeea','eeec']
    #for i in range(len(list)):
    #  p=list[i][0]
    #  q=list[i][1]
    #  r=list[i][2]
    #  s=list[i][3]
    #  name="W:"+p+q+r+s
    #  test=int2[inout[p][0]:inout[p][1],\
    #            inout[r][0]:inout[r][1],\
    #            inout[q][0]:inout[q][1],\
    #            inout[s][0]:inout[s][1]].transpose(0,2,1,3)
    #  print("Output: {:} Shape: {:}".format(name,test.shape))
    #  numpy.save(intfolder+name, numpy.asfortranarray(test))

    end = time.time()
    print('......savings of INGREDIENTS took %10.2f sec' %(end-start))
    print("")

    # Write "FCIDUMP_aaav0" and "FCIDUMP_aaac"
    if (not fully_ic):
      writeFCIDUMPs_NEVPT(mc,eris,eris_sp,aaavsplit,energy_core,energyE0,nfro)

    return norb, energyE0


def writeMRLCCIntegrals(mc, E1, E2, nfro, fully_ic=False, third_order=False):
    # Initializations
    ncor = mc.ncore
    nact = mc.ncas
    norb = mc.mo_coeff.shape[1]
    nvir = norb-ncor-nact
    nocc = ncor+nact
    mo   = mc.mo_coeff
    intfolder=mc.fcisolver.scratchDirectory+'/int/'



    # (Note: Integrals are in chemistry notation)
    start = time.time()
    print('Producing the integrals')
    eris={}
    eris_sp={}
    eris['pcpc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nfro:ncor], mo, mo[:,nfro:ncor]), compact=False)
    eris['ppcc'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nfro:ncor], mo[:,nfro:ncor]), compact=False)
    eris['ppee'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,nocc:], mo[:,nocc:]), compact=False)
    eris['pcpc'].shape=(norb, ncor-nfro, norb, ncor-nfro)
    eris['ppcc'].shape=(norb, norb, ncor-nfro, ncor-nfro)
    eris['ppee'].shape=(norb, norb, nvir, nvir)
    eris['papa'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,ncor:nocc], mo, mo[:,ncor:nocc]), compact=False)
    eris['ppaa'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo, mo[:,ncor:nocc], mo[:,ncor:nocc]), compact=False)
    eris['pepe'] = ao2mo.outcore.general_iofree(mc.mol, (mo, mo[:,nocc:], mo, mo[:,nocc:]), compact=False)
    eris['papa'].shape=(norb, nact, norb, nact)
    eris['ppaa'].shape=(norb, norb, nact, nact)
    eris['pepe'].shape=(norb, nvir, norb, nvir)
    end = time.time()
    print('......production of INT took %10.2f sec' %(end-start))
    print('')

    # h1eff
    hcore    = mc.get_hcore()
    dmcore   = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk   = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore  = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    eris_sp['h1eff'] = reduce(numpy.dot, (mo.T, hcore    , mo))+vhfcore
    numpy.save(intfolder+"int1eff",numpy.asfortranarray(eris_sp['h1eff'][nfro:,nfro:]))

    # int1
    hcore    = mc.get_hcore()
    dmcore   = numpy.dot(mo[:,:nfro], mo[:,:nfro].T)*2
    vj, vk   = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore  = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    int1     = reduce(numpy.dot, (mo.T, hcore    , mo))+vhfcore
    numpy.save(intfolder+"int1",   numpy.asfortranarray(int1[nfro:,nfro:]))

    # energy_core
    hcore  = mc.get_hcore()
    dmcore = numpy.dot(mo[:,:ncor], mo[:,:ncor].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    energy_core = numpy.einsum('ij,ji', dmcore, hcore) \
                + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5

    # energyE0
    energyE0 = 1.0*numpy.einsum('ij,ij',     E1, eris_sp['h1eff'][ncor:nocc,ncor:nocc])\
             + 0.5*numpy.einsum('ijkl,ijkl', E2, eris['ppaa'][ncor:nocc,ncor:nocc,:,:].transpose(0,2,1,3))
    energyE0 += energy_core
    energyE0 += mc.mol.energy_nuc()

    print("Energy_nuc  = %13.8f"%(mc.mol.energy_nuc()))
    print("Energy_core = %13.8f"%(energy_core))
    print("Energy      = %13.8f"%(energyE0))
    print("")

    # offdiagonal warning
    offdiagonal = 0.0
    for k in range(ncor):
      for l in range(ncor):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    for k in range(nocc, norb):
      for l in range(nocc,norb):
        if(k != l):
          offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    if (abs(offdiagonal) > 1e-6):
      print("WARNING: Have to use natural orbitals from CAASCF")
      print("         offdiagonal elements: %13.8f"%(offdiagonal))
      print("")

    # Write out ingredients to intfolder
    # 2 "C"
    start = time.time()
    print("Basic ingredients written to "+intfolder)
    numpy.save(intfolder+"W:ccae", numpy.asfortranarray(eris['pcpc'][ncor:nocc,     :    , nocc:    ,     :    ].transpose(1,3,0,2)))
    numpy.save(intfolder+"W:eecc", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    #
    numpy.save(intfolder+"W:ceec", numpy.asfortranarray(eris['pcpc'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:cccc", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
    numpy.save(intfolder+"W:cece", numpy.asfortranarray(eris['ppcc'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (third_order):
      numpy.save(intfolder+"W:ccca", numpy.asfortranarray(eris['ppcc'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:ccce", numpy.asfortranarray(eris['ppcc'][nfro:ncor, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:cace", numpy.asfortranarray(eris['ppcc'][ncor:nocc, nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    # 2 "A"
    numpy.save(intfolder+"W:caac", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,3,1,2)))
    numpy.save(intfolder+"W:ccaa", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nfro:ncor,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aeca", numpy.asfortranarray(eris['papa'][nfro:ncor,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:eeaa", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:aaaa", numpy.asfortranarray(eris['ppaa'][ncor:nocc, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:eaca", numpy.asfortranarray(eris['ppaa'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    numpy.save(intfolder+"W:caca", numpy.asfortranarray(eris['ppaa'][nfro:ncor, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
    #
    numpy.save(intfolder+"W:aeea", numpy.asfortranarray(eris['papa'][nocc:    ,     :    , nocc:    ,     :    ].transpose(1,2,0,3)))
    numpy.save(intfolder+"W:aeae", numpy.asfortranarray(eris['ppaa'][nocc:    , nocc:    ,     :    ,     :    ].transpose(2,0,3,1)))
    if (fully_ic):
      numpy.save(intfolder+"W:eaaa", numpy.asfortranarray(eris['ppaa'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:caaa", numpy.asfortranarray(eris['ppaa'][nfro:ncor, ncor:nocc,     :    ,     :    ].transpose(0,2,1,3)))
    # 2 "E"
    numpy.save(intfolder+"W:eeca", numpy.asfortranarray(eris['pepe'][nfro:ncor,     :    , ncor:nocc,     :    ].transpose(1,3,0,2)))
    if (third_order):
      numpy.save(intfolder+"W:aece", numpy.asfortranarray(eris['ppee'][ncor:nocc, nfro:ncor,     :    ,     :    ].transpose(0,2,1,3)))
      numpy.save(intfolder+"W:eeec", numpy.asfortranarray(eris['ppee'][nocc:    , nfro:ncor,     :    ,     :    ].transpose(2,0,3,1)))
      numpy.save(intfolder+"W:eeea", numpy.asfortranarray(eris['ppee'][nocc:    , ncor:nocc,     :    ,     :    ].transpose(2,0,3,1)))

    feri = h5py.File(intfolder+"int2eeee.hdf5", 'w')
    ao2mo.full(mc.mol, mo[:,nocc:], feri, compact=False)
    for o in range(nvir):
      int2eee = feri['eri_mo'][o*(norb-nocc):(o+1)*(norb-nocc),:]
      numpy.asfortranarray(int2eee).tofile(intfolder+"W:eeee%04d"%(o))

{'fraction_non_alphanumeric': 0.20062195087158713, 'fraction_numerical': 0.024393027303072665, 'mean_word_length': 3.3316856154226397, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4367479', 'n_tokens_mistral': 14855, 'n_tokens_neox': 14730, 'n_words': 2073}