1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      module chemical
9
10      use sys
11      use precision, only: dp
12
13      implicit none
14
15      private
16
17      public :: atomic_number
18      public :: number_of_species, species_label
19      public :: is_floating, is_bessel, is_synthetic
20      public :: read_chemical_types, print_chemical_type
21
22      ! Public due to Bcast routines
23      public :: chemical_types, chemical_list
24
25      ! Species information
26      type chemical_types
27         integer                    :: no_of_species
28         character(len=20), pointer :: spec_label(:)
29         integer, pointer           :: z(:)
30      end type chemical_types
31
32      type(chemical_types), save :: chemical_list
33
34
35      CONTAINS
36
37      subroutine check(i)
38      integer, intent(in) :: i
39      if (i.lt.0 .or. i.gt.chemical_list%no_of_species)
40     $     call die("Wrong species number requested")
41      end subroutine check
42
43      function number_of_species()
44      integer number_of_species
45      number_of_species = chemical_list%no_of_species
46      end function number_of_species
47
48      function species_label(i)
49      character(len=20) species_label
50      integer, intent(in)  :: i
51
52      call check(i)
53      species_label = chemical_list%spec_label(i)
54      end function species_label
55
56      function atomic_number(i)
57      integer atomic_number
58      integer, intent(in)  :: i
59
60      call check(i)
61      atomic_number = chemical_list%z(i)
62      end function atomic_number
63! -------
64      function is_floating(i)
65      logical is_floating
66      integer, intent(in)  :: i
67
68      call check(i)
69      is_floating = (chemical_list%z(i) .le. 0)
70      end function is_floating
71! -------
72
73      function is_bessel(i)
74      logical is_bessel
75      integer, intent(in)  :: i
76
77      call check(i)
78      is_bessel = (chemical_list%z(i) .eq. -100)
79      end function is_bessel
80! -------
81!     Checks whether we are dealing with a synthetic atom
82!
83      function is_synthetic(i)
84      logical is_synthetic
85      integer, intent(in)  :: i
86
87      call check(i)
88      ! Note that we could have a synthetic ghost atom, with
89      ! z <= -200
90      is_synthetic = (abs(chemical_list%z(i)) .gt. 200)
91      end function is_synthetic
92!---
93      subroutine read_chemical_types(silent)
94
95      use parallel,    only : Node
96      use fdf
97
98      logical, intent(in), optional :: silent ! default .false.
99
100      integer nsp, isp
101      integer ns_read
102
103      type(block_fdf)            :: bfdf
104      type(parsed_line), pointer :: pline
105
106      character(len=20) :: label
107      character(len=256) :: msg
108
109      integer :: z, is
110      logical :: found
111      logical :: lsilent
112
113      ! Determine whether we should be silent
114      lsilent = .false.
115      if ( present(silent) ) lsilent = silent
116      if ( Node /= 0 ) lsilent = .true.
117
118      ! Default to 0
119      nsp = fdf_integer('Number_of_species',0)
120
121      ! The most important thing to find is
122      ! the block containing the species
123      found = fdf_block('Chemical_species_label', bfdf)
124      if (.not. found )
125     $     call die("Block Chemical_species_label does not exist.")
126
127      if ( nsp == 0 ) then
128         ns_read = fdf_block_linecount('Chemical_species_label', 'iin')
129      else
130         ns_read = nsp
131      end if
132      ! If they are not equal we notify the user
133      if ( nsp /= ns_read ) then
134         nsp = ns_read
135      end if
136      if ( nsp == 0 ) call die("No species found!!!")
137
138      allocate(chemical_list%spec_label(nsp))
139      allocate(chemical_list%z(nsp))
140      chemical_list%no_of_species = nsp
141
142      ns_read = 0
143      do while( fdf_bline(bfdf,pline) )
144         if ( .not. fdf_bmatch(pline,'iin') ) cycle
145
146         ns_read = ns_read + 1
147
148         ! Get species information
149         isp = fdf_bintegers(pline,1)
150         label = fdf_bnames(pline,1)
151         z = fdf_bintegers(pline,2)
152
153         ! We cannot test label names in this
154         ! loop as isp may be non-linear
155         if ( isp < 1 .or. nsp < isp )
156     $ call die("Wrong specnum in Chemical_species_label")
157
158         chemical_list%z(isp) = z
159         chemical_list%spec_label(isp) = label
160
161      end do
162      if ( ns_read /= nsp )
163     &     call die("Not enough species in block")
164
165      if ( .not. lsilent ) then
166         ! Align output, always
167         do isp = 1 , nsp
168            call print_chemical_type(isp)
169         end do
170         write(*,*) ! new-line
171      end if
172
173      ! Check that none of the chemical species are the
174      ! same
175      if ( nsp > 1 ) then
176       do z = 1 , nsp - 1
177        do is = z+1, nsp
178         if (trim(species_label(z))==trim(species_label(is))) then
179          write(msg,'(2(a,i0,2a),2a)')
180     & "Specie index/label = ",
181     & z,'/',trim(species_label(z)), " has same label as ",
182     & is,'/',trim(species_label(is)),". ",
183     & " Use a different one for hygienic reasons."
184          call die(trim(msg))
185         end if
186        end do
187       end do
188      end if
189
190      end subroutine read_chemical_types
191
192      subroutine print_chemical_type(isp)
193      integer, intent(in)  :: isp
194
195      character(len=256) :: label
196      integer :: z
197      logical :: floating, bessel, synthetic
198
199      label = species_label(isp)
200      z = atomic_number(isp)
201
202      floating  =     z <= 0
203      bessel    =     z == -100
204      synthetic = abs(z) > 200
205
206      ! Bessel *must* be checked first!
207      if ( bessel ) then
208         write(6,'(a,i3,3a)') 'Species number: ',isp,
209     $        ' Label: ', trim(label),
210     $        ' (floating Bessel functions)'
211      else if ( floating ) then
212         write(6,'(a,i3,a,i4,3a)') 'Species number: ',isp,
213     $        ' Atomic number: ', z,
214     $        ' Label: ', trim(label),
215     $        ' (floating PAOs)'
216      else
217         write(6,'(a,i3,a,i4,2a)') 'Species number: ',isp,
218     $        ' Atomic number: ', z,
219     $        ' Label: ', trim(label)
220      end if
221      end subroutine print_chemical_type
222
223      end module chemical
224