| 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: 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 | 54 | subroutine create_images(this, max_bondlength) | |
| 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 | |||
| 46 | ! Local variables | ||
| 47 | integer :: is, ia, i, j, k | ||
| 48 | !! Loop indices. | ||
| 49 | integer :: amax, bmax, cmax | ||
| 50 | !! Maximum number of lattice vectors to consider. | ||
| 51 | real(real32), dimension(3) :: vtmp1 | ||
| 52 | !! Temporary vector for storing atom positions. | ||
| 53 |
13/20✓ Branch 0 taken 67 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 67 times.
✓ Branch 3 taken 54 times.
✓ Branch 4 taken 67 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 67 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 67 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 67 times.
✓ Branch 12 taken 67 times.
✓ Branch 13 taken 54 times.
✓ Branch 14 taken 67 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 67 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 67 times.
|
510 | type(species_type), dimension(this%nspec) :: image_species |
| 54 | !! Temporary store for the images. | ||
| 55 | |||
| 56 | |||
| 57 | !--------------------------------------------------------------------------- | ||
| 58 | ! get the maximum number of lattice vectors to consider | ||
| 59 | ! NOTE: this is not perfect | ||
| 60 | ! won't work for extremely acute/obtuse angle cells | ||
| 61 | ! (due to diagonal path being shorter than individual lattice vectors) | ||
| 62 | !--------------------------------------------------------------------------- | ||
| 63 |
7/8✓ Branch 0 taken 162 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 62 times.
✓ Branch 3 taken 100 times.
✓ Branch 4 taken 62 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 52 times.
|
216 | amax = ceiling(max_bondlength/norm2(this%lat(1,:))) |
| 64 |
7/8✓ Branch 0 taken 162 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 62 times.
✓ Branch 3 taken 100 times.
✓ Branch 4 taken 62 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 52 times.
|
216 | bmax = ceiling(max_bondlength/norm2(this%lat(2,:))) |
| 65 |
7/8✓ Branch 0 taken 162 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 54 times.
✓ Branch 3 taken 108 times.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 52 times.
|
216 | cmax = ceiling(max_bondlength/norm2(this%lat(3,:))) |
| 66 | |||
| 67 | |||
| 68 |
2/2✓ Branch 0 taken 67 times.
✓ Branch 1 taken 54 times.
|
121 | spec_loop: do is = 1, this%nspec |
| 69 | allocate( & | ||
| 70 | ✗ | image_species(is)%atom( & | |
| 71 | this%spec(is)%num*(2*amax+2)*(2*bmax+2)*(2*cmax+2), & | ||
| 72 | size(this%spec(is)%atom,2) & | ||
| 73 | ) & | ||
| 74 |
10/20✓ Branch 0 taken 67 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 67 times.
✓ Branch 6 taken 67 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 67 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 67 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 67 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 67 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 67 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 67 times.
|
67 | ) |
| 75 | 67 | image_species(is)%num = 0 | |
| 76 | 67 | image_species(is)%mass = this%spec(is)%mass | |
| 77 | 67 | image_species(is)%charge = this%spec(is)%charge | |
| 78 | 67 | image_species(is)%radius = this%spec(is)%radius | |
| 79 | 67 | image_species(is)%name = this%spec(is)%name | |
| 80 |
2/2✓ Branch 0 taken 411 times.
✓ Branch 1 taken 67 times.
|
532 | atom_loop: do ia = 1, this%spec(is)%num |
| 81 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 361 times.
|
411 | if(.not.this%spec(is)%atom_mask(ia)) cycle atom_loop |
| 82 |
2/2✓ Branch 0 taken 2170 times.
✓ Branch 1 taken 361 times.
|
2598 | do i=-amax,amax+1,1 |
| 83 | 2170 | vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32) | |
| 84 |
2/2✓ Branch 0 taken 13412 times.
✓ Branch 1 taken 2170 times.
|
15943 | do j=-bmax,bmax+1,1 |
| 85 | 13412 | vtmp1(2) = this%spec(is)%atom(ia,2) + real(j, real32) | |
| 86 |
2/2✓ Branch 0 taken 72320 times.
✓ Branch 1 taken 13412 times.
|
87902 | do k=-cmax,cmax+1,1 |
| 87 |
2/2✓ Branch 0 taken 361 times.
✓ Branch 1 taken 71959 times.
|
72320 | if( i .eq. 0 .and. j .eq. 0 .and. k .eq. 0 ) cycle |
| 88 | 71959 | vtmp1(3) = this%spec(is)%atom(ia,3) + real(k, real32) | |
| 89 |
2/2✓ Branch 1 taken 16791 times.
✓ Branch 2 taken 55168 times.
|
71959 | if( get_distance_from_unit_cell(vtmp1, this%lat) .le. & |
| 90 | 13412 | max_bondlength ) then | |
| 91 | ! add the image to the list | ||
| 92 | 16791 | image_species(is)%num = image_species(is)%num + 1 | |
| 93 |
2/2✓ Branch 0 taken 50373 times.
✓ Branch 1 taken 16791 times.
|
67164 | image_species(is)%atom(image_species(is)%num,:3) = vtmp1 |
| 94 | end if | ||
| 95 | end do | ||
| 96 | end do | ||
| 97 | end do | ||
| 98 | end do atom_loop | ||
| 99 | end do spec_loop | ||
| 100 | |||
| 101 | |||
| 102 |
13/24✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 54 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 54 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 54 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 54 times.
✓ Branch 17 taken 67 times.
✓ Branch 18 taken 54 times.
✓ Branch 19 taken 67 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 67 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 67 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 67 times.
|
121 | allocate(this%image_spec(this%nspec)) |
| 103 |
2/2✓ Branch 0 taken 67 times.
✓ Branch 1 taken 54 times.
|
121 | do is = 1, this%nspec |
| 104 | 67 | this%image_spec(is)%num = image_species(is)%num | |
| 105 | 67 | this%image_spec(is)%mass = image_species(is)%mass | |
| 106 | 67 | this%image_spec(is)%charge = image_species(is)%charge | |
| 107 | 67 | this%image_spec(is)%radius = image_species(is)%radius | |
| 108 | 67 | this%image_spec(is)%name = image_species(is)%name | |
| 109 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 67 times.
|
67 | if(image_species(is)%num .eq. 0) cycle |
| 110 | ✗ | allocate(this%image_spec(is)%atom( & | |
| 111 | image_species(is)%num, & | ||
| 112 | size(image_species(is)%atom,2) & | ||
| 113 |
9/18✓ Branch 0 taken 67 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 67 times.
✓ Branch 4 taken 67 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 67 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 67 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 67 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 67 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 67 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 67 times.
|
67 | ) ) |
| 114 | this%image_spec(is)%atom(:,:) = & | ||
| 115 |
4/4✓ Branch 0 taken 201 times.
✓ Branch 1 taken 67 times.
✓ Branch 2 taken 50373 times.
✓ Branch 3 taken 201 times.
|
50695 | image_species(is)%atom(:image_species(is)%num,:) |
| 116 | end do | ||
| 117 |
2/2✓ Branch 0 taken 67 times.
✓ Branch 1 taken 54 times.
|
121 | this%num_images = sum( this%image_spec(:)%num ) |
| 118 | |||
| 119 | 268 | end subroutine create_images | |
| 120 | !############################################################################### | ||
| 121 | |||
| 122 | |||
| 123 | !############################################################################### | ||
| 124 | 46 | subroutine update_images(this, max_bondlength, is, ia) | |
| 125 | !! Update the images for a specific atom | ||
| 126 | implicit none | ||
| 127 | |||
| 128 | ! Arguments | ||
| 129 | class(extended_basis_type), intent(inout) :: this | ||
| 130 | !! Parent of the procedure. Instance of the extended basis. | ||
| 131 | real(real32), intent(in) :: max_bondlength | ||
| 132 | !! Maximum distance to extend the basis set. | ||
| 133 | integer, intent(in) :: is, ia | ||
| 134 | !! Species and atom index to update. | ||
| 135 | |||
| 136 | |||
| 137 | ! Local variables | ||
| 138 | integer :: i, j, k, num_images, dim | ||
| 139 | !! Loop indices. | ||
| 140 | integer :: amax, bmax, cmax | ||
| 141 | !! Maximum number of lattice vectors to consider. | ||
| 142 | 46 | type(species_type) :: image_species | |
| 143 | !! Temporary store for the images. | ||
| 144 | real(real32), dimension(3) :: vtmp1 | ||
| 145 | !! Temporary vector for storing atom positions. | ||
| 146 | |||
| 147 | |||
| 148 | !--------------------------------------------------------------------------- | ||
| 149 | ! get the maximum number of lattice vectors to consider | ||
| 150 | ! NOTE: this is not perfect | ||
| 151 | ! won't work for extremely acute/obtuse angle cells | ||
| 152 | ! (due to diagonal path being shorter than individual lattice vectors) | ||
| 153 | !--------------------------------------------------------------------------- | ||
| 154 | 46 | num_images = this%image_spec(is)%num | |
| 155 |
7/8✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 46 times.
✓ Branch 3 taken 92 times.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 45 times.
|
184 | amax = ceiling(max_bondlength/norm2(this%lat(1,:))) |
| 156 |
7/8✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 46 times.
✓ Branch 3 taken 92 times.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 45 times.
|
184 | bmax = ceiling(max_bondlength/norm2(this%lat(2,:))) |
| 157 |
7/8✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 46 times.
✓ Branch 3 taken 92 times.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 45 times.
|
184 | cmax = ceiling(max_bondlength/norm2(this%lat(3,:))) |
| 158 | 46 | dim = 3 | |
| 159 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 46 times.
|
102 | do i = 1, this%nspec |
| 160 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
102 | if ( size(this%spec(i)%atom,2) .gt. dim) dim = size(this%spec(i)%atom,2) |
| 161 | end do | ||
| 162 | allocate( & | ||
| 163 | ✗ | image_species%atom( & | |
| 164 | num_images + & | ||
| 165 | (2*amax+2)*(2*bmax+2)*(2*cmax+2), & | ||
| 166 | dim & | ||
| 167 | ) & | ||
| 168 |
9/18✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 46 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 46 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 46 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 46 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 46 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 46 times.
|
46 | ) |
| 169 |
1/2✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
|
46 | if( num_images .ne. 0 ) then |
| 170 | image_species%atom(:num_images,:3) = & | ||
| 171 |
4/4✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 54015 times.
✓ Branch 3 taken 138 times.
|
54199 | this%image_spec(is)%atom(:num_images,:3) |
| 172 | end if | ||
| 173 | |||
| 174 | |||
| 175 |
2/2✓ Branch 0 taken 272 times.
✓ Branch 1 taken 46 times.
|
318 | do i=-amax,amax+1,1 |
| 176 | 272 | vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32) | |
| 177 |
2/2✓ Branch 0 taken 1624 times.
✓ Branch 1 taken 272 times.
|
1942 | do j=-bmax,bmax+1,1 |
| 178 | 1624 | vtmp1(2) = this%spec(is)%atom(ia,2) + real(j, real32) | |
| 179 |
2/2✓ Branch 0 taken 6488 times.
✓ Branch 1 taken 1624 times.
|
8384 | do k=-cmax,cmax+1,1 |
| 180 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 6442 times.
|
6488 | if( i .eq. 0 .and. j .eq. 0 .and. k .eq. 0 ) cycle |
| 181 | 6442 | vtmp1(3) = this%spec(is)%atom(ia,3) + real(k, real32) | |
| 182 |
2/2✓ Branch 1 taken 1555 times.
✓ Branch 2 taken 4887 times.
|
6442 | if( get_distance_from_unit_cell(vtmp1, this%lat) .le. & |
| 183 | 1624 | max_bondlength ) then | |
| 184 | ! add the image to the list | ||
| 185 | 1555 | num_images = num_images + 1 | |
| 186 |
2/2✓ Branch 0 taken 4665 times.
✓ Branch 1 taken 1555 times.
|
6220 | image_species%atom(num_images,:3) = vtmp1 |
| 187 | end if | ||
| 188 | end do | ||
| 189 | end do | ||
| 190 | end do | ||
| 191 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
|
46 | if( num_images .eq. this%image_spec(is)%num ) return |
| 192 | |||
| 193 | |||
| 194 | 46 | this%image_spec(is)%num = num_images | |
| 195 |
2/4✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
|
46 | if(allocated(this%image_spec(is)%atom)) deallocate(this%image_spec(is)%atom) |
| 196 | ✗ | allocate(this%image_spec(is)%atom( & | |
| 197 | num_images, & | ||
| 198 | size(image_species%atom,2) & | ||
| 199 |
9/18✓ Branch 0 taken 46 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 46 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 46 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 46 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 46 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 46 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 46 times.
|
46 | ) ) |
| 200 | this%image_spec(is)%atom(:,:) = & | ||
| 201 |
4/4✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 58680 times.
✓ Branch 3 taken 138 times.
|
58864 | image_species%atom(:num_images,:) |
| 202 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
|
46 | deallocate(image_species%atom) |
| 203 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 46 times.
|
102 | this%num_images = sum( this%image_spec(:)%num ) |
| 204 | |||
| 205 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 46 times.
|
46 | end subroutine update_images |
| 206 | !############################################################################### | ||
| 207 | |||
| 208 | |||
| 209 | !############################################################################### | ||
| 210 | 78401 | function get_distance_from_unit_cell( & | |
| 211 | point, lattice, closest_point, is_cartesian& | ||
| 212 | ) result(distance) | ||
| 213 | !! Get the distance of a point from the unit cell. | ||
| 214 | implicit none | ||
| 215 | |||
| 216 | ! Arguments | ||
| 217 | real(real32), intent(in) :: point(3) | ||
| 218 | !! Query point. | ||
| 219 | real(real32), intent(in) :: lattice(3,3) | ||
| 220 | !! Lattice vectors. | ||
| 221 | real(real32), intent(out), optional :: closest_point(3) | ||
| 222 | !! Closest point on the unit cell surface. | ||
| 223 | logical, optional, intent(in) :: is_cartesian | ||
| 224 | !! Boolean whether the point is in cartesian coordinates. | ||
| 225 | real(real32) :: distance | ||
| 226 | !! Distance of the point from the unit cell. | ||
| 227 | |||
| 228 | ! Local variables | ||
| 229 | integer :: i, j, k | ||
| 230 | !! Loop indices. | ||
| 231 | real(real32), dimension(3) :: point_ | ||
| 232 | !! Point in cartesian coordinates. | ||
| 233 | real(real32), dimension(3,3) :: inverse_lattice | ||
| 234 | !! Inverse of the lattice vectors. | ||
| 235 | real(real32), dimension(3) :: normal | ||
| 236 | !! Normal vector to the plane. | ||
| 237 | real(real32), dimension(3) :: plane_point | ||
| 238 | !! Point on the plane. | ||
| 239 | real(real32), dimension(3) :: projection, closest_point_ | ||
| 240 | !! Projection of the point onto the plane. | ||
| 241 | real(real32), dimension(3) :: inverse_projection | ||
| 242 | !! Inverse projection of the point onto the plane. | ||
| 243 | real(real32) :: min_distance | ||
| 244 | !! Minimum distance to the unit cell. | ||
| 245 | logical :: is_outside | ||
| 246 | !! Boolean whether the point is outside the unit cell. | ||
| 247 | integer, dimension(3) :: index_list | ||
| 248 | !! List of indices for the lattice vectors. | ||
| 249 | logical :: is_cartesian_ | ||
| 250 | !! Boolean whether the point is in cartesian coordinates. | ||
| 251 | |||
| 252 | |||
| 253 | !--------------------------------------------------------------------------- | ||
| 254 | ! check if the point is in cartesian coordinates | ||
| 255 | !--------------------------------------------------------------------------- | ||
| 256 | 78401 | is_outside = .false. | |
| 257 | 78401 | is_cartesian_ = .false. | |
| 258 | 78401 | index_list = [1, 2, 3] | |
| 259 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78401 times.
|
78401 | if(present(is_cartesian)) is_cartesian_ = is_cartesian |
| 260 | 78401 | inverse_lattice = inverse_3x3( lattice ) | |
| 261 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78401 times.
|
78401 | if(is_cartesian_) then |
| 262 | ! Convert point to fractional coordinates | ||
| 263 | ! point_ = matmul(LUinv(lattice), point) | ||
| 264 | ✗ | point_ = point | |
| 265 | else | ||
| 266 | 78401 | point_ = matmul(point, lattice) | |
| 267 | end if | ||
| 268 | |||
| 269 | 78401 | min_distance = huge(1._real32) | |
| 270 | |||
| 271 | ! get projection of point onto each face of the lattice | ||
| 272 | ! get the length of the projection vector | ||
| 273 | ! if negative, then the point is inside the unit cell | ||
| 274 | ! if positive, then the point is outside the unit cell | ||
| 275 | ! if the projection falls outside of the cell edges, use edge or corner | ||
| 276 | ! distances | ||
| 277 |
2/2✓ Branch 0 taken 235203 times.
✓ Branch 1 taken 78401 times.
|
313604 | face_loop: do i = 1, 3 |
| 278 |
2/2✓ Branch 1 taken 705609 times.
✓ Branch 2 taken 235203 times.
|
940812 | index_list = cshift(index_list, 1) |
| 279 | 235203 | plane_point = 0._real32 | |
| 280 |
2/2✓ Branch 0 taken 470406 times.
✓ Branch 1 taken 235203 times.
|
784010 | direction_loop: do j = 1, 2 |
| 281 | normal = (-1._real32)**j * cross( & | ||
| 282 | [ lattice(index_list(2),:3) ], & | ||
| 283 | [ lattice(index_list(3),:3) ] & | ||
| 284 |
10/10✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
✓ Branch 2 taken 1411218 times.
✓ Branch 3 taken 470406 times.
✓ Branch 4 taken 1411218 times.
✓ Branch 5 taken 470406 times.
✓ Branch 6 taken 1411218 times.
✓ Branch 7 taken 470406 times.
✓ Branch 9 taken 1411218 times.
✓ Branch 10 taken 470406 times.
|
7526496 | ) |
| 285 |
8/8✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
✓ Branch 2 taken 503046 times.
✓ Branch 3 taken 908172 times.
✓ Branch 4 taken 470406 times.
✓ Branch 5 taken 32640 times.
✓ Branch 6 taken 1411218 times.
✓ Branch 7 taken 470406 times.
|
3292842 | normal = normal / norm2(normal) |
| 286 | 470406 | projection = project_point_onto_plane(point_, plane_point, normal) | |
| 287 | |||
| 288 | ! check if point minus projection is negative | ||
| 289 | ! if so, it is on the wrong side of the plane and should be ignored | ||
| 290 |
4/4✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
✓ Branch 2 taken 263436 times.
✓ Branch 3 taken 206970 times.
|
1881624 | if( dot_product(point_ - projection, normal) .lt. 0._real32 )then |
| 291 |
2/2✓ Branch 0 taken 790308 times.
✓ Branch 1 taken 263436 times.
|
1053744 | plane_point = plane_point + lattice(index_list(1),:) |
| 292 | 263436 | cycle direction_loop | |
| 293 | end if | ||
| 294 | 206970 | is_outside = .true. | |
| 295 | |||
| 296 | ! check if projection is outside the surface | ||
| 297 | |||
| 298 | 206970 | inverse_projection = matmul(projection, inverse_lattice) | |
| 299 | if( & | ||
| 300 |
10/10✓ Branch 0 taken 485794 times.
✓ Branch 1 taken 91566 times.
✓ Branch 2 taken 115404 times.
✓ Branch 3 taken 370390 times.
✓ Branch 4 taken 451515 times.
✓ Branch 5 taken 63577 times.
✓ Branch 6 taken 143393 times.
✓ Branch 7 taken 308122 times.
✓ Branch 8 taken 196663 times.
✓ Branch 9 taken 10307 times.
|
885482 | any( inverse_projection .lt. 0._real32 ) .or. & |
| 301 | any( inverse_projection .gt. 1._real32 ) & | ||
| 302 | ) then | ||
| 303 | ! projection is outside the surface | ||
| 304 | ! check if the projection is outside the edges | ||
| 305 | ! if it is, then the closest point is the edge or corner | ||
| 306 | ! if it is not, then the closest point is the projection | ||
| 307 |
2/2✓ Branch 0 taken 589989 times.
✓ Branch 1 taken 196663 times.
|
786652 | do k = 1, 3 |
| 308 |
2/2✓ Branch 0 taken 139100 times.
✓ Branch 1 taken 450889 times.
|
786652 | if( inverse_projection(k) .lt. 0._real32 ) then |
| 309 | 139100 | inverse_projection(k) = 0._real32 | |
| 310 |
2/2✓ Branch 0 taken 185276 times.
✓ Branch 1 taken 265613 times.
|
450889 | else if( inverse_projection(k) .gt. 1._real32 ) then |
| 311 | 185276 | inverse_projection(k) = 1._real32 | |
| 312 | end if | ||
| 313 | end do | ||
| 314 | end if | ||
| 315 | 206970 | projection = matmul(inverse_projection, lattice) | |
| 316 |
6/6✓ Branch 0 taken 620910 times.
✓ Branch 1 taken 206970 times.
✓ Branch 2 taken 538482 times.
✓ Branch 3 taken 82428 times.
✓ Branch 4 taken 314137 times.
✓ Branch 5 taken 224345 times.
|
827880 | distance = norm2(point_ - projection) |
| 317 |
2/2✓ Branch 0 taken 81885 times.
✓ Branch 1 taken 125085 times.
|
206970 | if( distance .lt. min_distance ) then |
| 318 | 81885 | min_distance = distance | |
| 319 | 81885 | closest_point_ = projection | |
| 320 | end if | ||
| 321 | |||
| 322 | !! makes it apply to the next iteration | ||
| 323 |
2/2✓ Branch 0 taken 620910 times.
✓ Branch 1 taken 206970 times.
|
1063083 | plane_point = plane_point + lattice(index_list(1),:) |
| 324 | end do direction_loop | ||
| 325 | end do face_loop | ||
| 326 | |||
| 327 |
2/2✓ Branch 0 taken 78385 times.
✓ Branch 1 taken 16 times.
|
78401 | if( is_outside ) then |
| 328 | 78385 | distance = min_distance | |
| 329 | else | ||
| 330 | 16 | distance = 0._real32 | |
| 331 | end if | ||
| 332 | |||
| 333 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78401 times.
|
78401 | if( present(closest_point) ) then |
| 334 | ✗ | if(is_cartesian_) then | |
| 335 | ✗ | closest_point = closest_point_ | |
| 336 | else | ||
| 337 | ✗ | closest_point = matmul(closest_point_, inverse_lattice) | |
| 338 | end if | ||
| 339 | end if | ||
| 340 | |||
| 341 | 78401 | end function get_distance_from_unit_cell | |
| 342 | !############################################################################### | ||
| 343 | |||
| 344 | |||
| 345 | !############################################################################### | ||
| 346 |
1/2✓ Branch 0 taken 470406 times.
✗ Branch 1 not taken.
|
470406 | function project_point_onto_plane(point, plane_point, normal) result(output) |
| 347 | !! Project a point onto a plane. | ||
| 348 | implicit none | ||
| 349 | |||
| 350 | ! Arguments | ||
| 351 | real(real32), dimension(3), intent(in) :: point | ||
| 352 | !! Point to project. | ||
| 353 | real(real32), dimension(3), intent(in) :: plane_point | ||
| 354 | !! Point on the plane. | ||
| 355 | real(real32), dimension(3), intent(in) :: normal | ||
| 356 | !! Normal vector to the plane. | ||
| 357 | real(real32), dimension(3) :: output | ||
| 358 | !! Projected point. | ||
| 359 | |||
| 360 | ! Local variables | ||
| 361 | real(real32) :: distance | ||
| 362 | !! Distance of the point from the plane. | ||
| 363 | real(real32), dimension(3) :: vector_to_plane | ||
| 364 | !! Vector from the point to the plane. | ||
| 365 | |||
| 366 |
2/2✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
|
1881624 | vector_to_plane = point - plane_point |
| 367 | |||
| 368 | distance = & | ||
| 369 |
4/4✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
✓ Branch 2 taken 1411218 times.
✓ Branch 3 taken 470406 times.
|
3292842 | dot_product(vector_to_plane, normal) / dot_product(normal, normal) |
| 370 | |||
| 371 |
2/2✓ Branch 0 taken 1411218 times.
✓ Branch 1 taken 470406 times.
|
1881624 | output = point - distance * normal |
| 372 | |||
| 373 | 470406 | end function project_point_onto_plane | |
| 374 | !############################################################################### | ||
| 375 | |||
| 376 | ✗ | end module raffle__geom_extd | |
| 377 |