GCC Code Coverage Report


Directory: src/fortran/lib/
File: mod_evaluator.f90
Date: 2025-04-05 12:17:58
Exec Total Coverage
Lines: 104 107 97.2%
Functions: 0 0 -%
Branches: 192 248 77.4%

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
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 2125783 pure function evaluate_point( distribs_container, &
23
1/2
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
2125783 position, species, basis, atom_ignore_list, &
24
1/2
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
2125783 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 integer, dimension(:,:), intent(in) :: atom_ignore_list
43 !! List of atoms to ignore (i.e. indices of atoms not yet placed).
44 real(real32), dimension(:), intent(in) :: radius_list
45 !! List of radii for each pair of elements.
46 real(real32) :: output
47 !! Suitability of the test point.
48
49 ! Local variables
50 integer :: i, is, js, ia, ja
51 !! Loop counters.
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
59 2125783 integer, dimension(:,:), allocatable :: pair_index
60 !! Index of element pairs.
61
6/6
✓ Branch 0 taken 6377349 times.
✓ Branch 1 taken 2125783 times.
✓ Branch 2 taken 19132047 times.
✓ Branch 3 taken 6377349 times.
✓ Branch 4 taken 6377349 times.
✓ Branch 5 taken 2125783 times.
34012528 type(extended_basis_type) :: neighbour_basis
62 !! Basis of the neighbouring atoms.
63
64
65 ! Initialisation
66 2125783 output = 0._real32
67 2125783 viability_2body = 0._real32
68
69
70 !---------------------------------------------------------------------------
71 ! get list of element pair indices
72 ! (i.e. the index for bond_info for each element pair)
73 !---------------------------------------------------------------------------
74
13/22
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2125783 times.
✓ Branch 4 taken 2125783 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2125783 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2125783 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2125783 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2125783 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2125783 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2125783 times.
✓ Branch 21 taken 6021163 times.
✓ Branch 22 taken 2125783 times.
✓ Branch 23 taken 17707303 times.
✓ Branch 24 taken 6021163 times.
25854249 allocate(pair_index(basis%nspec, basis%nspec), source = 0)
75
2/2
✓ Branch 0 taken 6021163 times.
✓ Branch 1 taken 2125783 times.
8146946 do is = 1, basis%nspec
76
2/2
✓ Branch 0 taken 17707303 times.
✓ Branch 1 taken 6021163 times.
25854249 do js = 1, basis%nspec
77 pair_index(is, js) = distribs_container%get_pair_index( &
78 basis%spec(is)%name, basis%spec(js)%name &
79 23728466 )
80 end do
81 end do
82
83
84 !---------------------------------------------------------------------------
85 ! loop over all atoms in the system
86 !---------------------------------------------------------------------------
87
11/20
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2125783 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2125783 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2125783 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2125783 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2125783 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2125783 times.
✓ Branch 17 taken 6021163 times.
✓ Branch 18 taken 2125783 times.
✓ Branch 19 taken 6021163 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6021163 times.
8146946 allocate(neighbour_basis%spec(basis%nspec))
88
11/20
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2125783 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2125783 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2125783 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2125783 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2125783 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2125783 times.
✓ Branch 17 taken 6021163 times.
✓ Branch 18 taken 2125783 times.
✓ Branch 19 taken 6021163 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6021163 times.
8146946 allocate(neighbour_basis%image_spec(basis%nspec))
89 2125783 neighbour_basis%nspec = basis%nspec
90
4/4
✓ Branch 0 taken 6377349 times.
✓ Branch 1 taken 2125783 times.
✓ Branch 2 taken 19132047 times.
✓ Branch 3 taken 6377349 times.
27635179 neighbour_basis%lat = basis%lat
91 2125783 num_2body = 0
92
2/2
✓ Branch 0 taken 3247821 times.
✓ Branch 1 taken 604099 times.
3851920 species_loop: do is = 1, basis%nspec
93 allocate(neighbour_basis%spec(is)%atom( &
94 basis%spec(is)%num+basis%image_spec(is)%num, &
95 size(basis%spec(is)%atom,2) &
96
9/18
✓ Branch 0 taken 3247821 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3247821 times.
✓ Branch 4 taken 3247821 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3247821 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3247821 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3247821 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3247821 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3247821 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 3247821 times.
3247821 ) )
97 allocate(neighbour_basis%image_spec(is)%atom( &
98 basis%spec(is)%num+basis%image_spec(is)%num, &
99 size(basis%spec(is)%atom,2) &
100
9/18
✓ Branch 0 taken 3247821 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3247821 times.
✓ Branch 4 taken 3247821 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3247821 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3247821 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3247821 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3247821 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3247821 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 3247821 times.
3247821 ) )
101 3247821 neighbour_basis%spec(is)%num = 0
102 3247821 neighbour_basis%image_spec(is)%num = 0
103 !------------------------------------------------------------------------
104 ! 2-body map
105 ! check bondlength between test point and all other atoms
106 !------------------------------------------------------------------------
107
2/2
✓ Branch 0 taken 10830511 times.
✓ Branch 1 taken 2766102 times.
13596613 atom_loop: do ia = 1, basis%spec(is)%num
108 ! Check if the atom is in the ignore list
109 ! If it is, skip the atom.
110
2/2
✓ Branch 0 taken 41070813 times.
✓ Branch 1 taken 7454564 times.
48525377 do i = 1, size(atom_ignore_list,dim=2), 1
111
6/6
✓ Branch 0 taken 66349826 times.
✓ Branch 1 taken 3375947 times.
✓ Branch 2 taken 37694866 times.
✓ Branch 3 taken 28654960 times.
✓ Branch 4 taken 3375947 times.
✓ Branch 5 taken 37694866 times.
77180337 if(all(atom_ignore_list(:,i).eq.[is,ia])) cycle atom_loop
112 end do
113 2766102 associate( position_store => [ basis%spec(is)%atom(ia,1:3) ] )
114
2/2
✓ Branch 0 taken 22363692 times.
✓ Branch 1 taken 7454564 times.
29818256 bondlength = modu( matmul(position - position_store, basis%lat) )
115
2/2
✓ Branch 0 taken 1066237 times.
✓ Branch 1 taken 6388327 times.
7454564 if( bondlength .gt. distribs_container%cutoff_max(1) ) &
116 1066237 cycle atom_loop
117
2/2
✓ Branch 0 taken 481719 times.
✓ Branch 1 taken 5906608 times.
6388327 if( bondlength .lt. ( &
118 radius_list(pair_index(species,is)) * &
119 distribs_container%radius_distance_tol(1) &
120 ) )then
121 ! If the bond length is less than the minimum allowed bond,
122 ! return 0 (i.e. the point is not viable).
123 481719 return
124
2/2
✓ Branch 0 taken 1509098 times.
✓ Branch 1 taken 4397510 times.
5906608 elseif( bondlength .le. ( &
125 radius_list(pair_index(species,is)) * &
126 distribs_container%radius_distance_tol(2) &
127 ) )then
128 ! If the bond length is within the tolerance of the covalent
129 ! radius of the pair, add the atom to the list of
130 ! atoms to consider for 3-body interactions.
131 1509098 neighbour_basis%spec(is)%num = neighbour_basis%spec(is)%num + 1
132 neighbour_basis%spec(is)%atom( &
133 neighbour_basis%spec(is)%num,:3 &
134
2/2
✓ Branch 1 taken 4527294 times.
✓ Branch 2 taken 1509098 times.
6036392 ) = matmul(position_store, basis%lat)
135 end if
136
137 if( bondlength .ge. &
138 ( &
139 radius_list(pair_index(species,is)) * &
140 distribs_container%radius_distance_tol(3) &
141
2/2
✓ Branch 0 taken 2633638 times.
✓ Branch 1 taken 3272970 times.
5906608 ) .and. &
142 bondlength .le. min( &
143 distribs_container%cutoff_max(1), &
144 ( &
145 radius_list(pair_index(species,is)) * &
146 distribs_container%radius_distance_tol(4) &
147 ) &
148 ) &
149 )then
150 ! If the bond length is within the min and max allowed size,
151 ! add the atom to the list of atoms to consider for 4-body.
152 neighbour_basis%image_spec(is)%num = &
153 2633638 neighbour_basis%image_spec(is)%num + 1
154 neighbour_basis%image_spec(is)%atom( &
155 neighbour_basis%image_spec(is)%num,:3 &
156
2/2
✓ Branch 1 taken 7900914 times.
✓ Branch 2 taken 2633638 times.
10534552 ) = matmul(position_store, basis%lat)
157 end if
158
159 !------------------------------------------------------------------
160 ! Add the contribution of the bond length to the viability
161 !------------------------------------------------------------------
162 viability_2body = viability_2body + &
163 evaluate_2body_contributions( &
164 distribs_container, bondlength, pair_index(species,is) &
165 5906608 )
166
4/4
✓ Branch 0 taken 22363692 times.
✓ Branch 1 taken 7454564 times.
✓ Branch 2 taken 22363692 times.
✓ Branch 3 taken 7454564 times.
58088556 num_2body = num_2body + 1
167 end associate
168 end do atom_loop
169
170 !------------------------------------------------------------------------
171 ! Repeat the process for the image atoms.
172 ! i.e. atoms that are not in the unit cell but are within the cutoff
173 ! distance.
174 !------------------------------------------------------------------------
175
2/2
✓ Branch 0 taken 163329306 times.
✓ Branch 1 taken 1726137 times.
165659542 image_loop: do ia = 1, basis%image_spec(is)%num, 1
176 1726137 associate( position_store => [ basis%image_spec(is)%atom(ia,1:3) ] )
177
2/2
✓ Branch 0 taken 489987918 times.
✓ Branch 1 taken 163329306 times.
653317224 bondlength = modu( matmul(position - position_store, basis%lat) )
178
2/2
✓ Branch 0 taken 123164847 times.
✓ Branch 1 taken 40164459 times.
163329306 if( bondlength .gt. distribs_container%cutoff_max(1) ) &
179 123164847 cycle image_loop
180
2/2
✓ Branch 0 taken 1039965 times.
✓ Branch 1 taken 39124494 times.
40164459 if( bondlength .lt. ( &
181 radius_list(pair_index(species,is)) * &
182 distribs_container%radius_distance_tol(1) &
183 ) )then
184 1039965 return
185
2/2
✓ Branch 0 taken 3093439 times.
✓ Branch 1 taken 36031055 times.
39124494 elseif( bondlength .le. ( &
186 radius_list(pair_index(species,is)) * &
187 distribs_container%radius_distance_tol(2) &
188 ) )then
189 3093439 neighbour_basis%spec(is)%num = neighbour_basis%spec(is)%num + 1
190 neighbour_basis%spec(is)%atom( &
191 neighbour_basis%spec(is)%num,:3 &
192
2/2
✓ Branch 1 taken 9280317 times.
✓ Branch 2 taken 3093439 times.
12373756 ) = matmul(position_store, basis%lat)
193 end if
194
195 if( bondlength .ge. &
196 ( &
197 radius_list(pair_index(species,is)) * &
198 distribs_container%radius_distance_tol(3) &
199
2/2
✓ Branch 0 taken 16921599 times.
✓ Branch 1 taken 22202895 times.
39124494 ) .and. &
200 bondlength .le. min( &
201 distribs_container%cutoff_max(1), &
202 ( &
203 radius_list(pair_index(species,is)) * &
204 distribs_container%radius_distance_tol(4) &
205 ) &
206 ) &
207 )then
208 neighbour_basis%image_spec(is)%num = &
209 16921599 neighbour_basis%image_spec(is)%num + 1
210 neighbour_basis%image_spec(is)%atom( &
211 neighbour_basis%image_spec(is)%num,:3 &
212
2/2
✓ Branch 1 taken 50764797 times.
✓ Branch 2 taken 16921599 times.
67686396 ) = matmul(position_store, basis%lat)
213 end if
214
215 !------------------------------------------------------------------
216 ! Add the contribution of the bond length to the viability
217 !------------------------------------------------------------------
218 viability_2body = viability_2body + &
219 evaluate_2body_contributions( &
220 distribs_container, bondlength, pair_index(species,is) &
221 39124494 )
222
4/4
✓ Branch 0 taken 489987918 times.
✓ Branch 1 taken 163329306 times.
✓ Branch 2 taken 489987918 times.
✓ Branch 3 taken 163329306 times.
1182429636 num_2body = num_2body + 1
223 end associate
224 end do image_loop
225 end do species_loop
226
2/2
✓ Branch 0 taken 1499749 times.
✓ Branch 1 taken 604099 times.
2103848 neighbour_basis%natom = sum(neighbour_basis%spec(:)%num)
227
228
229 !---------------------------------------------------------------------------
230 ! Normalise the bond length viability
231 !---------------------------------------------------------------------------
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 604099 times.
604099 if(num_2body.eq.0)then
233 ! This does not matter as, if there are no 2-body bonds, the point is
234 ! not meant to be included in the viability set.
235 ! The evaluator cannot comment on the viability of the point.
236 viability_2body = 0.5_real32
237 else
238 604099 viability_2body = viability_2body / real( num_2body, real32 )
239 end if
240
241
242
243 ! store 3-body viable atoms in neighbour_basis%spec
244 ! store 4-body viable atoms in neighbour_basis%image_spec
245 ! for 3-bdoy, just cycle over neighbour_basis%spec
246 ! for 4-body, cycle over neighbour_basis%spec for first two atoms,
247 ! and neighbour_basis%image_spec for the third atom
248 604099 num_3body = 0
249 604099 num_4body = 0
250 604099 viability_3body = 1._real32
251 604099 viability_4body = 1._real32
252
2/2
✓ Branch 0 taken 1499749 times.
✓ Branch 1 taken 604099 times.
2103848 do is = 1, neighbour_basis%nspec
253
2/2
✓ Branch 0 taken 1964697 times.
✓ Branch 1 taken 1499749 times.
4068545 do ia = 1, neighbour_basis%spec(is)%num
254 !---------------------------------------------------------------------
255 ! 3-body map
256 ! check bondangle between test point and all other atoms
257 !---------------------------------------------------------------------
258 associate( &
259 position_1 => matmul(position, basis%lat), &
260 position_2 => [neighbour_basis%spec(is)%atom(ia,1:3)] &
261 1499749 )
262
4/4
✓ Branch 0 taken 4019201 times.
✓ Branch 1 taken 1964697 times.
✓ Branch 2 taken 739 times.
✓ Branch 3 taken 1963958 times.
5983898 if(sum(basis%spec(is:)%num).eq.ia) cycle
263 1963958 num_3body = num_3body + 1
264 viability_3body = viability_3body * &
265 evaluate_3body_contributions( distribs_container, &
266 position_1, &
267 position_2, &
268 neighbour_basis, species, [is, ia] &
269
2/2
✓ Branch 0 taken 3927916 times.
✓ Branch 1 taken 1963958 times.
5891874 )
270
6/6
✓ Branch 1 taken 5894091 times.
✓ Branch 2 taken 1964697 times.
✓ Branch 3 taken 5894091 times.
✓ Branch 4 taken 1964697 times.
✓ Branch 5 taken 4018462 times.
✓ Branch 6 taken 1963958 times.
19735299 do js = is, neighbour_basis%nspec
271
2/2
✓ Branch 0 taken 8157543 times.
✓ Branch 1 taken 4018462 times.
14139963 do ja = 1, neighbour_basis%spec(js)%num
272
2/2
✓ Branch 0 taken 4199012 times.
✓ Branch 1 taken 3958531 times.
8157543 if(js.eq.is .and. ja.le.ia) cycle
273
4/6
✓ Branch 0 taken 5229838 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3958531 times.
✓ Branch 3 taken 1271307 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3958531 times.
5229838 if(all(neighbour_basis%image_spec(:)%num.eq.0))cycle
274 3958531 num_4body = num_4body + 1
275 !------------------------------------------------------------
276 ! 4-body map
277 ! check improperdihedral angle between test point and all
278 ! other atoms
279 !------------------------------------------------------------
280 viability_4body = viability_4body * &
281 evaluate_4body_contributions( distribs_container, &
282 position_1, &
283 position_2, &
284 [neighbour_basis%spec(js)%atom(ja,1:3)], &
285 neighbour_basis, species &
286
4/4
✓ Branch 0 taken 11875593 times.
✓ Branch 1 taken 3958531 times.
✓ Branch 2 taken 11875593 times.
✓ Branch 3 taken 3958531 times.
31728179 )
287 end do
288 end do
289 end associate
290 end do
291 end do
292
293
294 !---------------------------------------------------------------------------
295 ! Normalise the angular viabilities
296 !---------------------------------------------------------------------------
297
2/2
✓ Branch 0 taken 63926 times.
✓ Branch 1 taken 540173 times.
604099 if(num_3body.eq.0)then
298 63926 viability_3body = distribs_container%viability_3body_default
299 else
300 viability_3body = viability_3body ** ( &
301 1._real32 / real(num_3body,real32) &
302 540173 )
303 end if
304
2/2
✓ Branch 0 taken 158551 times.
✓ Branch 1 taken 445548 times.
604099 if(num_4body.eq.0)then
305 158551 viability_4body = distribs_container%viability_4body_default
306 else
307 viability_4body = viability_4body ** ( &
308 1._real32 / real(num_4body,real32) &
309 445548 )
310 end if
311
312
313 !---------------------------------------------------------------------------
314 ! Combine the 2-, 3- and 4-body viabilities to get the overall viability
315 !---------------------------------------------------------------------------
316 604099 output = viability_2body * viability_3body * viability_4body
317
318
13/18
✓ Branch 0 taken 2125783 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2125783 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2125783 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6021163 times.
✓ Branch 7 taken 2125783 times.
✓ Branch 8 taken 3247821 times.
✓ Branch 9 taken 2773342 times.
✓ Branch 10 taken 2125783 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2125783 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6021163 times.
✓ Branch 15 taken 2125783 times.
✓ Branch 16 taken 3247821 times.
✓ Branch 17 taken 2773342 times.
14168109 end function evaluate_point
319 !###############################################################################
320
321
322 !###############################################################################
323 45031102 pure function evaluate_2body_contributions( distribs_container, &
324 bondlength, pair_index &
325 ) result(output)
326 !! Return the contribution to the viability from 2-body interactions
327 implicit none
328
329 ! Arguments
330 type(distribs_container_type), intent(in) :: distribs_container
331 !! Distribution function (gvector) container.
332 real(real32), intent(in) :: bondlength
333 !! Bond length.
334 integer, intent(in) :: pair_index
335 !! Index of the element pair.
336 real(real32) :: output
337 !! Contribution to the viability.
338
339
340 output = distribs_container%gdf%df_2body( &
341 distribs_container%get_bin(bondlength, dim = 1), &
342 pair_index &
343 45031102 )
344
345 45031102 end function evaluate_2body_contributions
346 !###############################################################################
347
348
349 !###############################################################################
350 1963958 pure function evaluate_3body_contributions( distribs_container, &
351 position_1, position_2, basis, species, current_idx &
352 ) result(output)
353 !! Return the contribution to the viability from 3-body interactions
354 implicit none
355
356 ! Arguments
357 type(distribs_container_type), intent(in) :: distribs_container
358 !! Distribution function (gvector) container.
359 real(real32), dimension(3), intent(in) :: position_1, position_2
360 !! Positions of the atoms.
361 type(extended_basis_type), intent(in) :: basis
362 !! Basis of the system.
363 integer, intent(in) :: species
364 !! Index of the query element.
365 integer, dimension(2), intent(in) :: current_idx
366 !! Index of the 1st-atom query element.
367 real(real32) :: output
368 !! Contribution to the viability.
369
370 ! Local variables
371 integer :: js, ja
372 !! Loop indices.
373 integer :: bin
374 !! Bin for the distribution function.
375 integer :: num_3body_local
376 !! Number of 3-body interactions local to the current atom pair.
377
378
379 1963958 output = 1._real32
380
2/2
✓ Branch 0 taken 4018462 times.
✓ Branch 1 taken 1963958 times.
5982420 num_3body_local = sum(basis%spec(current_idx(1):)%num) - current_idx(2)
381
2/2
✓ Branch 0 taken 539434 times.
✓ Branch 1 taken 1424524 times.
1963958 if(num_3body_local.eq.0) return
382
2/2
✓ Branch 0 taken 3155866 times.
✓ Branch 1 taken 1424524 times.
4580390 species_loop: do js = current_idx(1), basis%nspec, 1
383
2/2
✓ Branch 0 taken 7021774 times.
✓ Branch 1 taken 3155866 times.
11602164 atom_loop: do ja = 1, basis%spec(js)%num
384
2/2
✓ Branch 0 taken 3063243 times.
✓ Branch 1 taken 3958531 times.
7021774 if(js.eq.current_idx(1) .and. ja.le.current_idx(2))cycle
385 3155866 associate( position_store => [ basis%spec(js)%atom(ja,1:3) ] )
386 bin = distribs_container%get_bin( &
387 get_angle( &
388 position_2, &
389 position_1, &
390 position_store &
391 ), dim = 2 &
392 3958531 )
393 output = output * &
394 distribs_container%gdf%df_3body( &
395 bin, &
396 distribs_container%element_map(species) &
397
4/4
✓ Branch 0 taken 11875593 times.
✓ Branch 1 taken 3958531 times.
✓ Branch 2 taken 11875593 times.
✓ Branch 3 taken 3958531 times.
27709717 ) ** ( 1._real32 / real( num_3body_local, real32 ) )
398 end associate
399 end do atom_loop
400 end do species_loop
401
402 1424524 end function evaluate_3body_contributions
403 !###############################################################################
404
405
406 !###############################################################################
407 3958531 pure function evaluate_4body_contributions( distribs_container, &
408 position_1, position_2, position_3, basis, species ) result(output)
409 !! Return the contribution to the viability from 4-body interactions
410 implicit none
411
412 ! Arguments
413 type(distribs_container_type), intent(in) :: distribs_container
414 !! Distribution function (gvector) container.
415 real(real32), dimension(3), intent(in) :: position_1, position_2, position_3
416 !! Positions of the atoms.
417 type(extended_basis_type), intent(in) :: basis
418 !! Basis of the system.
419 integer, intent(in) :: species
420 !! Index of the query element.
421 real(real32) :: output
422 !! Contribution to the viability.
423
424 ! Local variables
425 integer :: ks, ka
426 !! Loop indices.
427 integer :: bin
428 !! Bin for the distribution function.
429 integer :: num_4body_local
430 !! Number of 4-body interactions local to the current atom triplet.
431
432
433 3958531 output = 1._real32
434
2/2
✓ Branch 0 taken 9859585 times.
✓ Branch 1 taken 3958531 times.
13818116 num_4body_local = sum(basis%image_spec(:)%num)
435
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3958531 times.
3958531 if(num_4body_local.eq.0) return
436
2/2
✓ Branch 0 taken 9859585 times.
✓ Branch 1 taken 3958531 times.
13818116 species_loop: do ks = 1, basis%nspec, 1
437
2/2
✓ Branch 0 taken 119170665 times.
✓ Branch 1 taken 9859585 times.
132988781 atom_loop: do ka = 1, basis%image_spec(ks)%num
438 9859585 associate( position_store => [ basis%image_spec(ks)%atom(ka,1:3) ] )
439 bin = distribs_container%get_bin( &
440 get_improper_dihedral_angle( &
441 position_1, &
442 position_2, &
443 position_3, &
444 position_store &
445 ), dim = 3 &
446 119170665 )
447 output = output * &
448 distribs_container%gdf%df_4body( &
449 bin, &
450 distribs_container%element_map(species) &
451
4/4
✓ Branch 0 taken 357511995 times.
✓ Branch 1 taken 119170665 times.
✓ Branch 2 taken 357511995 times.
✓ Branch 3 taken 119170665 times.
834194655 ) ** ( 1._real32 / real( num_4body_local, real32 ) )
452 end associate
453 end do atom_loop
454 end do species_loop
455
456 3958531 end function evaluate_4body_contributions
457 !###############################################################################
458
459 end module raffle__evaluator
460