Octopus
species_factory.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-2024 N. 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#include "global.h"
20
23 use debug_oct_m
25 use global_oct_m
26 use iihash_oct_m
27 use io_oct_m
31 use parser_oct_m
37 use unit_oct_m
39
40 implicit none
41
42 private
43 public :: &
47
48 type :: species_factory_t
49 logical :: initialized = .false.
50 integer :: default_pseudopotential_set_id
51 type(pseudo_set_t) :: default_pseudopotential_set
52 integer :: default_allelectron_type
53 real(real64) :: default_sigma
54 real(real64) :: default_anc_a
55 contains
56 procedure :: init => species_factory_init
57 procedure :: end => species_factory_end
58 procedure :: create_from_input => species_factory_create_from_input
59 procedure :: create_from_block => read_from_block
60 end type species_factory_t
61
62contains
63
64 ! ---------------------------------------------------------
65 subroutine species_factory_init(factory, namespace)
66 class(species_factory_t), intent(inout) :: factory
67 type(namespace_t), intent(in) :: namespace
68
69 integer :: ierr, default_val
70
71 if (factory%initialized) return
72
73 push_sub_with_profile(species_factory_init)
74
75 factory%initialized = .true.
76
77 call share_directory_set(conf%share)
78
79 !%Variable AllElectronType
80 !%Type integer
81 !%Default no
82 !%Section System::Species
83 !%Description
84 !% Selects the type of all-electron species that applies by default to all
85 !% atoms. This is not compatible with <tt>PseudopotentialSet</tt>, but it is
86 !% compatible with the <tt>Species</tt> block.
87 !%
88 !%Option no 0
89 !% Do not specify any default all-electron type of species. All species must be
90 !% specified in the Species block.
91 !%Option full_delta 1
92 !% All atoms are supposed to be by default of type <tt>species_full_delta</tt>.
93 !%Option full_gaussian 2
94 !% All atoms are supposed to be by default of type <tt>species_full_gaussian</tt>.
95 !%Option full_anc 3
96 !% All atoms are supposed to be by default of type <tt>species_full_anc</tt>.
97 !%End
98 call parse_variable(namespace, 'AllElectronType', option__allelectrontype__no, factory%default_allelectron_type)
99 call messages_print_var_option('AllElectronType', factory%default_allelectron_type, namespace=namespace)
100
101 !%Variable AllElectronSigma
102 !%Type integer
103 !%Default 0.6
104 !%Section System::Species
105 !%Description
106 !% Default value for the parameter <tt>gaussian_width</tt>. This is useful
107 !% for specifying multiple atoms without specifying the species block. The
108 !% default value is taken from the recommendation in
109 !% <i>Phys. Rev. B</i> <b>55</b>, 10289 (1997).
110 !%
111 !%End
112 call parse_variable(namespace, 'AllElectronSigma', 0.6_real64, factory%default_sigma)
113
114 !%Variable AllElectronANCParam
115 !%Type integer
116 !%Default 4
117 !%Section System::Species
118 !%Description
119 !% Default values for the parameter <tt>anc_a</tt>. This is usefull
120 !% for specifying multiple atoms without specifying the species block.
121 !%
122 !%End
123 call parse_variable(namespace, 'AllElectronANCParam', 4.0_real64, factory%default_anc_a)
124
125
126 !%Variable PseudopotentialSet
127 !%Type integer
128 !%Default standard
129 !%Section System::Species
130 !%Description
131 !% Selects the set of pseudopotentials used by default for species
132 !% not defined in the <tt>Species</tt> block.
133 !%
134 !% These sets of pseudopotentials come from different
135 !% sources. Octopus developers have not validated them. We include
136 !% them with the code for convenience of the users, but you are
137 !% expected to check the quality and suitability of the
138 !% pseudopotential for your application.
139 !%
140 !% Note that not all these pseudopotentials are compatible with spin-orbit coupling
141 !% Only hgh_lda, hgh_lda_sc and pseudodojo_pbe_fr provide SOC information.
142 !%
143 !%Option none 0
144 !% Do not load any pseudopotential by default. All species must be
145 !% specified in the Species block.
146 !%Option standard 1
147 !% The standard set of Octopus that provides LDA pseudopotentials
148 !% in the PSF format for some elements: H, Li, C, N, O, Na, Si, S, Ti, Se, Cd.
149 !%Option sg15 2
150 !% The set of Optimized Norm-Conserving Vanderbilt
151 !% PBE pseudopotentials. Ref: M. Schlipf and F. Gygi, <i>Comp. Phys. Commun.</i> <b>196</b>, 36 (2015).
152 !% This set provides pseudopotentials for elements up to Z = 83
153 !% (Bi), excluding Lanthanides.
154 !% Current version of the set is 1.2.
155 !%Option hgh_lda 3
156 !% The set of Hartwigsen-Goedecker-Hutter LDA pseudopotentials for elements from H to Rn.
157 !% Ref: C. Hartwigsen, S. Goedecker, and J. Hutter, <i>Phys. Rev. B</i> <b>58</b>, 3641 (1998).
158 !%Option hgh_lda_sc 31
159 !% The semicore set of Hartwigsen-Goedecker-Hutter LDA pseudopotentials.
160 !% Ref: C. Hartwigsen, S. Goedecker, and J. Hutter, <i>Phys. Rev. B</i> <b>58</b>, 3641 (1998).
161 !%Option hscv_lda 4
162 !% The set of Hamann-Schlueter-Chiang-Vanderbilt (HSCV) potentials
163 !% for LDA exchange and correlation downloaded from http://fpmd.ucdavis.edu/potentials/index.htm.
164 !% These pseudopotentials were originally intended for the QBox
165 !% code. They were generated using the method of Hamann, Schluter
166 !% and Chiang. Ref: D. Vanderbilt, <i>Phys. Rev. B</i> <b>32</b>, 8412 (1985).
167 !% Warning from the original site: The potentials provided in this
168 !% site are distributed without warranty. In most cases,
169 !% potentials were not tested. Potentials should be thoroughly
170 !% tested before being used in simulations.
171 !%Option hscv_pbe 5
172 !% PBE version of the HSCV pseudopotentials. Check the
173 !% documentation of the option <tt>hscv_lda</tt> for details and warnings.
174 !%Option pseudodojo_pbe 100
175 !% PBE version of the pseudopotentials of http://pseudo-dojo.org. Version 0.5
176 !%Option pseudodojo_lda 103
177 !% LDA pseudopotentials of http://pseudo-dojo.org. Version 0.41.
178 !%Option pseudodojo_pbesol 105
179 !% PBEsol version of the pseudopotentials of http://pseudo-dojo.org. Version 0.41.
180 !%Option pseudodojo_pbe_fr 106
181 !% Fully-relativistic PBE version of the pseudopotentials of http://pseudo-dojo.org. Version 0.4.
182 !%End
183 default_val = option__pseudopotentialset__standard
184 if(factory%default_allelectron_type /= option__allelectrontype__no) default_val = option__pseudopotentialset__none
185 call parse_variable(namespace, 'PseudopotentialSet', default_val, factory%default_pseudopotential_set_id)
186 call messages_print_var_option('PseudopotentialSet', factory%default_pseudopotential_set_id, namespace=namespace)
187 if (factory%default_pseudopotential_set_id /= option__pseudopotentialset__none) then
188 call pseudo_set_init(factory%default_pseudopotential_set, get_set_directory(factory%default_pseudopotential_set_id), ierr)
189 else
190 call pseudo_set_nullify(factory%default_pseudopotential_set)
191 end if
192
193 if (factory%default_pseudopotential_set_id /= option__pseudopotentialset__none &
194 .and. factory%default_allelectron_type /= option__allelectrontype__no) then
195 message(1) = "PseudopotentialSet /= none cannot be used with AllElectronType /= no."
196 call messages_fatal(1, namespace=namespace)
197 end if
198
199 pop_sub_with_profile(species_factory_init)
200 end subroutine species_factory_init
201
202 ! ---------------------------------------------------------
203
204 subroutine species_factory_end(factory)
205 class(species_factory_t), intent(inout) :: factory
206
207 push_sub(species_factory_end)
208
209 if (factory%initialized) then
210 call pseudo_set_end(factory%default_pseudopotential_set)
211 factory%initialized = .false.
212 end if
213
214 pop_sub(species_factory_end)
215 end subroutine species_factory_end
216
217
218 ! ---------------------------------------------------------
223 ! ---------------------------------------------------------
224 function species_factory_create_from_input(factory, namespace, label, index) result(spec)
225 class(species_factory_t), intent(in) :: factory
226 type(namespace_t), intent(in) :: namespace
227 character(len=*), intent(in) :: label
228 integer, intent(in) :: index
229 class(species_t), pointer :: spec
230
231 character(len=LABEL_LEN) :: lab
232 integer :: ib, row, n_spec_block, read_data
233 type(block_t) :: blk
234 type(element_t) :: element
235
237
238 read_data = 0
239
240 !%Variable Species
241 !%Type block
242 !%Section System::Species
243 !%Description
244 !% A species is by definition either an "ion" (nucleus + core electrons) described
245 !% through a pseudopotential, or a model potential.
246 !%
247 !% Note that some sets of pseudopotentials are distributed with
248 !% the code. To use these pseudopotentials, you do not need to define them
249 !% explicitly in the <tt>Species</tt> block, as default parameters
250 !% are provided.
251 !% You can select the set for default pseudopotentials using the
252 !% <tt>PseudopotentialSet</tt> variable.
253 !%
254 !% Supported norm-conserving pseudopotential formats are
255 !% detected by the file extension: UPF (<tt>.upf</tt>), PSF (SIESTA, <tt>.psf</tt>), FHI (ABINIT 6, <tt>.fhi</tt>),
256 !% CPI (Fritz-Haber, <tt>.cpi</tt>), QSO (quantum-simulation.org, for Qbox, <tt>.xml</tt>),
257 !% HGH (Hartwigsen-Goedecker-Hutter, <tt>.hgh</tt>).
258 !% PSPIO format can also be used via <tt>species_pspio</tt> if that library is linked.
259 !% Note: pseudopotentials may only be used in 3D.
260 !%
261 !% The format of this block is the following: The first field is a
262 !% string that defines the name of the species. The second field
263 !% defines the type of species (the valid options are detailed
264 !% below).
265 !%
266 !% Then a list of parameters follows. The parameters are specified
267 !% by a first field with the parameter name and the field that
268 !% follows with the value of the parameter. Some parameters are
269 !% specific to a certain species while others are accepted by all
270 !% species. These are <tt>mass</tt>, <tt>max_spacing</tt>, and <tt>min_radius</tt>.
271 !%
272 !% These are examples of possible species:
273 !%
274 !% <tt>%Species
275 !% <br>&nbsp;&nbsp;'O' | species_pseudo | file | 'O.psf' | lmax | 1 | lloc | 1
276 !% <br>&nbsp;&nbsp;'H' | species_pseudo | file | '../H.hgh'
277 !% <br>&nbsp;&nbsp;'Xe' | species_pseudo | set | pseudojo_pbe
278 !% <br>&nbsp;&nbsp;'C' | species_pseudo | file | "carbon.xml"
279 !% <br>&nbsp;&nbsp;'jlm' | species_jellium | jellium_radius | 5.0
280 !% <br>&nbsp;&nbsp;'rho' | species_charge_density | density_formula | "exp(-r/a)" | mass | 17.0 | valence | 6
281 !% <br>&nbsp;&nbsp;'udf' | species_user_defined | potential_formula | "1/2*r^2" | valence | 8
282 !% <br>&nbsp;&nbsp;'He_all' | species_full_delta
283 !% <br>&nbsp;&nbsp;'H_all' | species_full_gaussian | gaussian_width | 0.2
284 !% <br>&nbsp;&nbsp;'Li1D' | species_soft_coulomb | softening | 1.5 | valence | 3
285 !% <br>&nbsp;&nbsp;'H_all' | species_full_anc | anc_a | 4
286 !% <br>%</tt>
287 !%Option species_pseudo -7
288 !% The species is a pseudopotential. How to get the
289 !% pseudopotential can be specified by the <tt>file</tt> or
290 !% the <tt>set</tt> parameters. If both are missing, the
291 !% pseudopotential will be taken from the <tt>PseudopotentialSet</tt>
292 !% specified for the run, this is useful if you want to change
293 !% some parameters of the pseudo, like the <tt>mass</tt>.
294 !%
295 !% The optional parameters for this type of species are
296 !% <tt>lmax</tt>, that defines the maximum angular momentum
297 !% component to be used, and <tt>lloc</tt>, that defines the
298 !% angular momentum to be considered as local. When these
299 !% parameters are not set, the value for lmax is the maximum
300 !% angular component from the pseudopotential file. The default
301 !% value for <tt>lloc</tt> is taken from the pseudopotential if
302 !% available, if not, it is set to 0. Note that, depending on the
303 !% type of pseudopotential, it might not be possible to select
304 !% <tt>lmax</tt> and <tt>lloc</tt>, if that is the case the
305 !% parameters will be ignored.
306 !%
307 !%Option species_pspio -110
308 !% (experimental) Alternative method to read pseudopotentials
309 !% using the PSPIO library. This species uses the same parameters
310 !% as <tt>species_pseudo</tt>.
311 !%Option species_user_defined -123
312 !% Species with user-defined potential. The potential for the
313 !% species is defined by the formula given by the <tt>potential_formula</tt>
314 !% parameter.
315 !% The
316 !% <tt>valence</tt> parameter determines the number of electrons
317 !% associated with the species. By default, a valence of 0 is assumed.
318 !%Option species_charge_density -125
319 !% The potential for this species is created from the distribution
320 !% of charge given by the <tt>density_formula</tt> parameter.
321 !% The
322 !% <tt>valence</tt> parameter determines the number of electrons
323 !% associated with the species. By default, a valence of 0 is assumed.
324 !%Option species_point -3
325 !%Option species_jellium -3
326 !% Jellium sphere.
327 !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
328 !%Option species_jellium_slab -4
329 !% A slab of jellium that extends across the simulation box in the
330 !% <i>xy</i>-plane. The dimension along the <i>z</i> direction is
331 !% determined by the required parameter <tt>thickness</tt>.
332 !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
333 !%Option species_full_delta -127
334 !% Full atomic potential represented by a delta charge
335 !% distribution. The atom will be displaced to the nearest grid
336 !% point. The atomic number is determined from the name of the species.
337 !%Option species_full_anc -130
338 !% Analytical norm-conserving regulized Coulomb potential from
339 !% [Gygi J. Chem. Theory Comput. 2023, 19, 1300−1309].
340 !%Option species_full_gaussian -124
341 !% A full-potential atom is defined by a Gaussian accumulation of
342 !% positive charge (distorted if curvilinear coordinates are
343 !% used), in the form:
344 !%
345 !% <math>q(r) = z \beta \exp[ - (\vec{r}-\vec{r_0})^2 / (\sqrt{2} \delta \sigma) ] </math>
346 !%
347 !% <math>\beta</math> is chosen in order to maintain proper
348 !% normalization (the integral of <math>q</math> should sum up to
349 !% <math>z</math>). <math>\delta</math> is the grid spacing (the
350 !% grid spacing in the first dimension, to be precise).
351 !% <math>\vec{r_0}</math> is calculated in such a way that the the
352 !% first moment of <math>q(r)/z</math> is equal to the atomic
353 !% position. For a precise description, see N. A. Modine,
354 !% <i>Phys. Rev. B</i> <b>55</b>, 10289 (1997). The width of the
355 !% Gaussian is set by parameter <tt>gaussian_width</tt>. The
356 !% atomic number is determined from the name of the species.
357 !%Option species_from_file -126
358 !% The potential is read from a file. Accepted file formats, detected by extension: obf, ncdf and csv.
359 !% The
360 !% <tt>valence</tt> parameter determines the number of electrons
361 !% associated with the species. By default, a valence of 0 is assumed.
362 !%Option species_soft_coulomb -128
363 !% The potential is a soft-Coulomb function, <i>i.e.</i> a function in the form:
364 !%
365 !% <math>v(r) = - z_{val} / \sqrt{a^2 + r^2}</math>
366 !%
367 !% The value of <i>a</i> should be given by the mandatory <tt>softening</tt> parameter.
368 !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
369 !%Option species_jellium_charge_density -129
370 !% The parameter is the name of a volume block specifying the shape of the jellium.
371 !%Option lmax -10003
372 !% The maximum angular-momentum channel that will be used for the pseudopotential.
373 !%Option lloc -10004
374 !% The angular-momentum channel of the pseudopotential to be considered local.
375 !%Option mass -10005
376 !% The mass of the species in atomic mass units, <i>i.e.</i> the mass of a proton is
377 !% roughly one. It is set automatically for pseudopotentials from the
378 !% <a href=http://www.nist.gov/pml/data/comp.cfm>NIST values</a>.
379 !% For other species, the default is 1.0.
380 !%Option valence -10006
381 !% The number of electrons of the species. It is set automatically from the name of the species.
382 !% if it correspond to the name in the periodic table. If not specified and if the name
383 !% does not match an atom name, a value of 0 is assumed.
384 !%Option jellium_radius -10007
385 !% The radius of the sphere for <tt>species_jellium</tt>. If this value is not specified,
386 !% the default of 0.5 bohr is used.
387 !%Option set -10017
388 !% For a <tt>species_pseudo</tt>, get the pseudopotential from a
389 !% particular set. This flag must be followed with one of the
390 !% valid values for the variable <tt>PseudopotentialSet</tt>.
391 !%Option gaussian_width -10008
392 !% The width of the Gaussian (in units of spacing) used to represent
393 !% the nuclear charge for <tt>species_full_gaussian</tt>. If not present,
394 !% the default is 0.6.
395 !%Option softening -10009
396 !% The softening parameter <i>a</i> for <tt>species_soft_coulomb</tt> in units of length.
397 !%Option file -10010
398 !% The path for the file that describes the species.
399 !%Option db_file -10011
400 !% Obsolete. Use the <tt>set</tt> option of the <tt>PseudopotentialSet</tt> variable instead.
401 !%Option potential_formula -10012
402 !% Mathematical expression that defines the potential for <tt>species_user_defined</tt>. You can use
403 !% any of the <i>x</i>, <i>y</i>, <i>z</i> or <i>r</i> variables.
404 !%Option density_formula -10013
405 !% Mathematical expression that defines the charge density for <tt>species_charge_density</tt>. You can use
406 !% any of the <i>x</i>, <i>y</i>, <i>z</i> or <i>r</i> variables.
407 !%Option thickness -10014
408 !% The thickness of the slab for species_jellium_slab. Must be positive.
409 !%Option vdw_radius -10015
410 !% The van der Waals radius that will be used for this species.
411 !%Option volume -10016
412 !% Name of a volume block
413 !%Option hubbard_l -10018
414 !% The angular-momentum for which the effective U will be applied.
415 !%Option hubbard_u -10019
416 !% The effective U that will be used for the DFT+U calculations.
417 !%Option hubbard_j -10020
418 !% The value of j (hubbard_l-1/2 or hubbard_l+1/2) on which the effective U is applied.
419 !%Option hubbard_alpha -10021
420 !% The strength of the potential constraining the occupations of the localized subspace
421 !% as defined in PRB 71, 035105 (2005)
422 !%Option anc_a -10022
423 !% The value of the parameter a of the ANC potential, as defined in [Gygi, JCTC 2023, 19, 1300−1309].
424 !% This parameter has the unit of inverse length and determines the range of regularization.
425 !%End
426
427 call messages_obsolete_variable(namespace, 'SpecieAllElectronSigma', 'Species')
428 call messages_obsolete_variable(namespace, 'SpeciesAllElectronSigma', 'Species')
429
430 ! First, find out if there is a Species block.
431 n_spec_block = 0
432 if (parse_block(namespace, 'Species', blk) == 0) then
433 n_spec_block = parse_block_n(blk)
434 end if
435
436 ! Find out if the sought species is in the block
437 row = -1
438 block: do ib = 1, n_spec_block
439 call parse_block_string(blk, ib-1, 0, lab)
440 if (trim(lab) == trim(label)) then
441 row = ib - 1
442 exit block
443 end if
444 end do block
445
446 ! Read whatever may be read from the block
447 if (row >= 0) then
448 spec => factory%create_from_block(namespace, blk, row, label, index, read_data)
449 call parse_block_end(blk)
450
451 assert(read_data > 0)
452
454 return
455 end if
456
457 ! We get here if there is a Species block but it does not contain
458 ! the species we are looking for.
459 if (n_spec_block > 0) call parse_block_end(blk)
460
461 ! Initialize all electron species (except soft-Coulomb) from specified allelectron type
462 if(factory%default_allelectron_type /= 0) then
463 select case(factory%default_allelectron_type)
464 case(option__allelectrontype__full_delta)
465 spec => full_delta_t(label, index, factory%default_sigma)
466 case(option__allelectrontype__full_gaussian)
467 spec => full_gaussian_t(label, index, factory%default_sigma)
468 case(option__allelectrontype__full_anc)
469 spec => full_anc_t(label, index, factory%default_anc_a)
470 case default
471 assert(.false.)
472 end select
473
474 ! get the mass, vdw radius and atomic number for this element
475 call element_init(element, get_symbol(spec%get_label()))
476
477 assert(element_valid(element))
478
479 call spec%set_z(real(element_atomic_number(element), real64))
480 call spec%set_zval(spec%get_z())
481 call spec%set_mass(element_mass(element))
482 call spec%set_vdw_radius(element_vdw_radius(element))
483
484 call element_end(element)
485
486 else ! Pseudopotential from specified set or the default one
487 spec => pseudopotential_t(label, index)
488 select type(spec)
489 type is(pseudopotential_t)
490 call read_from_set(spec, factory%default_pseudopotential_set_id, factory%default_pseudopotential_set, read_data)
491
492 if (read_data == 0) then
493 call messages_write( 'Species '//trim(spec%get_label())//' not found in default pseudopotential set.', new_line=.true. )
494 call messages_write('( '//trim(get_set_directory(factory%default_pseudopotential_set_id))//' )')
495 call messages_fatal(namespace=namespace)
496 end if
497 end select
498 end if
499
502
503!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
504! Private procedures
505
506 ! ---------------------------------------------------------
508 function read_from_block(factory, namespace, blk, row, label, index, read_data) result(spec)
509 class(species_factory_t), intent(in) :: factory
510 type(namespace_t), intent(in) :: namespace
511 type(block_t), intent(in) :: blk
512 integer, intent(in) :: row
513 character(len=*), intent(in) :: label
514 integer, intent(in) :: index
515 integer, intent(out):: read_data
516 class(species_t), pointer :: spec
517
518 integer :: ncols, icol, flag, set_read_data, ierr, type
519 type(element_t) :: element
520 type(iihash_t) :: read_parameters
521 integer :: user_lmax, user_llocal, hubbard_l, pseudopotential_set_id
522 real(real64) :: hubbard_u, hubbard_j, hubbard_alpha, mass, z_val, jradius, jthick, vdw_radius, aa
523 real(real64) :: sigma, softening
524 character(len=MAX_PATH_LEN) :: filename
525
526 integer, parameter :: &
527 species_jellium = 3, & !< jellium sphere.
530 species_pseudo = 7, &
531 species_pspio = 110, &
532 species_usdef = 123, &
533 species_full_gaussian = 124, &
534 species_charge_density = 125, &
535 species_from_file = 126, &
536 species_full_delta = 127, &
537 species_soft_coulomb = 128, &
538 species_full_anc = 130
539
540 push_sub(read_from_block)
541
542 ncols = parse_block_cols(blk, row)
543 read_data = 0
544
545 call parse_block_integer(blk, row, 1, type)
546
547 ! To detect the old species block format, options are represented
548 ! as negative values. If we get a non-negative value we know we
549 ! are reading a mass.
550 if (type >= 0) then
551 call messages_write('Found a species with the old format. Please update', new_line = .true.)
552 call messages_write('the Species block to the new format, where the second', new_line = .true.)
553 call messages_write('column indicates the type of the species.')
554 call messages_fatal(namespace=namespace)
555 end if
556
557 ! now we convert back to positive
558 type = -type
559
560 read_data = 2
561
562
563 select case (type)
564
565 case (species_soft_coulomb)
566 spec => soft_coulomb_t(label, index)
567
568 case (species_usdef) ! user-defined
569 spec => species_user_defined_t(label, index)
570
571 case (species_from_file)
572 spec => species_from_file_t(label, index)
573
574 case (species_jellium)
575 spec => jellium_sphere_t(label, index)
576
578 spec => jellium_slab_t(label, index)
579
580 case (species_full_delta)
581 spec => full_delta_t(label, index, factory%default_sigma)
582
583 case (species_full_gaussian)
584 spec => full_gaussian_t(label, index, factory%default_sigma)
585
586 case (species_full_anc)
587 spec => full_anc_t(label, index, factory%default_anc_a)
588
590 spec => species_charge_density_t(label, index)
591
593 spec => jellium_charge_t(label, index)
594
595 case (species_pseudo)
596 spec => pseudopotential_t(label, index)
597
598 case (species_pspio) ! a pseudopotential file to be handled by the pspio library
599 spec => pseudopotential_t(label, index)
600
601 case default
602 call messages_input_error(namespace, 'Species', "Unknown type for species '"//trim(spec%get_label())//"'", row=row, column=1)
603 end select
604
605 call spec%set_mass(-m_one)
606 call spec%set_vdw_radius(-m_one)
607 call spec%set_zval(-m_one)
608
609 call iihash_init(read_parameters)
610
611 icol = read_data
612 do
613 if (icol >= ncols) exit
614
615 call parse_block_integer(blk, row, icol, flag)
616
617 select case (flag)
618
619 case (option__species__lmax)
620 call check_duplication(option__species__lmax)
621 call parse_block_integer(blk, row, icol + 1, user_lmax)
622
623 select type(spec)
624 class is(pseudopotential_t)
625 call spec%set_user_lmax(user_lmax)
626 class default
627 call messages_input_error(namespace, 'Species', &
628 "The 'lmax' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
629 row=row, column=icol+1)
630 end select
631
632 if (user_lmax < 0) then
633 call messages_input_error(namespace, 'Species', &
634 "The 'lmax' parameter in species "//trim(spec%get_label())//" cannot be negative", &
635 row=row, column=icol+1)
636 end if
637
638 case (option__species__lloc)
639 call check_duplication(option__species__lloc)
640 call parse_block_integer(blk, row, icol + 1, user_llocal)
641
642 select type(spec)
643 class is(pseudopotential_t)
644 call spec%set_user_lloc(user_llocal)
645 class default
646 call messages_input_error(namespace, 'Species', &
647 "The 'lloc' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
648 row=row, column=icol+1)
649 end select
650
651 if (user_llocal < 0) then
652 call messages_input_error(namespace, 'Species', &
653 "The 'lloc' parameter in species "//trim(spec%get_label())//" cannot be negative", row=row, column=icol+1)
654 end if
655
656 case (option__species__hubbard_l)
657 call check_duplication(option__species__hubbard_l)
658 call parse_block_integer(blk, row, icol + 1, hubbard_l)
659
660 select type(spec)
661 class is(pseudopotential_t)
662 call spec%set_hubbard_l(hubbard_l)
663 class default
664 call messages_input_error(namespace, 'Species', &
665 "The 'hubbard_l' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
666 row=row, column=icol+1)
667 end select
668
669 if (hubbard_l < 0) then
670 call messages_input_error(namespace, 'Species', &
671 "The 'hubbard_l' parameter in species "//trim(spec%get_label())//" cannot be negative", row=row, column=icol+1)
672 end if
673
674 case (option__species__hubbard_u)
675 call check_duplication(option__species__hubbard_u)
676 call parse_block_float(blk, row, icol + 1, hubbard_u, unit = units_inp%energy)
677 call spec%set_hubbard_u(hubbard_u)
678
679 case (option__species__hubbard_alpha)
680 call check_duplication(option__species__hubbard_alpha)
681 call parse_block_float(blk, row, icol + 1, hubbard_alpha, unit = units_inp%energy)
682 call spec%set_hubbard_u(hubbard_alpha)
683
684 case (option__species__hubbard_j)
685 call check_duplication(option__species__hubbard_j)
686 call parse_block_float(blk, row, icol + 1, hubbard_j)
687 call spec%set_hubbard_u(hubbard_alpha)
688
689 if (abs(abs(spec%get_hubbard_j()-spec%get_hubbard_l())-m_half) <= m_epsilon) then
690 call messages_input_error(namespace, 'Species', "The 'hubbard_j' parameter in species "// &
691 trim(spec%get_label())//" can only be hubbard_l +/- 1/2", row=row, column=icol+1)
692 end if
693
694 case (option__species__mass)
695 call check_duplication(option__species__mass)
696 call parse_block_float(blk, row, icol + 1, mass, unit = units_inp%mass)
697 call spec%set_mass(mass)
698
699 case (option__species__valence)
700 call check_duplication(option__species__valence)
701 call parse_block_float(blk, row, icol + 1, z_val)
702 call spec%set_zval(z_val)
703 call spec%set_z(z_val)
704
705 case (option__species__jellium_radius)
706 call check_duplication(option__species__jellium_radius)
707 call parse_block_float(blk, row, icol + 1, jradius)
708 if (jradius <= m_zero) call messages_input_error(namespace, 'Species', 'jellium_radius must be positive', &
709 row=row, column=icol+1)
710
711 select type(spec)
712 type is(jellium_sphere_t)
713 call spec%set_radius(jradius)
714 class default
715 call messages_input_error(namespace, 'Species', 'jellium_radius can only be used with species_jellium', &
716 row=row, column=icol+1)
717 end select
718
719 case (option__species__gaussian_width)
720 call check_duplication(option__species__gaussian_width)
721 call parse_block_float(blk, row, icol + 1, sigma)
722 if (sigma <= m_zero) call messages_input_error(namespace, 'Species', 'gaussian_width must be positive', &
723 row=row, column=icol+1)
724 select type(spec)
725 type is(full_gaussian_t)
726 call spec%set_sigma(sigma)
727 class default
728 call messages_input_error(namespace, 'Species', 'gaussian_width can only be used with species_full_gaussian', &
729 row=row, column=icol+1)
730 end select
731
732 case (option__species__anc_a)
733 call check_duplication(option__species__anc_a)
734 call parse_block_float(blk, row, icol + 1, aa)
735 if (aa <= m_zero) call messages_input_error(namespace, 'Species', 'anc_a must be positive', &
736 row=row, column=icol+1)
737
738 select type(spec)
739 type is(full_anc_t)
740 call spec%set_a(aa)
741 class default
742 call messages_input_error(namespace, 'Species', 'anc_a can only be used with species_full_anc', &
743 row=row, column=icol+1)
744 end select
745
746 case (option__species__softening)
747 call check_duplication(option__species__softening)
748 call parse_block_float(blk, row, icol + 1, softening)
749 softening = softening**2
750
751 select type(spec)
752 type is(soft_coulomb_t)
753 call spec%set_softening2(softening)
754 class default
755 call messages_input_error(namespace, 'Species', 'softening can only be used with species_soft_coulomb', &
756 row=row, column=icol+1)
757 end select
758
759 case (option__species__file)
760 call check_duplication(option__species__file)
761 call parse_block_string(blk, row, icol + 1, filename)
762 call spec%set_filename(filename)
763
764 case (option__species__db_file)
765 call messages_write("The 'db_file' option for 'Species' block is obsolete. Please use", new_line = .true.)
766 call messages_write("the option 'set' or the variable 'PseudopotentialSet' instead.")
767 call messages_fatal(namespace=namespace)
768
769 case (option__species__set)
770 call check_duplication(option__species__set)
771 call parse_block_integer(blk, row, icol + 1, pseudopotential_set_id)
772
773 select type(spec)
774 type is(pseudopotential_t)
775 spec%pseudopotential_set_initialized = .true.
776 spec%pseudopotential_set_id = pseudopotential_set_id
777 call pseudo_set_init(spec%pseudopotential_set, get_set_directory(spec%pseudopotential_set_id), ierr)
778 class default
779 call messages_input_error(namespace, 'Species', 'set can only be used with species_pseudo', &
780 row=row, column=icol+1)
781 end select
782
783 case (option__species__potential_formula)
784 call check_duplication(option__species__potential_formula)
785 select type(spec)
787 call parse_block_string(blk, row, icol + 1, spec%potential_formula, convert_to_c=.true.)
788 class default
789 call messages_input_error(namespace, 'Species', 'potential_formula can only be used with species_user_defined', &
790 row=row, column=icol+1)
791 end select
792
793 case (option__species__volume)
794 call check_duplication(option__species__volume)
795
796 select type(spec)
797 type is(jellium_charge_t)
798 call parse_block_string(blk, row, icol + 1, spec%density_formula, convert_to_c=.true.)
799
800 class default
801 call messages_input_error(namespace, 'Species', 'volume can only be used with species_jellium_charge_density', &
802 row=row, column=icol+1)
803 end select
804
805 case (option__species__density_formula)
806 call check_duplication(option__species__density_formula)
807
808 select type(spec)
810 call parse_block_string(blk, row, icol + 1, spec%density_formula, convert_to_c=.true.)
811
812 class default
813 call messages_input_error(namespace, 'Species', 'density_formula can only be used with species_charge_density', &
814 row=row, column=icol+1)
815 end select
816
817 case (option__species__thickness)
818 call check_duplication(option__species__thickness)
819 call parse_block_float(blk, row, icol + 1, jthick) ! thickness of the jellium slab
820
821 if (jthick <= m_zero) then
822 call messages_input_error(namespace, 'Species', 'the value of the thickness parameter in species '&
823 //trim(spec%get_label())//' must be positive.', row=row, column=icol+1)
824 end if
825
826 select type(spec)
827 type is(jellium_slab_t)
828 call spec%set_thickness(jthick)
829 class default
830 call messages_input_error(namespace, 'Species', 'thickness can only be used with species_jellium_slab', &
831 row=row, column=icol+1)
832 end select
833
834 case (option__species__vdw_radius)
835 call check_duplication(option__species__vdw_radius)
836 call parse_block_float(blk, row, icol + 1, vdw_radius, unit = units_inp%length)
837 call spec%set_vdw_radius(vdw_radius)
838
839 case default
840 call messages_input_error(namespace, 'Species', "Unknown parameter in species '"//trim(spec%get_label())//"'", &
841 row=row, column=icol)
842
843 end select
844
845 icol = icol + 2
846 end do
847 ! CHECK THAT WHAT WE PARSED MAKES SENSE
848
849 select type(spec)
850 type is(soft_coulomb_t)
851 if (.not. parameter_defined(option__species__softening)) then
852 call messages_input_error(namespace, 'Species', &
853 "The 'softening' parameter is missing for species "//trim(spec%get_label()))
854 end if
855
857 if (.not. parameter_defined(option__species__potential_formula)) then
858 call messages_input_error(namespace, 'Species', &
859 "The 'potential_formula' parameter is missing for species '"//trim(spec%get_label())//"'")
860 end if
861
863 if (.not. parameter_defined(option__species__density_formula)) then
864 call messages_input_error(namespace, 'Species', &
865 "The 'density_formula' parameter is missing for species '"//trim(spec%get_label())//"'")
866 end if
867
868 type is(species_from_file_t)
869 if( .not. (parameter_defined(option__species__file) .or. parameter_defined(option__species__db_file))) then
870 call messages_input_error(namespace, 'Species', &
871 "The 'file' or 'db_file' parameter is missing for species '"//trim(spec%get_label())//"'")
872 end if
873
874 type is(jellium_slab_t)
875 if (.not. parameter_defined(option__species__thickness)) then
876 call messages_input_error(namespace, 'Species', &
877 "The 'thickness' parameter is missing for species '"//trim(spec%get_label())//"'")
878 end if
879
880 type is(jellium_charge_t)
881 if (.not. parameter_defined(option__species__volume)) then
882 call messages_input_error(namespace, 'Species', &
883 "The 'volume' parameter is missing for species '"//trim(spec%get_label())//"'")
884 end if
885
886 type is(pseudopotential_t)
887 if (parameter_defined(option__species__lmax) .and. parameter_defined(option__species__lloc)) then
888 if (spec%get_user_lloc() > spec%get_user_lmax()) then
889 call messages_input_error(namespace, 'Species', &
890 "the 'lloc' parameter cannot be larger than the 'lmax' parameter in species "//trim(spec%get_label()))
891 end if
892 end if
893
894 if(.not. (parameter_defined(option__species__file) .or. parameter_defined(option__species__db_file))) then
895 ! we need to read the species from the pseudopotential set
896
897 !if the set was not defined, use the default set
898 if (.not. parameter_defined(option__species__set)) then
899 spec%pseudopotential_set_id = factory%default_pseudopotential_set_id
900 spec%pseudopotential_set = factory%default_pseudopotential_set
901 end if
902
903 call read_from_set(spec, spec%pseudopotential_set_id, spec%pseudopotential_set, set_read_data)
904
905 if (set_read_data == 0) then
906 call messages_write('Species '//trim(spec%get_label())//' is not defined in the requested pseudopotential set.', &
907 new_line=.true.)
908 call messages_write('( '//trim(get_set_directory(spec%pseudopotential_set_id))//' )')
909 call messages_fatal(namespace=namespace)
910 end if
911 end if
912
913 end select
914
915 select type (spec)
916 class is(pseudopotential_t)
918
919 class is(full_anc_t)
921
922 ! If z_val was not specified, we set it to be z
923 if (spec%get_zval() < m_zero) then
924 call spec%set_zval(spec%get_z())
925 end if
926
927 class is(full_gaussian_t)
929
930 ! If z_val was not specified, we set it to be z
931 if (spec%get_zval() < m_zero) then
932 call spec%set_zval(spec%get_z())
933 end if
934
935 class is(full_delta_t)
937
938 ! If z_val was not specified, we set it to be z
939 if (spec%get_zval() < m_zero) then
940 call spec%set_zval(spec%get_z())
941 end if
942
943 class default
944 if (.not. parameter_defined(option__species__mass)) then
945 call spec%set_mass(m_one)
946 call messages_write('Info: default mass for species '//trim(spec%get_label())//':')
947 call messages_write(spec%get_mass())
948 call messages_write(' amu.')
949 call messages_info(namespace=namespace)
950 end if
951
952 if (.not. parameter_defined(option__species__vdw_radius)) then
953 call spec%set_vdw_radius(m_zero)
954 call messages_write('Info: default vdW radius for species '//trim(spec%get_label())//':')
955 call messages_write(spec%get_vdw_radius())
956 call messages_write(' [b]')
957 call messages_info(namespace=namespace)
958 end if
959
960 if (.not. parameter_defined(option__species__valence)) then
961 if (spec%is_user_defined()) then
962 call spec%set_zval(m_zero)
963 else
964 call messages_input_error(namespace, 'Species', &
965 "The 'valence' parameter is missing for species '"//trim(spec%get_label())//"'")
966 end if
967 end if
968
969 end select
970
971 call iihash_end(read_parameters)
972
973 pop_sub(read_from_block)
974
975 contains
976
977 logical function parameter_defined(param) result(defined)
978 integer(int64), intent(in) :: param
979
980 integer :: tmp
981
983
984 tmp = iihash_lookup(read_parameters, int(-param), defined)
985
987 end function parameter_defined
988
989 !------------------------------------------------------
990 subroutine check_duplication(param)
991 integer(int64), intent(in) :: param
992
994
995 if (parameter_defined(param)) then
996 call messages_input_error(namespace, 'Species', "Duplicated parameter in species '"//trim(spec%get_label())//"'")
997 end if
998
999 call iihash_insert(read_parameters, int(-param), 1)
1000
1002 end subroutine check_duplication
1003
1004
1005 !------------------------------------------------------
1006 subroutine check_real_atom_species()
1007
1008 call element_init(element, get_symbol(spec%get_label()))
1009
1010 if (.not. element_valid(element)) then
1011 call messages_write('Cannot determine the element for species '//trim(spec%get_label())//'.')
1012 call messages_fatal(namespace=namespace)
1013 end if
1014
1015 call spec%set_z(real(element_atomic_number(element), real64))
1016
1017 if (spec%get_mass() < m_zero) then
1018 call spec%set_mass(element_mass(element))
1019 call messages_write('Info: default mass for species '//trim(spec%get_label())//':')
1020 call messages_write(spec%get_mass())
1021 call messages_write(' amu.')
1022 call messages_info(namespace=namespace)
1023 end if
1024
1025 if (spec%get_vdw_radius() < m_zero) then
1026 call spec%set_vdw_radius(element_vdw_radius(element))
1027 if (spec%get_vdw_radius() < m_zero) then
1028 call spec%set_vdw_radius(m_zero)
1029 call messages_write("The default vdW radius for species '"//trim(spec%get_label())//"' is not defined.", &
1030 new_line = .true.)
1031 call messages_write("You can specify the vdW radius in %Species block.")
1032 call messages_warning(namespace=namespace)
1033 end if
1034 call messages_write('Info: default vdW radius for species '//trim(spec%get_label())//':')
1035 call messages_write(spec%get_vdw_radius())
1036 call messages_write(' [b]')
1037 call messages_info(namespace=namespace)
1038 end if
1039
1040 call element_end(element)
1041
1042
1043 end subroutine check_real_atom_species
1044 end function read_from_block
1045
1046end module species_factory_oct_m
1047
1048!! Local Variables:
1049!! mode: f90
1050!! coding: utf-8
1051!! End:
subroutine check_duplication(param)
logical function parameter_defined(param)
Delete the C++ object.
Definition: pseudo_set.F90:156
Nullify the C++ pointer.
Definition: pseudo_set.F90:146
logical function, public element_valid(self)
Definition: element.F90:193
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_epsilon
Definition: global.F90:207
type(conf_t), public conf
Global instance of Octopus configuration.
Definition: global.F90:181
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
This module implements a simple hash table for non-negative integer keys and integer values.
Definition: iihash.F90:127
subroutine, public iihash_end(h)
Free a hash table.
Definition: iihash.F90:186
subroutine, public iihash_insert(h, key, val)
Insert a (key, val) pair into the hash table h.
Definition: iihash.F90:208
integer function, public iihash_lookup(h, key, found)
Look up a value in the hash table h. If found is present, it indicates if key could be found in the t...
Definition: iihash.F90:233
subroutine, public iihash_init(h)
Initialize a hash table h.
Definition: iihash.F90:163
Definition: io.F90:116
integer, parameter, public species_charge_density
user-defined function for charge density
Definition: jellium.F90:141
integer, parameter, public species_jellium_charge_density
jellium volume read from file
Definition: jellium.F90:141
integer, parameter, public species_jellium
jellium sphere.
Definition: jellium.F90:141
integer, parameter, public species_from_file
Definition: jellium.F90:141
integer, parameter, public species_usdef
user-defined function for local potential
Definition: jellium.F90:141
integer, parameter, public species_jellium_slab
jellium slab.
Definition: jellium.F90:141
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
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 messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
subroutine, public parse_block_string(blk, l, c, res, convert_to_c)
Definition: parser.F90:810
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public pseudo_set_init(pseudo_set, dirname, ierr)
Definition: pseudo_set.F90:191
character(len=max_path_len) function, public get_set_directory(set_id)
subroutine, public read_from_set(spec, set_id, set, read_data)
Creates a pseudopotential type from a set.
integer, parameter, public species_pseudo
pseudopotential
integer, parameter, public species_pspio
pseudopotential parsed by pspio library
subroutine, public share_directory_set(dir)
class(species_t) function, pointer species_factory_create_from_input(factory, namespace, label, index)
Reads the information (from the input file) about a species_t variable, initializing part of it (it h...
subroutine, public species_factory_init(factory, namespace)
class(species_t) function, pointer read_from_block(factory, namespace, blk, row, label, index, read_data)
Parses the species block for a given species.
subroutine, public species_factory_end(factory)
character(kind=c_char) function, dimension(label_len+1), public get_symbol(label)
Definition: species.F90:522
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
subroutine check_real_atom_species()
An abstract class for species. Derived classes include jellium, all electron, and pseudopotential spe...
Definition: species.F90:147
int true(void)