Strange results with m06-2x...

From NWChem

Viewed 1357 times, With a total of 9 Posts
Jump to: navigation, search

Clicked A Few Times
Threads 6
Posts 14
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

  • Bert Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Forum Vet
Threads 5
Posts 598
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
Threads 4
Posts 715
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
Threads 6
Posts 14
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.

  • Bert Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Forum Vet
Threads 5
Posts 598
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.

  • Huub Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Forum Regular
Threads 2
Posts 182
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
Threads 6
Posts 14
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

  • Huub Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Forum Regular
Threads 2
Posts 182
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
Threads 6
Posts 14
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"?

  • Huub Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Forum Regular
Threads 2
Posts 182
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



Who's here now Members 0 Guests 0 Bots/Crawler 1


AWC's: 2.5.10 MediaWiki - Stand Alone Forum Extension
Forum theme style by: AWC