LCOV - code coverage report
Current view: top level - src - qs_efield_berry.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 610 644 94.7 %
Date: 2024-03-28 07:31:50 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculates the energy contribution and the mo_derivative of
      10             : !>        a static periodic electric field
      11             : !> \par History
      12             : !>      none
      13             : !> \author fschiff (06.2010)
      14             : ! **************************************************************************************************
      15             : MODULE qs_efield_berry
      16             :    USE ai_moments,                      ONLY: cossin
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18             :                                               get_atomic_kind,&
      19             :                                               get_atomic_kind_set
      20             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      21             :                                               gto_basis_set_type
      22             :    USE block_p_types,                   ONLY: block_p_type
      23             :    USE cell_types,                      ONLY: cell_type,&
      24             :                                               pbc
      25             :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_scale_and_add_fm,&
      26             :                                               cp_cfm_solve
      27             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      28             :                                               cp_cfm_release,&
      29             :                                               cp_cfm_set_all,&
      30             :                                               cp_cfm_type
      31             :    USE cp_control_types,                ONLY: dft_control_type
      32             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      33             :                                               copy_fm_to_dbcsr,&
      34             :                                               cp_dbcsr_plus_fm_fm_t,&
      35             :                                               cp_dbcsr_sm_fm_multiply,&
      36             :                                               dbcsr_deallocate_matrix_set
      37             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      38             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      39             :                                               cp_fm_struct_release,&
      40             :                                               cp_fm_struct_type
      41             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      42             :                                               cp_fm_release,&
      43             :                                               cp_fm_set_all,&
      44             :                                               cp_fm_type
      45             :    USE dbcsr_api,                       ONLY: dbcsr_copy,&
      46             :                                               dbcsr_get_block_p,&
      47             :                                               dbcsr_p_type,&
      48             :                                               dbcsr_set,&
      49             :                                               dbcsr_type
      50             :    USE kinds,                           ONLY: dp
      51             :    USE mathconstants,                   ONLY: gaussi,&
      52             :                                               pi,&
      53             :                                               twopi,&
      54             :                                               z_one,&
      55             :                                               z_zero
      56             :    USE message_passing,                 ONLY: mp_para_env_type
      57             :    USE orbital_pointers,                ONLY: ncoset
      58             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      59             :    USE particle_types,                  ONLY: particle_type
      60             :    USE qs_energy_types,                 ONLY: qs_energy_type
      61             :    USE qs_environment_types,            ONLY: get_qs_env,&
      62             :                                               qs_environment_type,&
      63             :                                               set_qs_env
      64             :    USE qs_force_types,                  ONLY: qs_force_type
      65             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      66             :                                               get_qs_kind_set,&
      67             :                                               qs_kind_type
      68             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      69             :                                               mo_set_type
      70             :    USE qs_moments,                      ONLY: build_berry_moment_matrix
      71             :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      72             :                                               neighbor_list_iterate,&
      73             :                                               neighbor_list_iterator_create,&
      74             :                                               neighbor_list_iterator_p_type,&
      75             :                                               neighbor_list_iterator_release,&
      76             :                                               neighbor_list_set_p_type
      77             :    USE qs_period_efield_types,          ONLY: efield_berry_type,&
      78             :                                               init_efield_matrices,&
      79             :                                               set_efield_matrices
      80             :    USE virial_methods,                  ONLY: virial_pair_force
      81             :    USE virial_types,                    ONLY: virial_type
      82             : #include "./base/base_uses.f90"
      83             : 
      84             :    IMPLICIT NONE
      85             : 
      86             :    PRIVATE
      87             : 
      88             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_efield_berry'
      89             : 
      90             :    ! *** Public subroutines ***
      91             : 
      92             :    PUBLIC :: qs_efield_berry_phase
      93             : 
      94             : ! **************************************************************************************************
      95             : 
      96             : CONTAINS
      97             : 
      98             : ! **************************************************************************************************
      99             : 
     100             : ! **************************************************************************************************
     101             : !> \brief ...
     102             : !> \param qs_env ...
     103             : !> \param just_energy ...
     104             : !> \param calculate_forces ...
     105             : ! **************************************************************************************************
     106       96439 :    SUBROUTINE qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
     107             : 
     108             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     109             :       LOGICAL, INTENT(IN)                                :: just_energy, calculate_forces
     110             : 
     111             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_berry_phase'
     112             : 
     113             :       INTEGER                                            :: handle
     114             :       LOGICAL                                            :: s_mstruct_changed
     115             :       TYPE(dft_control_type), POINTER                    :: dft_control
     116             : 
     117       96439 :       CALL timeset(routineN, handle)
     118             : 
     119       96439 :       NULLIFY (dft_control)
     120             :       CALL get_qs_env(qs_env, s_mstruct_changed=s_mstruct_changed, &
     121       96439 :                       dft_control=dft_control)
     122             : 
     123       96439 :       IF (dft_control%apply_period_efield) THEN
     124        2058 :          IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env)
     125        2058 :          IF (dft_control%period_efield%displacement_field) THEN
     126         898 :             CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
     127             :          ELSE
     128        1160 :             CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces)
     129             :          END IF
     130             :       END IF
     131             : 
     132       96439 :       CALL timestop(handle)
     133             : 
     134       96439 :    END SUBROUTINE qs_efield_berry_phase
     135             : 
     136             : ! **************************************************************************************************
     137             : !> \brief ...
     138             : !> \param qs_env ...
     139             : ! **************************************************************************************************
     140         144 :    SUBROUTINE qs_efield_integrals(qs_env)
     141             : 
     142             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     143             : 
     144             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_integrals'
     145             : 
     146             :       INTEGER                                            :: handle, i
     147             :       REAL(dp), DIMENSION(3)                             :: kvec
     148             :       TYPE(cell_type), POINTER                           :: cell
     149         144 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: cosmat, matrix_s, sinmat
     150             :       TYPE(dft_control_type), POINTER                    :: dft_control
     151             :       TYPE(efield_berry_type), POINTER                   :: efield
     152             : 
     153         144 :       CALL timeset(routineN, handle)
     154         144 :       CPASSERT(ASSOCIATED(qs_env))
     155             : 
     156         144 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     157         144 :       NULLIFY (matrix_s)
     158         144 :       CALL get_qs_env(qs_env=qs_env, efield=efield, cell=cell, matrix_s=matrix_s)
     159         144 :       CALL init_efield_matrices(efield)
     160        1008 :       ALLOCATE (cosmat(3), sinmat(3))
     161         576 :       DO i = 1, 3
     162         432 :          ALLOCATE (cosmat(i)%matrix, sinmat(i)%matrix)
     163             : 
     164         432 :          CALL dbcsr_copy(cosmat(i)%matrix, matrix_s(1)%matrix, 'COS MAT')
     165         432 :          CALL dbcsr_copy(sinmat(i)%matrix, matrix_s(1)%matrix, 'SIN MAT')
     166         432 :          CALL dbcsr_set(cosmat(i)%matrix, 0.0_dp)
     167         432 :          CALL dbcsr_set(sinmat(i)%matrix, 0.0_dp)
     168             : 
     169        1728 :          kvec(:) = twopi*cell%h_inv(i, :)
     170         576 :          CALL build_berry_moment_matrix(qs_env, cosmat(i)%matrix, sinmat(i)%matrix, kvec)
     171             :       END DO
     172         144 :       CALL set_efield_matrices(efield=efield, cosmat=cosmat, sinmat=sinmat)
     173         144 :       CALL set_qs_env(qs_env=qs_env, efield=efield)
     174         144 :       CALL timestop(handle)
     175             : 
     176         144 :    END SUBROUTINE qs_efield_integrals
     177             : 
     178             : ! **************************************************************************************************
     179             : !> \brief ...
     180             : !> \param qs_env ...
     181             : !> \param just_energy ...
     182             : !> \param calculate_forces ...
     183             : ! **************************************************************************************************
     184        1160 :    SUBROUTINE qs_efield_derivatives(qs_env, just_energy, calculate_forces)
     185             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     186             :       LOGICAL, INTENT(IN)                                :: just_energy, calculate_forces
     187             : 
     188             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_derivatives'
     189             : 
     190             :       COMPLEX(dp)                                        :: zdet, zdeta, zi(3)
     191             :       INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, j, &
     192             :          jatom, jkind, jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, &
     193             :          nseta, nsetb, sgfa, sgfb
     194        1160 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     195        1160 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     196        1160 :                                                             npgfb, nsgfa, nsgfb
     197        1160 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     198             :       LOGICAL                                            :: found, uniform, use_virial
     199             :       REAL(dp)                                           :: charge, ci(3), cqi(3), dab, dd, &
     200             :                                                             ener_field, f0, fab, fieldpol(3), &
     201             :                                                             focc, fpolvec(3), hmat(3, 3), occ, &
     202             :                                                             qi(3), ti(3)
     203             :       REAL(dp), DIMENSION(3)                             :: forcea, forceb, kvec, ra, rab, rb, ria
     204        2320 :       REAL(dp), DIMENSION(:, :), POINTER                 :: cosab, iblock, rblock, sinab, work
     205        2320 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: dcosab, dsinab
     206        1160 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     207        1160 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
     208        1160 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     209       20880 :       TYPE(block_p_type), DIMENSION(3, 2)                :: dcost, dsint
     210             :       TYPE(cell_type), POINTER                           :: cell
     211        1160 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: eigrmat, inv_mat
     212             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     213        1160 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: mo_coeff_tmp, mo_derivs_tmp
     214        1160 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: inv_work, op_fm_set, opvec
     215             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     216        1160 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, mo_derivs
     217        1160 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: tempmat
     218             :       TYPE(dbcsr_type), POINTER                          :: cosmat, mo_coeff_b, sinmat
     219             :       TYPE(dft_control_type), POINTER                    :: dft_control
     220             :       TYPE(efield_berry_type), POINTER                   :: efield
     221        1160 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     222             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     223        1160 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     224             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     225             :       TYPE(neighbor_list_iterator_p_type), &
     226        1160 :          DIMENSION(:), POINTER                           :: nl_iterator
     227             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     228        1160 :          POINTER                                         :: sab_orb
     229        1160 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     230             :       TYPE(qs_energy_type), POINTER                      :: energy
     231        1160 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     232        1160 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     233             :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     234             :       TYPE(virial_type), POINTER                         :: virial
     235             : 
     236        1160 :       CALL timeset(routineN, handle)
     237             : 
     238        1160 :       NULLIFY (dft_control, cell, particle_set)
     239             :       CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
     240        1160 :                       particle_set=particle_set, virial=virial)
     241        1160 :       NULLIFY (qs_kind_set, efield, para_env, sab_orb)
     242             :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
     243        1160 :                       efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
     244             : 
     245             :       ! calculate stress only if forces requested also
     246        1160 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     247           0 :       use_virial = use_virial .AND. calculate_forces
     248             :       ! disable stress calculation
     249             :       IF (use_virial) THEN
     250           0 :          CPABORT("Stress tensor for periodic E-field not implemented")
     251             :       END IF
     252             : 
     253        4640 :       fieldpol = dft_control%period_efield%polarisation
     254        8120 :       fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
     255        4640 :       fieldpol = -fieldpol*dft_control%period_efield%strength
     256       15080 :       hmat = cell%hmat(:, :)/twopi
     257        4640 :       DO idir = 1, 3
     258        4640 :          fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
     259             :       END DO
     260             : 
     261             :       ! nuclear contribution
     262        1160 :       natom = SIZE(particle_set)
     263        1160 :       IF (calculate_forces) THEN
     264          28 :          CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
     265          28 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
     266             :       END IF
     267        4640 :       zi(:) = CMPLX(1._dp, 0._dp, dp)
     268        3960 :       DO ia = 1, natom
     269        2800 :          CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
     270        2800 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
     271       11200 :          ria = particle_set(ia)%r
     272       11200 :          ria = pbc(ria, cell)
     273       11200 :          DO idir = 1, 3
     274       33600 :             kvec(:) = twopi*cell%h_inv(idir, :)
     275       33600 :             dd = SUM(kvec(:)*ria(:))
     276        8400 :             zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
     277       11200 :             zi(idir) = zi(idir)*zdeta
     278             :          END DO
     279        2800 :          IF (calculate_forces) THEN
     280          74 :             IF (para_env%mepos == 0) THEN
     281          37 :                iatom = atom_of_kind(ia)
     282         148 :                forcea(:) = fieldpol(:)*charge
     283         148 :                force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + forcea(:)
     284             :             END IF
     285             :          END IF
     286        6760 :          IF (use_virial) THEN
     287           0 :             IF (para_env%mepos == 0) &
     288           0 :                CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria)
     289             :          END IF
     290             :       END DO
     291        4640 :       qi = AIMAG(LOG(zi))
     292             : 
     293             :       ! check uniform occupation
     294        1160 :       NULLIFY (mos)
     295        1160 :       CALL get_qs_env(qs_env=qs_env, mos=mos)
     296        2380 :       DO ispin = 1, dft_control%nspins
     297        1220 :          CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
     298        2380 :          IF (.NOT. uniform) THEN
     299           0 :             CPABORT("Berry phase moments for non uniform MOs' occupation numbers not implemented")
     300             :          END IF
     301             :       END DO
     302             : 
     303        1160 :       NULLIFY (mo_derivs)
     304        1160 :       CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
     305             :       ! initialize all work matrices needed
     306        7140 :       ALLOCATE (op_fm_set(2, dft_control%nspins))
     307        7140 :       ALLOCATE (opvec(2, dft_control%nspins))
     308        4700 :       ALLOCATE (eigrmat(dft_control%nspins))
     309        4700 :       ALLOCATE (inv_mat(dft_control%nspins))
     310        7140 :       ALLOCATE (inv_work(2, dft_control%nspins))
     311        4700 :       ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs)))
     312        4700 :       ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
     313             : 
     314             :       ! Allocate temp matrices for the wavefunction derivatives
     315        2380 :       DO ispin = 1, dft_control%nspins
     316        1220 :          NULLIFY (tmp_fm_struct, mo_coeff)
     317        1220 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
     318             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
     319        1220 :                                   ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
     320        1220 :          CALL cp_fm_create(mo_derivs_tmp(ispin), mo_coeff%matrix_struct)
     321        1220 :          CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
     322        1220 :          CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_tmp(ispin))
     323        3660 :          DO i = 1, SIZE(op_fm_set, 1)
     324        2440 :             CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
     325        2440 :             CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
     326        3660 :             CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
     327             :          END DO
     328        1220 :          CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
     329        1220 :          CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
     330        3600 :          CALL cp_fm_struct_release(tmp_fm_struct)
     331             :       END DO
     332             :       ! temp matrices for force calculation
     333        1160 :       IF (calculate_forces) THEN
     334          28 :          NULLIFY (matrix_s)
     335          28 :          CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
     336         192 :          ALLOCATE (tempmat(2, dft_control%nspins))
     337          64 :          DO ispin = 1, dft_control%nspins
     338          36 :             ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
     339          36 :             CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
     340          36 :             CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
     341          36 :             CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
     342          64 :             CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
     343             :          END DO
     344             :          ! integration
     345          28 :          CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
     346         280 :          ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
     347         252 :          ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
     348          28 :          lsab = MAX(ldab, lsab)
     349         140 :          DO i = 1, 3
     350         588 :             ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
     351         616 :             ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
     352             :          END DO
     353             :       END IF
     354             : 
     355             :       !Start the MO derivative calculation
     356             :       !loop over all cell vectors
     357        4640 :       DO idir = 1, 3
     358        3480 :          ci(idir) = 0.0_dp
     359             :          zi(idir) = z_zero
     360        4640 :          IF (ABS(fpolvec(idir)) .GT. 1.0E-12_dp) THEN
     361        2442 :             cosmat => efield%cosmat(idir)%matrix
     362        2442 :             sinmat => efield%sinmat(idir)%matrix
     363             :             !evaluate the expression needed for the derivative (S_berry * C  and [C^T S_berry C]^-1)
     364             :             !first step S_berry * C  and C^T S_berry C
     365        4944 :             DO ispin = 1, dft_control%nspins ! spin
     366        2502 :                IF (mos(ispin)%use_mo_coeff_b) THEN
     367        2502 :                   CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
     368        2502 :                   CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
     369             :                ELSE
     370           0 :                   CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
     371           0 :                   mo_coeff_tmp(ispin) = mo_coeff
     372             :                END IF
     373        2502 :                CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
     374             :                CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
     375        2502 :                                   op_fm_set(1, ispin))
     376        2502 :                CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
     377             :                CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
     378        4944 :                                   op_fm_set(2, ispin))
     379             :             END DO
     380             :             !second step invert C^T S_berry C
     381        2442 :             zdet = z_one
     382        4944 :             DO ispin = 1, dft_control%nspins
     383        2502 :                CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
     384        2502 :                CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
     385        2502 :                CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
     386        2502 :                CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
     387        4944 :                zdet = zdet*zdeta
     388             :             END DO
     389        2442 :             zi(idir) = zdet**occ
     390        2442 :             ci(idir) = AIMAG(LOG(zdet**occ))
     391             : 
     392        2442 :             IF (.NOT. just_energy) THEN
     393             :                !compute the orbital derivative
     394        2390 :                focc = fpolvec(idir)
     395        4818 :                DO ispin = 1, dft_control%nspins
     396       31564 :                   inv_work(1, ispin)%local_data(:, :) = REAL(inv_mat(ispin)%local_data(:, :), dp)
     397       31564 :                   inv_work(2, ispin)%local_data(:, :) = AIMAG(inv_mat(ispin)%local_data(:, :))
     398        2428 :                   CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
     399             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
     400        2428 :                                      1.0_dp, mo_derivs_tmp(ispin))
     401             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
     402        7246 :                                      1.0_dp, mo_derivs_tmp(ispin))
     403             :                END DO
     404             :             END IF
     405             : 
     406             :             !compute nuclear forces
     407        2442 :             IF (calculate_forces) THEN
     408          46 :                nkind = SIZE(qs_kind_set)
     409          46 :                natom = SIZE(particle_set)
     410         184 :                kvec(:) = twopi*cell%h_inv(idir, :)
     411             : 
     412             :                ! calculate: C [C^T S_berry C]^(-1) C^T
     413             :                ! Store this matrix in DBCSR form (only S overlap blocks)
     414         100 :                DO ispin = 1, dft_control%nspins
     415          54 :                   CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
     416          54 :                   CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
     417          54 :                   CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
     418             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
     419          54 :                                      opvec(1, ispin))
     420             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
     421          54 :                                      opvec(2, ispin))
     422             :                   CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
     423          54 :                                              matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
     424             :                   CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
     425         154 :                                              matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
     426             :                END DO
     427             : 
     428             :                ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
     429         230 :                ALLOCATE (basis_set_list(nkind))
     430         138 :                DO ikind = 1, nkind
     431          92 :                   qs_kind => qs_kind_set(ikind)
     432          92 :                   CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
     433         138 :                   IF (ASSOCIATED(basis_set_a)) THEN
     434          92 :                      basis_set_list(ikind)%gto_basis_set => basis_set_a
     435             :                   ELSE
     436           0 :                      NULLIFY (basis_set_list(ikind)%gto_basis_set)
     437             :                   END IF
     438             :                END DO
     439             :                !
     440          46 :                CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     441        1472 :                DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     442             :                   CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     443        1426 :                                          iatom=iatom, jatom=jatom, r=rab)
     444        1426 :                   basis_set_a => basis_set_list(ikind)%gto_basis_set
     445        1426 :                   IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     446        1426 :                   basis_set_b => basis_set_list(jkind)%gto_basis_set
     447        1426 :                   IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     448             :                   ! basis ikind
     449        1426 :                   first_sgfa => basis_set_a%first_sgf
     450        1426 :                   la_max => basis_set_a%lmax
     451        1426 :                   la_min => basis_set_a%lmin
     452        1426 :                   npgfa => basis_set_a%npgf
     453        1426 :                   nseta = basis_set_a%nset
     454        1426 :                   nsgfa => basis_set_a%nsgf_set
     455        1426 :                   rpgfa => basis_set_a%pgf_radius
     456        1426 :                   set_radius_a => basis_set_a%set_radius
     457        1426 :                   sphi_a => basis_set_a%sphi
     458        1426 :                   zeta => basis_set_a%zet
     459             :                   ! basis jkind
     460        1426 :                   first_sgfb => basis_set_b%first_sgf
     461        1426 :                   lb_max => basis_set_b%lmax
     462        1426 :                   lb_min => basis_set_b%lmin
     463        1426 :                   npgfb => basis_set_b%npgf
     464        1426 :                   nsetb = basis_set_b%nset
     465        1426 :                   nsgfb => basis_set_b%nsgf_set
     466        1426 :                   rpgfb => basis_set_b%pgf_radius
     467        1426 :                   set_radius_b => basis_set_b%set_radius
     468        1426 :                   sphi_b => basis_set_b%sphi
     469        1426 :                   zetb => basis_set_b%zet
     470             : 
     471        1426 :                   atom_a = atom_of_kind(iatom)
     472        1426 :                   atom_b = atom_of_kind(jatom)
     473             : 
     474        1426 :                   ldsa = SIZE(sphi_a, 1)
     475        1426 :                   ldsb = SIZE(sphi_b, 1)
     476        1426 :                   ra(:) = pbc(particle_set(iatom)%r(:), cell)
     477        5704 :                   rb(:) = ra + rab
     478        1426 :                   dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
     479             : 
     480        1426 :                   IF (iatom <= jatom) THEN
     481         960 :                      irow = iatom
     482         960 :                      icol = jatom
     483             :                   ELSE
     484         466 :                      irow = jatom
     485         466 :                      icol = iatom
     486             :                   END IF
     487             : 
     488        1426 :                   IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
     489             :                      fab = 1.0_dp*occ
     490             :                   ELSE
     491        1372 :                      fab = 2.0_dp*occ
     492             :                   END IF
     493             : 
     494        5704 :                   DO i = 1, 3
     495     1167894 :                      dcost(i, 1)%block = 0.0_dp
     496     1167894 :                      dsint(i, 1)%block = 0.0_dp
     497     1167894 :                      dcost(i, 2)%block = 0.0_dp
     498     1169320 :                      dsint(i, 2)%block = 0.0_dp
     499             :                   END DO
     500             : 
     501        4278 :                   DO iset = 1, nseta
     502        2852 :                      ncoa = npgfa(iset)*ncoset(la_max(iset))
     503        2852 :                      sgfa = first_sgfa(1, iset)
     504        9982 :                      DO jset = 1, nsetb
     505        5704 :                         IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     506        2744 :                         ncob = npgfb(jset)*ncoset(lb_max(jset))
     507        2744 :                         sgfb = first_sgfb(1, jset)
     508             :                         ! Calculate the primitive integrals (da|b)
     509             :                         CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     510             :                                     lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     511        2744 :                                     ra, rb, kvec, cosab, sinab, dcosab, dsinab)
     512       10976 :                         DO i = 1, 3
     513             :                            CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
     514             :                                              ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
     515             :                                              ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
     516       10976 :                                              dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
     517             :                         END DO
     518             :                         ! Calculate the primitive integrals (a|db)
     519             :                         CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     520             :                                     la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     521        2744 :                                     rb, ra, kvec, cosab, sinab, dcosab, dsinab)
     522       13828 :                         DO i = 1, 3
     523     1264188 :                            dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i))
     524     1264188 :                            dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i))
     525             :                            CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
     526             :                                              ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
     527             :                                              ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
     528       13936 :                                              dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
     529             :                         END DO
     530             :                      END DO
     531             :                   END DO
     532        1426 :                   forcea = 0.0_dp
     533        1426 :                   forceb = 0.0_dp
     534        3252 :                   DO ispin = 1, dft_control%nspins
     535        1826 :                      NULLIFY (rblock, iblock)
     536             :                      CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
     537        1826 :                                             row=irow, col=icol, BLOCK=rblock, found=found)
     538        1826 :                      CPASSERT(found)
     539             :                      CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
     540        1826 :                                             row=irow, col=icol, BLOCK=iblock, found=found)
     541        1826 :                      CPASSERT(found)
     542        1826 :                      n1 = SIZE(rblock, 1)
     543        1826 :                      n2 = SIZE(rblock, 2)
     544        1826 :                      CPASSERT(SIZE(iblock, 1) == n1)
     545        1826 :                      CPASSERT(SIZE(iblock, 2) == n2)
     546        1826 :                      CPASSERT(lsab >= n1)
     547        1826 :                      CPASSERT(lsab >= n2)
     548        6904 :                      IF (iatom <= jatom) THEN
     549        4976 :                         DO i = 1, 3
     550             :                            forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
     551      572148 :                                        - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
     552             :                            forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
     553      573392 :                                        - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
     554             :                         END DO
     555             :                      ELSE
     556        2328 :                         DO i = 1, 3
     557             :                            forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
     558      211626 :                                        - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
     559             :                            forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
     560      212208 :                                        - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
     561             :                         END DO
     562             :                      END IF
     563             :                   END DO
     564        5704 :                   force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fab*fpolvec(idir)*forcea(1:3)
     565        5704 :                   force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fab*fpolvec(idir)*forceb(1:3)
     566        1472 :                   IF (use_virial) THEN
     567           0 :                      f0 = -fab*fpolvec(idir)
     568           0 :                      CALL virial_pair_force(virial%pv_virial, f0, forcea, ra)
     569           0 :                      CALL virial_pair_force(virial%pv_virial, f0, forceb, rb)
     570             :                   END IF
     571             : 
     572             :                END DO
     573          46 :                CALL neighbor_list_iterator_release(nl_iterator)
     574          46 :                DEALLOCATE (basis_set_list)
     575             : 
     576             :             END IF
     577             :          END IF
     578             :       END DO
     579             : 
     580             :       ! Energy
     581        4640 :       ener_field = 0.0_dp
     582             :       ti = 0.0_dp
     583        4640 :       DO idir = 1, 3
     584             :          ! make sure the total normalized polarization is within [-1:1]
     585        3480 :          cqi(idir) = qi(idir) + ci(idir)
     586        3480 :          IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
     587        3480 :          IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
     588             :          ! now check for log branch
     589        3480 :          IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
     590           0 :             ti(idir) = (efield%polarisation(idir) - cqi(idir))/pi
     591           0 :             DO i = 1, 10
     592           0 :                cqi(idir) = cqi(idir) + SIGN(1.0_dp, ti(idir))*twopi
     593           0 :                IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
     594             :             END DO
     595             :          END IF
     596        4640 :          ener_field = ener_field + fpolvec(idir)*cqi(idir)
     597             :       END DO
     598             : 
     599             :       ! update the references
     600        1160 :       IF (calculate_forces) THEN
     601             :          ! check for smoothness of energy surface
     602         112 :          IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(fpolvec))) THEN
     603           0 :             CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
     604             :          END IF
     605          28 :          efield%field_energy = ener_field
     606         112 :          efield%polarisation(:) = cqi(:)
     607             :       END IF
     608        1160 :       energy%efield = ener_field
     609             : 
     610        1160 :       IF (.NOT. just_energy) THEN
     611             :          ! Add the result to mo_derivativs
     612        2254 :          DO ispin = 1, dft_control%nspins
     613        2254 :             CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin), mo_derivs(ispin)%matrix)
     614             :          END DO
     615        1108 :          IF (use_virial) THEN
     616           0 :             ti = 0.0_dp
     617           0 :             DO i = 1, 3
     618           0 :                DO j = 1, 3
     619           0 :                   ti(j) = ti(j) + hmat(j, i)*cqi(i)
     620             :                END DO
     621             :             END DO
     622           0 :             DO i = 1, 3
     623           0 :                DO j = 1, 3
     624           0 :                   virial%pv_virial(i, j) = virial%pv_virial(i, j) - fieldpol(i)*ti(j)
     625             :                END DO
     626             :             END DO
     627             :          END IF
     628             :       END IF
     629             : 
     630        2380 :       DO ispin = 1, dft_control%nspins
     631        1220 :          CALL cp_cfm_release(eigrmat(ispin))
     632        1220 :          CALL cp_cfm_release(inv_mat(ispin))
     633        1220 :          CALL cp_fm_release(mo_derivs_tmp(ispin))
     634        1220 :          IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
     635        4820 :          DO i = 1, SIZE(op_fm_set, 1)
     636        2440 :             CALL cp_fm_release(opvec(i, ispin))
     637        2440 :             CALL cp_fm_release(op_fm_set(i, ispin))
     638        3660 :             CALL cp_fm_release(inv_work(i, ispin))
     639             :          END DO
     640             :       END DO
     641        1160 :       DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
     642        1160 :       DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
     643             : 
     644        1160 :       IF (calculate_forces) THEN
     645          84 :          DO ikind = 1, SIZE(atomic_kind_set)
     646         676 :             CALL para_env%sum(force(ikind)%efield)
     647             :          END DO
     648          28 :          DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
     649         112 :          DO i = 1, 3
     650          84 :             DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
     651         112 :             DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
     652             :          END DO
     653          28 :          CALL dbcsr_deallocate_matrix_set(tempmat)
     654             :       END IF
     655        1160 :       CALL timestop(handle)
     656             : 
     657        4640 :    END SUBROUTINE qs_efield_derivatives
     658             : 
     659             : ! **************************************************************************************************
     660             : !> \brief ...
     661             : !> \param qs_env ...
     662             : !> \param just_energy ...
     663             : !> \param calculate_forces ...
     664             : ! **************************************************************************************************
     665         898 :    SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
     666             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     667             :       LOGICAL, INTENT(IN)                                :: just_energy, calculate_forces
     668             : 
     669             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_dispfield_derivatives'
     670             : 
     671             :       COMPLEX(dp)                                        :: zdet, zdeta, zi(3)
     672             :       INTEGER :: handle, i, ia, iatom, icol, idir, ikind, iodeb, irow, iset, ispin, jatom, jkind, &
     673             :          jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, &
     674             :          sgfa, sgfb
     675         898 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     676         898 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     677         898 :                                                             npgfb, nsgfa, nsgfb
     678         898 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     679             :       LOGICAL                                            :: found, uniform, use_virial
     680             :       REAL(dp)                                           :: charge, ci(3), cqi(3), dab, dd, di(3), &
     681             :                                                             ener_field, fab, fieldpol(3), focc, &
     682             :                                                             hmat(3, 3), occ, omega, qi(3), &
     683             :                                                             rlog(3), zlog(3)
     684             :       REAL(dp), DIMENSION(3)                             :: dfilter, forcea, forceb, kvec, ra, rab, &
     685             :                                                             rb, ria
     686        1796 :       REAL(dp), DIMENSION(:, :), POINTER                 :: cosab, iblock, rblock, sinab, work
     687        2694 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: dcosab, dsinab, force_tmp
     688         898 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     689         898 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
     690         898 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     691       16164 :       TYPE(block_p_type), DIMENSION(3, 2)                :: dcost, dsint
     692             :       TYPE(cell_type), POINTER                           :: cell
     693         898 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: eigrmat, inv_mat
     694             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     695         898 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: mo_coeff_tmp
     696         898 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: inv_work, mo_derivs_tmp, op_fm_set, opvec
     697             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     698         898 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, mo_derivs
     699         898 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: tempmat
     700             :       TYPE(dbcsr_type), POINTER                          :: cosmat, mo_coeff_b, sinmat
     701             :       TYPE(dft_control_type), POINTER                    :: dft_control
     702             :       TYPE(efield_berry_type), POINTER                   :: efield
     703         898 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     704             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     705         898 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     706             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     707             :       TYPE(neighbor_list_iterator_p_type), &
     708         898 :          DIMENSION(:), POINTER                           :: nl_iterator
     709             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     710         898 :          POINTER                                         :: sab_orb
     711         898 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     712             :       TYPE(qs_energy_type), POINTER                      :: energy
     713         898 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     714         898 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     715             :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     716             :       TYPE(virial_type), POINTER                         :: virial
     717             : 
     718         898 :       CALL timeset(routineN, handle)
     719             : 
     720         898 :       NULLIFY (dft_control, cell, particle_set)
     721             :       CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
     722         898 :                       particle_set=particle_set, virial=virial)
     723         898 :       NULLIFY (qs_kind_set, efield, para_env, sab_orb)
     724             :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
     725         898 :                       efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
     726             : 
     727             :       ! calculate stress only if forces requested also
     728         898 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     729           0 :       use_virial = use_virial .AND. calculate_forces
     730             :       ! disable stress calculation
     731             :       IF (use_virial) THEN
     732           0 :          CPABORT("Stress tensor for periodic D-field not implemented")
     733             :       END IF
     734             : 
     735        3592 :       dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
     736             : 
     737        3592 :       fieldpol = dft_control%period_efield%polarisation
     738        6286 :       fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
     739        3592 :       fieldpol = fieldpol*dft_control%period_efield%strength
     740             : 
     741         898 :       omega = cell%deth
     742       11674 :       hmat = cell%hmat(:, :)/(twopi*omega)
     743             : 
     744             :       ! nuclear contribution to polarization
     745         898 :       natom = SIZE(particle_set)
     746         898 :       IF (calculate_forces) THEN
     747          10 :          CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
     748          10 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
     749          40 :          ALLOCATE (force_tmp(natom, 3, 3))
     750         310 :          force_tmp = 0.0_dp
     751             :       END IF
     752        3592 :       zi(:) = CMPLX(1._dp, 0._dp, dp)
     753        2694 :       DO ia = 1, natom
     754        1796 :          CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
     755        1796 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
     756        7184 :          ria = particle_set(ia)%r
     757        7184 :          ria = pbc(ria, cell)
     758        7184 :          DO idir = 1, 3
     759       21552 :             kvec(:) = twopi*cell%h_inv(idir, :)
     760       21552 :             dd = SUM(kvec(:)*ria(:))
     761        5388 :             zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
     762        7184 :             zi(idir) = zi(idir)*zdeta
     763             :          END DO
     764        4490 :          IF (calculate_forces) THEN
     765          20 :             IF (para_env%mepos == 0) THEN
     766          40 :                DO i = 1, 3
     767          40 :                   force_tmp(ia, i, i) = force_tmp(ia, i, i) + charge/omega
     768             :                END DO
     769             :             END IF
     770             :          END IF
     771             :       END DO
     772        3592 :       rlog = AIMAG(LOG(zi))
     773             : 
     774             :       ! check uniform occupation
     775         898 :       NULLIFY (mos)
     776         898 :       CALL get_qs_env(qs_env=qs_env, mos=mos)
     777        1796 :       DO ispin = 1, dft_control%nspins
     778         898 :          CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
     779        1796 :          IF (.NOT. uniform) THEN
     780           0 :             CPABORT("Berry phase moments for non uniform MO occupation numbers not implemented")
     781             :          END IF
     782             :       END DO
     783             : 
     784             :       ! initialize all work matrices needed
     785         898 :       NULLIFY (mo_derivs)
     786         898 :       CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
     787        5388 :       ALLOCATE (op_fm_set(2, dft_control%nspins))
     788        5388 :       ALLOCATE (opvec(2, dft_control%nspins))
     789        3592 :       ALLOCATE (eigrmat(dft_control%nspins))
     790        3592 :       ALLOCATE (inv_mat(dft_control%nspins))
     791        5388 :       ALLOCATE (inv_work(2, dft_control%nspins))
     792        6286 :       ALLOCATE (mo_derivs_tmp(3, SIZE(mo_derivs)))
     793        3592 :       ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
     794             : 
     795             :       ! Allocate temp matrices for the wavefunction derivatives
     796        1796 :       DO ispin = 1, dft_control%nspins
     797         898 :          NULLIFY (tmp_fm_struct, mo_coeff)
     798         898 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
     799             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
     800         898 :                                   ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
     801         898 :          CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
     802        3592 :          DO i = 1, 3
     803        2694 :             CALL cp_fm_create(mo_derivs_tmp(i, ispin), mo_coeff%matrix_struct)
     804        3592 :             CALL cp_fm_set_all(matrix=mo_derivs_tmp(i, ispin), alpha=0.0_dp)
     805             :          END DO
     806        2694 :          DO i = 1, SIZE(op_fm_set, 1)
     807        1796 :             CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
     808        1796 :             CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
     809        2694 :             CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
     810             :          END DO
     811         898 :          CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
     812         898 :          CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
     813        2694 :          CALL cp_fm_struct_release(tmp_fm_struct)
     814             :       END DO
     815             :       ! temp matrices for force calculation
     816         898 :       IF (calculate_forces) THEN
     817          10 :          NULLIFY (matrix_s)
     818          10 :          CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
     819          60 :          ALLOCATE (tempmat(2, dft_control%nspins))
     820          20 :          DO ispin = 1, dft_control%nspins
     821          10 :             ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
     822          10 :             CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
     823          10 :             CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
     824          10 :             CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
     825          20 :             CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
     826             :          END DO
     827             :          ! integration
     828          10 :          CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
     829         100 :          ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
     830          90 :          ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
     831          10 :          lsab = MAX(lsab, ldab)
     832          50 :          DO i = 1, 3
     833         210 :             ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
     834         220 :             ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
     835             :          END DO
     836             :       END IF
     837             : 
     838             :       !Start the MO derivative calculation
     839             :       !loop over all cell vectors
     840        3592 :       DO idir = 1, 3
     841        2694 :          zi(idir) = z_zero
     842        2694 :          cosmat => efield%cosmat(idir)%matrix
     843        2694 :          sinmat => efield%sinmat(idir)%matrix
     844             :          !evaluate the expression needed for the derivative (S_berry * C  and [C^T S_berry C]^-1)
     845             :          !first step S_berry * C  and C^T S_berry C
     846        5388 :          DO ispin = 1, dft_control%nspins ! spin
     847        2694 :             IF (mos(ispin)%use_mo_coeff_b) THEN
     848        2694 :                CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
     849        2694 :                CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
     850             :             ELSE
     851           0 :                CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
     852           0 :                mo_coeff_tmp(ispin) = mo_coeff
     853             :             END IF
     854        2694 :             CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
     855             :             CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
     856        2694 :                                op_fm_set(1, ispin))
     857        2694 :             CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
     858             :             CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
     859        5388 :                                op_fm_set(2, ispin))
     860             :          END DO
     861             :          !second step invert C^T S_berry C
     862        2694 :          zdet = z_one
     863        5388 :          DO ispin = 1, dft_control%nspins
     864        2694 :             CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
     865        2694 :             CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
     866        2694 :             CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
     867        2694 :             CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
     868        5388 :             zdet = zdet*zdeta
     869             :          END DO
     870        2694 :          zi(idir) = zdet**occ
     871        2694 :          zlog(idir) = AIMAG(LOG(zi(idir)))
     872             : 
     873        2694 :          IF (.NOT. just_energy) THEN
     874             :             !compute the orbital derivative
     875        5388 :             DO ispin = 1, dft_control%nspins
     876       35022 :                inv_work(1, ispin)%local_data(:, :) = REAL(inv_mat(ispin)%local_data(:, :), dp)
     877       35022 :                inv_work(2, ispin)%local_data(:, :) = AIMAG(inv_mat(ispin)%local_data(:, :))
     878        2694 :                CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
     879       13470 :                DO i = 1, 3
     880        8082 :                   focc = hmat(idir, i)
     881             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
     882        8082 :                                      1.0_dp, mo_derivs_tmp(idir, ispin))
     883             :                   CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
     884       10776 :                                      1.0_dp, mo_derivs_tmp(idir, ispin))
     885             :                END DO
     886             :             END DO
     887             :          END IF
     888             : 
     889             :          !compute nuclear forces
     890        3592 :          IF (calculate_forces) THEN
     891          30 :             nkind = SIZE(qs_kind_set)
     892          30 :             natom = SIZE(particle_set)
     893         120 :             kvec(:) = twopi*cell%h_inv(idir, :)
     894             : 
     895             :             ! calculate: C [C^T S_berry C]^(-1) C^T
     896             :             ! Store this matrix in DBCSR form (only S overlap blocks)
     897          60 :             DO ispin = 1, dft_control%nspins
     898          30 :                CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
     899          30 :                CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
     900          30 :                CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
     901             :                CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
     902          30 :                                   opvec(1, ispin))
     903             :                CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
     904          30 :                                   opvec(2, ispin))
     905             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
     906          30 :                                           matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
     907             :                CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
     908          90 :                                           matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
     909             :             END DO
     910             : 
     911             :             ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
     912         150 :             ALLOCATE (basis_set_list(nkind))
     913          90 :             DO ikind = 1, nkind
     914          60 :                qs_kind => qs_kind_set(ikind)
     915          60 :                CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
     916          90 :                IF (ASSOCIATED(basis_set_a)) THEN
     917          60 :                   basis_set_list(ikind)%gto_basis_set => basis_set_a
     918             :                ELSE
     919           0 :                   NULLIFY (basis_set_list(ikind)%gto_basis_set)
     920             :                END IF
     921             :             END DO
     922             :             !
     923          30 :             CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     924         585 :             DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     925             :                CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     926         555 :                                       iatom=iatom, jatom=jatom, r=rab)
     927         555 :                basis_set_a => basis_set_list(ikind)%gto_basis_set
     928         555 :                IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     929         555 :                basis_set_b => basis_set_list(jkind)%gto_basis_set
     930         555 :                IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     931             :                ! basis ikind
     932         555 :                first_sgfa => basis_set_a%first_sgf
     933         555 :                la_max => basis_set_a%lmax
     934         555 :                la_min => basis_set_a%lmin
     935         555 :                npgfa => basis_set_a%npgf
     936         555 :                nseta = basis_set_a%nset
     937         555 :                nsgfa => basis_set_a%nsgf_set
     938         555 :                rpgfa => basis_set_a%pgf_radius
     939         555 :                set_radius_a => basis_set_a%set_radius
     940         555 :                sphi_a => basis_set_a%sphi
     941         555 :                zeta => basis_set_a%zet
     942             :                ! basis jkind
     943         555 :                first_sgfb => basis_set_b%first_sgf
     944         555 :                lb_max => basis_set_b%lmax
     945         555 :                lb_min => basis_set_b%lmin
     946         555 :                npgfb => basis_set_b%npgf
     947         555 :                nsetb = basis_set_b%nset
     948         555 :                nsgfb => basis_set_b%nsgf_set
     949         555 :                rpgfb => basis_set_b%pgf_radius
     950         555 :                set_radius_b => basis_set_b%set_radius
     951         555 :                sphi_b => basis_set_b%sphi
     952         555 :                zetb => basis_set_b%zet
     953             : 
     954         555 :                ldsa = SIZE(sphi_a, 1)
     955         555 :                ldsb = SIZE(sphi_b, 1)
     956         555 :                ra(:) = pbc(particle_set(iatom)%r(:), cell)
     957        2220 :                rb(:) = ra + rab
     958         555 :                dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
     959             : 
     960         555 :                IF (iatom <= jatom) THEN
     961         354 :                   irow = iatom
     962         354 :                   icol = jatom
     963             :                ELSE
     964         201 :                   irow = jatom
     965         201 :                   icol = iatom
     966             :                END IF
     967             : 
     968         555 :                IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
     969             :                   fab = 1.0_dp*occ
     970             :                ELSE
     971         525 :                   fab = 2.0_dp*occ
     972             :                END IF
     973             : 
     974        2220 :                DO i = 1, 3
     975      454545 :                   dcost(i, 1)%block = 0.0_dp
     976      454545 :                   dsint(i, 1)%block = 0.0_dp
     977      454545 :                   dcost(i, 2)%block = 0.0_dp
     978      455100 :                   dsint(i, 2)%block = 0.0_dp
     979             :                END DO
     980             : 
     981        1665 :                DO iset = 1, nseta
     982        1110 :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     983        1110 :                   sgfa = first_sgfa(1, iset)
     984        3885 :                   DO jset = 1, nsetb
     985        2220 :                      IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     986        1140 :                      ncob = npgfb(jset)*ncoset(lb_max(jset))
     987        1140 :                      sgfb = first_sgfb(1, jset)
     988             :                      ! Calculate the primitive integrals (da|b)
     989             :                      CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     990             :                                  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     991        1140 :                                  ra, rb, kvec, cosab, sinab, dcosab, dsinab)
     992        4560 :                      DO i = 1, 3
     993             :                         CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
     994             :                                           ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
     995             :                                           ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
     996        4560 :                                           dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
     997             :                      END DO
     998             :                      ! Calculate the primitive integrals (a|db)
     999             :                      CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
    1000             :                                  la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
    1001        1140 :                                  rb, ra, kvec, cosab, sinab, dcosab, dsinab)
    1002        5670 :                      DO i = 1, 3
    1003      560556 :                         dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i))
    1004      560556 :                         dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i))
    1005             :                         CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
    1006             :                                           ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
    1007             :                                           ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
    1008        5640 :                                           dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
    1009             :                      END DO
    1010             :                   END DO
    1011             :                END DO
    1012         555 :                forcea = 0.0_dp
    1013         555 :                forceb = 0.0_dp
    1014        1110 :                DO ispin = 1, dft_control%nspins
    1015         555 :                   NULLIFY (rblock, iblock)
    1016             :                   CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
    1017         555 :                                          row=irow, col=icol, BLOCK=rblock, found=found)
    1018         555 :                   CPASSERT(found)
    1019             :                   CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
    1020         555 :                                          row=irow, col=icol, BLOCK=iblock, found=found)
    1021         555 :                   CPASSERT(found)
    1022         555 :                   n1 = SIZE(rblock, 1)
    1023         555 :                   n2 = SIZE(rblock, 2)
    1024         555 :                   CPASSERT(SIZE(iblock, 1) == n1)
    1025         555 :                   CPASSERT(SIZE(iblock, 2) == n2)
    1026         555 :                   CPASSERT(lsab >= n1)
    1027         555 :                   CPASSERT(lsab >= n2)
    1028        2220 :                   IF (iatom <= jatom) THEN
    1029        1416 :                      DO i = 1, 3
    1030             :                         forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
    1031      160542 :                                     - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
    1032             :                         forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
    1033      160896 :                                     - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
    1034             :                      END DO
    1035             :                   ELSE
    1036         804 :                      DO i = 1, 3
    1037             :                         forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
    1038       85023 :                                     - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
    1039             :                         forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
    1040       85224 :                                     - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
    1041             :                      END DO
    1042             :                   END IF
    1043             :                END DO
    1044        2250 :                DO i = 1, 3
    1045        6660 :                   force_tmp(iatom, :, i) = force_tmp(iatom, :, i) - fab*hmat(i, idir)*forcea(:)
    1046        7215 :                   force_tmp(jatom, :, i) = force_tmp(jatom, :, i) - fab*hmat(i, idir)*forceb(:)
    1047             :                END DO
    1048             :             END DO
    1049          30 :             CALL neighbor_list_iterator_release(nl_iterator)
    1050          30 :             DEALLOCATE (basis_set_list)
    1051             :          END IF
    1052             :       END DO
    1053             : 
    1054             :       ! make sure the total normalized polarization is within [-1:1]
    1055        3592 :       DO idir = 1, 3
    1056        2694 :          cqi(idir) = rlog(idir) + zlog(idir)
    1057        2694 :          IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
    1058        2694 :          IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
    1059             :          ! now check for log branch
    1060        3592 :          IF (calculate_forces) THEN
    1061          30 :             IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
    1062           0 :                di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
    1063           0 :                DO i = 1, 10
    1064           0 :                   cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi
    1065           0 :                   IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
    1066             :                END DO
    1067             :             END IF
    1068             :          END IF
    1069             :       END DO
    1070        3592 :       DO idir = 1, 3
    1071        2694 :          qi(idir) = 0.0_dp
    1072        2694 :          ci(idir) = 0.0_dp
    1073       11674 :          DO i = 1, 3
    1074       10776 :             ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
    1075             :          END DO
    1076             :       END DO
    1077             : 
    1078             :       ! update the references
    1079         898 :       IF (calculate_forces) THEN
    1080          40 :          ener_field = SUM(ci)
    1081             :          ! check for smoothness of energy surface
    1082         130 :          IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN
    1083           0 :             CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
    1084             :          END IF
    1085          10 :          efield%field_energy = ener_field
    1086          40 :          efield%polarisation(:) = cqi(:)
    1087             :       END IF
    1088             : 
    1089             :       ! Energy
    1090         898 :       ener_field = 0.0_dp
    1091        3592 :       DO i = 1, 3
    1092        3592 :          ener_field = ener_field + dfilter(i)*(fieldpol(i) - 2._dp*twopi*ci(i))**2
    1093             :       END DO
    1094         898 :       energy%efield = 0.25_dp*omega/twopi*ener_field
    1095             : 
    1096             :       ! debugging output
    1097         898 :       IF (para_env%is_source()) THEN
    1098             :          iodeb = -1
    1099             :          IF (iodeb > 0) THEN
    1100             :             WRITE (iodeb, '(A,T61,F20.10)') "  Polarisation Quantum:  ", 2._dp*twopi*twopi*hmat(3, 3)
    1101             :             WRITE (iodeb, '(A,T21,3F20.10)') "  Polarisation: ", 2._dp*twopi*ci(1:3)
    1102             :             WRITE (iodeb, '(A,T21,3F20.10)') "  Displacement: ", fieldpol(1:3)
    1103             :             WRITE (iodeb, '(A,T21,3F20.10)') "  E-Field:      ", ((fieldpol(i) - 2._dp*twopi*ci(i)), i=1, 3)
    1104             :             WRITE (iodeb, '(A,T61,F20.10)') "  Disp Free Energy:", energy%efield
    1105             :          END IF
    1106             :       END IF
    1107             : 
    1108         898 :       IF (.NOT. just_energy) THEN
    1109        3592 :          DO i = 1, 3
    1110        3592 :             di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i)
    1111             :          END DO
    1112             :          ! Add the result to mo_derivativs
    1113        1796 :          DO ispin = 1, dft_control%nspins
    1114         898 :             CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_coeff_tmp(ispin))
    1115        4490 :             DO idir = 1, 3
    1116             :                CALL cp_fm_scale_and_add(1.0_dp, mo_coeff_tmp(ispin), di(idir), &
    1117        3592 :                                         mo_derivs_tmp(idir, ispin))
    1118             :             END DO
    1119             :          END DO
    1120        1796 :          DO ispin = 1, dft_control%nspins
    1121        1796 :             CALL copy_fm_to_dbcsr(mo_coeff_tmp(ispin), mo_derivs(ispin)%matrix)
    1122             :          END DO
    1123             :       END IF
    1124             : 
    1125         898 :       IF (calculate_forces) THEN
    1126          40 :          DO i = 1, 3
    1127         100 :             DO ia = 1, natom
    1128          60 :                CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
    1129          60 :                iatom = atom_of_kind(ia)
    1130         450 :                force(ikind)%efield(1:3, iatom) = force(ikind)%efield(1:3, iatom) + di(i)*force_tmp(ia, 1:3, i)
    1131             :             END DO
    1132             :          END DO
    1133             :       END IF
    1134             : 
    1135        1796 :       DO ispin = 1, dft_control%nspins
    1136         898 :          CALL cp_cfm_release(eigrmat(ispin))
    1137         898 :          CALL cp_cfm_release(inv_mat(ispin))
    1138         898 :          IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
    1139        3592 :          DO i = 1, 3
    1140        3592 :             CALL cp_fm_release(mo_derivs_tmp(i, ispin))
    1141             :          END DO
    1142        3592 :          DO i = 1, SIZE(op_fm_set, 1)
    1143        1796 :             CALL cp_fm_release(opvec(i, ispin))
    1144        1796 :             CALL cp_fm_release(op_fm_set(i, ispin))
    1145        2694 :             CALL cp_fm_release(inv_work(i, ispin))
    1146             :          END DO
    1147             :       END DO
    1148         898 :       DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
    1149         898 :       DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
    1150             : 
    1151         898 :       IF (calculate_forces) THEN
    1152          30 :          DO ikind = 1, SIZE(atomic_kind_set)
    1153         190 :             CALL para_env%sum(force(ikind)%efield)
    1154             :          END DO
    1155          10 :          DEALLOCATE (force_tmp)
    1156          10 :          DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
    1157          40 :          DO i = 1, 3
    1158          30 :             DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
    1159          40 :             DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
    1160             :          END DO
    1161          10 :          CALL dbcsr_deallocate_matrix_set(tempmat)
    1162             :       END IF
    1163         898 :       CALL timestop(handle)
    1164             : 
    1165        3592 :    END SUBROUTINE qs_dispfield_derivatives
    1166             : 
    1167             : ! **************************************************************************************************
    1168             : !> \brief ...
    1169             : !> \param cos_block ...
    1170             : !> \param sin_block ...
    1171             : !> \param ncoa ...
    1172             : !> \param nsgfa ...
    1173             : !> \param sgfa ...
    1174             : !> \param sphi_a ...
    1175             : !> \param ldsa ...
    1176             : !> \param ncob ...
    1177             : !> \param nsgfb ...
    1178             : !> \param sgfb ...
    1179             : !> \param sphi_b ...
    1180             : !> \param ldsb ...
    1181             : !> \param cosab ...
    1182             : !> \param sinab ...
    1183             : !> \param ldab ...
    1184             : !> \param work ...
    1185             : !> \param ldwork ...
    1186             : ! **************************************************************************************************
    1187       23304 :    SUBROUTINE contract_all(cos_block, sin_block, &
    1188       46608 :                            ncoa, nsgfa, sgfa, sphi_a, ldsa, &
    1189       46608 :                            ncob, nsgfb, sgfb, sphi_b, ldsb, &
    1190       23304 :                            cosab, sinab, ldab, work, ldwork)
    1191             : 
    1192             :       REAL(dp), DIMENSION(:, :), POINTER                 :: cos_block, sin_block
    1193             :       INTEGER, INTENT(IN)                                :: ncoa, nsgfa, sgfa
    1194             :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: sphi_a
    1195             :       INTEGER, INTENT(IN)                                :: ldsa, ncob, nsgfb, sgfb
    1196             :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: sphi_b
    1197             :       INTEGER, INTENT(IN)                                :: ldsb
    1198             :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: cosab, sinab
    1199             :       INTEGER, INTENT(IN)                                :: ldab
    1200             :       REAL(dp), DIMENSION(:, :)                          :: work
    1201             :       INTEGER, INTENT(IN)                                :: ldwork
    1202             : 
    1203             : ! Calculate cosine
    1204             : 
    1205             :       CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, cosab(1, 1), ldab, &
    1206       23304 :                  sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
    1207             : 
    1208             :       CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
    1209       23304 :                  work(1, 1), ldwork, 1.0_dp, cos_block(sgfa, sgfb), SIZE(cos_block, 1))
    1210             : 
    1211             :       ! Calculate sine
    1212             :       CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sinab(1, 1), ldab, &
    1213       23304 :                  sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
    1214             : 
    1215             :       CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
    1216       23304 :                  work(1, 1), ldwork, 1.0_dp, sin_block(sgfa, sgfb), SIZE(sin_block, 1))
    1217             : 
    1218       23304 :    END SUBROUTINE contract_all
    1219             : 
    1220             : END MODULE qs_efield_berry

Generated by: LCOV version 1.15