Skip to content
Permalink
Browse files

gCP potentials and Short range bond correction for vdW D3 method (ge…

…neral formula) (#417)

* Short range bond correction for vdW D3 method, parameterize general formula

* gCP potentials input

* gCP energy/forces

* Debug forces

* gCP method

* Update to latest trunk version

* Adjust input interface

* Copyright banner and remove replicated file
  • Loading branch information...
juerghutter committed Jun 17, 2019
1 parent d22e34f commit ec01bc2fbd6ca57f2730fe5a4690e2cc169df67b

Large diffs are not rendered by default.

@@ -0,0 +1,62 @@
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-12
&END
&END
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3(BJ)
D3BJ_SCALING 1.0000 0.4171 0.8777 2.9149
PARAMETER_FILE_NAME dftd3.dat
SHORT_RANGE_CORRECTION TRUE
SHORT_RANGE_CORRECTION_PARAMETERS 0.03 0.70 1.50 0.75
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
&gcp_potential
GLOBAL_PARAMETERS 0.1290 1.1549 1.1763 1.1526

DELTA_ENERGY H 0.04240
DELTA_ENERGY He 0.02832

DELTA_ENERGY Li 0.17787
DELTA_ENERGY Be 0.17160
DELTA_ENERGY B 0.22424
DELTA_ENERGY C 0.27995
DELTA_ENERGY N 0.35791
DELTA_ENERGY O 0.47901
DELTA_ENERGY F 0.63852
DELTA_ENERGY Ne 0.83235

DELTA_ENERGY Na 1.11411
DELTA_ENERGY Mg 1.27115
DELTA_ENERGY Al 1.44695
DELTA_ENERGY Si 1.61098
DELTA_ENERGY P 1.76661
DELTA_ENERGY S 1.98823
DELTA_ENERGY Cl 2.22845
DELTA_ENERGY Ar 2.48796

DELTA_ENERGY K 0.37425
DELTA_ENERGY Ca 0.46097
DELTA_ENERGY Sc 0.44489
DELTA_ENERGY Ti 0.40499
DELTA_ENERGY V 0.37841
DELTA_ENERGY Cr 0.37344
DELTA_ENERGY Mn 0.36125
DELTA_ENERGY Fe 0.36001
DELTA_ENERGY Co 0.36293
DELTA_ENERGY Ni 0.24380
DELTA_ENERGY Cu 0.40530
DELTA_ENERGY Zn 0.39651
DELTA_ENERGY Ga 0.35002
DELTA_ENERGY Ge 0.34578
DELTA_ENERGY As 0.34953
DELTA_ENERGY Se 0.36731
DELTA_ENERGY Br 0.38201
DELTA_ENERGY Kr 0.39971
&END GCP_POTENTIAL
&END XC
@@ -69,7 +69,7 @@ MODULE bibliography
Putrino2002, Sebastiani2001, Weber2009, Tran2013, Golze2013, &
Golze2015, Golze2017a, Golze2017b, Tuckerman1992, Zhao1994, &
Tozer1996, Goedecker2004, &
Khaliullin2013, Hutter2014, Bengtsson1999, Kantorovich2008, &
Khaliullin2013, Kruse2012, Hutter2014, Bengtsson1999, Kantorovich2008, &
Kantorovich2008a, Wellendorff2012, Niklasson2014, Borstnik2014, &
Rayson2009, Grimme2011, Fattebert2002, Andreussi2012, &
Khaliullin2007, Khaliullin2008, Merlot2014, Lin2009, Lin2013, &
@@ -3983,6 +3983,23 @@ SUBROUTINE add_all_references()
"ER"), &
DOI="10.4208/cicp.OA-2018-0053")

CALL add_reference(key=Kruse2012, ISI_record=s2a( &
"AU Kruse,Holger", &
" Grimme,Stefan", &
"AF Kruse,Holger", &
" Grimme,Stefan", &
"TI A geometrical correction for the inter- and intra-molecular basis set superposition error", &
" in Hartree-Fock and density functional theory calculations for large systems", &
"SO The Journal of Chemical Physics", &
"SN 0021-9606", &
"PY 2012", &
"VL 136", &
"IS 15", &
"AR 154101", &
"DI 10.1063/1.3700154", &
"ER"), &
DOI="10.1063/1.3700154")

END SUBROUTINE add_all_references

END MODULE bibliography
@@ -12,7 +12,7 @@
MODULE input_cp2k_xc
USE bibliography, ONLY: &
Becke1988, Becke1997, BeckeRoussel1989, Goedecker1996, Grimme2006, Grimme2010, Grimme2011, &
Heyd2004, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Perdew1981, Perdew1996, &
Heyd2004, Kruse2012, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Perdew1981, Perdew1996, &
Perdew2008, Proynov2007, Tao2003, Tran2013, Vosko1980, Wellendorff2012, Zhang1998
USE cp_output_handling, ONLY: add_last_numeric,&
cp_print_key_section_create,&
@@ -961,8 +961,16 @@ SUBROUTINE create_vdw_potential_section(section)
usage="SHORT_RANGE_CORRECTION", default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
! KG molecular corrections
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SHORT_RANGE_CORRECTION_PARAMETERS", &
description="Parameters for the short-range bond correction to the DFT-D3 model. "// &
" s*(za*zb)^t1*EXP(-g*dr*r0ab^t2), parameters: s, g, t1, t2 "// &
" Defaults: s=0.08, g=10.0, t1=0.5, t2=-1.0 ", &
usage="SHORT_RANGE_CORRECTION_PARAMETRS", default_r_vals=(/0.08_dp, 10.0_dp, 0.5_dp, -1.0_dp/), &
n_var=4, type_of_var=real_t)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
! KG molecular corrections
CALL keyword_create(keyword, __LOCATION__, name="MOLECULE_CORRECTION", &
description="Calculate a intermolecular correction to the DFT-D3 model", &
usage="MOLECULE_CORRECTION", default_l_val=.FALSE., &
@@ -1075,6 +1083,58 @@ SUBROUTINE create_vdw_potential_section(section)

END SUBROUTINE create_vdw_potential_section

! **************************************************************************************************
!> \brief creates the structure of the section needed for gCP potentials
!> \param section the section that will be created
!> \author jgh
! **************************************************************************************************
SUBROUTINE create_gcp_potential_section(section)
TYPE(section_type), POINTER :: section

CHARACTER(len=*), PARAMETER :: routineN = 'create_gcp_potential_section', &
routineP = moduleN//':'//routineN

TYPE(keyword_type), POINTER :: keyword

CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="gcp_potential", &
description="This section combines geometrical counterpoise potentials."// &
" This is a simple empirical pair potential to correct for BSSE. ", &
citations=(/Kruse2012/), &
n_keywords=1, n_subsections=1, repeats=.FALSE.)

NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="PARAMETER_FILE_NAME", &
description="Name of the parameter file, may include a path", &
usage="PARAMETER_FILE_NAME <FILENAME>", &
default_lc_val="---")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="GLOBAL_PARAMETERS", &
description="Global parameters of the gCP method."// &
" Parameters are sigma, alpha, beta, eta from the original paper.", &
usage="GLOBAL_PARAMETERS 1.0 1.0 1.0 1.0", n_var=4, &
default_r_vals=(/0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DELTA_ENERGY", &
description="Specify the delta energy [Hartree] term for an atom kind", &
usage="DELTA_ENERGY type value", &
type_of_var=char_t, repeats=.TRUE., n_var=-1, default_c_vals=(/"XX ", "0.0"/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="VERBOSE", &
description="Verbose output for gCP calculation", &
usage="VERBOSE logical_value", &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

END SUBROUTINE create_gcp_potential_section

! **************************************************************************************************
!> \brief creates the input section for the xc part
!> \param section the section to create
@@ -1196,6 +1256,10 @@ SUBROUTINE create_xc_section(section)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

CALL create_gcp_potential_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

END SUBROUTINE create_xc_section

! **************************************************************************************************
@@ -776,6 +776,9 @@ SUBROUTINE ec_energy(qs_env, ec_env, unit_nr)
! dispersion through pairpotentials
CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, energy, .FALSE.)
ec_env%edispersion = ec_env%edispersion+energy
!deb ! gCP pairpotentials
!deb CALL calculate_gcp_pairpot(qs_env, ec_env%gcp_env, energy, .FALSE.)
!deb ec_env%edispersion = ec_env%edispersion+energy

SELECT CASE (ec_env%energy_functional)
CASE (kg_ec_functional_harris)
@@ -873,8 +873,8 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, dd, de6, de8, de91, &
de921, de922, dea, devdw, dfdab6, dfdab8, dfdabc, dfdmp, dr, drk, e6, e6tot, e8, e8tot, &
e9, e9tot, elrc6, elrc8, elrc9, eps_cn, er, esrb, evdw, f0ab, fac, fac0, fdab6, fdab8, &
fdabc, fdmp, gnorm, kgc8, nab, nabc, r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, rabc, &
rc2, rcc, rcut, rs6, rs8, s1, s2, s3, s6, s8, s8i, s9, srbe, xp
fdabc, fdmp, gnorm, gsrb, kgc8, nab, nabc, r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, &
rabc, rc2, rcc, rcut, rs6, rs8, s1, s2, s3, s6, s8, s8i, s9, srbe, ssrb, t1srb, t2srb, xp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atom2mol, c6d2, cnkind, cnumbers, &
cnumfix, radd2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
@@ -1035,6 +1035,11 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
idmp = 2
END IF
! SRB parameters
ssrb = dispersion_env%srb_params(1)
gsrb = dispersion_env%srb_params(2)
t1srb = dispersion_env%srb_params(3)
t2srb = dispersion_env%srb_params(4)

IF (unit_nr > 0) THEN
WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
@@ -1049,6 +1054,7 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
IF (dispersion_env%lrc) THEN
WRITE (unit_nr, *) " Apply a long range correction"
WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
END IF
IF (dispersion_env%srb) THEN
WRITE (unit_nr, *) " Apply a short range bond correction"
@@ -1239,7 +1245,7 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
CPABORT("Unknown DFT-D3 damping function:")
END IF
IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
srbe = 0.08_dp*SQRT(REAL((za*zb), KIND=dp))*EXP(-10.0_dp*dr/dispersion_env%r0ab(za, zb))
srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
esrb = esrb+srbe
evdw = evdw-srbe
ELSE
@@ -1267,7 +1273,7 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
END IF
fdij(:) = (de6+de8)*rij(:)/dr*fac
IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
fdij(:) = fdij(:)+10._dp*srbe*rij(:)/(dr*dispersion_env%r0ab(za, zb))
fdij(:) = fdij(:)+srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
END IF
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
@@ -50,6 +50,7 @@ MODULE qs_dispersion_types
LOGICAL :: c9cnst !use constant c9 terms
LOGICAL :: lrc !calculate a long range correction
LOGICAL :: srb !calculate a short range bond correction
REAL(KIND=dp), DIMENSION(4) :: srb_params ! parameters for SRB (s,g,t1,t2)
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: sab_vdw, sab_cn ! neighborlists for pair interactions
REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER &
@@ -113,6 +113,8 @@ SUBROUTINE qs_dispersion_env_set(dispersion_env, xc_section)
CALL section_vals_val_get(pp_section, "REFERENCE_C9_TERM", l_val=dispersion_env%c9cnst)
CALL section_vals_val_get(pp_section, "LONG_RANGE_CORRECTION", l_val=dispersion_env%lrc)
CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION", l_val=dispersion_env%srb)
CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION_PARAMETERS", r_vals=params)
dispersion_env%srb_params(1:4) = params(1:4)
! KG corrections
CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION", l_val=dispersion_env%domol)
CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION_C8", r_val=dispersion_env%kgc8)
@@ -48,6 +48,8 @@ MODULE qs_energy_init
USE qs_external_density, ONLY: external_read_density
USE qs_external_potential, ONLY: external_c_potential,&
external_e_potential
USE qs_gcp_method, ONLY: calculate_gcp_pairpot
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_ks_methods, ONLY: qs_ks_allocate_basics
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
@@ -261,6 +263,7 @@ SUBROUTINE qs_energies_init_hamiltonians(qs_env, calc_forces, molecule_only)
POINTER :: sab_nl
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(section_vals_type), POINTER :: input

CALL timeset(routineN, handle)
@@ -339,6 +342,12 @@ SUBROUTINE qs_energies_init_hamiltonians(qs_env, calc_forces, molecule_only)
! energy info at the end of the SCF
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env, energy=energy)
CALL calculate_dispersion_pairpot(qs_env, dispersion_env, energy%dispersion, calc_forces)
! Add possible pair potential gCP energy - Evaluate first so we can print
! energy info at the end of the SCF
CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env, energy=energy)
IF (ASSOCIATED(gcp_env)) THEN
CALL calculate_gcp_pairpot(qs_env, gcp_env, energy%gcp, calc_forces)
END IF

END IF
! Embedding potential
@@ -27,6 +27,7 @@ MODULE qs_energy_types
core_self, &
repulsive, &
dispersion, &
gcp, &
ex, &
exc, &
exc_aux_fit, &
@@ -154,6 +155,7 @@ SUBROUTINE init_qs_energy(qs_energy)
qs_energy%core_self = 0.0_dp
qs_energy%repulsive = 0.0_dp
qs_energy%dispersion = 0.0_dp
qs_energy%gcp = 0.0_dp
qs_energy%qmmm_el = 0.0_dp
qs_energy%qmmm_nu = 0.0_dp
qs_energy%ex = 0.0_dp
@@ -250,6 +250,7 @@ SUBROUTINE qs_energies_properties(qs_env)
CALL atprop_array_add(atprop%atener, atprop%atexc)
CALL atprop_array_add(atprop%atener, atprop%atecoul)
CALL atprop_array_add(atprop%atener, atprop%atevdw)
CALL atprop_array_add(atprop%atener, atprop%ategcp)
CALL atprop_array_add(atprop%atener, atprop%atecc)
CALL atprop_array_add(atprop%atener, atprop%ate1c)
END IF
@@ -142,6 +142,9 @@ MODULE qs_environment
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_gcp_utils, ONLY: qs_gcp_env_set,&
qs_gcp_init
USE qs_interactions, ONLY: init_interaction_radii,&
init_se_nlradius,&
write_core_charge_radii,&
@@ -516,6 +519,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(qs_ks_env_type), POINTER :: ks_env
@@ -1141,6 +1145,21 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
END IF

! Initialize possible geomertical counterpoise correction potential
IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
dft_control%qs_control%method_id == do_method_gapw .OR. &
dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
dft_control%qs_control%method_id == do_method_lrigpw .OR. &
dft_control%qs_control%method_id == do_method_rigpw .OR. &
dft_control%qs_control%method_id == do_method_ofgpw) THEN
ALLOCATE (gcp_env)
NULLIFY (xc_section)
xc_section => section_vals_get_subs_vals(dft_section, "XC")
CALL qs_gcp_env_set(gcp_env, xc_section)
CALL qs_gcp_init(qs_env, gcp_env)
CALL set_qs_env(qs_env, gcp_env=gcp_env)
END IF

! *** Allocate the MO data types ***
CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)

0 comments on commit ec01bc2

Please sign in to comment.
You can’t perform that action at this time.