Octopus
ps_xml.F90
Go to the documentation of this file.
1!! Copyright (C) 2015-2018 Xavier Andrade
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module ps_xml_oct_m
22 use debug_oct_m
23 use global_oct_m
24 use, intrinsic :: iso_fortran_env
28 use pseudo_oct_m
29 use string_oct_m
30
31 implicit none
32
33 private
34 public :: &
35 ps_xml_t, &
38
39 type ps_xml_t
40 ! Components are public by default
41 logical :: initialized
42 logical :: kleinman_bylander
43 logical :: nlcc
44 logical :: has_density
45 integer :: atomic_number
46 real(real64) :: mass
47 real(real64) :: valence_charge
48 integer :: lmax
49 integer :: llocal
50 integer :: nchannels
51 integer :: grid_size
52 integer :: nwavefunctions
53 real(real64), allocatable :: grid(:)
54 real(real64), allocatable :: weights(:)
55 real(real64), allocatable :: potential(:, :)
56 real(real64), allocatable :: wavefunction(:, :)
57 real(real64), allocatable :: projector(:, :, :)
58 real(real64), allocatable :: dij(:, :, :)
59 real(real64), allocatable :: nlcc_density(:)
60 real(real64), allocatable :: density(:)
61 integer, allocatable :: wf_n(:)
62 integer, allocatable :: wf_l(:)
63 real(real64), allocatable :: wf_occ(:)
64 type(pseudo_t) :: pseudo
65 end type ps_xml_t
66
67contains
68
69 ! ---------------------------------------------------------
70 subroutine ps_xml_init(this, namespace, filename, fmt, ierr)
71 type(ps_xml_t), intent(inout) :: this
72 type(namespace_t), intent(in) :: namespace
73 character(len=*), intent(in) :: filename
74 integer, intent(in) :: fmt
75 integer, intent(out) :: ierr
76
77 integer :: ll, ii, ic, jc
78 type(pseudo_t) :: pseudo
79
80 push_sub(ps_xml_init)
81
82 this%initialized = .false.
83
84 call pseudo_init(pseudo, string_f_to_c(filename), fmt, ierr)
85
86 if (ierr == pseudo_status_file_not_found) then
87 call messages_write("Pseudopotential file '" // trim(filename) // "' not found")
88 call messages_fatal(namespace=namespace)
89 end if
90
91 if (ierr == pseudo_status_unknown_format) then
92 call messages_write("Cannot determine the format for pseudopotential file '" // trim(filename) // "'")
93 call messages_fatal(namespace=namespace)
94 end if
95
97 call messages_write("Ultrasoft pseudopotential file '" // trim(filename) // "' not supported")
98 call messages_fatal(namespace=namespace)
99 end if
100
101 if (ierr == pseudo_status_unsupported_type_paw) then
102 call messages_write("PAW pseudopotential file '" // trim(filename) // "' not supported")
103 call messages_fatal(namespace=namespace)
104 end if
105
106 if (ierr == pseudo_status_unsupported_type) then
107 call messages_write("Pseudopotential file '" // trim(filename) // "' not supported")
108 call messages_fatal(namespace=namespace)
109 end if
110
111 if (ierr == pseudo_status_format_not_supported) then
112 call messages_write("Pseudopotential file '" // trim(filename) // "' not supported")
113 call messages_fatal(namespace=namespace)
114 end if
115
116 this%initialized = .true.
117 this%valence_charge = pseudo_valence_charge(pseudo)
118 this%mass = pseudo_mass(pseudo)
119 this%lmax = pseudo_lmax(pseudo)
120 this%llocal = pseudo_llocal(pseudo)
121 this%nchannels = pseudo_nchannels(pseudo)
122
123 this%kleinman_bylander = (pseudo_type(pseudo) == pseudo_type_kleinman_bylander)
124
125 this%grid_size = pseudo_mesh_size(pseudo)
126
127 safe_allocate(this%grid(1:this%grid_size))
128 safe_allocate(this%weights(1:this%grid_size))
129
130 call pseudo_grid(pseudo, this%grid(1))
131 call pseudo_grid_weights(pseudo, this%weights(1))
132
133 if (.not. this%kleinman_bylander) then
135 safe_allocate(this%potential(1:this%grid_size, 0:this%lmax))
136 safe_allocate(this%wavefunction(1:this%grid_size, 0:this%lmax))
138 do ll = 0, this%lmax
139 if (.not. pseudo_has_radial_function(pseudo, ll)) then
140 call messages_write("The pseudopotential file '"//trim(filename)//"' does not contain")
142 call messages_write("the wave functions. Octopus cannot use it.")
143 call messages_fatal(namespace=namespace)
144 end if
146 call pseudo_radial_function(pseudo, ll, this%wavefunction(1, ll))
147 call pseudo_radial_potential(pseudo, ll, this%potential(1, ll))
148 end do
150 else
152 safe_allocate(this%potential(1:this%grid_size, this%llocal:this%llocal))
154 call pseudo_local_potential(pseudo, this%potential(1, this%llocal))
156 safe_allocate(this%projector(1:this%grid_size, 0:this%lmax, 1:this%nchannels))
158 this%projector = 0.0_real64
160 do ll = 0, this%lmax
161 if (this%llocal == ll) cycle
162 do ic = 1, pseudo_nprojectors_per_l(pseudo, ll)
163 call pseudo_projector(pseudo, ll, ic, this%projector(1, ll, ic))
164 end do
165 end do
166
167 safe_allocate(this%dij(0:this%lmax, 1:this%nchannels, 1:this%nchannels))
168
169 this%dij = 0.0_real64
170 do ll = 0, this%lmax
171 do ic = 1, this%nchannels
172 do jc = 1, this%nchannels
173 this%dij(ll, ic, jc) = pseudo_dij(pseudo, ll, ic, jc)
174 end do
175 end do
176 end do
177
178 this%nwavefunctions = pseudo_nwavefunctions(pseudo)
179
180 safe_allocate(this%wavefunction(1:this%grid_size, 1:this%nwavefunctions))
181 safe_allocate(this%wf_n(1:this%nwavefunctions))
182 safe_allocate(this%wf_l(1:this%nwavefunctions))
183 safe_allocate(this%wf_occ(1:this%nwavefunctions))
184
185 do ii = 1, this%nwavefunctions
186 call pseudo_wavefunction(pseudo, ii, this%wf_n(ii), this%wf_l(ii), this%wf_occ(ii), this%wavefunction(1, ii))
187 end do
188
189 end if
190
191 this%has_density = pseudo_has_density(pseudo)
192
193 if (this%has_density) then
194 safe_allocate(this%density(1:this%grid_size))
195 call pseudo_density(pseudo, this%density(1))
196 end if
197
198 this%nlcc = pseudo_has_nlcc(pseudo)
199 if (this%nlcc) then
200 safe_allocate(this%nlcc_density(1:this%grid_size))
201 call pseudo_nlcc_density(pseudo, this%nlcc_density(1))
202 end if
203
204 if (.not. this%kleinman_bylander) call ps_xml_check_normalization(this, namespace)
205
206 this%pseudo = pseudo
207
208 pop_sub(ps_xml_init)
209 end subroutine ps_xml_init
210
211 ! ---------------------------------------------------------
213 subroutine ps_xml_check_normalization(this, namespace)
214 type(ps_xml_t), intent(in) :: this
215 type(namespace_t), intent(in) :: namespace
216
217 integer :: ll, ip
218 real(real64) :: nrm, rr
219
221
222 ! checking normalization of the wavefunctions
223 do ll = 0, this%lmax
224 nrm = m_zero
225 do ip = 1, this%grid_size
226 rr = this%grid(ip)
227 nrm = nrm + this%wavefunction(ip, ll)**2*this%weights(ip)*rr**2
228 end do
229 nrm = sqrt(nrm)
230
231 nrm = abs(nrm - m_one)
232 if (nrm > 1.0e-5_real64) then
233 write(message(1), '(a,i2,a)') "Eigenstate for l = ", ll, ' is not normalized'
234 write(message(2), '(a, f12.6,a)') '(abs(1 - norm) = ', nrm, ')'
235 call messages_warning(2, namespace=namespace)
236 end if
237
238 end do
239
241 end subroutine ps_xml_check_normalization
242
243 ! ---------------------------------------------------------
244 subroutine ps_xml_end(this)
245 type(ps_xml_t), intent(inout) :: this
246
247 push_sub(ps_xml_end)
248
249 call pseudo_end(this%pseudo)
250
251 safe_deallocate_a(this%grid)
252 safe_deallocate_a(this%weights)
253 safe_deallocate_a(this%potential)
254 safe_deallocate_a(this%wavefunction)
255 safe_deallocate_a(this%projector)
256 safe_deallocate_a(this%dij)
257 safe_deallocate_a(this%nlcc_density)
258 safe_deallocate_a(this%density)
259 safe_deallocate_a(this%wf_n)
260 safe_deallocate_a(this%wf_l)
261 safe_deallocate_a(this%wf_occ)
262
263 pop_sub(ps_xml_end)
264 end subroutine ps_xml_end
265
266end module ps_xml_oct_m
267
268!! Local Variables:
269!! mode: f90
270!! coding: utf-8
271!! End:
double sqrt(double __x) __attribute__((__nothrow__
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_one
Definition: global.F90:192
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_new_line()
Definition: messages.F90:1112
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public ps_xml_end(this)
Definition: ps_xml.F90:340
subroutine ps_xml_check_normalization(this, namespace)
checks normalization of the pseudo wavefunctions
Definition: ps_xml.F90:309
subroutine, public ps_xml_init(this, namespace, filename, fmt, ierr)
Definition: ps_xml.F90:166
integer, parameter, public pseudo_type_kleinman_bylander
Definition: pseudo.F90:160
logical function, public pseudo_has_radial_function(pseudo, l)
Definition: pseudo.F90:570
logical function, public pseudo_has_density(pseudo)
Definition: pseudo.F90:533
logical function, public pseudo_has_nlcc(pseudo)
Definition: pseudo.F90:515
integer, parameter, public pseudo_status_unsupported_type
Definition: pseudo.F90:166
integer, parameter, public pseudo_status_file_not_found
Definition: pseudo.F90:166
integer, parameter, public pseudo_status_unsupported_type_ultrasoft
Definition: pseudo.F90:166
integer, parameter, public pseudo_status_unknown_format
Definition: pseudo.F90:166
integer, parameter, public pseudo_status_unsupported_type_paw
Definition: pseudo.F90:166
integer, parameter, public pseudo_status_format_not_supported
Definition: pseudo.F90:166
character(kind=c_char, len=1) function, dimension(:), allocatable, public string_f_to_c(f_string)
convert a Fortran string to a C string
Definition: string.F90:273
int true(void)