1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
|
!
! Copyright (C) 2001-2009 Quantum ESPRESSO group
! Copyright (C) 2009 Brian Kolb, Timo Thonhauser - Wake Forest University
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
MODULE kernel_table
!! This module is used to read in the kernel table file
!! "vdW_kernel_table" and store some of the important parameters. The
!! top of the vdW_kernel_table file holds the number of q points, the
!! number of radial points (r) used in the kernel generation, the maximum
!! value of r used (r is the parameter in the kernel function d=q*r where
!! q is defined in DION equation 11), and the values of the q points
!! used. These parameters are stored as public parameters for use in
!! various routines. This routine also reads the tabulated values of the
!! Fourier transformed kernel function for each pair of q values (see
!! SOLER equations 3 and 8). Since these kernel functions need to be
!! interpolated using splines, the second derivatives of the Fourier
!! transformed kernel functions (phi_alpha_beta) are also tabulated in
!! the vdW_kernel_table and are read in here.
!! This is done in a module because there are quite a few subroutines in
!! xc_vdW_DF.f90 that require knowledge of the number (identity) of q
!! points, the maximum value of the radius, and, of course, the tabulated
!! kernel function and its second derivatives (for spline interpolation).
!! Putting this routine in a module meas that those routines can just use
!! kernel_table rather than passing variables around all over the place.
USE kinds, ONLY : dp
USE io_files, ONLY : pseudo_dir
USE constants, ONLY : pi
use wrappers, ONLY : md5_from_file
implicit none
private
save
!! Variables to be used by various routines in xc_vdW_DF.f90, declared
!! public so they can be seen from outside
! --------------------------------------------------------------------------
public :: Nqs, Nr_points, r_max, q_mesh, q_cut, q_min, dk
public :: kernel, d2phi_dk2
public :: initialize_kernel_table
public :: vdw_table_name
public :: vdw_kernel_md5_cksum
integer :: Nqs, Nr_points !! The number of q points and radial points
! !! used in generating the kernel phi(q1*r, q2*r)
! !! (see DION 14-16 and SOLER 3)
real(dp) :: r_max, q_cut, q_min, dk !! The maximum value of r, the maximum and minimum
! !! values of q and the k-space spacing of grid points.
! !! Note that, during a vdW run, values of q0 found
! !! larger than q_cut will be saturated (SOLER 5) to
! !! q_cut
real(dp), allocatable :: q_mesh(:) !! The values of all the q points used
real(dp), allocatable :: kernel(:,:,:) !! A matrix holding the Fourier transformed kernel function
! !! for each pair of q values. The ordering is
! !! kernel(k_point, q1_value, q2_value)
real(dp), allocatable :: d2phi_dk2(:,:,:) !! A matrix holding the second derivatives of the above
! !! kernel matrix at each of the q points. Stored as
! !! d2phi_dk2(k_point, q1_value, q2_value)
!
character(len=256) :: vdw_table_name !! If present from input use this name
CHARACTER(LEN=30) :: double_format = "(1p4e23.14)"
CHARACTER(len=32) :: vdw_kernel_md5_cksum = 'NOT SET'
!
INTEGER, EXTERNAL :: find_free_unit
! --------------------------------------------------------------------------
CONTAINS
!! ###################################################################################
!! | |
!! | INITIALIZE_KERNEL_TABLE |
!! |___________________________|
!! Subroutine that actually reads the kernel file and stores the parameters. This routine
!! is called only once, at the start of a vdW run.
subroutine initialize_kernel_table(inlc)
integer, INTENT(IN) :: inlc
integer :: q1_i, q2_i !! Indexing variables
integer :: kernel_file !! The unit number for the kernel file
logical :: file_exists !! A variable to say whether
! !! needed file exists
character(len=1000) :: kernel_file_name !! The path to the kernel file.
! !! Although this name must be
! !! "vdW_kernel_table", this variable
! !! is used to hold the entire path
! !! since we check 3 places for it.
!!write(*,*) "Reading kernel table ... "
!! Get the unit number for the kernel file
kernel_file = find_free_unit()
!!
!!
if (TRIM(vdw_table_name)==' ') then
if (inlc==3) then
vdw_table_name='rVV10_kernel_table'
else
vdw_table_name='vdW_kernel_table'
endif
endif
if (allocated(kernel)) then
return
end if
!! First we check the current directory for the vdW_kernel_table file
!! If it is not found there it is looked for in the pseudopotential
!! directory. If it's not there the default kernel file installed
!! in the PW directory of the PWSCF source is tried. If none of those
!! exist the code crashes.
kernel_file_name=vdw_table_name
inquire(file=kernel_file_name, exist=file_exists)
!! If the file is found in the current directory we use that one
!! ------------------------------------------------------------------------------------------
if (.not. file_exists) then
!! No "vdW_kernel_table" file in the current directory. Try the pseudopotential directory
!! -----------------------------------------------------------------------------------------
kernel_file_name = trim(pseudo_dir)//'/'//vdw_table_name
inquire(file=kernel_file_name, exist=file_exists)
if (.not. file_exists) then
!! Finally, try the default pw_dir/PW/vdW_kernel_table file
!! --------------------------------------------------------------------------------------
kernel_file_name = 'DEFAULT_KERNEL_TABLE_FILE'
inquire(file=kernel_file_name, exist=file_exists)
if (.not. file_exists) then
!! No "vdW_kernel_table" file could be found. Time to die.
call errore('read_kernel_table', 'No \"vdW_kernel_table\" file could be found',1)
end if
end if
end if
!! Generates the md5 file
CALL md5_from_file(kernel_file_name, vdw_kernel_md5_cksum)
!! Open the file to read
open(unit=kernel_file, file=kernel_file_name, status='old', form='formatted', action='read')
!! Read in the number of q points used for this kernel file, the
!! number of r points, and the maximum value of the r point
read(kernel_file, '(2i5)') Nqs, Nr_points
read(kernel_file, double_format) r_max
allocate( q_mesh(Nqs) )
allocate( kernel(0:Nr_points,Nqs,Nqs), d2phi_dk2(0:Nr_points,Nqs,Nqs) )
!! Read in the values of the q points used to generate this kernel
read(kernel_file, double_format) q_mesh
!! For each pair of q values, read in the function phi_q1_q2(k).
!! That is, the fourier transformed kernel function assuming q1 and q2
!! for all the values of r used.
!! ----------------------------------------------------------------------------------------------
do q1_i = 1, Nqs
do q2_i = 1, q1_i
read(kernel_file, double_format) kernel(0:Nr_points, q1_i, q2_i)
kernel(0:Nr_points, q2_i, q1_i) = kernel(0:Nr_points, q1_i, q2_i)
end do
end do
!! ----------------------------------------------------------------------------------------------
!! Again, for each pair of q values (q1 and q2), read in the value
!! of the second derivative of the above mentiond Fourier transformed
!! kernel function phi_alpha_beta(k). These are used for spline
!! interpolation of the Fourier transformed kernel.
!! -----------------------------------------------------------------------------------------------
do q1_i = 1, Nqs
do q2_i = 1, q1_i
read(kernel_file, double_format) d2phi_dk2(0:Nr_points, q1_i, q2_i)
d2phi_dk2(0:Nr_points, q2_i, q1_i) = d2phi_dk2(0:Nr_points, q1_i, q2_i)
end do
end do
!! -----------------------------------------------------------------------------------------------
close(kernel_file)
!! Define a few more vaiables useful to some of the subroutines in xc_vdW_DF.f90
!! ------------------------------------------------------------------------------------------------
q_cut = q_mesh(Nqs)
q_min = q_mesh(1)
dk = 2.0D0*pi/r_max
!! ------------------------------------------------------------------------------------------------
end subroutine initialize_kernel_table
!! #################################################################################################
end MODULE kernel_table
|
