LCOV - code coverage report
Current view: top level - src - qs_linres_nmr_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 93.1 % 130 121
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            Line data    Source code
       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 Chemical shift calculation by dfpt
      10              : !>      Initialization of the nmr_env, creation of the special neighbor lists
      11              : !>      Perturbation Hamiltonians by application of the p and rxp oprtators to  psi0
      12              : !>      Write output
      13              : !>      Deallocate everything
      14              : !> \note
      15              : !>      The psi0 should be localized
      16              : !>      the Sebastiani method works within the assumption that the orbitals are
      17              : !>      completely contained in the simulation box
      18              : !> \par History
      19              : !>       created 07-2005 [MI]
      20              : !> \author MI
      21              : ! **************************************************************************************************
      22              : MODULE qs_linres_nmr_utils
      23              : 
      24              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      25              :    USE cell_types,                      ONLY: cell_type
      26              :    USE cp_control_types,                ONLY: dft_control_type
      27              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      28              :                                               cp_logger_type
      29              :    USE cp_output_handling,              ONLY: cp_p_file,&
      30              :                                               cp_print_key_finished_output,&
      31              :                                               cp_print_key_should_output,&
      32              :                                               cp_print_key_unit_nr
      33              :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      34              :                                               parser_get_object
      35              :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      36              :                                               parser_create,&
      37              :                                               parser_release
      38              :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      39              :                                               cp_unit_to_cp2k
      40              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      41              :                                               section_vals_type,&
      42              :                                               section_vals_val_get
      43              :    USE kinds,                           ONLY: default_path_length,&
      44              :                                               dp
      45              :    USE memory_utilities,                ONLY: reallocate
      46              :    USE particle_types,                  ONLY: particle_type
      47              :    USE pw_env_types,                    ONLY: pw_env_type
      48              :    USE qs_environment_types,            ONLY: get_qs_env,&
      49              :                                               qs_environment_type
      50              :    USE qs_linres_types,                 ONLY: linres_control_type,&
      51              :                                               nmr_env_type
      52              :    USE qs_matrix_pools,                 ONLY: qs_matrix_pools_type
      53              :    USE scf_control_types,               ONLY: scf_control_type
      54              : #include "./base/base_uses.f90"
      55              : 
      56              :    IMPLICIT NONE
      57              : 
      58              :    PRIVATE
      59              :    PUBLIC :: nmr_env_cleanup, nmr_env_init
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_utils'
      62              : 
      63              : CONTAINS
      64              : 
      65              : ! **************************************************************************************************
      66              : !> \brief Initialize the nmr environment
      67              : !> \param nmr_env ...
      68              : !> \param qs_env ...
      69              : !> \par History
      70              : !>      07.2006 created [MI]
      71              : !> \author MI
      72              : ! **************************************************************************************************
      73          160 :    SUBROUTINE nmr_env_init(nmr_env, qs_env)
      74              :       !
      75              :       TYPE(nmr_env_type)                                 :: nmr_env
      76              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      77              : 
      78              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'nmr_env_init'
      79              : 
      80              :       CHARACTER(LEN=default_path_length)                 :: nics_file_name
      81              :       INTEGER                                            :: handle, ini, ir, j, n_mo(2), n_rep, nao, &
      82              :                                                             nat_print, natom, nmoloc, nspins, &
      83              :                                                             output_unit
      84          160 :       INTEGER, DIMENSION(:), POINTER                     :: bounds, list
      85              :       LOGICAL                                            :: gapw
      86          160 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      87              :       TYPE(cell_type), POINTER                           :: cell
      88              :       TYPE(cp_logger_type), POINTER                      :: logger
      89              :       TYPE(dft_control_type), POINTER                    :: dft_control
      90              :       TYPE(linres_control_type), POINTER                 :: linres_control
      91          160 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      92              :       TYPE(pw_env_type), POINTER                         :: pw_env
      93              :       TYPE(qs_matrix_pools_type), POINTER                :: mpools
      94              :       TYPE(scf_control_type), POINTER                    :: scf_control
      95              :       TYPE(section_vals_type), POINTER                   :: lr_section, nmr_section
      96              : 
      97              : !
      98              : 
      99          160 :       CALL timeset(routineN, handle)
     100              : 
     101          160 :       NULLIFY (atomic_kind_set, cell, dft_control, linres_control, scf_control)
     102          160 :       NULLIFY (logger, mpools, nmr_section, particle_set)
     103          160 :       NULLIFY (pw_env)
     104              : 
     105              :       n_mo(1:2) = 0
     106          160 :       nao = 0
     107          160 :       nmoloc = 0
     108              : 
     109          160 :       logger => cp_get_default_logger()
     110          160 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     111              : 
     112              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     113          160 :                                          extension=".linresLog")
     114              : 
     115          160 :       CALL nmr_env_cleanup(nmr_env)
     116              : 
     117          160 :       IF (output_unit > 0) THEN
     118           80 :          WRITE (output_unit, "(/,T20,A,/)") "*** Start NMR Chemical Shift Calculation ***"
     119           80 :          WRITE (output_unit, "(T10,A,/)") "Inizialization of the NMR environment"
     120              :       END IF
     121              : 
     122              :       ! If current_density or full_nmr different allocations are required
     123              :       nmr_section => section_vals_get_subs_vals(qs_env%input, &
     124          160 :            &                                    "PROPERTIES%LINRES%NMR")
     125          160 :       CALL section_vals_val_get(nmr_section, "INTERPOLATE_SHIFT", l_val=nmr_env%interpolate_shift)
     126          160 :       CALL section_vals_val_get(nmr_section, "SHIFT_GAPW_RADIUS", r_val=nmr_env%shift_gapw_radius)
     127          160 :       CALL section_vals_val_get(nmr_section, "NICS", l_val=nmr_env%do_nics)
     128          160 :       IF (nmr_env%do_nics) THEN
     129              :          CALL section_vals_val_get(nmr_section, "NICS_FILE_NAME", &
     130            6 :                                    c_val=nics_file_name)
     131              :          BLOCK
     132              :             CHARACTER(LEN=2)                                   :: label
     133              :             LOGICAL :: my_end
     134              :             TYPE(cp_parser_type) :: parser
     135            6 :             CALL parser_create(parser, nics_file_name)
     136            6 :             CALL parser_get_next_line(parser, 1)
     137            6 :             CALL parser_get_object(parser, nmr_env%n_nics)
     138           18 :             ALLOCATE (nmr_env%r_nics(3, nmr_env%n_nics))
     139            6 :             CALL parser_get_next_line(parser, 2)
     140           24 :             DO j = 1, nmr_env%n_nics
     141           24 :                CALL parser_get_object(parser, label)
     142           24 :                CALL parser_get_object(parser, nmr_env%r_nics(1, j))
     143           24 :                CALL parser_get_object(parser, nmr_env%r_nics(2, j))
     144           24 :                CALL parser_get_object(parser, nmr_env%r_nics(3, j))
     145           24 :                nmr_env%r_nics(1, j) = cp_unit_to_cp2k(nmr_env%r_nics(1, j), "angstrom")
     146           24 :                nmr_env%r_nics(2, j) = cp_unit_to_cp2k(nmr_env%r_nics(2, j), "angstrom")
     147           24 :                nmr_env%r_nics(3, j) = cp_unit_to_cp2k(nmr_env%r_nics(3, j), "angstrom")
     148           24 :                CALL parser_get_next_line(parser, 1, at_end=my_end)
     149           24 :                IF (my_end) EXIT
     150              :             END DO
     151           24 :             CALL parser_release(parser)
     152              :          END BLOCK
     153              :       END IF
     154              : 
     155              :       CALL get_qs_env(qs_env=qs_env, &
     156              :                       atomic_kind_set=atomic_kind_set, &
     157              :                       cell=cell, &
     158              :                       dft_control=dft_control, &
     159              :                       linres_control=linres_control, &
     160              :                       mpools=mpools, &
     161              :                       particle_set=particle_set, &
     162              :                       pw_env=pw_env, &
     163          160 :                       scf_control=scf_control)
     164              :       !
     165              :       ! Check if restat also psi0 should be restarted
     166              :       !IF(nmr_env%restart_nmr .AND. scf_control%density_guess/=restart_guess)THEN
     167              :       !   CPABORT("restart_nmr requires density_guess=restart")
     168              :       !ENDIF
     169              :       !
     170              :       ! check that the psi0 are localized and you have all the centers
     171          160 :       IF (.NOT. linres_control%localized_psi0) THEN
     172            0 :          CPWARN("To get NMR parameters within PBC you need localized zero order orbitals")
     173              :       END IF
     174          160 :       gapw = dft_control%qs_control%gapw
     175          160 :       nspins = dft_control%nspins
     176          160 :       natom = SIZE(particle_set, 1)
     177              : 
     178              :       !
     179              :       ! Conversion factors
     180              :       ! factor for the CHEMICAL SHIFTS: alpha^2 *  ppm.
     181          160 :       nmr_env%shift_factor = (1.0_dp/137.03602_dp)**2*1.0E+6_dp/cell%deth
     182              :       ! factor for the CHEMICAL SHIFTS: alpha^2 *  ppm.
     183          160 :       nmr_env%shift_factor_gapw = (1.0_dp/137.03602_dp)**2*1.0E+6_dp
     184              :       ! chi_factor =  1/4 * e^2/m * a_0 ^2
     185          160 :       nmr_env%chi_factor = 1.9727566E-29_dp/1.0E-30_dp ! -> displayed in 10^-30 J/T^2
     186              :       ! Factor to convert 10^-30 J/T^2 into ppm cgs = ppm cm^3/mol
     187              :       ! = 10^-30 * mu_0/4pi * N_A * 10^6 * 10^6  [one 10^6 for ppm, one for m^3 -> cm^3]
     188          160 :       nmr_env%chi_SI2ppmcgs = 6.022045_dp/1.0E+2_dp
     189              :       ! Chi to Shift: 10^-30  *  2/3  mu_0 / Omega  * 1/ppm
     190              :       nmr_env%chi_SI2shiftppm = 1.0E-30_dp*8.37758041E-7_dp/ &
     191          160 :                                 (cp_unit_from_cp2k(cell%deth, "angstrom^3")*1.0E-30_dp)*1.0E+6_dp
     192              : 
     193          160 :       IF (output_unit > 0) THEN
     194           80 :          WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift gapw radius (a.u.) ", nmr_env%shift_gapw_radius
     195           80 :          IF (nmr_env%do_nics) THEN
     196            3 :             WRITE (output_unit, "(T2,A,T50,I5,A)") "NMR| NICS computed in ", nmr_env%n_nics, " additional points"
     197            3 :             WRITE (output_unit, "(T2,A,T60,A)") "NMR| NICS coordinates read on file ", TRIM(nics_file_name)
     198              :          END IF
     199           80 :          WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift factor (ppm)", nmr_env%shift_factor
     200           80 :          IF (gapw) THEN
     201           43 :             WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift factor gapw (ppm)", nmr_env%shift_factor_gapw
     202              :          END IF
     203           80 :          WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Chi factor (SI)", nmr_env%chi_factor
     204           80 :          WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Conversion Chi (ppm/cgs)", nmr_env%chi_SI2ppmcgs
     205           80 :          WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Conversion Chi to Shift", nmr_env%chi_SI2shiftppm
     206              :       END IF
     207              : 
     208          480 :       ALLOCATE (nmr_env%do_calc_cs_atom(natom))
     209          654 :       nmr_env%do_calc_cs_atom = 0
     210              : 
     211          160 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section,&
     212              :            &    "PRINT%SHIELDING_TENSOR"), cp_p_file)) THEN
     213              : 
     214          160 :          NULLIFY (bounds, list)
     215          160 :          nat_print = 0
     216              :          CALL section_vals_val_get(nmr_section,&
     217              :               &                    "PRINT%SHIELDING_TENSOR%ATOMS_LU_BOUNDS", &
     218          160 :                                    i_vals=bounds)
     219          160 :          nat_print = bounds(2) - bounds(1) + 1
     220          160 :          IF (nat_print > 0) THEN
     221            0 :             ALLOCATE (nmr_env%cs_atom_list(nat_print))
     222            0 :             DO ir = 1, nat_print
     223            0 :                nmr_env%cs_atom_list(ir) = bounds(1) + (ir - 1)
     224            0 :                nmr_env%do_calc_cs_atom(bounds(1) + (ir - 1)) = 1
     225              :             END DO
     226              :          END IF
     227              : 
     228          160 :          IF (.NOT. ASSOCIATED(nmr_env%cs_atom_list)) THEN
     229              :             CALL section_vals_val_get(nmr_section, "PRINT%SHIELDING_TENSOR%ATOMS_LIST", &
     230          160 :                                       n_rep_val=n_rep)
     231          160 :             nat_print = 0
     232          164 :             DO ir = 1, n_rep
     233            4 :                NULLIFY (list)
     234              :                CALL section_vals_val_get(nmr_section, "PRINT%SHIELDING_TENSOR%ATOMS_LIST", &
     235            4 :                                          i_rep_val=ir, i_vals=list)
     236          164 :                IF (ASSOCIATED(list)) THEN
     237            4 :                   CALL reallocate(nmr_env%cs_atom_list, 1, nat_print + SIZE(list))
     238           10 :                   DO ini = 1, SIZE(list)
     239            6 :                      nmr_env%cs_atom_list(ini + nat_print) = list(ini)
     240           10 :                      nmr_env%do_calc_cs_atom(list(ini)) = 1
     241              :                   END DO
     242            4 :                   nat_print = nat_print + SIZE(list)
     243              :                END IF
     244              :             END DO ! ir
     245              :          END IF
     246              : 
     247          160 :          IF (.NOT. ASSOCIATED(nmr_env%cs_atom_list)) THEN
     248          312 :             ALLOCATE (nmr_env%cs_atom_list(natom))
     249          638 :             DO ir = 1, natom
     250          638 :                nmr_env%cs_atom_list(ir) = ir
     251              :             END DO
     252          638 :             nmr_env%do_calc_cs_atom = 1
     253              :          END IF
     254              :          !
     255              :          ! check the list
     256          160 :          CPASSERT(ASSOCIATED(nmr_env%cs_atom_list))
     257          648 :          DO ir = 1, SIZE(nmr_env%cs_atom_list, 1)
     258          488 :             IF (nmr_env%cs_atom_list(ir) .LT. 1 .OR. nmr_env%cs_atom_list(ir) .GT. natom) THEN
     259            0 :                CPABORT("Unknown atom(s)")
     260              :             END IF
     261         2452 :             DO j = 1, SIZE(nmr_env%cs_atom_list, 1)
     262         1804 :                IF (j .EQ. ir) CYCLE
     263         1804 :                IF (nmr_env%cs_atom_list(ir) .EQ. nmr_env%cs_atom_list(j)) THEN
     264            0 :                   CPABORT("Duplicate atoms")
     265              :                END IF
     266              :             END DO
     267              :          END DO
     268              :       ELSE
     269            0 :          NULLIFY (nmr_env%cs_atom_list)
     270              :       END IF
     271              : 
     272          160 :       IF (output_unit > 0) THEN
     273           80 :          IF (ASSOCIATED(nmr_env%cs_atom_list)) THEN
     274           80 :             WRITE (output_unit, "(T2,A,T69,I5,A)") "NMR| Shielding tensor computed for ", &
     275          160 :                SIZE(nmr_env%cs_atom_list, 1), " atoms"
     276              :          ELSE
     277              :             WRITE (output_unit, "(T2,A,T50)")&
     278            0 :                  & "NMR| Shielding tensor not computed at the atomic positions"
     279              :          END IF
     280              :       END IF
     281              :       !
     282              :       ! Initialize the chemical shift tensor
     283              :       ALLOCATE (nmr_env%chemical_shift(3, 3, natom), &
     284          640 :                 nmr_env%chemical_shift_loc(3, 3, natom))
     285         6582 :       nmr_env%chemical_shift = 0.0_dp
     286         6582 :       nmr_env%chemical_shift_loc = 0.0_dp
     287          160 :       IF (nmr_env%do_nics) THEN
     288              :          ALLOCATE (nmr_env%chemical_shift_nics_loc(3, 3, nmr_env%n_nics), &
     289           24 :                    nmr_env%chemical_shift_nics(3, 3, nmr_env%n_nics))
     290          318 :          nmr_env%chemical_shift_nics_loc = 0.0_dp
     291          318 :          nmr_env%chemical_shift_nics = 0.0_dp
     292              :       END IF
     293              : 
     294              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
     295          160 :            &                            "PRINT%PROGRAM_RUN_INFO")
     296              : 
     297          160 :       CALL timestop(handle)
     298              : 
     299          160 :    END SUBROUTINE nmr_env_init
     300              : 
     301              : ! **************************************************************************************************
     302              : !> \brief Deallocate the nmr environment
     303              : !> \param nmr_env ...
     304              : !> \par History
     305              : !>      07.2005 created [MI]
     306              : !> \author MI
     307              : ! **************************************************************************************************
     308          320 :    SUBROUTINE nmr_env_cleanup(nmr_env)
     309              : 
     310              :       TYPE(nmr_env_type)                                 :: nmr_env
     311              : 
     312          320 :       IF (ASSOCIATED(nmr_env%cs_atom_list)) THEN
     313          160 :          DEALLOCATE (nmr_env%cs_atom_list)
     314              :       END IF
     315          320 :       IF (ASSOCIATED(nmr_env%do_calc_cs_atom)) THEN
     316          160 :          DEALLOCATE (nmr_env%do_calc_cs_atom)
     317              :       END IF
     318              :       !chemical_shift
     319          320 :       IF (ASSOCIATED(nmr_env%chemical_shift)) THEN
     320          160 :          DEALLOCATE (nmr_env%chemical_shift)
     321              :       END IF
     322          320 :       IF (ASSOCIATED(nmr_env%chemical_shift_loc)) THEN
     323          160 :          DEALLOCATE (nmr_env%chemical_shift_loc)
     324              :       END IF
     325              :       ! nics
     326          320 :       IF (ASSOCIATED(nmr_env%r_nics)) THEN
     327            6 :          DEALLOCATE (nmr_env%r_nics)
     328              :       END IF
     329          320 :       IF (ASSOCIATED(nmr_env%chemical_shift_nics)) THEN
     330            6 :          DEALLOCATE (nmr_env%chemical_shift_nics)
     331              :       END IF
     332          320 :       IF (ASSOCIATED(nmr_env%chemical_shift_nics_loc)) THEN
     333            6 :          DEALLOCATE (nmr_env%chemical_shift_nics_loc)
     334              :       END IF
     335              : 
     336          320 :    END SUBROUTINE nmr_env_cleanup
     337              : 
     338              : END MODULE qs_linres_nmr_utils
        

Generated by: LCOV version 2.0-1