Implementing new functionals

[August Melcher]

Hello,

I would like to implement the ωB97X-V functional as presented in this paper: http://www.ncbi.nlm.nih.gov/pubmed/24430168

I am relatively new to CP2K so I apologize if this is a trivial question, but would someone mind assisting me in implementing this functional?

Thanks,

August

[hut...@chem.uzh.ch]

Hi

I had a quick look at this functional and it seems to me there
are three non-trivial things to do for the implementation in
CP2K.

1) add a new parametrization to the BECKE97 functional
follow the lead from other parametrizations
2) add a new potential type in the HF section
follow the lead from other mixed potential types in the code
3) check if the implemented rVV10 functional can be used instead
of the original VV10. Only rVV10 is implemented in CP2K

Your input will look something like this:

&XC
&XC_FUNCTIONAL
&BECKE97
PARAMETRIZATION XXXXX
SCALE_X 0.0
SCALE_C 1.0
&END PBE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-3
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE SHORTRANGE-LONGRANGE
OMEGA 0.11
&END
&MEMORY
MAX_MEMORY 10
&END
FRACTION 0.25
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL NON_LOCAL
&NON_LOCAL
TYPE RVV10
PARAMETERS 6.3 0.0093
VERBOSE_OUTPUT
KERNEL_FILE_NAME ../rVV10_kernel_table.dat
CUTOFF 80
&END NON_LOCAL
&END vdW_POTENTIAL
&END XC


best regards

Juerg Hutter

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com
Date: 03/04/2014 09:30PM
Subject: [CP2K:5031] Implementing new functionals

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/groups/opt_out.

[hut...@chem.uzh.ch]

Hi

I checked again and realized that by rearranging the HFX terms
you get the 'MIX_CL' potential type. With the following
HF section you should get the HFX part of the functional
(parameters Cx and omega). MIX_CL_TRUNC has to be used in periodic systems.

&HF
&SCREENING
EPS_SCHWARZ 1.0E-8
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE MIX_CL (MIX_CL_TRUNC)
SCALE_COULOMB 1-Cx
SCALE_LONGRANGE Cx
OMEGA XXXX

&END
&MEMORY
MAX_MEMORY 10
&END

FRACTION 1.00
&END

Seems that all that is left are the B97 parameters and to hope
that the rVV10 is not changing the performance of the functional.

regards

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----Jürg Hutter/at/UZH wrote: -----
To: cp...@googlegroups.com
From: Jürg Hutter/at/UZH
Date: 03/05/2014 09:52AM
Subject: Re: [CP2K:5031] Implementing new functionals

[August Melcher]

Thanks for the help, but I am still confused about how to add a new parameterization for the BECKE97 functional. You said I am supposed to follow the lead from other parameterizations, so I assume that means I should look at the files where the functional is parameterized and base it off of those. However I am not sure where to find these files and how to make sure that once I have created the new parameterization file the input file can recognize it. Can I specify a file path after PARAMETERIZATION?

[Ralph Koitz]

Hi,

Starting from an up-to-date cp2k source tree, you'd have to make a few modifications to xc_b97.F.

First, add your parameter values as a new array and give it a name like the other two that are already there (~line 54).

You will also need to give your parameter set a number, ~line 602 in input_constants.F, and "USE" that in xc_b97.F (~line 28).

Then you can go through xc_b97.F, searching for all occurrences of "xc_b97_orig", and add appropriate "CASE" instructions for your own parametrization. Mostly it should be enough to copy-paste and switch out the variable names, and perhaps the lit. reference.

Finally, in order to make your parametrization available from the input file, you will need to appropriately modify lines 580--584 in input_cp2k_xc.F, simply adding one more element to the enum_* arrays and adapting the "description". (also "USE" the new input constant at the beginning of the file.)

I hope that's all. Then recompile and hope for the best :)

Ralph

--

You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.

For more options, visit https://groups.google.com/d/optout.

[hut...@chem.uzh.ch]

Hi

attached are the files with the new parameters. I hope you have
an independent source for testing. The rVV10 is slightly different
from the VV10 and I don't know exactly how the parameters should be
translated.
Be careful when using this functional in periodic boundary condition.
You have to use the truncated long range exchange.

regards

Juerg

the input section:

&XC
&XC_FUNCTIONAL
&BECKE97
PARAMETRIZATION wB97X-V
SCALE_X 1.0
SCALE_C 1.0
&END
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-8 #to be determined
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE MIX_CL
SCALE_COULOMB 0.833
SCALE_LONGRANGE 0.167
OMEGA 0.30
&END
&MEMORY
MAX_MEMORY 10 #depends on your computer
&END

&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL NON_LOCAL
&NON_LOCAL
TYPE RVV10

PARAMETERS 6.0 0.01 #can they be taken as is?
VERBOSE_OUTPUT
KERNEL_FILE_NAME ../rVV10_kernel_table.dat
CUTOFF 80 #to be determined

&END NON_LOCAL
&END vdW_POTENTIAL
&END XC

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 03/11/2014 07:41PM
Subject: [CP2K:5058] Re: Implementing new functionals

[hut...@chem.uzh.ch]

Hi

just realized that I interchanged some parameters.
new file xc_b97.F

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com
Date: 03/11/2014 07:41PM
Subject: [CP2K:5058] Re: Implementing new functionals

[hut...@chem.uzh.ch]

Hi

any feedback on this?

thank you

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----Jürg Hutter/at/UZH wrote: -----
To: cp...@googlegroups.com
From: Jürg Hutter/at/UZH
Date: 03/12/2014 10:10AM

Subject: [CP2K:5058] Re: Implementing new functionals

Hi

just realized that I interchanged some parameters.
new file xc_b97.F

Juerg
 
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com
Date: 03/11/2014 07:41PM
Subject: [CP2K:5058] Re: Implementing new functionals

Thanks for the help, but I am still confused about how to add a new parameterization for the BECKE97 functional. You said I am supposed to follow the lead from other parameterizations, so I assume that means I should look at the files where the functional is parameterized and base it off of those. However I am not sure where to find these files and how to make sure that once I have created the new parameterization file the input file can recognize it. Can I specify a file path after PARAMETERIZATION?

On Tuesday, March 4, 2014 12:30:36 PM UTC-8, August Melcher wrote:

  --
 You received this message because you are subscribed to the Google Groups "cp2k" group.
 To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
 To post to this group, send email to cp...@googlegroups.com.
 Visit this group at http://groups.google.com/group/cp2k.
 For more options, visit https://groups.google.com/d/optout.
 

[attachment "xc_b97.F" removed by Jürg Hutter/at/UZH]

[August Melcher]

Right now I'm getting the error message below

*****************************************************************************

 *** 13:34:28 ERRORL2 in input_enumeration_types:enum_c2i :: invalid value ***

 *** for enumeration:wB97X-V                                               ***

 *****************************************************************************

 *****************************************************************************

 *** 13:34:28 ERRORL2 in input_enumeration_types:enum_c2i :: invalid value ***

 *** for enumeration2:wB97X-V                                              ***

 *****************************************************************************

 Looking for words in the input similar to the unknown:

   'WB97X-V'

              CP2K failed to parse the input file.

              A full description of the input for this CP2K version

              can be generated using:

              cp2k.sopt --html-manual

              The manual for the latest version of CP2K is online available:

              http://manual.cp2k.org/trunk

              If this input was an input of a previous version

              of CP2K, you can try to convert it with --permissive-echo.

              However, this will just ignore the unknown keywords ...

Rank 0 [Mon Apr  7 13:34:29 2014] [c9-2c0s0n2] application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:

 CP2K| Abnormal program termination, stopped by process number 0

#0  0x2405E5D in _gfortran_backtrace at backtrace.c:258

#1  0x23E2290 in _gfortrani_backtrace_handler at compile_options.c:129

#2  0x26EEE5F in system

#3  0x26EEE0B in raise at pt-raise.c:41

#4  0x27028B0 in abort at abort.c:92

#5  0x25CE2F1 in MPID_Abort

#6  0x259FE22 in MPI_Abort

#7  0x1E74E24 in mpi_abort__

#8  0x6089A8 in __message_passing_MOD_mp_abort

#9  0x428D0D in MAIN__ at cp2k.F:0

_pmiu_daemon(SIGCHLD): [NID 04766] [c9-2c0s0n2] [Mon Apr  7 13:34:29 2014] PE RANK 0 exit signal Aborted

[hut...@chem.uzh.ch]

Hi

it runs in my implementation. I have committed the code to the
trunk version. There is also an example input
tests/QS/regtest-hybrid-4/wB97X-V.inp

I am rather sure that the functional is not yet correct.
Do you have a plan on how to test the correctness of the functional?

regards

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 04/07/2014 10:41PM
Subject: Re: [CP2K:5145] Re: Implementing new functionals

[August Melcher]

Alright, it's working now. I'm going to talk to Narbe Mardirossian (who developed the functional) about the best way to test it.

On Tuesday, March 4, 2014 12:30:36 PM UTC-8, August Melcher wrote:

[August Melcher]

We're going to start by comparing potential energy curves for the argon dimer. In the original computations, Narbe used different integration grids for the local xc functional (250,590) and VV10 (75,302). The only place I know to specify a grid is under the KIND section and there it seems that I can only specify one grid fineness. How would I go about specifying different grids for the two different parts of the functional?

On Tuesday, March 4, 2014 12:30:36 PM UTC-8, August Melcher wrote:

[hut...@chem.uzh.ch]

Hi

remember that the grids in CP2K have no connection to the grids
used in other QC codes.
You have to increase the energy cutoff to get better results.
The grids you can set in the KINDS sections are for the GAPW
method only (for the integration of the one center terms).

regards

Juerg Hutter

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 04/10/2014 09:13PM
Subject: [CP2K:5163] Re: Implementing new functionals

[August Melcher]

I'm currently getting this error when I run with the aug-cc-pVQZ basis, but not with aug-cc-pVDZ. 

 ****************************************************************************

 *** 15:21:22 ERRORL2 in hfx_libint_wrapper:initialize_libderiv processor ***

 *** 0  :: err=-300 the angular momentum needed exceeds the value assumed ***

 *** when configuring libderiv                                            ***

 ****************************************************************************

My input file looks like this:

&FORCE_EVAL

  METHOD Quickstep

  &DFT

    BASIS_SET_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/EMSL_BASIS_SETS

    POTENTIAL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/POTENTIAL

    &MGRID

      CUTOFF 1000

      NGRIDS 5

      REL_CUTOFF 50

    &END MGRID

    &QS

      METHOD GAPW

      EPS_PGF_ORB 1.0E-12

      EPS_FILTER_MATRIX 0.0e0

    &END QS

    &POISSON

      PERIODIC NONE

      PSOLVER MT

    &END

    &SCF

      EPS_SCF 1.0E-8

      SCF_GUESS ATOMIC

      MAX_SCF 400

      &OT ON

      &END

    &END SCF

    &XC

      &XC_FUNCTIONAL

        &BECKE97

            PARAMETRIZATION  wB97X-V

            SCALE_X 1.0

            SCALE_C 1.0

        &END

      &END XC_FUNCTIONAL

      &HF

        &SCREENING

           EPS_SCHWARZ 1.0E-14

        &END

        &INTERACTION_POTENTIAL

           POTENTIAL_TYPE MIX_CL

           SCALE_COULOMB   0.833

           SCALE_LONGRANGE 0.167

           OMEGA 0.30

        &END

        &MEMORY

          MAX_MEMORY 512

        &END

        #FRACTION 0.167

      &END

      &vdW_POTENTIAL

         DISPERSION_FUNCTIONAL NON_LOCAL

         &NON_LOCAL

           TYPE RVV10

           PARAMETERS 6.3 0.0093

           VERBOSE_OUTPUT

           KERNEL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/rVV10_kernel_table.dat

           CUTOFF  80

         &END NON_LOCAL

      &END vdW_POTENTIAL

   &END XC

  &END DFT

  &SUBSYS

    &CELL

      ABC 15.0 15.0 15.0

      PERIODIC NONE

    &END CELL

    &COORD

Ar       1.7500 0.0000 0.0000

Ar       -1.7500 0.0000 0.0000

    &END COORD

    &KIND Ar

      BASIS_SET aug-cc-pVDZ

      LEBEDEV_GRID 590

      RADIAL_GRID 250

      POTENTIAL ALL

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

&GLOBAL

  PROJECT 1.5

  PRINT_LEVEL LOW

  RUN_TYPE ENERGY

  &TIMINGS

    THRESHOLD 0.000000001

  &END

&END GLOBAL

[Michael Banck]

Hi,

On Tue, Apr 22, 2014 at 03:26:55PM -0700, August Melcher wrote:
> I'm currently getting this error when I run with the aug-cc-pVQZ basis, but
> not with aug-cc-pVDZ.
>
> ****************************************************************************
> *** 15:21:22 ERRORL2 in hfx_libint_wrapper:initialize_libderiv processor
> ***
> *** 0 :: err=-300 the angular momentum needed exceeds the value assumed
> ***
> *** when configuring libderiv
> ***
> ****************************************************************************

I can reproduce this with stock Debian/Ubuntu packages and the attached reduced
input file, which required much less memory.

There are two things here:

1. the HFX module appears to initialize libderiv, which should only be
required for forces, but your job is a single-point calculation. I was
able to reproduce a failure with PSI3 when running a geometry
optimization of Ar2, but not for a single-point run, so I think CP2K
should be able to do it in principle as well.

2. The default --with-libderiv-max-am1 configure option for the libint
library is 3, and the aug-cc-pVQZ basis for Ar has an angular momentum
of 4. So you can recompile libint with higher --with-libderiv-max-am1
(and maybe --with-libderiv-max-am2, just to be sure) options and link
CP2K against it.


Michael

[August Melcher]

I'm working on producing a neon dimer potential energy curve, and I keep getting results that look like this. I don't seem to understand how coordinates work because I keep seeing the binding part of the curve at what I expected to be the largest separation. Here is a sample input of what I thought would be a 1.5A separation. The difference in the x coordinates is increased correspondingly for the rest of the points. See attached for what the curve looks like.

&FORCE_EVAL

  METHOD Quickstep

  &DFT

    BASIS_SET_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/EMSL_BASIS_SETS

    POTENTIAL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/POTENTIAL

    &MGRID

      CUTOFF 1000

      NGRIDS 5

      REL_CUTOFF 50

    &END MGRID

    &QS

      METHOD GAPW

      EPS_PGF_ORB 1.0E-12

      EPS_FILTER_MATRIX 0.0e0

    &END QS

    &POISSON

      PERIODIC NONE

      PSOLVER MT

    &END

    &SCF

      EPS_SCF 1.0E-8

      SCF_GUESS ATOMIC

      MAX_SCF 40

      &OUTER_SCF

        EPS_SCF 1.0E-8

        MAX_SCF 50

      &END OUTER_SCF

      &OT ON

      &END

    &END SCF

    &XC

      &XC_FUNCTIONAL

        &BECKE97

            PARAMETRIZATION  wB97X-V

            SCALE_X 1.0

            SCALE_C 1.0

        &END

      &END XC_FUNCTIONAL

      &HF

        &SCREENING

           EPS_SCHWARZ 1.0E-14

        &END

        &INTERACTION_POTENTIAL

           POTENTIAL_TYPE MIX_CL

           SCALE_COULOMB   0.833

           SCALE_LONGRANGE 0.167

           OMEGA 0.30

        &END

        &MEMORY

          MAX_MEMORY 128

        &END

        #FRACTION 0.167

      &END

      &vdW_POTENTIAL

         DISPERSION_FUNCTIONAL NON_LOCAL

         &NON_LOCAL

           TYPE RVV10

           PARAMETERS 6.3 0.0093

           VERBOSE_OUTPUT

           KERNEL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/rVV10_kernel_table.dat

           CUTOFF  80

         &END NON_LOCAL

      &END vdW_POTENTIAL

   &END XC

  &END DFT

  &SUBSYS

    &CELL

      ABC 20.0 15.0 15.0

      PERIODIC NONE

    &END CELL

    &COORD

Ne       13.25 0.0000 0.0000

Ne       11.75 0.0000 0.0000

    &END COORD

    &KIND Ne

      BASIS_SET aug-cc-pVQZ

      LEBEDEV_GRID 590

      RADIAL_GRID 250

      POTENTIAL ALL

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

&GLOBAL

  PROJECT wb97xv15

[Matt W]

Hi,

I think the probable cause is the MT Poisson solver - it needs box length 2x the maximum distance between electron density in any given direction - so for your box you can only safely go out to ~10 A. 

Probably your actual binding curve is at the smaller distance, but you can't see it for the huge noise from the spurious results at long distance.

Matt

[August Melcher]

So should I use a different solver and rerun? Just zooming in on the plot doesn't help.

[Matt W]

I think your setup looks OK for distances ~< 10 A. Why you appear to get no binding I'm not sure.

Did you check with a known to be working functional?

Matt

[August Melcher]

Yes, I'm getting the exact same problem for B97

&FORCE_EVAL

  METHOD Quickstep

  &DFT

    BASIS_SET_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/EMSL_BASIS_SETS

    POTENTIAL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/POTENTIAL

    &MGRID

      CUTOFF 1000

      NGRIDS 5

      REL_CUTOFF 50

    &END MGRID

    &QS

      METHOD GAPW

      EPS_PGF_ORB 1.0E-12

      EPS_FILTER_MATRIX 0.0e0

    &END QS

    &SCF

      EPS_SCF 1.0E-8

      SCF_GUESS ATOMIC

      MAX_SCF 40

      &OUTER_SCF

        EPS_SCF 1.0E-8

        MAX_SCF 50

      &END OUTER_SCF

      &OT ON

      &END

    &END SCF

    &XC

      &XC_FUNCTIONAL

        &BECKE97

            PARAMETRIZATION  ORIG

        &END

      &END XC_FUNCTIONAL

   &END XC

  &END DFT

  &SUBSYS

    &CELL

      ABC 20.0 15.0 15.0

      PERIODIC NONE

    &END CELL

    &COORD

Ne       13.25 0.0000 0.0000

Ne       11.75 0.0000 0.0000

    &END COORD

    &KIND Ne

      BASIS_SET aug-cc-pVQZ

      LEBEDEV_GRID 590

      RADIAL_GRID 250

      POTENTIAL ALL

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

&GLOBAL

  PROJECT b9715

  PRINT_LEVEL LOW

  RUN_TYPE ENERGY

  &TIMINGS

    THRESHOLD 0.000000001

  &END

&END GLOBAL

On Tuesday, March 4, 2014 12:30:36 PM UTC-8, August Melcher wrote:

[Matt W]

Where'd you get your reference data? I just quickly and a B97 curve with a popular commercial QC code, and got a purely repulsive profile for B97 (i.e. pretty similar to your CP2K result).

Matt

[August Melcher]

My reference data is from Narbe Mardirossian, who developed the wB97X-V functional. He used Q-Chem.

[Matt W]

Hi again.

Well, maybe somehow apples and oranges are being compared?

I attach the PES I get for Ne dimer using Guassian 09 and the B97 functional (keyword B98, so possibly a revision?), with aug-cc-pVTZ basis (no BSSE correction applied).

I'll leave it for other more knowledgeable people to comment further.

Matt

[hut...@chem.uzh.ch]

Hi

a few more points:
- if possible I would suggest to test the different parts of the
functional individually (GGA + HFX + vdW)
- test first for the total energy using a simple system (He or Ne)
- remember that the VV10 part is in CP2K rVV10. In addition, this is
only calculated for the 'soft' density in a GAPW calculation.
You will not be able to compare total energies, only differences.

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com

From: Matt W
Sent by: cp...@googlegroups.com
Date: 05/01/2014 07:11PM
Subject: [CP2K:5235] Re: Implementing new functionals

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/d/optout.

[attachment "B97_Ne2.png" removed by Jürg Hutter/at/UZH]

[August Melcher]

So am I getting correct behavior for B97? The fact that the B97 and wB97X-V curves look so similar made me think that there was something wrong with the input file other than the GGA/HFX/vdW parts.

[August Melcher]

For reference, here are both PE curves on the same plot (as well as the input files again).

[Michael Banck]

Hi,

On Thu, May 01, 2014 at 01:10:42AM -0700, August Melcher wrote:
> So should I use a different solver and rerun? Just zooming in on the
> plot doesn't help.

First of all, anything beyond 5-7 Angstrom is probably not very
interesting, and the latest plot you sent (b97s.jpg) is again unreadable
due to the artefacts beyond 15 Angstrom.

What are the units in the second plot (the binding one, image.png), and
how did you make it?

I ran a scan with the PSI4 (Beta5) package using the (probably similar?)
wB97X-D functional and a aug-cc-pVQZ basis. This is the potential curve
I get:

R (Ang) E (a.u.)
-----------------------------------------------------
2.50 -257.90034199
2.75 -257.90174824
3.00 -257.90232610
3.25 -257.90253327
3.50 -257.90257489
3.75 -257.90255286
4.00 -257.90252022
4.25 -257.90249566
4.50 -257.90247981
4.75 -257.90246941
5.00 -257.90246231
5.50 -257.90245369
6.00 -257.90244979
6.50 -257.90244770
7.00 -257.90244637
7.50 -257.90244526

This translates to a binding energy of 0.08 kcal/mol which (assuming
identical units) is slightly higher than what your second plot shows for
B97.


Michael

[hut...@chem.uzh.ch]

Hi

here some numbers for He atom with the aug-cc-pvTZ basis:

Reference CP2K
HF -2.861183 -2.861183
BLYP -2.906512 -2.906512
PBE0 -2.894685 -2.894688
B97D -2.916114 -2.916115

wB97X(no V!) ????????? -2.920003

maybe this helps.

Juerg

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com

From: Michael Banck
Sent by: cp...@googlegroups.com
Date: 05/02/2014 11:35PM
Subject: Re: [CP2K:5238] Re: Implementing new functionals

[August Melcher]

Could someone please clarify whether I should be seeing binding for the neon dimer using the B97 functional? Because I currently am not seeing any but according to the Q-Chem reference data, I should be.

[August Melcher]

This may be getting ahead of myself a bit, but what pseudopotential should I use with wB97X-V? Should I make one using the ATOM program or will one of the existing pseudopotentials work fine?

[hut...@chem.uzh.ch]

Hi

range separated hybrids are not implemented in the atomic code.
We can therefore not optimize pseudopotentials for such functionals.
From experience I would suggest to use pseudpotentials for PBE0
(1. choice) or even PBE (2. choice), but don't mix them.

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com

From: August Melcher
Sent by: cp...@googlegroups.com

Date: 05/16/2014 09:17PM
Subject: [CP2K:5290] Re: Implementing new functionals

This may be getting ahead of myself a bit, but what pseudopotential should I use with wB97X-V? Should I make one using the ATOM program or will one of the existing pseudopotentials work fine?

[August Melcher]

Alright, thanks. Are the PBE0 GTH pseudopotentials contained somewhere in CP2K or should I be looking elsewhere for them?

[hut...@chem.uzh.ch]

No, they are not in the distribution.
I have a few generated for our own use (O,H).

regards

Juerg


--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491

Institut für Chemie C              FAX   : ++41 44 635 6838

Universität Zürich                  E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 05/19/2014 09:35PM
Subject: [CP2K:5295] Re: Implementing new functionals

Alright, thanks. Are the PBE0 GTH pseudopotentials contained somewhere in CP2K or should I be looking elsewhere for them?

[August Melcher]

I'm sorry if this is gauche, but would you mind letting me test your PBE0 O and H pseudopotentials for comparison versus PBE pseudopotentionals to determine if it's worth generating more?

[hut...@chem.uzh.ch]

Hi

the two pseudos are attached. They have been used and published
(supplemental material) in

Bulk liquid water at ambient temperature and pressure from MP2 theory,
Mauro Del Ben, Mandes Schönherr, Jürg Hutter and Joost VandeVondele,
Journal of Physical Chemistry Letters 4, 3753-3759 (2013)
http://dx.doi.org/10.1021/jz401931f

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C              FAX   : ++41 44 635 6838
Universität Zürich                  E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 05/21/2014 08:36PM
Subject: [CP2K:5302] Re: Implementing new functionals

I'm sorry if this is gauche, but would you mind letting me test your PBE0 O and H pseudopotentials for comparison versus PBE pseudopotentionals to determine if it's worth generating more?

[August Melcher]

After computing the total electronic energy for a neon atom using wB97X-V and comparing it to QChem reference results, it seems that the DFT XC energy in CP2K is not being computed correctly. Does anyone have any idea what might be wrong here? Perhaps the parametrization?

CP2K Results:

  Overlap energy of the core charge distribution:               0.00000000000000

  Self energy of the core charge distribution:               -148.47094303888323

  Core Hamiltonian energy:                                      0.61509680993492

  Hartree energy:                                              14.01686743627672

  Exchange-correlation energy:                                 -4.11859163513813

  Hartree-Fock Exchange energy:                                -3.38266780588875

  Dispersion energy:                                            0.02876776737550

  GAPW| Exc from hard and soft atomic rho1:                    -6.47114106845743

  GAPW| local Eh = 1 center integrals:                         17.39257016010844

  Total energy:                                              -130.39004137467197

  outer SCF iter =    2 RMS gradient =   0.85E-08 energy =       -130.3900413747

  outer SCF loop converged in   2 iterations or  435 steps

 ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):             -130.390041374671966

QChem results (aug-cc-pVQZ/250,590):

 Nonlocal correlation =     0.0441808008

    7    -128.9329200166      6.16E-09 Convergence criterion met

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

 One-Electron    Energy =  -182.3744736608

 Total Coulomb   Energy =    65.8986384358

 Alpha Exchange  Energy =    -1.6882204462

 Beta  Exchange  Energy =    -1.6882204462

 DFT   Exchange  Energy =    -8.6509462294

 DFT Correlation Energy =    -0.4296976698

 Nuclear Repu.   Energy =     0.0000000000

 Nuclear Attr.   Energy =  -311.0040993334

 Kinetic         Energy =   128.6296256727

&FORCE_EVAL

  METHOD Quickstep

  &DFT

    BASIS_SET_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/EMSL_BASIS_SETS

    POTENTIAL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/POTENTIAL

    &MGRID

      CUTOFF 1000

      NGRIDS 5

      REL_CUTOFF 50

    &END MGRID

    &QS

      METHOD GAPW

      EPS_PGF_ORB 1.0E-12

      EPS_FILTER_MATRIX 0.0e0

    &END QS

    &POISSON

      PERIODIC NONE

      PSOLVER MT

    &END

    &SCF

      EPS_SCF 1.0E-8

      SCF_GUESS ATOMIC

      MAX_SCF 400

      &OUTER_SCF

        EPS_SCF 1.0E-8

        MAX_SCF 50

      &END OUTER_SCF

      &OT ON

      &END

    &END SCF

    &XC

      &XC_FUNCTIONAL

        &BECKE97

            PARAMETRIZATION  wB97X-V

            SCALE_X 1.0

            SCALE_C 1.0

        &END

      &END XC_FUNCTIONAL

      &HF

        &SCREENING

           EPS_SCHWARZ 1.0E-14

        &END

        &INTERACTION_POTENTIAL

           POTENTIAL_TYPE MIX_CL

           SCALE_COULOMB   0.167

           SCALE_LONGRANGE 0.833

           OMEGA 0.30

        &END

        &MEMORY

          MAX_MEMORY 256

        &END

      &END

      &vdW_POTENTIAL

         DISPERSION_FUNCTIONAL NON_LOCAL

         &NON_LOCAL

           TYPE RVV10

           PARAMETERS 6.3 0.0093

           VERBOSE_OUTPUT

           KERNEL_FILE_NAME /global/u1/a/amelcher/soft/cp2k/cp2k/tests/QS/rVV10_kernel_table.dat

           CUTOFF  150

         &END NON_LOCAL

      &END vdW_POTENTIAL

   &END XC

  &END DFT

  &SUBSYS

    &CELL

      ABC 10.0 10.0 10.0

      PERIODIC NONE

    &END CELL

    &COORD

Ne       0.0000 0.0000 0.0000

    &END COORD

    &KIND Ne

      BASIS_SET aug-cc-pVQZ

      LEBEDEV_GRID 590

      RADIAL_GRID 250

      POTENTIAL ALL

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

&GLOBAL

  PROJECT wb97xv25

[hut...@chem.uzh.ch]

Hi

I think you need an additional

FRACTION 0.167

in the &HF Section

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 06/03/2014 02:02AM
Subject: [CP2K:5358] Re: Implementing new functionals

[August Melcher]

That wouldn't seem to be the case, it now gives 16.7% of the correct HF exchange energy. The DFT XC energy is unaffected.

  Overlap energy of the core charge distribution:               0.00000000000000

  Self energy of the core charge distribution:               -148.47094303888323

  Core Hamiltonian energy:                                     -0.18026539801860

  Hartree energy:                                              14.53040907815531

  Exchange-correlation energy:                                 -4.05872942140104

  Hartree-Fock Exchange energy:                                -0.55962243700625

  Dispersion energy:                                            0.02879894951182

  GAPW| Exc from hard and soft atomic rho1:                    -6.39896663215167

  GAPW| local Eh = 1 center integrals:                         17.52403430840716

  Total energy:                                              -127.58528459138651

  outer SCF iter =    1 RMS gradient =   0.93E-08 energy =       -127.5852845914

  outer SCF loop converged in   1 iterations or  365 steps

 ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):             -127.585284591386511

[hut...@chem.uzh.ch]

Hi

let's establish a common baseline. Using the attached input
I get the following PBE0 energy for Ne(aug-cc-pVQZ)

-128.86805545187937

I find a referenz energy of -128.8680317 from other sources.

Can you reproduce those energies with your local implementation
of CP2K and with QChem?

Some other points to consider:
- CP2K and QChem have a totally different way of calculating some
of the energies. Except for the total energy it is not clear
which of the intermediate energies to compare. You can establish
this by comparing maybe HF and PBE calculations.
- There is no way the VV10 (nonlocal dispersion) energies will
ever agree. In addition to the different functional implementation
(rVV10 vs. VV10), CP2K only calculates this energy for the 'smooth'
density, which is in turn highly run parameter dependent.
It would be best to compare without this term!

regards

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 06/03/2014 11:25PM
Subject: [CP2K:5365] Re: Implementing new functionals

[August Melcher]

I got -128.868055451878774 using that input file. However I was looking at the xc_b97.F file with Narbe today and he noticed what is likely the cause of the problem we've been noticing. Currently, the exchange energy density is being evaluated without attenuation. The exchange energy density in range-separated functionals uses the erfc(w*r)/r operator instead of 1/r, so there is an additional attenuation function (F) that appears in the integrand. The spin-polarized equation for F can be found in the wB97X-V paper (Equation 9). The only modification that needs to be made to the existing maple code is the inclusion of the F functional between e_lsda_x_a and gx_a (i.e. for the energy evaluation only, it would look like:

exc = scale_x * (e_lsda_x_a * Fx_a * gx_a + e_lsda_x_b * Fx_b *gx_b) + scale_c * (e_lsda_c_ab * gc_ab + e_lsda_c_a * gc_a + e_lsda_c_b * gc_b)

Would you mind modifying the source code or else directing me to the b97 maple file so that I can do so?

[hut...@chem.uzh.ch]

Yes, that is correct. The scaling term is missing.
It seems that Fawzi (he did this functional) has not
submitted the Maple worksheet to the repository.
At least I cannot find it.
I don't know when I will find time to work on the code directly.

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 06/05/2014 12:23AM
Subject: [CP2K:5368] Re: Implementing new functionals

I got -128.868055451878774 using that input file. However I was looking at the xc_b97.F file with Narbe today and he noticed what is likely the cause of the problem we've been noticing. Currently, the exchange energy density is being evaluated without attenuation. The exchange energy density in range-separated functionals uses the erfc(w*r)/r operator instead of 1/r, so there is an additional attenuation function (F) that appears in the integrand. The spin-polarized equation for F can be found in the wB97X-V paper (Equation 9). The only modification that needs to be made to the existing maple code is the inclusion of the F functional between e_lsda_x_a and gx_a (i.e. for the energy evaluation only, it would look like:

exc = scale_x * (e_lsda_x_a * Fx_a * gx_a + e_lsda_x_b * Fx_b *gx_b) + scale_c * (e_lsda_c_ab * gc_ab + e_lsda_c_a * gc_a + e_lsda_c_b * gc_b)

Would you mind modifying the source code or else directing me to the b97 maple file so that I can do so?

[August Melcher]

Ok, I just emailed Fawzi to ask if he has the file or knows where to find it. Just to be sure that I have his current email address, could you verify that this is it? fawzi at gmx dot ch

[August Melcher]

Hello, I was just looking at the GTH_POTENTIALS file and I noticed that Oxygen seems to have an extra parameter compared to the Phys. Rev. B 1996, 54, 1703 paper. The value 0.21720070 does not appear in the paper at all. Could someone tell me what this value corresponds to what I should set it to/how I should choose what to set it to for PBE0?

O GTH-BLYP-q6

    2    4

     0.24342026    2   -16.99189235     2.56614206

    2

     0.22083140    1    18.38885102

     0.21720070    0

Thanks,

August

[hut...@chem.uzh.ch]

Hi

you can ignore that number. The 0 after it, means it it is not used.

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: August Melcher
Sent by: cp...@googlegroups.com

Date: 06/16/2014 10:51PM
Subject: [CP2K:5412] Re: Implementing new functionals