GCC Code Coverage Report


Directory: src/fortran/lib/
File: mod_geom_rw.f90
Date: 2025-06-15 07:27:34
Exec Total Coverage
Lines: 575 1004 57.3%
Functions: 0 0 -%
Branches: 591 2177 27.1%

Line Branch Exec Source
1 module raffle__geom_rw
2 !! Module to store, read and write geometry files
3 !!
4 !! This module contains the procedures to read and write geometry files.
5 !! It also contains the derived types used to store the geometry data.
6 use raffle__constants, only: pi,real32
7 use raffle__io_utils, only: stop_program, print_warning
8 use raffle__misc, only: to_upper, to_lower, jump, icount, strip_null
9 use raffle__misc_linalg, only: inverse_3x3
10 implicit none
11
12
13 private
14
15 public :: igeom_input, igeom_output
16 public :: basis_type, species_type
17 public :: geom_read, geom_write
18 public :: get_element_properties
19
20
21 integer :: igeom_input = 1
22 !! geometry input file format
23 !! 1 = VASP
24 ! 2 = CASTEP
25 !! 3 = Quantum Espresso
26 !! 4 = CRYSTAL
27 !! 5 = XYZ
28 !! 6 = extended XYZ
29 integer :: igeom_output = 1
30 !! geometry output file format
31
32 type :: species_type
33 !! Derived type to store information about a species/element.
34 integer, allocatable, dimension(:) :: atom_idx
35 !! The indices of the atoms of this species in the basis.
36 !! For ASE compatibility.
37 logical, dimension(:), allocatable :: atom_mask
38 !! The mask of the atoms of this species.
39 real(real32), allocatable ,dimension(:,:) :: atom
40 !! The atomic positions of the species.
41 real(real32) :: mass
42 !! The mass of the species.
43 real(real32) :: charge
44 !! The charge of the species.
45 real(real32) :: radius
46 !! The radius of the species.
47 character(len=3) :: name
48 !! The name of the species.
49 integer :: num
50 !! The number of atoms of this species.
51 end type species_type
52 type :: basis_type
53 !! Derived type to store information about a basis.
54 type(species_type), allocatable, dimension(:) :: spec
55 !! Information about each species in the basis.
56 integer :: nspec = 0
57 !! The number of species in the basis.
58 integer :: natom = 0
59 !! The number of atoms in the basis.
60 real(real32) :: energy = 0._real32
61 !! The energy of the basis.
62 real(real32) :: lat(3,3) = 0._real32
63 !! The lattice vectors of the basis.
64 logical :: lcart = .false.
65 !! Boolean whether the basis is in cartesian coordinates.
66 logical, dimension(3) :: pbc = .true.
67 !! Boolean whether the basis has periodic boundary conditions.
68 character(len=128) :: sysname = "default"
69 !! The name of the system.
70 contains
71 procedure, pass(this) :: allocate_species
72 !! Procedure to allocate the species in the basis.
73 procedure, pass(this) :: convert
74 !! Procedure to convert the basis to cartesian coordinates.
75 procedure, pass(this) :: copy
76 !! Procedure to copy the basis.
77 procedure, pass(this) :: set_atom_mask
78 !! Procedure to set the atom mask of the basis.
79 procedure, pass(this) :: get_lattice_constants
80 !! Procedure to get the lattice constants of the basis.
81 procedure, pass(this) :: add_atom
82 !! Procedure to add an atom to the basis.
83 procedure, pass(this) :: remove_atom
84 !! Procedure to remove an atom from the basis.
85 procedure, pass(this) :: remove_atoms
86 !! Procedure to remove atoms from the basis.
87 end type basis_type
88
89
90 interface basis_type
91 module function init_basis_type(basis) result(output)
92 !! Initialise the basis type.
93 type(basis_type), intent(in), optional :: basis
94 !! Optional. Basis to copy.
95 type(basis_type) :: output
96 !! The basis to initialise.
97 end function init_basis_type
98 end interface basis_type
99
100
101
102 contains
103
104 !###############################################################################
105 1 module function init_basis_type(basis) result(output)
106 !! Initialise the basis type.
107 implicit none
108
109 ! Arguments
110 type(basis_type), intent(in), optional :: basis
111 !! Optional. Basis to copy.
112 type(basis_type) :: output
113 !! The basis to initialise.
114
115 1 if(present(basis)) call output%copy(basis)
116
117
7/8
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
17 end function init_basis_type
118 !###############################################################################
119
120
121 !###############################################################################
122 4 subroutine allocate_species( &
123 this, num_species, &
124
13/24
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 3 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 4 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
4 species_symbols, species_count, atoms, atom_idx_list )
125 !! Allocate the species in the basis.
126 implicit none
127
128 ! Arguments
129 class(basis_type), intent(inout) :: this
130 !! Parent. The basis to allocate the species in.
131 integer, intent(in), optional :: num_species
132 !! Optional. The number of species in the basis.
133 character(3), dimension(:), intent(in), optional :: species_symbols
134 !! Optional. The symbols of the species.
135 integer, dimension(:), intent(in), optional :: species_count
136 !! Optional. The number of atoms of each species.
137 real(real32), dimension(:,:), intent(in), optional :: atoms
138 !! Optional. The atomic positions of the species.
139 integer, dimension(:), intent(in), optional :: atom_idx_list
140 !! Optional. The indices of the atoms of the species.
141
142 ! Local variables
143 integer :: i, j, istart, iend
144 !! Loop index.
145
146
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
4 if(present(num_species)) this%nspec = num_species
147
148
11/14
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 4 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
13 if(allocated(this%spec)) deallocate(this%spec)
149
13/24
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 4 times.
✓ Branch 17 taken 6 times.
✓ Branch 18 taken 4 times.
✓ Branch 19 taken 6 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 6 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 6 times.
10 allocate(this%spec(this%nspec))
150
151
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
4 species_check: if(present(species_symbols))then
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(size(species_symbols).ne.this%nspec) exit species_check
153
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
3 this%spec(:)%name = species_symbols
154 end if species_check
155
156
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
4 natom_check: if(present(species_count))then
157
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(size(species_count).ne.this%nspec) exit natom_check
158
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
3 this%spec(:)%num = species_count
159 1 istart = 1
160
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
3 do i = 1, this%nspec
161 2 iend = istart + this%spec(i)%num - 1
162
9/16
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 3 times.
✓ Branch 18 taken 2 times.
5 allocate(this%spec(i)%atom_mask(this%spec(i)%num), source = .true.)
163
7/14
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
2 allocate(this%spec(i)%atom_idx(this%spec(i)%num))
164
8/16
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
2 allocate(this%spec(i)%atom(this%spec(i)%num,3))
165
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if(present(atoms))then
166
7/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 6 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 9 times.
✓ Branch 11 taken 6 times.
19 this%spec(i)%atom = atoms(istart:iend,:3)
167 end if
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if(present(atom_idx_list))then
169 this%spec(i)%atom_idx = atom_idx_list(istart:iend)
170 else
171
10/16
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✓ Branch 16 taken 2 times.
11 this%spec(i)%atom_idx = [ ( j, j = istart, iend, 1 ) ]
172 end if
173 3 istart = iend + 1
174 end do
175 end if natom_check
176
177
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4 times.
10 do i = 1, this%nspec
178 call get_element_properties( &
179 this%spec(i)%name, &
180 mass = this%spec(i)%mass, &
181 charge = this%spec(i)%charge, &
182 10 radius = this%spec(i)%radius )
183 end do
184
185 4 end subroutine allocate_species
186 !###############################################################################
187
188
189 !###############################################################################
190
7/18
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 9 times.
✓ Branch 13 taken 3 times.
✓ Branch 14 taken 27 times.
✓ Branch 15 taken 9 times.
✓ Branch 16 taken 9 times.
✓ Branch 17 taken 3 times.
48 subroutine geom_read(UNIT, basis, length, iostat)
191 !! Read geometry from a file.
192 implicit none
193
194 ! Arguments
195 integer, intent(in) :: UNIT
196 !! The unit number of the file to read from.
197 type(basis_type), intent(out) :: basis
198 !! The basis to read the geometry into.
199 integer, optional, intent(in) :: length
200 !! Optional. The dimension of the basis atom positions.
201 integer, optional, intent(out) :: iostat
202 !! Optional. The I/O status of the read.
203
204 ! Local variables
205 integer :: i
206 !! Loop index.
207 integer :: length_
208 !! The dimension of the basis atom positions.
209 integer :: iostat_
210 !! The I/O status of the read.
211
212
213 3 length_ = 3
214 3 iostat_ = 0
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if(present(length)) length_=length
216
217 2 select case(igeom_input)
218 case(1)
219 2 call VASP_geom_read(UNIT, basis, length_, iostat_)
220 case(2)
221 call CASTEP_geom_read(UNIT, basis, length_)
222 case(3)
223 call QE_geom_read(UNIT, basis, length_)
224 case(4)
225 call stop_program("Not yet set up for CRYSTAL")
226 return
227 case(5)
228 call XYZ_geom_read(UNIT, basis, length_, iostat_)
229 call print_warning("XYZ file format does not contain lattice data")
230 case(6)
231
2/7
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
3 call extXYZ_geom_read(UNIT, basis, length_, iostat_)
232 end select
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if(iostat_.ne.0) then
234 if(present(iostat)) iostat = iostat_
235 return
236 else
237
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 if(present(iostat)) iostat = 0
238 end if
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if(length_.eq.4)then
240 do i = 1, basis%nspec
241 basis%spec(i)%atom(:,4)=1._real32
242 end do
243 end if
244
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 3 times.
6 do i = 1, basis%nspec
245 call get_element_properties( &
246 basis%spec(i)%name, &
247 mass = basis%spec(i)%mass, &
248 charge = basis%spec(i)%charge, &
249 6 radius = basis%spec(i)%radius )
250 end do
251
252 end subroutine geom_read
253 !###############################################################################
254
255
256 !###############################################################################
257 2 subroutine geom_write(UNIT, basis)
258 !! Write geometry to a file.
259 implicit none
260
261 ! Arguments
262 integer, intent(in) :: UNIT
263 !! The unit number of the file to write to.
264 class(basis_type), intent(in) :: basis
265 !! The basis to write the geometry from.
266
267 ! MAKE IT CHANGE HERE IF USER SPECIFIES LCART OR NOT
268 ! AND GIVE IT THE CASTEP AND QE OPTION OF LABC !
269
270 1 select case(igeom_output)
271 case(1)
272 1 call VASP_geom_write(UNIT,basis)
273 case(2)
274 call CASTEP_geom_write(UNIT,basis)
275 case(3)
276 call QE_geom_write(UNIT,basis)
277 case(4)
278 call stop_program("ERROR: Not yet set up for CRYSTAL")
279 return
280 case(5)
281 call XYZ_geom_write(UNIT,basis)
282 case(6)
283
2/7
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
2 call extXYZ_geom_write(UNIT,basis)
284 end select
285
286 end subroutine geom_write
287 !###############################################################################
288
289
290 !###############################################################################
291
7/18
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 18 times.
✓ Branch 15 taken 6 times.
✓ Branch 16 taken 6 times.
✓ Branch 17 taken 2 times.
32 subroutine VASP_geom_read(UNIT, basis, length, iostat)
292 !! Read the structure in vasp poscar style format.
293 implicit none
294
295 ! Arguments
296 integer, intent(in) :: UNIT
297 !! The unit number of the file to read from.
298 type(basis_type), intent(out) :: basis
299 !! The basis to read the geometry into.
300 integer, intent(in), optional :: length
301 !! Optional. The dimension of the basis atom positions.
302 integer, intent(out), optional :: iostat
303 !! Optional. The I/O status of the read.
304
305 integer :: Reason
306 !! The I/O status of the read.
307 integer :: pos, count, natom
308 !! Temporary integer variables.
309 real(real32) :: scal
310 !! The scaling factor of the lattice.
311 character(len=100) :: lspec
312 !! The species names and number of each atomic species.
313 character(len=1024) :: buffer
314 !! Temporary character variable.
315 integer :: i, j, k
316 !! Loop index.
317 integer :: length_
318 !! The dimension of the basis atom positions.
319 integer :: iostat_
320 !! The I/O status of the read.
321
322
323 2 length_ = 3
324 2 iostat_ = 0
325 !---------------------------------------------------------------------------
326 ! determine dimension of basis (include translation dimension for symmetry?)
327 !---------------------------------------------------------------------------
328
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if(present(length)) length_ = length
329
330
331 !---------------------------------------------------------------------------
332 ! read system name
333 !---------------------------------------------------------------------------
334 2 read(UNIT,'(A)',iostat=Reason) basis%sysname
335
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if(Reason.ne.0)then
336 write(0,'("ERROR: The file is not in POSCAR format.")')
337 write(0,*) "Expected system name, got: ",trim(basis%sysname)
338 iostat_ = 1
339 if(present(iostat)) iostat = iostat_
340 return
341 end if
342 2 read(UNIT,*) scal
343
344
345 !---------------------------------------------------------------------------
346 ! read lattice
347 !---------------------------------------------------------------------------
348
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
8 do i = 1, 3
349
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 18 times.
✓ Branch 4 taken 6 times.
26 read(UNIT,*) (basis%lat(i,j),j=1,3)
350 end do
351
4/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 6 times.
26 basis%lat=scal*basis%lat
352
353
354 !---------------------------------------------------------------------------
355 ! read species names and number of each atomic species
356 !---------------------------------------------------------------------------
357 2 read(UNIT,'(A)') lspec
358 2 basis%nspec = icount(lspec)
359
13/24
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 2 times.
4 allocate(basis%spec(basis%nspec))
360
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if(verify(lspec,' 0123456789').ne.0) then
361 2 count=0;pos=1
362 2 speccount: do
363 4 i=verify(lspec(pos:), ' ')
364
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
4 if (i.eq.0) exit speccount
365 2 count=count+1
366 2 pos=i+pos-1
367 2 i=scan(lspec(pos:), ' ')
368
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if (i.eq.0) exit speccount
369
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 basis%spec(count)%name=lspec(pos:pos+i-1)
370 2 pos=i+pos-1
371 end do speccount
372
373
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
4 read(UNIT,*) (basis%spec(j)%num,j=1,basis%nspec)
374 else !only numbers
375 do count = 1, basis%nspec
376 write(basis%spec(count)%name,'(I0)') count
377 end do
378 read(lspec,*) (basis%spec(j)%num,j=1,basis%nspec)
379 end if
380
381
382 !---------------------------------------------------------------------------
383 ! determines whether input basis is in direct or cartesian coordinates
384 !---------------------------------------------------------------------------
385 2 basis%lcart=.false.
386 2 read(UNIT,'(A)') buffer
387
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
2 buffer = to_lower(buffer)
388
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
2 if(verify(trim(buffer),'direct').eq.0) basis%lcart=.false.
389
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
2 if(verify(trim(buffer),'cartesian').eq.0) basis%lcart=.true.
390
391
392 !---------------------------------------------------------------------------
393 ! read basis
394 !---------------------------------------------------------------------------
395 2 natom = 0
396
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
4 do i = 1, basis%nspec
397
7/14
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
2 allocate(basis%spec(i)%atom_idx(basis%spec(i)%num))
398
9/16
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 16 times.
✓ Branch 18 taken 2 times.
18 allocate(basis%spec(i)%atom_mask(basis%spec(i)%num), source = .true.)
399
9/18
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
2 allocate(basis%spec(i)%atom(basis%spec(i)%num,length_))
400
4/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 6 times.
56 basis%spec(i)%atom(:,:)=0._real32
401
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
20 do j = 1, basis%spec(i)%num
402 16 natom = natom + 1
403 16 basis%spec(i)%atom_idx(j) = natom
404
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 48 times.
✓ Branch 4 taken 16 times.
66 read(UNIT,*) (basis%spec(i)%atom(j,k),k=1,3)
405 end do
406 end do
407
408
409 !---------------------------------------------------------------------------
410 ! convert basis if in cartesian coordinates
411 !---------------------------------------------------------------------------
412
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if(basis%lcart) call basis%convert()
413
414
415 !---------------------------------------------------------------------------
416 ! normalise basis to between 0 and 1 in direct coordinates
417 !---------------------------------------------------------------------------
418
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
4 do i = 1, basis%nspec
419
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
20 do j = 1, basis%spec(i)%num
420
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 16 times.
66 do k = 1, 3
421 basis%spec(i)%atom(j,k)=&
422
1/2
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
64 basis%spec(i)%atom(j,k)-floor(basis%spec(i)%atom(j,k))
423 end do
424 end do
425 end do
426
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
4 basis%natom=sum(basis%spec(:)%num)
427
428
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if(present(iostat)) iostat = iostat_
429
430 end subroutine VASP_geom_read
431 !###############################################################################
432
433
434 !###############################################################################
435 1 subroutine VASP_geom_write(UNIT, basis, cartesian)
436 !! Write the structure in vasp poscar style format.
437 implicit none
438
439 ! Arguments
440 integer, intent(in) :: UNIT
441 !! The unit number of the file to write to.
442 class(basis_type), intent(in) :: basis
443 !! The basis to write the geometry from.
444 logical, intent(in), optional :: cartesian
445 !! Optional. Whether to write the basis in cartesian coordinates.
446
447 ! Local variables
448 integer :: i,j
449 !! Loop index.
450 character(100) :: fmt
451 !! Format string.
452 character(10) :: string
453 !! String to determine whether to write in direct or cartesian coordinates.
454
455
456 1 string="Direct"
457
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(present(cartesian))then
458 if(cartesian) string="Cartesian"
459 end if
460
461
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
1 write(UNIT,'(A)') trim(adjustl(basis%sysname))
462 1 write(UNIT,'(F15.9)') 1._real32
463
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
4 do i = 1, 3
464 4 write(UNIT,'(3(F15.9))') basis%lat(i,:)
465 end do
466 1 write(fmt,'("(",I0,"(A,1X))")') basis%nspec
467
4/6
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
2 write(UNIT,trim(adjustl(fmt))) (adjustl(basis%spec(j)%name),j=1,basis%nspec)
468 1 write(fmt,'("(",I0,"(I0,5X))")') basis%nspec
469
4/6
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
2 write(UNIT,trim(adjustl(fmt))) (basis%spec(j)%num,j=1,basis%nspec)
470
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
1 write(UNIT,'(A)') trim(adjustl(string))
471
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
2 do i = 1, basis%nspec
472
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
10 do j = 1, basis%spec(i)%num
473
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 8 times.
33 write(UNIT,'(3(F15.9))') basis%spec(i)%atom(j,1:3)
474 end do
475 end do
476
477 1 end subroutine VASP_geom_write
478 !###############################################################################
479
480
481 !###############################################################################
482 subroutine QE_geom_read(UNIT,basis,length)
483 !! Read the structure in Quantum Espresso style format.
484 implicit none
485
486 ! Arguments
487 integer, intent(in) :: UNIT
488 !! The unit number of the file to read from.
489 type(basis_type), intent(out) :: basis
490 !! The basis to read the geometry into.
491 integer, intent(in), optional :: length
492 !! Optional. The dimension of the basis atom positions.
493
494 ! Local variables
495 integer :: Reason
496 !! The I/O status of the read.
497 integer :: i, j, k, iline
498 !! Loop index.
499 integer :: length_ = 3
500 !! The dimension of the basis atom positions.
501 integer, dimension(1000) :: tmp_natom
502 !! Temporary array to store the number of atoms of each species.
503 real(real32), dimension(3) :: tmpvec
504 !! Temporary array to store the atomic positions.
505 character(len=3) :: ctmp
506 !! Temporary character variable.
507 character(256) :: stop_msg
508 !! Error message.
509 character(len=3), dimension(1000) :: tmp_spec
510 !! Temporary array to store the species names.
511 character(len=1024) :: buffer, buffer2
512 !! Temporary character variables.
513
514
515 !---------------------------------------------------------------------------
516 ! determine dimension of basis (include translation dimension for symmetry?)
517 !---------------------------------------------------------------------------
518 if(present(length)) length_ = length
519 basis%lcart = .false.
520 basis%sysname = "Converted_from_geom_file"
521
522
523 !---------------------------------------------------------------------------
524 ! read lattice
525 !---------------------------------------------------------------------------
526 rewind UNIT
527 cellparam: do
528 read(UNIT,'(A)',iostat=Reason) buffer
529 if(Reason.ne.0)then
530 call stop_program( &
531 "An issue with the QE input file format has been encountered." &
532 )
533 return
534 end if
535 if(index(trim(buffer),"ibrav").ne.0)then
536 write(stop_msg,*) &
537 "Internal error in QE_geom_read" // &
538 achar(13) // achar(10) // &
539 " Subroutine not yet set up to read IBRAV lattices"
540 call stop_program(stop_msg)
541 return
542 end if
543 if(verify("CELL_PARAMETERS",buffer).eq.0) then
544 exit cellparam
545 end if
546 end do cellparam
547 do i = 1, 3
548 read(UNIT,*) (basis%lat(i,j),j=1,3)
549 end do
550
551
552 !---------------------------------------------------------------------------
553 ! determines whether input basis is in direct or cartesian coordinates
554 !---------------------------------------------------------------------------
555 iline=0
556 rewind UNIT
557 basfind: do
558 read(UNIT,'(A)',iostat=Reason) buffer
559 iline=iline+1
560 if(verify("ATOMIC_POSITIONS",buffer).eq.0)then
561 backspace(UNIT)
562 read(UNIT,*) buffer,buffer2
563 if(verify("crystal",buffer2).eq.0) basis%lcart = .false.
564 if(verify("angstrom",buffer2).eq.0) basis%lcart = .true.
565 exit basfind
566 end if
567 end do basfind
568
569
570 !---------------------------------------------------------------------------
571 ! read basis
572 !---------------------------------------------------------------------------
573 basis%natom = 0
574 basis%nspec = 0
575 tmp_natom = 1
576 basread: do
577 read(UNIT,'(A)',iostat=Reason) buffer
578 read(buffer,*) ctmp
579 if(Reason.ne.0) exit
580 if(trim(ctmp).eq.'') exit
581 if(verify(buffer,' 0123456789').eq.0) exit
582 basis%natom = basis%natom + 1
583 if(.not.any(tmp_spec(1:basis%nspec).eq.ctmp))then
584 basis%nspec = basis%nspec + 1
585 tmp_spec(basis%nspec) = ctmp
586 else
587 where(tmp_spec(1:basis%nspec).eq.ctmp)
588 tmp_natom(1:basis%nspec) = tmp_natom(1:basis%nspec) + 1
589 end where
590 end if
591 end do basread
592
593 allocate(basis%spec(basis%nspec))
594 basis%spec(1:basis%nspec)%name = tmp_spec(1:basis%nspec)
595 do i = 1, basis%nspec
596 basis%spec(i)%num = 0
597 allocate(basis%spec(i)%atom_idx(tmp_natom(i)))
598 allocate(basis%spec(i)%atom_mask(tmp_natom(i)), source = .true.)
599 allocate(basis%spec(i)%atom(tmp_natom(i),length_))
600 end do
601
602 call jump(UNIT,iline)
603 basread2: do i = 1, basis%natom
604 read(UNIT,*,iostat=Reason) ctmp,tmpvec(1:3)
605 do j = 1, basis%nspec
606 if(basis%spec(j)%name.eq.ctmp)then
607 basis%spec(j)%num = basis%spec(j)%num + 1
608 basis%spec(j)%atom_idx(basis%spec(j)%num) = i
609 basis%spec(j)%atom(basis%spec(j)%num,1:3) = tmpvec(1:3)
610 exit
611 end if
612 end do
613 end do basread2
614
615
616 !---------------------------------------------------------------------------
617 ! convert basis if in cartesian coordinates
618 !---------------------------------------------------------------------------
619 if(basis%lcart) call basis%convert()
620
621
622 !---------------------------------------------------------------------------
623 ! normalise basis to between 0 and 1 in direct coordinates
624 !---------------------------------------------------------------------------
625 do i = 1, basis%nspec
626 do j = 1, basis%spec(i)%num
627 do k = 1, 3
628 basis%spec(i)%atom(j,k) = &
629 basis%spec(i)%atom(j,k) - floor( basis%spec(i)%atom(j,k) )
630 end do
631 end do
632 end do
633 basis%natom=sum(basis%spec(:)%num)
634
635 end subroutine QE_geom_read
636 !###############################################################################
637
638
639 !###############################################################################
640 subroutine QE_geom_write(UNIT, basis, cartesian)
641 !! Write the structure in Quantum Espresso style format.
642 implicit none
643
644 ! Arguments
645 integer, intent(in) :: UNIT
646 !! The unit number of the file to write to.
647 class(basis_type), intent(in) :: basis
648 !! The basis to write the geometry from.
649 logical, intent(in), optional :: cartesian
650 !! Optional. Whether to write the basis in cartesian coordinates.
651
652 ! Local variables
653 integer :: i,j
654 !! Loop index.
655 character(10) :: string
656 !! String to determine whether to write in crystal or angstrom coordinates.
657
658
659 string="crystal"
660 if(present(cartesian))then
661 if(cartesian) string="angstrom"
662 end if
663
664
665 write(UNIT,'("CELL_PARAMETERS angstrom")')
666 do i = 1, 3
667 write(UNIT,'(3(F15.9))') basis%lat(i,:)
668 end do
669 write(UNIT,'("ATOMIC_SPECIES")')
670 do i = 1, basis%nspec
671 write(UNIT,'(A)') trim(adjustl(basis%spec(i)%name))
672 end do
673 write(UNIT,'("ATOMIC_POSITIONS",1X,A)') trim(adjustl(string))
674 do i = 1, basis%nspec
675 do j = 1, basis%spec(i)%num
676 write(UNIT,'(A5,1X,3(F15.9))') &
677 basis%spec(i)%name,basis%spec(i)%atom(j,1:3)
678 end do
679 end do
680
681 end subroutine QE_geom_write
682 !###############################################################################
683
684
685 !###############################################################################
686 subroutine CASTEP_geom_read(UNIT, basis, length)
687 !! Read the structure in CASTEP style format.
688 implicit none
689
690 ! Arguments
691 integer, intent(in) :: UNIT
692 !! The unit number of the file to read from.
693 type(basis_type), intent(out) :: basis
694 !! The basis to read the geometry into.
695 integer, intent(in), optional :: length
696 !! Optional. The dimension of the basis atom positions.
697
698 ! Local variables
699 integer :: Reason
700 !! The I/O status of the read.
701 integer :: i, j, k, iline
702 !! Loop index.
703 integer :: length_ = 3
704 !! The dimension of the basis atom positions.
705 integer :: itmp1
706 !! Temporary integer variable.
707 character(len=3) :: ctmp
708 !! Temporary character variable.
709 character(len=20) :: units
710 !! Units of the lattice vectors.
711 character(len=200) :: buffer, store
712 !! Temporary character variables.
713 logical :: labc
714 !! Logical variable to determine whether the lattice is in abc or
715 !! cartesian coordinates.
716 integer, dimension(1000) :: tmp_natom
717 !! Temporary array to store the number of atoms of each species.
718 real(real32), dimension(3) :: abc, angle, dvtmp1
719 !! Temporary arrays to store the lattice vectors.
720 character(len=3), dimension(1000) :: tmp_spec
721 !! Temporary array to store the species names.
722
723
724 !---------------------------------------------------------------------------
725 ! determine dimension of basis (include translation dimension for symmetry?)
726 !---------------------------------------------------------------------------
727 if(present(length)) length_ = length
728
729
730 !---------------------------------------------------------------------------
731 ! reading loop of file
732 !---------------------------------------------------------------------------
733 tmp_spec = ""
734 tmp_natom = 0
735 iline = 0
736 labc = .true.
737 basis%sysname = "from CASTEP"
738 rewind(UNIT)
739 readloop: do
740 iline=iline+1
741 read(UNIT,'(A)',iostat=Reason) buffer
742 if(Reason.ne.0) exit
743 buffer=to_upper(buffer)
744 if(scan(trim(adjustl(buffer)),'%').ne.1) cycle readloop
745 if(index(trim(adjustl(buffer)),'%END').eq.1) cycle readloop
746 read(buffer,*) store, buffer
747 if(trim(buffer).eq.'') cycle readloop
748 !------------------------------------------------------------------------
749 ! read lattice
750 !------------------------------------------------------------------------
751 lattice_if: if(index(trim(buffer),"LATTICE").eq.1)then
752 if(index(trim(buffer),"ABC").ne.0) labc = .true.
753 if(index(trim(buffer),"CART").ne.0) labc = .false.
754 store = ""
755 itmp1 = 0
756 lattice_loop: do
757 itmp1 = itmp1 + 1
758 read(UNIT,'(A)',iostat=Reason) buffer
759 if(Reason.ne.0) exit lattice_loop
760 if(scan(trim(adjustl(buffer)),'%').eq.1) exit lattice_loop
761 if(itmp1.eq.5)then
762 call stop_program( &
763 "Too many lines in LATTICE block of structure file" &
764 )
765 return
766 end if
767 store=trim(store)//" "//trim(buffer)
768 end do lattice_loop
769 iline=iline+itmp1
770
771 if(labc)then
772 read(store,*) units,(abc(i),i=1,3), (angle(j),j=1,3)
773 basis%lat = convert_abc_to_lat(abc,angle,.false.)
774 else
775 read(store,*) units,(basis%lat(i,:),i=1,3)
776 end if
777 cycle readloop
778 end if lattice_if
779
780 !------------------------------------------------------------------------
781 ! read basis
782 !------------------------------------------------------------------------
783 basis_if: if(index(trim(buffer),"POSITIONS").eq.1) then
784 if(index(trim(buffer),"ABS").ne.0) basis%lcart=.true.
785 if(index(trim(buffer),"FRAC").ne.0) basis%lcart=.false.
786 itmp1 = 0
787 basis_loop1: do
788 read(UNIT,'(A)',iostat=Reason) buffer
789 if(Reason.ne.0) exit basis_loop1
790 if(scan(trim(adjustl(buffer)),'%').eq.1) exit basis_loop1
791 read(buffer,*) ctmp
792 if(trim(ctmp).eq.'') exit
793 if(verify(buffer,' 0123456789').eq.0) exit
794 basis%natom = basis%natom + 1
795 if(.not.any(tmp_spec(1:basis%nspec).eq.ctmp))then
796 basis%nspec = basis%nspec+1
797 tmp_natom(basis%nspec) = 1
798 tmp_spec(basis%nspec) = ctmp
799 else
800 where(tmp_spec(1:basis%nspec).eq.ctmp)
801 tmp_natom(1:basis%nspec) = tmp_natom(1:basis%nspec) + 1
802 end where
803 end if
804 end do basis_loop1
805
806 allocate(basis%spec(basis%nspec))
807 basis%spec(1:basis%nspec)%name = tmp_spec(1:basis%nspec)
808 do i = 1, basis%nspec
809 basis%spec(i)%num = 0
810 allocate(basis%spec(i)%atom(tmp_natom(i),length_))
811 end do
812
813 call jump(UNIT,iline)
814 basis_loop2: do i = 1, basis%natom
815 read(UNIT,'(A)',iostat=Reason) buffer
816 if(Reason.ne.0)then
817 call stop_program("Internal error in assigning the basis")
818 return
819 end if
820 read(buffer,*) ctmp,dvtmp1(1:3)
821 species_loop: do j = 1, basis%nspec
822 if(basis%spec(j)%name.eq.ctmp)then
823 basis%spec(j)%num = basis%spec(j)%num + 1
824 basis%spec(j)%atom(basis%spec(j)%num,1:3) = dvtmp1(1:3)
825 exit species_loop
826 end if
827 end do species_loop
828 end do basis_loop2
829
830 end if basis_if
831 end do readloop
832
833
834 !---------------------------------------------------------------------------
835 ! convert basis if in cartesian coordinates
836 !---------------------------------------------------------------------------
837 if(basis%lcart) call basis%convert()
838
839
840 !---------------------------------------------------------------------------
841 ! normalise basis to between 0 and 1 in direct coordinates
842 !---------------------------------------------------------------------------
843 do i = 1, basis%nspec
844 do j = 1, basis%spec(i)%num
845 do k = 1, 3
846 basis%spec(i)%atom(j,k) = &
847 basis%spec(i)%atom(j,k) - floor( basis%spec(i)%atom(j,k) )
848 end do
849 end do
850 end do
851 basis%natom=sum(basis%spec(:)%num)
852
853 end subroutine CASTEP_geom_read
854 !###############################################################################
855
856
857 !###############################################################################
858 subroutine CASTEP_geom_write(UNIT, basis, labc, cartesian)
859 !! Write the structure in CASTEP style format.
860 implicit none
861
862 ! Arguments
863 integer :: UNIT
864 !! The unit number of the file to write to.
865 class(basis_type), intent(in) :: basis
866 !! The basis to write the geometry from.
867 logical, intent(in), optional :: labc
868 !! Optional. Boolean whether to write the lattice in abc format.
869 logical, intent(in), optional :: cartesian
870 !! Optional. Boolean whether to write basis in cartesian coordinates.
871
872 ! Local variables
873 integer :: i, j
874 !! Loop index.
875 real(real32), dimension(2,3) :: abc_angle
876 !! Temporary arrays to store the lattice vectors.
877 character(4) :: string_lat, string_bas
878 !! Strings specifying lattice and basis format
879 character(len=256) :: stop_msg
880 !! Error message.
881
882
883 string_lat="CART"
884 if(present(labc))then
885 if(labc) string_lat="ABC"
886 end if
887
888 string_bas="FRAC"
889 if(present(cartesian))then
890 if(cartesian)then
891 string_bas="ABS"
892 write(stop_msg,*) &
893 "Internal error in CASTEP_geom_write" // &
894 achar(13) // achar(10) // &
895 " Subroutine not yet set up to output cartesian coordinates"
896 call stop_program(stop_msg)
897 return
898 end if
899 end if
900
901 write(UNIT,'("%block LATTICE_",A)') trim(string_lat)
902 write(UNIT,'("ang")')
903 if(present(labc))then
904 if(labc)then
905 abc_angle = convert_lat_to_abc(basis%lat)
906 write(UNIT,'(3(F15.9))') abc_angle(1,:)
907 write(UNIT,'(3(F15.9))') abc_angle(2,:)
908 goto 10
909 end if
910 end if
911 do i = 1, 3
912 write(UNIT,'(3(F15.9))') basis%lat(i,:)
913 end do
914
915 10 write(UNIT,'("%endblock LATTICE_",A)') trim(string_lat)
916
917 write(UNIT,*)
918 write(UNIT,'("%block POSITIONS_",A)') trim(string_bas)
919 do i = 1, basis%nspec
920 do j = 1, basis%spec(i)%num
921 write(UNIT,'(A5,1X,3(F15.9))') &
922 basis%spec(i)%name,basis%spec(i)%atom(j,1:3)
923 end do
924 end do
925 write(UNIT,'("%endblock POSITIONS_",A)') trim(string_bas)
926
927 end subroutine CASTEP_geom_write
928 !###############################################################################
929
930
931 !###############################################################################
932 subroutine XYZ_geom_read(UNIT, basis, length, iostat)
933 !! Read the structure in xyz style format.
934 implicit none
935
936 ! Arguments
937 integer, intent(in) :: UNIT
938 !! The unit number of the file to read from.
939 type(basis_type), intent(out) :: basis
940 !! The basis to read the geometry into.
941 integer, intent(in), optional :: length
942 !! Optional. The dimension of the basis atom positions.
943 integer, intent(out), optional :: iostat
944 !! Optional. The I/O status of the read.
945
946 ! Local variables
947 integer :: Reason
948 !! The I/O status of the read.
949 integer :: i, j
950 !! Loop index.
951 integer, allocatable, dimension(:) :: tmp_num
952 !! Temporary array to store the number of atoms of each species.
953 real(real32), dimension(3) :: vec
954 !! Temporary array to store the atomic positions.
955 real(real32), allocatable, dimension(:,:,:) :: tmp_bas
956 !! Temporary array to store the atomic positions.
957 character(len=3) :: ctmp
958 !! Temporary character variable.
959 character(len=3), allocatable, dimension(:) :: tmp_spec
960 !! Temporary array to store the species names.
961 integer :: length_
962 !! The dimension of the basis atom positions.
963 integer :: iostat_
964 !! The I/O status of the read.
965
966
967 length_ = 3
968 iostat_ = 0
969 if(present(length)) length_ = length
970
971
972 read(UNIT,*,iostat=Reason) basis%natom
973 if(Reason.ne.0)then
974 write(0,'("ERROR: The file is not in xyz format.")')
975 iostat_ = 1
976 if(present(iostat)) iostat = iostat_
977 return
978 end if
979 read(UNIT,'(A)',iostat=Reason) basis%sysname
980
981
982 !---------------------------------------------------------------------------
983 ! read basis
984 !---------------------------------------------------------------------------
985 allocate(tmp_spec(basis%natom))
986 allocate(tmp_num(basis%natom))
987 allocate(tmp_bas(basis%natom,basis%natom,length_))
988 tmp_num(:) = 0
989 tmp_spec = ""
990 tmp_bas = 0
991 basis%nspec = 0
992 do i = 1, basis%natom
993 read(UNIT,*,iostat=Reason) ctmp,vec(1:3)
994 if(.not.any(tmp_spec(1:basis%nspec).eq.ctmp))then
995 basis%nspec = basis%nspec + 1
996 tmp_spec(basis%nspec) = ctmp
997 tmp_bas(basis%nspec,1,1:3) = vec(1:3)
998 tmp_num(basis%nspec) = 1
999 else
1000 checkspec: do j = 1, basis%nspec
1001 if(tmp_spec(j).eq.ctmp)then
1002 tmp_num(j) = tmp_num(j)+1
1003 tmp_bas(j,tmp_num(j),1:3) = vec(1:3)
1004 exit checkspec
1005 end if
1006 end do checkspec
1007 end if
1008 end do
1009
1010
1011 !---------------------------------------------------------------------------
1012 ! move basis from temporary basis to main basis.
1013 ! done to allow for correct allocation of number of and per species
1014 !---------------------------------------------------------------------------
1015 allocate(basis%spec(basis%nspec))
1016 do i = 1, basis%nspec
1017 basis%spec(i)%name = tmp_spec(i)
1018 basis%spec(i)%num = tmp_num(i)
1019 allocate(basis%spec(i)%atom(tmp_num(i),length_))
1020 basis%spec(i)%atom(:,:) = 0
1021 basis%spec(i)%atom(1:tmp_num(i),1:3) = tmp_bas(i,1:tmp_num(i),1:3)
1022 end do
1023
1024 if(present(iostat)) iostat = iostat_
1025
1026 end subroutine XYZ_geom_read
1027 !###############################################################################
1028
1029
1030 !###############################################################################
1031 subroutine XYZ_geom_write(UNIT,basis)
1032 !! Write the structure in xyz style format.
1033 implicit none
1034
1035 ! Arguments
1036 integer, intent(in) :: UNIT
1037 !! The unit number of the file to write to.
1038 class(basis_type), intent(in) :: basis
1039 !! The basis to write the geometry from.
1040
1041 ! Local variables
1042 integer :: i, j
1043 !! Loop index.
1044
1045
1046 write(UNIT,'("I0")') basis%natom
1047 write(UNIT,'("A")') basis%sysname
1048 do i = 1, basis%nspec
1049 do j = 1, basis%spec(i)%num
1050 write(UNIT,'(A5,1X,3(F15.9))') &
1051 basis%spec(i)%name,basis%spec(i)%atom(j,1:3)
1052 end do
1053 end do
1054
1055 end subroutine XYZ_geom_write
1056 !###############################################################################
1057
1058
1059 !###############################################################################
1060
6/16
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 9 times.
✓ Branch 13 taken 3 times.
✓ Branch 14 taken 3 times.
✓ Branch 15 taken 1 times.
16 subroutine extXYZ_geom_read(UNIT, basis, length, iostat)
1061 !! Read the structure in extended xyz style format.
1062 implicit none
1063
1064 ! Arguments
1065 integer, intent(in) :: UNIT
1066 !! The unit number of the file to read from.
1067 type(basis_type), intent(out) :: basis
1068 !! The basis to read the geometry into.
1069 integer, intent(in), optional :: length
1070 !! Optional. The dimension of the basis atom positions.
1071 integer, intent(out), optional :: iostat
1072 !! Optional. The I/O status of the read.
1073
1074 ! Local variables
1075 integer :: Reason
1076 !! The I/O status of the read.
1077 integer :: i, j
1078 !! Loop index.
1079 integer :: index1, index2
1080 !! Index variables.
1081 1 integer, allocatable, dimension(:) :: tmp_num
1082 !! Temporary array to store the number of atoms of each species.
1083 real(real32), dimension(3) :: vec
1084 !! Temporary array to store the atomic positions.
1085 1 real(real32), allocatable, dimension(:,:,:) :: tmp_bas
1086 !! Temporary array to store the atomic positions.
1087 character(len=3) :: ctmp
1088 !! Temporary character variable.
1089
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 character(len=3), allocatable, dimension(:) :: tmp_spec
1090 !! Temporary array to store the species names.
1091 character(len=1024) :: buffer
1092 !! Temporary character variable.
1093 integer :: length_ = 3
1094 !! The dimension of the basis atom positions.
1095 integer :: iostat_ = 0
1096 !! The I/O status of the read.
1097
1098
1099 1 basis%lcart=.true.
1100
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(present(length)) length_ = length
1101
1102
1103 !---------------------------------------------------------------------------
1104 ! read system information
1105 !---------------------------------------------------------------------------
1106 1 read(UNIT,*,iostat=Reason) basis%natom
1107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(Reason.ne.0)then
1108 write(0,'("ERROR: The file is not in xyz format.")')
1109 iostat_ = 1
1110 if(present(iostat)) iostat = iostat_
1111 return
1112 end if
1113 1 read(UNIT,'(A)',iostat=Reason) buffer
1114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(Reason.ne.0)then
1115 write(0,'("ERROR: The file is not in xyz format.")')
1116 iostat_ = 1
1117 if(present(iostat)) iostat = iostat_
1118 return
1119 end if
1120 1 index1 = index(buffer,'Lattice="') + 9
1121 1 index2 = index(buffer(index1:),'"') + index1 - 2
1122
6/8
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 9 times.
✓ Branch 8 taken 3 times.
13 read(buffer(index1:index2),*) ( ( basis%lat(i,j), j = 1, 3), i = 1, 3)
1123
1124 1 index1 = index(buffer,'free_energy=') + 12
1125 1 read(buffer(index1:),*) basis%energy
1126
1127
1128 !---------------------------------------------------------------------------
1129 ! read basis
1130 !---------------------------------------------------------------------------
1131
7/14
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
1 allocate(tmp_spec(basis%natom))
1132
7/14
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
1 allocate(tmp_num(basis%natom))
1133
11/22
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
1 allocate(tmp_bas(basis%natom,basis%natom,length_))
1134
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
9 tmp_num(:) = 0
1135
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
9 tmp_spec = ""
1136
6/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 192 times.
✓ Branch 5 taken 24 times.
220 tmp_bas = 0
1137 1 basis%nspec = 0
1138
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
9 do i = 1, basis%natom
1139 8 read(UNIT,*,iostat=Reason) ctmp, vec(1:3)
1140
5/6
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 7 times.
9 if(.not.any(tmp_spec(1:basis%nspec).eq.ctmp))then
1141 1 basis%nspec=basis%nspec+1
1142
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
1 tmp_spec(basis%nspec) = trim(adjustl(ctmp))
1143
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
4 tmp_bas(basis%nspec,1,1:3) = vec(1:3)
1144 1 tmp_num(basis%nspec) = 1
1145 else
1146
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
7 checkspec: do j = 1, basis%nspec
1147
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
7 if(tmp_spec(j).eq.ctmp)then
1148 7 tmp_num(j) = tmp_num(j) + 1
1149
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 7 times.
28 tmp_bas(j,tmp_num(j),1:3) = vec(1:3)
1150 7 exit checkspec
1151 end if
1152 end do checkspec
1153 end if
1154 end do
1155
1156
1157 !---------------------------------------------------------------------------
1158 ! move basis from temporary basis to main basis.
1159 ! done to allow for correct allocation of number of and per species
1160 !---------------------------------------------------------------------------
1161
13/24
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 1 times.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 1 times.
2 allocate(basis%spec(basis%nspec))
1162 1 basis%sysname = ""
1163
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
2 do i = 1, basis%nspec
1164 1 basis%spec(i)%name = tmp_spec(i)
1165 1 basis%spec(i)%num = tmp_num(i)
1166
9/18
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
1 allocate(basis%spec(i)%atom(tmp_num(i),length_))
1167
4/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 3 times.
28 basis%spec(i)%atom(:,:) = 0
1168
4/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 3 times.
28 basis%spec(i)%atom(1:tmp_num(i),1:3) = tmp_bas(i,1:tmp_num(i),1:3)
1169
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
1 write(buffer,'(I0,A)') basis%spec(i)%num,trim(basis%spec(i)%name)
1170
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
1 basis%sysname = basis%sysname//trim(buffer)
1171
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
2 if(i.lt.basis%nspec) basis%sysname = trim(adjustl(basis%sysname))//"_"
1172 end do
1173
1174
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(present(iostat)) iostat = iostat_
1175
1176
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
1 end subroutine extXYZ_geom_read
1177 !###############################################################################
1178
1179
1180 !###############################################################################
1181 1 subroutine extXYZ_geom_write(UNIT, basis)
1182 !! Write the structure in extended xyz style format.
1183 implicit none
1184
1185 ! Arguments
1186 integer, intent(in) :: UNIT
1187 !! The unit number of the file to write to.
1188 class(basis_type), intent(in) :: basis
1189 !! The basis to write the geometry from.
1190
1191 ! Local variables
1192 integer :: i, j
1193 !! Loop index.
1194
1195
1196 1 write(UNIT,'(I0)') basis%natom
1197 write(UNIT,'(A,8(F0.8,1X),F0.8,A)', advance="no") &
1198
6/8
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 12 times.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 3 times.
13 'Lattice="',((basis%lat(i,j),j=1,3),i=1,3),'"'
1199 1 write(UNIT,'(A,F0.8)', advance="no") ' free_energy=',basis%energy
1200 1 write(UNIT,'(A)', advance="no") ' pbc="T T T"'
1201
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(basis%lcart)then
1202 do i = 1, basis%nspec
1203 do j = 1, basis%spec(i)%num
1204 write(UNIT,'(A8,3(1X, F16.8))') &
1205 basis%spec(i)%name,basis%spec(i)%atom(j,1:3)
1206 end do
1207 end do
1208 else
1209
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
2 do i = 1, basis%nspec
1210
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
10 do j = 1, basis%spec(i)%num
1211 8 write(UNIT,'(A8,3(1X, F16.8))') basis%spec(i)%name, &
1212
2/2
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 8 times.
41 matmul(basis%spec(i)%atom(j,1:3),basis%lat)
1213 end do
1214 end do
1215 end if
1216
1217 1 end subroutine extXYZ_geom_write
1218 !###############################################################################
1219
1220
1221 !###############################################################################
1222 1 subroutine convert(this)
1223 !! Convert the basis between direct and cartesian coordinates.
1224 implicit none
1225
1226 ! Arguments
1227 class(basis_type), intent(inout) :: this
1228 !! Parent. The basis to convert.
1229
1230 ! Local variables
1231 integer :: is, ia
1232 !! Loop index.
1233 real(real32), dimension(3,3) :: lattice
1234 !! The reciprocal lattice vectors.
1235
1236
1237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(this%lcart)then
1238 lattice = inverse_3x3( this%lat )
1239 else
1240
4/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
13 lattice = this%lat
1241 end if
1242
1243 1 this%lcart = .not.this%lcart
1244
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
2 do is = 1, this%nspec
1245
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
10 do ia = 1, this%spec(is)%num
1246 this%spec(is)%atom(ia,1:3) = &
1247
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 8 times.
33 matmul( this%spec(is)%atom(ia,1:3), lattice )
1248 end do
1249 end do
1250
1251 1 end subroutine convert
1252 !###############################################################################
1253
1254
1255 !#############################################################################
1256 function convert_abc_to_lat(abc,angle,radians) result(lattice)
1257 !! Convert the lattice from abc and αβγ to lattice matrix.
1258 implicit none
1259
1260 ! Arguments
1261 real(real32), dimension(3), intent(in) :: abc, angle
1262 !! lattice constants
1263 logical, intent(in), optional :: radians
1264 !! Optional. Boolean whether angles are in radians.
1265 real(real32), dimension(3,3) :: lattice
1266 !! The lattice matrix.
1267
1268 ! Local variables
1269 real(real32), dimension(3) :: in_angle
1270 !! The lattice angles in radians.
1271
1272
1273
1274 in_angle = angle
1275 if(present(radians))then
1276 if(.not.radians) in_angle = angle*pi/180._real32
1277 end if
1278
1279 lattice=0._real32
1280
1281 lattice(1,1)=abc(1)
1282 lattice(2,:2)=(/abc(2)*cos(in_angle(3)),abc(2)*sin(in_angle(3))/)
1283
1284 lattice(3,1) = abc(3)*cos(in_angle(2))
1285 lattice(3,2) = abc(3)*(cos(in_angle(1)) - cos(in_angle(2))*&
1286 cos(in_angle(3)))/sin(in_angle(3))
1287 lattice(3,3) = sqrt(abc(3)**2._real32 - &
1288 lattice(3,1)**2._real32 - &
1289 lattice(3,2)**2._real32)
1290
1291 end function convert_abc_to_lat
1292 !###############################################################################
1293
1294
1295 !###############################################################################
1296 function convert_lat_to_abc(lattice, radians) result(abc_angle)
1297 !! Convert the lattice from lattice matrix to abc and αβγ.
1298 implicit none
1299
1300 ! Arguments
1301 real(real32), dimension(3,3), intent(in) :: lattice
1302 !! The lattice matrix.
1303 logical, intent(in), optional :: radians
1304 !! Optional. Boolean whether to return angles in radians.
1305 real(real32), dimension(2,3) :: abc_angle
1306 !! The lattice constants and angles.
1307
1308 ! Local variables
1309 integer :: i
1310 !! Loop index.
1311
1312
1313 do i = 1, 3
1314 abc_angle(1,i)=norm2(lattice(i,:))
1315 end do
1316 do i = 1, 3
1317 end do
1318 abc_angle(2,1)=acos(dot_product(lattice(2,:),lattice(3,:))/&
1319 (abc_angle(1,2)*abc_angle(1,3)))
1320 abc_angle(2,3)=acos(dot_product(lattice(1,:),lattice(3,:))/&
1321 (abc_angle(1,1)*abc_angle(1,3)))
1322 abc_angle(2,3)=acos(dot_product(lattice(1,:),lattice(2,:))/&
1323 (abc_angle(1,1)*abc_angle(1,2)))
1324
1325 if(present(radians))then
1326 if(.not.radians) abc_angle(2,:)=abc_angle(2,:)*180._real32/pi
1327 end if
1328
1329 end function convert_lat_to_abc
1330 !###############################################################################
1331
1332
1333 !###############################################################################
1334 function get_lattice_constants(this, radians) result(output)
1335 !! Convert the lattice from lattice matrix to abc and αβγ.
1336 implicit none
1337
1338 ! Arguments
1339 class(basis_type), intent(in) :: this
1340 !! Parent. The basis.
1341 logical, intent(in), optional :: radians
1342 !! Optional. Boolean whether to return angles in radians.
1343 real(real32), dimension(2,3) :: output
1344 !! The lattice constants and angles.
1345
1346 ! Local variables
1347 logical :: radians_
1348 !! Boolean whether to return angles in radians.
1349
1350
1351 radians_ = .true.
1352 if(present(radians)) radians_ = radians
1353
1354 output = convert_lat_to_abc(this%lat, radians_)
1355
1356 end function get_lattice_constants
1357 !###############################################################################
1358
1359
1360 !###############################################################################
1361 78 subroutine copy(this, basis, length)
1362 !! Copy the basis.
1363 implicit none
1364
1365 ! Arguments
1366 class(basis_type), intent(inout) :: this
1367 !! Parent. The basis to copy into.
1368 class(basis_type), intent(in) :: basis
1369 !! The basis to copy from.
1370 integer, intent(in), optional :: length
1371 !! The dimension of the basis atom positions.
1372
1373
1374 ! Local variables
1375 integer :: i, j
1376 !! Loop indices.
1377 integer :: length_, length_input
1378 !! The dimension of the basis atom positions.
1379
1380
1381 !---------------------------------------------------------------------------
1382 ! determines whether user wants output basis extra translational dimension
1383 !---------------------------------------------------------------------------
1384 78 length_input = size(basis%spec(1)%atom(1,:),dim=1)
1385
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
78 if(present(length))then
1386 length_ = length
1387 else
1388 78 length_ = length_input
1389 end if
1390
1391
1392 !---------------------------------------------------------------------------
1393 ! if already allocated, deallocates output basis
1394 !---------------------------------------------------------------------------
1395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
78 if(allocated(this%spec))then
1396 do i = 1, this%nspec
1397 if(allocated(this%spec(i)%atom_mask)) &
1398 deallocate(this%spec(i)%atom_mask)
1399 if(allocated(this%spec(i)%atom_idx)) deallocate(this%spec(i)%atom_idx)
1400 if(allocated(this%spec(i)%atom)) deallocate(this%spec(i)%atom)
1401 end do
1402 deallocate(this%spec)
1403 end if
1404
1405
1406 !---------------------------------------------------------------------------
1407 ! allocates output basis and clones data from input basis to output basis
1408 !---------------------------------------------------------------------------
1409
13/24
✓ Branch 0 taken 78 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 78 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 78 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 78 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 78 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 78 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 78 times.
✓ Branch 17 taken 95 times.
✓ Branch 18 taken 78 times.
✓ Branch 19 taken 95 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 95 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 95 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 95 times.
173 allocate(this%spec(basis%nspec))
1410
2/2
✓ Branch 0 taken 95 times.
✓ Branch 1 taken 78 times.
173 do i = 1, basis%nspec
1411
9/16
✓ Branch 0 taken 95 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 95 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 95 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 95 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 95 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 95 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 95 times.
✓ Branch 17 taken 609 times.
✓ Branch 18 taken 95 times.
704 allocate(this%spec(i)%atom_mask(basis%spec(i)%num), source = .true.)
1412
7/14
✓ Branch 0 taken 95 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 95 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 95 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 95 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 95 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 95 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 95 times.
95 allocate(this%spec(i)%atom_idx(basis%spec(i)%num))
1413
9/18
✓ Branch 0 taken 95 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 95 times.
✓ Branch 4 taken 95 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 95 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 95 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 95 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 95 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 95 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 95 times.
95 allocate(this%spec(i)%atom(basis%spec(i)%num,length_))
1414
1415
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 40 times.
95 if(allocated(basis%spec(i)%atom_mask)) &
1416
4/12
✓ Branch 0 taken 55 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 55 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 402 times.
✓ Branch 11 taken 55 times.
512 this%spec(i)%atom_mask = basis%spec(i)%atom_mask
1417
1418
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 44 times.
95 if(allocated(basis%spec(i)%atom_idx))then
1419
4/12
✓ Branch 0 taken 51 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 376 times.
✓ Branch 11 taken 51 times.
478 this%spec(i)%atom_idx = basis%spec(i)%atom_idx
1420 else
1421 this%spec(i)%atom_idx = [ ( j, j = sum(basis%spec(1:i-1:1)%num) + 1, &
1422
14/20
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 53 times.
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 44 times.
✓ Branch 7 taken 233 times.
✓ Branch 8 taken 44 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 44 times.
✓ Branch 11 taken 233 times.
✓ Branch 12 taken 44 times.
✓ Branch 13 taken 44 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 44 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 233 times.
✓ Branch 20 taken 44 times.
805 sum(basis%spec(1:i)%num) ) ]
1423 end if
1424
1/2
✓ Branch 0 taken 95 times.
✗ Branch 1 not taken.
95 if(length_input.eq.length_)then
1425
4/4
✓ Branch 0 taken 285 times.
✓ Branch 1 taken 95 times.
✓ Branch 2 taken 1827 times.
✓ Branch 3 taken 285 times.
2207 this%spec(i)%atom(:,:length_) = basis%spec(i)%atom(:,:length_)
1426 elseif(length_input.gt.length_)then
1427 this%spec(i)%atom(:,:3) = basis%spec(i)%atom(:,:3)
1428 elseif(length_input.lt.length_)then
1429 this%spec(i)%atom(:,:3) = basis%spec(i)%atom(:,:3)
1430 this%spec(i)%atom(:,4) = 1._real32
1431 end if
1432 95 this%spec(i)%num = basis%spec(i)%num
1433 95 this%spec(i)%name = strip_null(basis%spec(i)%name)
1434
1435 95 this%spec(i)%mass = basis%spec(i)%mass
1436 95 this%spec(i)%charge = basis%spec(i)%charge
1437 173 this%spec(i)%radius = basis%spec(i)%radius
1438 end do
1439 78 this%nspec = basis%nspec
1440 78 this%natom = basis%natom
1441 78 this%lcart = basis%lcart
1442 78 this%sysname = basis%sysname
1443 78 this%energy = basis%energy
1444
4/4
✓ Branch 0 taken 234 times.
✓ Branch 1 taken 78 times.
✓ Branch 2 taken 702 times.
✓ Branch 3 taken 234 times.
1014 this%lat = basis%lat
1445
2/2
✓ Branch 0 taken 234 times.
✓ Branch 1 taken 78 times.
312 this%pbc = basis%pbc
1446
1447 78 end subroutine copy
1448 !###############################################################################
1449
1450
1451 !###############################################################################
1452
4/6
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
22 subroutine set_atom_mask(this, index_list)
1453 !! Set the mask for the atoms in the basis.
1454 implicit none
1455
1456 ! Arguments
1457 class(basis_type), intent(inout) :: this
1458 !! Parent. The basis.
1459 integer, dimension(:,:), intent(in), optional :: index_list
1460 !! The list of indices to set the mask for.
1461
1462 ! Local variables
1463 integer :: i
1464 !! Loop index.
1465
1466
1467
2/2
✓ Branch 0 taken 35 times.
✓ Branch 1 taken 22 times.
57 do i = 1, this%nspec
1468
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 28 times.
57 if(.not.allocated(this%spec(i)%atom_mask))then
1469 allocate( &
1470 this%spec(i)%atom_mask(this%spec(i)%num), source = .true. &
1471
9/16
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 7 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 7 times.
✓ Branch 16 taken 30 times.
✓ Branch 17 taken 7 times.
37 )
1472 end if
1473 end do
1474
1475
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 2 times.
22 if(present(index_list))then
1476
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 20 times.
51 do i = 1, size(index_list,2)
1477 this%spec(index_list(1,i))%atom_mask(index_list(2,i)) = &
1478 51 .not.this%spec(index_list(1,i))%atom_mask(index_list(2,i))
1479 end do
1480 end if
1481
1482 22 end subroutine set_atom_mask
1483 !###############################################################################
1484
1485
1486 !###############################################################################
1487 subroutine add_atom(this, species, position, is_cartesian, mask)
1488 !! Add an atom to the basis.
1489 implicit none
1490
1491 ! Arguments
1492 class(basis_type), intent(inout) :: this
1493 !! Parent. The basis.
1494 character(len=3), intent(in) :: species
1495 !! The species of the atom to add.
1496 real(real32), dimension(3), intent(in) :: position
1497 !! The position of the atom to add.
1498 logical, intent(in), optional :: is_cartesian
1499 !! Optional. Boolean whether the position is in cartesian coordinates.
1500 !! NOT YET IMPLEMENTED.
1501 logical, intent(in), optional :: mask
1502 !! Optional. Boolean whether to add a mask for the atom.
1503
1504 ! Local variables
1505 integer :: j
1506 !! Loop index.
1507 integer :: idx
1508 !! The index of the species in the basis.
1509 integer :: length
1510 !! The dimension of the basis atom positions.
1511 logical :: mask_
1512 !! Boolean mask for the atom.
1513 integer, dimension(:), allocatable :: atom_idx
1514 !! Temporary array.
1515 logical, dimension(:), allocatable :: atom_mask
1516 !! Temporary array.
1517 real(real32), dimension(:,:), allocatable :: positions
1518 !! Temporary array.
1519 type(species_type), dimension(:), allocatable :: species_list
1520 !! Temporary array.
1521
1522
1523 mask_ = .true.
1524 if(present(mask)) mask_ = mask
1525
1526 this%natom = this%natom + 1
1527 length = size(this%spec(1)%atom,dim=2)
1528 idx = findloc(this%spec(:)%name, strip_null(species), dim=1)
1529 if(idx.eq.0)then
1530 this%nspec = this%nspec + 1
1531 allocate(species_list(this%nspec))
1532 species_list(1:this%nspec-1) = this%spec(1:this%nspec-1)
1533 deallocate(this%spec)
1534 species_list(this%nspec)%name = strip_null(species)
1535 species_list(this%nspec)%num = 1
1536 call get_element_properties(species_list(this%nspec)%name, &
1537 species_list(this%nspec)%mass, &
1538 species_list(this%nspec)%charge, &
1539 species_list(this%nspec)%radius &
1540 )
1541 allocate(species_list(this%nspec)%atom_idx(1))
1542 allocate(species_list(this%nspec)%atom_mask(1), source = mask_)
1543 species_list(this%nspec)%atom_idx(1) = this%natom
1544 allocate(species_list(this%nspec)%atom(1,length))
1545 species_list(this%nspec)%atom(1,:) = 0._real32
1546 species_list(this%nspec)%atom(1,:3) = position
1547 this%spec = species_list
1548 deallocate(species_list)
1549 else
1550 allocate(atom_mask(this%spec(idx)%num+1), source = .true.)
1551 if(allocated(this%spec(idx)%atom_mask))then
1552 atom_mask(1:this%spec(idx)%num) = this%spec(idx)%atom_mask
1553 end if
1554 atom_mask(this%spec(idx)%num+1) = mask_
1555 allocate(atom_idx(this%spec(idx)%num+1))
1556 if(allocated(this%spec(idx)%atom_idx))then
1557 atom_idx(1:this%spec(idx)%num) = this%spec(idx)%atom_idx
1558 else
1559 atom_idx(1:this%spec(idx)%num) = [ ( j, j = 1, this%spec(idx)%num ) ]
1560 end if
1561 atom_idx(this%spec(idx)%num+1) = this%natom
1562 allocate(positions(this%spec(idx)%num+1,length))
1563 positions = 0._real32
1564 positions(1:this%spec(idx)%num,:) = this%spec(idx)%atom
1565 positions(this%spec(idx)%num+1,:3) = position
1566 this%spec(idx)%num = this%spec(idx)%num + 1
1567 this%spec(idx)%atom_mask = atom_mask
1568 this%spec(idx)%atom_idx = atom_idx
1569 this%spec(idx)%atom = positions
1570 deallocate(atom_mask)
1571 deallocate(atom_idx)
1572 deallocate(positions)
1573 end if
1574
1575 end subroutine add_atom
1576 !###############################################################################
1577
1578
1579 !###############################################################################
1580 subroutine remove_atom(this, ispec, iatom)
1581 !! Remove an atom from the basis.
1582 implicit none
1583
1584 ! Arguments
1585 class(basis_type), intent(inout) :: this
1586 !! Parent. The basis.
1587 integer, intent(in) :: ispec, iatom
1588 !! The species and atom to remove.
1589
1590 ! Local variables
1591 integer :: i
1592 !! Loop index.
1593 integer :: remove_idx
1594 !! The index associated with the atom to remove.
1595 integer, dimension(:), allocatable :: atom_idx
1596 !! Temporary array to store the atomic indices.
1597 logical, dimension(:), allocatable :: atom_mask
1598 !! Temporary array to store the atomic masks.
1599 real(real32), dimension(:,:), allocatable :: atom
1600 !! Temporary array to store the atomic positions.
1601
1602
1603 !---------------------------------------------------------------------------
1604 ! remove atom from basis
1605 !---------------------------------------------------------------------------
1606 remove_idx = this%spec(ispec)%atom_idx(iatom)
1607 do i = 1, this%nspec
1608 if(i.eq.ispec)then
1609 if(iatom.gt.this%spec(i)%num)then
1610 call stop_program("Atom to remove does not exist")
1611 return
1612 end if
1613 allocate(atom_mask(this%spec(i)%num-1), source = .true.)
1614 allocate(atom_idx(this%spec(i)%num-1))
1615 allocate(atom(this%spec(i)%num-1,size(this%spec(i)%atom,2)))
1616 if(iatom.eq.1)then
1617 atom_mask(1:this%spec(i)%num-1) = &
1618 this%spec(i)%atom_mask(2:this%spec(i)%num:1)
1619 atom_idx(1:this%spec(i)%num-1) = &
1620 this%spec(i)%atom_idx(2:this%spec(i)%num:1)
1621 atom(1:this%spec(i)%num-1:1,:) = &
1622 this%spec(i)%atom(2:this%spec(i)%num:1,:)
1623 elseif(iatom.eq.this%spec(i)%num)then
1624 atom_mask(1:this%spec(i)%num-1) = &
1625 this%spec(i)%atom_mask(1:this%spec(i)%num-1:1)
1626 atom_idx(1:this%spec(i)%num-1) = &
1627 this%spec(i)%atom_idx(1:this%spec(i)%num-1:1)
1628 atom(1:this%spec(i)%num-1:1,:) = &
1629 this%spec(i)%atom(1:this%spec(i)%num-1:1,:)
1630 else
1631 atom_mask(1:iatom-1:1) = this%spec(i)%atom_mask(1:iatom-1:1)
1632 atom_idx(1:iatom-1:1) = this%spec(i)%atom_idx(1:iatom-1:1)
1633 atom_idx(iatom:this%spec(i)%num-1:1) = &
1634 this%spec(i)%atom_idx(iatom+1:this%spec(i)%num:1)
1635 atom(1:iatom-1:1,:) = this%spec(i)%atom(1:iatom-1:1,:)
1636 atom(iatom:this%spec(i)%num-1:1,:) = &
1637 this%spec(i)%atom(iatom+1:this%spec(i)%num:1,:)
1638 end if
1639 where(atom_idx(1:this%spec(i)%num-1:1).gt.remove_idx)
1640 atom_idx(1:this%spec(i)%num-1:1) = &
1641 atom_idx(1:this%spec(i)%num-1:1) - 1
1642 end where
1643 this%spec(i)%atom_mask = atom_mask
1644 this%spec(i)%atom_idx = atom_idx
1645 this%spec(i)%atom = atom
1646 deallocate(atom_mask)
1647 deallocate(atom_idx)
1648 deallocate(atom)
1649 this%spec(i)%num = this%spec(i)%num - 1
1650 this%natom = this%natom - 1
1651 if(this%spec(i)%num.eq.0)then
1652 deallocate(this%spec(i)%atom)
1653 if(this%nspec.eq.0)then
1654 deallocate(this%spec)
1655 this%lcart = .true.
1656 this%sysname = ""
1657 this%energy = 0._real32
1658 this%lat = 0._real32
1659 this%pbc = .true.
1660 end if
1661 end if
1662 end if
1663 end do
1664
1665 end subroutine remove_atom
1666 !###############################################################################
1667
1668
1669 !###############################################################################
1670 subroutine remove_atoms(this, atoms)
1671 !! Remove atoms from the basis.
1672 use raffle__misc, only: swap
1673 implicit none
1674
1675 ! Arguments
1676 class(basis_type), intent(inout) :: this
1677 !! Parent. The basis.
1678 integer, dimension(:,:), intent(in) :: atoms
1679 !! The atoms to remove (2, number of atoms to remove)
1680 !! 1st value of 1st dimension is the species number
1681 !! 2nd value of 1st dimension is the atom number
1682 !! 2nd dimension is the number of atoms to remove
1683
1684 ! Local variables
1685 integer :: is, ia, i
1686 !! Loop index.
1687 integer :: n, m, start_idx, end_idx, loc
1688 !! Index variables.
1689 integer :: num_species
1690 !! The number of species.
1691 integer, dimension(:,:), allocatable :: atoms_ordered
1692 !! The atoms to remove ordered by species and atom
1693 real(real32), dimension(:,:), allocatable :: atom
1694 !! Temporary array to store the atomic positions.
1695
1696
1697 !---------------------------------------------------------------------------
1698 ! reorder atoms to remove
1699 !---------------------------------------------------------------------------
1700 allocate(atoms_ordered, source=atoms)
1701 n = size(atoms_ordered, 1)
1702 m = size(atoms_ordered, 2)
1703
1704 do i = 1, m
1705 loc = maxloc(atoms_ordered(1, i:n), dim=1) + i - 1
1706 if (loc .ne. i) then
1707 call swap(atoms_ordered(1, i), atoms_ordered(1, loc))
1708 call swap(atoms_ordered(2, i), atoms_ordered(2, loc))
1709 end if
1710 end do
1711 num_species = this%nspec
1712 do is = 1, num_species
1713 start_idx = findloc(atoms_ordered(1, :), is, dim=1)
1714 end_idx = findloc(atoms_ordered(1, :), is, dim=1, back=.true.)
1715 if (start_idx .eq. 0) cycle
1716 do ia = start_idx, end_idx, 1
1717 loc = maxloc( &
1718 atoms_ordered(2, ia:end_idx), &
1719 dim=1 &
1720 ) + ia - 1
1721 if (loc .ne. ia) then
1722 call swap(atoms_ordered(1, ia), atoms_ordered(1, loc))
1723 call swap(atoms_ordered(2, ia), atoms_ordered(2, loc))
1724 end if
1725 end do
1726 end do
1727
1728
1729 !---------------------------------------------------------------------------
1730 ! remove atoms from basis
1731 !---------------------------------------------------------------------------
1732 do i = 1, size(atoms_ordered, 2)
1733 call this%remove_atom(atoms_ordered(1, i), atoms_ordered(2, i))
1734 end do
1735
1736 do is = 1, this%nspec
1737 if (this%spec(is)%num .eq. 0) then
1738 this%spec = [ this%spec(1:is-1), this%spec(is+1:) ]
1739 this%nspec = this%nspec - 1
1740 end if
1741 end do
1742
1743 end subroutine remove_atoms
1744 !###############################################################################
1745
1746
1747 !###############################################################################
1748 256 subroutine get_element_properties(element, charge, mass, radius)
1749 !! Set the mass and charge of the element
1750 implicit none
1751
1752 ! Arguments
1753 character(len=3), intent(in) :: element
1754 !! Element name.
1755 real(real32), intent(out), optional :: charge
1756 !! Charge of the element.
1757 real(real32), intent(out), optional :: mass
1758 !! Mass of the element.
1759 real(real32), intent(out), optional :: radius
1760 !! Radius of the element.
1761
1762 ! Local variables
1763 real(real32) :: mass_, charge_, radius_
1764 !! Mass, charge and radius of the element.
1765
1766 2 select case(element)
1767 case('H')
1768 2 mass_ = 1.00784_real32
1769 2 charge_ = 1.0_real32
1770 2 radius_ = 0.31_real32
1771 case('He')
1772 2 mass_ = 4.0026_real32
1773 2 charge_ = 2.0_real32
1774 2 radius_ = 0.28_real32
1775 case('Li')
1776 2 mass_ = 6.94_real32
1777 2 charge_ = 3.0_real32
1778 2 radius_ = 1.28_real32
1779 case('Be')
1780 2 mass_ = 9.0122_real32
1781 2 charge_ = 4.0_real32
1782 2 radius_ = 0.96_real32
1783 case('B')
1784 2 mass_ = 10.81_real32
1785 2 charge_ = 5.0_real32
1786 2 radius_ = 0.84_real32
1787 case('C')
1788 7 mass_ = 12.011_real32
1789 7 charge_ = 6.0_real32
1790 7 radius_ = 0.76_real32
1791 case('N')
1792 2 mass_ = 14.007_real32
1793 2 charge_ = 7.0_real32
1794 2 radius_ = 0.71_real32
1795 case('O')
1796 5 mass_ = 15.999_real32
1797 5 charge_ = 8.0_real32
1798 5 radius_ = 0.66_real32
1799 case('F')
1800 2 mass_ = 18.998_real32
1801 2 charge_ = 9.0_real32
1802 2 radius_ = 0.57_real32
1803 case('Ne')
1804 2 mass_ = 20.180_real32
1805 2 charge_ = 10.0_real32
1806 2 radius_ = 0.58_real32
1807 case('Na')
1808 2 mass_ = 22.989_real32
1809 2 charge_ = 11.0_real32
1810 2 radius_ = 1.66_real32
1811 case('Mg')
1812 3 mass_ = 24.305_real32
1813 3 charge_ = 12.0_real32
1814 3 radius_ = 1.41_real32
1815 case('Al')
1816 2 mass_ = 26.982_real32
1817 2 charge_ = 13.0_real32
1818 2 radius_ = 1.21_real32
1819 case('Si')
1820 6 mass_ = 28.085_real32
1821 6 charge_ = 14.0_real32
1822 6 radius_ = 1.11_real32
1823 case('P')
1824 2 mass_ = 30.974_real32
1825 2 charge_ = 15.0_real32
1826 2 radius_ = 1.07_real32
1827 case('S')
1828 2 mass_ = 32.06_real32
1829 2 charge_ = 16.0_real32
1830 2 radius_ = 1.05_real32
1831 case('Cl')
1832 2 mass_ = 35.453_real32
1833 2 charge_ = 17.0_real32
1834 2 radius_ = 1.02_real32
1835 case('Ar')
1836 2 mass_ = 39.948_real32
1837 2 charge_ = 18.0_real32
1838 2 radius_ = 1.06_real32
1839 case('K')
1840 2 mass_ = 39.098_real32
1841 2 charge_ = 19.0_real32
1842 2 radius_ = 2.03_real32
1843 case('Ca')
1844 2 mass_ = 40.078_real32
1845 2 charge_ = 20.0_real32
1846 2 radius_ = 1.74_real32
1847 case('Sc')
1848 2 mass_ = 44.956_real32
1849 2 charge_ = 21.0_real32
1850 2 radius_ = 1.44_real32
1851 case('Ti')
1852 3 mass_ = 47.867_real32
1853 3 charge_ = 22.0_real32
1854 3 radius_ = 1.32_real32
1855 case('V')
1856 2 mass_ = 50.942_real32
1857 2 charge_ = 23.0_real32
1858 2 radius_ = 1.22_real32
1859 case('Cr')
1860 2 mass_ = 51.996_real32
1861 2 charge_ = 24.0_real32
1862 2 radius_ = 1.18_real32
1863 case('Mn')
1864 2 mass_ = 54.938_real32
1865 2 charge_ = 25.0_real32
1866 2 radius_ = 1.17_real32
1867 case('Fe')
1868 2 mass_ = 55.845_real32
1869 2 charge_ = 26.0_real32
1870 2 radius_ = 1.17_real32
1871 case('Co')
1872 2 mass_ = 58.933_real32
1873 2 charge_ = 27.0_real32
1874 2 radius_ = 1.16_real32
1875 case('Ni')
1876 2 mass_ = 58.693_real32
1877 2 charge_ = 28.0_real32
1878 2 radius_ = 1.15_real32
1879 case('Cu')
1880 2 mass_ = 63.546_real32
1881 2 charge_ = 29.0_real32
1882 2 radius_ = 1.17_real32
1883 case('Zn')
1884 2 mass_ = 65.38_real32
1885 2 charge_ = 30.0_real32
1886 2 radius_ = 1.25_real32
1887 case('Ga')
1888 2 mass_ = 69.723_real32
1889 2 charge_ = 31.0_real32
1890 2 radius_ = 1.26_real32
1891 case('Ge')
1892 2 mass_ = 72.63_real32
1893 2 charge_ = 32.0_real32
1894 2 radius_ = 1.22_real32
1895 case('As')
1896 2 mass_ = 74.922_real32
1897 2 charge_ = 33.0_real32
1898 2 radius_ = 1.19_real32
1899 case('Se')
1900 2 mass_ = 78.971_real32
1901 2 charge_ = 34.0_real32
1902 2 radius_ = 1.16_real32
1903 case('Br')
1904 2 mass_ = 79.904_real32
1905 2 charge_ = 35.0_real32
1906 2 radius_ = 1.14_real32
1907 case('Kr')
1908 2 mass_ = 83.798_real32
1909 2 charge_ = 36.0_real32
1910 2 radius_ = 1.12_real32
1911 case('Rb')
1912 2 mass_ = 85.468_real32
1913 2 charge_ = 37.0_real32
1914 2 radius_ = 2.16_real32
1915 case('Sr')
1916 2 mass_ = 87.62_real32
1917 2 charge_ = 38.0_real32
1918 2 radius_ = 1.91_real32
1919 case('Y')
1920 2 mass_ = 88.906_real32
1921 2 charge_ = 39.0_real32
1922 2 radius_ = 1.62_real32
1923 case('Zr')
1924 2 mass_ = 91.224_real32
1925 2 charge_ = 40.0_real32
1926 2 radius_ = 1.45_real32
1927 case('Nb')
1928 2 mass_ = 92.906_real32
1929 2 charge_ = 41.0_real32
1930 2 radius_ = 1.34_real32
1931 case('Mo')
1932 2 mass_ = 95.95_real32
1933 2 charge_ = 42.0_real32
1934 2 radius_ = 1.3_real32
1935 case('Tc')
1936 2 mass_ = 98.0_real32
1937 2 charge_ = 43.0_real32
1938 2 radius_ = 1.27_real32
1939 case('Ru')
1940 2 mass_ = 101.07_real32
1941 2 charge_ = 44.0_real32
1942 2 radius_ = 1.25_real32
1943 case('Rh')
1944 2 mass_ = 102.91_real32
1945 2 charge_ = 45.0_real32
1946 2 radius_ = 1.25_real32
1947 case('Pd')
1948 2 mass_ = 106.42_real32
1949 2 charge_ = 46.0_real32
1950 2 radius_ = 1.28_real32
1951 case('Ag')
1952 2 mass_ = 107.87_real32
1953 2 charge_ = 47.0_real32
1954 2 radius_ = 1.34_real32
1955 case('Cd')
1956 2 mass_ = 112.41_real32
1957 2 charge_ = 48.0_real32
1958 2 radius_ = 1.48_real32
1959 case('In')
1960 2 mass_ = 114.82_real32
1961 2 charge_ = 49.0_real32
1962 2 radius_ = 1.44_real32
1963 case('Sn')
1964 2 mass_ = 118.71_real32
1965 2 charge_ = 50.0_real32
1966 2 radius_ = 1.41_real32
1967 case('Sb')
1968 2 mass_ = 121.76_real32
1969 2 charge_ = 51.0_real32
1970 2 radius_ = 1.38_real32
1971 case('Te')
1972 2 mass_ = 127.6_real32
1973 2 charge_ = 52.0_real32
1974 2 radius_ = 1.35_real32
1975 case('I')
1976 2 mass_ = 126.9_real32
1977 2 charge_ = 53.0_real32
1978 2 radius_ = 1.33_real32
1979 case('Xe')
1980 2 mass_ = 131.29_real32
1981 2 charge_ = 54.0_real32
1982 2 radius_ = 1.31_real32
1983 case('Cs')
1984 2 mass_ = 132.91_real32
1985 2 charge_ = 55.0_real32
1986 2 radius_ = 2.35_real32
1987 case('Ba')
1988 3 mass_ = 137.33_real32
1989 3 charge_ = 56.0_real32
1990 3 radius_ = 1.98_real32
1991 case('La')
1992 2 mass_ = 138.91_real32
1993 2 charge_ = 57.0_real32
1994 2 radius_ = 1.69_real32
1995 case('Ce')
1996 2 mass_ = 140.12_real32
1997 2 charge_ = 58.0_real32
1998 2 radius_ = 1.65_real32
1999 case('Pr')
2000 2 mass_ = 140.91_real32
2001 2 charge_ = 59.0_real32
2002 2 radius_ = 1.65_real32
2003 case('Nd')
2004 2 mass_ = 144.24_real32
2005 2 charge_ = 60.0_real32
2006 2 radius_ = 1.64_real32
2007 case('Pm')
2008 2 mass_ = 145.0_real32
2009 2 charge_ = 61.0_real32
2010 2 radius_ = 1.63_real32
2011 case('Sm')
2012 2 mass_ = 150.36_real32
2013 2 charge_ = 62.0_real32
2014 2 radius_ = 1.62_real32
2015 case('Eu')
2016 2 mass_ = 152.0_real32
2017 2 charge_ = 63.0_real32
2018 2 radius_ = 1.85_real32
2019 case('Gd')
2020 2 mass_ = 157.25_real32
2021 2 charge_ = 64.0_real32
2022 2 radius_ = 1.61_real32
2023 case('Tb')
2024 2 mass_ = 158.93_real32
2025 2 charge_ = 65.0_real32
2026 2 radius_ = 1.59_real32
2027 case('Dy')
2028 2 mass_ = 162.5_real32
2029 2 charge_ = 66.0_real32
2030 2 radius_ = 1.59_real32
2031 case('Ho')
2032 2 mass_ = 164.93_real32
2033 2 charge_ = 67.0_real32
2034 2 radius_ = 1.58_real32
2035 case('Er')
2036 2 mass_ = 167.26_real32
2037 2 charge_ = 68.0_real32
2038 2 radius_ = 1.57_real32
2039 case('Tm')
2040 2 mass_ = 168.93_real32
2041 2 charge_ = 69.0_real32
2042 2 radius_ = 1.56_real32
2043 case('Yb')
2044 2 mass_ = 173.05_real32
2045 2 charge_ = 70.0_real32
2046 2 radius_ = 1.74_real32
2047 case('Lu')
2048 2 mass_ = 174.97_real32
2049 2 charge_ = 71.0_real32
2050 2 radius_ = 1.56_real32
2051 case('Hf')
2052 2 mass_ = 178.49_real32
2053 2 charge_ = 72.0_real32
2054 2 radius_ = 1.44_real32
2055 case('Ta')
2056 2 mass_ = 180.95_real32
2057 2 charge_ = 73.0_real32
2058 2 radius_ = 1.34_real32
2059 case('W')
2060 2 mass_ = 183.84_real32
2061 2 charge_ = 74.0_real32
2062 2 radius_ = 1.3_real32
2063 case('Re')
2064 2 mass_ = 186.21_real32
2065 2 charge_ = 75.0_real32
2066 2 radius_ = 1.28_real32
2067 case('Os')
2068 2 mass_ = 190.23_real32
2069 2 charge_ = 76.0_real32
2070 2 radius_ = 1.26_real32
2071 case('Ir')
2072 2 mass_ = 192.22_real32
2073 2 charge_ = 77.0_real32
2074 2 radius_ = 1.27_real32
2075 case('Pt')
2076 2 mass_ = 195.08_real32
2077 2 charge_ = 78.0_real32
2078 2 radius_ = 1.3_real32
2079 case('Au')
2080 2 mass_ = 196.97_real32
2081 2 charge_ = 79.0_real32
2082 2 radius_ = 1.34_real32
2083 case('Hg')
2084 2 mass_ = 200.59_real32
2085 2 charge_ = 80.0_real32
2086 2 radius_ = 1.49_real32
2087 case('Tl')
2088 2 mass_ = 204.38_real32
2089 2 charge_ = 81.0_real32
2090 2 radius_ = 1.48_real32
2091 case('Pb')
2092 2 mass_ = 207.2_real32
2093 2 charge_ = 82.0_real32
2094 2 radius_ = 1.47_real32
2095 case('Bi')
2096 2 mass_ = 208.98_real32
2097 2 charge_ = 83.0_real32
2098 2 radius_ = 1.46_real32
2099 case('Po')
2100 2 mass_ = 209.0_real32
2101 2 charge_ = 84.0_real32
2102 2 radius_ = 1.45_real32
2103 case('At')
2104 2 mass_ = 210.0_real32
2105 2 charge_ = 85.0_real32
2106 2 radius_ = 1.44_real32
2107 case('Rn')
2108 2 mass_ = 222.0_real32
2109 2 charge_ = 86.0_real32
2110 2 radius_ = 1.43_real32
2111 case('Fr')
2112 2 mass_ = 223.0_real32
2113 2 charge_ = 87.0_real32
2114 2 radius_ = 2.6_real32
2115 case('Ra')
2116 2 mass_ = 226.0_real32
2117 2 charge_ = 88.0_real32
2118 2 radius_ = 2.21_real32
2119 case('Ac')
2120 2 mass_ = 227.0_real32
2121 2 charge_ = 89.0_real32
2122 2 radius_ = 1.86_real32
2123 case('Th')
2124 2 mass_ = 232.04_real32
2125 2 charge_ = 90.0_real32
2126 2 radius_ = 1.75_real32
2127 case('Pa')
2128 2 mass_ = 231.04_real32
2129 2 charge_ = 91.0_real32
2130 2 radius_ = 1.61_real32
2131 case('U')
2132 2 mass_ = 238.03_real32
2133 2 charge_ = 92.0_real32
2134 2 radius_ = 1.58_real32
2135 case('Np')
2136 2 mass_ = 237.0_real32
2137 2 charge_ = 93.0_real32
2138 2 radius_ = 1.55_real32
2139 case('Pu')
2140 2 mass_ = 244.0_real32
2141 2 charge_ = 94.0_real32
2142 2 radius_ = 1.53_real32
2143 case('Am')
2144 2 mass_ = 243.0_real32
2145 2 charge_ = 95.0_real32
2146 2 radius_ = 1.51_real32
2147 case('Cm')
2148 2 mass_ = 247.0_real32
2149 2 charge_ = 96.0_real32
2150 2 radius_ = 1.69_real32
2151 case('Bk')
2152 2 mass_ = 247.0_real32
2153 2 charge_ = 97.0_real32
2154 2 radius_ = 1.48_real32
2155 case('Cf')
2156 2 mass_ = 251.0_real32
2157 2 charge_ = 98.0_real32
2158 2 radius_ = 1.47_real32
2159 case('Es')
2160 2 mass_ = 252.0_real32
2161 2 charge_ = 99.0_real32
2162 2 radius_ = 1.46_real32
2163 case('Fm')
2164 2 mass_ = 257.0_real32
2165 2 charge_ = 100.0_real32
2166 2 radius_ = 1.45_real32
2167 case('Md')
2168 2 mass_ = 258.0_real32
2169 2 charge_ = 101.0_real32
2170 2 radius_ = 1.44_real32
2171 case('No')
2172 2 mass_ = 259.0_real32
2173 2 charge_ = 102.0_real32
2174 2 radius_ = 1.43_real32
2175 case('Lr')
2176 2 mass_ = 262.0_real32
2177 2 charge_ = 103.0_real32
2178 2 radius_ = 1.62_real32
2179 case('Rf')
2180 2 mass_ = 267.0_real32
2181 2 charge_ = 104.0_real32
2182 2 radius_ = 1.57_real32
2183 case('Db')
2184 2 mass_ = 270.0_real32
2185 2 charge_ = 105.0_real32
2186 2 radius_ = 1.49_real32
2187 case('Sg')
2188 2 mass_ = 271.0_real32
2189 2 charge_ = 106.0_real32
2190 2 radius_ = 1.43_real32
2191 case('Bh')
2192 2 mass_ = 270.0_real32
2193 2 charge_ = 107.0_real32
2194 2 radius_ = 1.41_real32
2195 case('Hs')
2196 2 mass_ = 277.0_real32
2197 2 charge_ = 108.0_real32
2198 2 radius_ = 1.34_real32
2199 case('Mt')
2200 2 mass_ = 276.0_real32
2201 2 charge_ = 109.0_real32
2202 2 radius_ = 1.29_real32
2203 case('Ds')
2204 2 mass_ = 281.0_real32
2205 2 charge_ = 110.0_real32
2206 2 radius_ = 1.28_real32
2207 case('Rg')
2208 2 mass_ = 280.0_real32
2209 2 charge_ = 111.0_real32
2210 2 radius_ = 1.21_real32
2211 case('Cn')
2212 2 mass_ = 285.0_real32
2213 2 charge_ = 112.0_real32
2214 2 radius_ = 1.22_real32
2215 case('Nh')
2216 2 mass_ = 284.0_real32
2217 2 charge_ = 113.0_real32
2218 2 radius_ = 1.21_real32
2219 case('Fl')
2220 2 mass_ = 289.0_real32
2221 2 charge_ = 114.0_real32
2222 2 radius_ = 1.21_real32
2223 case('Mc')
2224 2 mass_ = 288.0_real32
2225 2 charge_ = 115.0_real32
2226 2 radius_ = 1.21_real32
2227 case('Lv')
2228 2 mass_ = 293.0_real32
2229 2 charge_ = 116.0_real32
2230 2 radius_ = 1.21_real32
2231 case('Ts')
2232 2 mass_ = 294.0_real32
2233 2 charge_ = 117.0_real32
2234 2 radius_ = 1.21_real32
2235 case('Og')
2236 2 mass_ = 294.0_real32
2237 2 charge_ = 118.0_real32
2238 2 radius_ = 1.21_real32
2239 case default
2240 ! handle unknown element
2241 5 mass_ = 0.0_real32
2242 5 charge_ = 0.0_real32
2243
119/119
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 7 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 5 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 3 times.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 6 times.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 2 times.
✓ Branch 20 taken 2 times.
✓ Branch 21 taken 3 times.
✓ Branch 22 taken 2 times.
✓ Branch 23 taken 2 times.
✓ Branch 24 taken 2 times.
✓ Branch 25 taken 2 times.
✓ Branch 26 taken 2 times.
✓ Branch 27 taken 2 times.
✓ Branch 28 taken 2 times.
✓ Branch 29 taken 2 times.
✓ Branch 30 taken 2 times.
✓ Branch 31 taken 2 times.
✓ Branch 32 taken 2 times.
✓ Branch 33 taken 2 times.
✓ Branch 34 taken 2 times.
✓ Branch 35 taken 2 times.
✓ Branch 36 taken 2 times.
✓ Branch 37 taken 2 times.
✓ Branch 38 taken 2 times.
✓ Branch 39 taken 2 times.
✓ Branch 40 taken 2 times.
✓ Branch 41 taken 2 times.
✓ Branch 42 taken 2 times.
✓ Branch 43 taken 2 times.
✓ Branch 44 taken 2 times.
✓ Branch 45 taken 2 times.
✓ Branch 46 taken 2 times.
✓ Branch 47 taken 2 times.
✓ Branch 48 taken 2 times.
✓ Branch 49 taken 2 times.
✓ Branch 50 taken 2 times.
✓ Branch 51 taken 2 times.
✓ Branch 52 taken 2 times.
✓ Branch 53 taken 2 times.
✓ Branch 54 taken 2 times.
✓ Branch 55 taken 3 times.
✓ Branch 56 taken 2 times.
✓ Branch 57 taken 2 times.
✓ Branch 58 taken 2 times.
✓ Branch 59 taken 2 times.
✓ Branch 60 taken 2 times.
✓ Branch 61 taken 2 times.
✓ Branch 62 taken 2 times.
✓ Branch 63 taken 2 times.
✓ Branch 64 taken 2 times.
✓ Branch 65 taken 2 times.
✓ Branch 66 taken 2 times.
✓ Branch 67 taken 2 times.
✓ Branch 68 taken 2 times.
✓ Branch 69 taken 2 times.
✓ Branch 70 taken 2 times.
✓ Branch 71 taken 2 times.
✓ Branch 72 taken 2 times.
✓ Branch 73 taken 2 times.
✓ Branch 74 taken 2 times.
✓ Branch 75 taken 2 times.
✓ Branch 76 taken 2 times.
✓ Branch 77 taken 2 times.
✓ Branch 78 taken 2 times.
✓ Branch 79 taken 2 times.
✓ Branch 80 taken 2 times.
✓ Branch 81 taken 2 times.
✓ Branch 82 taken 2 times.
✓ Branch 83 taken 2 times.
✓ Branch 84 taken 2 times.
✓ Branch 85 taken 2 times.
✓ Branch 86 taken 2 times.
✓ Branch 87 taken 2 times.
✓ Branch 88 taken 2 times.
✓ Branch 89 taken 2 times.
✓ Branch 90 taken 2 times.
✓ Branch 91 taken 2 times.
✓ Branch 92 taken 2 times.
✓ Branch 93 taken 2 times.
✓ Branch 94 taken 2 times.
✓ Branch 95 taken 2 times.
✓ Branch 96 taken 2 times.
✓ Branch 97 taken 2 times.
✓ Branch 98 taken 2 times.
✓ Branch 99 taken 2 times.
✓ Branch 100 taken 2 times.
✓ Branch 101 taken 2 times.
✓ Branch 102 taken 2 times.
✓ Branch 103 taken 2 times.
✓ Branch 104 taken 2 times.
✓ Branch 105 taken 2 times.
✓ Branch 106 taken 2 times.
✓ Branch 107 taken 2 times.
✓ Branch 108 taken 2 times.
✓ Branch 109 taken 2 times.
✓ Branch 110 taken 2 times.
✓ Branch 111 taken 2 times.
✓ Branch 112 taken 2 times.
✓ Branch 113 taken 2 times.
✓ Branch 114 taken 2 times.
✓ Branch 115 taken 2 times.
✓ Branch 116 taken 2 times.
✓ Branch 117 taken 2 times.
✓ Branch 118 taken 5 times.
256 radius_ = 0.0_real32
2244 end select
2245
2246 !---------------------------------------------------------------------------
2247 ! Return the values
2248 !---------------------------------------------------------------------------
2249
2/2
✓ Branch 0 taken 246 times.
✓ Branch 1 taken 10 times.
256 if(present(mass)) mass = mass_
2250
2/2
✓ Branch 0 taken 246 times.
✓ Branch 1 taken 10 times.
256 if(present(charge)) charge = charge_
2251
1/2
✓ Branch 0 taken 256 times.
✗ Branch 1 not taken.
256 if(present(radius)) radius = radius_
2252
2253 256 end subroutine get_element_properties
2254 !###############################################################################
2255
2256
7/24
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
4 end module raffle__geom_rw
2257