CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
howto:xas_tdp

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
Next revisionBoth sides next revision
howto:xas_tdp [2021/07/29 15:06] – [The GW2X input subsection] abussyhowto:xas_tdp [2021/08/16 14:48] – [Tetrahedral NaAlO$_2$ (K-edge, periodic)] abussy
Line 318: Line 318:
 </code> </code>
  
-There are many performance oriented keywords and subsection in the above input. Most importantly, only one atom is treated at the all-electron level (the one atom from which the excitation takes place), all other are described using pseudopotentials. Also quite important is the usage of the ADMM method. This allows for very efficient evaluation of the HFX energy in the ground state calculation. Finally, the OT iterative solver is used. Since only a handful of eigenvalues are calculated (those within 20.0 eV of the first excitation energy), this scales much better than a full digonalization.+There are many performance oriented keywords and subsection in the above input. Most importantly, only one atom is treated at the all-electron level (the one atom from which the excitation takes place), all other are described using pseudopotentials. Also quite important is the usage of the ADMM method. This allows for very efficient evaluation of the HFX energy in the ground state calculation. Finally, the OT iterative solver is used. Since only a handful of eigenvalues are calculated (those within 20.0 eV of the first excitation energy), this scales much better than a full digonalization. Note that the ''RI_REGION'' keyword is absent (it is set to 0 by default). Since the neighbors of the excited Al atom are described with pseudopotentials, there is no need for extra RI basis function for the projection of the density.
  
 This input file would generate a spectrum such as the one visible on figure 4 of the [[ https://pubs.rsc.org/en/Content/ArticleLanding/2021/CP/D0CP06164F#!divAbstract | reference work]]. This is a much larger calculation than the first example though and would require a few hours on 20-30 processors (mostly to converge the SCF). In you are interested in reproducing this result, input, geometry and pcseg-2/admm-2 basis sets are available {{ :howto:sodal.zip | here}}.  This input file would generate a spectrum such as the one visible on figure 4 of the [[ https://pubs.rsc.org/en/Content/ArticleLanding/2021/CP/D0CP06164F#!divAbstract | reference work]]. This is a much larger calculation than the first example though and would require a few hours on 20-30 processors (mostly to converge the SCF). In you are interested in reproducing this result, input, geometry and pcseg-2/admm-2 basis sets are available {{ :howto:sodal.zip | here}}. 
Line 516: Line 516:
 $$ $$
  
-where $\varepsilon_a$ is the orbital energy of a virtual MO and $\varepsilon_I$ the energy of the donor core MO. Under Koopman's condition, these energies are interpreted as the electron affinity and and the ionization potential (IP). However, DFT is notoriously bad at prediction accurate absolute orbital eigenvalues. Therefore, and because $|\varepsilon_I| >> |\varepsilon_a|$, excitation energies are expected to be widely improved if the DFT energy $\varepsilon_I$ were to be replace by an accurate value of the IP.+where $\varepsilon_a$ is the orbital energy of a virtual MO and $\varepsilon_I$ the energy of the donor core MO. Under Koopman's condition, these energies are interpreted as the electron affinity and and the ionization potential (IP). However, DFT is notoriously bad at predicting accurate absolute orbital eigenvalues. Therefore, and because $|\varepsilon_I| >> |\varepsilon_a|$, excitation energies are expected to be widely improved if the DFT energy $\varepsilon_I$ were to be replace by an accurate value of the IP.
  
 The IP can be accurately calculated using the second-order electron propagator equation: The IP can be accurately calculated using the second-order electron propagator equation:
Line 532: Line 532:
 ==== Simple examples ==== ==== Simple examples ====
  
-=== SO$_2$ molecule (L-edge + SOC) ===+=== OCS molecule (L-edge + SOC) === 
 + 
 +This example covers GW2X corrected L-edge spectroscopy with spin-orbit coupling. 
 + 
 +<code - OCS.inp> 
 + 
 +&GLOBAL 
 +  PROJECT OCS 
 +  PRINT_LEVEL MEDIUM 
 +  RUN_TYPE ENERGY 
 +&END GLOBAL 
 +&FORCE_EVAL 
 +  METHOD Quickstep 
 +  &DFT 
 +    BASIS_SET_FILE_NAME BASIS_GW2X 
 +    POTENTIAL_FILE_NAME POTENTIAL 
 +    AUTO_BASIS RI_XAS MEDIUM 
 + 
 +    &MGRID 
 +      CUTOFF 800 
 +      REL_CUTOFF 50 
 +      NGRIDS 5 
 +    &END MGRID 
 +    &QS 
 +      METHOD GAPW 
 +    &END QS 
 + 
 +    &POISSON 
 +      PERIODIC NONE 
 +      PSOLVER MT  
 +    &END 
 + 
 +    &SCF 
 +      EPS_SCF 1.0E-8 
 +      MAX_SCF 50 
 +    &END SCF 
 + 
 +    &XC 
 +      &XC_FUNCTIONAL                ! The PBEh(45%) functional  
 +         &LIBXC  
 +            FUNCTIONAL GGA_C_PBE  
 +         &END LIBXC  
 +         &LIBXC  
 +            FUNCTIONAL GGA_X_PBE  
 +            SCALE 0.55  
 +         &END LIBXC  
 +      &END XC_FUNCTIONAL 
 + 
 +      &HF 
 +         FRACTION 0.45 
 +      &END HF 
 +    &END XC 
 + 
 +    &XAS_TDP 
 +      &DONOR_STATES 
 +         DEFINE_EXCITED BY_KIND 
 +         KIND_LIST S 
 +         STATE_TYPES 2p          ! Need to look for the S 2p states within the 7 MOs with lowest energy; 
 +         N_SEARCH 7              ! one S 1s, one S 2s, three S 2s , one C 1s and one O 1s 
 +         LOCALIZE                ! Localization is required 
 +      &END DONOR_STATES 
 + 
 +      EXCITATIONS RCS_SINGLET 
 +      EXCITATIONS RCS_TRIPLET 
 +      SOC 
 + 
 +      GRID S  300 500 
 + 
 +      N_EXCITED 150 
 +      TAMM_DANCOFF 
 + 
 +      &GW2X                      ! This is the only difference in the input file with respect to a 
 +      &END GW2X                  ! standard XAS_TDP calculation (defaults parameters are used) 
 + 
 +      &KERNEL 
 +         RI_REGION 3.0 
 +         &XC_FUNCTIONAL  
 +            &LIBXC  
 +               FUNCTIONAL GGA_C_PBE  
 +            &END LIBXC  
 +            &LIBXC  
 +               FUNCTIONAL GGA_X_PBE  
 +               SCALE 0.55  
 +            &END LIBXC  
 +         &END XC_FUNCTIONAL 
 +         &EXACT_EXCHANGE 
 +            FRACTION 0.45 
 +         &END EXACT_EXCHANGE 
 +      &END KERNEL 
 + 
 +    &END XAS_TDP 
 +  &END DFT 
 +  &SUBSYS 
 +    &CELL 
 +      ABC 10.0 10.0 10.0 
 +      PERIODIC NONE 
 +    &END CELL 
 +    &COORD 
 +      C         5.0000000209        4.9999999724        5.2021372095 
 +      O         5.0000000094        5.0000000290        6.3579624316 
 +      S         5.0000000207        5.0000000007        3.6399034216 
 +    &END COORD 
 +    &KIND C 
 +      BASIS_SET aug-pcX-2 
 +      POTENTIAL ALL 
 +    &END KIND 
 +    &KIND O 
 +      BASIS_SET aug-pcX-2 
 +      POTENTIAL ALL 
 +    &END KIND 
 +    &KIND S 
 +      BASIS_SET aug-pcX-2 
 +      POTENTIAL ALL 
 +    &END KIND 
 +  &END SUBSYS 
 +&END FORCE_EVAL 
 + 
 +</code> 
 + 
 +The only difference between the above input file and that of a standard XAS LR-TDDFT calculation is the addition of the ''&GW2X'' subsection. In this case, only default parameters are used, which corresponds to the original GW2X scheme with a convergence threshold of 0.01 eV. Note that the core specific all-electron aug-pcX-2 basis set is used (triple zeta quality). This inputs corresponds to an entry of table II in the [[https://doi.org/10.1063/5.0058124|reference paper]], although slacker parameters are used here (in order to make this tutorial cheap and easy to run, this particular calculations takes ~2 minutes on 4 cores). 
 + 
 +In the output file, the correction for each S $2p$ is displayed. Note that the correction amounts to a shift of 1.9 eV compared to standard XAS LR-TDDFT, leading to a first singlet excitation energy of 164.4 eV (at the L$_3$ edge). This fits [[https://doi.org/10.1016/s0301-0104(97)00111-0|experimental results]] within 0.1 eV. thus clearly improving the XAS LR-TDDFT result. Note that the core IPs, including spin-orbit coupling effects, are also provided. These can be directly used to produce a XPS spectrum. The content of the ''OCS.spectrum'' file yields the corrected spectrum directly. 
 + 
 +<code> 
 +    - GW2X correction for donor MO with spin  1 and MO index    5:                                    
 +                             iteration                convergence (eV)                                
 +                                     1                       10.047536                                
 +                                     2                        1.237503                                
 +                                     3                        0.014416                                
 +                                     4                       -0.000000                                
 +                                                                                                      
 +      Final GW2X shift for this donor MO (eV):   1.927146                                             
 +                                                                                                      
 +                                                                                                      
 +    - GW2X correction for donor MO with spin  1 and MO index    6:                                    
 +                             iteration                convergence (eV)                                
 +                                     1                        6.197650                                
 +                                     2                        4.963008                                
 +                                     3                        0.241838                                
 +                                     4                        0.000439                                
 +                                                                                                      
 +      Final GW2X shift for this donor MO (eV):   1.907648                                             
 +                                                                                                      
 +                                                                                                      
 +    - GW2X correction for donor MO with spin  1 and MO index    7:                                    
 +                             iteration                convergence (eV)                                
 +                                     1                        6.197650                                
 +                                     2                        4.963008                                
 +                                     3                        0.241838                                
 +                                     4                        0.000439                                
 +                                                                                                      
 +      Final GW2X shift for this donor MO (eV):   1.907648                                             
 +                                                                                                      
 +                                                                                                      
 +    Calculations done:                                                                                
 +                                                                                                      
 +    First singlet XAS excitation energy (eV):                165.014087                               
 +    First triplet XAS excitation energy (eV):                164.681850                               
 +    First SOC XAS excitation energy (eV):                    164.396537                               
 +                                                                                                      
 +    Ionization potentials for XPS (GW2X + SOC):              170.602279                               
 +                                                             169.457339                               
 +                                                             169.367465                        
 + 
 +</code>
  
 === Solid NH$_3$ (K-edge, periodic) === === Solid NH$_3$ (K-edge, periodic) ===
 +
 +This is a much larger example of a periodic system, namely solid ammonia. This example is much heavier to run (~45 minutes on 24 cores).
 +
 +<code NH3.inp>
 +
 +&GLOBAL
 +  PROJECT NH3
 +  RUN_TYPE ENERGY
 +  PRINT_LEVEL MEDIUM
 +&END GLOBAL
 +&FORCE_EVAL
 +  METHOD QS
 +  &DFT
 +    BASIS_SET_FILE_NAME BASIS_GW2X
 +    BASIS_SET_FILE_NAME BASIS_ADMM
 +    BASIS_SET_FILE_NAME BASIS_MOLOPT
 +    POTENTIAL_FILE_NAME POTENTIAL
 +    AUTO_BASIS RI_XAS MEDIUM
 +
 +    &QS
 +      METHOD GAPW
 +    &END QS
 +
 +    &MGRID
 +      CUTOFF 600
 +      REL_CUTOFF 50
 +      NGRIDS 5
 +    &END MGRID
 +
 +    &SCF
 +      SCF_GUESS RESTART
 +      EPS_SCF 1.0E-8
 +      MAX_SCF 30
 +
 +      &OT
 +         MINIMIZER CG
 +         PRECONDITIONER FULL_ALL
 +      &END OT
 +
 +      &OUTER_SCF
 +         MAX_SCF 6
 +         EPS_SCF 1.0E-8
 +      &END OUTER_SCF
 +
 +    &END SCF
 +
 +    &AUXILIARY_DENSITY_MATRIX_METHOD
 +      ADMM_PURIFICATION_METHOD NONE
 +    &END AUXILIARY_DENSITY_MATRIX_METHOD
 +
 +    &XC
 +      &XC_FUNCTIONAL
 +         &LIBXC
 +            FUNCTIONAL GGA_X_PBE
 +            SCALE 0.55 
 +         &END
 +         &LIBXC
 +            FUNCTIONAL GGA_C_PBE
 +         &END
 +      &END XC_FUNCTIONAL
 +      &HF
 +         FRACTION 0.45
 +         &INTERACTION_POTENTIAL
 +            POTENTIAL_TYPE TRUNCATED
 +            CUTOFF_RADIUS 5.0
 +         &END INTERACTION_POTENTIAL
 +      &END HF
 +    &END XC
 +
 +    &XAS_TDP
 +      &DONOR_STATES
 +         DEFINE_EXCITED BY_KIND
 +         KIND_LIST Nx
 +         STATE_TYPES 1s
 +         N_SEARCH 1
 +         LOCALIZE
 +      &END DONOR_STATES
 +
 +      TAMM_DANCOFF
 +      GRID Nx 300 500
 +      E_RANGE 30.0
 +
 +      &GW2X
 +      &END
 +
 +      &KERNEL
 +         &XC_FUNCTIONAL
 +            &LIBXC
 +               FUNCTIONAL GGA_X_PBE
 +               SCALE 0.55 
 +            &END
 +            &LIBXC
 +               FUNCTIONAL GGA_C_PBE
 +            &END
 +         &END XC_FUNCTIONAL
 +         &EXACT_EXCHANGE
 +            OPERATOR TRUNCATED
 +            CUTOFF_RADIUS 5.0
 +            FRACTION 0.45
 +         &END EXACT_EXCHANGE
 +      &END KERNEL
 +
 +    &END XAS_TDP
 +  &END DFT
 +  &SUBSYS
 +    &CELL
 +      ABC   10.016118  10.016118  10.016118
 +    &END CELL
 +    &TOPOLOGY
 +      COORD_FILE_FORMAT XYZ
 +      COORD_FILE_NAME NH3.xyz
 +    &END TOPOLOGY
 +    &KIND H
 +      BASIS_SET DZVP-MOLOPT-SR-GTH
 +      BASIS_SET AUX_FIT FIT3
 +      POTENTIAL GTH-PBE
 +    &END KIND
 +    &KIND N
 +      BASIS_SET DZVP-MOLOPT-SR-GTH
 +      BASIS_SET AUX_FIT FIT3
 +      POTENTIAL GTH-PBE
 +    &END KIND
 +    &KIND Nx
 +      ELEMENT N
 +      BASIS_SET aug-pcseg-2
 +      BASIS_SET AUX_FIT aug-admm-2
 +      POTENTIAL ALL
 +    &END KIND
 +  &END SUBSYS
 +&END FORCE_EVAL
 +
 +</code>
 +
 +Again, the only difference with respect to a standard XAS-LRTDDFT input file is the ''&GW2X'' subsection. This input file corresponds exactly to figure 3 a) of the [[ https://doi.org/10.1063/5.0058124 |reference papaer]]. In this case, the GW2X correction amounts to a blue shift of 3.7 eV, aligning the calculated spectrum to the experimental one remarkably well. The NH3 structure file as well as the necessary basis set file (also for the OCS example) are available {{ :howto:gw2x.zip |here}}.
 +
  
 ==== FAQ ==== ==== FAQ ====
  
 === How can I make the GW2X correction run faster ?=== === How can I make the GW2X correction run faster ?===
 +
 +The GW2X correction scheme scales cubically with the number of MOs in the system. Therefore, the best way to improve performance is to reduce that number. Because an accurate description of the core region is only necessary for the exited atoms, all other atoms can be described with pseudopotentials. This drastically reduces the number of MOs since only valence states are kept. In the solid NH3 example above, all nitrogen atoms are equivalent under symmetry. Therefore, their individual contribution to the XAS spectrum is bound to be the same. This allows for the description of a single nitrogen atom at the all-electron level, while all others (and the hydrogens) use pseudopotentials. Note that the ADMM approximation is also utilized. This greatly reduces the cost of the underlying hybrid DFT calculation, as well as the evaluation of the generalized Fock matrix as required by GW2X.
  
 === Why don't I get the absolute core IP in periodic systems ? === === Why don't I get the absolute core IP in periodic systems ? ===
 +
 +For molecules in non-periodic boundary conditions, the potential is such that it is zero far away. In the periodic case, the zero is ill defined. As a consequence, all Kohn-Sham eigenvalues end up shifted by some unknown, constant amount. Therefore, their absolute values and that of the calculated IP cannot be interpreted in a physical manner. However, the correction scheme depends on the difference $|\varepsilon_a-\varepsilon_I|$, where the shift cancels out. 
 +
 +=== Why is the LOCALIZE keyword required ? ===
 +In order to efficiently evaluate the antisymmetric integrals of the type $\langle Ia || jk \rangle$, the same local RI scheme as XAS_TDP is used. Therefore, the core state $I$ needs to be local in space. However, the rotation required to get the pseudocanonical orbitals needed for the original GW2X scheme may break this localization, provided that there are other equivalent atoms in the system. To prevent that from happening, all core states localized on other atoms are ignored for the rotation and the subsequent IP calculation. This has negligible impact since core states belonging to different atoms only weakly interact. It is however important to keep the value of the LOCALIZE keyword to a minimum to insure that only core states are ignored.
howto/xas_tdp.txt · Last modified: 2024/02/24 10:01 by oschuett