Strange results with m06-2x...
Clicked A Few Times
11:17:07 AM PST - Wed, Nov 14th 2012
Hello. I'm trying to get NWChem 6.0 to reproduce some G09 results so
that I can use some features of NWChem that don't exist on G09 (namely
nudged elastic band) and am
running into some issues.

When I include the line
grid lebedev 90 14 ssf euler
in the input file, I can can get good agreement with G09 (0.004 eV) on
the total energy using PBE. However, when I use m06-2X, I get very
different results for the total energy (54 Hartrees lower in NWChem)
and the geometry that it's relaxing to is very odd (the bond lengths
of all atoms are much larger). It seems that I am calling for m06-2x
properly ( xc xm06-2x cm06-2x), as it's properly specified in the
output file

             XC Information

             --------------

                       M06-2X Exchange Functional  1.000

                     M06-2X Correlation Potential  1.000 non-local

But it also includes the line:
cm06-2x uses PW91LDA as defaults.These defaults cannot be changed
above. I don't know if this is correct or not, but perhaps this
indicates the origin of the issue.

I've attached the input file and output file below. Thanks for any
insight into what's going on here and how to get proper m06-2x
results.


Adam

start
scratch_dir /mnt/lustre_scratch_2012/pelze1aCMICH9
echo
memory 1200 mb
geometry bridge1 noautoz noautosym units angstroms
C -1.269808 0.784902 -1.932051
C -2.62542 0.489323 -2.090349
C -3.291588 -0.22736 -1.114172
C -2.604401 -0.65443 0.037477
C -1.261942 -0.35329 0.194067
C -0.586074 0.370308 -0.798368
H -0.750592 1.343133 -2.702751
H -0.724218 -0.672298 1.078888
O -4.610106 -0.518672 -1.259002
O -3.376948 -1.347046 0.916929
C -2.76287 -1.809241 2.11268
H -3.536266 -2.337048 2.664224
H -1.94063 -2.491621 1.884016
H -2.394065 -0.968454 2.705452
H -4.900412 -1.010022 -0.478869
C 0.877029 0.659564 -0.627351
C 3.685554 0.606039 -0.387715
C 5.073129 0.508699 -0.268154
C 2.947666 -0.531954 -0.71439
C 5.726983 -0.698946 -0.470821
H 5.639971 1.396804 -0.013435
C 3.600436 -1.751878 -0.921283
C 4.978637 -1.831026 -0.799637
H 6.803933 -0.762214 -0.375785
H 5.473009 -2.78182 -0.962085
O 1.594321 -0.555607 -0.85702
H -3.172949 0.80934 -2.968799
H 3.205102 1.561539 -0.22685
H 1.087862 1.012898 0.387837
H 3.006689 -2.621801 -1.176405
H 1.20427 1.421195 -1.341739
end

BASIS "ao basis"  spherical PRINT
H library 6-311++G**
C library 6-311++G**
O library 6-311++G**
END

driver
maxiter 1000
end
set geometry bridge1
charge 0

property
nbofile
end

dft
mult 1
xc xm06-2x cm06-2x
grid xfine
  1. convergence damp 95 ncydp 600
vectors input atomic output CS_singlet.movecs
convergence energy 1.000000E-07
convergence density 1.000000E-05
convergence gradient 1E-05
iterations 1600
mulliken
end
task dft optimize
  1. task dft property

______

I'm not sure how to attach the full output file as on the list-serve, but here is the geometry after five optimization steps. You will see that it is indeed very strange. 

         --------

         Step   5

         --------





                        Geometry "geometry" -> "bridge1"

                        --------------------------------



Output coordinates in angstroms (scale by  1.889725989 to convert to a.u.)



 No.       Tag          Charge          X              Y              Z

---- ---------------- ---------- -------------- -------------- --------------

   1 C                    6.0000    -1.52251588     1.38924494    -1.71773538

   2 C                    6.0000    -3.04535544     1.03402489    -1.96525149

   3 C                    6.0000    -3.79905369     0.22636978    -0.83028863

   4 C                    6.0000    -3.02805455    -0.30600565     0.54448526

   5 C                    6.0000    -1.51720012     0.03215947     0.85997062

   6 C                    6.0000    -0.94412141     1.02800021    -0.26886746

   7 H                    1.0000    -0.75530924     2.04148024    -2.58560511

   8 H                    1.0000    -0.74234131    -0.34691887     1.89395730

   9 O                    8.0000    -5.36548871    -0.16212587    -0.90283064

  10 O                    8.0000    -4.01972372    -1.10845730     1.51869590

  11 C                    6.0000    -3.05197476    -1.53691917     2.86399246

  12 H                    1.0000    -3.99546604    -2.14550941     3.46852296

  13 H                    1.0000    -2.05387530    -2.24196059     2.35892594

  14 H                    1.0000    -2.60090108    -0.39448789     3.35516578

  15 H                    1.0000    -5.44219033    -0.71242023     0.09910977

  16 C                    6.0000     0.62404803     1.31674605    -0.09946475

  17 C                    6.0000     3.40686357     1.22948228     0.12730032

  18 C                    6.0000     4.99921540     1.12357963     0.26764245

  19 C                    6.0000     2.67177416    -0.20729126    -0.27667005

  20 C                    6.0000     5.74865423    -0.28996436     0.02893383

  21 H                    1.0000     5.64196706     2.21489604     0.57512181

  22 C                    6.0000     3.36027785    -1.61851010    -0.51070624

  23 C                    6.0000     4.93477155    -1.66474885    -0.35976522

  24 H                    1.0000     7.03513727    -0.36162716     0.14266787

  25 H                    1.0000     5.48573620    -2.80959551    -0.56132056

  26 O                    8.0000     1.11718864    -0.31008067    -0.46710013

  27 H                    1.0000    -3.73778650     1.40992221    -2.98488063

  28 H                    1.0000     2.80849014     2.37489588     0.31529663

  29 H                    1.0000     0.94902010     1.55660438     1.15648684

  30 H                    1.0000     2.50467796    -2.57392023    -0.81245292

  31 H                    1.0000     1.09952147     2.09238246    -1.05205279

Forum Vet
2:29:52 PM PST - Wed, Nov 14th 2012
Have you looked at the orbitals from the first DFT calculation? I wonder if they are in the wrong order due to a bad intial guess.

One good way to get a good initial guess is to do the following:

scf

 vectors input atomic output CS_singlet.movecs

 maxiter 100

end

task scf ignore



dft

mult 1

xc xm06-2x cm06-2x

grid xfine

convergence damp 95 ncydp 600

vectors input CS_singlet.movecs output CS_singlet.movecs

convergence energy 1.000000E-07

convergence density 1.000000E-05

convergence gradient 1E-05

iterations 1600  

mulliken 

end

task dft optimize

Bert


Quote:Adampelzer Nov 14th 6:17 pm
Hello. I'm trying to get NWChem 6.0 to reproduce some G09 results so
that I can use some features of NWChem that don't exist on G09 (namely
nudged elastic band) and am
running into some issues.

When I include the line
grid lebedev 90 14 ssf euler
in the input file, I can can get good agreement with G09 (0.004 eV) on
the total energy using PBE. However, when I use m06-2X, I get very
different results for the total energy (54 Hartrees lower in NWChem)
and the geometry that it's relaxing to is very odd (the bond lengths
of all atoms are much larger). It seems that I am calling for m06-2x
properly ( xc xm06-2x cm06-2x), as it's properly specified in the
output file

             XC Information

             --------------

                       M06-2X Exchange Functional  1.000

                     M06-2X Correlation Potential  1.000 non-local

But it also includes the line:
cm06-2x uses PW91LDA as defaults.These defaults cannot be changed
above. I don't know if this is correct or not, but perhaps this
indicates the origin of the issue.

I've attached the input file and output file below. Thanks for any
insight into what's going on here and how to get proper m06-2x
results.


Adam

start
scratch_dir /mnt/lustre_scratch_2012/pelze1aCMICH9
echo
memory 1200 mb
geometry bridge1 noautoz noautosym units angstroms
C -1.269808 0.784902 -1.932051
C -2.62542 0.489323 -2.090349
C -3.291588 -0.22736 -1.114172
C -2.604401 -0.65443 0.037477
C -1.261942 -0.35329 0.194067
C -0.586074 0.370308 -0.798368
H -0.750592 1.343133 -2.702751
H -0.724218 -0.672298 1.078888
O -4.610106 -0.518672 -1.259002
O -3.376948 -1.347046 0.916929
C -2.76287 -1.809241 2.11268
H -3.536266 -2.337048 2.664224
H -1.94063 -2.491621 1.884016
H -2.394065 -0.968454 2.705452
H -4.900412 -1.010022 -0.478869
C 0.877029 0.659564 -0.627351
C 3.685554 0.606039 -0.387715
C 5.073129 0.508699 -0.268154
C 2.947666 -0.531954 -0.71439
C 5.726983 -0.698946 -0.470821
H 5.639971 1.396804 -0.013435
C 3.600436 -1.751878 -0.921283
C 4.978637 -1.831026 -0.799637
H 6.803933 -0.762214 -0.375785
H 5.473009 -2.78182 -0.962085
O 1.594321 -0.555607 -0.85702
H -3.172949 0.80934 -2.968799
H 3.205102 1.561539 -0.22685
H 1.087862 1.012898 0.387837
H 3.006689 -2.621801 -1.176405
H 1.20427 1.421195 -1.341739
end

BASIS "ao basis"  spherical PRINT
H library 6-311++G**
C library 6-311++G**
O library 6-311++G**
END

driver
maxiter 1000
end
set geometry bridge1
charge 0

property
nbofile
end

dft
mult 1
xc xm06-2x cm06-2x
grid xfine
  1. convergence damp 95 ncydp 600
vectors input atomic output CS_singlet.movecs
convergence energy 1.000000E-07
convergence density 1.000000E-05
convergence gradient 1E-05
iterations 1600
mulliken
end
task dft optimize
  1. task dft property

______

I'm not sure how to attach the full output file as on the list-serve, but here is the geometry after five optimization steps. You will see that it is indeed very strange. 

         --------

         Step   5

         --------





                        Geometry "geometry" -> "bridge1"

                        --------------------------------



Output coordinates in angstroms (scale by  1.889725989 to convert to a.u.)



 No.       Tag          Charge          X              Y              Z

---- ---------------- ---------- -------------- -------------- --------------

   1 C                    6.0000    -1.52251588     1.38924494    -1.71773538

   2 C                    6.0000    -3.04535544     1.03402489    -1.96525149

   3 C                    6.0000    -3.79905369     0.22636978    -0.83028863

   4 C                    6.0000    -3.02805455    -0.30600565     0.54448526

   5 C                    6.0000    -1.51720012     0.03215947     0.85997062

   6 C                    6.0000    -0.94412141     1.02800021    -0.26886746

   7 H                    1.0000    -0.75530924     2.04148024    -2.58560511

   8 H                    1.0000    -0.74234131    -0.34691887     1.89395730

   9 O                    8.0000    -5.36548871    -0.16212587    -0.90283064

  10 O                    8.0000    -4.01972372    -1.10845730     1.51869590

  11 C                    6.0000    -3.05197476    -1.53691917     2.86399246

  12 H                    1.0000    -3.99546604    -2.14550941     3.46852296

  13 H                    1.0000    -2.05387530    -2.24196059     2.35892594

  14 H                    1.0000    -2.60090108    -0.39448789     3.35516578

  15 H                    1.0000    -5.44219033    -0.71242023     0.09910977

  16 C                    6.0000     0.62404803     1.31674605    -0.09946475

  17 C                    6.0000     3.40686357     1.22948228     0.12730032

  18 C                    6.0000     4.99921540     1.12357963     0.26764245

  19 C                    6.0000     2.67177416    -0.20729126    -0.27667005

  20 C                    6.0000     5.74865423    -0.28996436     0.02893383

  21 H                    1.0000     5.64196706     2.21489604     0.57512181

  22 C                    6.0000     3.36027785    -1.61851010    -0.51070624

  23 C                    6.0000     4.93477155    -1.66474885    -0.35976522

  24 H                    1.0000     7.03513727    -0.36162716     0.14266787

  25 H                    1.0000     5.48573620    -2.80959551    -0.56132056

  26 O                    8.0000     1.11718864    -0.31008067    -0.46710013

  27 H                    1.0000    -3.73778650     1.40992221    -2.98488063

  28 H                    1.0000     2.80849014     2.37489588     0.31529663

  29 H                    1.0000     0.94902010     1.55660438     1.15648684

  30 H                    1.0000     2.50467796    -2.57392023    -0.81245292

  31 H                    1.0000     1.09952147     2.09238246    -1.05205279[/quote]

Forum Vet
5:43:59 PM PST - Wed, Nov 14th 2012
Quote:Adampelzer Nov 14th 10:17 am

of all atoms are much larger). It seems that I am calling for m06-2x
properly ( xc xm06-2x cm06-2x), as it's properly specified in the
output file

             XC Information

             --------------

                       M06-2X Exchange Functional  1.000

                     M06-2X Correlation Potential  1.000 non-local



Adam
The complete definition of M06-2X contains 54% of HF exchange,
therefore your input line for the xc functional should change to
xc cm06-2x xm06-2x hfexch 0.54

However, a simpler solution would be to use the
combined m06-2x keyword (that properly contains hf exchange)
xc  m06-2x


Cheers, Edo

Clicked A Few Times
7:18:23 PM PST - Wed, Nov 14th 2012
I checked the orbital energies from the NWChem M06-2x calculation and there was no "hole" such that an unoccupied orbital was lower in energy than an occupied one. I've seen that happen in these calculations before, and they've never resulted in a 54 Hartree difference.

I also did as you suggested and did a HF calculation to get a better set of guess orbitals. The total energy for this came out to -763.057301433098 Hartrees, which is nearly identical to the G09 results of -763.057667241 Hartrees. However, using these orbitals as a guess for the M06-2X calculation led to an energy of -713.1746899070 Hartrees after the first SCF iteration. This is very close to the energy of the initial SCF calculation with M06-2x where I didn't use HF orbitals as a guess: -713.17469465. This once again leads me to believe that the issue is with the implementation of the M06-2x functional.

Forum Vet
10:04:41 AM PST - Thu, Nov 15th 2012
See the note from Edo on the proper definition of the M06-2X functional.

Bert


Quote:Adampelzer Nov 15th 2:18 am
I checked the orbital energies from the NWChem M06-2x calculation and there was no "hole" such that an unoccupied orbital was lower in energy than an occupied one. I've seen that happen in these calculations before, and they've never resulted in a 54 Hartree difference.

I also did as you suggested and did a HF calculation to get a better set of guess orbitals. The total energy for this came out to -763.057301433098 Hartrees, which is nearly identical to the G09 results of -763.057667241 Hartrees. However, using these orbitals as a guess for the M06-2X calculation led to an energy of -713.1746899070 Hartrees after the first SCF iteration. This is very close to the energy of the initial SCF calculation with M06-2x where I didn't use HF orbitals as a guess: -713.17469465. This once again leads me to believe that the issue is with the implementation of the M06-2x functional.

Forum Regular
4:30:39 PM PST - Thu, Nov 15th 2012
Hi Adam,
I can confirm that the right setting is as Edo suggested:

  xc m06-2x

The first few iterations of your job with this functional specification give energies of

                iter       energy          gnorm     gmax       time

            ----- ------------------- --------- --------- --------

                1     -767.2664469117  1.94D+00  2.15D-01    528.7

                2     -767.4315112635  7.80D-01  8.62D-02   1217.0

I trust that this is more what you were looking for than the -713 Hartree you got before.

Best wishes, Huub

Clicked A Few Times
5:20:52 AM PST - Fri, Nov 16th 2012
Quote:Edoapra Nov 15th 12:43 am
Quote:Adampelzer Nov 14th 10:17 am

of all atoms are much larger). It seems that I am calling for m06-2x
properly ( xc xm06-2x cm06-2x), as it's properly specified in the
output file

             XC Information

             --------------

                       M06-2X Exchange Functional  1.000

                     M06-2X Correlation Potential  1.000 non-local



Adam
The complete definition of M06-2X contains 54% of HF exchange,
therefore your input line for the xc functional should change to
xc cm06-2x xm06-2x hfexch 0.54

However, a simpler solution would be to use the
combined m06-2x keyword (that properly contains hf exchange)
xc  m06-2x


Cheers, Edo


Thanks, this seems to fix things. Strangely though, I'm getting a 0.02eV difference in the total energy, which is an order of magnitude larger than I got for PBE. If anyone has any ideas about why this might be, I'd appreciate it.

Thanks for the help thus far everyone.

Adam

Forum Regular
1:08:06 PM PST - Fri, Nov 16th 2012
Hi Adam,
The M06-2X functional is a meta-GGA. As far as I know meta-GGAs tend to generate quantities that are harder to integrate accurately on a grid. Hence these calculations tend to be much more sensitive to the quadrature technology being used. This might explain the 10 times larger deviations between codes for meta-GGA functionals relative to GGA.
Best wishes,
Huub

Clicked A Few Times
10:24:02 PM PST - Sun, Nov 18th 2012
Quote:Huub Nov 16th 8:08 pm
Hi Adam,
The M06-2X functional is a meta-GGA. As far as I know meta-GGAs tend to generate quantities that are harder to integrate accurately on a grid. Hence these calculations tend to be much more sensitive to the quadrature technology being used. This might explain the 10 times larger deviations between codes for meta-GGA functionals relative to GGA.
Best wishes,
Huub


Is this the sort of thing that I should expect to be improved by using "grid xfine"?

Forum Regular
11:01:06 AM PST - Mon, Nov 19th 2012
Hi Adam,
I had a brief chat with Edo and the message is that with Meta-GGAs you need to work harder to get good answers. So, yes, the results should improve with going to the finest grid available, but you may have to do more than that. For example you may have to tighten the SCF convergence as well.
Best wishes,
Huub


Forum >> NWChem's corner >> Running NWChem