Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief all routins needed for a nonperiodic  electric field
      10             : ! **************************************************************************************************
      11             : 
      12             : MODULE efield_utils
      13             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14             :                                               get_atomic_kind
      15             :    USE cell_types,                      ONLY: cell_type,&
      16             :                                               pbc
      17             :    USE cp_control_types,                ONLY: dft_control_type,&
      18             :                                               efield_type
      19             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      20             :                                               dbcsr_deallocate_matrix_set
      21             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      22             :                                               dbcsr_copy,&
      23             :                                               dbcsr_p_type,&
      24             :                                               dbcsr_set
      25             :    USE input_constants,                 ONLY: constant_env,&
      26             :                                               custom_env,&
      27             :                                               gaussian_env,&
      28             :                                               ramp_env
      29             :    USE kinds,                           ONLY: dp
      30             :    USE mathconstants,                   ONLY: pi
      31             :    USE particle_types,                  ONLY: particle_type
      32             :    USE qs_energy_types,                 ONLY: qs_energy_type
      33             :    USE qs_environment_types,            ONLY: get_qs_env,&
      34             :                                               qs_environment_type
      35             :    USE qs_force_types,                  ONLY: qs_force_type
      36             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      37             :                                               qs_kind_type
      38             :    USE qs_moments,                      ONLY: build_local_moment_matrix
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_utils'
      46             : 
      47             : ! *** Public subroutines ***
      48             : 
      49             :    PUBLIC :: efield_potential_lengh_gauge, &
      50             :              calculate_ecore_efield, &
      51             :              make_field
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! **************************************************************************************************
      56             : !> \brief Replace the original implementation of the electric-electronic
      57             : !>        interaction in the length gauge. This calculation is no longer done in
      58             : !>        the grid but using matrices to match the velocity gauge implementation.
      59             : !>        Note: The energy is stored in energy%core and computed later on.
      60             : !> \param qs_env ...
      61             : !> \author Guillaume Le Breton (02.23)
      62             : ! **************************************************************************************************
      63             : 
      64         214 :    SUBROUTINE efield_potential_lengh_gauge(qs_env)
      65             : 
      66             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      67             : 
      68             :       CHARACTER(len=*), PARAMETER :: routineN = 'efield_potential_lengh_gauge'
      69             : 
      70             :       INTEGER                                            :: handle, i, image
      71             :       REAL(kind=dp)                                      :: field(3)
      72         214 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, moments
      73         214 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h
      74             :       TYPE(dft_control_type), POINTER                    :: dft_control
      75             : 
      76         214 :       NULLIFY (dft_control)
      77         214 :       CALL timeset(routineN, handle)
      78             : 
      79             :       CALL get_qs_env(qs_env, &
      80             :                       dft_control=dft_control, &
      81             :                       matrix_h_kp=matrix_h, &
      82         214 :                       matrix_s=matrix_s)
      83             : 
      84         214 :       NULLIFY (moments)
      85         214 :       CALL dbcsr_allocate_matrix_set(moments, 3)
      86         856 :       DO i = 1, 3
      87         642 :          ALLOCATE (moments(i)%matrix)
      88         642 :          CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
      89         856 :          CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
      90             :       END DO
      91             : 
      92         214 :       CALL build_local_moment_matrix(qs_env, moments, 1)
      93             : 
      94         214 :       CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
      95             : 
      96         856 :       DO i = 1, 3
      97        1498 :          DO image = 1, dft_control%nimages
      98        1284 :             CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i))
      99             :          END DO
     100             :       END DO
     101             : 
     102         214 :       CALL dbcsr_deallocate_matrix_set(moments)
     103             : 
     104         214 :       CALL timestop(handle)
     105             : 
     106         214 :    END SUBROUTINE efield_potential_lengh_gauge
     107             : 
     108             : ! **************************************************************************************************
     109             : !> \brief computes the amplitude of the efield within a given envelop
     110             : !> \param dft_control ...
     111             : !> \param field ...
     112             : !> \param sim_step ...
     113             : !> \param sim_time ...
     114             : !> \author Florian Schiffmann (02.09)
     115             : ! **************************************************************************************************
     116             : 
     117         566 :    SUBROUTINE make_field(dft_control, field, sim_step, sim_time)
     118             :       TYPE(dft_control_type), INTENT(IN)                 :: dft_control
     119             :       REAL(dp), INTENT(OUT)                              :: field(3)
     120             :       INTEGER, INTENT(IN)                                :: sim_step
     121             :       REAL(KIND=dp), INTENT(IN)                          :: sim_time
     122             : 
     123             :       INTEGER                                            :: i, lower, nfield, upper
     124             :       REAL(dp)                                           :: c, env, nu, pol(3), strength
     125             :       REAL(KIND=dp)                                      :: dt
     126             :       TYPE(efield_type), POINTER                         :: efield
     127             : 
     128         566 :       c = 137.03599962875_dp
     129         566 :       field = 0._dp
     130         566 :       nfield = SIZE(dft_control%efield_fields)
     131        1132 :       DO i = 1, nfield
     132         566 :          efield => dft_control%efield_fields(i)%efield
     133         566 :          IF (.NOT. efield%envelop_id == custom_env) nu = c/(efield%wavelength) !in case of a custom efield we do not need nu
     134         566 :          strength = SQRT(efield%strength/(3.50944_dp*10.0_dp**16))
     135        2264 :          IF (DOT_PRODUCT(efield%polarisation, efield%polarisation) == 0) THEN
     136           0 :             pol(:) = 1.0_dp/3.0_dp
     137             :          ELSE
     138        3962 :             pol(:) = efield%polarisation(:)/(SQRT(DOT_PRODUCT(efield%polarisation, efield%polarisation)))
     139             :          END IF
     140        1132 :          IF (efield%envelop_id == constant_env) THEN
     141         212 :             IF (sim_step .GE. efield%envelop_i_vars(1) .AND. &
     142             :                 (sim_step .LE. efield%envelop_i_vars(2) .OR. efield%envelop_i_vars(2) .LT. 0)) THEN
     143             :                field = field + strength*COS(sim_time*nu*2.0_dp*pi + &
     144         514 :                                             efield%phase_offset*pi)*pol(:)
     145             :             END IF
     146         354 :          ELSE IF (efield%envelop_id == ramp_env) THEN
     147         118 :             IF (sim_step .GE. efield%envelop_i_vars(1) .AND. sim_step .LE. efield%envelop_i_vars(2)) &
     148          52 :                strength = strength*(sim_step - efield%envelop_i_vars(1))/(efield%envelop_i_vars(2) - efield%envelop_i_vars(1))
     149         118 :             IF (sim_step .GE. efield%envelop_i_vars(3) .AND. sim_step .LE. efield%envelop_i_vars(4)) &
     150          60 :                strength = strength*(efield%envelop_i_vars(4) - sim_step)/(efield%envelop_i_vars(4) - efield%envelop_i_vars(3))
     151         118 :             IF (sim_step .GT. efield%envelop_i_vars(4) .AND. efield%envelop_i_vars(4) .GT. 0) strength = 0.0_dp
     152         118 :             IF (sim_step .LE. efield%envelop_i_vars(1)) strength = 0.0_dp
     153             :             field = field + strength*COS(sim_time*nu*2.0_dp*pi + &
     154         472 :                                          efield%phase_offset*pi)*pol(:)
     155         236 :          ELSE IF (efield%envelop_id == gaussian_env) THEN
     156         122 :             env = EXP(-0.5_dp*((sim_time - efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
     157             :             field = field + strength*env*COS(sim_time*nu*2.0_dp*pi + &
     158         488 :                                              efield%phase_offset*pi)*pol(:)
     159         114 :          ELSE IF (efield%envelop_id == custom_env) THEN
     160         114 :             dt = efield%envelop_r_vars(1)
     161         114 :             IF (sim_time .LT. (SIZE(efield%envelop_r_vars) - 2)*dt) THEN
     162             :                !make a linear interpolation between the two next points
     163          62 :                lower = FLOOR(sim_time/dt)
     164          62 :                upper = lower + 1
     165          62 :      strength = (efield%envelop_r_vars(lower + 2)*(upper*dt - sim_time) + efield%envelop_r_vars(upper + 2)*(sim_time - lower*dt))/dt
     166             :             ELSE
     167             :                strength = 0.0_dp
     168             :             END IF
     169         456 :             field = field + strength*pol(:)
     170             :          END IF
     171             :       END DO
     172             : 
     173         566 :    END SUBROUTINE make_field
     174             : 
     175             : ! **************************************************************************************************
     176             : !> \brief Computes the force and the energy due to a efield on the cores
     177             : !>        Note: In the velocity gauge, the energy term is not added because
     178             : !>        it would lead to an unbalanced energy (center of negative charge not
     179             : !>        involved in the electric energy in this gauge).
     180             : !> \param qs_env ...
     181             : !> \param calculate_forces ...
     182             : !> \author Florian Schiffmann (02.09)
     183             : ! **************************************************************************************************
     184             : 
     185       16239 :    SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
     186             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     187             :       LOGICAL, OPTIONAL                                  :: calculate_forces
     188             : 
     189             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_efield'
     190             : 
     191             :       INTEGER                                            :: atom_a, handle, iatom, ikind, natom, &
     192             :                                                             nkind
     193       16239 :       INTEGER, DIMENSION(:), POINTER                     :: list
     194             :       LOGICAL                                            :: my_force
     195             :       REAL(KIND=dp)                                      :: efield_ener, zeff
     196             :       REAL(KIND=dp), DIMENSION(3)                        :: field, r
     197       16239 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     198             :       TYPE(cell_type), POINTER                           :: cell
     199             :       TYPE(dft_control_type), POINTER                    :: dft_control
     200       16239 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     201             :       TYPE(qs_energy_type), POINTER                      :: energy
     202       16239 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     203       16239 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     204             : 
     205       16239 :       NULLIFY (dft_control)
     206       16239 :       CALL timeset(routineN, handle)
     207       16239 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     208       16239 :       IF (dft_control%apply_efield_field .OR. dft_control%apply_vector_potential) THEN
     209         316 :          my_force = .FALSE.
     210         316 :          IF (PRESENT(calculate_forces)) my_force = calculate_forces
     211             : 
     212             :          CALL get_qs_env(qs_env=qs_env, &
     213             :                          atomic_kind_set=atomic_kind_set, &
     214             :                          qs_kind_set=qs_kind_set, &
     215             :                          energy=energy, &
     216             :                          particle_set=particle_set, &
     217         316 :                          cell=cell)
     218         316 :          efield_ener = 0.0_dp
     219         316 :          nkind = SIZE(atomic_kind_set)
     220         316 :          CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
     221             : 
     222         710 :          DO ikind = 1, SIZE(atomic_kind_set)
     223         394 :             CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
     224         394 :             CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     225             : 
     226         394 :             natom = SIZE(list)
     227        1420 :             DO iatom = 1, natom
     228         710 :                IF (dft_control%apply_efield_field) THEN
     229         498 :                   atom_a = list(iatom)
     230         498 :                   r(:) = pbc(particle_set(atom_a)%r(:), cell)
     231        1992 :                   efield_ener = efield_ener - zeff*DOT_PRODUCT(r, field)
     232             :                END IF
     233        1104 :                IF (my_force) THEN
     234         408 :                   CALL get_qs_env(qs_env=qs_env, force=force)
     235        1632 :                   force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
     236             :                END IF
     237             : !               END IF
     238             :             END DO
     239             : 
     240             :          END DO
     241         316 :          IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
     242             : !          energy%efield_core = efield_ener
     243             :       END IF
     244       16239 :       CALL timestop(handle)
     245       16239 :    END SUBROUTINE calculate_ecore_efield
     246             : END MODULE efield_utils