| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module raffle__geom_extd | ||
| 2 | !! Module to extend the basis set to include images of atoms. | ||
| 3 | !! | ||
| 4 | !! This module is designed to extend the basis set to include images of atoms | ||
| 5 | !! within a specified distance of the unit cell. This is useful for | ||
| 6 | !! calculating interactions between atoms that are not within the unit cell. | ||
| 7 | use raffle__constants, only: real32, pi | ||
| 8 | use raffle__misc_linalg, only: modu, cross, inverse_3x3 | ||
| 9 | use raffle__geom_rw, only: basis_type, species_type | ||
| 10 | implicit none | ||
| 11 | |||
| 12 | |||
| 13 | private | ||
| 14 | |||
| 15 | public :: extended_basis_type | ||
| 16 | |||
| 17 | |||
| 18 | type, extends(basis_type) :: extended_basis_type | ||
| 19 | !! Extended basis set type | ||
| 20 | real(real32) :: max_extension | ||
| 21 | !! Maximum distance to extend the basis set | ||
| 22 | integer :: num_images | ||
| 23 | !! Number of images in the extended basis set | ||
| 24 | type(species_type), dimension(:), allocatable :: image_spec | ||
| 25 | !! Species type for the images | ||
| 26 | contains | ||
| 27 | procedure, pass(this) :: create_images | ||
| 28 | !! Create the images for the basis set | ||
| 29 | procedure, pass(this) :: update_images | ||
| 30 | !! Update the images for a specific atom | ||
| 31 | end type extended_basis_type | ||
| 32 | |||
| 33 | contains | ||
| 34 | |||
| 35 | !############################################################################### | ||
| 36 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 29 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
|
39 | subroutine create_images(this, max_bondlength, atom_ignore_list) |
| 37 | !! Create the images for the basis set. | ||
| 38 | implicit none | ||
| 39 | |||
| 40 | ! Arguments | ||
| 41 | class(extended_basis_type), intent(inout) :: this | ||
| 42 | !! Parent of the procedure. Instance of the extended basis. | ||
| 43 | real(real32), intent(in) :: max_bondlength | ||
| 44 | !! Maximum distance to extend the basis set. | ||
| 45 | integer, dimension(:,:), intent(in), optional :: atom_ignore_list | ||
| 46 | !! List of atoms to ignore when creating images. | ||
| 47 | |||
| 48 | ! Local variables | ||
| 49 | integer :: is, ia, i, j, k | ||
| 50 | !! Loop indices. | ||
| 51 | integer :: amax, bmax, cmax | ||
| 52 | !! Maximum number of lattice vectors to consider. | ||
| 53 | real(real32), dimension(3) :: vtmp1 | ||
| 54 | !! Temporary vector for storing atom positions. | ||
| 55 |
9/12✓ Branch 0 taken 48 times.
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 39 times.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 48 times.
✓ Branch 8 taken 48 times.
✓ Branch 9 taken 39 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 48 times.
|
270 | type(species_type), dimension(this%nspec) :: image_species |
| 56 | !! Temporary store for the images. | ||
| 57 | 39 | integer, dimension(:,:), allocatable :: atom_ignore_list_ | |
| 58 | !! List of atoms to ignore when creating images. | ||
| 59 | |||
| 60 | |||
| 61 | !--------------------------------------------------------------------------- | ||
| 62 | ! initialise the atom_ignore_list | ||
| 63 | !--------------------------------------------------------------------------- | ||
| 64 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 29 times.
|
39 | if(present(atom_ignore_list)) then |
| 65 |
7/14✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 50 times.
✓ Branch 11 taken 10 times.
✓ Branch 12 taken 100 times.
✓ Branch 13 taken 50 times.
|
160 | atom_ignore_list_ = atom_ignore_list(:,:) |
| 66 | else | ||
| 67 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 29 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 29 times.
|
29 | allocate(atom_ignore_list_(0,0)) |
| 68 | end if | ||
| 69 | |||
| 70 | |||
| 71 | !--------------------------------------------------------------------------- | ||
| 72 | ! get the maximum number of lattice vectors to consider | ||
| 73 | ! NOTE: this is not perfect | ||
| 74 | ! won't work for extremely acute/obtuse angle cells | ||
| 75 | ! (due to diagonal path being shorter than individual lattice vectors) | ||
| 76 | !--------------------------------------------------------------------------- | ||
| 77 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 37 times.
|
39 | amax = ceiling(max_bondlength/modu(this%lat(1,:))) |
| 78 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 37 times.
|
39 | bmax = ceiling(max_bondlength/modu(this%lat(2,:))) |
| 79 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 37 times.
|
39 | cmax = ceiling(max_bondlength/modu(this%lat(3,:))) |
| 80 | |||
| 81 | |||
| 82 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 39 times.
|
87 | spec_loop: do is = 1, this%nspec |
| 83 | allocate( & | ||
| 84 | ✗ | image_species(is)%atom( & | |
| 85 | this%spec(is)%num*(2*amax+2)*(2*bmax+2)*(2*cmax+2), & | ||
| 86 | size(this%spec(is)%atom,2) & | ||
| 87 | ) & | ||
| 88 |
10/20✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 48 times.
✓ Branch 6 taken 48 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 48 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 48 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 48 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 48 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 48 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 48 times.
|
48 | ) |
| 89 | 48 | image_species(is)%num = 0 | |
| 90 | 48 | image_species(is)%mass = this%spec(is)%mass | |
| 91 | 48 | image_species(is)%charge = this%spec(is)%charge | |
| 92 | 48 | image_species(is)%radius = this%spec(is)%radius | |
| 93 | 48 | image_species(is)%name = this%spec(is)%name | |
| 94 |
2/2✓ Branch 0 taken 307 times.
✓ Branch 1 taken 48 times.
|
394 | atom_loop: do ia = 1, this%spec(is)%num |
| 95 |
2/2✓ Branch 0 taken 604 times.
✓ Branch 1 taken 257 times.
|
861 | do i = 1, size(atom_ignore_list_,2), 1 |
| 96 |
6/6✓ Branch 0 taken 1187 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 554 times.
✓ Branch 3 taken 633 times.
✓ Branch 4 taken 50 times.
✓ Branch 5 taken 554 times.
|
1494 | if(all(atom_ignore_list_(:,i).eq.[is,ia])) cycle atom_loop |
| 97 | end do | ||
| 98 |
2/2✓ Branch 0 taken 1514 times.
✓ Branch 1 taken 257 times.
|
1819 | do i=-amax,amax+1,1 |
| 99 | 1514 | vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32) | |
| 100 |
2/2✓ Branch 0 taken 9220 times.
✓ Branch 1 taken 1514 times.
|
10991 | do j=-bmax,bmax+1,1 |
| 101 | 9220 | vtmp1(2) = this%spec(is)%atom(ia,2) + real(j, real32) | |
| 102 |
2/2✓ Branch 0 taken 49216 times.
✓ Branch 1 taken 9220 times.
|
59950 | do k=-cmax,cmax+1,1 |
| 103 |
2/2✓ Branch 0 taken 257 times.
✓ Branch 1 taken 48959 times.
|
49216 | if( i .eq. 0 .and. j .eq. 0 .and. k .eq. 0 ) cycle |
| 104 | 48959 | vtmp1(3) = this%spec(is)%atom(ia,3) + real(k, real32) | |
| 105 |
2/2✓ Branch 1 taken 11447 times.
✓ Branch 2 taken 37512 times.
|
48959 | if( get_distance_from_unit_cell(vtmp1, this%lat) .le. & |
| 106 | 9220 | max_bondlength ) then | |
| 107 | ! add the image to the list | ||
| 108 | 11447 | image_species(is)%num = image_species(is)%num + 1 | |
| 109 |
2/2✓ Branch 0 taken 34341 times.
✓ Branch 1 taken 11447 times.
|
45788 | image_species(is)%atom(image_species(is)%num,:3) = vtmp1 |
| 110 | end if | ||
| 111 | end do | ||
| 112 | end do | ||
| 113 | end do | ||
| 114 | end do atom_loop | ||
| 115 | end do spec_loop | ||
| 116 | |||
| 117 | |||
| 118 |
11/20✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 39 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 39 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 39 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 39 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 39 times.
✓ Branch 17 taken 48 times.
✓ Branch 18 taken 39 times.
✓ Branch 19 taken 48 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 48 times.
|
87 | allocate(this%image_spec(this%nspec)) |
| 119 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 39 times.
|
87 | do is = 1, this%nspec |
| 120 | 48 | this%image_spec(is)%num = image_species(is)%num | |
| 121 | 48 | this%image_spec(is)%mass = image_species(is)%mass | |
| 122 | 48 | this%image_spec(is)%charge = image_species(is)%charge | |
| 123 | 48 | this%image_spec(is)%radius = image_species(is)%radius | |
| 124 | 48 | this%image_spec(is)%name = image_species(is)%name | |
| 125 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
|
48 | if(image_species(is)%num .eq. 0) cycle |
| 126 | ✗ | allocate(this%image_spec(is)%atom( & | |
| 127 | image_species(is)%num, & | ||
| 128 | size(image_species(is)%atom,2) & | ||
| 129 |
9/18✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 48 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 48 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 48 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 48 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 48 times.
|
48 | ) ) |
| 130 | this%image_spec(is)%atom(:,:) = & | ||
| 131 |
4/4✓ Branch 0 taken 144 times.
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 34341 times.
✓ Branch 3 taken 144 times.
|
34572 | image_species(is)%atom(:image_species(is)%num,:) |
| 132 | end do | ||
| 133 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 39 times.
|
87 | this%num_images = sum( this%image_spec(:)%num ) |
| 134 | |||
| 135 |
1/2✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
|
135 | end subroutine create_images |
| 136 | !############################################################################### | ||
| 137 | |||
| 138 | |||
| 139 | !############################################################################### | ||
| 140 | 38 | subroutine update_images(this, max_bondlength, is, ia) | |
| 141 | !! Update the images for a specific atom | ||
| 142 | implicit none | ||
| 143 | |||
| 144 | ! Arguments | ||
| 145 | class(extended_basis_type), intent(inout) :: this | ||
| 146 | !! Parent of the procedure. Instance of the extended basis. | ||
| 147 | real(real32), intent(in) :: max_bondlength | ||
| 148 | !! Maximum distance to extend the basis set. | ||
| 149 | integer, intent(in) :: is, ia | ||
| 150 | !! Species and atom index to update. | ||
| 151 | |||
| 152 | |||
| 153 | ! Local variables | ||
| 154 | integer :: i, j, k, num_images, dim | ||
| 155 | !! Loop indices. | ||
| 156 | integer :: amax, bmax, cmax | ||
| 157 | !! Maximum number of lattice vectors to consider. | ||
| 158 | 38 | type(species_type) :: image_species | |
| 159 | !! Temporary store for the images. | ||
| 160 | real(real32), dimension(3) :: vtmp1 | ||
| 161 | !! Temporary vector for storing atom positions. | ||
| 162 | |||
| 163 | |||
| 164 | !--------------------------------------------------------------------------- | ||
| 165 | ! get the maximum number of lattice vectors to consider | ||
| 166 | ! NOTE: this is not perfect | ||
| 167 | ! won't work for extremely acute/obtuse angle cells | ||
| 168 | ! (due to diagonal path being shorter than individual lattice vectors) | ||
| 169 | !--------------------------------------------------------------------------- | ||
| 170 | 38 | num_images = this%image_spec(is)%num | |
| 171 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 37 times.
|
38 | amax = ceiling(max_bondlength/modu(this%lat(1,:))) |
| 172 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 37 times.
|
38 | bmax = ceiling(max_bondlength/modu(this%lat(2,:))) |
| 173 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 37 times.
|
38 | cmax = ceiling(max_bondlength/modu(this%lat(3,:))) |
| 174 | 38 | dim = 3 | |
| 175 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 38 times.
|
86 | do i = 1, this%nspec |
| 176 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
|
86 | if ( size(this%spec(i)%atom,2) .gt. dim) dim = size(this%spec(i)%atom,2) |
| 177 | end do | ||
| 178 | allocate( & | ||
| 179 | ✗ | image_species%atom( & | |
| 180 | num_images + & | ||
| 181 | (2*amax+2)*(2*bmax+2)*(2*cmax+2), & | ||
| 182 | dim & | ||
| 183 | ) & | ||
| 184 |
9/18✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 38 times.
✓ Branch 4 taken 38 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 38 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 38 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 38 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 38 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 38 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 38 times.
|
38 | ) |
| 185 |
1/2✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
|
38 | if( num_images .ne. 0 ) then |
| 186 | image_species%atom(:num_images,:3) = & | ||
| 187 |
4/4✓ Branch 0 taken 114 times.
✓ Branch 1 taken 38 times.
✓ Branch 2 taken 40212 times.
✓ Branch 3 taken 114 times.
|
40364 | this%image_spec(is)%atom(:num_images,:3) |
| 188 | end if | ||
| 189 | |||
| 190 | |||
| 191 |
2/2✓ Branch 0 taken 224 times.
✓ Branch 1 taken 38 times.
|
262 | do i=-amax,amax+1,1 |
| 192 | 224 | vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32) | |
| 193 |
2/2✓ Branch 0 taken 1336 times.
✓ Branch 1 taken 224 times.
|
1598 | do j=-bmax,bmax+1,1 |
| 194 | 1336 | vtmp1(2) = this%spec(is)%atom(ia,2) + real(j, real32) | |
| 195 |
2/2✓ Branch 0 taken 5336 times.
✓ Branch 1 taken 1336 times.
|
6896 | do k=-cmax,cmax+1,1 |
| 196 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 5298 times.
|
5336 | if( i .eq. 0 .and. j .eq. 0 .and. k .eq. 0 ) cycle |
| 197 | 5298 | vtmp1(3) = this%spec(is)%atom(ia,3) + real(k, real32) | |
| 198 |
2/2✓ Branch 1 taken 1280 times.
✓ Branch 2 taken 4018 times.
|
5298 | if( get_distance_from_unit_cell(vtmp1, this%lat) .le. & |
| 199 | 1336 | max_bondlength ) then | |
| 200 | ! add the image to the list | ||
| 201 | 1280 | num_images = num_images + 1 | |
| 202 |
2/2✓ Branch 0 taken 3840 times.
✓ Branch 1 taken 1280 times.
|
5120 | image_species%atom(num_images,:3) = vtmp1 |
| 203 | end if | ||
| 204 | end do | ||
| 205 | end do | ||
| 206 | end do | ||
| 207 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 38 times.
|
38 | if( num_images .eq. this%image_spec(is)%num ) return |
| 208 | |||
| 209 | |||
| 210 | 38 | this%image_spec(is)%num = num_images | |
| 211 |
2/4✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 38 times.
|
38 | if(allocated(this%image_spec(is)%atom)) deallocate(this%image_spec(is)%atom) |
| 212 | ✗ | allocate(this%image_spec(is)%atom( & | |
| 213 | num_images, & | ||
| 214 | size(image_species%atom,2) & | ||
| 215 |
9/18✓ Branch 0 taken 38 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 38 times.
✓ Branch 4 taken 38 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 38 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 38 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 38 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 38 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 38 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 38 times.
|
38 | ) ) |
| 216 | this%image_spec(is)%atom(:,:) = & | ||
| 217 |
4/4✓ Branch 0 taken 114 times.
✓ Branch 1 taken 38 times.
✓ Branch 2 taken 44052 times.
✓ Branch 3 taken 114 times.
|
44204 | image_species%atom(:num_images,:) |
| 218 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 38 times.
|
38 | deallocate(image_species%atom) |
| 219 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 38 times.
|
86 | this%num_images = sum( this%image_spec(:)%num ) |
| 220 | |||
| 221 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 38 times.
|
38 | end subroutine update_images |
| 222 | !############################################################################### | ||
| 223 | |||
| 224 | |||
| 225 | !############################################################################### | ||
| 226 | 54257 | function get_distance_from_unit_cell( & | |
| 227 | point, lattice, closest_point, is_cartesian& | ||
| 228 | ) result(distance) | ||
| 229 | !! Get the distance of a point from the unit cell. | ||
| 230 | implicit none | ||
| 231 | |||
| 232 | ! Arguments | ||
| 233 | real(real32), intent(in) :: point(3) | ||
| 234 | !! Query point. | ||
| 235 | real(real32), intent(in) :: lattice(3,3) | ||
| 236 | !! Lattice vectors. | ||
| 237 | real(real32), intent(out), optional :: closest_point(3) | ||
| 238 | !! Closest point on the unit cell surface. | ||
| 239 | logical, optional, intent(in) :: is_cartesian | ||
| 240 | !! Boolean whether the point is in cartesian coordinates. | ||
| 241 | real(real32) :: distance | ||
| 242 | !! Distance of the point from the unit cell. | ||
| 243 | |||
| 244 | ! Local variables | ||
| 245 | integer :: i, j, k | ||
| 246 | !! Loop indices. | ||
| 247 | real(real32), dimension(3) :: point_ | ||
| 248 | !! Point in cartesian coordinates. | ||
| 249 | real(real32), dimension(3,3) :: inverse_lattice | ||
| 250 | !! Inverse of the lattice vectors. | ||
| 251 | real(real32), dimension(3) :: normal | ||
| 252 | !! Normal vector to the plane. | ||
| 253 | real(real32), dimension(3) :: plane_point | ||
| 254 | !! Point on the plane. | ||
| 255 | real(real32), dimension(3) :: projection, closest_point_ | ||
| 256 | !! Projection of the point onto the plane. | ||
| 257 | real(real32), dimension(3) :: inverse_projection | ||
| 258 | !! Inverse projection of the point onto the plane. | ||
| 259 | real(real32) :: min_distance | ||
| 260 | !! Minimum distance to the unit cell. | ||
| 261 | logical :: is_outside | ||
| 262 | !! Boolean whether the point is outside the unit cell. | ||
| 263 | integer, dimension(3) :: index_list | ||
| 264 | !! List of indices for the lattice vectors. | ||
| 265 | logical :: is_cartesian_ | ||
| 266 | !! Boolean whether the point is in cartesian coordinates. | ||
| 267 | |||
| 268 | |||
| 269 | !--------------------------------------------------------------------------- | ||
| 270 | ! check if the point is in cartesian coordinates | ||
| 271 | !--------------------------------------------------------------------------- | ||
| 272 | 54257 | is_outside = .false. | |
| 273 | 54257 | is_cartesian_ = .false. | |
| 274 | 54257 | index_list = [1, 2, 3] | |
| 275 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54257 times.
|
54257 | if(present(is_cartesian)) is_cartesian_ = is_cartesian |
| 276 | 54257 | inverse_lattice = inverse_3x3( lattice ) | |
| 277 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54257 times.
|
54257 | if(is_cartesian_) then |
| 278 | ! Convert point to fractional coordinates | ||
| 279 | ! point_ = matmul(LUinv(lattice), point) | ||
| 280 | ✗ | point_ = point | |
| 281 | else | ||
| 282 | 54257 | point_ = matmul(point, lattice) | |
| 283 | end if | ||
| 284 | |||
| 285 | 54257 | min_distance = huge(1._real32) | |
| 286 | |||
| 287 | ! get projection of point onto each face of the lattice | ||
| 288 | ! get the length of the projection vector | ||
| 289 | ! if negative, then the point is inside the unit cell | ||
| 290 | ! if positive, then the point is outside the unit cell | ||
| 291 | ! if the projection falls outside of the cell edges, use edge or corner | ||
| 292 | ! distances | ||
| 293 |
2/2✓ Branch 0 taken 162771 times.
✓ Branch 1 taken 54257 times.
|
217028 | face_loop: do i = 1, 3 |
| 294 |
2/2✓ Branch 1 taken 488313 times.
✓ Branch 2 taken 162771 times.
|
651084 | index_list = cshift(index_list, 1) |
| 295 | 162771 | plane_point = 0._real32 | |
| 296 |
2/2✓ Branch 0 taken 325542 times.
✓ Branch 1 taken 162771 times.
|
542570 | direction_loop: do j = 1, 2 |
| 297 | normal = (-1._real32)**j * cross( & | ||
| 298 | [ lattice(index_list(2),:3) ], & | ||
| 299 | [ lattice(index_list(3),:3) ] & | ||
| 300 |
10/10✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
✓ Branch 2 taken 976626 times.
✓ Branch 3 taken 325542 times.
✓ Branch 4 taken 976626 times.
✓ Branch 5 taken 325542 times.
✓ Branch 6 taken 976626 times.
✓ Branch 7 taken 325542 times.
✓ Branch 9 taken 976626 times.
✓ Branch 10 taken 325542 times.
|
5208672 | ) |
| 301 |
8/8✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
✓ Branch 2 taken 341862 times.
✓ Branch 3 taken 634764 times.
✓ Branch 4 taken 325542 times.
✓ Branch 5 taken 16320 times.
✓ Branch 6 taken 976626 times.
✓ Branch 7 taken 325542 times.
|
2278794 | normal = normal / norm2(normal) |
| 302 | 325542 | projection = project_point_onto_plane(point_, plane_point, normal) | |
| 303 | |||
| 304 | ! check if point minus projection is negative | ||
| 305 | ! if so, it is on the wrong side of the plane and should be ignored | ||
| 306 |
4/4✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
✓ Branch 2 taken 182608 times.
✓ Branch 3 taken 142934 times.
|
1302168 | if( dot_product(point_ - projection, normal) .lt. 0._real32 )then |
| 307 |
2/2✓ Branch 0 taken 547824 times.
✓ Branch 1 taken 182608 times.
|
730432 | plane_point = plane_point + lattice(index_list(1),:) |
| 308 | 182608 | cycle direction_loop | |
| 309 | end if | ||
| 310 | 142934 | is_outside = .true. | |
| 311 | |||
| 312 | ! check if projection is outside the surface | ||
| 313 | |||
| 314 | 142934 | inverse_projection = matmul(projection, inverse_lattice) | |
| 315 | if( & | ||
| 316 |
10/10✓ Branch 0 taken 336566 times.
✓ Branch 1 taken 63914 times.
✓ Branch 2 taken 79020 times.
✓ Branch 3 taken 257546 times.
✓ Branch 4 taken 311500 times.
✓ Branch 5 taken 43832 times.
✓ Branch 6 taken 99102 times.
✓ Branch 7 taken 212398 times.
✓ Branch 8 taken 135654 times.
✓ Branch 9 taken 7280 times.
|
612878 | any( inverse_projection .lt. 0._real32 ) .or. & |
| 317 | any( inverse_projection .gt. 1._real32 ) & | ||
| 318 | ) then | ||
| 319 | ! projection is outside the surface | ||
| 320 | ! check if the projection is outside the edges | ||
| 321 | ! if it is, then the closest point is the edge or corner | ||
| 322 | ! if it is not, then the closest point is the projection | ||
| 323 |
2/2✓ Branch 0 taken 406962 times.
✓ Branch 1 taken 135654 times.
|
542616 | do k = 1, 3 |
| 324 |
2/2✓ Branch 0 taken 94960 times.
✓ Branch 1 taken 312002 times.
|
542616 | if( inverse_projection(k) .lt. 0._real32 ) then |
| 325 | 94960 | inverse_projection(k) = 0._real32 | |
| 326 |
2/2✓ Branch 0 taken 128108 times.
✓ Branch 1 taken 183894 times.
|
312002 | else if( inverse_projection(k) .gt. 1._real32 ) then |
| 327 | 128108 | inverse_projection(k) = 1._real32 | |
| 328 | end if | ||
| 329 | end do | ||
| 330 | end if | ||
| 331 | 142934 | projection = matmul(inverse_projection, lattice) | |
| 332 |
6/6✓ Branch 0 taken 428802 times.
✓ Branch 1 taken 142934 times.
✓ Branch 2 taken 369398 times.
✓ Branch 3 taken 59404 times.
✓ Branch 4 taken 215103 times.
✓ Branch 5 taken 154295 times.
|
571736 | distance = norm2(point_ - projection) |
| 333 |
2/2✓ Branch 0 taken 56077 times.
✓ Branch 1 taken 86857 times.
|
142934 | if( distance .lt. min_distance ) then |
| 334 | 56077 | min_distance = distance | |
| 335 | 56077 | closest_point_ = projection | |
| 336 | end if | ||
| 337 | |||
| 338 | !! makes it apply to the next iteration | ||
| 339 |
2/2✓ Branch 0 taken 428802 times.
✓ Branch 1 taken 142934 times.
|
734507 | plane_point = plane_point + lattice(index_list(1),:) |
| 340 | end do direction_loop | ||
| 341 | end do face_loop | ||
| 342 | |||
| 343 |
2/2✓ Branch 0 taken 54249 times.
✓ Branch 1 taken 8 times.
|
54257 | if( is_outside ) then |
| 344 | 54249 | distance = min_distance | |
| 345 | else | ||
| 346 | 8 | distance = 0._real32 | |
| 347 | end if | ||
| 348 | |||
| 349 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54257 times.
|
54257 | if( present(closest_point) ) then |
| 350 | ✗ | if(is_cartesian_) then | |
| 351 | ✗ | closest_point = closest_point_ | |
| 352 | else | ||
| 353 | ✗ | closest_point = matmul(closest_point_, inverse_lattice) | |
| 354 | end if | ||
| 355 | end if | ||
| 356 | |||
| 357 | 54257 | end function get_distance_from_unit_cell | |
| 358 | !############################################################################### | ||
| 359 | |||
| 360 | |||
| 361 | !############################################################################### | ||
| 362 |
1/2✓ Branch 0 taken 325542 times.
✗ Branch 1 not taken.
|
325542 | function project_point_onto_plane(point, plane_point, normal) result(output) |
| 363 | !! Project a point onto a plane. | ||
| 364 | implicit none | ||
| 365 | |||
| 366 | ! Arguments | ||
| 367 | real(real32), dimension(3), intent(in) :: point | ||
| 368 | !! Point to project. | ||
| 369 | real(real32), dimension(3), intent(in) :: plane_point | ||
| 370 | !! Point on the plane. | ||
| 371 | real(real32), dimension(3), intent(in) :: normal | ||
| 372 | !! Normal vector to the plane. | ||
| 373 | real(real32), dimension(3) :: output | ||
| 374 | !! Projected point. | ||
| 375 | |||
| 376 | ! Local variables | ||
| 377 | real(real32) :: distance | ||
| 378 | !! Distance of the point from the plane. | ||
| 379 | real(real32), dimension(3) :: vector_to_plane | ||
| 380 | !! Vector from the point to the plane. | ||
| 381 | |||
| 382 |
2/2✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
|
1302168 | vector_to_plane = point - plane_point |
| 383 | |||
| 384 | distance = & | ||
| 385 |
4/4✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
✓ Branch 2 taken 976626 times.
✓ Branch 3 taken 325542 times.
|
2278794 | dot_product(vector_to_plane, normal) / dot_product(normal, normal) |
| 386 | |||
| 387 |
2/2✓ Branch 0 taken 976626 times.
✓ Branch 1 taken 325542 times.
|
1302168 | output = point - distance * normal |
| 388 | |||
| 389 | 325542 | end function project_point_onto_plane | |
| 390 | !############################################################################### | ||
| 391 | |||
| 392 | ✗ | end module raffle__geom_extd | |
| 393 |