GCC Code Coverage Report


Directory: src/fortran/lib/
File: mod_evaluator.f90
Date: 2025-06-15 07:27:34
Exec Total Coverage
Lines: 137 144 95.1%
Functions: 0 0 -%
Branches: 215 280 76.8%

Line Branch Exec Source
1 module raffle__evaluator
2 !! Module to build viability map of a structure
3 !!
4 !! This module handles the viability map for a structure, which is a map of
5 !! the system with each point in the map representing the suitability of
6 !! that point for a new atom. The map is built by checking the bond lengths,
7 !! bond angles and dihedral angles between the test point and all atoms.
8 use raffle__constants, only: real32, tau, pi
9 use raffle__misc_linalg, only: modu, get_angle, get_improper_dihedral_angle
10 use raffle__geom_extd, only: extended_basis_type
11 use raffle__distribs_container, only: distribs_container_type
12 implicit none
13
14
15 private
16 public :: evaluate_point
17
18
19 contains
20
21 !###############################################################################
22 2305359 pure function evaluate_point( distribs_container, &
23 position, species, basis, &
24
1/2
✓ Branch 0 taken 2305359 times.
✗ Branch 1 not taken.
2305359 radius_list &
25 ) result(output)
26 !! Return the viability of a point in a basis for a specified species
27 !!
28 !! This function evaluates the viability of a point in a basis for a
29 !! specified species. The viability is determined by the bond lengths,
30 !! bond angles and dihedral angles between the test point and all atoms.
31 implicit none
32
33 ! Arguments
34 integer, intent(in) :: species
35 !! Index of the query element.
36 type(distribs_container_type), intent(in) :: distribs_container
37 !! Distribution function (gvector) container.
38 type(extended_basis_type), intent(in) :: basis
39 !! Basis of the system.
40 real(real32), dimension(3), intent(in) :: position
41 !! Position of the test point.
42 real(real32), dimension(:), intent(in) :: radius_list
43 !! List of radii for each pair of elements.
44 real(real32) :: output
45 !! Suitability of the test point.
46
47 ! Local variables
48 integer :: i, is, js, ia, ja, ia_end, ja_start, is_end
49 !! Loop counters.
50 integer :: element_idx
51 !! Index of the query element.
52 integer :: num_2body, num_3body, num_4body
53 !! Number of 2-, 3- and 4-body interactions.
54 real(real32) :: viability_2body
55 !! Viability of the test point for 2-body interactions.
56 real(real32) :: viability_3body, viability_4body
57 !! Viability of the test point for 3- and 4-body interactions.
58 real(real32) :: bondlength, rtmp1, min_distance, min_from_max_cutoff
59 !! Temporary variables.
60 logical :: has_4body
61 !! Boolean whether the system has 4-body interactions.
62 real(real32) :: cos_scale1, cos_scale2
63 !! Cosine scales for the 3- and 4-body interactions.
64 real(real32), dimension(3) :: position_1
65 !! Cartesian coordinates of the test point.
66 real(real32), dimension(4) :: position_2
67 !! Cartesian coordinates of the second atom and its cutoff weighting.
68 real(real32), dimension(4) :: tolerances
69 !! Tolerance for the distance between atoms for 3- and 4-body.
70 2305359 integer, dimension(:,:), allocatable :: pair_index
71 !! Index of element pairs.
72
6/6
✓ Branch 0 taken 6916077 times.
✓ Branch 1 taken 2305359 times.
✓ Branch 2 taken 20748231 times.
✓ Branch 3 taken 6916077 times.
✓ Branch 4 taken 6916077 times.
✓ Branch 5 taken 2305359 times.
36885744 type(extended_basis_type) :: neighbour_basis
73 !! Basis of the neighbouring atoms.
74
75
76 ! Initialisation
77 2305359 output = 0._real32
78 2305359 viability_2body = 0._real32
79 2305359 min_distance = distribs_container%cutoff_max(1)
80 2305359 min_from_max_cutoff = 0._real32
81
82
83 !---------------------------------------------------------------------------
84 ! get list of element pair indices
85 ! (i.e. the index for bond_info for each element pair)
86 !---------------------------------------------------------------------------
87
13/22
✓ Branch 0 taken 2305359 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2305359 times.
✓ Branch 4 taken 2305359 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2305359 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2305359 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2305359 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2305359 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2305359 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2305359 times.
✓ Branch 21 taken 6200739 times.
✓ Branch 22 taken 2305359 times.
✓ Branch 23 taken 17886879 times.
✓ Branch 24 taken 6200739 times.
26392977 allocate(pair_index(basis%nspec, basis%nspec), source = 0)
88
2/2
✓ Branch 0 taken 6200739 times.
✓ Branch 1 taken 2305359 times.
8506098 do is = 1, basis%nspec
89
2/2
✓ Branch 0 taken 17886879 times.
✓ Branch 1 taken 6200739 times.
26392977 do js = 1, basis%nspec
90 pair_index(is, js) = distribs_container%get_pair_index( &
91 basis%spec(is)%name, basis%spec(js)%name &
92 24087618 )
93 end do
94 end do
95
96
97 !---------------------------------------------------------------------------
98 ! loop over all atoms in the system
99 !---------------------------------------------------------------------------
100
13/24
✓ Branch 0 taken 2305359 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2305359 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2305359 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2305359 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2305359 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2305359 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2305359 times.
✓ Branch 17 taken 6200739 times.
✓ Branch 18 taken 2305359 times.
✓ Branch 19 taken 6200739 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6200739 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 6200739 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 6200739 times.
8506098 allocate(neighbour_basis%spec(basis%nspec))
101
13/24
✓ Branch 0 taken 2305359 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2305359 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2305359 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2305359 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2305359 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2305359 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2305359 times.
✓ Branch 17 taken 6200739 times.
✓ Branch 18 taken 2305359 times.
✓ Branch 19 taken 6200739 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6200739 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 6200739 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 6200739 times.
8506098 allocate(neighbour_basis%image_spec(basis%nspec))
102 2305359 neighbour_basis%nspec = basis%nspec
103
4/4
✓ Branch 0 taken 6916077 times.
✓ Branch 1 taken 2305359 times.
✓ Branch 2 taken 20748231 times.
✓ Branch 3 taken 6916077 times.
29969667 neighbour_basis%lat = basis%lat
104 2305359 num_2body = 0
105
2/2
✓ Branch 0 taken 3427397 times.
✓ Branch 1 taken 784705 times.
4212102 species_loop: do is = 1, basis%nspec
106 allocate(neighbour_basis%spec(is)%atom( &
107 basis%spec(is)%num+basis%image_spec(is)%num, 4 &
108
8/16
✓ Branch 0 taken 3427397 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3427397 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3427397 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3427397 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3427397 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3427397 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3427397 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3427397 times.
3427397 ) )
109 allocate(neighbour_basis%image_spec(is)%atom( &
110 basis%spec(is)%num+basis%image_spec(is)%num, 4 &
111
8/16
✓ Branch 0 taken 3427397 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3427397 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3427397 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3427397 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3427397 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3427397 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3427397 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3427397 times.
3427397 ) )
112 3427397 neighbour_basis%spec(is)%num = 0
113 3427397 neighbour_basis%image_spec(is)%num = 0
114 tolerances = distribs_container%radius_distance_tol * &
115
2/2
✓ Branch 0 taken 13709588 times.
✓ Branch 1 taken 3427397 times.
17136985 radius_list(pair_index(species,is))
116 3427397 tolerances(1) = max( distribs_container%cutoff_min(1), tolerances(1) )
117 3427397 tolerances(3) = max( distribs_container%cutoff_min(1), tolerances(3) )
118 3427397 tolerances(2) = min( distribs_container%cutoff_max(1), tolerances(2) )
119 3427397 tolerances(4) = min( distribs_container%cutoff_max(1), tolerances(4) )
120 3427397 cos_scale1 = tau / ( tolerances(2) - tolerances(1) )
121 3427397 cos_scale2 = tau / ( tolerances(4) - tolerances(3) )
122 !------------------------------------------------------------------------
123 ! 2-body map
124 ! check bondlength between test point and all other atoms
125 !------------------------------------------------------------------------
126
2/2
✓ Branch 0 taken 14237097 times.
✓ Branch 1 taken 2946558 times.
17183655 atom_loop: do ia = 1, basis%spec(is)%num
127 ! Check if the atom is in the ignore list
128 ! If it is, skip the atom.
129
2/2
✓ Branch 0 taken 4602854 times.
✓ Branch 1 taken 9634243 times.
15388600 if(.not.basis%spec(is)%atom_mask(ia)) cycle atom_loop
130 2946558 associate( position_store => [ basis%spec(is)%atom(ia,1:3) ] )
131
8/8
✓ Branch 0 taken 28902729 times.
✓ Branch 1 taken 9634243 times.
✓ Branch 3 taken 28902729 times.
✓ Branch 4 taken 9634243 times.
✓ Branch 5 taken 28175063 times.
✓ Branch 6 taken 727666 times.
✓ Branch 7 taken 15784631 times.
✓ Branch 8 taken 12390432 times.
67439701 bondlength = norm2( matmul(position - position_store, basis%lat) )
132
2/2
✓ Branch 0 taken 1151503 times.
✓ Branch 1 taken 8482740 times.
9634243 if( bondlength .gt. distribs_container%cutoff_max(1) ) &
133 1151503 cycle atom_loop
134
2/2
✓ Branch 0 taken 480839 times.
✓ Branch 1 taken 8001901 times.
8482740 if( bondlength .lt. tolerances(1) )then
135 ! If the bond length is less than the minimum allowed bond,
136 ! return 0 (i.e. the point is not viable).
137 480839 return
138
2/2
✓ Branch 0 taken 1890874 times.
✓ Branch 1 taken 6111027 times.
8001901 elseif( bondlength .le. tolerances(2) )then
139 ! If the bond length is within the tolerance of the covalent
140 ! radius of the pair, add the atom to the list of
141 ! atoms to consider for 3-body interactions.
142 1890874 neighbour_basis%spec(is)%num = neighbour_basis%spec(is)%num + 1
143 neighbour_basis%spec(is)%atom( &
144 neighbour_basis%spec(is)%num,:3 &
145
2/2
✓ Branch 1 taken 5672622 times.
✓ Branch 2 taken 1890874 times.
7563496 ) = matmul(position_store, basis%lat)
146 neighbour_basis%spec(is)%atom( &
147 neighbour_basis%spec(is)%num,4 &
148 ) = 0.5_real32 * abs( 1._real32 - &
149 1890874 cos( cos_scale1 * ( bondlength - tolerances(1) ) ) )
150 end if
151
152
2/2
✓ Branch 0 taken 3846287 times.
✓ Branch 1 taken 4155614 times.
8001901 if( bondlength .ge. tolerances(3) .and. &
153 bondlength .le. tolerances(4) &
154 )then
155 ! If the bond length is within the min and max allowed size,
156 ! add the atom to the list of atoms to consider for 4-body.
157 neighbour_basis%image_spec(is)%num = &
158 3846287 neighbour_basis%image_spec(is)%num + 1
159 neighbour_basis%image_spec(is)%atom( &
160 neighbour_basis%image_spec(is)%num,:3 &
161
2/2
✓ Branch 1 taken 11538861 times.
✓ Branch 2 taken 3846287 times.
15385148 ) = matmul(position_store, basis%lat)
162 neighbour_basis%image_spec(is)%atom( &
163 neighbour_basis%image_spec(is)%num,4 &
164 ) = 0.5_real32 * abs( 1._real32 - &
165 3846287 cos( cos_scale2 * ( bondlength - tolerances(3) ) ) )
166 end if
167
168 !------------------------------------------------------------------
169 ! Add the contribution of the bond length to the viability
170 !------------------------------------------------------------------
171
2/2
✓ Branch 0 taken 3733478 times.
✓ Branch 1 taken 4268423 times.
8001901 if(bondlength - tolerances(1) .lt. min_distance)then
172 3733478 min_distance = bondlength - tolerances(1)
173 end if
174
2/2
✓ Branch 0 taken 3951421 times.
✓ Branch 1 taken 4050480 times.
8001901 if(distribs_container%cutoff_max(1) - bondlength .gt. &
175 min_from_max_cutoff)then
176 min_from_max_cutoff = distribs_container%cutoff_max(1) - &
177 3951421 bondlength
178 end if
179 viability_2body = viability_2body + &
180 evaluate_2body_contributions( &
181 distribs_container, bondlength, pair_index(species,is) &
182 8001901 )
183
4/4
✓ Branch 0 taken 28902729 times.
✓ Branch 1 taken 9634243 times.
✓ Branch 2 taken 28902729 times.
✓ Branch 3 taken 9634243 times.
75441602 num_2body = num_2body + 1
184 end associate
185 end do atom_loop
186
187 !------------------------------------------------------------------------
188 ! Repeat the process for the image atoms.
189 ! i.e. atoms that are not in the unit cell but are within the cutoff
190 ! distance.
191 !------------------------------------------------------------------------
192
2/2
✓ Branch 0 taken 237757591 times.
✓ Branch 1 taken 1906743 times.
239664334 image_loop: do ia = 1, basis%image_spec(is)%num, 1
193 1906743 associate( position_store => [ basis%image_spec(is)%atom(ia,1:3) ] )
194
8/8
✓ Branch 0 taken 713272773 times.
✓ Branch 1 taken 237757591 times.
✓ Branch 3 taken 713272773 times.
✓ Branch 4 taken 237757591 times.
✓ Branch 5 taken 708301084 times.
✓ Branch 6 taken 4971689 times.
✓ Branch 7 taken 415043930 times.
✓ Branch 8 taken 293257154 times.
1664303137 bondlength = norm2( matmul(position - position_store, basis%lat) )
195
2/2
✓ Branch 0 taken 177290505 times.
✓ Branch 1 taken 60467086 times.
237757591 if( bondlength .gt. distribs_container%cutoff_max(1) ) &
196 177290505 cycle image_loop
197
2/2
✓ Branch 0 taken 1039815 times.
✓ Branch 1 taken 59427271 times.
60467086 if( bondlength .lt. tolerances(1) )then
198 1039815 return
199
2/2
✓ Branch 0 taken 3435746 times.
✓ Branch 1 taken 55991525 times.
59427271 elseif( bondlength .le. tolerances(2) )then
200 3435746 neighbour_basis%spec(is)%num = neighbour_basis%spec(is)%num + 1
201 neighbour_basis%spec(is)%atom( &
202 neighbour_basis%spec(is)%num,:3 &
203
2/2
✓ Branch 1 taken 10307238 times.
✓ Branch 2 taken 3435746 times.
13742984 ) = matmul(position_store, basis%lat)
204 neighbour_basis%spec(is)%atom( &
205 neighbour_basis%spec(is)%num,4 &
206 ) = 0.5_real32 * ( 1._real32 - &
207 3435746 cos( cos_scale1 * ( bondlength - tolerances(1) ) ) )
208 end if
209
210
2/2
✓ Branch 0 taken 24142692 times.
✓ Branch 1 taken 35284579 times.
59427271 if( bondlength .ge. tolerances(3) .and. &
211 bondlength .le. tolerances(4) &
212 )then
213 neighbour_basis%image_spec(is)%num = &
214 24142692 neighbour_basis%image_spec(is)%num + 1
215 neighbour_basis%image_spec(is)%atom( &
216 neighbour_basis%image_spec(is)%num,:3 &
217
2/2
✓ Branch 1 taken 72428076 times.
✓ Branch 2 taken 24142692 times.
96570768 ) = matmul(position_store, basis%lat)
218 neighbour_basis%image_spec(is)%atom( &
219 neighbour_basis%image_spec(is)%num,4 &
220 ) = 0.5_real32 * abs( 1._real32 - &
221 24142692 cos( cos_scale2 * ( bondlength - tolerances(3) ) ) )
222 end if
223
224 !------------------------------------------------------------------
225 ! Add the contribution of the bond length to the viability
226 !------------------------------------------------------------------
227
2/2
✓ Branch 0 taken 2399324 times.
✓ Branch 1 taken 57027947 times.
59427271 if(bondlength - tolerances(1) .lt. min_distance)then
228 2399324 min_distance = bondlength - tolerances(1)
229 end if
230
2/2
✓ Branch 0 taken 2550898 times.
✓ Branch 1 taken 56876373 times.
59427271 if(distribs_container%cutoff_max(1) - bondlength .gt. &
231 min_from_max_cutoff)then
232 min_from_max_cutoff = distribs_container%cutoff_max(1) - &
233 2550898 bondlength
234 end if
235 viability_2body = viability_2body + &
236 evaluate_2body_contributions( &
237 distribs_container, bondlength, pair_index(species,is) &
238 59427271 )
239
4/4
✓ Branch 0 taken 713272773 times.
✓ Branch 1 taken 237757591 times.
✓ Branch 2 taken 713272773 times.
✓ Branch 3 taken 237757591 times.
1723730408 num_2body = num_2body + 1
240 end associate
241 end do image_loop
242 !------------------------------------------------------------------------
243 ! DEVELOPER TOOL: This conditional allows the user to retrieve
244 ! results closer to first arXiv paper release.
245 ! It is not recommended to use this option for production runs and is
246 ! only here for testing purposes.
247 !------------------------------------------------------------------------
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1906743 times.
2691448 if(.not.distribs_container%smooth_viability)then
249 neighbour_basis%spec(is)%atom(1:neighbour_basis%spec(is)%num,4) = &
250 1._real32
251 neighbour_basis%image_spec(is)%atom( &
252 1:neighbour_basis%image_spec(is)%num,4 &
253 ) = 1._real32
254 end if
255 end do species_loop
256
2/2
✓ Branch 0 taken 1680355 times.
✓ Branch 1 taken 784705 times.
2465060 neighbour_basis%natom = sum(neighbour_basis%spec(:)%num)
257
258
259 !---------------------------------------------------------------------------
260 ! Normalise the bond length viability
261 !---------------------------------------------------------------------------
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 784705 times.
784705 if(num_2body.eq.0)then
263 ! This does not matter as, if there are no 2-body bonds, the point is
264 ! not meant to be included in the viability set.
265 ! The evaluator cannot comment on the viability of the point.
266 viability_2body = distribs_container%viability_2body_default
267 else
268 784705 viability_2body = viability_2body / real( num_2body, real32 )
269
2/2
✓ Branch 0 taken 361362 times.
✓ Branch 1 taken 423343 times.
784705 if(min_distance .lt. 0.25_real32 )then
270 viability_2body = viability_2body * 0.5_real32 * &
271 361362 ( 1._real32 - cos( pi * min_distance / 0.25_real32 ) )
272 end if
273
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 784705 times.
784705 if( min_from_max_cutoff .lt. 0.5_real32 )then
274 rtmp1 = 0.5_real32 * ( 1._real32 - &
275 cos( pi * min_from_max_cutoff / 0.5_real32 ) )
276 viability_2body = &
277 distribs_container%viability_2body_default * &
278 abs( 1._real32 - rtmp1 ) + &
279 rtmp1 * viability_2body
280 end if
281 end if
282
283
284 ! store 3-body viable atoms in neighbour_basis%spec
285 ! store 4-body viable atoms in neighbour_basis%image_spec
286 ! for 3-bdoy, just cycle over neighbour_basis%spec
287 ! for 4-body, cycle over neighbour_basis%spec for first two atoms,
288 ! and neighbour_basis%image_spec for the third atom
289 784705 num_3body = 0
290 784705 num_4body = 0
291 784705 viability_3body = 1._real32
292 784705 viability_4body = 1._real32
293 784705 position_1 = matmul(position, basis%lat)
294 784705 element_idx = distribs_container%element_map(species)
295
4/4
✓ Branch 0 taken 840232 times.
✓ Branch 1 taken 687 times.
✓ Branch 2 taken 784018 times.
✓ Branch 3 taken 56214 times.
840919 has_4body = any(neighbour_basis%image_spec(:)%num .gt. 0)
296 784705 is_end = 0
297
2/2
✓ Branch 0 taken 1680355 times.
✓ Branch 1 taken 784705 times.
2465060 do is = 1, neighbour_basis%nspec, 1
298
2/2
✓ Branch 0 taken 1074058 times.
✓ Branch 1 taken 606297 times.
2465060 if(neighbour_basis%spec(is)%num .gt. 0) is_end = is
299 end do
300
2/2
✓ Branch 0 taken 1193924 times.
✓ Branch 1 taken 784705 times.
1978629 do is = 1, is_end, 1
301
2/2
✓ Branch 0 taken 704492 times.
✓ Branch 1 taken 489432 times.
1193924 ia_end = neighbour_basis%spec(is)%num - merge( 1, 0, is .eq. is_end )
302
2/2
✓ Branch 0 taken 1985842 times.
✓ Branch 1 taken 1193924 times.
3964471 do ia = 1, ia_end, 1
303
2/2
✓ Branch 0 taken 7943368 times.
✓ Branch 1 taken 1985842 times.
9929210 position_2 = neighbour_basis%spec(is)%atom(ia,1:4)
304 !---------------------------------------------------------------------
305 ! 3-body map
306 ! check bondangle between test point and all other atoms
307 !---------------------------------------------------------------------
308 viability_3body = viability_3body * max( &
309 0._real32, &
310 evaluate_3body_contributions( distribs_container, &
311 position_1, &
312 position_2, &
313 neighbour_basis, element_idx, [is, ia] &
314
2/2
✓ Branch 0 taken 3971684 times.
✓ Branch 1 taken 1985842 times.
5957526 ) )
315 1985842 num_3body = num_3body + 1
316
1/2
✓ Branch 0 taken 1985842 times.
✗ Branch 1 not taken.
3179766 if( has_4body )then
317
2/2
✓ Branch 0 taken 3717184 times.
✓ Branch 1 taken 1985842 times.
5703026 do js = is, neighbour_basis%nspec
318
2/2
✓ Branch 0 taken 1985842 times.
✓ Branch 1 taken 1731342 times.
3717184 ja_start = merge(ia + 1, 1, js .eq. is)
319
2/2
✓ Branch 0 taken 5694882 times.
✓ Branch 1 taken 3717184 times.
11397908 do ja = ja_start, neighbour_basis%spec(js)%num, 1
320 !------------------------------------------------------------
321 ! 4-body map
322 ! check improperdihedral angle between test point and all
323 ! other atoms
324 !------------------------------------------------------------
325 rtmp1 = max( &
326 0._real32, &
327 evaluate_4body_contributions( distribs_container, &
328 position_1, &
329 position_2, &
330 [ neighbour_basis%spec(js)%atom(ja,1:4) ], &
331 neighbour_basis, element_idx &
332
4/4
✓ Branch 0 taken 22779528 times.
✓ Branch 1 taken 5694882 times.
✓ Branch 2 taken 22779528 times.
✓ Branch 3 taken 5694882 times.
51253938 ) )
333 5694882 viability_4body = viability_4body * rtmp1
334 9412066 num_4body = num_4body + 1
335 end do
336 end do
337 end if
338 end do
339 end do
340
341
342 !---------------------------------------------------------------------------
343 ! Normalise the angular viabilities
344 !---------------------------------------------------------------------------
345
2/2
✓ Branch 0 taken 581532 times.
✓ Branch 1 taken 203173 times.
784705 if(num_3body.gt.0)then
346 viability_3body = viability_3body ** ( &
347 1._real32 / real(num_3body,real32) &
348 581532 )
349 else
350 203173 viability_3body = distribs_container%viability_3body_default
351 end if
352
2/2
✓ Branch 0 taken 581532 times.
✓ Branch 1 taken 203173 times.
784705 if(num_4body.gt.0)then
353 viability_4body = viability_4body ** ( &
354 1._real32 / real(num_4body,real32) &
355 581532 )
356 else
357 203173 viability_4body = distribs_container%viability_4body_default
358 end if
359
360
361 !---------------------------------------------------------------------------
362 ! Combine the 2-, 3- and 4-body viabilities to get the overall viability
363 !---------------------------------------------------------------------------
364 784705 output = viability_2body * viability_3body * viability_4body
365
366
17/26
✓ Branch 0 taken 2305359 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2305359 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2305359 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6200739 times.
✓ Branch 7 taken 2305359 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 6200739 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 6200739 times.
✓ Branch 12 taken 3427397 times.
✓ Branch 13 taken 2773342 times.
✓ Branch 14 taken 2305359 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2305359 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 6200739 times.
✓ Branch 19 taken 2305359 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 6200739 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 6200739 times.
✓ Branch 24 taken 3427397 times.
✓ Branch 25 taken 2773342 times.
14706837 end function evaluate_point
367 !###############################################################################
368
369
370 !###############################################################################
371 67429172 pure function evaluate_2body_contributions( distribs_container, &
372 bondlength, pair_index &
373 ) result(output)
374 !! Return the contribution to the viability from 2-body interactions
375 implicit none
376
377 ! Arguments
378 type(distribs_container_type), intent(in) :: distribs_container
379 !! Distribution function (gvector) container.
380 real(real32), intent(in) :: bondlength
381 !! Bond length.
382 integer, intent(in) :: pair_index
383 !! Index of the element pair.
384 real(real32) :: output
385 !! Contribution to the viability.
386
387
388 output = distribs_container%gdf%df_2body( &
389 distribs_container%get_bin(bondlength, dim = 1), &
390 pair_index &
391 67429172 )
392
393 67429172 end function evaluate_2body_contributions
394 !###############################################################################
395
396
397 !###############################################################################
398 1985842 pure function evaluate_3body_contributions( distribs_container, &
399 position_1, position_2, basis, element_idx, current_idx &
400 ) result(output)
401 !! Return the contribution to the viability from 3-body interactions
402 implicit none
403
404 ! Arguments
405 type(distribs_container_type), intent(in) :: distribs_container
406 !! Distribution function (gvector) container.
407 real(real32), dimension(3), intent(in) :: position_1
408 !! Positions of the atom.
409 real(real32), dimension(4), intent(in) :: position_2
410 !! Positions of the second atom and its cutoff weighting.
411 type(extended_basis_type), intent(in) :: basis
412 !! Basis of the system.
413 integer, intent(in) :: element_idx
414 !! Index of the query element.
415 integer, dimension(2), intent(in) :: current_idx
416 !! Index of the 1st-atom query element.
417 real(real32) :: output
418 !! Contribution to the viability.
419
420 ! Local variables
421 integer :: js, ja
422 !! Loop indices.
423 integer :: bin
424 !! Bin for the distribution function.
425 integer :: num_3body_local
426 !! Number of 3-body interactions local to the current atom pair.
427 real(real32) :: power
428 !! Power for the contribution to the viability.
429 integer :: start_idx
430 !! Start index for the loop over the atoms.
431 real(real32) :: rtmp1, weight_2
432 !! Temporary variables.
433
434
435
2/2
✓ Branch 0 taken 3717184 times.
✓ Branch 1 taken 1985842 times.
5703026 num_3body_local = sum(basis%spec(current_idx(1):)%num) - current_idx(2)
436 1985842 output = 1._real32
437 1985842 power = 1._real32 / real( num_3body_local, real32 )
438 1985842 weight_2 = sqrt( position_2(4) )
439
2/2
✓ Branch 0 taken 3717184 times.
✓ Branch 1 taken 1985842 times.
5703026 species_loop: do js = current_idx(1), basis%nspec, 1
440
2/2
✓ Branch 0 taken 1985842 times.
✓ Branch 1 taken 1731342 times.
3717184 start_idx = merge(current_idx(2) + 1, 1, js .eq. current_idx(1))
441
2/2
✓ Branch 0 taken 5694882 times.
✓ Branch 1 taken 3717184 times.
11397908 atom_loop: do ja = start_idx, basis%spec(js)%num
442 bin = distribs_container%get_bin( &
443 get_angle( &
444 position_2(1:3), &
445 position_1, &
446 [ basis%spec(js)%atom(ja,1:3) ] &
447 ), dim = 2 &
448
4/4
✓ Branch 0 taken 17084646 times.
✓ Branch 1 taken 5694882 times.
✓ Branch 2 taken 17084646 times.
✓ Branch 3 taken 5694882 times.
39864174 )
449 5694882 rtmp1 = weight_2 * sqrt( basis%spec(js)%atom(ja,4) )
450 output = output * ( &
451 distribs_container%viability_3body_default * &
452 abs( 1._real32 - rtmp1 ) + &
453 rtmp1 * distribs_container%gdf%df_3body( bin, element_idx ) &
454 9412066 ) ** power
455 end do atom_loop
456 end do species_loop
457
458 1985842 end function evaluate_3body_contributions
459 !###############################################################################
460
461
462 !###############################################################################
463 5694882 pure function evaluate_4body_contributions( distribs_container, &
464 position_1, position_2, position_3, basis, element_idx ) result(output)
465 !! Return the contribution to the viability from 4-body interactions
466 implicit none
467
468 ! Arguments
469 type(distribs_container_type), intent(in) :: distribs_container
470 !! Distribution function (gvector) container.
471 real(real32), dimension(3), intent(in) :: position_1
472 !! Positions of the atoms.
473 real(real32), dimension(4), intent(in) :: position_2, position_3
474 type(extended_basis_type), intent(in) :: basis
475 !! Basis of the system.
476 integer, intent(in) :: element_idx
477 !! Index of the query element.
478 real(real32) :: output
479 !! Contribution to the viability.
480
481 ! Local variables
482 integer :: ks, ka
483 !! Loop indices.
484 integer :: bin
485 !! Bin for the distribution function.
486 integer :: num_4body_local
487 !! Number of 4-body interactions local to the current atom triplet.
488 real(real32) :: power
489 !! Power for the contribution to the viability.
490 real(real32) :: rtmp1, third, weight_2_3
491 !! Temporary variables.
492
493
494
2/2
✓ Branch 0 taken 11595936 times.
✓ Branch 1 taken 5694882 times.
17290818 num_4body_local = sum(basis%image_spec(:)%num)
495 5694882 output = 1._real32
496 5694882 power = 1._real32 / real( num_4body_local, real32 )
497 5694882 third = 1._real32 / 3._real32
498 5694882 weight_2_3 = ( position_2(4) * position_3(4) ) ** third
499
2/2
✓ Branch 0 taken 11595936 times.
✓ Branch 1 taken 5694882 times.
17290818 species_loop: do ks = 1, basis%nspec, 1
500
2/2
✓ Branch 0 taken 216393559 times.
✓ Branch 1 taken 11595936 times.
233684377 atom_loop: do ka = 1, basis%image_spec(ks)%num
501 bin = distribs_container%get_bin( &
502 get_improper_dihedral_angle( &
503 position_1, &
504 position_2(1:3), &
505 position_3(1:3), &
506 [ basis%image_spec(ks)%atom(ka,1:3) ] &
507 ), dim = 3 &
508
4/4
✓ Branch 0 taken 649180677 times.
✓ Branch 1 taken 216393559 times.
✓ Branch 2 taken 649180677 times.
✓ Branch 3 taken 216393559 times.
1514754913 )
509 216393559 rtmp1 = weight_2_3 * ( basis%image_spec(ks)%atom(ka,4) ) ** third
510 output = output * ( &
511 distribs_container%viability_4body_default * &
512 abs( 1._real32 - rtmp1 ) + &
513 rtmp1 * distribs_container%gdf%df_4body( bin, element_idx ) &
514 227989495 ) ** power
515 end do atom_loop
516 end do species_loop
517
518 5694882 end function evaluate_4body_contributions
519 !###############################################################################
520
521 end module raffle__evaluator
522