18 USE dbcsr_api,
ONLY: dbcsr_iterator_blocks_left,&
19 dbcsr_iterator_next_block,&
20 dbcsr_iterator_start,&
54#include "./base/base_uses.f90"
60 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pao_ml'
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
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
91 TYPE(training_list_type),
ALLOCATABLE, &
92 DIMENSION(:) :: training_lists
94 IF (
SIZE(pao%ml_training_set) == 0)
RETURN
96 IF (pao%iw > 0)
WRITE (pao%iw, *)
'PAO|ML| Initializing maschine learning...'
99 cpabort(
"PAO maschine learning requires ROTINV parametrization")
101 CALL get_qs_env(qs_env, para_env=para_env, atomic_kind_set=atomic_kind_set)
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)
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)
115 CALL sanity_check(qs_env, training_lists)
118 CALL training_list2matrix(training_lists, pao%ml_training_matrices, para_env)
121 CALL pao_ml_substract_prior(pao%ml_prior, pao%ml_training_matrices)
124 CALL pao_ml_print(pao, pao%ml_training_matrices)
127 CALL pao_ml_train(pao)
138 SUBROUTINE add_to_training_list(pao, qs_env, training_lists, filename)
141 TYPE(training_list_type),
DIMENSION(:) :: training_lists
142 CHARACTER(LEN=default_path_length) :: filename
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
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
158 NULLIFY (new_point, cell)
160 IF (pao%iw > 0)
WRITE (pao%iw,
'(A,A)')
" PAO|ML| Reading training frame from file: ", trim(filename)
165 IF (para_env%is_source())
THEN
166 CALL pao_read_raw(filename, param, hmat,
kinds, atom2kind, positions, xblocks, ml_range)
169 IF (trim(param) .NE. trim(adjustl(
id2str(pao%parameterization)))) &
170 cpabort(
"Restart PAO parametrization does not match")
173 CALL match_kinds(pao, qs_env,
kinds, kindsmap)
174 nkinds =
SIZE(kindsmap)
175 natoms =
SIZE(positions, 1)
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))
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)
193 IF (ml_range(1) /= 1 .OR. ml_range(2) /= natoms) &
194 cpwarn(
"Skipping some atoms for PAO-ML training.")
200 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
201 ALLOCATE (my_particle_set(natoms))
203 ikind = kindsmap(atom2kind(iatom))
204 my_particle_set(iatom)%atomic_kind => atomic_kind_set(ikind)
205 my_particle_set(iatom)%r = positions(iatom, :)
213 IF (iatom < ml_range(1) .OR. ml_range(2) < iatom) cycle
217 IF (mod(iatom - 1, para_env%num_pe) == para_env%mepos)
THEN
223 descriptor=new_point%input)
227 IF (para_env%is_source())
THEN
228 nparams =
SIZE(xblocks(iatom)%p, 1)
229 ALLOCATE (new_point%output(nparams))
230 new_point%output(:) = xblocks(iatom)%p(:, 1)
234 ikind = kindsmap(atom2kind(iatom))
235 training_lists(ikind)%npoints = training_lists(ikind)%npoints + 1
236 new_point%next => training_lists(ikind)%head
237 training_lists(ikind)%head => new_point
240 DEALLOCATE (cell, my_particle_set, hmat, kindsmap, positions, atom2kind)
242 END SUBROUTINE add_to_training_list
251 SUBROUTINE match_kinds(pao, qs_env, kinds, kindsmap)
255 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: kindsmap
257 CHARACTER(LEN=default_string_length) :: name
258 INTEGER :: ikind, jkind
261 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
263 cpassert(.NOT.
ALLOCATED(kindsmap))
264 ALLOCATE (kindsmap(
SIZE(
kinds)))
267 DO ikind = 1,
SIZE(
kinds)
268 DO jkind = 1,
SIZE(atomic_kind_set)
271 IF (trim(
kinds(ikind)%name) .EQ. trim(name))
THEN
273 kindsmap(ikind) = jkind
279 IF (any(kindsmap < 1)) &
280 cpabort(
"PAO: Could not match all kinds from training set")
281 END SUBROUTINE match_kinds
288 SUBROUTINE sanity_check(qs_env, training_lists)
290 TYPE(training_list_type),
DIMENSION(:),
TARGET :: training_lists
292 INTEGER :: ikind, pao_basis_size
294 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
295 TYPE(training_list_type),
POINTER :: training_list
297 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
299 DO ikind = 1,
SIZE(training_lists)
300 training_list => training_lists(ikind)
301 IF (training_list%npoints > 0) cycle
302 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_basis_size=pao_basis_size)
303 IF (pao_basis_size /= basis_set%nsgf) &
304 cpabort(
"Found no training-points for kind: "//trim(training_list%kindname))
307 END SUBROUTINE sanity_check
315 SUBROUTINE training_list2matrix(training_lists, training_matrices, para_env)
316 TYPE(training_list_type),
ALLOCATABLE, &
317 DIMENSION(:),
TARGET :: training_lists
319 DIMENSION(:),
TARGET :: training_matrices
322 INTEGER :: i, ikind, inp_size, ninputs, noutputs, &
324 TYPE(training_list_type),
POINTER :: training_list
326 TYPE(training_point_type),
POINTER :: cur_point, prev_point
328 cpassert(
ALLOCATED(training_lists) .AND. .NOT.
ALLOCATED(training_matrices))
330 ALLOCATE (training_matrices(
SIZE(training_lists)))
332 DO ikind = 1,
SIZE(training_lists)
333 training_list => training_lists(ikind)
334 training_matrix => training_matrices(ikind)
335 training_matrix%kindname = training_list%kindname
336 npoints = training_list%npoints
337 IF (npoints == 0)
THEN
338 ALLOCATE (training_matrix%inputs(0, 0))
339 ALLOCATE (training_matrix%outputs(0, 0))
344 inp_size = 0; out_size = 0
345 IF (
ALLOCATED(training_list%head%input)) &
346 inp_size =
SIZE(training_list%head%input)
347 IF (
ALLOCATED(training_list%head%output)) &
348 out_size =
SIZE(training_list%head%output)
349 CALL para_env%sum(inp_size)
350 CALL para_env%sum(out_size)
353 ALLOCATE (training_matrix%inputs(inp_size, npoints))
354 ALLOCATE (training_matrix%outputs(out_size, npoints))
355 training_matrix%inputs(:, :) = 0.0_dp
356 training_matrix%outputs(:, :) = 0.0_dp
359 ninputs = 0; noutputs = 0
360 cur_point => training_list%head
361 NULLIFY (training_list%head)
363 IF (
ALLOCATED(cur_point%input))
THEN
364 training_matrix%inputs(:, i) = cur_point%input(:)
365 ninputs = ninputs + 1
367 IF (
ALLOCATED(cur_point%output))
THEN
368 training_matrix%outputs(:, i) = cur_point%output(:)
369 noutputs = noutputs + 1
372 prev_point => cur_point
373 cur_point => cur_point%next
374 DEALLOCATE (prev_point)
376 training_list%npoints = 0
379 CALL para_env%sum(training_matrix%inputs)
380 CALL para_env%sum(training_matrix%outputs)
383 CALL para_env%sum(noutputs)
384 CALL para_env%sum(ninputs)
385 cpassert(noutputs == npoints .AND. ninputs == npoints)
388 END SUBROUTINE training_list2matrix
395 SUBROUTINE pao_ml_substract_prior(ml_prior, training_matrices)
396 INTEGER,
INTENT(IN) :: ml_prior
399 INTEGER :: i, ikind, npoints, out_size
402 DO ikind = 1,
SIZE(training_matrices)
403 training_matrix => training_matrices(ikind)
404 out_size =
SIZE(training_matrix%outputs, 1)
405 npoints =
SIZE(training_matrix%outputs, 2)
406 IF (npoints == 0) cycle
407 ALLOCATE (training_matrix%prior(out_size))
410 SELECT CASE (ml_prior)
412 training_matrix%prior(:) = 0.0_dp
414 training_matrix%prior(:) = sum(training_matrix%outputs, 2)/real(npoints,
dp)
416 cpabort(
"PAO: unknown prior")
421 training_matrix%outputs(:, i) = training_matrix%outputs(:, i) - training_matrix%prior
425 END SUBROUTINE pao_ml_substract_prior
432 SUBROUTINE pao_ml_print(pao, training_matrices)
436 INTEGER :: i, ikind, n, npoints
440 IF (pao%iw_mldata > 0)
THEN
441 DO ikind = 1,
SIZE(training_matrices)
442 training_matrix => training_matrices(ikind)
443 npoints =
SIZE(training_matrix%outputs, 2)
445 WRITE (pao%iw_mldata, *)
"PAO|ML| training-point kind: ", trim(training_matrix%kindname), &
446 " point:", i,
" in:", training_matrix%inputs(:, i), &
447 " out:", training_matrix%outputs(:, i)
455 DO ikind = 1,
SIZE(training_matrices)
456 training_matrix => training_matrices(ikind)
457 n =
SIZE(training_matrix%inputs)
459 WRITE (pao%iw,
"(A,I3,A,E10.1,1X,E10.1,1X,E10.1)")
" PAO|ML| Descriptor for kind: "// &
460 trim(training_matrix%kindname)//
" size: ", &
461 SIZE(training_matrix%inputs, 1),
" min/mean/max: ", &
462 minval(training_matrix%inputs), &
463 sum(training_matrix%inputs)/real(n,
dp), &
464 maxval(training_matrix%inputs)
468 END SUBROUTINE pao_ml_print
474 SUBROUTINE pao_ml_train(pao)
477 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_ml_train'
481 CALL timeset(routinen, handle)
483 SELECT CASE (pao%ml_method)
491 cpabort(
"PAO: unknown machine learning scheme")
494 CALL timestop(handle)
496 END SUBROUTINE pao_ml_train
507 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_ml_predict'
509 INTEGER :: acol, arow, handle, iatom, ikind, natoms
510 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: descriptor, variances
511 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_x
514 TYPE(dbcsr_iterator_type) :: iter
517 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
519 CALL timeset(routinen, handle)
524 particle_set=particle_set, &
525 atomic_kind_set=atomic_kind_set, &
526 qs_kind_set=qs_kind_set, &
530 ALLOCATE (variances(natoms))
531 variances(:) = 0.0_dp
534 CALL dbcsr_iterator_start(iter, pao%matrix_X)
535 DO WHILE (dbcsr_iterator_blocks_left(iter))
536 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
537 iatom = arow; cpassert(arow == acol)
538 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
539 IF (
SIZE(block_x) == 0) cycle
550 CALL pao_ml_predict_low(pao, ikind=ikind, &
551 descriptor=descriptor, &
552 output=block_x(:, 1), &
553 variance=variances(iatom))
555 DEALLOCATE (descriptor)
558 block_x(:, 1) = block_x(:, 1) + pao%ml_training_matrices(ikind)%prior
560 CALL dbcsr_iterator_stop(iter)
564 CALL para_env%sum(variances)
565 IF (pao%iw_mlvar > 0)
THEN
567 WRITE (pao%iw_mlvar, *)
"PAO|ML| atom:", iatom,
" prediction variance:", variances(iatom)
573 IF (pao%iw > 0)
WRITE (pao%iw,
"(A,E20.10,A,T71,I10)")
" PAO|ML| max prediction variance:", &
574 maxval(variances),
" for atom:", maxloc(variances)
576 IF (maxval(variances) > pao%ml_tolerance) &
577 cpabort(
"Variance of prediction above ML_TOLERANCE.")
579 DEALLOCATE (variances)
581 CALL timestop(handle)
593 SUBROUTINE pao_ml_predict_low(pao, ikind, descriptor, output, variance)
595 INTEGER,
INTENT(IN) :: ikind
596 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: descriptor
597 REAL(
dp),
DIMENSION(:),
INTENT(OUT) :: output
598 REAL(
dp),
INTENT(OUT) :: variance
600 SELECT CASE (pao%ml_method)
609 cpabort(
"PAO: unknown machine learning scheme")
612 END SUBROUTINE pao_ml_predict_low
624 TYPE(dbcsr_type) :: matrix_g
625 REAL(
dp),
DIMENSION(:, :),
INTENT(INOUT) :: forces
627 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_ml_forces'
629 INTEGER :: acol, arow, handle, iatom, ikind
630 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: descr_grad, descriptor
631 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_g
634 TYPE(dbcsr_iterator_type) :: iter
637 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
639 CALL timeset(routinen, handle)
644 particle_set=particle_set, &
645 atomic_kind_set=atomic_kind_set, &
646 qs_kind_set=qs_kind_set)
651 CALL dbcsr_iterator_start(iter, matrix_g)
652 DO WHILE (dbcsr_iterator_blocks_left(iter))
653 CALL dbcsr_iterator_next_block(iter, arow, acol, block_g)
654 iatom = arow; cpassert(arow == acol)
655 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
656 IF (
SIZE(block_g) == 0) cycle
664 descriptor=descriptor)
667 CALL pao_ml_gradient_low(pao, ikind=ikind, &
668 descriptor=descriptor, &
669 outer_deriv=block_g(:, 1), &
678 descr_grad=descr_grad, &
681 DEALLOCATE (descriptor, descr_grad)
683 CALL dbcsr_iterator_stop(iter)
686 CALL timestop(handle)
698 SUBROUTINE pao_ml_gradient_low(pao, ikind, descriptor, outer_deriv, gradient)
700 INTEGER,
INTENT(IN) :: ikind
701 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: descriptor, outer_deriv
702 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: gradient
704 ALLOCATE (gradient(
SIZE(descriptor)))
706 SELECT CASE (pao%ml_method)
714 cpabort(
"PAO: unknown machine learning scheme")
717 END SUBROUTINE pao_ml_gradient_low
Define the atomic kind types and their sub types.
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.
Handles all functions related to the CELL.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Interface to the message passing library MPI.
Routines for reading and writing restart files.
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.
subroutine, public pao_read_raw(filename, param, hmat, kinds, atom2kind, positions, xblocks, ml_range)
Reads a restart file into temporary datastructures.
Feature vectors for describing chemical environments in a rotationally invariant fashion.
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.
Gaussian Process implementation.
subroutine, public pao_ml_gp_gradient(pao, ikind, descriptor, outer_deriv, gradient)
Calculate gradient of Gaussian process.
subroutine, public pao_ml_gp_train(pao)
Builds the covariance matrix.
subroutine, public pao_ml_gp_predict(pao, ikind, descriptor, output, variance)
Uses covariance matrix to make prediction.
Neural Network implementation.
subroutine, public pao_ml_nn_gradient(pao, ikind, descriptor, outer_deriv, gradient)
Calculate gradient of neural network.
subroutine, public pao_ml_nn_train(pao)
Trains the neural network on given training points.
subroutine, public pao_ml_nn_predict(pao, ikind, descriptor, output, variance)
Uses neural network to make a prediction.
Main module for PAO Machine Learning.
subroutine, public pao_ml_forces(pao, qs_env, matrix_g, forces)
Calculate forces contributed by machine learning.
subroutine, public pao_ml_init(pao, qs_env)
Initializes the learning machinery.
subroutine, public pao_ml_predict(pao, qs_env)
Fills paomatrix_X based on machine learning predictions.
Types used by the PAO machinery.
Define the data structure for the particle information.
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_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, 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, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.