Octopus
xc_vdw.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2023 . Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
23module xc_vdw_oct_m
24 use atom_oct_m
25 use debug_oct_m
26 use global_oct_m
27 use grid_oct_m
28 use ions_oct_m
29 use, intrinsic :: iso_fortran_env
35 use parser_oct_m
37 use space_oct_m
41 use vdw_ts_oct_m
42 use xc_oct_m
45
46 ! from the dftd3 library
47 use dftd3_api
48
49 implicit none
50
51 private
52 public :: &
53 xc_vdw_t, &
55
56 type xc_vdw_t
57 private
58 integer, public :: vdw_correction = -1
59 logical :: self_consistent = .false.
60 type(vdw_ts_t), public :: vdw_ts
61 type(dftd3_calc) :: vdw_d3
62 real(real64), public, allocatable :: forces(:, :)
63
64 real(real64), public :: stress(3, 3)
65
66 contains
67 procedure :: init => xc_vdw_init
68 procedure :: calc => xc_vdw_calc
69 procedure :: d3_version => xc_vdw_octopus_input_to_version
70 procedure :: end => xc_vdw_end
71 end type xc_vdw_t
72
74 ! No reason for this to be int64, but options are appended with this
75 ! precision when parsed.
76 integer(int64), parameter :: d3_lib_options(5) = [&
77 option__vdwcorrection__vdw_d2, &
78 option__vdwcorrection__vdw_d3_zero_damping, &
79 option__vdwcorrection__vdw_d3_bj, &
80 option__vdwcorrection__vdw_d3m_zero_damping, &
81 option__vdwcorrection__vdw_d3m_bj &
82 ]
83
84contains
85
87 subroutine xc_vdw_init(this, namespace, space, gr, xc, ions, x_id, c_id)
88 class(xc_vdw_t), intent(inout) :: this
89 type(namespace_t), intent(in) :: namespace
90 class(space_t), intent(in) :: space
91 type(grid_t), intent(in) :: gr
92 type(xc_t), intent(inout) :: xc
93 type(ions_t), intent(in) :: ions
94 integer, intent(in) :: x_id, c_id
95
96 integer :: iatom
97 type(dftd3_input) :: d3_input
98 character(len=20) :: d3func_def, d3func
99
100 push_sub(xc_vdw_init)
101
102 if (in_family(xc%family, [xc_family_libvdwxc])) then
103 call libvdwxc_set_geometry(xc%functional(func_c,1)%libvdwxc, namespace, space, gr)
104 end if
105
106 !%Variable VDWCorrection
107 !%Type integer
108 !%Default no
109 !%Section Hamiltonian::XC
110 !%Description
111 !% (Experimental) This variable selects which van der Waals
112 !% correction to apply to the correlation functional.
113 !%Option none 0
114 !% No correction is applied.
115 !%Option vdw_ts 1
116 !% The scheme of Tkatchenko and Scheffler, Phys. Rev. Lett. 102
117 !% 073005 (2009).
118 !%Option vdw_d2 2
119 !% DFT-D2 corrections, equivalent to specifying version = 2 in the DFTD3 library.
120 !%Option vdw_d3_zero_damping 3
121 !% DFT-D3 corrections, damping the dispersion contribution to zero for short ranges.
122 !% This is equivalent to specifying version = 3 in the DFTD3 library.
123 !% For more details, see [J. Chem. Phys. 132, 154104 (2010)](https://doi.org/10.1063/1.3382344).
124 !%Option vdw_d3_bj 4
125 !% DFT-D3 with Becke-Johnson finite-damping, variant 2. In Octopus 15.0 and older,
126 !% his was only option with the DFT-D3 library, and was
127 !% specified with "vdw_d3".
128 !% This is equivalent to specifying version = 4 in the DFTD3 library.
129 !% For more details, see
130 !% [Effect of the damping function in dispersion corrected density functional theory](https://doi.org/10.1002/jcc.21759).
131 !%Option vdw_d3m_zero_damping 5
132 !% DFT-D3 corrections, damping the dispersion contribution to zero for short ranges.
133 !% This uses the modified damping function and parameters
134 !% of Sherrill and coworkers. For more details see
135 !% [J. Phys. Chem. Lett. 2016, 7, 12, 2197–2203](https://doi.org/10.1021/acs.jpclett.6b00780)
136 !% This is equivalent to specifying version = 5 in the DFTD3 library.
137 !%Option vdw_d3m_bj 6
138 !% DFT-D3 with Becke-Johnson finite-damping, and the modified parameters of Sherrill and coworkers.
139 !% For more details see [J. Phys. Chem. Lett. 2016, 7, 12, 2197–2203](https://doi.org/10.1021/acs.jpclett.6b00780)
140 !% This is equivalent to specifying version = 6 in the DFTD3 library.
141 !%End
142 call parse_variable(namespace, 'VDWCorrection', option__vdwcorrection__none, this%vdw_correction)
143
144 if (this%vdw_correction == option__vdwcorrection__none) then
145 this%self_consistent = .false.
146 pop_sub(xc_vdw_init)
147 return
148 endif
150 call messages_experimental('VDWCorrection')
152 if( this%vdw_correction == option__vdwcorrection__vdw_ts) then
153 !%Variable VDWSelfConsistent
154 !%Type logical
155 !%Default yes
156 !%Section Hamiltonian::XC
157 !%Description
158 !% This variable controls whether the VDW correction is applied
159 !% self-consistently, the default, or just as a correction to
160 !% the total energy. This option only works with vdw_ts.
161 !%End
162 call parse_variable(namespace, 'VDWSelfConsistent', .true., this%self_consistent)
164 call vdw_ts_init(this%vdw_ts, namespace, ions)
165
166 ! DFTD3 library
167 else
168 this%self_consistent = .false.
170 if (space%dim /= 3) then
171 call messages_write('xc_vdw_d3 can only be used in 3-dimensional systems')
172 call messages_fatal(namespace=namespace)
173 end if
174
175 do iatom = 1, ions%natoms
176 if (.not. ions%atom(iatom)%species%represents_real_atom()) then
177 call messages_write('xc_vdw_d3 is not implemented when non-atomic species are present')
178 call messages_fatal(namespace=namespace)
179 end if
180 end do
181
182 d3func_def = ''
183
184 ! The list of valid values can be found in 'external_libs/dftd3/core.f90'.
185 ! For the moment I include the most common ones.
186 if (x_id == option__xcfunctional__gga_x_b88 .and. c_id*libxc_c_index == option__xcfunctional__gga_c_lyp) then
187 d3func_def = 'b-lyp'
188 end if
189 if (x_id == option__xcfunctional__gga_x_pbe .and. c_id*libxc_c_index == option__xcfunctional__gga_c_pbe) then
190 d3func_def = 'pbe'
191 end if
192 if (x_id == option__xcfunctional__gga_x_pbe_sol .and. c_id*libxc_c_index == option__xcfunctional__gga_c_pbe_sol) then
193 d3func_def = 'pbesol'
194 end if
195 if (c_id*libxc_c_index == option__xcfunctional__hyb_gga_xc_b3lyp) then
196 d3func_def = 'b3-lyp'
197 end if
198 if (c_id*libxc_c_index == option__xcfunctional__hyb_gga_xc_pbeh) then
199 d3func_def = 'pbe0'
200 end if
201
202 !%Variable VDWD3Functional
203 !%Type string
204 !%Section Hamiltonian::XC
205 !%Description
206 !% (Experimental) You can use this variable to override the
207 !% parametrization used by the DFT-D3 van deer Waals
208 !% correction. Normally you need not set this variable, as the
209 !% proper value will be selected by Octopus (if available).
210 !%
211 !% This variable takes a string value, the valid values can
212 !% be found in the source file 'external_libs/dftd3/core.f90'.
213 !% For example you can use:
214 !%
215 !% VDWD3Functional = 'pbe'
216 !%
217 !%End
218 if (parse_is_defined(namespace, 'VDWD3Functional')) call messages_experimental('VDWD3Functional')
219 call parse_variable(namespace, 'VDWD3Functional', d3func_def, d3func)
220
221 if (d3func == '') then
222 call messages_write('Cannot find a matching parametrization of DFT-D3 for the current')
223 call messages_new_line()
224 call messages_write('XCFunctional. Please select a different XCFunctional, or select a')
225 call messages_new_line()
226 call messages_write('functional for DFT-D3 using the <tt>VDWD3Functional</tt> variable.')
227 call messages_fatal(namespace=namespace)
228 end if
229
230 if (space%periodic_dim /= 0 .and. space%periodic_dim /= 3) then
231 call messages_write('For partially periodic systems, the xc_vdw_d3 interaction is assumed')
232 call messages_new_line()
233 call messages_write('to be periodic in three dimensions.')
234 call messages_warning(namespace=namespace)
235 end if
236
237 call dftd3_init(this%vdw_d3, d3_input, trim(conf%share)//'/dftd3/pars.dat')
238 call dftd3_set_functional(this%vdw_d3, func = d3func, version = this%d3_version(), tz = .false.)
239 endif
240
241 pop_sub(xc_vdw_init)
242 end subroutine xc_vdw_init
243
244
246 function xc_vdw_octopus_input_to_version(this) result(d3_version)
247 class(xc_vdw_t), intent(in) :: this
248
249 integer :: d3_version
250
252
253 select case (this%vdw_correction)
254 case (option__vdwcorrection__vdw_d2)
255 d3_version = 2
256 case (option__vdwcorrection__vdw_d3_zero_damping)
257 d3_version = 3
258 case (option__vdwcorrection__vdw_d3_bj)
259 d3_version = 4
260 case (option__vdwcorrection__vdw_d3m_zero_damping)
261 d3_version = 5
262 case (option__vdwcorrection__vdw_d3m_bj)
263 d3_version = 6
264 case default
265 assert(.false.)
266 end select
267
269
271
272
273 ! ---------------------------------------------------------
274 subroutine xc_vdw_end(this)
275 class(xc_vdw_t), intent(inout) :: this
276
277 push_sub(xc_vdw_end)
278
279 select case (this%vdw_correction)
280 case (option__vdwcorrection__vdw_ts)
281 call vdw_ts_end(this%vdw_ts)
282 end select
283
284 pop_sub(xc_vdw_end)
285 end subroutine xc_vdw_end
286
287 ! ---------------------------------------------------------
288 subroutine xc_vdw_calc(this, namespace, space, latt, atom, natoms, pos, gr, st, energy, vxc)
289 class(xc_vdw_t), intent(inout) :: this
290 type(namespace_t), intent(in) :: namespace
291 type(space_t), intent(in) :: space
292 type(lattice_vectors_t), intent(in) :: latt
293 type(atom_t), intent(in) :: atom(:)
294 integer, intent(in) :: natoms
295 real(real64), intent(in) :: pos(1:space%dim,1:natoms)
296 type(grid_t), intent(in) :: gr
297 type(states_elec_t), intent(inout) :: st
298 real(real64), intent(out) :: energy
299 real(real64), contiguous, intent(inout) :: vxc(:,:)
300
301 integer :: ispin, iatom, idir
302 real(real64), allocatable :: vxc_vdw(:)
303 integer, allocatable :: atnum(:)
304 real(real64) :: rlattice(3, 3)
305 ! For non-periodic directions, we make the length of the box very large
306 real(real64), parameter :: aperiodic_scaling = 1000.0_real64
307
308 energy = m_zero
309 if (this%vdw_correction == option__vdwcorrection__none) return
310
311 push_sub(xc_vdw_calc)
312
313 assert(space%dim == 3)
314
315 safe_allocate(vxc_vdw(1:gr%np))
316 safe_allocate(this%forces(1:space%dim, 1:natoms))
317
318 if(this%vdw_correction == option__vdwcorrection__vdw_ts) then
319 vxc_vdw = m_zero
320 call vdw_ts_calculate(this%vdw_ts, namespace, space, latt, atom, natoms, pos, &
321 gr, st%d%nspin, st%rho, energy, vxc_vdw, this%forces)
322 this%stress = m_zero
323
324 elseif(any(this%vdw_correction == d3_lib_options)) then
325
326 safe_allocate(atnum(1:natoms))
327
328 atnum = [(nint(atom(iatom)%species%get_z()), iatom=1, natoms)]
329
330 if (space%is_periodic()) then
331 ! DFTD3 treats interactions as 3D periodic. In the case of
332 ! mixed-periodicity, we set the length of the periodic lattice along
333 ! the non-periodic dimensions to be much larger than the system
334 ! size, such that for all practical purposes the system is treated
335 ! as isolated.
336 rlattice = latt%rlattice(1:3, 1:3)
337 if (space%periodic_dim < 3) then
338 message(1) = "Info: Using DFT-D3 van der Walls corrections for a system with mixed-periodicity,"
339 message(2) = " but DFTD3 treats system as fully periodic. Octopus will set the DFT-D3"
340 message(3) = " lattice vectors along non-periodic dimensions to a suitably large value."
341 call messages_info(3, namespace=namespace)
342
343 do idir = space%periodic_dim + 1, 3
344 rlattice(idir,idir) = (maxval(abs(pos(idir,:))) + m_one)*aperiodic_scaling
345 end do
346 end if
347 call dftd3_pbc_dispersion(this%vdw_d3, pos, atnum, rlattice, energy, this%forces, &
348 this%stress)
349 else
350 call dftd3_dispersion(this%vdw_d3, pos, atnum, energy, this%forces)
351 this%stress = m_zero
352 end if
353
354 safe_deallocate_a(atnum)
355
356 else
357 assert(.false.)
358 endif
359
360 if (this%self_consistent) then
361 do ispin = 1, st%d%nspin
362 call lalg_axpy(gr%np, m_one, vxc_vdw, vxc(:, ispin))
363 end do
364 end if
365
366 safe_deallocate_a(vxc_vdw)
368 pop_sub(xc_vdw_calc)
369 end subroutine xc_vdw_calc
370
371end module xc_vdw_oct_m
372
373!! Local Variables:
374!! mode: f90
375!! coding: utf-8
376!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:171
real(real64) function func(r1, rn, n, a)
Definition: logrid.F90:353
real(real64), parameter, public m_zero
Definition: global.F90:188
type(conf_t), public conf
Global instance of Octopus configuration.
Definition: global.F90:178
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public libvdwxc_set_geometry(this, namespace, space, mesh)
Definition: libvdwxc.F90:229
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
integer, parameter, private libxc_c_index
Definition: species.F90:276
Tkatchenko-Scheffler pairwise method for van der Waals (vdW, dispersion) interactions.
Definition: vdw_ts.F90:119
subroutine, public vdw_ts_init(this, namespace, ions)
Definition: vdw_ts.F90:173
subroutine, public vdw_ts_calculate(this, namespace, space, latt, atom, natoms, pos, mesh, nspin, density, energy, potential, force)
Calculate the potential for the Tatchenko-Scheffler vdW correction as described in Phys....
Definition: vdw_ts.F90:267
subroutine, public vdw_ts_end(this)
Definition: vdw_ts.F90:249
integer, parameter, public xc_family_libvdwxc
integer, parameter, public func_c
Definition: xc.F90:114
pure logical function, public in_family(family, xc_families)
Definition: xc.F90:623
A module that takes care of xc contribution from vdW interactions.
Definition: xc_vdw.F90:116
subroutine xc_vdw_end(this)
Definition: xc_vdw.F90:368
subroutine xc_vdw_init(this, namespace, space, gr, xc, ions, x_id, c_id)
Initialize a xc_vdw_t instance.
Definition: xc_vdw.F90:181
integer function xc_vdw_octopus_input_to_version(this)
Map Octopus''s VDWCORRECTION input option to the DFTD3 library''s input option "version".
Definition: xc_vdw.F90:340
integer(int64), dimension(5), parameter, public d3_lib_options
VDWCORRECTION options that correspond to the DFT-D3 library.
Definition: xc_vdw.F90:169
subroutine xc_vdw_calc(this, namespace, space, latt, atom, natoms, pos, gr, st, energy, vxc)
Definition: xc_vdw.F90:382
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:169
The states_elec_t class contains all electronic wave functions.
int true(void)