(git:ed6f26b)
Loading...
Searching...
No Matches
pao_ml.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Main module for PAO Machine Learning
10!> \author Ole Schuett
11! **************************************************************************************************
12MODULE pao_ml
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind
15 USE basis_set_types, ONLY: gto_basis_set_type
16 USE cell_methods, ONLY: cell_create
17 USE cell_types, ONLY: cell_type
18 USE cp_dbcsr_api, ONLY: dbcsr_iterator_blocks_left,&
19 dbcsr_iterator_next_block,&
20 dbcsr_iterator_start,&
21 dbcsr_iterator_stop,&
22 dbcsr_iterator_type,&
23 dbcsr_type
24 USE kinds, ONLY: default_path_length,&
25 default_string_length,&
26 dp
27 USE machine, ONLY: m_flush
28 USE message_passing, ONLY: mp_para_env_type
29 USE pao_input, ONLY: id2str,&
30 pao_ml_gp,&
31 pao_ml_lazy,&
32 pao_ml_nn,&
33 pao_ml_prior_mean,&
34 pao_ml_prior_zero,&
35 pao_rotinv_param
36 USE pao_io, ONLY: pao_ioblock_type,&
37 pao_iokind_type,&
38 pao_kinds_ensure_equal,&
39 pao_read_raw
40 USE pao_ml_descriptor, ONLY: pao_ml_calc_descriptor
41 USE pao_ml_gaussprocess, ONLY: pao_ml_gp_gradient,&
42 pao_ml_gp_predict,&
43 pao_ml_gp_train
44 USE pao_ml_neuralnet, ONLY: pao_ml_nn_gradient,&
45 pao_ml_nn_predict,&
46 pao_ml_nn_train
47 USE pao_types, ONLY: pao_env_type,&
48 training_matrix_type
49 USE particle_types, ONLY: particle_type
50 USE qs_environment_types, ONLY: get_qs_env,&
51 qs_environment_type
52 USE qs_kind_types, ONLY: get_qs_kind,&
53 qs_kind_type
54#include "./base/base_uses.f90"
55
56 IMPLICIT NONE
57
58 PRIVATE
59
60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_ml'
61
62 PUBLIC :: pao_ml_init, pao_ml_predict, pao_ml_forces
63
64 ! linked list used to group training points by kind
65 TYPE training_point_type
66 TYPE(training_point_type), POINTER :: next => null()
67 REAL(dp), DIMENSION(:), ALLOCATABLE :: input
68 REAL(dp), DIMENSION(:), ALLOCATABLE :: output
69 END TYPE training_point_type
70
71 TYPE training_list_type
72 CHARACTER(LEN=default_string_length) :: kindname = ""
73 TYPE(training_point_type), POINTER :: head => null()
74 INTEGER :: npoints = 0
75 END TYPE training_list_type
76
77CONTAINS
78
79! **************************************************************************************************
80!> \brief Initializes the learning machinery
81!> \param pao ...
82!> \param qs_env ...
83! **************************************************************************************************
84 SUBROUTINE pao_ml_init(pao, qs_env)
85 TYPE(pao_env_type), POINTER :: pao
86 TYPE(qs_environment_type), POINTER :: qs_env
87
88 INTEGER :: i
89 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
90 TYPE(mp_para_env_type), POINTER :: para_env
91 TYPE(training_list_type), ALLOCATABLE, &
92 DIMENSION(:) :: training_lists
93
94 IF (SIZE(pao%ml_training_set) == 0) RETURN
95
96 IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO|ML| Initializing maschine learning...'
97
98 IF (pao%parameterization /= pao_rotinv_param) &
99 cpabort("PAO maschine learning requires ROTINV parametrization")
100
101 CALL get_qs_env(qs_env, para_env=para_env, atomic_kind_set=atomic_kind_set)
102
103 ! create training-set data-structure
104 ALLOCATE (training_lists(SIZE(atomic_kind_set)))
105 DO i = 1, SIZE(training_lists)
106 CALL get_atomic_kind(atomic_kind_set(i), name=training_lists(i)%kindname)
107 END DO
108
109 ! parses training files, calculates descriptors and stores all training-points as linked lists
110 DO i = 1, SIZE(pao%ml_training_set)
111 CALL add_to_training_list(pao, qs_env, training_lists, filename=pao%ml_training_set(i)%fn)
112 END DO
113
114 ! ensure there there are training points for all kinds that use pao
115 CALL sanity_check(qs_env, training_lists)
116
117 ! turns linked lists into matrices and syncs them across ranks
118 CALL training_list2matrix(training_lists, pao%ml_training_matrices, para_env)
119
120 ! calculate and subtract prior
121 CALL pao_ml_substract_prior(pao%ml_prior, pao%ml_training_matrices)
122
123 ! print some statistics about the training set and dump it upon request
124 CALL pao_ml_print(pao, pao%ml_training_matrices)
125
126 ! use training-set to train model
127 CALL pao_ml_train(pao)
128
129 END SUBROUTINE pao_ml_init
130
131! **************************************************************************************************
132!> \brief Reads the given file and adds its training points to linked lists.
133!> \param pao ...
134!> \param qs_env ...
135!> \param training_lists ...
136!> \param filename ...
137! **************************************************************************************************
138 SUBROUTINE add_to_training_list(pao, qs_env, training_lists, filename)
139 TYPE(pao_env_type), POINTER :: pao
140 TYPE(qs_environment_type), POINTER :: qs_env
141 TYPE(training_list_type), DIMENSION(:) :: training_lists
142 CHARACTER(LEN=default_path_length) :: filename
143
144 CHARACTER(LEN=default_string_length) :: param
145 INTEGER :: iatom, ikind, natoms, nkinds, nparams
146 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2kind, kindsmap
147 INTEGER, DIMENSION(2) :: ml_range
148 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hmat, positions
149 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
150 TYPE(cell_type), POINTER :: cell
151 TYPE(mp_para_env_type), POINTER :: para_env
152 TYPE(pao_ioblock_type), ALLOCATABLE, DIMENSION(:) :: xblocks
153 TYPE(pao_iokind_type), ALLOCATABLE, DIMENSION(:) :: kinds
154 TYPE(particle_type), DIMENSION(:), POINTER :: my_particle_set
155 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
156 TYPE(training_point_type), POINTER :: new_point
157
158 NULLIFY (new_point, cell)
159
160 IF (pao%iw > 0) WRITE (pao%iw, '(A,A)') " PAO|ML| Reading training frame from file: ", trim(filename)
161
162 CALL get_qs_env(qs_env, para_env=para_env)
163
164 ! parse training data on first rank
165 IF (para_env%is_source()) THEN
166 CALL pao_read_raw(filename, param, hmat, kinds, atom2kind, positions, xblocks, ml_range)
167
168 ! check parametrization
169 IF (trim(param) .NE. trim(adjustl(id2str(pao%parameterization)))) &
170 cpabort("Restart PAO parametrization does not match")
171
172 ! map read-in kinds onto kinds of this run
173 CALL match_kinds(pao, qs_env, kinds, kindsmap)
174 nkinds = SIZE(kindsmap)
175 natoms = SIZE(positions, 1)
176 END IF
177
178 ! broadcast parsed raw training data
179 CALL para_env%bcast(nkinds)
180 CALL para_env%bcast(natoms)
181 IF (.NOT. para_env%is_source()) THEN
182 ALLOCATE (hmat(3, 3))
183 ALLOCATE (kindsmap(nkinds))
184 ALLOCATE (positions(natoms, 3))
185 ALLOCATE (atom2kind(natoms))
186 END IF
187 CALL para_env%bcast(hmat)
188 CALL para_env%bcast(kindsmap)
189 CALL para_env%bcast(atom2kind)
190 CALL para_env%bcast(positions)
191 CALL para_env%bcast(ml_range)
192
193 IF (ml_range(1) /= 1 .OR. ml_range(2) /= natoms) THEN
194 cpwarn("Skipping some atoms for PAO-ML training.")
195 END IF
196
197 ! create cell from read-in h-matrix
198 CALL cell_create(cell, hmat)
199
200 ! create a particle_set based on read-in positions and refere to kinds of this run
201 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
202 ALLOCATE (my_particle_set(natoms))
203 DO iatom = 1, natoms
204 ikind = kindsmap(atom2kind(iatom))
205 my_particle_set(iatom)%atomic_kind => atomic_kind_set(ikind)
206 my_particle_set(iatom)%r = positions(iatom, :)
207 END DO
208
209 ! fill linked list with training points
210 ! Afterwards all ranks will have lists with the same number of entries,
211 ! however the input and output arrays will only be allocated on one rank per entry.
212 ! We farm out the expensive calculation of the descriptor across ranks.
213 DO iatom = 1, natoms
214 IF (iatom < ml_range(1) .OR. ml_range(2) < iatom) cycle
215 ALLOCATE (new_point)
216
217 ! training-point input, calculate descriptor only on one rank
218 IF (mod(iatom - 1, para_env%num_pe) == para_env%mepos) THEN
219 CALL pao_ml_calc_descriptor(pao, &
220 my_particle_set, &
221 qs_kind_set, &
222 cell, &
223 iatom=iatom, &
224 descriptor=new_point%input)
225 END IF
226
227 ! copy training-point output on first rank
228 IF (para_env%is_source()) THEN
229 nparams = SIZE(xblocks(iatom)%p, 1)
230 ALLOCATE (new_point%output(nparams))
231 new_point%output(:) = xblocks(iatom)%p(:, 1)
232 END IF
233
234 ! add to linked list
235 ikind = kindsmap(atom2kind(iatom))
236 training_lists(ikind)%npoints = training_lists(ikind)%npoints + 1
237 new_point%next => training_lists(ikind)%head
238 training_lists(ikind)%head => new_point
239 END DO
240
241 DEALLOCATE (cell, my_particle_set, hmat, kindsmap, positions, atom2kind)
242
243 END SUBROUTINE add_to_training_list
244
245! **************************************************************************************************
246!> \brief Make read-in kinds on to atomic-kinds of this run
247!> \param pao ...
248!> \param qs_env ...
249!> \param kinds ...
250!> \param kindsmap ...
251! **************************************************************************************************
252 SUBROUTINE match_kinds(pao, qs_env, kinds, kindsmap)
253 TYPE(pao_env_type), POINTER :: pao
254 TYPE(qs_environment_type), POINTER :: qs_env
255 TYPE(pao_iokind_type), DIMENSION(:) :: kinds
256 INTEGER, ALLOCATABLE, DIMENSION(:) :: kindsmap
257
258 CHARACTER(LEN=default_string_length) :: name
259 INTEGER :: ikind, jkind
260 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
261
262 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
263
264 cpassert(.NOT. ALLOCATED(kindsmap))
265 ALLOCATE (kindsmap(SIZE(kinds)))
266 kindsmap(:) = -1
267
268 DO ikind = 1, SIZE(kinds)
269 DO jkind = 1, SIZE(atomic_kind_set)
270 CALL get_atomic_kind(atomic_kind_set(jkind), name=name)
271 ! match kinds via their name
272 IF (trim(kinds(ikind)%name) .EQ. trim(name)) THEN
273 CALL pao_kinds_ensure_equal(pao, qs_env, jkind, kinds(ikind))
274 kindsmap(ikind) = jkind
275 EXIT
276 END IF
277 END DO
278 END DO
279
280 IF (any(kindsmap < 1)) &
281 cpabort("PAO: Could not match all kinds from training set")
282 END SUBROUTINE match_kinds
283
284! **************************************************************************************************
285!> \brief Checks that there is at least one training point per pao-enabled kind
286!> \param qs_env ...
287!> \param training_lists ...
288! **************************************************************************************************
289 SUBROUTINE sanity_check(qs_env, training_lists)
290 TYPE(qs_environment_type), POINTER :: qs_env
291 TYPE(training_list_type), DIMENSION(:), TARGET :: training_lists
292
293 INTEGER :: ikind, pao_basis_size
294 TYPE(gto_basis_set_type), POINTER :: basis_set
295 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
296 TYPE(training_list_type), POINTER :: training_list
297
298 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
299
300 DO ikind = 1, SIZE(training_lists)
301 training_list => training_lists(ikind)
302 IF (training_list%npoints > 0) cycle ! it's ok
303 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_basis_size=pao_basis_size)
304 IF (pao_basis_size /= basis_set%nsgf) & ! if this kind has pao enabled...
305 cpabort("Found no training-points for kind: "//trim(training_list%kindname))
306 END DO
307
308 END SUBROUTINE sanity_check
309
310! **************************************************************************************************
311!> \brief Turns the linked lists of training points into matrices
312!> \param training_lists ...
313!> \param training_matrices ...
314!> \param para_env ...
315! **************************************************************************************************
316 SUBROUTINE training_list2matrix(training_lists, training_matrices, para_env)
317 TYPE(training_list_type), ALLOCATABLE, &
318 DIMENSION(:), TARGET :: training_lists
319 TYPE(training_matrix_type), ALLOCATABLE, &
320 DIMENSION(:), TARGET :: training_matrices
321 TYPE(mp_para_env_type), POINTER :: para_env
322
323 INTEGER :: i, ikind, inp_size, ninputs, noutputs, &
324 npoints, out_size
325 TYPE(training_list_type), POINTER :: training_list
326 TYPE(training_matrix_type), POINTER :: training_matrix
327 TYPE(training_point_type), POINTER :: cur_point, prev_point
328
329 cpassert(ALLOCATED(training_lists) .AND. .NOT. ALLOCATED(training_matrices))
330
331 ALLOCATE (training_matrices(SIZE(training_lists)))
332
333 DO ikind = 1, SIZE(training_lists)
334 training_list => training_lists(ikind)
335 training_matrix => training_matrices(ikind)
336 training_matrix%kindname = training_list%kindname ! copy kindname
337 npoints = training_list%npoints ! number of points
338 IF (npoints == 0) THEN
339 ALLOCATE (training_matrix%inputs(0, 0))
340 ALLOCATE (training_matrix%outputs(0, 0))
341 cycle
342 END IF
343
344 ! figure out size of input and output
345 inp_size = 0; out_size = 0
346 IF (ALLOCATED(training_list%head%input)) &
347 inp_size = SIZE(training_list%head%input)
348 IF (ALLOCATED(training_list%head%output)) &
349 out_size = SIZE(training_list%head%output)
350 CALL para_env%sum(inp_size)
351 CALL para_env%sum(out_size)
352
353 ! allocate matices to hold all training points
354 ALLOCATE (training_matrix%inputs(inp_size, npoints))
355 ALLOCATE (training_matrix%outputs(out_size, npoints))
356 training_matrix%inputs(:, :) = 0.0_dp
357 training_matrix%outputs(:, :) = 0.0_dp
358
359 ! loop over all training points, consume linked-list in the process
360 ninputs = 0; noutputs = 0
361 cur_point => training_list%head
362 NULLIFY (training_list%head)
363 DO i = 1, npoints
364 IF (ALLOCATED(cur_point%input)) THEN
365 training_matrix%inputs(:, i) = cur_point%input(:)
366 ninputs = ninputs + 1
367 END IF
368 IF (ALLOCATED(cur_point%output)) THEN
369 training_matrix%outputs(:, i) = cur_point%output(:)
370 noutputs = noutputs + 1
371 END IF
372 ! advance to next entry and deallocate the current one
373 prev_point => cur_point
374 cur_point => cur_point%next
375 DEALLOCATE (prev_point)
376 END DO
377 training_list%npoints = 0 ! list is now empty
378
379 ! sync training_matrix across ranks
380 CALL para_env%sum(training_matrix%inputs)
381 CALL para_env%sum(training_matrix%outputs)
382
383 ! sanity check
384 CALL para_env%sum(noutputs)
385 CALL para_env%sum(ninputs)
386 cpassert(noutputs == npoints .AND. ninputs == npoints)
387 END DO
388
389 END SUBROUTINE training_list2matrix
390
391! **************************************************************************************************
392!> \brief TODO
393!> \param ml_prior ...
394!> \param training_matrices ...
395! **************************************************************************************************
396 SUBROUTINE pao_ml_substract_prior(ml_prior, training_matrices)
397 INTEGER, INTENT(IN) :: ml_prior
398 TYPE(training_matrix_type), DIMENSION(:), TARGET :: training_matrices
399
400 INTEGER :: i, ikind, npoints, out_size
401 TYPE(training_matrix_type), POINTER :: training_matrix
402
403 DO ikind = 1, SIZE(training_matrices)
404 training_matrix => training_matrices(ikind)
405 out_size = SIZE(training_matrix%outputs, 1)
406 npoints = SIZE(training_matrix%outputs, 2)
407 IF (npoints == 0) cycle
408 ALLOCATE (training_matrix%prior(out_size))
409
410 ! calculate prior
411 SELECT CASE (ml_prior)
412 CASE (pao_ml_prior_zero)
413 training_matrix%prior(:) = 0.0_dp
414 CASE (pao_ml_prior_mean)
415 training_matrix%prior(:) = sum(training_matrix%outputs, 2)/real(npoints, dp)
416 CASE DEFAULT
417 cpabort("PAO: unknown prior")
418 END SELECT
419
420 ! subtract prior from all training points
421 DO i = 1, npoints
422 training_matrix%outputs(:, i) = training_matrix%outputs(:, i) - training_matrix%prior
423 END DO
424 END DO
425
426 END SUBROUTINE pao_ml_substract_prior
427
428! **************************************************************************************************
429!> \brief Print some statistics about the training set and dump it upon request
430!> \param pao ...
431!> \param training_matrices ...
432! **************************************************************************************************
433 SUBROUTINE pao_ml_print(pao, training_matrices)
434 TYPE(pao_env_type), POINTER :: pao
435 TYPE(training_matrix_type), DIMENSION(:), TARGET :: training_matrices
436
437 INTEGER :: i, ikind, n, npoints
438 TYPE(training_matrix_type), POINTER :: training_matrix
439
440 ! dump training data
441 IF (pao%iw_mldata > 0) THEN
442 DO ikind = 1, SIZE(training_matrices)
443 training_matrix => training_matrices(ikind)
444 npoints = SIZE(training_matrix%outputs, 2)
445 DO i = 1, npoints
446 WRITE (pao%iw_mldata, *) "PAO|ML| training-point kind: ", trim(training_matrix%kindname), &
447 " point:", i, " in:", training_matrix%inputs(:, i), &
448 " out:", training_matrix%outputs(:, i)
449 END DO
450 END DO
451 CALL m_flush(pao%iw_mldata)
452 END IF
453
454 ! print stats
455 IF (pao%iw > 0) THEN
456 DO ikind = 1, SIZE(training_matrices)
457 training_matrix => training_matrices(ikind)
458 n = SIZE(training_matrix%inputs)
459 IF (n == 0) cycle
460 WRITE (pao%iw, "(A,I3,A,E10.1,1X,E10.1,1X,E10.1)") " PAO|ML| Descriptor for kind: "// &
461 trim(training_matrix%kindname)//" size: ", &
462 SIZE(training_matrix%inputs, 1), " min/mean/max: ", &
463 minval(training_matrix%inputs), &
464 sum(training_matrix%inputs)/real(n, dp), &
465 maxval(training_matrix%inputs)
466 END DO
467 END IF
468
469 END SUBROUTINE pao_ml_print
470
471! **************************************************************************************************
472!> \brief Calls the actual learning algorthim to traing on the given matrices
473!> \param pao ...
474! **************************************************************************************************
475 SUBROUTINE pao_ml_train(pao)
476 TYPE(pao_env_type), POINTER :: pao
477
478 CHARACTER(len=*), PARAMETER :: routinen = 'pao_ml_train'
479
480 INTEGER :: handle
481
482 CALL timeset(routinen, handle)
483
484 SELECT CASE (pao%ml_method)
485 CASE (pao_ml_gp)
486 CALL pao_ml_gp_train(pao)
487 CASE (pao_ml_nn)
488 CALL pao_ml_nn_train(pao)
489 CASE (pao_ml_lazy)
490 ! nothing to do
491 CASE DEFAULT
492 cpabort("PAO: unknown machine learning scheme")
493 END SELECT
494
495 CALL timestop(handle)
496
497 END SUBROUTINE pao_ml_train
498
499! **************************************************************************************************
500!> \brief Fills pao%matrix_X based on machine learning predictions
501!> \param pao ...
502!> \param qs_env ...
503! **************************************************************************************************
504 SUBROUTINE pao_ml_predict(pao, qs_env)
505 TYPE(pao_env_type), POINTER :: pao
506 TYPE(qs_environment_type), POINTER :: qs_env
507
508 CHARACTER(len=*), PARAMETER :: routinen = 'pao_ml_predict'
509
510 INTEGER :: acol, arow, handle, iatom, ikind, natoms
511 REAL(dp), ALLOCATABLE, DIMENSION(:) :: descriptor, variances
512 REAL(dp), DIMENSION(:, :), POINTER :: block_x
513 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
514 TYPE(cell_type), POINTER :: cell
515 TYPE(dbcsr_iterator_type) :: iter
516 TYPE(mp_para_env_type), POINTER :: para_env
517 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
518 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
519
520 CALL timeset(routinen, handle)
521
522 CALL get_qs_env(qs_env, &
523 para_env=para_env, &
524 cell=cell, &
525 particle_set=particle_set, &
526 atomic_kind_set=atomic_kind_set, &
527 qs_kind_set=qs_kind_set, &
528 natom=natoms)
529
530 ! fill matrix_X
531 ALLOCATE (variances(natoms))
532 variances(:) = 0.0_dp
533!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,particle_set,qs_kind_set,cell,variances) &
534!$OMP PRIVATE(iter,arow,acol,iatom,ikind,descriptor,block_X)
535 CALL dbcsr_iterator_start(iter, pao%matrix_X)
536 DO WHILE (dbcsr_iterator_blocks_left(iter))
537 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
538 iatom = arow; cpassert(arow == acol)
539 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
540 IF (SIZE(block_x) == 0) cycle ! pao disabled for iatom
541
542 ! calculate descriptor
543 CALL pao_ml_calc_descriptor(pao, &
544 particle_set, &
545 qs_kind_set, &
546 cell, &
547 iatom, &
548 descriptor)
549
550 ! call actual machine learning for prediction
551 CALL pao_ml_predict_low(pao, ikind=ikind, &
552 descriptor=descriptor, &
553 output=block_x(:, 1), &
554 variance=variances(iatom))
555
556 DEALLOCATE (descriptor)
557
558 !add prior
559 block_x(:, 1) = block_x(:, 1) + pao%ml_training_matrices(ikind)%prior
560 END DO
561 CALL dbcsr_iterator_stop(iter)
562!$OMP END PARALLEL
563
564 ! print variances
565 CALL para_env%sum(variances)
566 IF (pao%iw_mlvar > 0) THEN
567 DO iatom = 1, natoms
568 WRITE (pao%iw_mlvar, *) "PAO|ML| atom:", iatom, " prediction variance:", variances(iatom)
569 END DO
570 CALL m_flush(pao%iw_mlvar)
571 END IF
572
573 ! one-line summary
574 IF (pao%iw > 0) WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO|ML| max prediction variance:", &
575 maxval(variances), " for atom:", maxloc(variances)
576
577 IF (maxval(variances) > pao%ml_tolerance) &
578 cpabort("Variance of prediction above ML_TOLERANCE.")
579
580 DEALLOCATE (variances)
581
582 CALL timestop(handle)
583
584 END SUBROUTINE pao_ml_predict
585
586! **************************************************************************************************
587!> \brief Queries the actual learning algorthim to make a prediction
588!> \param pao ...
589!> \param ikind ...
590!> \param descriptor ...
591!> \param output ...
592!> \param variance ...
593! **************************************************************************************************
594 SUBROUTINE pao_ml_predict_low(pao, ikind, descriptor, output, variance)
595 TYPE(pao_env_type), POINTER :: pao
596 INTEGER, INTENT(IN) :: ikind
597 REAL(dp), DIMENSION(:), INTENT(IN) :: descriptor
598 REAL(dp), DIMENSION(:), INTENT(OUT) :: output
599 REAL(dp), INTENT(OUT) :: variance
600
601 SELECT CASE (pao%ml_method)
602 CASE (pao_ml_gp)
603 CALL pao_ml_gp_predict(pao, ikind, descriptor, output, variance)
604 CASE (pao_ml_nn)
605 CALL pao_ml_nn_predict(pao, ikind, descriptor, output, variance)
606 CASE (pao_ml_lazy)
607 output = 0.0_dp ! let's be really lazy and just rely on the prior
608 variance = 0
609 CASE DEFAULT
610 cpabort("PAO: unknown machine learning scheme")
611 END SELECT
612
613 END SUBROUTINE pao_ml_predict_low
614
615! **************************************************************************************************
616!> \brief Calculate forces contributed by machine learning
617!> \param pao ...
618!> \param qs_env ...
619!> \param matrix_G ...
620!> \param forces ...
621! **************************************************************************************************
622 SUBROUTINE pao_ml_forces(pao, qs_env, matrix_G, forces)
623 TYPE(pao_env_type), POINTER :: pao
624 TYPE(qs_environment_type), POINTER :: qs_env
625 TYPE(dbcsr_type) :: matrix_g
626 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
627
628 CHARACTER(len=*), PARAMETER :: routinen = 'pao_ml_forces'
629
630 INTEGER :: acol, arow, handle, iatom, ikind
631 REAL(dp), ALLOCATABLE, DIMENSION(:) :: descr_grad, descriptor
632 REAL(dp), DIMENSION(:, :), POINTER :: block_g
633 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
634 TYPE(cell_type), POINTER :: cell
635 TYPE(dbcsr_iterator_type) :: iter
636 TYPE(mp_para_env_type), POINTER :: para_env
637 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
638 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
639
640 CALL timeset(routinen, handle)
641
642 CALL get_qs_env(qs_env, &
643 para_env=para_env, &
644 cell=cell, &
645 particle_set=particle_set, &
646 atomic_kind_set=atomic_kind_set, &
647 qs_kind_set=qs_kind_set)
648
649!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_G,particle_set,qs_kind_set,cell) &
650!$OMP REDUCTION(+:forces) &
651!$OMP PRIVATE(iter,arow,acol,iatom,ikind,block_G,descriptor,descr_grad)
652 CALL dbcsr_iterator_start(iter, matrix_g)
653 DO WHILE (dbcsr_iterator_blocks_left(iter))
654 CALL dbcsr_iterator_next_block(iter, arow, acol, block_g)
655 iatom = arow; cpassert(arow == acol)
656 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
657 IF (SIZE(block_g) == 0) cycle ! pao disabled for iatom
658
659 ! calculate descriptor
660 CALL pao_ml_calc_descriptor(pao, &
661 particle_set, &
662 qs_kind_set, &
663 cell, &
664 iatom=iatom, &
665 descriptor=descriptor)
666
667 ! calcaulte derivate of machine learning prediction
668 CALL pao_ml_gradient_low(pao, ikind=ikind, &
669 descriptor=descriptor, &
670 outer_deriv=block_g(:, 1), &
671 gradient=descr_grad)
672
673 ! calculate force contributions from descriptor
674 CALL pao_ml_calc_descriptor(pao, &
675 particle_set, &
676 qs_kind_set, &
677 cell, &
678 iatom=iatom, &
679 descr_grad=descr_grad, &
680 forces=forces)
681
682 DEALLOCATE (descriptor, descr_grad)
683 END DO
684 CALL dbcsr_iterator_stop(iter)
685!$OMP END PARALLEL
686
687 CALL timestop(handle)
688
689 END SUBROUTINE pao_ml_forces
690
691! **************************************************************************************************
692!> \brief Calculate gradient of machine learning algorithm
693!> \param pao ...
694!> \param ikind ...
695!> \param descriptor ...
696!> \param outer_deriv ...
697!> \param gradient ...
698! **************************************************************************************************
699 SUBROUTINE pao_ml_gradient_low(pao, ikind, descriptor, outer_deriv, gradient)
700 TYPE(pao_env_type), POINTER :: pao
701 INTEGER, INTENT(IN) :: ikind
702 REAL(dp), DIMENSION(:), INTENT(IN) :: descriptor, outer_deriv
703 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gradient
704
705 ALLOCATE (gradient(SIZE(descriptor)))
706
707 SELECT CASE (pao%ml_method)
708 CASE (pao_ml_gp)
709 CALL pao_ml_gp_gradient(pao, ikind, descriptor, outer_deriv, gradient)
710 CASE (pao_ml_nn)
711 CALL pao_ml_nn_gradient(pao, ikind, descriptor, outer_deriv, gradient)
712 CASE (pao_ml_lazy)
713 gradient = 0.0_dp
714 CASE DEFAULT
715 cpabort("PAO: unknown machine learning scheme")
716 END SELECT
717
718 END SUBROUTINE pao_ml_gradient_low
719
720END MODULE pao_ml
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:85
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pao_input
Definition pao_input.F:8
pao_input::id2str
character(len=11) function, public id2str(id)
Helper routine.
Definition pao_input.F:216
pao_input::pao_ml_gp
integer, parameter, public pao_ml_gp
Definition pao_input.F:45
pao_input::pao_ml_prior_mean
integer, parameter, public pao_ml_prior_mean
Definition pao_input.F:45
pao_input::pao_ml_lazy
integer, parameter, public pao_ml_lazy
Definition pao_input.F:45
pao_input::pao_rotinv_param
integer, parameter, public pao_rotinv_param
Definition pao_input.F:45
pao_input::pao_ml_prior_zero
integer, parameter, public pao_ml_prior_zero
Definition pao_input.F:45
pao_input::pao_ml_nn
integer, parameter, public pao_ml_nn
Definition pao_input.F:45
pao_io
Routines for reading and writing restart files.
Definition pao_io.F:12
pao_io::pao_kinds_ensure_equal
subroutine, public pao_kinds_ensure_equal(pao, qs_env, ikind, pao_kind)
Ensure that the kind read from the restart is equal to the kind curretly in use.
Definition pao_io.F:326
pao_io::pao_read_raw
subroutine, public pao_read_raw(filename, param, hmat, kinds, atom2kind, positions, xblocks, ml_range)
Reads a restart file into temporary datastructures.
Definition pao_io.F:183
pao_ml_descriptor
Feature vectors for describing chemical environments in a rotationally invariant fashion.
Definition pao_ml_descriptor.F:12
pao_ml_descriptor::pao_ml_calc_descriptor
subroutine, public pao_ml_calc_descriptor(pao, particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
Calculates a descriptor for chemical environment of given atom.
Definition pao_ml_descriptor.F:55
pao_ml_gaussprocess
Gaussian Process implementation.
Definition pao_ml_gaussprocess.F:12
pao_ml_gaussprocess::pao_ml_gp_gradient
subroutine, public pao_ml_gp_gradient(pao, ikind, descriptor, outer_deriv, gradient)
Calculate gradient of Gaussian process.
Definition pao_ml_gaussprocess.F:130
pao_ml_gaussprocess::pao_ml_gp_train
subroutine, public pao_ml_gp_train(pao)
Builds the covariance matrix.
Definition pao_ml_gaussprocess.F:33
pao_ml_gaussprocess::pao_ml_gp_predict
subroutine, public pao_ml_gp_predict(pao, ikind, descriptor, output, variance)
Uses covariance matrix to make prediction.
Definition pao_ml_gaussprocess.F:81
pao_ml_neuralnet
Neural Network implementation.
Definition pao_ml_neuralnet.F:12
pao_ml_neuralnet::pao_ml_nn_gradient
subroutine, public pao_ml_nn_gradient(pao, ikind, descriptor, outer_deriv, gradient)
Calculate gradient of neural network.
Definition pao_ml_neuralnet.F:68
pao_ml_neuralnet::pao_ml_nn_train
subroutine, public pao_ml_nn_train(pao)
Trains the neural network on given training points.
Definition pao_ml_neuralnet.F:121
pao_ml_neuralnet::pao_ml_nn_predict
subroutine, public pao_ml_nn_predict(pao, ikind, descriptor, output, variance)
Uses neural network to make a prediction.
Definition pao_ml_neuralnet.F:44
pao_ml
Main module for PAO Machine Learning.
Definition pao_ml.F:12
pao_ml::pao_ml_forces
subroutine, public pao_ml_forces(pao, qs_env, matrix_g, forces)
Calculate forces contributed by machine learning.
Definition pao_ml.F:623
pao_ml::pao_ml_init
subroutine, public pao_ml_init(pao, qs_env)
Initializes the learning machinery.
Definition pao_ml.F:85
pao_ml::pao_ml_predict
subroutine, public pao_ml_predict(pao, qs_env)
Fills paomatrix_X based on machine learning predictions.
Definition pao_ml.F:505
pao_types
Types used by the PAO machinery.
Definition pao_types.F:12
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pao_io::pao_ioblock_type
Definition pao_io.F:63
pao_io::pao_iokind_type
Definition pao_io.F:67
pao_types::pao_env_type
Definition pao_types.F:147
pao_types::training_matrix_type
Definition pao_types.F:48
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171