How to generate wfn file with ECP

I want to get wfn file (for QTAIM analysis) from DFT calculations using ECP. I use
The resulting file looks as follows:

GAUSSIAN 113 MOL ORBITALS 1922 PRIMITIVES 58 NUCLEI
END DATA
THE HF ENERGY = -2332.856626577337 THE VIRIAL(-V/T)= 0.00000000

Without ECP, the correct wfn file is generated.

Can you please post both a working (no ECP) and a failing (with ECP) input file? The following input file seems to work fine with the latest Q-Chem trunk

$molecule
-1 1
        I         0.6539056222    0.0300979939    0.0000000000
        O        -2.8593081518    0.0835488510    0.0000000000
        H        -1.8947008157    0.2596704547    0.0000000000
        H        -2.8641221666   -0.8778240298    0.0000000000
$end

$rem
method               hf
basis                gen
purecart             111
print_orbitals       true
print_general_basis  true
write_wfn            aimecp
ecp                  def2-ecp
def
$end

$basis
i   0
def2-tzvp
****
O   0
6-31G
****
H   0
6-31G
****
$end

The job runs to completion and a .wfn file is generated, which seems to be complete although I don’t have a way to test it. ECP is definitely turned on because there are only 36 electrons, not 64.

Thank you, John.
This simple input works for me:

$molecule
0 1
C 0. 0. 0.
O 0. 0. 1.1314
$end

$rem
METHOD PBE0
BASIS def2-TZVPP
MAX_SCF_CYCLES 200
WRITE_WFN co
IQMOL_FCHK = TRUE
MOLDEN_FORMAT TRUE
XC_GRID 3
$end

It produces co.fchk and CO.wfn, but no molden file. Then I can use CO.wfn with AIMALL program.

The second file looks as follows:
$molecule
1 1
C 1.219637 0.825015 0.614786
C 1.241431 2.195814 0.813241
C 0.080034 2.916204 1.051608
C -1.130170 2.266752 0.873943
C -1.197290 0.896002 0.664667
C -0.008819 0.087689 0.763446
H 2.184662 2.724232 0.774701
H 0.117517 3.975693 1.261592
H -2.044001 2.848463 0.885725
P -2.733264 0.199523 0.063802
P 2.694047 -0.010194 0.033185
S 2.149921 -1.777090 -0.794615
S -2.483534 -1.767023 -0.299905
C -4.058966 0.494514 1.278706
H -4.119064 1.584778 1.371020
C 3.853276 -0.340177 1.411155
H 3.328942 -1.144846 1.939654
C 3.503182 0.963666 -1.301279
H 4.187447 0.205044 -1.701282
C -3.197806 1.022156 -1.514096
H -4.046525 0.406086 -1.833435
Rh -0.160358 -1.637214 -0.531989
C -0.010359 -1.227648 1.522862
H -0.046319 -2.340641 0.971571
H 0.913360 -1.365210 2.077926
H -0.885368 -1.324227 2.155985
C 5.168087 -0.905794 0.887429
H 5.767845 -1.260781 1.726705
H 5.007230 -1.749411 0.214391
H 5.753676 -0.149618 0.362998
C 4.033774 0.825850 2.373798
H 3.081083 1.181452 2.766335
H 4.639464 0.494486 3.219068
H 4.551473 1.667797 1.915004
C 4.324739 2.177929 -0.887643
H 5.073533 1.949845 -0.131827
H 4.853417 2.548571 -1.767906
H 3.702354 2.996418 -0.527500
C 2.489228 1.310224 -2.384128
H 1.926277 0.434162 -2.705701
H 1.787267 2.069347 -2.035285
H 3.015678 1.711398 -3.251433
C -3.662802 -0.094112 2.625635
H -3.568297 -1.179526 2.561762
H -2.721070 0.320436 2.988981
H -4.434127 0.132355 3.363151
C -5.390833 -0.049331 0.777365
H -5.333267 -1.126511 0.609428
H -6.161073 0.128982 1.529082
H -5.716072 0.425930 -0.148574
C -2.085798 0.882318 -2.542633
H -1.200434 1.446885 -2.243185
H -1.804276 -0.161266 -2.692510
H -2.427383 1.273475 -3.502265
C -3.656405 2.467674 -1.373785
H -2.823970 3.129769 -1.134699
H -4.065708 2.800610 -2.329197
H -4.436942 2.597853 -0.623712
$end

$rem
METHOD PBE0
BASIS def2-TZVP
ECP def2-ecp
SCF_PRINT TRUE
WRITE_WFN ag5-tmp
PRINT_ORBITALS TRUE
print_general_basis true
$end

If I use def2-TZVPP basis set, it produced the 3-line wfn file that I send before. In this input file, I use def2-TZVP and get wfn file that looks correct, but AIMALL says: “Invalid traditional wfn file! Could be improperly formatted 2nd line or other line.” So, my first problem was in G-functions on the transition metal and not ECP. Can you say what is the problem now?

I ran those two inputs using the current Q-Chem trunk (pre-release v. 6.1). For the first (little) one, it does indeed produce a MolDen input file:

======= MOLDEN-FORMATTED INPUT FILE FOLLOWS =======
[Molden Format]
[Atoms] (Angs)
  C      1    6      0.00000000      0.00000000     -0.64651429
  O      2    8      0.00000000      0.00000000      0.48488571

======= END OF MOLDEN-FORMATTED INPUT FILE =======

This not a separate file but is contained at the end of the Q-Chem output file. If you want it to contain orbitals you need to stipulate PRINT_ORBITALS = TRUE in $rem. Using IQMOL_FCHK=TRUE to obtain .fchk file, and loading that into IQmol, is another way to visualize orbitals.

For the second job, Q-Chem does generate a file AG5-TMP.wfn. I don’t run QTAIM calculations so I can’t say for sure that it’s a working .wfn file, but it does look complete. Starts with


GAUSSIAN            113 MOL ORBITALS    1591 PRIMITIVES       58 NUCLEI
  C    1    (CENTRE  1)   2.30477994  1.55905244  1.16177704  CHARGE =  6.0

and ends with

END DATA
THE  HF ENERGY =   -2332.902175761614 THE VIRIAL(-V/T)=    0.00000000

(If you send me an email I will send you the file so that you can verify whether it loads in AIMALL.) What version of Q-Chem are you using?

It appears that I had two problems:

  1. If basis set contains g-function, def2-TZVPP in my case, wfn file is not generated at all. Using def2-TZVP, I get wfn file that cannot be used for topological analysis by AIMALL or Multiwfn.
  2. AIMALL can open fchk file from QChem calculations. But wavefunction in this file is wrong if ECP is used, and AIMALL finds non-nuclear attractor critical points on the atom with ECP.

(1) Are you certain that these programs (AIMALL and Multiwfn) support g functions?

(2) A wild guess might be that somehow the number of electrons or the atomic numbers are not adjusted properly (both should be reduced to account for the electrons that are replaced by the ECP). If I look at the .fchk file for the iodide-water input above (input contains a small typo at the end of $rem, apologies), then I find

Charge                                     I               -1
Multiplicity                               I                1
Number of electrons                        I               36
Number of alpha electrons                  I               18
Number of beta electrons                   I               18
Atomic numbers                             I   N=           4
          53           8           1           1

Here, the number of electrons is reduced by iodine still has its full atomic number. It’s not clear to me that this is a bug (IQmol may read it properly), but perhaps an incompatibility. I wonder what would happen if you were to manually change iodine’s atomic number to 25? (That is, reduce it by the 28, corresponding to the number of electrons that are replaced by the ECP). Does AIMALL produce sensible results in that case?

  1. Yes, I used def2-TZVPP basis set with g-functions on Rh in Gaussian, and the resulting wfn/wfx files were smoothly opened by AIMALL.
  2. It does not help. I tried it and get the same error message. The ECP story is not so simple for topological analysis because the spatial distribution of the electron density near the nucleus is incorrect.