231 TYPE(
nnp_type),
INTENT(INOUT),
POINTER :: nnp_env
232 CHARACTER(LEN=*),
INTENT(IN) :: printtag
234 CHARACTER(len=*),
PARAMETER :: routinen =
'nnp_init_model'
235 INTEGER,
PARAMETER :: def_str_len = 256
237 CHARACTER(len=1),
ALLOCATABLE,
DIMENSION(:) :: cactfnct
238 CHARACTER(len=2) :: ele
239 CHARACTER(len=def_str_len) :: dummy, line
240 CHARACTER(len=default_path_length) :: file_name
241 INTEGER :: handle, i, io, j, k, nele, nuc_ele, &
242 symfnct_type, unit_nr
243 LOGICAL :: atom_e_found, explicit, first, found
244 REAL(kind=
dp) :: energy
245 REAL(kind=
dp),
DIMENSION(7) :: test_array
246 REAL(kind=
dp),
DIMENSION(:),
POINTER :: work
251 CALL timeset(routinen, handle)
258 IF (unit_nr > 0)
THEN
259 WRITE (unit_nr, *)
""
260 WRITE (unit_nr, *) trim(printtag)//
"| Neural Network Potential Force Environment"
265 ALLOCATE (nnp_env%atomic_energy(nnp_env%num_atoms, nnp_env%n_committee))
266 ALLOCATE (nnp_env%committee_energy(nnp_env%n_committee))
267 ALLOCATE (nnp_env%myforce(3, nnp_env%num_atoms, nnp_env%n_committee))
268 ALLOCATE (nnp_env%committee_forces(3, nnp_env%num_atoms, nnp_env%n_committee))
269 ALLOCATE (nnp_env%committee_stress(3, 3, nnp_env%n_committee))
273 IF (nnp_env%rad_spline_n < 2)
THEN
274 cpabort(
"NNP| RAD_SPLINE_N must be >= 2.")
277 CALL parser_create(parser, file_name, para_env=logger%para_env)
280 nnp_env%scale_acsf = .false.
281 nnp_env%scale_sigma_acsf = .false.
283 nnp_env%scmin = 0.0_dp
284 nnp_env%scmax = 1.0_dp
285 nnp_env%center_acsf = .false.
286 nnp_env%normnodes = .false.
289 IF (unit_nr > 0)
THEN
290 WRITE (unit_nr, *) trim(printtag)//
"| Reading NNP input from file: ", trim(file_name)
294 search_from_begin_of_file=.true.)
296 READ (line, *) dummy, nnp_env%n_ele
298 CALL cp_abort(__location__, trim(printtag)// &
299 "| number of elements missing in NNP_INPUT_FILE")
303 search_from_begin_of_file=.true.)
304 nnp_env%scale_sigma_acsf = found
307 search_from_begin_of_file=.true.)
308 nnp_env%scale_acsf = found
314 IF (found .AND. nnp_env%scale_sigma_acsf)
THEN
315 cpwarn(
"Ignoring sigma ACSF scaling; both keywords were set.")
316 nnp_env%scale_sigma_acsf = .false.
317 ELSE IF (.NOT. found .AND. nnp_env%scale_sigma_acsf)
THEN
318 nnp_env%scale_acsf = .false.
322 search_from_begin_of_file=.true.)
323 IF (found)
READ (line, *) dummy, nnp_env%scmin
326 search_from_begin_of_file=.true.)
327 IF (found)
READ (line, *) dummy, nnp_env%scmax
330 search_from_begin_of_file=.true.)
331 nnp_env%center_acsf = found
333 IF (nnp_env%scale_sigma_acsf .AND. nnp_env%center_acsf)
THEN
334 nnp_env%scale_sigma_acsf = .false.
337 IF (nnp_env%center_acsf .AND. nnp_env%scale_acsf)
THEN
338 IF ((abs(nnp_env%scmin) > epsilon(0.0_dp)*1.0e+4_dp) .OR. (abs(nnp_env%scmax - 1.0_dp) > epsilon(0.0_dp)*1.0e+4_dp))
THEN
339 CALL cp_warn(__location__, &
340 "Centering and scaling of symmetry functions requested while scale_min_short_atomic != 0 and/or "// &
341 "scale_max_short_atomic != 1. Make sure that scaling and centering of symmetry functions in CP2K "// &
342 "is consistent with your training code. "// &
343 "In CP2K: G* = (G - ave(G)) / (max(G) - min(G)) * (Smax - Smin) + Smin")
348 search_from_begin_of_file=.true.)
349 nnp_env%normnodes = found
352 search_from_begin_of_file=.true.)
354 READ (line, *) dummy, nnp_env%cut_type
356 CALL cp_abort(__location__, trim(printtag)// &
357 "| no cutoff type specified in NNP_INPUT_FILE")
361 search_from_begin_of_file=.true.)
363 READ (line, *) dummy, nnp_env%n_hlayer
365 CALL cp_abort(__location__, trim(printtag)// &
366 "| number of hidden layers missing in NNP_INPUT_FILE")
368 nnp_env%n_layer = nnp_env%n_hlayer + 2
371 ALLOCATE (nnp_env%rad(nele))
372 ALLOCATE (nnp_env%ang(nele))
373 ALLOCATE (nnp_env%n_rad(nele))
374 ALLOCATE (nnp_env%n_ang(nele))
375 ALLOCATE (nnp_env%actfnct(nnp_env%n_hlayer + 1))
376 ALLOCATE (cactfnct(nnp_env%n_hlayer + 1))
377 ALLOCATE (nnp_env%ele(nele))
378 ALLOCATE (nnp_env%nuc_ele(nele))
379 ALLOCATE (nnp_env%arc(nele))
381 ALLOCATE (nnp_env%arc(i)%layer(nnp_env%n_layer))
382 ALLOCATE (nnp_env%arc(i)%n_nodes(nnp_env%n_layer))
384 ALLOCATE (nnp_env%n_hnodes(nnp_env%n_hlayer))
385 ALLOCATE (nnp_env%atom_energies(nele))
386 nnp_env%atom_energies = 0.0_dp
394 IF (trim(adjustl(dummy)) ==
"elements")
THEN
395 READ (line, *) dummy, nnp_env%ele(:)
400 CALL cp_abort(__location__, trim(printtag)// &
401 "| elements not specified in NNP_INPUT_FILE")
406 search_from_begin_of_file=.true.)
408 IF (atom_e_found)
THEN
414 READ (line, *) dummy, ele, energy
416 IF (nnp_env%ele(j) == trim(ele))
THEN
418 nnp_env%atom_energies(j) = energy
423 CALL cp_abort(__location__, trim(printtag)// &
424 "| atom energies are not specified")
430 search_from_begin_of_file=.true.)
432 READ (line, *) dummy, nnp_env%n_hnodes(:)
434 CALL cp_abort(__location__, trim(printtag)// &
435 "NNP| global_nodes_short not specified in NNP_INPUT_FILE")
439 search_from_begin_of_file=.true.)
441 READ (line, *) dummy, cactfnct(:)
443 CALL cp_abort(__location__, trim(printtag)// &
444 "| global_activation_short not specified in NNP_INPUT_FILE")
447 DO i = 1, nnp_env%n_hlayer + 1
448 SELECT CASE (cactfnct(i))
468 CALL cp_abort(__location__, trim(printtag)// &
469 "| Activation function unkown")
485 READ (line, *) dummy, ele, symfnct_type
487 IF (trim(ele) == nnp_env%ele(i))
THEN
488 IF (symfnct_type == 2)
THEN
489 nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
490 ELSE IF (symfnct_type == 3)
THEN
491 nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
493 CALL cp_abort(__location__, trim(printtag)// &
494 "| Symmetry function type not supported")
500 IF (first)
CALL cp_abort(__location__, trim(printtag)// &
501 "| no symfunction_short specified in NNP_INPUT_FILE")
508 ALLOCATE (nnp_env%rad(i)%y(nnp_env%n_rad(i)))
509 ALLOCATE (nnp_env%rad(i)%funccut(nnp_env%n_rad(i)))
510 ALLOCATE (nnp_env%rad(i)%eta(nnp_env%n_rad(i)))
511 ALLOCATE (nnp_env%rad(i)%rs(nnp_env%n_rad(i)))
512 ALLOCATE (nnp_env%rad(i)%loc_min(nnp_env%n_rad(i)))
513 ALLOCATE (nnp_env%rad(i)%loc_max(nnp_env%n_rad(i)))
514 ALLOCATE (nnp_env%rad(i)%loc_av(nnp_env%n_rad(i)))
515 ALLOCATE (nnp_env%rad(i)%sigma(nnp_env%n_rad(i)))
516 ALLOCATE (nnp_env%rad(i)%ele(nnp_env%n_rad(i)))
517 ALLOCATE (nnp_env%rad(i)%nuc_ele(nnp_env%n_rad(i)))
518 nnp_env%rad(i)%funccut = 0.0_dp
519 nnp_env%rad(i)%eta = 0.0_dp
520 nnp_env%rad(i)%rs = 0.0_dp
521 nnp_env%rad(i)%ele =
'X'
522 nnp_env%rad(i)%nuc_ele = 0
524 ALLOCATE (nnp_env%ang(i)%y(nnp_env%n_ang(i)))
525 ALLOCATE (nnp_env%ang(i)%funccut(nnp_env%n_ang(i)))
526 ALLOCATE (nnp_env%ang(i)%eta(nnp_env%n_ang(i)))
527 ALLOCATE (nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
528 ALLOCATE (nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)))
529 ALLOCATE (nnp_env%ang(i)%lam(nnp_env%n_ang(i)))
530 ALLOCATE (nnp_env%ang(i)%loc_min(nnp_env%n_ang(i)))
531 ALLOCATE (nnp_env%ang(i)%loc_max(nnp_env%n_ang(i)))
532 ALLOCATE (nnp_env%ang(i)%loc_av(nnp_env%n_ang(i)))
533 ALLOCATE (nnp_env%ang(i)%sigma(nnp_env%n_ang(i)))
534 ALLOCATE (nnp_env%ang(i)%ele1(nnp_env%n_ang(i)))
535 ALLOCATE (nnp_env%ang(i)%ele2(nnp_env%n_ang(i)))
536 ALLOCATE (nnp_env%ang(i)%nuc_ele1(nnp_env%n_ang(i)))
537 ALLOCATE (nnp_env%ang(i)%nuc_ele2(nnp_env%n_ang(i)))
538 nnp_env%ang(i)%funccut = 0.0_dp
539 nnp_env%ang(i)%eta = 0.0_dp
540 nnp_env%ang(i)%zeta = 0.0_dp
541 nnp_env%ang(i)%prefzeta = 1.0_dp
542 nnp_env%ang(i)%lam = 0.0_dp
543 nnp_env%ang(i)%ele1 =
'X'
544 nnp_env%ang(i)%ele2 =
'X'
545 nnp_env%ang(i)%nuc_ele1 = 0
546 nnp_env%ang(i)%nuc_ele2 = 0
549 nnp_env%arc(i)%n_nodes(1) = nnp_env%n_rad(i) + nnp_env%n_ang(i)
550 nnp_env%arc(i)%n_nodes(2:nnp_env%n_layer - 1) = nnp_env%n_hnodes
551 nnp_env%arc(i)%n_nodes(nnp_env%n_layer) = 1
552 DO j = 1, nnp_env%n_layer
553 ALLOCATE (nnp_env%arc(i)%layer(j)%node(nnp_env%arc(i)%n_nodes(j)))
554 ALLOCATE (nnp_env%arc(i)%layer(j)%node_grad(nnp_env%arc(i)%n_nodes(j)))
555 ALLOCATE (nnp_env%arc(i)%layer(j)%tmp_der(nnp_env%arc(i)%n_nodes(1), nnp_env%arc(i)%n_nodes(j)))
566 nnp_env%max_cut = 0.0_dp
570 READ (line, *) dummy, ele, symfnct_type
572 IF (trim(ele) == nnp_env%ele(i))
THEN
573 IF (symfnct_type == 2)
THEN
574 nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
575 READ (line, *) dummy, ele, symfnct_type, &
576 nnp_env%rad(i)%ele(nnp_env%n_rad(i)), &
577 nnp_env%rad(i)%eta(nnp_env%n_rad(i)), &
578 nnp_env%rad(i)%rs(nnp_env%n_rad(i)), &
579 nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
580 IF (nnp_env%max_cut < nnp_env%rad(i)%funccut(nnp_env%n_rad(i)))
THEN
581 nnp_env%max_cut = nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
583 ELSE IF (symfnct_type == 3)
THEN
584 nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
585 READ (line, *) dummy, ele, symfnct_type, &
586 nnp_env%ang(i)%ele1(nnp_env%n_ang(i)), &
587 nnp_env%ang(i)%ele2(nnp_env%n_ang(i)), &
588 nnp_env%ang(i)%eta(nnp_env%n_ang(i)), &
589 nnp_env%ang(i)%lam(nnp_env%n_ang(i)), &
590 nnp_env%ang(i)%zeta(nnp_env%n_ang(i)), &
591 nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
592 nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)) = &
593 2.0_dp**(1.0_dp - nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
594 IF (nnp_env%max_cut < nnp_env%ang(i)%funccut(nnp_env%n_ang(i)))
THEN
595 nnp_env%max_cut = nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
598 CALL cp_abort(__location__, trim(printtag)// &
599 "| Symmetry function type not supported")
605 IF (first)
CALL cp_abort(__location__, trim(printtag)// &
606 "| no symfunction_short specified in NNP_INPUT_FILE")
613 DO j = 1, nnp_env%n_rad(i)
614 CALL get_ptable_info(nnp_env%rad(i)%ele(j), number=nnp_env%rad(i)%nuc_ele(j))
616 DO j = 1, nnp_env%n_ang(i)
617 CALL get_ptable_info(nnp_env%ang(i)%ele1(j), number=nnp_env%ang(i)%nuc_ele1(j))
618 CALL get_ptable_info(nnp_env%ang(i)%ele2(j), number=nnp_env%ang(i)%nuc_ele2(j))
620 IF (nnp_env%ang(i)%nuc_ele1(j) > nnp_env%ang(i)%nuc_ele2(j))
THEN
621 ele = nnp_env%ang(i)%ele1(j)
622 nnp_env%ang(i)%ele1(j) = nnp_env%ang(i)%ele2(j)
623 nnp_env%ang(i)%ele2(j) = ele
624 nuc_ele = nnp_env%ang(i)%nuc_ele1(j)
625 nnp_env%ang(i)%nuc_ele1(j) = nnp_env%ang(i)%nuc_ele2(j)
626 nnp_env%ang(i)%nuc_ele2(j) = nuc_ele
639 IF (nnp_env%scale_acsf .OR. nnp_env%center_acsf .OR. nnp_env%scale_sigma_acsf)
THEN
640 IF (unit_nr > 0)
THEN
641 WRITE (unit_nr, *) trim(printtag)//
"| Reading scaling information from file: ", trim(file_name)
645 CALL parser_create(parser, file_name, para_env=logger%para_env)
651 READ (parser%input_line, *, iostat=io) test_array(1:k)
657 IF (k == 5 .AND. nnp_env%scale_sigma_acsf)
THEN
658 cpabort(
"Sigma scaling requested, but scaling.data does not contain sigma.")
662 DO i = 1, nnp_env%n_ele
663 DO j = 1, nnp_env%n_rad(i)
665 IF (nnp_env%scale_sigma_acsf)
THEN
666 READ (parser%input_line, *) dummy, dummy, &
667 nnp_env%rad(i)%loc_min(j), &
668 nnp_env%rad(i)%loc_max(j), &
669 nnp_env%rad(i)%loc_av(j), &
670 nnp_env%rad(i)%sigma(j)
672 READ (parser%input_line, *) dummy, dummy, &
673 nnp_env%rad(i)%loc_min(j), &
674 nnp_env%rad(i)%loc_max(j), &
675 nnp_env%rad(i)%loc_av(j)
678 DO j = 1, nnp_env%n_ang(i)
680 IF (nnp_env%scale_sigma_acsf)
THEN
681 READ (parser%input_line, *) dummy, dummy, &
682 nnp_env%ang(i)%loc_min(j), &
683 nnp_env%ang(i)%loc_max(j), &
684 nnp_env%ang(i)%loc_av(j), &
685 nnp_env%ang(i)%sigma(j)
687 READ (parser%input_line, *) dummy, dummy, &
688 nnp_env%ang(i)%loc_min(j), &
689 nnp_env%ang(i)%loc_max(j), &
690 nnp_env%ang(i)%loc_av(j)
699 IF (nnp_env%scale_acsf)
THEN
700 DO i = 1, nnp_env%n_ele
701 DO j = 1, nnp_env%n_rad(i)
702 IF (nnp_env%rad(i)%loc_max(j) <= nnp_env%rad(i)%loc_min(j))
THEN
703 WRITE (line,
'(A,I0,A,I0,A,2(1X,ES13.6))') &
704 "scaling.data: radial sf range non-positive for element ", i, &
705 " sf ", j,
" (loc_min, loc_max) =", &
706 nnp_env%rad(i)%loc_min(j), nnp_env%rad(i)%loc_max(j)
710 DO j = 1, nnp_env%n_ang(i)
711 IF (nnp_env%ang(i)%loc_max(j) <= nnp_env%ang(i)%loc_min(j))
THEN
712 WRITE (line,
'(A,I0,A,I0,A,2(1X,ES13.6))') &
713 "scaling.data: angular sf range non-positive for element ", i, &
714 " sf ", j,
" (loc_min, loc_max) =", &
715 nnp_env%ang(i)%loc_min(j), nnp_env%ang(i)%loc_max(j)
721 IF (nnp_env%scale_sigma_acsf)
THEN
722 DO i = 1, nnp_env%n_ele
723 DO j = 1, nnp_env%n_rad(i)
724 IF (nnp_env%rad(i)%sigma(j) <= 0.0_dp)
THEN
725 WRITE (line,
'(A,I0,A,I0,A,1X,ES13.6)') &
726 "scaling.data: radial sf sigma non-positive for element ", i, &
727 " sf ", j,
" sigma =", nnp_env%rad(i)%sigma(j)
731 DO j = 1, nnp_env%n_ang(i)
732 IF (nnp_env%ang(i)%sigma(j) <= 0.0_dp)
THEN
733 WRITE (line,
'(A,I0,A,I0,A,1X,ES13.6)') &
734 "scaling.data: angular sf sigma non-positive for element ", i, &
735 " sf ", j,
" sigma =", nnp_env%ang(i)%sigma(j)
748 DO i = 1, nnp_env%n_ele
749 DO j = 2, nnp_env%n_layer
750 ALLOCATE (nnp_env%arc(i)%layer(j)%weights(nnp_env%arc(i)%n_nodes(j - 1), &
751 nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
752 ALLOCATE (nnp_env%arc(i)%layer(j)%bweights(nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
755 CALL nnp_read_committee_weights(nnp_env, model_section, printtag)
760 NULLIFY (bias_section)
764 IF (
ASSOCIATED(bias_section))
CALL section_vals_get(bias_section, explicit=explicit)
765 nnp_env%bias = .false.
767 IF (nnp_env%n_committee > 1)
THEN
768 IF (unit_nr > 0)
THEN
769 WRITE (unit_nr, *)
"NNP| Biasing of committee disagreement enabled"
771 nnp_env%bias = .true.
772 ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
773 ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
776 nnp_env%bias_e_avrg(:) = 0.0_dp
778 nnp_env%bias_align = explicit
782 IF (
SIZE(work) /= nnp_env%n_committee)
THEN
783 cpabort(
"ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
785 nnp_env%bias_e_avrg(:) = work
786 IF (unit_nr > 0)
THEN
787 WRITE (unit_nr, *) trim(printtag)//
"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
791 cpwarn(
"NNP committee size is 1, BIAS section is ignored.")
795 IF (unit_nr > 0)
THEN
796 WRITE (unit_nr, *) trim(printtag)//
"| NNP force environment initialized"
799 CALL timestop(handle)