Open SourceMolecular Dynamics

Sidebar

For Developers

exercises:2017_uzh_cp2k-tutorial:wfc

This is an old revision of the document!

In order to go beyond GGA and hybrid DFT, one option is to use wave function correlation methods. Recently, second-order Møller-Plesset perturbation theory (MP2) and random phase approximation (RPA) have been added to CP2K . The implementations are aimed at condensed phase calculations, and in the case of MP2 provides energies, forces and stress, as well as electron (spin) densities at O(N^5) cost, while RPA provides energies at O(N^4) or O(N^3) cost.

However, significant computational resources are needed for the condensed phase calculations (for the overview see the corresponding lecture). Therefore, we will perform the gas phase calculations in this tutorial, even though RI-GPW is not very efficient in this case.

Some references:

MP2 and RPA energies implementations in CP2K:10.1016/j.cpc.2014.10.021, 10.1021/ct4002202, 10.1021/ct300531w

MP2 forces: 10.1063/1.4919238 and 10.1021/acs.jctc.6b00015

Cubic scaling RPA implementation: 10.1021/acs.jctc.6b00840

Since the correlated wave function calculations are expensive, please use 4 cores for execution (with one OMP thread):

mpirun -np 4 -x OMP_NUM_THREADS=1  cp2k.psmp -i input -o output

1. Task: Benzene dimer MP2 binding energy

Employ the provided input file to compute the benzene dimer binding energy. The provided dimer geometry is optimized already. To obtain the energy of the monomer, geometry optimization is in principle necessary, although the geometry of a benzene monomer must be close to the one in the dimer. This can be done by removing one of the benzene molecules editing the coordinates in the input. Perform a few optimization steps to see how much the energy changes and to calculate electron density (saved as .cube files) at the MP2 level. If energy does not change much, this means that the monomer structure is already close to the optimal. Compare the times needed for one geometry optimization step (energy+forces) for a monomer and energy for a dimer, as well as the relative timing for energy and forces evaluation for a monomer.

During the optimization of benzene, one will calculate gradient which, in turn, requires density matrices. Hence, one can calculate electronic densities. Add the following to the &GLOBAL section:

&PRINT        MEDIUM
&END

and the following lines to the &DFT section

&PRINT
&E_DENSITY_CUBE MEDIUM
&END
&END

Importantly, during the force calculations one will have to solve the coupled-perturbed equations invoking exact exchange calculations. If there is enough memory we can reuse the integrals from the HF calculation by setting the following keyword in the &RI_MP2 section:

FREE_HFX_BUFFER .FALSE.

Perform two optimizations setting FREE_HFX_BUFFER to .TRUE. and .FALSE.. Compare the overall timings and especially the times for performing Hartree-Fock exchange calculations:

integrate_four_center               71 12.3    2.261    4.012  108.805  109.179

2. Task: Benzene monomer RPA energy: frequency integration

RPA is HFX+RPA correlation. It can be performed with HFX orbitals and eigenvalues, but also based on e.g. GGA or hybrid orbitals. Two advantages over MP2 are : 1) scaling is O(N^4) (for RI-dRPA with frequency integration). 2) it is applicable for systems with a small gap. A O(N^3)implementation is also available in CP2K, although with a very large prefactor, which makes it worth using for very large systems.

Here, we look at the convergence of the RPA energy as a function of the number of integration points (RPA_NUM_QUAD_POINTS). Change this parameter in the range 4-16. For large systems the time needed for the calculations is proportional to the number of integration points.

Use the following section to change to RPA:

      ! with WF correlation
&WF_CORRELATION
! use the RI-GPW approach
METHOD  RI_RPA_GPW
&WFC_GPW
&END
&RI_RPA
! number of quadrature points, essential for accurate energies.
! small gap systems need more points
! essentially always use minimax
MINIMAX
&HF
FRACTION 1.0000000
&SCREENING
EPS_SCHWARZ 1.0E-8
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&END HF
&END RI_RPA
MEMORY    1800
NUMBER_PROC  1
&END

It is strongly recommended to use the MINIMAX integration scheme.

Required files

Input files and a wave function restart file can be found in the archive: mp2_rpa.tar.gz input file for an RI-MP2 calculation on a benzene dimer

benzene_dimer.inp
&GLOBAL
PROJECT      benzene_dimer
RUN_TYPE     ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD         Quickstep
&DFT
! specification of basis and potential files (cp2k/data)
BASIS_SET_FILE_NAME    HFX_BASIS
POTENTIAL_FILE_NAME    HF_POTENTIALS
! Wave function restart file to accelerate HF convergence
WFN_RESTART_FILE_NAME  benzene_dimer.wfn
&MGRID
CUTOFF    400
&END MGRID

&QS
METHOD GPW
EPS_DEFAULT  1.0E-10
EPS_PGF_ORB  1.0E-8
&END QS

! standard OT
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-6
MAX_SCF 40
&OT
MINIMIZER       CG
PRECONDITIONER  FULL_SINGLE_INVERSE
&END
&OUTER_SCF
EPS_SCF 1.0E-6
MAX_SCF 10
&END
&END SCF

! Non periodic calculation needs Poisson solver: use wavelet solver
&POISSON
PERIODIC NONE
PSOLVER  WAVELET
&END POISSON

&XC
! no XC functional
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
! and 100%HFX
&HF
FRACTION    1.0
&SCREENING
EPS_SCHWARZ          1.0E-9
&END SCREENING
&MEMORY
MAX_MEMORY  1800
&END
&END HF
! with WF correlation
&WF_CORRELATION
! use the RI-GPW approach
METHOD  RI_MP2_GPW
&WFC_GPW
CUTOFF      200
REL_CUTOFF  50
&END
! block sizes and memory affect parallel efficiency
&RI_MP2
BLOCK_SIZE       1
&END
MEMORY    1800
NUMBER_PROC  1
&END
&END XC
&END DFT

&SUBSYS
! sufficiently large non-periodic unit cell
&CELL
ABC [angstrom]  12 12 12
PERIODIC NONE           ! Non periodic calculation.
&END CELL

! specification of the system with coordinates
&COORD
H         8.5709951714        6.1617188657        9.1769626364
C         8.0387778483        5.9757981379        8.2528397815
C         6.6427474511        5.9561428224        8.2430496957
H         6.0937543445        6.1249650636        9.1610309894
C         5.9546039990        5.7314555716        7.0475649921
H         4.8717398728        5.7196290593        7.0404902271
C         6.6635660065        5.5279574532        5.8611612413
H         6.1300650069        5.3536860275        4.9350551393
C         8.0615723757        5.5446390494        5.8709776417
H         8.6121087088        5.3836499080        4.9523278088
C         8.7479484762        5.7683467907        7.0666774749
H         9.8306982448        5.7893185194        7.0736980974
H         6.4284687494        8.8405192314        5.8228620343
H         8.9057373702        8.8752743844        5.8387264727
C         6.9607954533        9.0255100555        6.7471184900
C         8.3568440363        9.0439429763        6.7567676586
C         6.2520232187        9.2329626290        7.9335080936
C         9.0454365692        9.2673522387        7.9522880615
H         5.1692478273        9.2131682112        7.9263504668
H        10.1283042427        9.2780238064        7.9594286415
C         6.9387312641        9.4553058890        9.1292884396
C         8.3367433685        9.4704785560        9.1388664843
H         6.3883177411        9.6163935534       10.0479844468
H         8.8703567259        9.6435467092       10.0651664724
&END

! keep atoms away from box borders,
! a requirement for the wavelet Poisson solver
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY

! MP2 needs correlation consistent basis set
! RI-MP2 needs an auxiliary basis set
! We employ HF pseudo potentials
&KIND H
BASIS_SET         cc-TZV2P-GTH
RI_AUX_BASIS_SET  TZ_fiodena_opt
POTENTIAL         GTH-HF-q1
&END KIND
&KIND C
BASIS_SET         cc-TZV2P-GTH
RI_AUX_BASIS_SET  TZ_fiodena_opt
POTENTIAL         GTH-HF-q4
&END KIND
&END SUBSYS
&END FORCE_EVAL

! how to propagate the system, selection via RUN_TYPE in the &GLOBAL section
&MOTION
&GEO_OPT
OPTIMIZER BFGS ! Good choice for 'small' systems (use LBFGS for large systems)
MAX_ITER  100
MAX_DR    [bohr] 0.003 ! adjust target as needed
&BFGS
&END
&END
&END