Geometry convergance problems for organometallic compound

Just Got Here

2:13:00 AM PDT - Tue, Apr 30th 2019

Dear NWChem developers and users:

I have been using NWChem 6.6 to perform computations on organometallic compounds using B3LYP/Lanl2DZ/6-311+G* for a couple of months now but I keep running in geometry convergence issues, i.e. I end up with one or two imaginary frequencies. I know this topic has been discussed numerous times in this forum but all solutions that seemed to work for the other users did not work for me. First of all here is my input file

start ProductCompound+Pd_ActivNewTry

title "ProductCompound+Pd_Activ in 6-311+g* basis set"

scratch_dir ../../../../var/scratch/user

memory total 8000 mb

charge 1

geometry units angstroms print xyz autosym

  C       -0.84169       -4.67838       -0.78686

  C       -1.65397       -3.38660       -0.96403

  H       -1.71396       -3.11420       -2.02404

  H       -2.68067       -3.53152       -0.62818

  C       -0.81108       -5.06104        0.70418

  H       -1.82691       -5.25146        1.06680

  H       -0.24681       -5.98894        0.84871

  C        0.59057       -4.45649       -1.29716

  H        1.18915       -5.36596       -1.17793

  H        0.58450       -4.22342       -2.36829

  C       -1.03586       -2.19919       -0.15858

  C       -0.16400       -3.91760        1.50747

  H       -0.13525       -4.18115        2.57003

  C        1.23527       -3.30783       -0.50083

  H        2.25590       -3.13864       -0.86420

  C       -1.00137       -2.63356        1.33562

  C        0.40263       -1.99950       -0.73850

  C        1.26997       -3.69073        0.99150

  H       -2.02947       -2.81257        1.66299

  H        1.85682       -4.60796        1.10790

  H        1.78896       -2.93546        1.59039

  H        0.26814       -1.91138       -1.82192

  H       -0.62254       -1.83117        1.97726

  H       -1.32011       -5.47621       -1.36413

  C       -2.03374       -1.02677       -0.30303

  N       -1.65887        0.28979       -0.39970

  O       -3.23714       -1.27466       -0.35105

  C       -2.68492        1.24806        0.05181

 Pd        0.18915        1.02023       -0.40189

  C       -0.19606        3.93530        0.45343

  C       -0.93663        5.05841        0.77117

  N       -0.77952        2.74485        0.21089

  H       -3.49750        1.29078       -0.68187

  H       -3.14971        0.89072        0.97802

  C       -2.32714        4.94388        0.83541

  C       -2.11951        2.61368        0.27843

  C       -2.92018        3.71491        0.58740

  H       -0.43566        5.99866        0.96547

  H       -2.93777        5.80587        1.07986

  H       -3.99642        3.59503        0.63422

  H        0.88556        3.96331        0.39216

  C        1.23672       -0.78444       -0.35938

  C        1.99822       -0.17545       -1.42220

  C        1.56525       -0.36838        0.96422

  C        2.90352        0.83525       -1.21601

  C        2.42046        0.72670        1.18037

  C        3.09822        1.33314        0.10341

  H        1.83795       -0.53558       -2.43275

  H        1.14869       -0.86934        1.82439

  H        3.44415        1.25444       -2.05408

  H        2.63506        1.06956        2.18525

  O        3.91464        2.33607        0.41531

  C        4.68995        2.99342       -0.60493

  H        5.37414        2.28919       -1.08103

  H        4.03828        3.45909       -1.34701

 H        5.25736        3.75752       -0.08147

end

basis

  H library 6-311+G*

  O library 6-311+G*

  N library 6-311+G*

  C library 6-311+G*

  Pd library lanl2dz_ecp

end

ecp

  Pd library lanl2dz_ecp

end

stepper

  maxiter 100

end

scf

thresh 1.0 d-8

end

dft

 xc b3lyp

 grid fine

 tolerances tight

 convergence energy 1.0d-8

 iterations 100

end

task dft optimize

task dft freq

This was my last attempt to reach convergence, without success. Before I also tried

1. the "default" driver module with settings tight
2. did computations with and without "autosym"
3. changing the structure randomly to approach the minimum from different directions
4. use the normal vectors of the imaginary frequencies to change the structure along these normal vectors

All attempts without success and I always end up with imaginary frequencies (about -60 cm^-1). Does someone has another idea what to try and could this also come from an incorrect compilation of the NWChem code and corresponding libraries?

Best regards,

StarkEffect

Forum Vet

7:41:33 PM PDT - Tue, Apr 30th 2019
The following reference from Tang, et al. may be helpful but it might be unsuitable in your case. The Electrolyte Genome project: A big data approach in battery materials discovery. Computational Materials Science. 103(2015) 56-67 "The results, shown in Fig. 4, demonstrate that if the frequency is very low (<39.0 cm-1), a high-accuracy frequency calculation can remove up to 42% of the imaginary frequencies. However, high accuracy re-optimization is less helpful for larger imaginary frequencies. In contrast, the two molecular geometry perturbation strategies work extremely well: they remove at least 21 out of the 24 imaginary frequencies for the whole range of imaginary frequencies."

Just Got Here

3:29:19 AM PDT - Thu, Jun 13th 2019
Dear Xiongyan21 and NWChem developers: Thanks for your reply. Unfortunetely, I was unable to resolve the problem posted above. Additionally, I tested other quantum chemistry codes with their implementation of the same level of theory but in these codes the geometry optimization and subsequent frequency computations revealed that the structure has converged and is a local minimum. To double check I asked a friend to run the same computation on another machine with the newest NWChem version. The computation also failed. Therefore, my follow-up question: is anything known about numerical inaccuracies for some special cases? Are there ways to check for probelms or even resolve these? Best regards, StarkEffect

Forum Vet

6:00:24 PM PST - Sat, Nov 23rd 2019

Somebody recommended the decrease of the symmetry, e.g., using C1, to eliminate it. I have just tried thiabendazole, which also can be considered as to be planar even in C1, using GAMESS with 6-311++G(2d,2p) and b3lyp, and the approximate 55 cm-1 imaginary frequency has been removed during HF vibrational anaysis with the warning meaning it is not a true stationary point even if the convergence criteria are met, thus HESSIAN analysis is not rigorious.
The obtained THz results with GAMESS exhibits two peaks in the range of 0-2.0 THz, assuming 1THZ =33.33 cm-1.

MODE FREQ(CM**-1)  SYMMETRY  RED. MASS  IR INTENS.

    1       0.026    A        ...               0.000000

    2       0.004    A        ...               0.000000

    3       0.034    A        ...               0.000000

    4      33.966    A       ...               0.000687

    5      50.002    A       ...               0.010410

...
If the imaginary frequency has not been eliminated, it gives

MODE FREQ(CM**-1) SYMMETRY RED. MASS IR INTENS.

    1      55.684    A"     ...              0.105512

    2       1.891    A"      ...              0.090056

    3       0.114    A"      ...               0.001076

    4       0.052    A'       ...               0.000036

    5       0.034    A'       ...               0.000001

    6       0.814    A'       ...               0.010790

    7       1.182    A"      ...               0.000087

    8      66.596    A"     ...               0.068779

...

                           1           2           3           4           5

      FREQUENCY:        55.68 I      1.89        0.11        0.05        0.03  

       SYMMETRY:         A"        A"         A"         A'          A'  

   REDUCED MASS:     ...

   IR INTENSITY:      0.10551     0.09006     0.00108     0.00004     0.00000

...

NWCHEM SCF and 6-31+G(d) frequency analysis gives
Normal Eigenvalue || Projected Infra Red Intensities

 Mode   [cm**-1]  || [atomic units] [(debye/angs)**2] [(KM/mol)] [arbitrary]

------ ---------- || -------------- ----------------- ---------- -----------

   1       -0.000 ||    0.000331           0.008         0.323       0.196

   2       -0.000 ||    0.002410           0.056         2.350       1.424

   3        0.000 ||    0.000866           0.020         0.844       0.512

   4        0.000 ||    0.000963           0.022         0.939       0.569

   5        0.000 ||    0.000105           0.002         0.102       0.062

   6        0.000 ||    0.000074           0.002         0.073       0.044

   7       92.624 ||    0.000461           0.011         0.449       0.272

   8       99.156 ||    0.000909           0.021         0.886       0.537

  ...

Forum Vet

9:11:04 PM PST - Sat, Nov 30th 2019
With GAMESS, the 6-311++G(2d,2p) optimized Carbendazim can give the following peaks, when doing vibrational analysis using the same basis set and HF, even having the warning. MODE FREQ(CM-1) SYMMETRY RED. MASS IR INTENS. 1 0.013 A ... 0.000000 2 0.008 A ... 0.000000 3 0.006 A ... 0.000000 4 37.138 A ... 0.002125 5 55.879 A ... 0.007225 6 90.644 A ... 0.029414 Please note, after an mini basis set Monte Carlo search** first, the imaginary frequency of this molecule in later vibrational analysis after an geometry optimization has been eliminated. ...

Forum Vet

1:13:34 AM PST - Wed, Dec 25th 2019

The 6-311G++(2d,2p) vibrational analysis of the same functional, i.e., b3lyp, same basis set optimized thiabendazole gives the following using GAMESS.

    1      37.366    A      

    2       7.082    A       

    3       0.027    A       

    4       0.003    A       

    5       0.035    A       

    6       6.201    A       

    7      15.991    A      

    8      65.706    A

Please not there is an imaginary frequency existing, and the frequencies mismatching.

The previously done b3lyp 6-31G+(d,p) optimized Cs thiabendazole gives the following during the RHF Hessian analysis using the same basis set even using Cs symmetry and the same functional
MODE FREQ(CM**-1)

    1       0.046    A'    

    2       0.020    A 

    3       0.008    A'    

    4      36.068    A  

    5      52.861    A'   ,

but 6-31G+(d.p) optimized thiabendazole gives the following during the RHF b3lyp Hessian analysis using the same basis set with C1
MODE FREQ(CM**-1)

    1      33.580    A       

    2       4.732    A        

    3       0.055    A        

    4       0.047    A        

    5       0.022    A        

    6       5.641    A        

    7      14.790    A        

    8      67.181   A

Forum Vet

5:39:18 PM PST - Thu, Dec 26th 2019
The 6-311G++(2d,2p) b3lyp optimized carbendazim whose initial geometry was obtained by an mini basis set Monte Carlo search first, gives no imaginary frequency when doing vibrational analysis using the same functional and same basis set, but the number and positions of peaks changes when compared with the HF result. All the DFT and HF vibrational analyses employs analytical methods.

Forum Vet

9:23:38 PM PST - Sun, Jan 5th 2020

It seems only HF Hessian analysis using GAMESS with the default solver setup can be compared with the experimental data in peak numbers and positions in the THz region.
For example, in the case of alanine, the RHF semi-numerical and analytical methods with 6-311G++(2d,2p) of the geometry obtained by mp2 optimization with the almost same basis set, having a warning meaning the analysis is not very rigorous because of the geometry obtained by optimization is not a true stationary point. The analytical method gives
MODE FREQ(CM**-1) ... IR INTENS.

    1       0.015    A          0.000000

    2       0.008    A          0.000000

    3       0.016    A          0.000000

    4      61.236    A         0.032632

    5      79.238    A         0.013808

which can be compared with the data of D or L alanine in Terahertz absorption spectra of L-,D-,and DL-alanine and their application to determination of enantiomeric composition, in peak numbers and positions, but the mp2 semi-numerical
Hessian analysis with the same basis set and geometry, with no such warning, renders
MODE FREQ(CM**-1) ... IR INTENS.

    1       2.723    A           0.057955

    2       2.307    A           0.101703

    3       1.393    A           0.010289

    4       0.149    A           0.000367

    5       0.340    A           0.001699

    6       2.174    A           0.005833

    7      56.927    A          0.082500

The b3lyp analysis with the default analytical hessian input setup gives
MODE FREQ(CM**-1) ... IR INTENS.

    1      64.741    A          0.012959

    2      38.942    A          0.034262

    3       0.098    A           0.000032

    4       0.007    A           0.000000

    5       0.021    A           0.000000

    6       4.677    A           0.072961

    7      71.778    A          0.191482

Forum Vet

9:28:18 PM PST - Sun, Jan 5th 2020
In the case of Asprin, with rich Thz peaks, it seems also RHF frequency analysis perhaps can be compared with the experimental NIST THZ spectrum, whereas the one using B3lyp mismatches. Toward precise assignment of THz spectra by DFT suggests correction methods to improve the gaps between the theoretically calculatedted and the experimentally obtained. Since the gaps between the values of the theoretically calculated and the experimentally obtained are obvious with the default GAMESS input setup, according to the above article, all those mentioned here are only for educational usage, and I am not interested in any other molecule's THz quantum chemical calculations.

Forum Vet

3:15:21 AM PDT - Fri, Mar 13th 2020
I am now trying the fluorescence and phosphorescence of pesticides, tring to use the same optimization method, if necessary.

Forum >> NWChem's corner >> NWChem functionality