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

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

Generated by: LCOV version 2.0-1