[Debichem-devel] Bug#695540: nwchem: Grave numerical errors in MCSCF analytical gradients and DFT semidirect calculations

[Michael Banck]

Package: nwchem
Version: 6.1-3
Severity: serious
Justification: Renders package unfit for release

Version 6.1 has two grave regressions compared to 6.0:

1. Analytic MCSCF gradients are totally wrong, resulting in non
converging geometry optimizations

2. Semidirect DFT calculations result in qualitatively wrong energies

Ad 1:

This is the input file (QA/tests/mcscf_ch2/mcscf_ch2.nw in the source
tree):

--8<--
echo
start ch2

title "CASSCF for CH2 1A1 - 3B1 splitting at the 3B1 geometry"

geometry units au
  C    0   0     0
  H    0  1.87  -0.82
  symmetry c2v
end

basis
  C library 6-31g**
  H library 6-31g**
end

task scf

mcscf
  active  6
  actelec 6
  state   3B1
end

task mcscf optimize

mcscf
  state  1A1
end

task mcscf
--8<--

This is the output from 6.0-2, piped through "grep ^@ | cut -c -68":

@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax  
@ ---- ---------------- -------- -------- -------- -------- --------
@    0     -38.95855140  0.0D+00  0.00740  0.00610  0.00000  0.00000
@    1     -38.95874023 -1.9D-04  0.00103  0.00059  0.00865  0.01516
@    2     -38.95876004 -2.0D-05  0.00046  0.00037  0.00637  0.01345
@    3     -38.95876065 -6.1D-07  0.00000  0.00000  0.00063  0.00132
@    3     -38.95876065 -6.1D-07  0.00000  0.00000  0.00063  0.00132

This is the output from 6.1-3:

@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax  
@ ---- ---------------- -------- -------- -------- -------- --------
@    0     -38.95855140  0.0D+00  0.07784  0.07395  0.00000  0.00000
@    1     -38.95696026  1.6D-03  0.05334  0.04885  0.04900  0.10113
@    2     -38.95570226  1.3D-03  0.05502  0.04832  0.01535  0.02930
@    3     -38.95665326 -9.5D-04  0.06689  0.05471  0.01439  0.02852
@    4     -38.95717107 -5.2D-04  0.07566  0.06209  0.01026  0.01901
@    5     -38.95708206  8.9D-05  0.05443  0.05175  0.00229  0.00432
@    6     -38.95701518  6.7D-05  0.07294  0.05985  0.00118  0.00248
@    7     -38.95693884  7.6D-05  0.05269  0.05043  0.00121  0.00256
@    8     -38.95691119  2.8D-05  0.06577  0.05379  0.00042  0.00086
@    9     -38.95691174 -5.5D-07  0.18549  0.15361  0.00002  0.00004
@   10     -38.95697624 -6.4D-05  0.04261  0.03892  0.00115  0.00234
@   11     -38.95700216 -2.6D-05  0.04244  0.03878  0.00044  0.00075
@   12     -38.95842296 -1.4D-03  0.09093  0.07566  0.03564  0.06462
@   13     -38.95819225  2.3D-04  0.09130  0.07554  0.00882  0.01549
@   14     -38.95767335  5.2D-04  0.05886  0.05604  0.01328  0.02317
@   15     -38.95757527  9.8D-05  0.05789  0.05516  0.00206  0.00364
@   16     -38.95736754  2.1D-04  0.05598  0.05343  0.00407  0.00715
@   17     -38.95714460  2.2D-04  0.07376  0.06058  0.00401  0.00704
@   18     -38.95711506  3.0D-05  0.18058  0.14965  0.00049  0.00094
@   19     -38.95711712 -2.1D-06  0.05094  0.04951  0.00004  0.00006
@   20     -38.95711525  1.9D-06  0.04207  0.03840  0.00010  0.00025
@   20     -38.95711525  1.9D-06  0.04207  0.03840  0.00010  0.00025

Note the random values for Gmax (maximum gradient) and Grms (median
gradient).  The test suite calculation aborted after the default maximum
amount of 20 iterations without convergence.

Ad 2:

This is the input file (adopted from
QA/tests/dft_semidirect/dft_semidirect.nw in the source tree, which
shows the same behaviour, but is more expensive to run):

--8<--
echo

start dft_semidirect

memory stack 300 mb heap 150 mb global 250 mb

# do calculation with direct for comparison

set int:txs:limxmem 30000000

geometry units angstrom noautoz nocenter noautosym
O         0.00000000000000      0.00000000000000     -0.26192472620000
H         0.00000000000000      0.77162768440000      0.34261631290000
H         0.00000000000000     -0.77162768440000      0.34261631290000
end

basis "ao basis" spherical
  * library cc-pvdz
end

tddft
  nroots 2
end

dft
  xc b3lyp
  convergence energy 1e-8 density 1e-7 gradient 1e-7
  vectors input atomic
  direct
end

task tddft energy

# test new method of semidirect input

dft
  xc b3lyp
  convergence energy 1e-8 density 1e-7 gradient 1e-7
  vectors input atomic
  semidirect memsize 51200000 filesize 0 # use 400 mb of memory and no disk
end

task tddft energy

# test old method of semidirect input

scf
  semidirect memsize 51200000 filesize 0 # use 400 mb of memory and no disk
end

dft
  xc b3lyp
  convergence energy 1e-8 density 1e-7 gradient 1e-7
  vectors input atomic
end

task tddft energy
--8<--

This is the output from 6.0-2, piped through "egrep -e
'(Total.DFT|Excitation)'":

         Total DFT energy =      -76.420328702970
        Excitation energy =      0.272382539701
        Excitation energy =      0.245024871199
         Total DFT energy =      -76.420328702972
        Excitation energy =      0.272382539704
        Excitation energy =      0.245024871205
         Total DFT energy =      -76.420328702972
        Excitation energy =      0.272382539704
        Excitation energy =      0.245024871205

This is the output from 6.1-3:

         Total DFT energy =      -76.420328702972
        Excitation energy =      0.272382539701
        Excitation energy =      0.245024871199
         Total DFT energy =      -84.163008005292
        Excitation energy =      0.432147825828
        Excitation energy =      0.392466494410
         Total DFT energy =      -92.285951773846
        Excitation energy =      0.586049142135
        Excitation energy =      0.532145253480

Note how the total energies are different for 6.1-3 between the three
(physically identical) runs, same for the excitation energies.


Cheers,

Michael