Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Set disperson types for DFT calculations
10 : !> \author JGH (04.2014)
11 : ! **************************************************************************************************
12 : MODULE qs_dispersion_utils
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE input_constants, ONLY: &
21 : vdw_nl_DRSLL, vdw_nl_LMKLL, vdw_nl_RVV10, vdw_pairpot_dftd2, vdw_pairpot_dftd3, &
22 : vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
23 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
24 : section_vals_type,&
25 : section_vals_val_get
26 : USE kinds, ONLY: dp
27 : USE physcon, ONLY: bohr,&
28 : kjmol
29 : USE qs_dispersion_pairpot, ONLY: qs_scaling_dftd3,&
30 : qs_scaling_dftd3bj,&
31 : qs_scaling_init
32 : USE qs_dispersion_types, ONLY: qs_atom_dispersion_type,&
33 : qs_dispersion_type
34 : USE qs_environment_types, ONLY: get_qs_env,&
35 : qs_environment_type
36 : USE qs_kind_types, ONLY: get_qs_kind,&
37 : qs_kind_type
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_utils'
45 :
46 : PUBLIC :: qs_dispersion_env_set, qs_write_dispersion
47 :
48 : ! **************************************************************************************************
49 : CONTAINS
50 : ! **************************************************************************************************
51 : !> \brief ...
52 : !> \param dispersion_env ...
53 : !> \param xc_section ...
54 : ! **************************************************************************************************
55 5350 : SUBROUTINE qs_dispersion_env_set(dispersion_env, xc_section)
56 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
57 : TYPE(section_vals_type), POINTER :: xc_section
58 :
59 : LOGICAL :: exfun, explicit
60 5350 : REAL(dp), POINTER :: params(:), scal(:)
61 : TYPE(section_vals_type), POINTER :: nl_section, pp_section, vdw_section, &
62 : xc_fun_section
63 :
64 0 : CPASSERT(ASSOCIATED(dispersion_env))
65 :
66 : ! set general defaults
67 5350 : dispersion_env%doabc = .FALSE.
68 5350 : dispersion_env%c9cnst = .FALSE.
69 5350 : dispersion_env%lrc = .FALSE.
70 5350 : dispersion_env%srb = .FALSE.
71 5350 : dispersion_env%verbose = .FALSE.
72 5350 : dispersion_env%nd3_exclude_pair = 0
73 5350 : NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
74 5350 : dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
75 5350 : dispersion_env%d3_exclude_pair)
76 5350 : NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
77 5350 : dispersion_env%d2y_dx2)
78 5350 : NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
79 5350 : NULLIFY (dispersion_env%dftd_section)
80 5350 : NULLIFY (vdw_section, xc_fun_section)
81 5350 : vdw_section => section_vals_get_subs_vals(xc_section, "vdw_potential")
82 5350 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
83 5350 : CALL section_vals_val_get(vdw_section, "POTENTIAL_TYPE", i_val=dispersion_env%type)
84 5350 : IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
85 80 : NULLIFY (pp_section)
86 80 : pp_section => section_vals_get_subs_vals(vdw_section, "PAIR_POTENTIAL")
87 80 : CALL section_vals_val_get(pp_section, "VERBOSE_OUTPUT", l_val=dispersion_env%verbose)
88 80 : CALL section_vals_val_get(pp_section, "TYPE", i_val=dispersion_env%pp_type)
89 80 : IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
90 : ! functional parameters for Grimme D2 type
91 20 : CALL section_vals_val_get(pp_section, "EXP_PRE", r_val=dispersion_env%exp_pre)
92 20 : CALL section_vals_val_get(pp_section, "SCALING", explicit=explicit)
93 20 : IF (.NOT. explicit) THEN
94 8 : CALL section_vals_val_get(pp_section, "REFERENCE_FUNCTIONAL", explicit=exfun)
95 8 : CPASSERT(exfun)
96 8 : CALL qs_scaling_init(dispersion_env%scaling, vdw_section)
97 : ELSE
98 12 : CALL section_vals_val_get(pp_section, "SCALING", r_val=dispersion_env%scaling)
99 : END IF
100 : ELSE
101 60 : dispersion_env%exp_pre = 0._dp
102 60 : dispersion_env%scaling = 0._dp
103 : END IF
104 80 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
105 : dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
106 : ! functional parameters for Grimme DFT-D3 type
107 54 : CALL section_vals_val_get(pp_section, "EPS_CN", r_val=dispersion_env%eps_cn)
108 54 : CALL section_vals_val_get(pp_section, "CALCULATE_C9_TERM", l_val=dispersion_env%doabc)
109 54 : CALL section_vals_val_get(pp_section, "REFERENCE_C9_TERM", l_val=dispersion_env%c9cnst)
110 54 : CALL section_vals_val_get(pp_section, "LONG_RANGE_CORRECTION", l_val=dispersion_env%lrc)
111 54 : CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION", l_val=dispersion_env%srb)
112 54 : CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION_PARAMETERS", r_vals=params)
113 540 : dispersion_env%srb_params(1:4) = params(1:4)
114 : ! KG corrections
115 54 : CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION", l_val=dispersion_env%domol)
116 54 : CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION_C8", r_val=dispersion_env%kgc8)
117 54 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
118 42 : CALL section_vals_val_get(pp_section, "D3_SCALING", explicit=explicit)
119 12 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
120 12 : CALL section_vals_val_get(pp_section, "D3BJ_SCALING", explicit=explicit)
121 : END IF
122 54 : IF (.NOT. explicit) THEN
123 48 : CALL section_vals_val_get(pp_section, "REFERENCE_FUNCTIONAL", explicit=exfun)
124 48 : CPASSERT(exfun)
125 48 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
126 38 : CALL qs_scaling_dftd3(dispersion_env%s6, dispersion_env%sr6, dispersion_env%s8, vdw_section)
127 10 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
128 : CALL qs_scaling_dftd3bj(dispersion_env%s6, dispersion_env%a1, dispersion_env%s8, &
129 10 : dispersion_env%a2, vdw_section)
130 : END IF
131 : ELSE
132 6 : IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
133 : ! zero damping
134 4 : CALL section_vals_val_get(pp_section, "D3_SCALING", r_vals=scal)
135 4 : dispersion_env%s6 = scal(1)
136 4 : dispersion_env%sr6 = scal(2)
137 4 : dispersion_env%s8 = scal(3)
138 4 : dispersion_env%a1 = 0.0_dp
139 4 : dispersion_env%a2 = 0.0_dp
140 2 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
141 : ! BJ damping
142 2 : CALL section_vals_val_get(pp_section, "D3BJ_SCALING", r_vals=scal)
143 2 : dispersion_env%s6 = scal(1)
144 2 : dispersion_env%a1 = scal(2)
145 2 : dispersion_env%s8 = scal(3)
146 2 : dispersion_env%a2 = scal(4)
147 2 : dispersion_env%sr6 = 0.0_dp
148 : END IF
149 : END IF
150 : ELSE
151 26 : dispersion_env%s6 = 0._dp
152 26 : dispersion_env%sr6 = 0._dp
153 26 : dispersion_env%s8 = 0._dp
154 26 : dispersion_env%a1 = 0._dp
155 26 : dispersion_env%a2 = 0._dp
156 26 : dispersion_env%eps_cn = 0._dp
157 : END IF
158 80 : IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
159 6 : CALL section_vals_val_get(pp_section, "REFERENCE_FUNCTIONAL", explicit=exfun)
160 6 : CPASSERT(exfun)
161 : CALL section_vals_val_get(vdw_section, "PAIR_POTENTIAL%REFERENCE_FUNCTIONAL", &
162 6 : c_val=dispersion_env%ref_functional)
163 6 : CALL section_vals_val_get(pp_section, "D4_CUTOFF", r_val=dispersion_env%rc_d4)
164 6 : CALL section_vals_val_get(pp_section, "D4_CN_CUTOFF", r_val=dispersion_env%rc_cn)
165 6 : CALL section_vals_val_get(pp_section, "FACTOR_S9_TERM", r_val=dispersion_env%s9)
166 : END IF
167 80 : CALL section_vals_val_get(pp_section, "R_CUTOFF", r_val=dispersion_env%rc_disp)
168 : CALL section_vals_val_get(pp_section, "PARAMETER_FILE_NAME", &
169 80 : c_val=dispersion_env%parameter_file_name)
170 : ! set DFTD section for output handling
171 80 : dispersion_env%dftd_section => pp_section
172 5270 : ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
173 46 : NULLIFY (nl_section)
174 46 : nl_section => section_vals_get_subs_vals(vdw_section, "NON_LOCAL")
175 46 : CALL section_vals_val_get(nl_section, "VERBOSE_OUTPUT", l_val=dispersion_env%verbose)
176 : CALL section_vals_val_get(nl_section, "KERNEL_FILE_NAME", &
177 46 : c_val=dispersion_env%kernel_file_name)
178 46 : CALL section_vals_val_get(nl_section, "TYPE", i_val=dispersion_env%nl_type)
179 46 : CALL section_vals_val_get(nl_section, "CUTOFF", r_val=dispersion_env%pw_cutoff)
180 46 : CALL section_vals_val_get(nl_section, "PARAMETERS", r_vals=params)
181 46 : CALL section_vals_val_get(nl_section, "SCALE", r_val=dispersion_env%scale_rvv10)
182 46 : dispersion_env%b_value = params(1)
183 46 : dispersion_env%c_value = params(2)
184 : END IF
185 5350 : END SUBROUTINE qs_dispersion_env_set
186 :
187 : ! **************************************************************************************************
188 :
189 : ! **************************************************************************************************
190 : !> \brief ...
191 : !> \param qs_env ...
192 : !> \param dispersion_env ...
193 : !> \param ounit ...
194 : ! **************************************************************************************************
195 5586 : SUBROUTINE qs_write_dispersion(qs_env, dispersion_env, ounit)
196 : TYPE(qs_environment_type), POINTER :: qs_env
197 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
198 : INTEGER, INTENT(in), OPTIONAL :: ounit
199 :
200 : CHARACTER(LEN=2) :: symbol
201 : INTEGER :: i, ikind, nkind, output_unit
202 5586 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
203 : TYPE(cp_logger_type), POINTER :: logger
204 : TYPE(qs_atom_dispersion_type), POINTER :: disp
205 5586 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
206 : TYPE(section_vals_type), POINTER :: dft_section
207 :
208 5586 : IF (PRESENT(ounit)) THEN
209 0 : output_unit = ounit
210 : ELSE
211 5586 : NULLIFY (logger)
212 5586 : logger => cp_get_default_logger()
213 :
214 5586 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
215 : output_unit = cp_print_key_unit_nr(logger, dft_section, &
216 5586 : "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
217 : END IF
218 :
219 5586 : IF (output_unit > 0) THEN
220 : ! vdW type specific output
221 1394 : IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
222 58 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T67,'Pair Potential')")
223 : ! Pair potentials
224 58 : IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
225 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'DFT-D2')")
226 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Potential Form: S. Grimme, JCC 27: 1787 (2006)')")
227 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp
228 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Scaling Factor:',T73,F8.4)") dispersion_env%scaling
229 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Exp Prefactor for Damping:',T73,F8.1)") dispersion_env%exp_pre
230 8 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
231 8 : nkind = SIZE(atomic_kind_set)
232 19 : DO ikind = 1, nkind
233 11 : CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol)
234 11 : CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp)
235 19 : IF (disp%defined) THEN
236 : WRITE (output_unit, fmt="(' vdW PARAMETER| ',T18,'Atom=',A2, "// &
237 : "T28,'C6[J*nm^6*mol^-1]=',F8.4,T63,'r(vdW)[A]=',F8.4)") &
238 11 : symbol, disp%c6/(1000._dp*bohr**6/kjmol), disp%vdw_radii/bohr
239 : ELSE
240 0 : WRITE (output_unit, fmt="(' vdW PARAMETER| ',T20,'Atom=',A2,T70,'not defined')")
241 : END IF
242 : END DO
243 50 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
244 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D3 (Version 3.1)')")
245 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Potential Form: S. Grimme et al, JCP 132: 154104 (2010)')")
246 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Zero Damping')")
247 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp
248 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s6 Scaling Factor:',T73,F8.4)") dispersion_env%s6
249 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'sr6 Scaling Factor:',T73,F8.4)") dispersion_env%sr6
250 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8
251 8 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn
252 8 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
253 0 : DO i = 1, dispersion_env%nd3_exclude_pair
254 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Excluded Pairs: ',T76,I2,' ',I2)") &
255 0 : dispersion_env%d3_exclude_pair(i, :)
256 : END DO
257 : END IF
258 42 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
259 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D3 (Version 3.1)')")
260 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Potential Form: S. Grimme et al, JCP 132: 154104 (2010)')")
261 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'BJ Damping: S. Grimme et al, JCC 32: 1456 (2011)')")
262 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp
263 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s6 Scaling Factor:',T73,F8.4)") dispersion_env%s6
264 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'a1 Damping Factor:',T73,F8.4)") dispersion_env%a1
265 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8
266 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'a2 Damping Factor:',T73,F8.4)") dispersion_env%a2
267 39 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn
268 39 : IF (dispersion_env%nd3_exclude_pair > 0) THEN
269 0 : DO i = 1, dispersion_env%nd3_exclude_pair
270 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Excluded Kind Pairs: ',T76,I2,' ',I2)") &
271 0 : dispersion_env%d3_exclude_pair(i, :)
272 : END DO
273 : END IF
274 3 : ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
275 3 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D4(Version 3.6.0)')")
276 3 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'see https://github.com/dftd4/dftd4')")
277 3 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'E. Caldeweyher et al, PCCP 22: 8499 (2020)')")
278 3 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'E. Caldeweyher et al, JCP 150: 154122 (2019)')")
279 3 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'E. Caldeweyher et al, JCP 147: 034112 (2017)')")
280 : END IF
281 1336 : ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
282 13 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T61,'Non-local Functional')")
283 : WRITE (output_unit, &
284 13 : fmt="(' vdW POTENTIAL| ','Implementation: G. Roman-Perez, J. Soler, PRL 103: 096102 (2009)')")
285 : WRITE (output_unit, &
286 13 : fmt="(' vdW POTENTIAL| ',T38,' T. Thonhauser et al, PRB 76: 125112 (2007)')")
287 : WRITE (output_unit, &
288 13 : fmt="(' vdW POTENTIAL| ',T22,' R. Sabatini et al, J.Phys:Condens Matter 24: 424209 (2012)')")
289 : WRITE (output_unit, &
290 13 : fmt="(' vdW POTENTIAL| ',T16,' Based on QE implementation by Brian Kolb, Timo Thonhauser (2009)')")
291 13 : SELECT CASE (dispersion_env%nl_type)
292 : CASE DEFAULT
293 : ! unknown functional
294 0 : CPABORT("")
295 : CASE (vdw_nl_DRSLL)
296 : WRITE (output_unit, &
297 8 : fmt="(' vdW POTENTIAL| ','DRSLL Functional: M. Dion et al, PRL 92: 246401 (2004)')")
298 : CASE (vdw_nl_LMKLL)
299 : WRITE (output_unit, &
300 3 : fmt="(' vdW POTENTIAL| ','LMKLL Functional: K. Lee et al, PRB 82: 081101 (2010)')")
301 : CASE (vdw_nl_RVV10)
302 : WRITE (output_unit, &
303 13 : fmt="(' vdW POTENTIAL| ','RVV10 Functional: R. Sabatini et al, PRB 87: 041108(R) (2013)')")
304 : END SELECT
305 13 : IF (dispersion_env%verbose) THEN
306 : WRITE (output_unit, &
307 12 : fmt="(' vdW POTENTIAL| ',' Carrying out vdW-DF run using the following parameters:')")
308 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ','Nqs =',I8,' Nr_points =',I8,' r_max =',F10.3)") &
309 12 : dispersion_env%nqs, dispersion_env%nr_points, dispersion_env%r_max
310 12 : WRITE (output_unit, fmt="(' vdW POTENTIAL| ','q_mesh =')")
311 252 : WRITE (output_unit, fmt="(8X,4F18.8)") (dispersion_env%q_mesh(i), i=1, dispersion_env%nqs)
312 : WRITE (output_unit, &
313 : fmt="(' vdW POTENTIAL| ','Density cutoff for convolution [a.u.]:',T71,F10.1)") &
314 12 : dispersion_env%pw_cutoff
315 : END IF
316 : END IF
317 : END IF
318 5586 : IF (.NOT. PRESENT(ounit)) THEN
319 : CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
320 5586 : "PRINT%DFT_CONTROL_PARAMETERS")
321 : END IF
322 :
323 5586 : END SUBROUTINE qs_write_dispersion
324 :
325 : ! **************************************************************************************************
326 :
327 : END MODULE qs_dispersion_utils
328 :
|