GCC Code Coverage Report


Directory: src/fortran/lib/
File: mod_geom_rw.f90
Date: 2025-04-05 12:17:58
Exec Total Coverage
Lines: 591 906 65.2%
Functions: 0 0 -%
Branches: 556 1463 38.0%

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