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

            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              : !> \par History
      10              : !>       created 06-2006 [RD]
      11              : !> \author RD
      12              : ! **************************************************************************************************
      13              : MODULE qs_linres_epr_ownutils
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15              :                                               get_atomic_kind
      16              :    USE cell_types,                      ONLY: cell_type
      17              :    USE cp_control_types,                ONLY: dft_control_type
      18              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      19              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      20              :                                               cp_logger_type
      21              :    USE cp_output_handling,              ONLY: cp_p_file,&
      22              :                                               cp_print_key_finished_output,&
      23              :                                               cp_print_key_should_output,&
      24              :                                               cp_print_key_unit_nr
      25              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      26              :                                               section_vals_type,&
      27              :                                               section_vals_val_get
      28              :    USE kinds,                           ONLY: default_string_length,&
      29              :                                               dp
      30              :    USE mathlib,                         ONLY: diamat_all
      31              :    USE message_passing,                 ONLY: mp_para_env_type
      32              :    USE particle_types,                  ONLY: particle_type
      33              :    USE pw_env_types,                    ONLY: pw_env_get,&
      34              :                                               pw_env_type
      35              :    USE pw_methods,                      ONLY: pw_axpy,&
      36              :                                               pw_integral_ab,&
      37              :                                               pw_scale,&
      38              :                                               pw_transfer,&
      39              :                                               pw_zero
      40              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      41              :                                               pw_pool_type
      42              :    USE pw_spline_utils,                 ONLY: Eval_Interp_Spl3_pbc,&
      43              :                                               find_coeffs,&
      44              :                                               pw_spline_do_precond,&
      45              :                                               pw_spline_precond_create,&
      46              :                                               pw_spline_precond_release,&
      47              :                                               pw_spline_precond_set_kind,&
      48              :                                               pw_spline_precond_type,&
      49              :                                               spl3_pbc
      50              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      51              :                                               pw_r3d_rs_type
      52              :    USE qs_core_energies,                ONLY: calculate_ptrace
      53              :    USE qs_environment_types,            ONLY: get_qs_env,&
      54              :                                               qs_environment_type
      55              :    USE qs_grid_atom,                    ONLY: grid_atom_type
      56              :    USE qs_harmonics_atom,               ONLY: harmonics_atom_type
      57              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      58              :                                               qs_kind_type
      59              :    USE qs_linres_nmr_epr_common_utils,  ONLY: mult_G_ov_G2_grid
      60              :    USE qs_linres_op,                    ONLY: fac_vecp,&
      61              :                                               set_vecp,&
      62              :                                               set_vecp_rev
      63              :    USE qs_linres_types,                 ONLY: current_env_type,&
      64              :                                               epr_env_type,&
      65              :                                               get_current_env,&
      66              :                                               get_epr_env,&
      67              :                                               jrho_atom_type,&
      68              :                                               nablavks_atom_type
      69              :    USE qs_rho_atom_types,               ONLY: get_rho_atom,&
      70              :                                               rho_atom_coeff,&
      71              :                                               rho_atom_type
      72              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      73              :                                               qs_rho_p_type,&
      74              :                                               qs_rho_type
      75              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_type
      76              :    USE util,                            ONLY: get_limit
      77              : #include "./base/base_uses.f90"
      78              : 
      79              :    IMPLICIT NONE
      80              : 
      81              :    PRIVATE
      82              :    PUBLIC :: epr_g_print, epr_g_zke, epr_g_so, epr_g_soo, epr_ind_magnetic_field
      83              : 
      84              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_ownutils'
      85              : 
      86              : CONTAINS
      87              : 
      88              : ! **************************************************************************************************
      89              : !> \brief Prints the g tensor
      90              : !> \param epr_env ...
      91              : !> \param qs_env ...
      92              : !> \par History
      93              : !>      06.2006 created [RD]
      94              : !> \author RD
      95              : ! **************************************************************************************************
      96           14 :    SUBROUTINE epr_g_print(epr_env, qs_env)
      97              : 
      98              :       TYPE(epr_env_type)                                 :: epr_env
      99              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     100              : 
     101              :       CHARACTER(LEN=default_string_length)               :: title
     102              :       INTEGER                                            :: idir1, idir2, output_unit, unit_nr
     103              :       REAL(KIND=dp)                                      :: eigenv_g(3), g_sym(3, 3), gsum
     104              :       TYPE(cp_logger_type), POINTER                      :: logger
     105              :       TYPE(section_vals_type), POINTER                   :: lr_section
     106              : 
     107           14 :       NULLIFY (logger, lr_section)
     108              : 
     109           14 :       logger => cp_get_default_logger()
     110           14 :       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           14 :                                          extension=".linresLog")
     114              : 
     115           14 :       gsum = 0.0_dp
     116              : 
     117           56 :       DO idir1 = 1, 3
     118          182 :          DO idir2 = 1, 3
     119          168 :             gsum = gsum + epr_env%g_total(idir1, idir2)
     120              :          END DO
     121              :       END DO
     122              : 
     123           14 :       IF (output_unit > 0) THEN
     124              :          WRITE (UNIT=output_unit, FMT="(T2,A,T56,E14.6)") &
     125            7 :             "epr|TOT:checksum", gsum
     126              :       END IF
     127              : 
     128              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     129           14 :                                         "PRINT%PROGRAM_RUN_INFO")
     130              : 
     131           14 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, lr_section, &
     132              :                                            "EPR%PRINT%G_TENSOR"), cp_p_file)) THEN
     133              : 
     134              :          unit_nr = cp_print_key_unit_nr(logger, lr_section, "EPR%PRINT%G_TENSOR", &
     135              :                                         extension=".data", middle_name="GTENSOR", &
     136           14 :                                         log_filename=.FALSE.)
     137              : 
     138           14 :          IF (unit_nr > 0) THEN
     139              : 
     140            7 :             WRITE (title, "(A)") "G tensor "
     141            7 :             WRITE (unit_nr, "(T2,A)") title
     142              : 
     143            7 :             WRITE (unit_nr, "(T2,A)") "gmatrix_zke"
     144            7 :             WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_zke, &
     145           14 :                " XY=", 0.0_dp, " XZ=", 0.0_dp
     146            7 :             WRITE (unit_nr, "(3(A,f15.10))") " YX=", 0.0_dp, &
     147           14 :                " YY=", epr_env%g_zke, " YZ=", 0.0_dp
     148            7 :             WRITE (unit_nr, "(3(A,f15.10))") " ZX=", 0.0_dp, &
     149           14 :                " ZY=", 0.0_dp, " ZZ=", epr_env%g_zke
     150              : 
     151            7 :             WRITE (unit_nr, "(T2,A)") "gmatrix_so"
     152            7 :             WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_so(1, 1), &
     153           14 :                " XY=", epr_env%g_so(1, 2), " XZ=", epr_env%g_so(1, 3)
     154            7 :             WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_so(2, 1), &
     155           14 :                " YY=", epr_env%g_so(2, 2), " YZ=", epr_env%g_so(2, 3)
     156            7 :             WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_so(3, 1), &
     157           14 :                " ZY=", epr_env%g_so(3, 2), " ZZ=", epr_env%g_so(3, 3)
     158              : 
     159            7 :             WRITE (unit_nr, "(T2,A)") "gmatrix_soo"
     160            7 :             WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_soo(1, 1), &
     161           14 :                " XY=", epr_env%g_soo(1, 2), " XZ=", epr_env%g_soo(1, 3)
     162            7 :             WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_soo(2, 1), &
     163           14 :                " YY=", epr_env%g_soo(2, 2), " YZ=", epr_env%g_soo(2, 3)
     164            7 :             WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_soo(3, 1), &
     165           14 :                " ZY=", epr_env%g_soo(3, 2), " ZZ=", epr_env%g_soo(3, 3)
     166              : 
     167            7 :             WRITE (unit_nr, "(T2,A)") "gmatrix_total"
     168            7 :             WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_total(1, 1) + epr_env%g_free_factor, &
     169           14 :                " XY=", epr_env%g_total(1, 2), " XZ=", epr_env%g_total(1, 3)
     170            7 :             WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_total(2, 1), &
     171           14 :                " YY=", epr_env%g_total(2, 2) + epr_env%g_free_factor, " YZ=", epr_env%g_total(2, 3)
     172            7 :             WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_total(3, 1), &
     173           14 :                " ZY=", epr_env%g_total(3, 2), " ZZ=", epr_env%g_total(3, 3) + epr_env%g_free_factor
     174              : 
     175           28 :             DO idir1 = 1, 3
     176           91 :                DO idir2 = 1, 3
     177              :                   g_sym(idir1, idir2) = (epr_env%g_total(idir1, idir2) + &
     178           84 :                                          epr_env%g_total(idir2, idir1))/2.0_dp
     179              :                END DO
     180              :             END DO
     181              : 
     182            7 :             WRITE (unit_nr, "(T2,A)") "gtensor_total"
     183            7 :             WRITE (unit_nr, "(3(A,f15.10))") " XX=", g_sym(1, 1) + epr_env%g_free_factor, &
     184           14 :                " XY=", g_sym(1, 2), " XZ=", g_sym(1, 3)
     185            7 :             WRITE (unit_nr, "(3(A,f15.10))") " YX=", g_sym(2, 1), &
     186           14 :                " YY=", g_sym(2, 2) + epr_env%g_free_factor, " YZ=", g_sym(2, 3)
     187            7 :             WRITE (unit_nr, "(3(A,f15.10))") " ZX=", g_sym(3, 1), &
     188           14 :                " ZY=", g_sym(3, 2), " ZZ=", g_sym(3, 3) + epr_env%g_free_factor
     189              : 
     190            7 :             CALL diamat_all(g_sym, eigenv_g)
     191           28 :             eigenv_g(:) = eigenv_g(:)*1.0e6_dp
     192              : 
     193            7 :             WRITE (unit_nr, "(T2,A)") "delta_g principal values in ppm"
     194            7 :             WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(1), " X=", g_sym(1, 1), &
     195           14 :                " Y=", g_sym(2, 1), " Z=", g_sym(3, 1)
     196            7 :             WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(2), " X=", g_sym(1, 2), &
     197           14 :                " Y=", g_sym(2, 2), " Z=", g_sym(3, 2)
     198            7 :             WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(3), " X=", g_sym(1, 3), &
     199           14 :                " Y=", g_sym(2, 3), " Z=", g_sym(3, 3)
     200              : 
     201              :          END IF
     202              : 
     203              :          CALL cp_print_key_finished_output(unit_nr, logger, lr_section,&
     204           14 :               &                            "EPR%PRINT%G_TENSOR")
     205              : 
     206              :       END IF
     207              : 
     208           14 :    END SUBROUTINE epr_g_print
     209              : 
     210              : ! **************************************************************************************************
     211              : !> \brief Calculate zke part of the g tensor
     212              : !> \param epr_env ...
     213              : !> \param qs_env ...
     214              : !> \par History
     215              : !>      06.2006 created [RD]
     216              : !> \author RD
     217              : ! **************************************************************************************************
     218           14 :    SUBROUTINE epr_g_zke(epr_env, qs_env)
     219              : 
     220              :       TYPE(epr_env_type)                                 :: epr_env
     221              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     222              : 
     223              :       INTEGER                                            :: i1, ispin, output_unit
     224              :       REAL(KIND=dp)                                      :: epr_g_zke_temp(2)
     225              :       TYPE(cp_logger_type), POINTER                      :: logger
     226           14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: kinetic, rho_ao
     227              :       TYPE(dft_control_type), POINTER                    :: dft_control
     228              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     229              :       TYPE(qs_rho_type), POINTER                         :: rho
     230              :       TYPE(section_vals_type), POINTER                   :: lr_section
     231              : 
     232           14 :       NULLIFY (dft_control, logger, lr_section, rho, kinetic, para_env, rho_ao)
     233              : 
     234           28 :       logger => cp_get_default_logger()
     235           14 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     236              : 
     237              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     238           14 :                                          extension=".linresLog")
     239              : 
     240              :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     241           14 :                       kinetic=kinetic, rho=rho, para_env=para_env)
     242              : 
     243           14 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     244              : 
     245           42 :       DO ispin = 1, dft_control%nspins
     246              :          CALL calculate_ptrace(kinetic(1)%matrix, rho_ao(ispin)%matrix, &
     247           42 :                                ecore=epr_g_zke_temp(ispin))
     248              :       END DO
     249              : 
     250           14 :       epr_env%g_zke = epr_env%g_zke_factor*(epr_g_zke_temp(1) - epr_g_zke_temp(2))
     251           56 :       DO i1 = 1, 3
     252           56 :          epr_env%g_total(i1, i1) = epr_env%g_total(i1, i1) + epr_env%g_zke
     253              :       END DO
     254              : 
     255           14 :       IF (output_unit > 0) THEN
     256              :          WRITE (UNIT=output_unit, FMT="(T2,A,T56,E24.16)") &
     257            7 :             "epr|ZKE:g_zke", epr_env%g_zke
     258              :       END IF
     259              : 
     260              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     261           14 :                                         "PRINT%PROGRAM_RUN_INFO")
     262              : 
     263           14 :    END SUBROUTINE epr_g_zke
     264              : 
     265              : ! **************************************************************************************************
     266              : !> \brief Calculates g_so
     267              : !> \param epr_env ...
     268              : !> \param current_env ...
     269              : !> \param qs_env ...
     270              : !> \param iB ...
     271              : !> \par History
     272              : !>      06.2006 created [RD]
     273              : !> \author RD
     274              : ! **************************************************************************************************
     275           42 :    SUBROUTINE epr_g_so(epr_env, current_env, qs_env, iB)
     276              : 
     277              :       TYPE(epr_env_type)                                 :: epr_env
     278              :       TYPE(current_env_type)                             :: current_env
     279              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     280              :       INTEGER, INTENT(IN)                                :: iB
     281              : 
     282              :       INTEGER                                            :: aint_precond, ia, iat, iatom, idir1, &
     283              :                                                             idir2, idir3, ikind, ir, ispin, &
     284              :                                                             max_iter, natom, nkind, nspins, &
     285              :                                                             output_unit, precond_kind
     286              :       INTEGER, DIMENSION(2)                              :: bo
     287           42 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
     288              :       LOGICAL                                            :: gapw, paw_atom, success
     289              :       REAL(dp)                                           :: eps_r, eps_x, hard_radius, temp_so_soft, &
     290              :                                                             vks_ra_idir2, vks_ra_idir3
     291              :       REAL(dp), DIMENSION(3, 3)                          :: temp_so_gapw
     292           42 :       REAL(dp), DIMENSION(:, :), POINTER                 :: g_so, g_total
     293              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     294           42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     295              :       TYPE(cp_logger_type), POINTER                      :: logger
     296              :       TYPE(dft_control_type), POINTER                    :: dft_control
     297              :       TYPE(grid_atom_type), POINTER                      :: grid_atom
     298              :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
     299           42 :       TYPE(jrho_atom_type), DIMENSION(:), POINTER        :: jrho1_atom_set
     300              :       TYPE(jrho_atom_type), POINTER                      :: jrho1_atom
     301              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     302           42 :       TYPE(nablavks_atom_type), DIMENSION(:), POINTER    :: nablavks_atom_set
     303              :       TYPE(nablavks_atom_type), POINTER                  :: nablavks_atom
     304           42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     305              :       TYPE(pw_env_type), POINTER                         :: pw_env
     306              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     307           42 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: vks_pw_spline
     308           42 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: jrho2_r, jrho3_r, nrho1_r, nrho2_r, &
     309           42 :                                                             nrho3_r
     310              :       TYPE(pw_spline_precond_type)                       :: precond
     311           42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     312           42 :       TYPE(qs_rho_p_type), DIMENSION(:), POINTER         :: jrho1_set
     313           42 :       TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER      :: nablavks_set
     314              :       TYPE(section_vals_type), POINTER                   :: interp_section, lr_section
     315              : 
     316           42 :       NULLIFY (atomic_kind_set, qs_kind_set, atom_list, dft_control, &
     317           42 :                grid_atom, g_so, g_total, harmonics, interp_section, jrho1_atom_set, &
     318           42 :                jrho1_set, logger, lr_section, nablavks_atom, nablavks_atom_set, &
     319           42 :                nablavks_set, para_env, particle_set, jrho2_r, jrho3_r, nrho1_r, nrho2_r, nrho3_r)
     320              : 
     321           84 :       logger => cp_get_default_logger()
     322           42 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     323              : 
     324              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     325           42 :                                          extension=".linresLog")
     326              : 
     327              :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     328              :                       atomic_kind_set=atomic_kind_set, &
     329              :                       qs_kind_set=qs_kind_set, &
     330              :                       para_env=para_env, pw_env=pw_env, &
     331           42 :                       particle_set=particle_set)
     332              : 
     333              :       CALL get_epr_env(epr_env=epr_env, &
     334              :                        nablavks_set=nablavks_set, &
     335              :                        nablavks_atom_set=nablavks_atom_set, &
     336           42 :                        g_total=g_total, g_so=g_so)
     337              : 
     338              :       CALL get_current_env(current_env=current_env, &
     339           42 :                            jrho1_set=jrho1_set, jrho1_atom_set=jrho1_atom_set)
     340              : 
     341           42 :       gapw = dft_control%qs_control%gapw
     342           42 :       nkind = SIZE(qs_kind_set, 1)
     343           42 :       nspins = dft_control%nspins
     344              : 
     345          168 :       DO idir1 = 1, 3
     346          126 :          CALL set_vecp(idir1, idir2, idir3)
     347              :          ! j_pw x nabla_vks_pw
     348          126 :          temp_so_soft = 0.0_dp
     349          378 :          DO ispin = 1, nspins
     350          252 :             CALL qs_rho_get(jrho1_set(idir2)%rho, rho_r=jrho2_r)
     351          252 :             CALL qs_rho_get(jrho1_set(idir3)%rho, rho_r=jrho3_r)
     352          252 :             CALL qs_rho_get(nablavks_set(idir2, ispin)%rho, rho_r=nrho2_r)
     353          252 :             CALL qs_rho_get(nablavks_set(idir3, ispin)%rho, rho_r=nrho3_r)
     354              :             temp_so_soft = temp_so_soft + (-1.0_dp)**(1 + ispin)*( &
     355              :                            pw_integral_ab(jrho2_r(ispin), nrho3_r(1)) - &
     356          378 :                            pw_integral_ab(jrho3_r(ispin), nrho2_r(1)))
     357              :          END DO
     358          126 :          temp_so_soft = -1.0_dp*epr_env%g_so_factor*temp_so_soft
     359          126 :          IF (output_unit > 0) THEN
     360              :             WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
     361           63 :                "epr|SOX:soft", iB, idir1, temp_so_soft
     362              :          END IF
     363          294 :          g_so(iB, idir1) = temp_so_soft
     364              :       END DO !idir1
     365              : 
     366           42 :       IF (gapw) THEN
     367              : 
     368           30 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     369          330 :          ALLOCATE (vks_pw_spline(3, nspins))
     370              : 
     371              :          interp_section => section_vals_get_subs_vals(lr_section, &
     372           30 :                                                       "EPR%INTERPOLATOR")
     373              :          CALL section_vals_val_get(interp_section, "aint_precond", &
     374           30 :                                    i_val=aint_precond)
     375           30 :          CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
     376           30 :          CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
     377           30 :          CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
     378           30 :          CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
     379              : 
     380           90 :          DO ispin = 1, nspins
     381          270 :             DO idir1 = 1, 3
     382          180 :                CALL auxbas_pw_pool%create_pw(vks_pw_spline(idir1, ispin))
     383              :                ! calculate spline coefficients
     384              :                CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
     385          180 :                                              pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
     386          180 :                CALL qs_rho_get(nablavks_set(idir1, ispin)%rho, rho_r=nrho1_r)
     387              :                CALL pw_spline_do_precond(precond, nrho1_r(1), &
     388          180 :                                          vks_pw_spline(idir1, ispin))
     389          180 :                CALL pw_spline_precond_set_kind(precond, precond_kind)
     390              :                success = find_coeffs(values=nrho1_r(1), &
     391              :                                      coeffs=vks_pw_spline(idir1, ispin), linOp=spl3_pbc, &
     392              :                                      preconditioner=precond, pool=auxbas_pw_pool, &
     393          180 :                                      eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
     394          180 :                CPASSERT(success)
     395          240 :                CALL pw_spline_precond_release(precond)
     396              :             END DO ! idir1
     397              :          END DO ! ispin
     398              : 
     399           30 :          temp_so_gapw = 0.0_dp
     400              : 
     401           90 :          DO ikind = 1, nkind
     402           60 :             NULLIFY (atom_list, grid_atom, harmonics)
     403           60 :             CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
     404              :             CALL get_qs_kind(qs_kind_set(ikind), &
     405              :                              hard_radius=hard_radius, &
     406              :                              grid_atom=grid_atom, &
     407              :                              harmonics=harmonics, &
     408           60 :                              paw_atom=paw_atom)
     409              : 
     410           60 :             IF (.NOT. paw_atom) CYCLE
     411              : 
     412              :             ! Distribute the atoms of this kind
     413              : 
     414           60 :             bo = get_limit(natom, para_env%num_pe, para_env%mepos)
     415              : 
     416          240 :             DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation
     417              :                !                         ! routines (i.e. waiting for parallel sum)
     418           90 :                iatom = atom_list(iat)
     419           90 :                NULLIFY (jrho1_atom, nablavks_atom)
     420           90 :                jrho1_atom => jrho1_atom_set(iatom)
     421           90 :                nablavks_atom => nablavks_atom_set(iatom)
     422          420 :                DO idir1 = 1, 3
     423          270 :                   CALL set_vecp(idir1, idir2, idir3)
     424          900 :                   DO ispin = 1, nspins
     425        27810 :                      DO ir = 1, grid_atom%nr
     426              : 
     427        27000 :                         IF (grid_atom%rad(ir) >= hard_radius) CYCLE
     428              : 
     429       771660 :                         DO ia = 1, grid_atom%ng_sphere
     430              : 
     431      3024000 :                            ra = particle_set(iatom)%r
     432      3024000 :                            ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia)
     433              :                            vks_ra_idir2 = Eval_Interp_Spl3_pbc(ra, &
     434       756000 :                                                                vks_pw_spline(idir2, ispin))
     435              :                            vks_ra_idir3 = Eval_Interp_Spl3_pbc(ra, &
     436       756000 :                                                                vks_pw_spline(idir3, ispin))
     437              : 
     438       756000 :                            IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty:
     439              :                            !                                    !here take care of the partition
     440              : 
     441              :                            ! + sum_A j_loc_h_A x nabla_vks_s_A
     442              :                            temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + &
     443              :                                                      (-1.0_dp)**(1 + ispin)*( &
     444              :                                                      jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
     445              :                                                      vks_ra_idir3 - &
     446              :                                                      jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
     447              :                                                      vks_ra_idir2 &
     448       378000 :                                                      )*grid_atom%wr(ir)*grid_atom%wa(ia)
     449              : 
     450              :                            ! - sum_A j_loc_s_A x nabla_vks_s_A
     451              :                            temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - &
     452              :                                                      (-1.0_dp)**(1 + ispin)*( &
     453              :                                                      jrho1_atom%jrho_vec_rad_s(idir2, ispin)%r_coef(ir, ia)* &
     454              :                                                      vks_ra_idir3 - &
     455              :                                                      jrho1_atom%jrho_vec_rad_s(idir3, ispin)%r_coef(ir, ia)* &
     456              :                                                      vks_ra_idir2 &
     457       378000 :                                                      )*grid_atom%wr(ir)*grid_atom%wa(ia)
     458              : 
     459              :                            ! + sum_A j_loc_h_A x nabla_vks_loc_h_A
     460              :                            temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + &
     461              :                                                      (-1.0_dp)**(1 + ispin)*( &
     462              :                                                      jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
     463              :                                                      nablavks_atom%nablavks_vec_rad_h(idir3, ispin)%r_coef(ir, ia) - &
     464              :                                                      jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
     465              :                                                      nablavks_atom%nablavks_vec_rad_h(idir2, ispin)%r_coef(ir, ia) &
     466       378000 :                                                      )*grid_atom%wr(ir)*grid_atom%wa(ia)
     467              : 
     468              :                            ! - sum_A j_loc_h_A x nabla_vks_loc_s_A
     469              :                            temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - &
     470              :                                                      (-1.0_dp)**(1 + ispin)*( &
     471              :                                                      jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
     472              :                                                      nablavks_atom%nablavks_vec_rad_s(idir3, ispin)%r_coef(ir, ia) - &
     473              :                                                      jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
     474              :                                                      nablavks_atom%nablavks_vec_rad_s(idir2, ispin)%r_coef(ir, ia) &
     475       783000 :                                                      )*grid_atom%wr(ir)*grid_atom%wa(ia)
     476              : 
     477              : ! ORIGINAL
     478              : !                            ! + sum_A j_loc_h_A x nabla_vks_loc_h_A
     479              : !                            temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + &
     480              : !                               (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( &
     481              : !                               jrho1_atom%jrho_vec_rad_h(idir2,iB,ispin)%r_coef(ir,ia) * &
     482              : !                               nablavks_atom%nablavks_vec_rad_h(idir3,ispin)%r_coef(ir,ia) - &
     483              : !                               jrho1_atom%jrho_vec_rad_h(idir3,iB,ispin)%r_coef(ir,ia) * &
     484              : !                               nablavks_atom%nablavks_vec_rad_h(idir2,ispin)%r_coef(ir,ia) &
     485              : !                               ) * grid_atom%wr(ir)*grid_atom%wa(ia)
     486              : !                            ! - sum_A j_loc_s_A x nabla_vks_loc_s_A
     487              : !                            temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - &
     488              : !                               (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( &
     489              : !                               jrho1_atom%jrho_vec_rad_s(idir2,iB,ispin)%r_coef(ir,ia) * &
     490              : !                               nablavks_atom%nablavks_vec_rad_s(idir3,ispin)%r_coef(ir,ia) - &
     491              : !                               jrho1_atom%jrho_vec_rad_s(idir3,iB,ispin)%r_coef(ir,ia) * &
     492              : !                               nablavks_atom%nablavks_vec_rad_s(idir2,ispin)%r_coef(ir,ia) &
     493              : !                               ) * grid_atom%wr(ir)*grid_atom%wa(ia)
     494              :                         END DO !ia
     495              :                      END DO !ir
     496              :                   END DO !ispin
     497              :                END DO !idir1
     498              :             END DO !iat
     499              :          END DO !ikind
     500              : 
     501           30 :          CALL para_env%sum(temp_so_gapw)
     502          390 :          temp_so_gapw(:, :) = -1.0_dp*epr_env%g_so_factor_gapw*temp_so_gapw(:, :)
     503              : 
     504           30 :          IF (output_unit > 0) THEN
     505           60 :             DO idir1 = 1, 3
     506              :                WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
     507           60 :                   "epr|SOX:gapw", iB, idir1, temp_so_gapw(iB, idir1)
     508              :             END DO
     509              :          END IF
     510              : 
     511          120 :          g_so(iB, :) = g_so(iB, :) + temp_so_gapw(iB, :)
     512              : 
     513           90 :          DO ispin = 1, nspins
     514          270 :             DO idir1 = 1, 3
     515          240 :                CALL auxbas_pw_pool%give_back_pw(vks_pw_spline(idir1, ispin))
     516              :             END DO
     517              :          END DO
     518           60 :          DEALLOCATE (vks_pw_spline)
     519              : 
     520              :       END IF ! gapw
     521              : 
     522          294 :       g_total(iB, :) = g_total(iB, :) + g_so(iB, :)
     523              : 
     524              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     525           42 :                                         "PRINT%PROGRAM_RUN_INFO")
     526              : 
     527          336 :    END SUBROUTINE epr_g_so
     528              : 
     529              : ! **************************************************************************************************
     530              : !> \brief Calculates g_soo (soft part only for now)
     531              : !> \param epr_env ...
     532              : !> \param current_env ...
     533              : !> \param qs_env ...
     534              : !> \param iB ...
     535              : !> \par History
     536              : !>      06.2006 created [RD]
     537              : !> \author RD
     538              : ! **************************************************************************************************
     539           42 :    SUBROUTINE epr_g_soo(epr_env, current_env, qs_env, iB)
     540              : 
     541              :       TYPE(epr_env_type)                                 :: epr_env
     542              :       TYPE(current_env_type)                             :: current_env
     543              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     544              :       INTEGER, INTENT(IN)                                :: iB
     545              : 
     546              :       INTEGER                                            :: aint_precond, ia, iat, iatom, idir1, &
     547              :                                                             ikind, ir, iso, ispin, max_iter, &
     548              :                                                             natom, nkind, nspins, output_unit, &
     549              :                                                             precond_kind
     550              :       INTEGER, DIMENSION(2)                              :: bo
     551           42 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
     552              :       LOGICAL                                            :: gapw, paw_atom, soo_rho_hard, success
     553              :       REAL(dp)                                           :: bind_ra_idir1, chi_tensor(3, 3, 2), &
     554              :                                                             eps_r, eps_x, hard_radius, rho_spin, &
     555              :                                                             temp_soo_soft
     556              :       REAL(dp), DIMENSION(3, 3)                          :: temp_soo_gapw
     557           42 :       REAL(dp), DIMENSION(:, :), POINTER                 :: g_soo, g_total
     558              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     559           42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     560              :       TYPE(cp_logger_type), POINTER                      :: logger
     561              :       TYPE(dft_control_type), POINTER                    :: dft_control
     562              :       TYPE(grid_atom_type), POINTER                      :: grid_atom
     563              :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
     564              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     565           42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     566              :       TYPE(pw_env_type), POINTER                         :: pw_env
     567              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     568           42 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: bind_pw_spline
     569           84 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: brho1_r, rho_r
     570              :       TYPE(pw_spline_precond_type)                       :: precond
     571           42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     572           42 :       TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER      :: bind_set
     573              :       TYPE(qs_rho_type), POINTER                         :: rho
     574           42 :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: rho_rad_h, rho_rad_s
     575           42 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
     576              :       TYPE(rho_atom_type), POINTER                       :: rho_atom
     577              :       TYPE(section_vals_type), POINTER                   :: g_section, interp_section, lr_section
     578              : 
     579           42 :       NULLIFY (atomic_kind_set, qs_kind_set, atom_list, bind_set, dft_control, &
     580           42 :                grid_atom, g_section, g_soo, g_total, harmonics, interp_section, &
     581           42 :                logger, lr_section, para_env, particle_set, rho, rho_atom, &
     582           42 :                rho_atom_set, rho_r, brho1_r)
     583              : 
     584           84 :       logger => cp_get_default_logger()
     585           42 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     586              : 
     587              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     588           42 :                                          extension=".linresLog")
     589              : 
     590              :       g_section => section_vals_get_subs_vals(lr_section, &
     591           42 :                                               "EPR%PRINT%G_TENSOR")
     592              : 
     593           42 :       CALL section_vals_val_get(g_section, "soo_rho_hard", l_val=soo_rho_hard)
     594              : 
     595              :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
     596              :                       dft_control=dft_control, para_env=para_env, particle_set=particle_set, &
     597           42 :                       pw_env=pw_env, rho=rho, rho_atom_set=rho_atom_set)
     598              : 
     599              :       CALL get_epr_env(epr_env=epr_env, bind_set=bind_set, &
     600           42 :                        g_soo=g_soo, g_total=g_total)
     601              : 
     602              :       CALL get_current_env(current_env=current_env, &
     603           42 :                            chi_tensor=chi_tensor)
     604           42 :       CALL qs_rho_get(rho, rho_r=rho_r)
     605              : 
     606           42 :       gapw = dft_control%qs_control%gapw
     607           42 :       nkind = SIZE(qs_kind_set, 1)
     608           42 :       nspins = dft_control%nspins
     609              : 
     610          168 :       DO idir1 = 1, 3
     611          126 :          temp_soo_soft = 0.0_dp
     612          378 :          DO ispin = 1, nspins
     613          252 :             CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r)
     614              :             temp_soo_soft = temp_soo_soft + (-1.0_dp)**(1 + ispin)* &
     615          378 :                             pw_integral_ab(brho1_r(1), rho_r(ispin))
     616              :          END DO
     617          126 :          temp_soo_soft = 1.0_dp*epr_env%g_soo_factor*temp_soo_soft
     618          126 :          IF (output_unit > 0) THEN
     619              :             WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") &
     620           63 :                "epr|SOO:soft", iB, idir1, temp_soo_soft
     621              :          END IF
     622          168 :          g_soo(iB, idir1) = temp_soo_soft
     623              :       END DO
     624              : 
     625          168 :       DO idir1 = 1, 3
     626              :          temp_soo_soft = 1.0_dp*epr_env%g_soo_chicorr_factor*chi_tensor(idir1, iB, 2)* &
     627          126 :                          (REAL(dft_control%multiplicity, KIND=dp) - 1.0_dp)
     628          126 :          IF (output_unit > 0) THEN
     629              :             WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") &
     630           63 :                "epr|SOO:soft_g0", iB, idir1, temp_soo_soft
     631              :          END IF
     632          168 :          g_soo(iB, idir1) = g_soo(iB, idir1) + temp_soo_soft
     633              :       END DO
     634              : 
     635           42 :       IF (gapw .AND. soo_rho_hard) THEN
     636              : 
     637           12 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     638          156 :          ALLOCATE (bind_pw_spline(3, 3))
     639              : 
     640              :          interp_section => section_vals_get_subs_vals(lr_section, &
     641           12 :                                                       "EPR%INTERPOLATOR")
     642              :          CALL section_vals_val_get(interp_section, "aint_precond", &
     643           12 :                                    i_val=aint_precond)
     644           12 :          CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
     645           12 :          CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
     646           12 :          CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
     647           12 :          CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
     648              : 
     649           48 :          DO idir1 = 1, 3
     650           36 :             CALL auxbas_pw_pool%create_pw(bind_pw_spline(idir1, iB))
     651              :             ! calculate spline coefficients
     652              :             CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
     653           36 :                                           pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
     654           36 :             CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r)
     655              :             CALL pw_spline_do_precond(precond, brho1_r(1), &
     656           36 :                                       bind_pw_spline(idir1, iB))
     657           36 :             CALL pw_spline_precond_set_kind(precond, precond_kind)
     658              :             success = find_coeffs(values=brho1_r(1), &
     659              :                                   coeffs=bind_pw_spline(idir1, iB), linOp=spl3_pbc, &
     660              :                                   preconditioner=precond, pool=auxbas_pw_pool, &
     661           36 :                                   eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
     662           36 :             CPASSERT(success)
     663           48 :             CALL pw_spline_precond_release(precond)
     664              :          END DO ! idir1
     665              : 
     666           12 :          temp_soo_gapw = 0.0_dp
     667              : 
     668           36 :          DO ikind = 1, nkind
     669           24 :             NULLIFY (atom_list, grid_atom, harmonics)
     670           24 :             CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
     671              :             CALL get_qs_kind(qs_kind_set(ikind), &
     672              :                              hard_radius=hard_radius, &
     673              :                              grid_atom=grid_atom, &
     674              :                              harmonics=harmonics, &
     675           24 :                              paw_atom=paw_atom)
     676              : 
     677           24 :             IF (.NOT. paw_atom) CYCLE
     678              : 
     679              :             ! Distribute the atoms of this kind
     680              : 
     681           24 :             bo = get_limit(natom, para_env%num_pe, para_env%mepos)
     682              : 
     683           96 :             DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation
     684              :                !                         ! routines (i.e. waiting for parallel sum)
     685           36 :                iatom = atom_list(iat)
     686           36 :                rho_atom => rho_atom_set(iatom)
     687           36 :                NULLIFY (rho_rad_h, rho_rad_s)
     688              :                CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, &
     689           36 :                                  rho_rad_s=rho_rad_s)
     690          168 :                DO idir1 = 1, 3
     691          360 :                   DO ispin = 1, nspins
     692        11124 :                      DO ir = 1, grid_atom%nr
     693              : 
     694        10800 :                         IF (grid_atom%rad(ir) >= hard_radius) CYCLE
     695              : 
     696       308664 :                         DO ia = 1, grid_atom%ng_sphere
     697              : 
     698      1209600 :                            ra = particle_set(iatom)%r
     699      1209600 :                            ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia)
     700              :                            bind_ra_idir1 = Eval_Interp_Spl3_pbc(ra, &
     701       302400 :                                                                 bind_pw_spline(idir1, iB))
     702              : 
     703       302400 :                            IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty:
     704              :                            !                                    !here take care of the partition
     705              : 
     706       151200 :                            rho_spin = 0.0_dp
     707              : 
     708      2721600 :                            DO iso = 1, harmonics%max_iso_not0
     709              :                               rho_spin = rho_spin + &
     710              :                                          (rho_rad_h(ispin)%r_coef(ir, iso) - &
     711              :                                           rho_rad_s(ispin)%r_coef(ir, iso))* &
     712      2721600 :                                          harmonics%slm(ia, iso)
     713              :                            END DO
     714              : 
     715              :                            temp_soo_gapw(iB, idir1) = temp_soo_gapw(iB, idir1) + &
     716              :                                                       (-1.0_dp)**(1 + ispin)*( &
     717              :                                                       bind_ra_idir1*rho_spin &
     718       313200 :                                                       )*grid_atom%wr(ir)*grid_atom%wa(ia)
     719              : 
     720              :                         END DO !ia
     721              :                      END DO !ir
     722              :                   END DO ! ispin
     723              :                END DO !idir1
     724              :             END DO !iat
     725              :          END DO !ikind
     726              : 
     727           12 :          CALL para_env%sum(temp_soo_gapw)
     728          156 :          temp_soo_gapw(:, :) = 1.0_dp*epr_env%g_soo_factor*temp_soo_gapw(:, :)
     729              : 
     730           12 :          IF (output_unit > 0) THEN
     731           24 :             DO idir1 = 1, 3
     732              :                WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
     733           24 :                   "epr|SOO:gapw", iB, idir1, temp_soo_gapw(iB, idir1)
     734              :             END DO
     735              :          END IF
     736              : 
     737           48 :          g_soo(iB, :) = g_soo(iB, :) + temp_soo_gapw(iB, :)
     738              : 
     739           48 :          DO idir1 = 1, 3
     740           48 :             CALL auxbas_pw_pool%give_back_pw(bind_pw_spline(idir1, iB))
     741              :          END DO
     742           24 :          DEALLOCATE (bind_pw_spline)
     743              : 
     744              :       END IF ! gapw
     745              : 
     746          294 :       g_total(iB, :) = g_total(iB, :) + g_soo(iB, :)
     747              : 
     748              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     749           42 :                                         "PRINT%PROGRAM_RUN_INFO")
     750              : 
     751          336 :    END SUBROUTINE epr_g_soo
     752              : 
     753              : ! **************************************************************************************************
     754              : !> \brief ...
     755              : !> \param epr_env ...
     756              : !> \param current_env ...
     757              : !> \param qs_env ...
     758              : !> \param iB ...
     759              : ! **************************************************************************************************
     760           42 :    SUBROUTINE epr_ind_magnetic_field(epr_env, current_env, qs_env, iB)
     761              : 
     762              :       TYPE(epr_env_type)                                 :: epr_env
     763              :       TYPE(current_env_type)                             :: current_env
     764              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     765              :       INTEGER, INTENT(IN)                                :: iB
     766              : 
     767              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_ind_magnetic_field'
     768              : 
     769              :       INTEGER                                            :: handle, idir, idir2, idir3, iiB, iiiB, &
     770              :                                                             ispin, natom, nspins
     771              :       LOGICAL                                            :: gapw
     772              :       REAL(dp)                                           :: scale_fac
     773              :       TYPE(cell_type), POINTER                           :: cell
     774              :       TYPE(dft_control_type), POINTER                    :: dft_control
     775           42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     776              :       TYPE(pw_c1d_gs_type)                               :: pw_gspace_work
     777           42 :       TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:, :) :: shift_pw_gspace
     778           42 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: jrho1_g
     779              :       TYPE(pw_env_type), POINTER                         :: pw_env
     780           42 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     781              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     782              :       TYPE(pw_r3d_rs_type)                               :: shift_pw_rspace
     783           42 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: epr_rho_r
     784              :       TYPE(realspace_grid_desc_type), POINTER            :: auxbas_rs_desc
     785              : 
     786           42 :       CALL timeset(routineN, handle)
     787              : 
     788           42 :       NULLIFY (cell, dft_control, pw_env, auxbas_rs_desc, auxbas_pw_pool, &
     789           42 :                pw_pools, particle_set, jrho1_g, epr_rho_r)
     790              : 
     791              :       CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, &
     792           42 :                       particle_set=particle_set)
     793              : 
     794           42 :       gapw = dft_control%qs_control%gapw
     795           42 :       natom = SIZE(particle_set, 1)
     796           42 :       nspins = dft_control%nspins
     797              : 
     798           42 :       CALL get_epr_env(epr_env=epr_env)
     799              : 
     800           42 :       CALL get_current_env(current_env=current_env)
     801              : 
     802           42 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     803              :       CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
     804           42 :                       auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
     805              :       !
     806              :       ! Initialize
     807              :       ! Allocate grids for the calculation of jrho and the shift
     808          462 :       ALLOCATE (shift_pw_gspace(3, nspins))
     809          126 :       DO ispin = 1, nspins
     810          378 :          DO idir = 1, 3
     811          252 :             CALL auxbas_pw_pool%create_pw(shift_pw_gspace(idir, ispin))
     812          336 :             CALL pw_zero(shift_pw_gspace(idir, ispin))
     813              :          END DO
     814              :       END DO
     815           42 :       CALL auxbas_pw_pool%create_pw(shift_pw_rspace)
     816           42 :       CALL pw_zero(shift_pw_rspace)
     817           42 :       CALL auxbas_pw_pool%create_pw(pw_gspace_work)
     818           42 :       CALL pw_zero(pw_gspace_work)
     819              :       !
     820           42 :       CALL set_vecp(iB, iiB, iiiB)
     821              :       !
     822          126 :       DO ispin = 1, nspins
     823              :          !
     824          336 :          DO idir3 = 1, 3
     825              :             ! set to zero for the calculation of the shift
     826          336 :             CALL pw_zero(shift_pw_gspace(idir3, ispin))
     827              :          END DO
     828          378 :          DO idir = 1, 3
     829          252 :             CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g)
     830              :             ! Field gradient
     831              :             ! loop over the Gvec  components: x,y,z
     832         1092 :             DO idir2 = 1, 3
     833         1008 :                IF (idir /= idir2) THEN
     834              :                   ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i))
     835              :                   CALL mult_G_ov_G2_grid(auxbas_pw_pool, jrho1_g(ispin), &
     836          504 :                                          pw_gspace_work, idir2, 0.0_dp)
     837              :                   !
     838              :                   ! scale and add to the correct component of the shift column
     839          504 :                   CALL set_vecp_rev(idir, idir2, idir3)
     840          504 :                   scale_fac = fac_vecp(idir3, idir2, idir)
     841          504 :                   CALL pw_scale(pw_gspace_work, scale_fac)
     842          504 :                   CALL pw_axpy(pw_gspace_work, shift_pw_gspace(idir3, ispin))
     843              :                END IF
     844              :             END DO
     845              :             !
     846              :          END DO ! idir
     847              :       END DO ! ispin
     848              : 
     849              :       ! Store the total soft induced magnetic field (corrected for sic)
     850           42 :       IF (dft_control%nspins == 2) THEN
     851          168 :          DO idir = 1, 3
     852          126 :             CALL qs_rho_get(epr_env%bind_set(idir, iB)%rho, rho_r=epr_rho_r)
     853          168 :             CALL pw_transfer(shift_pw_gspace(idir, 2), epr_rho_r(1))
     854              :          END DO
     855              :       END IF
     856              :       !
     857              :       ! Dellocate grids for the calculation of jrho and the shift
     858           42 :       CALL auxbas_pw_pool%give_back_pw(pw_gspace_work)
     859          126 :       DO ispin = 1, dft_control%nspins
     860          378 :          DO idir = 1, 3
     861          336 :             CALL auxbas_pw_pool%give_back_pw(shift_pw_gspace(idir, ispin))
     862              :          END DO
     863              :       END DO
     864           42 :       DEALLOCATE (shift_pw_gspace)
     865           42 :       CALL auxbas_pw_pool%give_back_pw(shift_pw_rspace)
     866              :       !
     867              :       ! Finalize
     868           42 :       CALL timestop(handle)
     869              :       !
     870           84 :    END SUBROUTINE epr_ind_magnetic_field
     871              : 
     872              : END MODULE qs_linres_epr_ownutils
     873              : 
        

Generated by: LCOV version 2.0-1