Dear Tian,
I would like to report a reproducible problem that we have encountered in ORCA 6.1.1 and would be very grateful for your advice.
Our study concerns a series of Lewis acid–base complexes formed between Group 13 Lewis acids and the nitrogen atom of a six-membered nitrogen-containing ring (CA). The systems under investigation are
• BF₃–CA
• AlF₃–CA
• GaF₃–CA
• InF₃–CA
• TlF₃–CA
where "CA" denotes the nitrogen-containing six-membered ring used throughout our study.
The final goal of this project is to calculate adiabatic electron affinities for these complexes. Therefore, diffuse basis functions are essential for obtaining reliable energies of both the neutral and anionic species.
To construct the augmented basis, we did not use arbitrary diffuse exponents. Instead, we generated the additional diffuse functions using the even-tempered procedure, and only one additional diffuse S shell and one additional diffuse P shell were added to each element. The same augmentation strategy was applied consistently to all atoms.
The computational protocol was as follows:
1. The lowest-energy conformer was first located using the xTB GOAT conformer search.
2. The lowest-energy conformer obtained from GOAT was then fully optimized using ORCA.
3. Vibrational frequency calculations were subsequently performed.
This protocol works perfectly for BF₃–CA, AlF₃–CA, GaF₃–CA, and InF₃–CA.
However, the TlF₃–CA complex behaves differently.
When the additional diffuse functions are included, geometry optimization converges normally, but the subsequent frequency calculation always produces a single low-frequency imaginary mode (approximately −20 to −30 cm⁻¹).
If the additional diffuse functions are completely removed, the imaginary frequency disappears and the optimized structure becomes a true minimum.
Unfortunately, removing the diffuse functions is not an acceptable solution for our project because the final objective is the calculation of adiabatic electron affinities, where diffuse functions are required.
During the past several days we have carefully investigated this problem and tested many different approaches.
The following attempts were performed:
• Performing a full TightOpt + Freq calculation.
• Performing geometry optimization and frequency calculation as two completely separate jobs.
• Starting from the optimized geometry written in the final XYZ file.
• Using analytical Hessians during optimization.
• Recalculating the Hessian repeatedly during optimization.
• Increasing the maximum optimization iterations.
• Using tighter SCF convergence.
• Enforcing SCF convergence (ConvForced true).
• Using DEFGRID3 integration grids.
• Carefully checking SCF convergence and optimization convergence.
• Repeating the calculations from scratch.
• Verifying that the optimized geometry is reproducible.
In addition, we carefully investigated whether the imaginary frequency originated from an unfavorable initial geometry.
We manually displaced the structure several times along the imaginary mode, generated modified geometries, reoptimized these structures, and repeated the frequency calculations. We also introduced several reasonable structural perturbations and repeated the complete optimization procedure.
Nevertheless, after every optimization the same low-frequency imaginary mode reappeared.
Therefore, the problem does not appear to originate from an unfortunate starting geometry or from being trapped in an obvious local minimum.
At present, this behavior appears to be specific to the TlF₃–CA complex, since all other Lewis acid complexes in the series behave normally under exactly the same computational protocol.
Because the issue only appears when diffuse functions are added to the Tl-containing system, we are wondering whether this could be related to one or more of the following:
• the treatment of very diffuse basis functions on heavy elements,
• the X2C relativistic Hamiltonian,
• Hessian construction,
• numerical integration,
• optimization algorithms,
• or another aspect of ORCA's implementation.
For your convenience, we have attached a ZIP archive containing all relevant files, including:
• input files,
• output files,
• optimized geometries,
• Hessian files,
• GBW files,
• and any other files that may help reproduce and investigate the problem.
If this behavior is expected, we would greatly appreciate an explanation.
If it is not expected, we would be grateful if you could indicate whether this may represent a numerical issue or an implementation issue in ORCA.
Finally, if there is any ORCA-specific solution to this problem—such as using additional keywords, SCF options, geometry optimization settings, Hessian-related options, integration grid settings, or any other input blocks recommended by the ORCA developers—we would greatly appreciate it if you could provide a complete working ORCA 6.1.1 input file (not just isolated keyword suggestions).
Since this computational protocol will be applied consistently to the entire series of Lewis acid complexes (BF₃, AlF₃, GaF₃, InF₃, and TlF₃), we would prefer to use an input that follows the officially recommended ORCA methodology for such systems. If any modification of our current input is advisable, please rewrite the entire input file accordingly.
Thank you very much for your time and your assistance.
We look forward to your advice.

Kind regards,

Saeed

Hello,

I am not sure what you mean by G09 being unable to compute RESP charges. Both the ORCA and Gaussian09 manuals provide the necessary input to compute them.

For Gaussian09, the input is as follows:

```
#P HF/6-31G* Pop=MK IOp(6/33=2) IOp(6/41=15) IOp(6/42=15) IOp(6/50=1)
# Test Units(Ang,Deg)
```

I am also using the MK method with Multiwfn, yet there is still a large discrepancy compared to ORCA, which uses COSMO I believe

Please carefully look at Kollman's original paper of RESP charge (J. Phys. Chem., 97, 10269 (1993)), and relevant part of Multiwfn manual (Section 3.9.16). RESP charge is not a variant of ESP fitting charge just with additional constraints, but it defined specific fitting scheme. The so-called RESP charge in ORCA and Gaussian only enable users to manually set simple constraints, while it doesn't directly follow the full procedure of RESP charge fitting in Kollman's paper, so the results are not comparable.

Even for MK charges, which is one of the simplest ESP fitting charge, the results of Multiwfn and ORCA/Gaussian can also be detectably different even based one exactly the same wavefunction, because the spatial distribution of fitting points are not completely identical in these programs. ESP fitting charge is known to be sensitive to distribution and density of fitting points. In Multiwfn, the default density of fitting points is fine enough and can lead to basically converged result, and the point distribution strictly follows the original MK definition.

Sincerely yours,
Saeed

1 delta-SCF itself is an iterative process, reaction field is intrinsically updated during the SCF iteration, so external iteration is never needed for consideration of solvent effect.

2 ORCA doesn't have this feature, this is one of evident limitations of ORCA (especially for studying CT excitation in solvent environment as linear response model performs poorly in this case). The corrected linear-response model that available in G16 is also not supported by ORCA.

Best regards

R.

In Multiwfn, all indices, including orbital indices, always starts from 1 (PS: I quite dislike starting from 0)

I used the phenanthrene molecule from section 4.14.3 in the Multiwfn manual and compared what I believe is the same calculation in ORCA. For others to see, here is the input file for the Gaussian calculation:

# b3lyp/3-21g pop=nboread

b3lyp/3-21g opted

0 1
C                  0.00000000    3.56061700   -0.29722900
C                  0.00000000    2.83932500    0.87979300
C                  0.00000000    1.42361400    0.86771500
C                  0.00000000    0.72986200   -0.38070000
C                  0.00000000    1.49924400   -1.56931800
C                  0.00000000    2.88151800   -1.53149000
C                  0.00000000    0.67926500    2.09717700
C                  0.00000000   -0.72986200   -0.38070000
C                  0.00000000   -1.42361400    0.86771500
C                  0.00000000   -0.67926500    2.09717700
C                  0.00000000   -2.83932500    0.87979300
H                  0.00000000   -3.34927800    1.83744900
C                  0.00000000   -3.56061700   -0.29722900
C                  0.00000000   -2.88151800   -1.53149000
C                  0.00000000   -1.49924400   -1.56931800
H                  0.00000000    1.23356400    3.02968400
H                  0.00000000    4.64400700   -0.27588200
H                  0.00000000    3.34927800    1.83744900
H                  0.00000000    1.00203800   -2.53051300
H                  0.00000000    3.44643700   -2.45642700
H                  0.00000000   -1.23356400    3.02968400
H                  0.00000000   -4.64400700   -0.27588200
H                  0.00000000   -3.44643700   -2.45642700
H                  0.00000000   -1.00203800   -2.53051300

$nbo dmnao aonao $end

And here is the input for the ORCA calculation:
! B3LYP 3-21G NBO

%nbo
nbokeylist="$nbo aonao dmnao $end"
end

%MaxCore 7500
%pal nprocs 16
           end

* xyz 0 1
C                  0.00000000    3.56061700   -0.29722900
C                  0.00000000    2.83932500    0.87979300
C                  0.00000000    1.42361400    0.86771500
C                  0.00000000    0.72986200   -0.38070000
C                  0.00000000    1.49924400   -1.56931800
C                  0.00000000    2.88151800   -1.53149000
C                  0.00000000    0.67926500    2.09717700
C                  0.00000000   -0.72986200   -0.38070000
C                  0.00000000   -1.42361400    0.86771500
C                  0.00000000   -0.67926500    2.09717700
C                  0.00000000   -2.83932500    0.87979300
H                  0.00000000   -3.34927800    1.83744900
C                  0.00000000   -3.56061700   -0.29722900
C                  0.00000000   -2.88151800   -1.53149000
C                  0.00000000   -1.49924400   -1.56931800
H                  0.00000000    1.23356400    3.02968400
H                  0.00000000    4.64400700   -0.27588200
H                  0.00000000    3.34927800    1.83744900
H                  0.00000000    1.00203800   -2.53051300
H                  0.00000000    3.44643700   -2.45642700
H                  0.00000000   -1.23356400    3.02968400
H                  0.00000000   -4.64400700   -0.27588200
H                  0.00000000   -3.44643700   -2.45642700
H                  0.00000000   -1.00203800   -2.53051300
*

In the Multiwfn manual, the residual density distribution after accepting the first 27 orbitals is:
1 C : 1.0250     2C : 1.0370      3C : 1.0280     4C : 1.0414
5 C : 1.0339      6C : 1.0262      7C : 0.1322      8C : 1.0414  
9 C : 1.0280      10C : 0.1322      11C : 1.0370      12H : 0.0117  
13C : 1.0250      14C : 1.0262      15C : 1.0339      16H : 0.0121  
17H : 0.0113      18H : 0.0117      19H : 0.0126     20H : 0.0111  
21H : 0.0121     22H : 0.0113      23H : 0.0111      24H : 0.0126  

In my ORCA calculations, I receive the following:
1C :  1.0234    2C :  1.0295    3C :  1.0243    4C :  1.0381
5C :  1.0246    6C :  1.0246    7C :  0.1302    8C :  1.0473
9C :  1.0167    10C :  0.1341    11C :  1.0345    12H :  0.0112
13C :  1.0234    14C :  1.0246    15C :  1.0217    16H :  0.0115
17H :  0.0106    18H :  0.0112    19H :  0.0120    20H :  0.0104
21H :  0.0115    22H :  0.0106    23H :  0.0104    24H :  0.0120

Later, the analysis of orbital 31 shows the following in the Multiwfn manual:
NAO# Center    Label    Type            Composition
67      8(C )    px            Val( 2p)     2.247%
76      9(C )    px            Val( 2p)     1.929%
94    11(C )    px            Val( 2p)    13.276%
105    13(C )    px            Val( 2p)    32.846%
114    14(C )    px            Val( 2p)    34.484%
123    15(C )    px            Val( 2p)    15.204%

In my ORCA calculations, I receive the following:
NAO# Center    Label    Type            Composition
67      8(C )    px            Val( 2p)     2.217%
76      9(C )    px            Val( 2p)     1.942%
94    11(C )    px            Val( 2p)    13.387%
105    13(C )    px            Val( 2p)    32.956%
114    14(C )    px            Val( 2p)    34.401%
123    15(C )    px            Val( 2p)    15.082%

As you can see, the differences between the two calculations in not very significant (>0.2 %); however, the visualization of the orbitals is still weird.
These are the orbitals from the Multiwfn manual

And here are the orbitals from my calculation (using a contour value of 0.03)

Saeed

I suggest checking what is the nature of the low-energy transitions (>1000 nm) using hole-electron analysis in Multiwfn, which can provide valuable physical insight.

Also, don't forget to confirm that the reference state wavefunction is stable.