efield_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief all routins needed for a nonperiodic  electric field
! **************************************************************************************************
 
MODULE efield_utils
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind,&
                                              get_atomic_kind_set
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE cp_control_types,                ONLY: dft_control_type,&
                                              efield_type
   USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
                                              dbcsr_copy,&
                                              dbcsr_p_type,&
                                              dbcsr_set
   USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
                                              dbcsr_deallocate_matrix_set
   USE input_constants,                 ONLY: constant_env,&
                                              custom_env,&
                                              gaussian_env,&
                                              ramp_env
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi,&
                                              twopi
   USE message_passing,                 ONLY: mp_para_env_type
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: c_light_au
   USE qs_efield_local,                 ONLY: qs_efield_dipole_force
   USE qs_energy_types,                 ONLY: qs_energy_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_force_types,                  ONLY: qs_force_type
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              qs_kind_type
   USE qs_moments,                      ONLY: build_local_moment_matrix
   USE qs_rho_types,                    ONLY: qs_rho_get,&
                                              qs_rho_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_utils'
 
! *** Public subroutines ***
 
   PUBLIC :: efield_potential_lengh_gauge, &
             calculate_ecore_efield, &
             make_field
 
CONTAINS
 
! **************************************************************************************************
!> \brief Replace the original implementation of the electric-electronic
!>        interaction in the length gauge. This calculation is no longer done in
!>        the grid but using matrices to match the velocity gauge implementation.
!>        Note: The energy is stored in energy%core and computed later on.
!> \param qs_env ...
!> \param calculate_forces add the electronic force from the derivative of the dipole integrals
!> \author Guillaume Le Breton (02.23)
! **************************************************************************************************
 
   SUBROUTINE efield_potential_lengh_gauge(qs_env, calculate_forces)
 
      TYPE(qs_environment_type), POINTER                 :: qs_env
      LOGICAL, INTENT(IN), OPTIONAL                      :: calculate_forces
 
      CHARACTER(len=*), PARAMETER :: routinen = 'efield_potential_lengh_gauge'
 
      INTEGER                                            :: handle, i, iatom, image, natom
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
      LOGICAL                                            :: my_force
      REAL(kind=dp)                                      :: field(3)
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: fdip
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p, matrix_s, moments
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_rho_type), POINTER                         :: rho
 
      NULLIFY (dft_control)
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, &
                      dft_control=dft_control, &
                      matrix_h_kp=matrix_h, &
                      matrix_s=matrix_s)
 
      NULLIFY (moments)
      CALL dbcsr_allocate_matrix_set(moments, 3)
      DO i = 1, 3
         ALLOCATE (moments(i)%matrix)
         CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
         CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
      END DO
 
      CALL build_local_moment_matrix(qs_env, moments, 1)
 
      CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
 
      DO i = 1, 3
         DO image = 1, dft_control%nimages
            CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i))
         END DO
      END DO
 
      CALL dbcsr_deallocate_matrix_set(moments)
 
      my_force = .false.
      IF (PRESENT(calculate_forces)) my_force = calculate_forces
      IF (my_force) THEN
         ! The dipole integrals are real and symmetric: only the real part of P contributes
         CALL get_qs_env(qs_env, rho=rho, force=force, atomic_kind_set=atomic_kind_set, &
                         natom=natom, para_env=para_env)
         CALL qs_rho_get(rho, rho_ao=matrix_p)
         ALLOCATE (fdip(3, natom))
         fdip = 0.0_dp
         CALL qs_efield_dipole_force(qs_env, matrix_p, field, [0.0_dp, 0.0_dp, 0.0_dp], fdip)
         CALL para_env%sum(fdip)
         CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
         DO iatom = 1, natom
            force(kind_of(iatom))%efield(:, atom_of_kind(iatom)) = &
               force(kind_of(iatom))%efield(:, atom_of_kind(iatom)) + fdip(:, iatom)
         END DO
         DEALLOCATE (fdip)
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE efield_potential_lengh_gauge
 
! **************************************************************************************************
!> \brief computes the amplitude of the efield within a given envelop
!> \param dft_control ...
!> \param field ...
!> \param sim_step ...
!> \param sim_time ...
!> \author Florian Schiffmann (02.09)
! **************************************************************************************************
 
   SUBROUTINE make_field(dft_control, field, sim_step, sim_time)
      TYPE(dft_control_type), INTENT(IN)                 :: dft_control
      REAL(dp), INTENT(OUT)                              :: field(3)
      INTEGER, INTENT(IN)                                :: sim_step
      REAL(kind=dp), INTENT(IN)                          :: sim_time
 
      INTEGER                                            :: i, lower, nfield, upper
      REAL(dp)                                           :: e_0, env, nu, pol(3), strength
      REAL(kind=dp)                                      :: dt
      TYPE(efield_type), POINTER                         :: efield
 
      field = 0._dp
      nfield = SIZE(dft_control%efield_fields)
      DO i = 1, nfield
         efield => dft_control%efield_fields(i)%efield
         nu = 0.0_dp
         IF (.NOT. efield%envelop_id == custom_env .AND. efield%wavelength > epsilon(0.0_dp)) THEN
            nu = c_light_au/(efield%wavelength) !in case of a custom efield we do not need nu
         END IF
         e_0 = efield%amplitude
 
         IF (dot_product(efield%polarisation, efield%polarisation) == 0) THEN
            pol(:) = 1.0_dp/3.0_dp
         ELSE
            pol(:) = efield%polarisation(:)/norm2(efield%polarisation)
         END IF
 
         SELECT CASE (efield%envelop_id)
         CASE (constant_env)
            IF (sim_step >= efield%envelop_i_vars(1) .AND. &
                (sim_step <= efield%envelop_i_vars(2) .OR. efield%envelop_i_vars(2) < 0)) THEN
               field = field + e_0*cos(sim_time*nu*twopi + &
                                       efield%phase_offset*pi)*pol(:)
            END IF
         CASE (ramp_env)
            IF (sim_step >= efield%envelop_i_vars(1) .AND. sim_step <= efield%envelop_i_vars(2)) THEN
               strength = e_0*(sim_step - efield%envelop_i_vars(1))/(efield%envelop_i_vars(2) - efield%envelop_i_vars(1))
            END IF
            IF (sim_step >= efield%envelop_i_vars(3) .AND. sim_step <= efield%envelop_i_vars(4)) THEN
               strength = e_0*(efield%envelop_i_vars(4) - sim_step)/(efield%envelop_i_vars(4) - efield%envelop_i_vars(3))
            END IF
            IF (sim_step > efield%envelop_i_vars(4) .AND. efield%envelop_i_vars(4) > 0) strength = 0.0_dp
            IF (sim_step <= efield%envelop_i_vars(1)) strength = 0.0_dp
            field = field + strength*cos(sim_time*nu*twopi + &
                                         efield%phase_offset*pi)*pol(:)
         CASE (gaussian_env)
            env = exp(-0.5_dp*((sim_time - efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
            field = field + e_0*env*cos(sim_time*nu*twopi + &
                                        efield%phase_offset*pi)*pol(:)
         CASE (custom_env)
            dt = efield%envelop_r_vars(1)
            IF (sim_time < (SIZE(efield%envelop_r_vars) - 2)*dt) THEN
               !make a linear interpolation between the two next points
               lower = floor(sim_time/dt)
               upper = lower + 1
               strength = (efield%envelop_r_vars(lower + 2)*(upper*dt - sim_time) &
                           + efield%envelop_r_vars(upper + 2)*(sim_time - lower*dt))/dt
            ELSE
               strength = 0.0_dp
            END IF
            field = field + strength*pol(:)
         CASE default
         END SELECT
      END DO
 
   END SUBROUTINE make_field
 
! **************************************************************************************************
!> \brief Computes the force and the energy due to a efield on the cores
!>        Note: In the velocity gauge, the energy term is not added because
!>        it would lead to an unbalanced energy (center of negative charge not
!>        involved in the electric energy in this gauge).
!> \param qs_env ...
!> \param calculate_forces ...
!> \author Florian Schiffmann (02.09)
! **************************************************************************************************
 
   SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      LOGICAL, OPTIONAL                                  :: calculate_forces
 
      CHARACTER(len=*), PARAMETER :: routinen = 'calculate_ecore_efield'
 
      INTEGER                                            :: atom_a, handle, iatom, ikind, natom, &
                                                            nkind
      INTEGER, DIMENSION(:), POINTER                     :: list
      LOGICAL                                            :: my_force
      REAL(kind=dp)                                      :: efield_ener, zeff
      REAL(kind=dp), DIMENSION(3)                        :: field, r
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_energy_type), POINTER                      :: energy
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      NULLIFY (dft_control)
      CALL timeset(routinen, handle)
      CALL get_qs_env(qs_env, dft_control=dft_control)
      IF (dft_control%apply_efield_field .OR. dft_control%apply_vector_potential) THEN
         my_force = .false.
         IF (PRESENT(calculate_forces)) my_force = calculate_forces
 
         CALL get_qs_env(qs_env=qs_env, &
                         atomic_kind_set=atomic_kind_set, &
                         qs_kind_set=qs_kind_set, &
                         energy=energy, &
                         particle_set=particle_set, &
                         cell=cell)
         efield_ener = 0.0_dp
         nkind = SIZE(atomic_kind_set)
         CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
 
         DO ikind = 1, SIZE(atomic_kind_set)
            CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
            CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
 
            natom = SIZE(list)
            DO iatom = 1, natom
               IF (dft_control%apply_efield_field) THEN
                  atom_a = list(iatom)
                  r(:) = pbc(particle_set(atom_a)%r(:), cell)
                  efield_ener = efield_ener - zeff*dot_product(r, field)
               END IF
               IF (my_force) THEN
                  CALL get_qs_env(qs_env=qs_env, force=force)
                  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
               END IF
            END DO
 
         END DO
         IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
      END IF
      CALL timestop(handle)
   END SUBROUTINE calculate_ecore_efield
END MODULE efield_utils