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

            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 Contains the setup for  the calculation of properties by linear response
      10              : !>      by the application of second order density functional perturbation theory.
      11              : !>      The knowledge of the ground state energy, density and wavefunctions is assumed.
      12              : !>      Uses the self consistent approach.
      13              : !>      Properties that can be calculated : none
      14              : !> \par History
      15              : !>       created 06-2005 [MI]
      16              : !> \author MI
      17              : ! **************************************************************************************************
      18              : MODULE qs_linres_module
      19              :    USE bibliography,                    ONLY: Ditler2021,&
      20              :                                               Ditler2022,&
      21              :                                               Weber2009,&
      22              :                                               cite_reference
      23              :    USE cp_control_types,                ONLY: dft_control_type
      24              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      25              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26              :                                               cp_logger_type
      27              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      28              :                                               cp_print_key_unit_nr
      29              :    USE force_env_types,                 ONLY: force_env_get,&
      30              :                                               force_env_type,&
      31              :                                               use_qmmm,&
      32              :                                               use_qs_force
      33              :    USE input_constants,                 ONLY: lr_current,&
      34              :                                               lr_none,&
      35              :                                               ot_precond_full_all,&
      36              :                                               ot_precond_full_kinetic,&
      37              :                                               ot_precond_full_single,&
      38              :                                               ot_precond_full_single_inverse,&
      39              :                                               ot_precond_none,&
      40              :                                               ot_precond_s_inverse
      41              :    USE input_section_types,             ONLY: section_vals_get,&
      42              :                                               section_vals_get_subs_vals,&
      43              :                                               section_vals_type,&
      44              :                                               section_vals_val_get
      45              :    USE kinds,                           ONLY: dp
      46              :    USE qs_dcdr,                         ONLY: apt_dR,&
      47              :                                               apt_dR_localization,&
      48              :                                               dcdr_build_op_dR,&
      49              :                                               dcdr_response_dR,&
      50              :                                               prepare_per_atom
      51              :    USE qs_dcdr_utils,                   ONLY: dcdr_env_cleanup,&
      52              :                                               dcdr_env_init,&
      53              :                                               dcdr_print
      54              :    USE qs_density_matrices,             ONLY: calculate_density_matrix
      55              :    USE qs_environment_types,            ONLY: get_qs_env,&
      56              :                                               qs_environment_type,&
      57              :                                               set_qs_env
      58              :    USE qs_linres_current,               ONLY: current_build_chi,&
      59              :                                               current_build_current
      60              :    USE qs_linres_current_utils,         ONLY: current_env_cleanup,&
      61              :                                               current_env_init,&
      62              :                                               current_response
      63              :    USE qs_linres_epr_nablavks,          ONLY: epr_nablavks
      64              :    USE qs_linres_epr_ownutils,          ONLY: epr_g_print,&
      65              :                                               epr_g_so,&
      66              :                                               epr_g_soo,&
      67              :                                               epr_g_zke,&
      68              :                                               epr_ind_magnetic_field
      69              :    USE qs_linres_epr_utils,             ONLY: epr_env_cleanup,&
      70              :                                               epr_env_init
      71              :    USE qs_linres_issc_utils,            ONLY: issc_env_cleanup,&
      72              :                                               issc_env_init,&
      73              :                                               issc_issc,&
      74              :                                               issc_print,&
      75              :                                               issc_response
      76              :    USE qs_linres_methods,               ONLY: linres_localize
      77              :    USE qs_linres_nmr_shift,             ONLY: nmr_shift,&
      78              :                                               nmr_shift_print
      79              :    USE qs_linres_nmr_utils,             ONLY: nmr_env_cleanup,&
      80              :                                               nmr_env_init
      81              :    USE qs_linres_op,                    ONLY: current_operators,&
      82              :                                               issc_operators,&
      83              :                                               polar_operators,&
      84              :                                               polar_operators_local,&
      85              :                                               polar_operators_local_wannier
      86              :    USE qs_linres_polar_utils,           ONLY: polar_env_init,&
      87              :                                               polar_polar,&
      88              :                                               polar_print,&
      89              :                                               polar_response
      90              :    USE qs_linres_types,                 ONLY: &
      91              :         current_env_type, dcdr_env_type, epr_env_type, get_polar_env, issc_env_type, &
      92              :         linres_control_type, nmr_env_type, polar_env_type, vcd_env_type
      93              :    USE qs_mfp,                          ONLY: mfp_aat,&
      94              :                                               mfp_build_operator_gauge_dependent,&
      95              :                                               mfp_build_operator_gauge_independent,&
      96              :                                               mfp_response
      97              :    USE qs_mo_types,                     ONLY: mo_set_type
      98              :    USE qs_p_env_methods,                ONLY: p_env_create,&
      99              :                                               p_env_psi0_changed
     100              :    USE qs_p_env_types,                  ONLY: p_env_release,&
     101              :                                               qs_p_env_type
     102              :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho
     103              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     104              :                                               qs_rho_type
     105              :    USE qs_vcd,                          ONLY: aat_dV,&
     106              :                                               apt_dV,&
     107              :                                               prepare_per_atom_vcd,&
     108              :                                               vcd_build_op_dV,&
     109              :                                               vcd_response_dV
     110              :    USE qs_vcd_utils,                    ONLY: vcd_env_cleanup,&
     111              :                                               vcd_env_init,&
     112              :                                               vcd_print
     113              : #include "./base/base_uses.f90"
     114              : 
     115              :    IMPLICIT NONE
     116              : 
     117              :    PRIVATE
     118              :    PUBLIC :: linres_calculation, linres_calculation_low
     119              : 
     120              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_module'
     121              : 
     122              : CONTAINS
     123              : ! *****************************************************************************
     124              : !> \brief Calculates the derivatives of the MO coefficients dC/dV^lambda_beta
     125              : !>         wrt to nuclear velocities. The derivative is indexed by `beta`, the
     126              : !>         electric dipole operator by `alpha`.
     127              : !>        Calculates the APT and AAT in velocity form
     128              : !>               P^lambda_alpha,beta = d< mu_alpha >/dV^lambda_beta
     129              : !>               M^lambda_alpha,beta = d< m_alpha >/dV^lambda_beta
     130              : !> \param qs_env ...
     131              : !> \param p_env ...
     132              : !> \author Edward Ditler
     133              : ! **************************************************************************************************
     134            2 :    SUBROUTINE vcd_linres(qs_env, p_env)
     135              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     136              :       TYPE(qs_p_env_type)                                :: p_env
     137              : 
     138              :       INTEGER                                            :: beta, i, latom
     139              :       LOGICAL                                            :: mfp_is_done, mfp_repeat
     140           60 :       TYPE(vcd_env_type)                                 :: vcd_env
     141              : 
     142            2 :       CALL cite_reference(Ditler2022)
     143              : 
     144              :       ! We need the position perturbation for the velocity perturbation operator
     145            2 :       CALL vcd_env_init(vcd_env, qs_env)
     146              : 
     147            2 :       mfp_repeat = vcd_env%distributed_origin
     148            2 :       mfp_is_done = .FALSE.
     149              : 
     150            2 :       qs_env%linres_control%linres_restart = .TRUE.
     151              : 
     152              :       ! Iterate over the list of atoms for which we want to calculate the APTs/AATs
     153              :       !  default is all atoms.
     154            8 :       DO latom = 1, SIZE(vcd_env%dcdr_env%list_of_atoms)
     155            6 :          vcd_env%dcdr_env%lambda = vcd_env%dcdr_env%list_of_atoms(latom)
     156              : 
     157            6 :          CALL prepare_per_atom(vcd_env%dcdr_env, qs_env)
     158            6 :          CALL prepare_per_atom_vcd(vcd_env, qs_env)
     159              : 
     160           24 :          DO beta = 1, 3                   ! in every direction
     161              : 
     162           18 :             vcd_env%dcdr_env%beta = beta
     163           18 :             vcd_env%dcdr_env%deltaR(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) = 1._dp
     164              : 
     165              :             ! Since we do the heavy lifting anyways, we might also calculate the length form APTs here
     166           18 :             CALL dcdr_build_op_dR(vcd_env%dcdr_env, qs_env)
     167           18 :             CALL dcdr_response_dR(vcd_env%dcdr_env, p_env, qs_env)
     168           18 :             CALL apt_dR(qs_env, vcd_env%dcdr_env)
     169              : 
     170              :             ! And with the position perturbation ready, we can calculate the NVP
     171           18 :             CALL vcd_build_op_dV(vcd_env, qs_env)
     172           18 :             CALL vcd_response_dV(vcd_env, p_env, qs_env)
     173              : 
     174           18 :             CALL apt_dV(vcd_env, qs_env)
     175           18 :             CALL aat_dV(vcd_env, qs_env)
     176              : 
     177           24 :             IF (vcd_env%do_mfp) THEN
     178              :                ! Since we came so far, we might as well calculate the MFP AATs
     179              :                ! If we use a distributed origin we need to compute the MFP response again for each
     180              :                !   atom, because the reference point changes.
     181            0 :                IF (.NOT. mfp_is_done .OR. mfp_repeat) THEN
     182            0 :                   DO i = 1, 3
     183            0 :                      IF (vcd_env%origin_dependent_op_mfp) THEN
     184            0 :                         CPWARN("Using the origin dependent MFP operator")
     185            0 :                         CALL mfp_build_operator_gauge_dependent(vcd_env, qs_env, i)
     186              :                      ELSE
     187            0 :                         CALL mfp_build_operator_gauge_independent(vcd_env, qs_env, i)
     188              :                      END IF
     189            0 :                      CALL mfp_response(vcd_env, p_env, qs_env, i)
     190              :                   END DO
     191              :                   mfp_is_done = .TRUE.
     192              :                END IF
     193              : 
     194            0 :                CALL mfp_aat(vcd_env, qs_env)
     195              :             END IF
     196              :          END DO ! beta
     197              : 
     198              :          vcd_env%dcdr_env%apt_total_dcdr(:, :, vcd_env%dcdr_env%lambda) = &
     199              :             vcd_env%dcdr_env%apt_el_dcdr(:, :, vcd_env%dcdr_env%lambda) &
     200           78 :             + vcd_env%dcdr_env%apt_nuc_dcdr(:, :, vcd_env%dcdr_env%lambda)
     201              : 
     202              :          vcd_env%apt_total_nvpt(:, :, vcd_env%dcdr_env%lambda) = &
     203           78 :             vcd_env%apt_el_nvpt(:, :, vcd_env%dcdr_env%lambda) + vcd_env%apt_nuc_nvpt(:, :, vcd_env%dcdr_env%lambda)
     204              : 
     205            6 :          IF (vcd_env%do_mfp) &
     206            2 :             vcd_env%aat_atom_mfp(:, :, vcd_env%dcdr_env%lambda) = vcd_env%aat_atom_mfp(:, :, vcd_env%dcdr_env%lambda)*4._dp
     207              : 
     208              :       END DO !lambda
     209              : 
     210            2 :       CALL vcd_print(vcd_env, qs_env)
     211            2 :       CALL vcd_env_cleanup(qs_env, vcd_env)
     212              : 
     213            2 :    END SUBROUTINE vcd_linres
     214              : 
     215              : ! **************************************************************************************************
     216              : !> \brief Calculates the derivatives of the MO coefficients dC/dR^lambda_beta
     217              : !>         wrt to nuclear coordinates. The derivative is index by `beta`, the
     218              : !>         electric dipole operator by `alpha`.
     219              : !>        Also calculates the APT
     220              : !>               P^lambda_alpha,beta = d< mu_alpha >/dR^lambda_beta
     221              : !>        and calculates the sum rules for the APT elements.
     222              : !> \param qs_env ...
     223              : !> \param p_env ...
     224              : ! **************************************************************************************************
     225           22 :    SUBROUTINE dcdr_linres(qs_env, p_env)
     226              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     227              :       TYPE(qs_p_env_type)                                :: p_env
     228              : 
     229              :       INTEGER                                            :: beta, latom
     230          308 :       TYPE(dcdr_env_type)                                :: dcdr_env
     231              :       TYPE(polar_env_type), POINTER                      :: polar_env
     232              : 
     233           22 :       CALL cite_reference(Ditler2021)
     234           22 :       CALL dcdr_env_init(dcdr_env, qs_env)
     235              : 
     236           22 :       IF (.NOT. dcdr_env%z_matrix_method) THEN
     237              : 
     238           72 :          DO latom = 1, SIZE(dcdr_env%list_of_atoms)
     239           54 :             dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
     240           54 :             CALL prepare_per_atom(dcdr_env, qs_env)
     241              : 
     242          216 :             DO beta = 1, 3                   ! in every direction
     243          162 :                dcdr_env%beta = beta
     244          162 :                dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
     245              : 
     246          162 :                CALL dcdr_build_op_dR(dcdr_env, qs_env)
     247          162 :                CALL dcdr_response_dR(dcdr_env, p_env, qs_env)
     248              : 
     249          216 :                IF (.NOT. dcdr_env%localized_psi0) THEN
     250          126 :                   CALL apt_dR(qs_env, dcdr_env)
     251              :                ELSE IF (dcdr_env%localized_psi0) THEN
     252           36 :                   CALL apt_dR_localization(qs_env, dcdr_env)
     253              :                END IF
     254              : 
     255              :             END DO !beta
     256              : 
     257              :             dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
     258          720 :                dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
     259              :          END DO !lambda
     260              : 
     261              :       ELSE
     262              : 
     263            4 :          CALL polar_env_init(qs_env)
     264            4 :          CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
     265            4 :          CALL get_polar_env(polar_env=polar_env)
     266              : 
     267            4 :          IF (.NOT. dcdr_env%localized_psi0) THEN
     268            4 :             CALL polar_operators_local(qs_env)
     269              :          ELSE
     270            0 :             CALL polar_operators_local_wannier(qs_env, dcdr_env)
     271              :          END IF
     272              : 
     273            4 :          polar_env%do_periodic = .FALSE.
     274            4 :          CALL polar_response(p_env, qs_env)
     275              : 
     276           16 :          DO latom = 1, SIZE(dcdr_env%list_of_atoms)
     277           12 :             dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
     278           12 :             CALL prepare_per_atom(dcdr_env, qs_env)
     279              : 
     280           48 :             DO beta = 1, 3                   ! in every direction
     281           36 :                dcdr_env%beta = beta
     282           36 :                dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
     283              : 
     284           36 :                CALL dcdr_build_op_dR(dcdr_env, qs_env)
     285           48 :                IF (.NOT. dcdr_env%localized_psi0) THEN
     286           36 :                   CALL apt_dR(qs_env, dcdr_env)
     287              :                ELSE
     288            0 :                   CALL apt_dR_localization(qs_env, dcdr_env)
     289              :                END IF
     290              :             END DO !beta
     291              : 
     292              :             dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
     293          160 :                dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
     294              :          END DO !lambda
     295              : 
     296              :       END IF
     297              : 
     298           22 :       CALL dcdr_print(dcdr_env, qs_env)
     299           22 :       CALL dcdr_env_cleanup(qs_env, dcdr_env)
     300           22 :    END SUBROUTINE dcdr_linres
     301              : 
     302              : ! **************************************************************************************************
     303              : !> \brief Driver for the linear response calculatios
     304              : !> \param force_env ...
     305              : !> \par History
     306              : !>      06.2005 created [MI]
     307              : !> \author MI
     308              : ! **************************************************************************************************
     309          188 :    SUBROUTINE linres_calculation(force_env)
     310              : 
     311              :       TYPE(force_env_type), POINTER                      :: force_env
     312              : 
     313              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation'
     314              : 
     315              :       INTEGER                                            :: handle
     316              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     317              : 
     318          188 :       CALL timeset(routineN, handle)
     319              : 
     320          188 :       NULLIFY (qs_env)
     321              : 
     322          188 :       CPASSERT(ASSOCIATED(force_env))
     323          188 :       CPASSERT(force_env%ref_count > 0)
     324              : 
     325          370 :       SELECT CASE (force_env%in_use)
     326              :       CASE (use_qs_force)
     327          182 :          CALL force_env_get(force_env, qs_env=qs_env)
     328              :       CASE (use_qmmm)
     329            6 :          qs_env => force_env%qmmm_env%qs_env
     330              :       CASE DEFAULT
     331          188 :          CPABORT("Does not recognize this force_env")
     332              :       END SELECT
     333              : 
     334          188 :       qs_env%linres_run = .TRUE.
     335              : 
     336          188 :       CALL linres_calculation_low(qs_env)
     337              : 
     338          188 :       CALL timestop(handle)
     339              : 
     340          188 :    END SUBROUTINE linres_calculation
     341              : 
     342              : ! **************************************************************************************************
     343              : !> \brief Linear response can be called as run type or as post scf calculation
     344              : !>      Initialize the perturbation environment
     345              : !>      Define which properties is to be calculated
     346              : !>      Start up the optimization of the response density and wfn
     347              : !> \param qs_env ...
     348              : !> \par History
     349              : !>      06.2005 created [MI]
     350              : !>      02.2013 added polarizability section [SL]
     351              : !> \author MI
     352              : ! **************************************************************************************************
     353        21389 :    SUBROUTINE linres_calculation_low(qs_env)
     354              : 
     355              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     356              : 
     357              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation_low'
     358              : 
     359              :       INTEGER                                            :: every_n_step, handle, iounit
     360              :       LOGICAL                                            :: dcdr_present, do_apt_fd, epr_present, &
     361              :                                                             issc_present, lr_calculation, &
     362              :                                                             nmr_present, polar_present, vcd_present
     363              :       TYPE(cp_logger_type), POINTER                      :: logger
     364              :       TYPE(dft_control_type), POINTER                    :: dft_control
     365              :       TYPE(linres_control_type), POINTER                 :: linres_control
     366              :       TYPE(qs_p_env_type)                                :: p_env
     367              :       TYPE(section_vals_type), POINTER                   :: lr_section, prop_section
     368              : 
     369        21389 :       CALL timeset(routineN, handle)
     370              : 
     371              :       lr_calculation = .FALSE.
     372        21389 :       nmr_present = .FALSE.
     373        21389 :       epr_present = .FALSE.
     374        21389 :       issc_present = .FALSE.
     375        21389 :       polar_present = .FALSE.
     376        21389 :       dcdr_present = .FALSE.
     377              :       do_apt_fd = .FALSE.
     378              : 
     379        21389 :       NULLIFY (dft_control, linres_control, logger, prop_section, lr_section)
     380              : 
     381        21389 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     382        21389 :       CALL section_vals_get(lr_section, explicit=lr_calculation)
     383              : 
     384        21389 :       CALL section_vals_val_get(lr_section, "DCDR%APT_FD", explicit=do_apt_fd)
     385        21389 :       IF (do_apt_fd) THEN
     386           14 :          CALL timestop(handle)
     387           14 :          RETURN
     388              :       END IF
     389              : 
     390        21375 :       logger => cp_get_default_logger()
     391              : 
     392        21375 :       CALL section_vals_val_get(lr_section, "EVERY_N_STEP", i_val=every_n_step)
     393              : 
     394        21375 :       IF (lr_calculation .AND. MODULO(qs_env%sim_step, every_n_step) == 0) THEN
     395          318 :          CALL linres_init(lr_section, p_env, qs_env)
     396              :          iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     397          318 :                                        extension=".linresLog")
     398              :          CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     399          318 :                          linres_control=linres_control)
     400              : 
     401              :          ! The type of perturbation has not been defined yet
     402          318 :          linres_control%property = lr_none
     403              : 
     404              :          ! We do NMR or EPR, then compute the current response
     405          318 :          prop_section => section_vals_get_subs_vals(lr_section, "NMR")
     406          318 :          CALL section_vals_get(prop_section, explicit=nmr_present)
     407          318 :          prop_section => section_vals_get_subs_vals(lr_section, "EPR")
     408          318 :          CALL section_vals_get(prop_section, explicit=epr_present)
     409              : 
     410          318 :          IF (nmr_present .OR. epr_present) THEN
     411              :             CALL nmr_epr_linres(linres_control, qs_env, p_env, dft_control, &
     412          174 :                                 nmr_present, epr_present, iounit)
     413              :          END IF
     414              : 
     415              :          ! We do the indirect spin-spin coupling calculation
     416          318 :          prop_section => section_vals_get_subs_vals(lr_section, "SPINSPIN")
     417          318 :          CALL section_vals_get(prop_section, explicit=issc_present)
     418              : 
     419          318 :          IF (issc_present) THEN
     420           12 :             CALL issc_linres(linres_control, qs_env, p_env, dft_control)
     421              :          END IF
     422              : 
     423              :          ! We do the polarizability calculation
     424          318 :          prop_section => section_vals_get_subs_vals(lr_section, "POLAR")
     425          318 :          CALL section_vals_get(prop_section, explicit=polar_present)
     426          318 :          IF (polar_present) THEN
     427          108 :             CALL polar_linres(qs_env, p_env)
     428              :          END IF
     429              : 
     430              :          ! Nuclear Position Perturbation
     431          318 :          prop_section => section_vals_get_subs_vals(lr_section, "dcdr")
     432          318 :          CALL section_vals_get(prop_section, explicit=dcdr_present)
     433              : 
     434          318 :          IF (dcdr_present) THEN
     435           22 :             CALL dcdr_linres(qs_env, p_env)
     436              :          END IF
     437              : 
     438              :          ! VCD
     439          318 :          prop_section => section_vals_get_subs_vals(lr_section, "VCD")
     440          318 :          CALL section_vals_get(prop_section, explicit=vcd_present)
     441              : 
     442          318 :          IF (vcd_present) THEN
     443            2 :             CALL vcd_linres(qs_env, p_env)
     444              :          END IF
     445              : 
     446              :          ! Other possible LR calculations can be introduced here
     447              : 
     448          318 :          CALL p_env_release(p_env)
     449              : 
     450          318 :          IF (iounit > 0) THEN
     451              :             WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") &
     452          159 :                REPEAT("=", 79), &
     453          159 :                "ENDED LINRES CALCULATION", &
     454          318 :                REPEAT("=", 79)
     455              :          END IF
     456              :          CALL cp_print_key_finished_output(iounit, logger, lr_section, &
     457          318 :                                            "PRINT%PROGRAM_RUN_INFO")
     458              :       END IF
     459              : 
     460        21375 :       CALL timestop(handle)
     461              : 
     462       106945 :    END SUBROUTINE linres_calculation_low
     463              : 
     464              : ! **************************************************************************************************
     465              : !> \brief Initialize some general settings like the p_env
     466              : !>      Localize the psi0 if required
     467              : !> \param lr_section ...
     468              : !> \param p_env ...
     469              : !> \param qs_env ...
     470              : !> \par History
     471              : !>      06.2005 created [MI]
     472              : !> \author MI
     473              : !> \note
     474              : !>      - The localization should probably be always for all the occupied states
     475              : ! **************************************************************************************************
     476         1908 :    SUBROUTINE linres_init(lr_section, p_env, qs_env)
     477              : 
     478              :       TYPE(section_vals_type), POINTER                   :: lr_section
     479              :       TYPE(qs_p_env_type), INTENT(OUT)                   :: p_env
     480              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     481              : 
     482              :       INTEGER                                            :: iounit, ispin
     483              :       LOGICAL                                            :: do_it
     484              :       TYPE(cp_logger_type), POINTER                      :: logger
     485          318 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, rho_ao
     486              :       TYPE(dft_control_type), POINTER                    :: dft_control
     487              :       TYPE(linres_control_type), POINTER                 :: linres_control
     488          318 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     489              :       TYPE(qs_rho_type), POINTER                         :: rho
     490              :       TYPE(section_vals_type), POINTER                   :: loc_section
     491              : 
     492          318 :       NULLIFY (logger)
     493          318 :       logger => cp_get_default_logger()
     494              :       iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     495          318 :                                     extension=".linresLog")
     496          318 :       NULLIFY (dft_control, linres_control, loc_section, rho, mos, matrix_ks, rho_ao)
     497              : 
     498          318 :       ALLOCATE (linres_control)
     499          318 :       CALL set_qs_env(qs_env=qs_env, linres_control=linres_control)
     500              :       CALL get_qs_env(qs_env=qs_env, &
     501          318 :                       dft_control=dft_control, matrix_ks=matrix_ks, mos=mos, rho=rho)
     502          318 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     503              : 
     504              :       ! Localized Psi0 are required when the position operator has to be defined (nmr)
     505          318 :       loc_section => section_vals_get_subs_vals(lr_section, "LOCALIZE")
     506              :       CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", &
     507          318 :                                 l_val=linres_control%localized_psi0)
     508          318 :       IF (linres_control%localized_psi0) THEN
     509          190 :          IF (iounit > 0) THEN
     510              :             WRITE (UNIT=iounit, FMT="(/,T3,A,A)") &
     511           95 :                "Localization of ground state orbitals", &
     512          190 :                " before starting linear response calculation"
     513              :          END IF
     514              : 
     515          190 :          CALL linres_localize(qs_env, linres_control, dft_control%nspins)
     516              : 
     517          458 :          DO ispin = 1, dft_control%nspins
     518          458 :             CALL calculate_density_matrix(mos(ispin), rho_ao(ispin)%matrix)
     519              :          END DO
     520              :          ! ** update qs_env%rho
     521          190 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     522              :       END IF
     523              : 
     524          318 :       CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart)
     525          318 :       CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
     526          318 :       CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps)
     527          318 :       CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter)
     528          318 :       CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every)
     529          318 :       CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
     530          318 :       CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
     531              : 
     532          318 :       IF (iounit > 0) THEN
     533              :          WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") &
     534          159 :             REPEAT("=", 79), &
     535          159 :             "START LINRES CALCULATION", &
     536          318 :             REPEAT("=", 79)
     537              : 
     538              :          WRITE (UNIT=iounit, FMT="(T2,A)") &
     539          159 :             "LINRES| Properties to be calculated:"
     540          159 :          CALL section_vals_val_get(lr_section, "NMR%_SECTION_PARAMETERS_", l_val=do_it)
     541          159 :          IF (do_it) WRITE (UNIT=iounit, FMT="(T62,A)") "NMR Chemical Shift"
     542          159 :          CALL section_vals_val_get(lr_section, "EPR%_SECTION_PARAMETERS_", l_val=do_it)
     543          159 :          IF (do_it) WRITE (UNIT=iounit, FMT="(T68,A)") "EPR g Tensor"
     544          159 :          CALL section_vals_val_get(lr_section, "SPINSPIN%_SECTION_PARAMETERS_", l_val=do_it)
     545          159 :          IF (do_it) WRITE (UNIT=iounit, FMT="(T43,A)") "Indirect spin-spin coupling constants"
     546          159 :          CALL section_vals_val_get(lr_section, "POLAR%_SECTION_PARAMETERS_", l_val=do_it)
     547          159 :          IF (do_it) WRITE (UNIT=iounit, FMT="(T57,A)") "Electric Polarizability"
     548              : 
     549          159 :          IF (linres_control%localized_psi0) WRITE (UNIT=iounit, FMT="(T2,A,T65,A)") &
     550           95 :             "LINRES|", " LOCALIZED PSI0"
     551              : 
     552              :          WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     553          159 :             "LINRES| Optimization algorithm", " Conjugate Gradients"
     554              : 
     555          160 :          SELECT CASE (linres_control%preconditioner_type)
     556              :          CASE (ot_precond_none)
     557              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     558            1 :                "LINRES| Preconditioner", "                NONE"
     559              :          CASE (ot_precond_full_single)
     560              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     561            2 :                "LINRES| Preconditioner", "         FULL_SINGLE"
     562              :          CASE (ot_precond_full_kinetic)
     563              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     564            3 :                "LINRES| Preconditioner", "        FULL_KINETIC"
     565              :          CASE (ot_precond_s_inverse)
     566              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     567           12 :                "LINRES| Preconditioner", "      FULL_S_INVERSE"
     568              :          CASE (ot_precond_full_single_inverse)
     569              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     570           32 :                "LINRES| Preconditioner", " FULL_SINGLE_INVERSE"
     571              :          CASE (ot_precond_full_all)
     572              :             WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
     573          109 :                "LINRES| Preconditioner", "            FULL_ALL"
     574              :          CASE DEFAULT
     575          159 :             CPABORT("Preconditioner NYI")
     576              :          END SELECT
     577              : 
     578              :          WRITE (UNIT=iounit, FMT="(T2,A,T72,ES8.1)") &
     579          159 :             "LINRES| EPS", linres_control%eps
     580              :          WRITE (UNIT=iounit, FMT="(T2,A,T72,I8)") &
     581          159 :             "LINRES| MAX_ITER", linres_control%max_iter
     582              :       END IF
     583              : 
     584              :       !------------------!
     585              :       ! create the p_env !
     586              :       !------------------!
     587          318 :       CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE., linres_control=linres_control)
     588              : 
     589              :       ! update the m_epsilon matrix
     590          318 :       CALL p_env_psi0_changed(p_env, qs_env)
     591              : 
     592          318 :       p_env%new_preconditioner = .TRUE.
     593              :       CALL cp_print_key_finished_output(iounit, logger, lr_section, &
     594          318 :                                         "PRINT%PROGRAM_RUN_INFO")
     595              : 
     596          318 :    END SUBROUTINE linres_init
     597              : 
     598              : ! **************************************************************************************************
     599              : !> \brief ...
     600              : !> \param linres_control ...
     601              : !> \param qs_env ...
     602              : !> \param p_env ...
     603              : !> \param dft_control ...
     604              : !> \param nmr_present ...
     605              : !> \param epr_present ...
     606              : !> \param iounit ...
     607              : ! **************************************************************************************************
     608          174 :    SUBROUTINE nmr_epr_linres(linres_control, qs_env, p_env, dft_control, nmr_present, epr_present, iounit)
     609              : 
     610              :       TYPE(linres_control_type), POINTER                 :: linres_control
     611              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     612              :       TYPE(qs_p_env_type)                                :: p_env
     613              :       TYPE(dft_control_type), POINTER                    :: dft_control
     614              :       LOGICAL                                            :: nmr_present, epr_present
     615              :       INTEGER                                            :: iounit
     616              : 
     617              :       INTEGER                                            :: iB
     618              :       LOGICAL                                            :: do_qmmm
     619              :       TYPE(current_env_type)                             :: current_env
     620              :       TYPE(epr_env_type)                                 :: epr_env
     621              :       TYPE(nmr_env_type)                                 :: nmr_env
     622              : 
     623          174 :       linres_control%property = lr_current
     624              : 
     625          174 :       CALL cite_reference(Weber2009)
     626              : 
     627          174 :       IF (.NOT. linres_control%localized_psi0) THEN
     628              :          CALL cp_abort(__LOCATION__, &
     629              :                        "Are you sure that you want to calculate the chemical "// &
     630            0 :                        "shift without localized psi0?")
     631              :          CALL linres_localize(qs_env, linres_control, &
     632            0 :                               dft_control%nspins, centers_only=.TRUE.)
     633              :       END IF
     634          174 :       IF (dft_control%nspins /= 2 .AND. epr_present) THEN
     635            0 :          CPABORT("LSD is needed to perform a g tensor calculation!")
     636              :       END IF
     637              :       !
     638              :       !Initialize the current environment
     639          174 :       do_qmmm = .FALSE.
     640          174 :       IF (qs_env%qmmm) do_qmmm = .TRUE.
     641          174 :       current_env%do_qmmm = do_qmmm
     642              :       !current_env%prop='nmr'
     643          174 :       CALL current_env_init(current_env, qs_env)
     644          174 :       CALL current_operators(current_env, qs_env)
     645          174 :       CALL current_response(current_env, p_env, qs_env)
     646              :       !
     647          174 :       IF (current_env%all_pert_op_done) THEN
     648              :          !Initialize the nmr environment
     649          174 :          IF (nmr_present) THEN
     650          160 :             CALL nmr_env_init(nmr_env, qs_env)
     651              :          END IF
     652              :          !
     653              :          !Initialize the epr environment
     654          174 :          IF (epr_present) THEN
     655           14 :             CALL epr_env_init(epr_env, qs_env)
     656           14 :             CALL epr_g_zke(epr_env, qs_env)
     657           14 :             CALL epr_nablavks(epr_env, qs_env)
     658              :          END IF
     659              :          !
     660              :          ! Build the rs_gauge if needed
     661              :          !CALL current_set_gauge(current_env,qs_env)
     662              :          !
     663              :          ! Loop over field direction
     664          696 :          DO iB = 1, 3
     665              :             !
     666              :             ! Build current response and succeptibility
     667          522 :             CALL current_build_current(current_env, qs_env, iB)
     668          522 :             CALL current_build_chi(current_env, qs_env, iB)
     669              :             !
     670              :             ! Compute NMR shift
     671          522 :             IF (nmr_present) THEN
     672          480 :                CALL nmr_shift(nmr_env, current_env, qs_env, iB)
     673              :             END IF
     674              :             !
     675              :             ! Compute EPR
     676          696 :             IF (epr_present) THEN
     677           42 :                CALL epr_ind_magnetic_field(epr_env, current_env, qs_env, iB)
     678           42 :                CALL epr_g_so(epr_env, current_env, qs_env, iB)
     679           42 :                CALL epr_g_soo(epr_env, current_env, qs_env, iB)
     680              :             END IF
     681              :          END DO
     682              :          !
     683              :          ! Finalized the nmr environment
     684          174 :          IF (nmr_present) THEN
     685          160 :             CALL nmr_shift_print(nmr_env, current_env, qs_env)
     686          160 :             CALL nmr_env_cleanup(nmr_env)
     687              :          END IF
     688              :          !
     689              :          ! Finalized the epr environment
     690          174 :          IF (epr_present) THEN
     691           14 :             CALL epr_g_print(epr_env, qs_env)
     692           14 :             CALL epr_env_cleanup(epr_env)
     693              :          END IF
     694              :          !
     695              :       ELSE
     696            0 :          IF (iounit > 0) THEN
     697              :             WRITE (iounit, "(T10,A,/T20,A,/)") &
     698            0 :                "CURRENT: Not all responses to perturbation operators could be calculated.", &
     699            0 :                " Hence: NO nmr and NO epr possible."
     700              :          END IF
     701              :       END IF
     702              :       ! Finalized the current environment
     703          174 :       CALL current_env_cleanup(current_env)
     704              : 
     705        12702 :    END SUBROUTINE nmr_epr_linres
     706              : 
     707              : ! **************************************************************************************************
     708              : !> \brief ...
     709              : !> \param linres_control ...
     710              : !> \param qs_env ...
     711              : !> \param p_env ...
     712              : !> \param dft_control ...
     713              : ! **************************************************************************************************
     714           12 :    SUBROUTINE issc_linres(linres_control, qs_env, p_env, dft_control)
     715              : 
     716              :       TYPE(linres_control_type), POINTER                 :: linres_control
     717              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     718              :       TYPE(qs_p_env_type)                                :: p_env
     719              :       TYPE(dft_control_type), POINTER                    :: dft_control
     720              : 
     721              :       INTEGER                                            :: iatom
     722              :       LOGICAL                                            :: do_qmmm
     723              :       TYPE(current_env_type)                             :: current_env
     724              :       TYPE(issc_env_type)                                :: issc_env
     725              : 
     726           12 :       linres_control%property = lr_current
     727           12 :       IF (.NOT. linres_control%localized_psi0) THEN
     728              :          CALL cp_abort(__LOCATION__, &
     729              :                        "Are you sure that you want to calculate the chemical "// &
     730            0 :                        "shift without localized psi0?")
     731              :          CALL linres_localize(qs_env, linres_control, &
     732            0 :                               dft_control%nspins, centers_only=.TRUE.)
     733              :       END IF
     734              :       !
     735              :       !Initialize the current environment
     736           12 :       do_qmmm = .FALSE.
     737           12 :       IF (qs_env%qmmm) do_qmmm = .TRUE.
     738           12 :       current_env%do_qmmm = do_qmmm
     739              :       !current_env%prop='issc'
     740              :       !CALL current_env_init(current_env,qs_env)
     741              :       !CALL current_response(current_env,p_env,qs_env)
     742              :       !
     743              :       !Initialize the issc environment
     744           12 :       CALL issc_env_init(issc_env, qs_env)
     745              :       !
     746              :       ! Loop over atoms
     747           56 :       DO iatom = 1, issc_env%issc_natms
     748           44 :          CALL issc_operators(issc_env, qs_env, iatom)
     749           44 :          CALL issc_response(issc_env, p_env, qs_env)
     750           56 :          CALL issc_issc(issc_env, qs_env, iatom)
     751              :       END DO
     752              :       !
     753              :       ! Finalized the issc environment
     754           12 :       CALL issc_print(issc_env, qs_env)
     755           12 :       CALL issc_env_cleanup(issc_env)
     756              : 
     757          888 :    END SUBROUTINE issc_linres
     758              : 
     759              : ! **************************************************************************************************
     760              : !> \brief ...
     761              : !> \param qs_env ...
     762              : !> \param p_env ...
     763              : !> \par History
     764              : !>      06.2018 polar_env integrated into qs_env (MK)
     765              : ! **************************************************************************************************
     766          108 :    SUBROUTINE polar_linres(qs_env, p_env)
     767              : 
     768              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     769              :       TYPE(qs_p_env_type)                                :: p_env
     770              : 
     771          108 :       CALL polar_env_init(qs_env)
     772          108 :       CALL polar_operators(qs_env)
     773          108 :       CALL polar_response(p_env, qs_env)
     774          108 :       CALL polar_polar(qs_env)
     775          108 :       CALL polar_print(qs_env)
     776              : 
     777          108 :    END SUBROUTINE polar_linres
     778              : 
     779              : END MODULE qs_linres_module
        

Generated by: LCOV version 2.0-1