LCOV - code coverage report
Current view: top level - src - qs_vcd_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 88.6 % 446 395
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              : MODULE qs_vcd_utils
       9              :    USE cell_types,                      ONLY: cell_type
      10              :    USE commutator_rpnl,                 ONLY: build_com_mom_nl
      11              :    USE cp_control_types,                ONLY: dft_control_type
      12              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      13              :                                               dbcsr_create,&
      14              :                                               dbcsr_init_p,&
      15              :                                               dbcsr_p_type,&
      16              :                                               dbcsr_set,&
      17              :                                               dbcsr_type_antisymmetric,&
      18              :                                               dbcsr_type_no_symmetry
      19              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      20              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      21              :                                               dbcsr_deallocate_matrix_set
      22              :    USE cp_files,                        ONLY: close_file,&
      23              :                                               open_file
      24              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25              :                                               cp_fm_get_info,&
      26              :                                               cp_fm_get_submatrix,&
      27              :                                               cp_fm_release,&
      28              :                                               cp_fm_set_submatrix,&
      29              :                                               cp_fm_type
      30              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      31              :                                               cp_logger_get_default_io_unit,&
      32              :                                               cp_logger_type,&
      33              :                                               cp_to_string
      34              :    USE cp_output_handling,              ONLY: cp_p_file,&
      35              :                                               cp_print_key_finished_output,&
      36              :                                               cp_print_key_generate_filename,&
      37              :                                               cp_print_key_should_output,&
      38              :                                               cp_print_key_unit_nr
      39              :    USE cp_result_methods,               ONLY: get_results
      40              :    USE cp_result_types,                 ONLY: cp_result_type
      41              :    USE input_constants,                 ONLY: use_mom_ref_user
      42              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      43              :                                               section_vals_type,&
      44              :                                               section_vals_val_get
      45              :    USE kinds,                           ONLY: default_path_length,&
      46              :                                               default_string_length,&
      47              :                                               dp
      48              :    USE message_passing,                 ONLY: mp_para_env_type
      49              :    USE molecule_types,                  ONLY: molecule_type
      50              :    USE moments_utils,                   ONLY: get_reference_point
      51              :    USE orbital_pointers,                ONLY: init_orbital_pointers
      52              :    USE particle_types,                  ONLY: particle_type
      53              :    USE qs_dcdr_utils,                   ONLY: dcdr_env_cleanup,&
      54              :                                               dcdr_env_init
      55              :    USE qs_environment_types,            ONLY: get_qs_env,&
      56              :                                               qs_environment_type
      57              :    USE qs_kind_types,                   ONLY: qs_kind_type
      58              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      59              :    USE qs_linres_types,                 ONLY: vcd_env_type
      60              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      61              :                                               mo_set_type
      62              :    USE qs_moments,                      ONLY: build_local_moments_der_matrix
      63              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      64              :    USE qs_operators_ao,                 ONLY: build_lin_mom_matrix
      65              :    USE qs_vcd_ao,                       ONLY: build_com_rpnl_r,&
      66              :                                               build_matrix_r_vhxc,&
      67              :                                               build_rcore_matrix,&
      68              :                                               build_rpnl_matrix,&
      69              :                                               build_tr_matrix
      70              :    USE string_utilities,                ONLY: xstring
      71              : #include "./base/base_uses.f90"
      72              : 
      73              :    IMPLICIT NONE
      74              : 
      75              :    PRIVATE
      76              :    PUBLIC :: vcd_env_cleanup, vcd_env_init
      77              :    PUBLIC :: vcd_read_restart, vcd_write_restart
      78              :    PUBLIC :: vcd_print
      79              : 
      80              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd_utils'
      81              : 
      82              :    REAL(dp), DIMENSION(3, 3, 3), PARAMETER  :: Levi_Civita = RESHAPE((/ &
      83              :                                                           0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
      84              :                                                           0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
      85              :                                              0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), (/3, 3, 3/))
      86              : 
      87              : CONTAINS
      88              : 
      89              : ! *****************************************************************************
      90              : !> \brief Initialize the vcd environment
      91              : !> \param vcd_env ...
      92              : !> \param qs_env ...
      93              : !> \author Edward Ditler
      94              : ! **************************************************************************************************
      95            2 :    SUBROUTINE vcd_env_init(vcd_env, qs_env)
      96              :       TYPE(vcd_env_type), TARGET                         :: vcd_env
      97              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      98              : 
      99              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_env_init'
     100              : 
     101              :       INTEGER                                            :: handle, i, idir, ispin, j, natom, &
     102              :                                                             nspins, output_unit, reference, &
     103              :                                                             unit_number
     104              :       LOGICAL                                            :: explicit
     105            2 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
     106              :       TYPE(cell_type), POINTER                           :: cell
     107              :       TYPE(cp_logger_type), POINTER                      :: logger
     108            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, my_matrix_hr_1d
     109              :       TYPE(dft_control_type), POINTER                    :: dft_control
     110              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     111            2 :          POINTER                                         :: sab_all, sab_orb, sap_ppnl
     112            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     113            2 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     114              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     115              :       TYPE(section_vals_type), POINTER                   :: lr_section, vcd_section
     116              : 
     117            2 :       CALL timeset(routineN, handle)
     118            2 :       vcd_env%do_mfp = .FALSE.
     119              : 
     120              :       ! Set up the logger
     121            2 :       NULLIFY (logger, vcd_section, lr_section)
     122            2 :       logger => cp_get_default_logger()
     123            2 :       vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
     124              :       vcd_env%output_unit = cp_print_key_unit_nr(logger, vcd_section, "PRINT%VCD", &
     125              :                                                  extension=".data", middle_name="vcd", log_filename=.FALSE., &
     126            2 :                                                  file_position="REWIND", file_status="REPLACE")
     127              : 
     128            2 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     129              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     130            2 :                                          extension=".linresLog")
     131            2 :       unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
     132              : 
     133              :       ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
     134            2 :       CALL dcdr_env_init(vcd_env%dcdr_env, qs_env)
     135              :       ! vcd_env%dcdr_env%output_unit = vcd_env%output_unit
     136              : 
     137            2 :       IF (output_unit > 0) THEN
     138            1 :          WRITE (output_unit, "(/,T20,A,/)") "*** Start NVPT/MFPT calculation ***"
     139              :       END IF
     140              : 
     141              :       ! Just to make sure. The memory requirements are tiny.
     142            2 :       CALL init_orbital_pointers(12)
     143              : 
     144            2 :       CALL section_vals_val_get(vcd_section, "DISTRIBUTED_ORIGIN", l_val=vcd_env%distributed_origin)
     145            2 :       CALL section_vals_val_get(vcd_section, "ORIGIN_DEPENDENT_MFP", l_val=vcd_env%origin_dependent_op_mfp)
     146              : 
     147              :       ! Reference point
     148            8 :       vcd_env%magnetic_origin = 0._dp
     149            8 :       vcd_env%spatial_origin = 0._dp
     150              :       ! Get the magnetic origin from the input
     151            2 :       CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN", i_val=reference)
     152            2 :       CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", explicit=explicit)
     153            2 :       IF (explicit) THEN
     154            0 :          CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", r_vals=ref_point)
     155              :       ELSE
     156            2 :          IF (reference == use_mom_ref_user) &
     157            0 :             CPABORT("User-defined reference point should be given explicitly")
     158              :       END IF
     159              : 
     160              :       CALL get_reference_point(rpoint=vcd_env%magnetic_origin, qs_env=qs_env, &
     161              :                                reference=reference, &
     162            2 :                                ref_point=ref_point)
     163              : 
     164              :       ! Get the spatial origin from the input
     165            2 :       CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN", i_val=reference)
     166            2 :       CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", explicit=explicit)
     167            2 :       IF (explicit) THEN
     168            0 :          CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", r_vals=ref_point)
     169              :       ELSE
     170            2 :          IF (reference == use_mom_ref_user) &
     171            0 :             CPABORT("User-defined reference point should be given explicitly")
     172              :       END IF
     173              : 
     174              :       CALL get_reference_point(rpoint=vcd_env%spatial_origin, qs_env=qs_env, &
     175              :                                reference=reference, &
     176            2 :                                ref_point=ref_point)
     177              : 
     178            8 :       IF (vcd_env%distributed_origin .AND. ANY(vcd_env%magnetic_origin /= vcd_env%spatial_origin)) THEN
     179            0 :          CPWARN("The magnetic and spatial origins don't match")
     180              :          ! This is fine for NVP but will give unphysical results for MFP.
     181              :       END IF
     182              : 
     183            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
     184            1 :          'The reference point is', vcd_env%dcdr_env%ref_point
     185            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
     186            1 :          'The magnetic origin is', vcd_env%magnetic_origin
     187            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
     188            1 :          'The velocity origin is', vcd_env%spatial_origin
     189              : 
     190            8 :       vcd_env%magnetic_origin_atom = vcd_env%magnetic_origin
     191            8 :       vcd_env%spatial_origin_atom = vcd_env%spatial_origin
     192              : 
     193              :       CALL get_qs_env(qs_env=qs_env, &
     194              :                       ks_env=ks_env, &
     195              :                       dft_control=dft_control, &
     196              :                       sab_orb=sab_orb, &
     197              :                       sab_all=sab_all, &
     198              :                       sap_ppnl=sap_ppnl, &
     199              :                       particle_set=particle_set, &
     200              :                       matrix_ks=matrix_ks, &
     201              :                       cell=cell, &
     202            2 :                       qs_kind_set=qs_kind_set)
     203              : 
     204            2 :       natom = SIZE(particle_set)
     205            2 :       nspins = dft_control%nspins
     206              : 
     207            6 :       ALLOCATE (vcd_env%apt_el_nvpt(3, 3, natom))
     208            4 :       ALLOCATE (vcd_env%apt_nuc_nvpt(3, 3, natom))
     209            4 :       ALLOCATE (vcd_env%apt_total_nvpt(3, 3, natom))
     210            4 :       ALLOCATE (vcd_env%aat_atom_nvpt(3, 3, natom))
     211            4 :       ALLOCATE (vcd_env%aat_atom_mfp(3, 3, natom))
     212           80 :       vcd_env%apt_el_nvpt = 0._dp
     213           80 :       vcd_env%apt_nuc_nvpt = 0._dp
     214           80 :       vcd_env%apt_total_nvpt = 0._dp
     215           80 :       vcd_env%aat_atom_nvpt = 0._dp
     216           80 :       vcd_env%aat_atom_mfp = 0._dp
     217              : 
     218            8 :       ALLOCATE (vcd_env%dCV(nspins))
     219            6 :       ALLOCATE (vcd_env%dCV_prime(nspins))
     220            6 :       ALLOCATE (vcd_env%op_dV(nspins))
     221            6 :       ALLOCATE (vcd_env%op_dB(nspins))
     222            4 :       DO ispin = 1, nspins
     223            2 :          CALL cp_fm_create(vcd_env%dCV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     224            2 :          CALL cp_fm_create(vcd_env%dCV_prime(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     225            2 :          CALL cp_fm_create(vcd_env%op_dV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     226            4 :          CALL cp_fm_create(vcd_env%op_dB(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     227              :       END DO
     228              : 
     229            8 :       ALLOCATE (vcd_env%dCB(3))
     230            8 :       ALLOCATE (vcd_env%dCB_prime(3))
     231            8 :       DO i = 1, 3
     232            6 :          CALL cp_fm_create(vcd_env%dCB(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     233            8 :          CALL cp_fm_create(vcd_env%dCB_prime(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
     234              :       END DO
     235              : 
     236              :       ! DBCSR matrices
     237            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%moments_der, 9, 3)
     238            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_right, 9, 3)
     239            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_left, 9, 3)
     240            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_difdip2, 3, 3)
     241            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdV, 3)
     242            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdB, 3)
     243            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hxc_dsdv, nspins)
     244              : 
     245            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%hcom, 3)
     246            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rcomr, 3, 3)
     247            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rrcom, 3, 3)
     248            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dcom, 3, 3)
     249            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hr, nspins, 3)
     250            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rh, nspins, 3)
     251            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_drpnl, 3)
     252            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao, 3)
     253            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao_delta, 3)
     254            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxrv, 3)
     255            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_rxvr, 3, 3)
     256            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxvr_r, 3, 3)
     257            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_doublecom, 3, 3)
     258              : 
     259            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp_33, 3, 3)
     260            2 :       CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp2_33, 3, 3)
     261           20 :       DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
     262           74 :          DO idir = 1, 3 ! d/dx, d/dy, d/dz
     263           54 :             CALL dbcsr_init_p(vcd_env%moments_der(i, idir)%matrix)
     264           54 :             CALL dbcsr_init_p(vcd_env%moments_der_right(i, idir)%matrix)
     265           54 :             CALL dbcsr_init_p(vcd_env%moments_der_left(i, idir)%matrix)
     266              : 
     267              :             CALL dbcsr_create(vcd_env%moments_der(i, idir)%matrix, template=matrix_ks(1)%matrix, &
     268           54 :                               matrix_type=dbcsr_type_antisymmetric)
     269           54 :             CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%moments_der(i, idir)%matrix, sab_orb)
     270           54 :             CALL dbcsr_set(vcd_env%moments_der(i, idir)%matrix, 0.0_dp)
     271              : 
     272              :             ! And the ones which will be multiplied by delta_(mu/nu)
     273              :             CALL dbcsr_copy(vcd_env%moments_der_right(i, idir)%matrix, &
     274           54 :                             vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     275              :             CALL dbcsr_copy(vcd_env%moments_der_left(i, idir)%matrix, &
     276           72 :                             vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     277              :          END DO
     278              :       END DO
     279              : 
     280            8 :       DO i = 1, 3
     281           24 :          DO j = 1, 3
     282           18 :             CALL dbcsr_init_p(vcd_env%matrix_difdip2(i, j)%matrix)
     283           18 :             CALL dbcsr_copy(vcd_env%matrix_difdip2(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     284           18 :             CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0.0_dp)
     285              : 
     286           18 :             CALL dbcsr_init_p(vcd_env%matrix_nosym_temp_33(i, j)%matrix)
     287              :             CALL dbcsr_create(vcd_env%matrix_nosym_temp_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
     288           18 :                               matrix_type=dbcsr_type_no_symmetry)
     289           18 :             CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp_33(i, j)%matrix, sab_all)
     290           18 :             CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(i, j)%matrix, 0._dp)
     291              : 
     292           18 :             CALL dbcsr_init_p(vcd_env%matrix_nosym_temp2_33(i, j)%matrix)
     293              :             CALL dbcsr_create(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
     294           18 :                               matrix_type=dbcsr_type_no_symmetry)
     295           18 :             CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, sab_all)
     296           24 :             CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, 0._dp)
     297              : 
     298              :          END DO
     299            6 :          CALL dbcsr_init_p(vcd_env%matrix_dSdV(i)%matrix)
     300            6 :          CALL dbcsr_copy(vcd_env%matrix_dSdV(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     301            6 :          CALL dbcsr_init_p(vcd_env%matrix_dSdB(i)%matrix)
     302            8 :          CALL dbcsr_copy(vcd_env%matrix_dSdB(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     303              :       END DO
     304              : 
     305            4 :       DO ispin = 1, nspins
     306            2 :          CALL dbcsr_init_p(vcd_env%matrix_hxc_dsdv(ispin)%matrix)
     307            4 :          CALL dbcsr_copy(vcd_env%matrix_hxc_dsdv(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     308              :       END DO
     309              : 
     310              :       ! Things for op_dV
     311              :       ! lin_mom
     312            8 :       DO i = 1, 3
     313            6 :          CALL dbcsr_init_p(vcd_env%dipvel_ao(i)%matrix)
     314            6 :          CALL dbcsr_copy(vcd_env%dipvel_ao(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     315              : 
     316            6 :          CALL dbcsr_init_p(vcd_env%dipvel_ao_delta(i)%matrix)
     317            8 :          CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     318              :       END DO
     319              : 
     320              :       ! [V, r]
     321            8 :       DO i = 1, 3
     322            6 :          CALL dbcsr_init_p(vcd_env%hcom(i)%matrix)
     323              :          CALL dbcsr_create(vcd_env%hcom(i)%matrix, template=matrix_ks(1)%matrix, &
     324            6 :                            matrix_type=dbcsr_type_antisymmetric)
     325            6 :          CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%hcom(i)%matrix, sab_orb)
     326              : 
     327            6 :          CALL dbcsr_init_p(vcd_env%matrix_rxrv(i)%matrix)
     328              :          CALL dbcsr_create(vcd_env%matrix_rxrv(i)%matrix, template=matrix_ks(1)%matrix, &
     329            6 :                            matrix_type=dbcsr_type_antisymmetric)
     330            6 :          CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_rxrv(i)%matrix, sab_orb)
     331              : 
     332           26 :          DO j = 1, 3
     333           18 :             CALL dbcsr_init_p(vcd_env%matrix_rcomr(i, j)%matrix)
     334           18 :             CALL dbcsr_copy(vcd_env%matrix_rcomr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     335           18 :             CALL dbcsr_init_p(vcd_env%matrix_rrcom(i, j)%matrix)
     336           18 :             CALL dbcsr_copy(vcd_env%matrix_rrcom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     337           18 :             CALL dbcsr_init_p(vcd_env%matrix_dcom(i, j)%matrix)
     338           18 :             CALL dbcsr_copy(vcd_env%matrix_dcom(i, j)%matrix, matrix_ks(1)%matrix)
     339           18 :             CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0._dp)
     340              : 
     341           18 :             CALL dbcsr_init_p(vcd_env%matrix_r_rxvr(i, j)%matrix)
     342           18 :             CALL dbcsr_copy(vcd_env%matrix_r_rxvr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     343           18 :             CALL dbcsr_set(vcd_env%matrix_r_rxvr(i, j)%matrix, 0._dp)
     344              : 
     345           18 :             CALL dbcsr_init_p(vcd_env%matrix_rxvr_r(i, j)%matrix)
     346           18 :             CALL dbcsr_copy(vcd_env%matrix_rxvr_r(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     347           18 :             CALL dbcsr_set(vcd_env%matrix_rxvr_r(i, j)%matrix, 0._dp)
     348              : 
     349           18 :             CALL dbcsr_init_p(vcd_env%matrix_r_doublecom(i, j)%matrix)
     350           18 :             CALL dbcsr_copy(vcd_env%matrix_r_doublecom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     351           24 :             CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
     352              :          END DO
     353              :       END DO
     354              : 
     355              :       ! matrix_hr: nonsymmetric dbcsr matrix
     356            4 :       DO ispin = 1, nspins
     357           10 :          DO i = 1, 3
     358            6 :             CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
     359            6 :             CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     360              : 
     361            6 :             CALL dbcsr_init_p(vcd_env%matrix_rh(ispin, i)%matrix)
     362            8 :             CALL dbcsr_copy(vcd_env%matrix_rh(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     363              :          END DO
     364              :       END DO
     365              : 
     366              :       ! drpnl for the operator
     367            8 :       DO i = 1, 3
     368            6 :          CALL dbcsr_init_p(vcd_env%matrix_drpnl(i)%matrix)
     369            8 :          CALL dbcsr_copy(vcd_env%matrix_drpnl(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     370              :       END DO
     371              : 
     372              :       ! NVP matrices
     373              :       ! hr matrices
     374            2 :       my_matrix_hr_1d => vcd_env%matrix_hr(1, 1:3)
     375              :       CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
     376              :                              dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
     377            2 :                              direction_Or=.TRUE.)
     378              :       CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
     379            2 :                            direction_Or=.TRUE., rc=[0._dp, 0._dp, 0._dp])
     380            2 :       CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
     381            2 :       CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, [0._dp, 0._dp, 0._dp])
     382              : 
     383            2 :       my_matrix_hr_1d => vcd_env%matrix_rh(1, 1:3)
     384              :       CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
     385              :                              dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
     386            2 :                              direction_Or=.FALSE.)
     387              :       CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
     388            2 :                            direction_Or=.FALSE., rc=[0._dp, 0._dp, 0._dp])
     389            2 :       CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
     390            2 :       CALL build_matrix_r_vhxc(vcd_env%matrix_rh, qs_env, [0._dp, 0._dp, 0._dp])
     391              : 
     392              :       ! commutator terms
     393              :       ! - [V, r]
     394              :       CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
     395            2 :                             particle_set, cell=cell, matrix_rv=vcd_env%hcom)
     396              :       ! <[V, r] * r> and <r * [V, r]>
     397              :       CALL build_com_rpnl_r(vcd_env%matrix_rcomr, qs_kind_set, sab_all, sap_ppnl, &
     398            2 :                             dft_control%qs_control%eps_ppnl, particle_set, cell, .TRUE.)
     399              :       CALL build_com_rpnl_r(vcd_env%matrix_rrcom, qs_kind_set, sab_all, sap_ppnl, &
     400            2 :                             dft_control%qs_control%eps_ppnl, particle_set, cell, .FALSE.)
     401              : 
     402              :       ! lin_mom
     403            2 :       CALL build_lin_mom_matrix(qs_env, vcd_env%dipvel_ao)
     404              : 
     405              :       ! AAT
     406              :       ! The moments are set to zero and then recomputed in the routine.
     407              :       CALL build_local_moments_der_matrix(qs_env, moments_der=vcd_env%moments_der, &
     408            2 :                                           nmoments_der=2, nmoments=0, ref_point=[0._dp, 0._dp, 0._dp])
     409              : 
     410              :       ! PP terms
     411              :       CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
     412              :                             particle_set, matrix_rxrv=vcd_env%matrix_rxrv, ref_point=[0._dp, 0._dp, 0._dp], &
     413            2 :                             cell=cell)
     414              : 
     415              :       CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
     416              :                             particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
     417            2 :                             matrix_r_rxvr=vcd_env%matrix_r_rxvr)
     418              : 
     419              :       CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
     420              :                             particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
     421            2 :                             matrix_rxvr_r=vcd_env%matrix_rxvr_r)
     422              : 
     423              :       ! Done with NVP matrices
     424              : 
     425              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     426            2 :                                         "PRINT%PROGRAM_RUN_INFO")
     427              : 
     428            2 :       CALL timestop(handle)
     429              : 
     430            6 :    END SUBROUTINE vcd_env_init
     431              : 
     432              : ! *****************************************************************************
     433              : !> \brief Deallocate the vcd environment
     434              : !> \param qs_env ...
     435              : !> \param vcd_env ...
     436              : !> \author Edward Ditler
     437              : ! **************************************************************************************************
     438            2 :    SUBROUTINE vcd_env_cleanup(qs_env, vcd_env)
     439              : 
     440              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     441              :       TYPE(vcd_env_type)                                 :: vcd_env
     442              : 
     443              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_env_cleanup'
     444              : 
     445              :       INTEGER                                            :: handle
     446              : 
     447            2 :       CALL timeset(routineN, handle)
     448              : 
     449              :       ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
     450            2 :       CALL dcdr_env_cleanup(qs_env, vcd_env%dcdr_env)
     451              : 
     452            2 :       DEALLOCATE (vcd_env%apt_el_nvpt)
     453            2 :       DEALLOCATE (vcd_env%apt_nuc_nvpt)
     454            2 :       DEALLOCATE (vcd_env%apt_total_nvpt)
     455            2 :       DEALLOCATE (vcd_env%aat_atom_nvpt)
     456            2 :       DEALLOCATE (vcd_env%aat_atom_mfp)
     457              : 
     458            2 :       CALL cp_fm_release(vcd_env%dCV)
     459            2 :       CALL cp_fm_release(vcd_env%dCV_prime)
     460            2 :       CALL cp_fm_release(vcd_env%op_dV)
     461            2 :       CALL cp_fm_release(vcd_env%op_dB)
     462              : 
     463            2 :       CALL cp_fm_release(vcd_env%dCB)
     464            2 :       CALL cp_fm_release(vcd_env%dCB_prime)
     465              : 
     466              :       ! DBCSR matrices
     467              :       ! Probably, the memory requirements could be reduced by quite a bit
     468              :       ! by not storing each term in its own set of matrices.
     469              :       ! On the other hand, the memory bottleneck is usually the numerical
     470              :       ! integration grid.
     471            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der)
     472            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_right)
     473            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_left)
     474            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_difdip2)
     475            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdV)
     476            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdB)
     477            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hxc_dsdv)
     478            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%hcom)
     479            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rcomr)
     480            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rrcom)
     481            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dcom)
     482            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hr)
     483            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rh)
     484            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_drpnl)
     485            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao)
     486            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao_delta)
     487            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxrv)
     488            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_rxvr)
     489            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxvr_r)
     490            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_doublecom)
     491            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp_33)
     492            2 :       CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp2_33)
     493            2 :       CALL timestop(handle)
     494              : 
     495            2 :    END SUBROUTINE vcd_env_cleanup
     496              : 
     497              : ! **************************************************************************************************
     498              : !> \brief Copied from linres_read_restart
     499              : !> \param qs_env ...
     500              : !> \param linres_section ...
     501              : !> \param vec ...
     502              : !> \param lambda ...
     503              : !> \param beta ...
     504              : !> \param tag ...
     505              : !> \author Edward Ditler
     506              : ! **************************************************************************************************
     507           18 :    SUBROUTINE vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
     508              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     509              :       TYPE(section_vals_type), POINTER                   :: linres_section
     510              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: vec
     511              :       INTEGER, INTENT(IN)                                :: lambda, beta
     512              :       CHARACTER(LEN=*)                                   :: tag
     513              : 
     514              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_read_restart'
     515              : 
     516              :       CHARACTER(LEN=default_path_length)                 :: filename
     517              :       CHARACTER(LEN=default_string_length)               :: my_middle
     518              :       INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
     519              :          max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
     520              :       LOGICAL                                            :: file_exists
     521           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     522              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     523              :       TYPE(cp_logger_type), POINTER                      :: logger
     524           18 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     525              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     526              :       TYPE(section_vals_type), POINTER                   :: print_key
     527              : 
     528           18 :       file_exists = .FALSE.
     529              : 
     530           18 :       CALL timeset(routineN, handle)
     531              : 
     532           18 :       NULLIFY (mos, para_env, logger, print_key, vecbuffer)
     533           18 :       logger => cp_get_default_logger()
     534              : 
     535              :       iounit = cp_print_key_unit_nr(logger, linres_section, &
     536           18 :                                     "PRINT%PROGRAM_RUN_INFO", extension=".Log")
     537              : 
     538              :       CALL get_qs_env(qs_env=qs_env, &
     539              :                       para_env=para_env, &
     540           18 :                       mos=mos)
     541              : 
     542           18 :       nspins = SIZE(mos)
     543              : 
     544           18 :       rst_unit = -1
     545           18 :       IF (para_env%is_source()) THEN
     546              :          CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
     547            9 :                                    n_rep_val=n_rep_val)
     548              : 
     549            9 :          CALL XSTRING(tag, ia, ie)
     550              :          my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
     551            9 :                      //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
     552              : 
     553            9 :          IF (n_rep_val > 0) THEN
     554            0 :             CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
     555            0 :             CALL xstring(filename, ia, ie)
     556            0 :             filename = filename(ia:ie)//TRIM(my_middle)//".lr"
     557              :          ELSE
     558              :             ! try to read from the filename that is generated automatically from the printkey
     559            9 :             print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
     560              :             filename = cp_print_key_generate_filename(logger, print_key, &
     561            9 :                                                       extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
     562              :          END IF
     563            9 :          INQUIRE (FILE=filename, exist=file_exists)
     564              :          !
     565              :          ! open file
     566            9 :          IF (file_exists) THEN
     567              :             CALL open_file(file_name=TRIM(filename), &
     568              :                            file_action="READ", &
     569              :                            file_form="UNFORMATTED", &
     570              :                            file_position="REWIND", &
     571              :                            file_status="OLD", &
     572            0 :                            unit_number=rst_unit)
     573              : 
     574            0 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     575            0 :                "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
     576              :          ELSE
     577            9 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     578            9 :                "LINRES| Restart file  <"//TRIM(ADJUSTL(filename))//"> not found"
     579              :          END IF
     580              :       END IF
     581              : 
     582           18 :       CALL para_env%bcast(file_exists)
     583              : 
     584           18 :       IF (file_exists) THEN
     585              : 
     586            0 :          CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
     587            0 :          CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
     588              : 
     589            0 :          ALLOCATE (vecbuffer(nao, max_block))
     590              :          !
     591              :          ! read headers
     592            0 :          IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
     593            0 :          CALL para_env%bcast(iostat)
     594              : 
     595            0 :          CALL para_env%bcast(beta_tmp)
     596            0 :          CALL para_env%bcast(lambda_tmp)
     597            0 :          CALL para_env%bcast(nspins_tmp)
     598            0 :          CALL para_env%bcast(nao_tmp)
     599              : 
     600              :          ! check that the number nao, nmo and nspins are
     601              :          ! the same as in the current mos
     602            0 :          IF (nspins_tmp .NE. nspins) THEN
     603            0 :             CPABORT("nspins not consistent")
     604              :          END IF
     605            0 :          IF (nao_tmp .NE. nao) CPABORT("nao not consistent")
     606              :          ! check that it's the right file
     607              :          ! the same as in the current mos
     608            0 :          IF (lambda_tmp .NE. lambda) CPABORT("lambda not consistent")
     609            0 :          IF (beta_tmp .NE. beta) CPABORT("beta not consistent")
     610              :          !
     611            0 :          DO ispin = 1, nspins
     612            0 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
     613            0 :             CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
     614              :             !
     615            0 :             IF (rst_unit > 0) READ (rst_unit) nmo_tmp
     616            0 :             CALL para_env%bcast(nmo_tmp)
     617            0 :             IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent")
     618              :             !
     619              :             ! read the response
     620            0 :             DO i = 1, nmo, MAX(max_block, 1)
     621            0 :                i_block = MIN(max_block, nmo - i + 1)
     622            0 :                DO j = 1, i_block
     623            0 :                   IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
     624              :                END DO
     625            0 :                CALL para_env%bcast(vecbuffer)
     626            0 :                CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
     627              :             END DO
     628              :          END DO
     629              : 
     630            0 :          IF (iostat /= 0) THEN
     631            0 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     632            0 :                "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
     633              :          END IF
     634              : 
     635            0 :          DEALLOCATE (vecbuffer)
     636              : 
     637              :       END IF
     638              : 
     639           18 :       IF (para_env%is_source()) THEN
     640            9 :          IF (file_exists) CALL close_file(unit_number=rst_unit)
     641              :       END IF
     642              : 
     643           18 :       CALL timestop(handle)
     644              : 
     645           18 :    END SUBROUTINE vcd_read_restart
     646              : 
     647              : ! **************************************************************************************************
     648              : !> \brief Copied from linres_write_restart
     649              : !> \param qs_env ...
     650              : !> \param linres_section ...
     651              : !> \param vec ...
     652              : !> \param lambda ...
     653              : !> \param beta ...
     654              : !> \param tag ...
     655              : !> \author Edward Ditler
     656              : ! **************************************************************************************************
     657           18 :    SUBROUTINE vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
     658              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     659              :       TYPE(section_vals_type), POINTER                   :: linres_section
     660              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: vec
     661              :       INTEGER, INTENT(IN)                                :: lambda, beta
     662              :       CHARACTER(LEN=*)                                   :: tag
     663              : 
     664              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_write_restart'
     665              : 
     666              :       CHARACTER(LEN=default_path_length)                 :: filename
     667              :       CHARACTER(LEN=default_string_length)               :: my_middle, my_pos, my_status
     668              :       INTEGER                                            :: handle, i, i_block, ia, ie, iounit, &
     669              :                                                             ispin, j, max_block, nao, nmo, nspins, &
     670              :                                                             rst_unit
     671           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     672              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     673              :       TYPE(cp_logger_type), POINTER                      :: logger
     674           18 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     675              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     676              :       TYPE(section_vals_type), POINTER                   :: print_key
     677              : 
     678           18 :       NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
     679              : 
     680           18 :       CALL timeset(routineN, handle)
     681              : 
     682           18 :       logger => cp_get_default_logger()
     683              : 
     684           18 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
     685              :                                            used_print_key=print_key), &
     686              :                 cp_p_file)) THEN
     687              : 
     688              :          iounit = cp_print_key_unit_nr(logger, linres_section, &
     689           18 :                                        "PRINT%PROGRAM_RUN_INFO", extension=".Log")
     690              : 
     691              :          CALL get_qs_env(qs_env=qs_env, &
     692              :                          mos=mos, &
     693           18 :                          para_env=para_env)
     694              : 
     695           18 :          nspins = SIZE(mos)
     696              : 
     697           18 :          my_status = "REPLACE"
     698           18 :          my_pos = "REWIND"
     699           18 :          CALL XSTRING(tag, ia, ie)
     700              :          my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
     701           18 :                      //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
     702              :          rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
     703              :                                          extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), &
     704           18 :                                          file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
     705              : 
     706              :          filename = cp_print_key_generate_filename(logger, print_key, &
     707           18 :                                                    extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
     708              : 
     709           18 :          IF (iounit > 0) THEN
     710              :             WRITE (UNIT=iounit, FMT="(T2,A)") &
     711            9 :                "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">"
     712              :          END IF
     713              : 
     714              :          !
     715              :          ! write data to file
     716              :          ! use the scalapack block size as a default for buffering columns
     717           18 :          CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
     718           18 :          CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
     719           72 :          ALLOCATE (vecbuffer(nao, max_block))
     720              : 
     721           18 :          IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
     722              : 
     723           36 :          DO ispin = 1, nspins
     724           18 :             CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
     725              : 
     726           18 :             IF (rst_unit > 0) WRITE (rst_unit) nmo
     727              : 
     728           72 :             DO i = 1, nmo, MAX(max_block, 1)
     729           18 :                i_block = MIN(max_block, nmo - i + 1)
     730           18 :                CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
     731              :                ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
     732              :                ! to old ones, and in cases where max_block is different between runs, as might happen during
     733              :                ! restarts with a different number of CPUs
     734          108 :                DO j = 1, i_block
     735           90 :                   IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
     736              :                END DO
     737              :             END DO
     738              :          END DO
     739              : 
     740           18 :          DEALLOCATE (vecbuffer)
     741              : 
     742              :          CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
     743           36 :                                            "PRINT%RESTART")
     744              :       END IF
     745              : 
     746           18 :       CALL timestop(handle)
     747              : 
     748           18 :    END SUBROUTINE vcd_write_restart
     749              : 
     750              : ! **************************************************************************************************
     751              : !> \brief Print the APTs, AATs, and sum rules
     752              : !> \param vcd_env ...
     753              : !> \param qs_env ...
     754              : !> \author Edward Ditler
     755              : ! **************************************************************************************************
     756            2 :    SUBROUTINE vcd_print(vcd_env, qs_env)
     757              :       TYPE(vcd_env_type)                                 :: vcd_env
     758              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     759              : 
     760              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'vcd_print'
     761              : 
     762              :       CHARACTER(LEN=default_string_length)               :: description
     763              :       INTEGER                                            :: alpha, beta, delta, gamma, handle, i, l, &
     764              :                                                             lambda, natom, nsubset, output_unit
     765              :       REAL(dp)                                           :: mean, standard_deviation, &
     766              :                                                             standard_deviation_sum
     767            2 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: apt_el_dcdr, apt_el_nvpt, apt_nuc_dcdr, &
     768            2 :                                                             apt_nuc_nvpt, apt_total_dcdr, &
     769            2 :                                                             apt_total_nvpt
     770            2 :       REAL(dp), DIMENSION(:, :, :, :), POINTER           :: apt_center_dcdr, apt_subset_dcdr
     771              :       REAL(kind=dp), DIMENSION(3, 3)                     :: sum_rule_0, sum_rule_0_second, &
     772              :                                                             sum_rule_1, sum_rule_2, &
     773              :                                                             sum_rule_2_second, sum_rule_3_mfp, &
     774              :                                                             sum_rule_3_second
     775              :       TYPE(cp_logger_type), POINTER                      :: logger
     776              :       TYPE(cp_result_type), POINTER                      :: results
     777            2 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     778            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     779              :       TYPE(section_vals_type), POINTER                   :: vcd_section
     780              : 
     781            2 :       CALL timeset(routineN, handle)
     782              : 
     783            2 :       NULLIFY (logger)
     784              : 
     785            2 :       logger => cp_get_default_logger()
     786            2 :       output_unit = cp_logger_get_default_io_unit(logger)
     787              : 
     788            2 :       vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
     789              : 
     790            2 :       NULLIFY (particle_set)
     791            2 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
     792            2 :       natom = SIZE(particle_set)
     793            2 :       nsubset = SIZE(molecule_set)
     794              : 
     795            2 :       apt_el_dcdr => vcd_env%dcdr_env%apt_el_dcdr
     796            2 :       apt_nuc_dcdr => vcd_env%dcdr_env%apt_nuc_dcdr
     797            2 :       apt_total_dcdr => vcd_env%dcdr_env%apt_total_dcdr
     798            2 :       apt_subset_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_subset
     799            2 :       apt_center_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_center
     800              : 
     801            2 :       apt_el_nvpt => vcd_env%apt_el_nvpt
     802            2 :       apt_nuc_nvpt => vcd_env%apt_nuc_nvpt
     803            2 :       apt_total_nvpt => vcd_env%apt_total_nvpt
     804              : 
     805            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
     806            1 :          'APT | Write the final APT matrix per atom (Position perturbation)'
     807            8 :       DO l = 1, natom
     808            6 :          IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
     809            3 :             'APT | Atom', l, ' - GAPT ', &
     810              :             (apt_total_dcdr(1, 1, l) &
     811              :              + apt_total_dcdr(2, 2, l) &
     812            6 :              + apt_total_dcdr(3, 3, l))/3._dp
     813           26 :          DO i = 1, 3
     814           24 :             IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
     815              :          END DO
     816              :       END DO
     817              : 
     818            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
     819            1 :          'NVP | Write the final APT matrix per atom (Velocity perturbation)'
     820            8 :       DO l = 1, natom
     821            6 :          IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
     822            3 :             'NVP | Atom', l, ' - GAPT ', &
     823              :             (apt_total_nvpt(1, 1, l) &
     824              :              + apt_total_nvpt(2, 2, l) &
     825            6 :              + apt_total_nvpt(3, 3, l))/3._dp
     826           26 :          DO i = 1, 3
     827           18 :             IF (vcd_env%output_unit > 0) &
     828              :                WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
     829           15 :                "NVP | ", apt_total_nvpt(i, :, l)
     830              :          END DO
     831              :       END DO
     832              : 
     833            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
     834            1 :          'NVP | Write the final AAT matrix per atom (Velocity perturbation)'
     835            8 :       DO l = 1, natom
     836            6 :          IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
     837            3 :             'NVP | Atom', l
     838           26 :          DO i = 1, 3
     839           18 :             IF (vcd_env%output_unit > 0) &
     840              :                WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
     841           15 :                "NVP | ", vcd_env%aat_atom_nvpt(i, :, l)
     842              :          END DO
     843              :       END DO
     844              : 
     845            2 :       IF (vcd_env%do_mfp) THEN
     846            0 :          IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
     847            0 :             'MFP | Write the final AAT matrix per atom (Magnetic Field perturbation)'
     848            0 :          DO l = 1, natom
     849            0 :             IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
     850            0 :                'MFP | Atom', l
     851            0 :             DO i = 1, 3
     852            0 :                IF (vcd_env%output_unit > 0) &
     853              :                   WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
     854            0 :                   "MFP | ", vcd_env%aat_atom_mfp(i, :, l)
     855              :             END DO
     856              :          END DO
     857              :       END IF
     858              : 
     859              :       ! Get the dipole
     860            2 :       CALL get_qs_env(qs_env, results=results)
     861            2 :       description = "[DIPOLE]"
     862            2 :       CALL get_results(results=results, description=description, values=vcd_env%dcdr_env%dipole_pos(1:3))
     863              : 
     864              :       ! Sum rules [for all alpha, beta]
     865            2 :       sum_rule_0 = 0._dp
     866            2 :       sum_rule_1 = 0._dp
     867            2 :       sum_rule_2 = 0._dp
     868            2 :       sum_rule_0_second = 0._dp
     869            2 :       sum_rule_2_second = 0._dp
     870            2 :       sum_rule_3_second = 0._dp
     871            2 :       sum_rule_3_mfp = 0._dp
     872            2 :       standard_deviation = 0._dp
     873            2 :       standard_deviation_sum = 0._dp
     874              : 
     875            8 :       DO alpha = 1, 3
     876           26 :          DO beta = 1, 3
     877              :             ! 0: sum_lambda apt(alpha, beta, lambda)
     878           72 :             DO lambda = 1, natom
     879              :                sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
     880           54 :                                          + apt_total_dcdr(alpha, beta, lambda)
     881              :                sum_rule_0_second(alpha, beta) = sum_rule_0_second(alpha, beta) &
     882           72 :                                                 + apt_total_nvpt(alpha, beta, lambda)
     883              :             END DO
     884              : 
     885              :             ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
     886           72 :             DO gamma = 1, 3
     887              :                sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
     888           72 :                                          + Levi_Civita(alpha, beta, gamma)*vcd_env%dcdr_env%dipole_pos(gamma)
     889              :             END DO
     890              : 
     891              :             ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
     892           72 :             DO lambda = 1, natom
     893          234 :                DO gamma = 1, 3
     894          702 :                   DO delta = 1, 3
     895              :                      sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
     896              :                                                + Levi_Civita(beta, gamma, delta) &
     897              :                                                *particle_set(lambda)%r(gamma) &
     898          486 :                                                *apt_total_dcdr(delta, alpha, lambda)
     899              :                      sum_rule_2_second(alpha, beta) = sum_rule_2_second(alpha, beta) &
     900              :                                                       + Levi_Civita(beta, gamma, delta) &
     901              :                                                       *particle_set(lambda)%r(gamma) &
     902          648 :                                                       *apt_total_nvpt(delta, alpha, lambda)
     903              :                   END DO
     904              :                END DO
     905              :             END DO
     906              : 
     907              :             ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
     908           72 :             DO lambda = 1, natom
     909              :                sum_rule_3_second(alpha, beta) = sum_rule_3_second(alpha, beta) &
     910           72 :                                                 + vcd_env%aat_atom_nvpt(alpha, beta, lambda)
     911              :                ! + 2._dp*c_light_au*vcd_env%aat_atom_nvpt(alpha, beta, lambda)
     912              :             END DO
     913              : 
     914           24 :             IF (vcd_env%do_mfp) THEN
     915              :                ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
     916            0 :                DO lambda = 1, natom
     917              :                   sum_rule_3_mfp(alpha, beta) = sum_rule_3_mfp(alpha, beta) &
     918            0 :                                                 + vcd_env%aat_atom_mfp(alpha, beta, lambda)
     919              :                END DO
     920              :             END IF
     921              : 
     922              :          END DO ! beta
     923              :       END DO   ! alpha
     924              : 
     925            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") "APT | Position perturbation sum rules"
     926            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") &
     927            1 :          "APT |", " Total APT", "Dipole", "R * APT", "AAT"
     928              :       standard_deviation_sum = 0._dp
     929            8 :       DO alpha = 1, 3
     930           26 :          DO beta = 1, 3
     931           18 :             mean = (sum_rule_1(alpha, beta) + sum_rule_2(alpha, beta) + sum_rule_3_mfp(alpha, beta))/3
     932              :             standard_deviation = &
     933              :                SQRT((sum_rule_1(alpha, beta)**2 + sum_rule_2(alpha, beta)**2 + sum_rule_3_mfp(alpha, beta)**2)/3 &
     934           18 :                     - mean**2)
     935           18 :             standard_deviation_sum = standard_deviation_sum + standard_deviation
     936              : 
     937           18 :             IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
     938              :                                                 "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
     939            9 :                "APT | ", &
     940            9 :                alpha, beta, &
     941            9 :                sum_rule_0(alpha, beta), &
     942            9 :                sum_rule_1(alpha, beta), &
     943            9 :                sum_rule_2(alpha, beta), &
     944            9 :                sum_rule_3_mfp(alpha, beta), &
     945           24 :                standard_deviation
     946              :          END DO
     947              :       END DO
     948            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
     949              : 
     950            2 :       IF (vcd_env%output_unit > 0) THEN
     951            1 :          WRITE (vcd_env%output_unit, "(A)") "NVP | Velocity perturbation sum rules"
     952            1 :          WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") "NVP |", " Total APT", "Dipole", "R * APT", "AAT"
     953              :       END IF
     954              : 
     955            2 :       standard_deviation_sum = 0._dp
     956            8 :       DO alpha = 1, 3
     957           26 :          DO beta = 1, 3
     958           18 :             mean = (sum_rule_1(alpha, beta) + sum_rule_2_second(alpha, beta) + sum_rule_3_second(alpha, beta))/3
     959              :             standard_deviation = &
     960              :                SQRT((sum_rule_1(alpha, beta)**2 + sum_rule_2_second(alpha, beta)**2 + sum_rule_3_second(alpha, beta)**2)/3 &
     961           18 :                     - mean**2)
     962           18 :             standard_deviation_sum = standard_deviation_sum + standard_deviation
     963           18 :             IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
     964              :                                                 "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
     965            9 :                "NVP | ", &
     966            9 :                alpha, &
     967            9 :                beta, &
     968            9 :                sum_rule_0_second(alpha, beta), &
     969            9 :                sum_rule_1(alpha, beta), &
     970            9 :                sum_rule_2_second(alpha, beta), &
     971            9 :                sum_rule_3_second(alpha, beta), &
     972           24 :                standard_deviation
     973              :          END DO
     974              :       END DO
     975            2 :       IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
     976              : 
     977            2 :       CALL timestop(handle)
     978            2 :    END SUBROUTINE vcd_print
     979              : 
     980              : END MODULE qs_vcd_utils
        

Generated by: LCOV version 2.0-1