LCOV - code coverage report
Current view: top level - src - qs_linres_nmr_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 121 130 93.1 %
Date: 2024-03-29 07:50:05 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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) &
     172           0 :          CPWARN(' To get NMR parameters within PBC you need localized zero order orbitals ')
     173             : 
     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         468 :             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         800 :                 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          30 :                    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 1.15