TDDFPT calculations with CP2K7.1 using the range-separated functionals

[liu xiangyang]

Dear All,

I have tried to use the TDDFPT method implemented in CP2K7.1 to do excited state calculations with the range-separated functionals such as wB97XD. 

However, after several tests with a small molecule, namely H2Pc, I found that the first two excitation energies are greatly underestimated in comparison with the LR-TDDFT results obtained in  GAUSSIAN09 (ca. 0.51 and 0.60 eV (CP2K7.1) vs. 2.18 and 2.26 eV (GAUSSIAN09)).

I wonder whether there are some mistakes with my input file or there are some problem of TDDFPT for such calculations?

The input file used in my calculations is written as follows:

&GLOBAL

  PROJECT tddfpt

  RUN_TYPE energy

  PRINT_LEVEL medium

&END GLOBAL

&FORCE_EVAL

  METHOD Quickstep

  &PROPERTIES

    &TDDFPT

       NSTATES  5            # number of excited states

    &END TDDFPT

  &END PROPERTIES

&DFT

    BASIS_SET_FILE_NAME GTH_BASIS_SETS

    POTENTIAL_FILE_NAME POTENTIAL

    BASIS_SET_FILE_NAME BASIS_ADMM_MOLOPT

    BASIS_SET_FILE_NAME BASIS_ADMM

    &AUXILIARY_DENSITY_MATRIX_METHOD

      METHOD BASIS_PROJECTION

      ADMM_PURIFICATION_METHOD NONE 

      EXCH_CORRECTION_FUNC BECKE88X

    &END AUXILIARY_DENSITY_MATRIX_METHOD

    CHARGE 0 

    &MGRID

      CUTOFF 400

    &END MGRID

    &QS

      METHOD gpw

      EPS_PGF_ORB 1e-12

    &END QS

    &SCF

      MAX_SCF   100

      EPS_SCF   1e-5

      SCF_GUESS atomic

      &DIAGONALIZATION

        ALGORITHM STANDARD

      &END DIAGONALIZATION

      &MIXING T

        ALPHA 0.5

        METHOD PULAY_MIXING

        NPULAY 5

      &END MIXING

    &END SCF

    &POISSON

       PERIODIC XYZ

    &END POISSON

    &XC

      &XC_FUNCTIONAL

        &LIBXC

          FUNCTIONAL HYB_GGA_XC_WB97X_D

        &END

      &END XC_FUNCTIONAL

      &HF

       &SCREENING

         EPS_SCHWARZ 1.0E-6

       &END

       &MEMORY

         MAX_MEMORY 100

       &END

       &INTERACTION_POTENTIAL

         POTENTIAL_TYPE MIX_CL_TRUNC

         OMEGA 0.2

         SCALE_LONGRANGE 0.777964

         SCALE_COULOMB 0.222036

         CUTOFF_RADIUS 5.0

         T_C_G_DATA t_c_g.dat

       &END

      &END

      &XC_GRID

        XC_DERIV SPLINE2_SMOOTH      # this is needed for the 2nd derivatives of the XC functional

      &END XC_GRID

    &END XC

  &END DFT

  &SUBSYS

    &TOPOLOGY

       &CENTER_COORDINATES

       &END CENTER_COORDINATES

    &END TOPOLOGY

    &CELL

      ABC 30.0 30.0 30.0

      PERIODIC XYZ

    &END CELL

    &COORD

 N                  0.01468800    2.00544200   -0.29683900

 N                 -2.37562900    2.41155000   -0.13186300

 N                 -2.41200400   -2.37377100   -0.17491600

 N                  2.41112400    2.37299900   -0.12658700

 N                 -0.01307800   -2.00591400   -0.34375500

 N                  2.37788800   -2.41418900   -0.17065800

 C                 -4.17222800    0.73645800   -0.00309900

 C                 -0.67412900    4.18239500    0.01934000

 C                 -4.18380100   -0.67562700   -0.02014000

 C                  0.73808000    4.17070300    0.02178600

 C                 -2.78102300    1.14403300   -0.18692600

 C                 -1.10192300    2.79691600   -0.16605600

 C                 -2.79906300   -1.10173400   -0.21220500

 C                  1.14375300    2.77826000   -0.16260500

 C                 -0.73890800   -4.16630900   -0.00019300

 C                  4.18311400    0.67096100    0.00115100

 C                  0.67355900   -4.17905700    0.00226800

 C                  4.17314700   -0.74140800   -0.01277100

 C                 -1.14329800   -2.77792600   -0.20621200

 C                  2.79688300    1.09912600   -0.17512500

 C                  1.10326800   -2.79777700   -0.20275300

 C                  2.78178900   -1.14729500   -0.19912900

 C                 -5.34524900    1.46617200    0.20504700

 C                 -1.38645100    5.36877700    0.20954200

 C                 -5.36863500   -1.39063500    0.17089500

 C                  1.47037200    5.34443400    0.21437700

 C                 -6.53049700    0.75063900    0.38667600

 C                 -0.65342400    6.54324100    0.39406800

 C                 -6.54209900   -0.66018000    0.36971000

 C                  0.75719100    6.53149500    0.39645600

 C                 -1.47251300   -5.33425600    0.22345700

 C                  5.36943600    1.38205000    0.19903500

 C                  1.38473700   -5.36039200    0.22884400

 C                  5.34798200   -1.47459500    0.17084800

 C                 -0.76071300   -6.51537100    0.44153000

 C                  6.54461400    0.64821000    0.37381100

 C                  0.65044900   -6.52827600    0.44412700

 C                  6.53422100   -0.76263800    0.36007100

 H                 -5.31982100    2.55648400    0.22741200

 H                 -2.47726000    5.36199000    0.21499900

 H                 -5.36047700   -2.48141400    0.16783300

 H                  2.56090800    5.31878300    0.22274500

 H                 -7.46690300    1.28881400    0.54987000

 H                 -1.17840000    7.48942700    0.54335800

 H                 -7.48727000   -1.18662900    0.52012100

 H                  1.29740200    7.46875700    0.54746900

 H                 -2.56294400   -5.30668600    0.23110100

 H                  5.36158400    2.47267500    0.21857200

 H                  2.47541900   -5.35281000    0.24119000

 H                  5.32277600   -2.56515200    0.16738300

 H                 -1.30192600   -7.44735900    0.61932200

 H                  7.49034300    1.17231500    0.52877000

 H                  1.17373400   -7.47001200    0.62407400

 H                  7.47217800   -1.30348900    0.50426500

 H                 -0.95410508    0.01210298   -0.52371821

 H                  0.94993522   -0.00596112   -0.48792422

 N                 -2.00779500    0.01729500   -0.33767900

 N                  2.00585900   -0.01636500   -0.31524900

    &END COORD

    &KIND N

      BASIS_SET DZVP-GTH

      POTENTIAL GTH-PBE

      BASIS_SET AUX_FIT cFIT3

    &END KIND

    &KIND C

      BASIS_SET DZVP-GTH

      POTENTIAL GTH-PBE

      BASIS_SET AUX_FIT cFIT3

    &END KIND

    &KIND H

      BASIS_SET DZVP-GTH

      POTENTIAL GTH-PBE

      BASIS_SET AUX_FIT cFIT3

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

Best wishes, 

Xiang-Yang Liu

[Vladimir Rybkin]

Dear Xiang-Yang Liu,

most importantly: you are using different basis sets in Gaussian and CP2K (CP2K also use pseudopotentials). With this difference in mind you differences are within reasonable. Generally, difference below 0.1 eV for TDDFT implementations is not "great".

Yours,

Vladimir

понедельник, 11 января 2021 г. в 14:49:47 UTC+1, lxyl...@gmail.com:

[Lucas Lodeiro]

Hi Liu,

I did not run TDDFT calculations, but I did some tests between CP2K and other programs like G09. As Vladimir says, your basis sets are not the same, and some difference could appear due to this reason. But in your case the differences are big, 1.6 eV approx. I found that some default settings of convergences criterium are differents, for example the EPS_SCF which is 1E-8 in G09, you could tight your convergence criterion, EPS_SCF, EPS_DEFAULT, EPS_SCHWARZ to -8, -12 and -8 to get results with similar convergences in both programs. Also, you are using a PBC calculation in a big cell, but maybe it is no sufficient to mimic the isolated molecule as in G09... and yout cutoff radius for HF is a little bit short, if you run a non-periodic calculation, you can use just the long range potential without the truncation.

Finally, just to speed up, you can use OT instead of diagonalization method, with it you can use ADMM_PURIFICATION_METHOD MO_DIAG.

In order to have the same basis sets, as vladimir says, you could explore  to use the same basis sets in both programs, you can get basis sets to both programs from: https://www.basissetexchange.org/

And the Auxiliary basis set is the minimum one, you could explore FIT and pFIT basis set to check if the result is sensitive.

Regards

Lucas Lodeiro

--
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 view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/4897d4d1-a074-44f6-89ff-3b2c613eab83n%40googlegroups.com.

[Vladimir Rybkin]

Sorry, I really read the differences incorrectly. Thank you, Lukas, for correcting me. 

понедельник, 11 января 2021 г. в 21:02:30 UTC+1, Lucas Lodeiro:

[liu xiangyang]

Hi  Vladimir  & Lucas，

Thanks a lot for your kindly responses. Actually, I have also tested small molecule such as H2CO using present settings and the results obtained by CP2K and G09 are quanlitatively in good agreement with each other and THAT IS THE FACT THAT PUZZLES ME MOST.

According to your suggestions, I have made several changes to my input file, including removal of the periodic conditions, removal of the ADMM approximations to avoid possible mistakes, and  lower the convergence in G09 since the computational efforts of CP2K is much larger. The basis sets are not changed since the TDDFPT in CP2K is implemented only in GPW method, which is unable to do calculations with all electron basis sets such as Pople basis sets. However, I believe the difference of basis set (DZVP-GTH vs. 6-31G*) should not be so large. Unfortunately, the results of CP2K7.1 is still about 1 eV lower than that obtained in G09 (1.19 and 1.25 eV in CP2K vs. 2.18 and 2.26 eV in G09). 

PS. I have also run GAPW calculations in CP2K using the WB97XD/6-31G** and the resulted ground state energy are nearly the same as that obtained in G09 (differenceless than 0.002 hartree). 

The modified DFT parts of CP2K  input file is attatched also:

&DFT

    BASIS_SET_FILE_NAME BASIS_MOLOPT

    POTENTIAL_FILE_NAME POTENTIAL

    CHARGE 0 

    &MGRID

      CUTOFF 400

    &END MGRID

    &QS

      METHOD gpw

      EPS_PGF_ORB 1e-12

    &END QS

    &SCF

      MAX_SCF   100

      EPS_SCF   1e-5

      SCF_GUESS atomic

      &DIAGONALIZATION

        ALGORITHM STANDARD

      &END DIAGONALIZATION

      &MIXING T

        ALPHA 0.5

        METHOD PULAY_MIXING

        NPULAY 5

      &END MIXING

    &END SCF

    &POISSON

       PERIODIC NONE

       PSOLVER MT

    &END POISSON

    &XC

      &XC_FUNCTIONAL

        &LIBXC

          FUNCTIONAL HYB_GGA_XC_WB97X_D

        &END

      &END XC_FUNCTIONAL

      &HF

       &SCREENING

         EPS_SCHWARZ 1.0E-6

       &END

       &MEMORY

         MAX_MEMORY 100

       &END

       &INTERACTION_POTENTIAL

         POTENTIAL_TYPE MIX_CL

         OMEGA 0.2

         SCALE_LONGRANGE 0.777964

         SCALE_COULOMB 0.222036

       &END

      &END

    &END XC

  &END DFT

Best wishes,

Xiang-Yang Liu

[Lucas Lodeiro]

Hi Liu,

As you mention it is a little weird if other molecules work fine with the same input. I guess the next test is to tighten the convergence criteria in CP2K if it is possible.

Another thing, remember that for non periodic calculations, you have to set periodic none in the poisson and cell section.

    &CELL

      ABC 30.0 30.0 30.0

      PERIODIC NONE

    &END CELL

I see the transition energies change a little bit, it is at least comfortable, showing that one or more changes you did affect the result. Maybe as the last chance there will be  useful to know which change was the one that caused the variation a the energies.

Sorry, it does not seem a simple problem.

Regards

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/64940b8a-9bef-4f42-9833-55b0ddaa7326n%40googlegroups.com.

[Fangyong Yan]

Hi, XIangyan,

I used you structure, and tested both Gaussian and cp2k, and they have good agreement. 

So Gaussian and cp2k can both give reasonable results for the excitation energies of your molecule. 

Gaussian 09, blyp/6-31G*, 

Excited State   1:      Singlet-A      1.9414 eV  638.63 nm  f=0.2542  <S**2>=0.000

Excited State   2:      Singlet-A      1.9724 eV  628.58 nm  f=0.3319  <S**2>=0.000

Excited State   3:      Singlet-A      2.5113 eV  493.71 nm  f=0.0038  <S**2>=0.000

Excited State   4:      Singlet-A      2.5205 eV  491.90 nm  f=0.0341  <S**2>=0.000

Excited State   5:      Singlet-A      2.5764 eV  481.23 nm  f=0.0001  <S**2>=0.000

cp2k, blyp/TZVP, with GTH pseudopotential, 

I get

         State    Excitation        Transition dipole (a.u.)        Oscillator

         number   energy (eV)       x           y           z     strength (a.u.)

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

 TDDFPT|      1       2.08760   2.4868E+00 -8.6129E-03  6.4653E-03   3.16297E-01

 TDDFPT|      2       2.16739  -1.4224E-02 -2.9924E+00 -1.0393E-02   4.75501E-01

 TDDFPT|      3       2.52956  -1.3916E-01 -1.4062E-03  1.0520E-04   1.20031E-03

 TDDFPT|      4       2.56128   8.4944E-01 -4.4270E-03  1.5135E-04   4.52790E-02

 TDDFPT|      5       2.59590  -1.3595E-04  4.4904E-02 -4.4546E-02   2.54439E-04

 TDDFPT|      6       2.60688  -3.1839E-03  4.1336E-02  2.5674E-02   1.51876E-04

 TDDFPT|      7       2.76485  -3.0533E-03 -2.4437E-01  1.0913E-03   4.04581E-03

 TDDFPT|      8       2.80552  -2.7107E-01 -4.8867E-03 -9.2586E-04   5.05222E-03

 TDDFPT|      9       2.88205  -1.2938E-02 -1.3849E-04  1.2856E-03   1.19366E-05

 TDDFPT|     10       2.92954   2.8231E-05  6.4011E-03 -8.3886E-02   5.07991E-04

****************** Below is the input. 

  &SUBSYS

    &CELL

      ABC  30. 30. 6.

      PERIODIC none

    &END CELL

    &COORD

    &KIND C

      BASIS_SET TZVP-GTH

      POTENTIAL GTH-PBE-q4

    &END KIND

    &KIND N

      BASIS_SET TZVP-GTH

      POTENTIAL GTH-PBE-q5

    &END KIND

    &KIND H

      BASIS_SET TZVP-GTH

      POTENTIAL GTH-PBE-q1

    &END KIND

    &QS

      METHOD GPW

      MAP_CONSISTENT  T

      EPS_DEFAULT 1.0E-9

    &END QS

    &XC

      &XC_FUNCTIONAL BLYP

      &END XC_FUNCTIONAL

      &XC_GRID

        XC_DERIV SPLINE2_SMOOTH      # this is needed for the 2nd derivatives of the XC functional

      &END XC_GRID

    &END XC

    &SCF

      MAX_SCF   250

      EPS_SCF   1e-8

      SCF_GUESS atomic

    &END SCF

    &MGRID

      CUTOFF 600

      NGRIDS 4

      REL_CUTOFF 60

    &END MGRID

    &POISSON

       PERIODIC none

       POISSON_SOLVER MT

    &END POISSON

  &PROPERTIES

    &TDDFPT

       NSTATES     10             # number of excited states

       MAX_ITER    100             # maximum number of Davidson iterations

       CONVERGENCE [eV] 1.0e-3    # convergence on maximum energy change between iterations

       &MGRID

          CUTOFF 300  # separate cutoff for TDDFPT calc

       &END

    &END TDDFPT

  &END PROPERTIES

Regares,

Fangyong

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAOFT4PJ1bAsK-iK5Lnd7YikGQjJYgDxb47dNh5bNH_FeM1U%2BrQ%40mail.gmail.com.

[Fangyong Yan]

Hi, XIangyan,

I changed the functional type for atoms, from GTH-PBE-q4 => GTH-BLYP-q4, since I am using BLYP. (Sorry, I just copied your basis set and functionals and did not realize the difference.)

Here is the result, and is in good agreement with Gaussian 09.

         State    Excitation        Transition dipole (a.u.)        Oscillator

         number   energy (eV)       x           y           z     strength (a.u.)

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

 TDDFPT|      1       2.08325  -2.4674E+00  8.6733E-03 -6.5509E-03   3.10723E-01

 TDDFPT|      2       2.16407   1.4338E-02  2.9818E+00  1.0596E-02   4.71397E-01

 TDDFPT|      3       2.52106   1.3816E-01  1.4082E-03 -1.0486E-04   1.17904E-03

 TDDFPT|      4       2.55291  -8.4802E-01  4.4219E-03 -1.5266E-04   4.49799E-02

 TDDFPT|      5       2.58935   1.4471E-04 -4.8114E-02  4.4750E-02   2.73894E-04

I will test your functional, since it is hybrid, it will take some time. 

Regards,

Fangyong

[Fangyong Yan]

Hi, Xiangyang,

I tested your functional, and got the similar result as yours. 

However, I have tested both blyp and pbe, both are in good agreement with Gaussian. So the discrepancy may only happen for hybrid functionals, and I dont know the reason.

Regards,

Fangyong

[Fangyong Yan]

Hi, XIangyang,

One possibility is that your structure has been optimized using Gaussian, which uses all-electron calculation (I assume your structure is optimized). While cp2k uses pseudopotential. So the Gaussian optimized structure may not be the cp2k optimized structure. 

While this optimization issue may not be a problem with PBE/BLYP, it may be a problem for hybrid functional, especially for your structure. 

So maybe you can optimize your structure using pseudopotential in cp2k, once you get your optimized structure, you do tddfpt calculation. 

Regards,

Fangyong

[liu xiangyang]

Hi, Fangyong,

Thanks a lot for your kindness and patience for doing so much testings and Sorry for the delay to reply so late since I have not login this site these days.

I agree with your results since I have also test the results using pure functionals GGA such as PBE functional previously and the results agree perferctly well when using G09 and CP2K respectively, just as the BLYP functional you tested in the first place. However, as you have found, the results obtained using the range separated functional wB97XD are distinct from each other. 

I have noticed your suggestion that the optimization might be the origin of the discrepancy. However, this structure is modified from a previously optimized structure of ZnPc and I have not optimized it at all.

Since TDDFT is a single point calculation and I believe the single point calculations using different programs should obtain similar results, which is the spirit of Science.

I tend to believe that there are some differences between the TDDFPT method implemented in CP2K and TDDFT in G09, which makes the TDDFPT unsuitable for the calculation using the range-separated functionals.

However, I have not found a efficient that could fix such differences yet and I believe this issue is important since the GGA functional might not be accurate enough for describing charge transfer states, which always need range separated functionals.

Sincerely Yours,

Xiang-Yang Liu

[Fangyong Yan]

Hi, Xiang-Yang,

Regards your comment:

"Since TDDFT is a single point calculation and I believe the single point calculations using different programs should obtain similar results, which is the spirit of Science."

I partially agree with you:

1) The cp2k uses pseudopotential for TDDFT, and Gaussian uses all-electron, so they should not give the exactly same result;

2) the main difference only happens for hybrid functional, in Gaussian the potential operator is 1/r; in cp2k, you use mix_CL_truncated, so it can result in difference. Also, I dont know how the parameters are set up, such as omega, scaled coulomb, scaled long range, etc. 

Regards,

Fangyong

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/6540e59d-ff9f-4c8c-b550-2247e0d60cd1n%40googlegroups.com.

[Fangyong Yan]

Hi, Xiang-Yang,

I tested B3LYP some time ago using your geometry, which is reasonable:

         State    Excitation        Transition dipole (a.u.)        Oscillator

         number   energy (eV)       x           y           z     strength (a.u.)

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

 TDDFPT|      1       1.86120  -2.0023E+00  8.1993E-03 -3.3226E-03   1.82825E-01

 TDDFPT|      2       1.92870   9.9269E-03  2.2110E+00  4.4788E-03   2.31001E-01

 TDDFPT|      3       3.01758   8.6432E-02  1.8994E-04  7.7692E-06   5.52295E-04

 TDDFPT|      4       3.04539  -6.9340E-01  4.6654E-03 -6.3691E-05   3.58749E-02

 TDDFPT|      5       3.25891  -1.5458E-03 -7.2118E-02 -3.4034E-02   5.07928E-04

Below is the input (part):

      ABC  30. 30. 16.

      PERIODIC none

    &END CELL

    &COORD

 N                  0.01468800    2.00544200   -0.29683900

 N                 -2.37562900    2.41155000   -0.13186300

....

        &INTERACTION_POTENTIAL

          ! for condensed phase systems

          POTENTIAL_TYPE TRUNCATED

          ! should be less than halve the cell

          CUTOFF_RADIUS 6.0

          ! data file needed with the truncated operator

          T_C_G_DATA ./t_c_g.dat

        &END

Regards,

Fangyong

[Vladimir Rybkin]

Actually, the CUTOFF_RADIUS should be a bit less than the half of the smallest cell dimension, i.e. <16/2. I would take it as 7.9. Should improve the result.

среда, 20 января 2021 г. в 19:53:51 UTC+1, fyya...@gmail.com: