LCOV - code coverage report
Current view: top level - src - molecular_dipoles.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 105 105 100.0 %
Date: 2024-03-29 07:50:05 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Set of routines handling the localization for molecular properties
      10             : ! **************************************************************************************************
      11             : MODULE molecular_dipoles
      12             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      13             :                                               get_atomic_kind
      14             :    USE cell_types,                      ONLY: cell_type,&
      15             :                                               pbc
      16             :    USE cp_control_types,                ONLY: dft_control_type
      17             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      18             :                                               cp_logger_type
      19             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      20             :                                               cp_print_key_unit_nr
      21             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      22             :    USE input_section_types,             ONLY: section_get_ival,&
      23             :                                               section_vals_type,&
      24             :                                               section_vals_val_get
      25             :    USE kinds,                           ONLY: dp
      26             :    USE mathconstants,                   ONLY: twopi
      27             :    USE message_passing,                 ONLY: mp_para_env_type
      28             :    USE molecule_kind_types,             ONLY: get_molecule_kind,&
      29             :                                               molecule_kind_type
      30             :    USE molecule_types,                  ONLY: molecule_type
      31             :    USE moments_utils,                   ONLY: get_reference_point
      32             :    USE particle_types,                  ONLY: particle_type
      33             :    USE physcon,                         ONLY: debye
      34             :    USE qs_environment_types,            ONLY: get_qs_env,&
      35             :                                               qs_environment_type
      36             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      37             :                                               qs_kind_type
      38             :    USE qs_loc_types,                    ONLY: qs_loc_env_type
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    ! *** Public ***
      46             :    PUBLIC :: calculate_molecular_dipole
      47             : 
      48             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_dipoles'
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief maps wfc's to molecules and also prints molecular dipoles
      54             : !> \param qs_env the qs_env in which the qs_env lives
      55             : !> \param qs_loc_env ...
      56             : !> \param loc_print_key ...
      57             : !> \param molecule_set ...
      58             : ! **************************************************************************************************
      59          16 :    SUBROUTINE calculate_molecular_dipole(qs_env, qs_loc_env, loc_print_key, molecule_set)
      60             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      61             :       TYPE(qs_loc_env_type), INTENT(IN)                  :: qs_loc_env
      62             :       TYPE(section_vals_type), POINTER                   :: loc_print_key
      63             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
      64             : 
      65             :       COMPLEX(KIND=dp)                                   :: zeta
      66             :       COMPLEX(KIND=dp), DIMENSION(3)                     :: ggamma, zphase
      67             :       INTEGER                                            :: akind, first_atom, i, iatom, ikind, &
      68             :                                                             imol, imol_now, iounit, ispin, istate, &
      69             :                                                             j, natom, nkind, nmol, nspins, nstate, &
      70             :                                                             reference
      71             :       LOGICAL                                            :: do_berry, floating, ghost
      72             :       REAL(KIND=dp)                                      :: charge_tot, dipole(3), ria(3), theta, &
      73             :                                                             zeff, zwfc
      74             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: charge_set
      75             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dipole_set
      76             :       REAL(KIND=dp), DIMENSION(3)                        :: ci, gvec, rcc
      77          16 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
      78          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: center(:, :)
      79             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      80             :       TYPE(cell_type), POINTER                           :: cell
      81             :       TYPE(cp_logger_type), POINTER                      :: logger
      82             :       TYPE(dft_control_type), POINTER                    :: dft_control
      83             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
      84             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      85             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      86          16 :       TYPE(particle_type), POINTER                       :: particle_set(:)
      87          16 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      88             : 
      89          32 :       logger => cp_get_default_logger()
      90             : 
      91          16 :       CALL get_qs_env(qs_env, dft_control=dft_control)
      92          16 :       nspins = dft_control%nspins
      93             : 
      94             :       ! Setup reference point
      95          16 :       reference = section_get_ival(loc_print_key, keyword_name="MOLECULAR_DIPOLES%REFERENCE")
      96          16 :       CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%REF_POINT", r_vals=ref_point)
      97          16 :       CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%PERIODIC", l_val=do_berry)
      98             : 
      99          16 :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell)
     100          16 :       particle_set => qs_loc_env%particle_set
     101          16 :       para_env => qs_loc_env%para_env
     102          16 :       local_molecules => qs_loc_env%local_molecules
     103          16 :       nkind = SIZE(local_molecules%n_el)
     104          16 :       zwfc = 3.0_dp - REAL(nspins, KIND=dp)
     105             : 
     106          48 :       ALLOCATE (dipole_set(3, SIZE(molecule_set)))
     107          48 :       ALLOCATE (charge_set(SIZE(molecule_set)))
     108         168 :       dipole_set = 0.0_dp
     109          54 :       charge_set = 0.0_dp
     110             : 
     111          36 :       DO ispin = 1, nspins
     112          20 :          center => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
     113          20 :          nstate = SIZE(center, 2)
     114          82 :          DO ikind = 1, nkind ! loop over different molecules
     115          46 :             nmol = SIZE(local_molecules%list(ikind)%array)
     116          89 :             DO imol = 1, nmol ! all the molecules of the kind
     117          23 :                imol_now = local_molecules%list(ikind)%array(imol) ! index in the global array
     118          23 :                IF (.NOT. ASSOCIATED(molecule_set(imol_now)%lmi(ispin)%states)) CYCLE
     119          23 :                molecule_kind => molecule_set(imol_now)%molecule_kind
     120          23 :                first_atom = molecule_set(imol_now)%first_atom
     121          23 :                CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom)
     122             : 
     123             :                ! Get reference point for this molecule
     124             :                CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, &
     125             :                                         ref_point=ref_point, ifirst=first_atom, &
     126          23 :                                         ilast=first_atom + natom - 1)
     127             : 
     128          23 :                dipole = 0.0_dp
     129          23 :                IF (do_berry) THEN
     130          84 :                   rcc = pbc(rcc, cell)
     131             :                   ! Find out the total charge of the molecule
     132          67 :                   DO iatom = 1, natom
     133          46 :                      i = first_atom + iatom - 1
     134          46 :                      atomic_kind => particle_set(i)%atomic_kind
     135          46 :                      CALL get_atomic_kind(atomic_kind, kind_number=akind)
     136          46 :                      CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
     137         113 :                      IF (.NOT. ghost .AND. .NOT. floating) THEN
     138          46 :                         CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
     139          46 :                         charge_set(imol_now) = charge_set(imol_now) + zeff
     140             :                      END IF
     141             :                   END DO
     142             :                   ! Charges of the wfc involved
     143          83 :                   DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
     144          83 :                      charge_set(imol_now) = charge_set(imol_now) - zwfc
     145             :                   END DO
     146             : 
     147          21 :                   charge_tot = charge_set(imol_now)
     148         336 :                   ria = twopi*MATMUL(cell%h_inv, rcc)
     149          84 :                   zphase = CMPLX(COS(ria), SIN(ria), KIND=dp)**charge_tot
     150          84 :                   ggamma = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
     151             : 
     152             :                   ! Nuclear charges
     153          21 :                   IF (ispin == 1) THEN
     154          51 :                   DO iatom = 1, natom
     155          34 :                      i = first_atom + iatom - 1
     156          34 :                      atomic_kind => particle_set(i)%atomic_kind
     157          34 :                      CALL get_atomic_kind(atomic_kind, kind_number=akind)
     158          34 :                      CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
     159          85 :                      IF (.NOT. ghost .AND. .NOT. floating) THEN
     160          34 :                         CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
     161          34 :                         ria = pbc(particle_set(i)%r, cell)
     162         136 :                         DO j = 1, 3
     163         408 :                            gvec = twopi*cell%h_inv(j, :)
     164         408 :                            theta = SUM(ria(:)*gvec(:))
     165         102 :                            zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(zeff)
     166         136 :                            ggamma(j) = ggamma(j)*zeta
     167             :                         END DO
     168             :                      END IF
     169             :                   END DO
     170             :                   END IF
     171             : 
     172             :                   ! Charges of the wfc involved
     173          83 :                   DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
     174          62 :                      i = molecule_set(imol_now)%lmi(ispin)%states(istate)
     175          62 :                      ria = pbc(center(1:3, i), cell)
     176         269 :                      DO j = 1, 3
     177         744 :                         gvec = twopi*cell%h_inv(j, :)
     178         744 :                         theta = SUM(ria(:)*gvec(:))
     179         186 :                         zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-zwfc)
     180         248 :                         ggamma(j) = ggamma(j)*zeta
     181             :                      END DO
     182             :                   END DO
     183             : 
     184          84 :                   ggamma = ggamma*zphase
     185          84 :                   ci = AIMAG(LOG(ggamma))/twopi
     186         273 :                   dipole = MATMUL(cell%hmat, ci)
     187             :                ELSE
     188           2 :                   IF (ispin == 1) THEN
     189             :                      ! Nuclear charges
     190           8 :                      DO iatom = 1, natom
     191           6 :                         i = first_atom + iatom - 1
     192           6 :                         atomic_kind => particle_set(i)%atomic_kind
     193           6 :                         CALL get_atomic_kind(atomic_kind, kind_number=akind)
     194           6 :                         CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
     195          14 :                         IF (.NOT. ghost .AND. .NOT. floating) THEN
     196           6 :                            CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
     197          30 :                            ria = pbc(particle_set(i)%r, cell) - rcc
     198          24 :                            dipole = dipole + zeff*(ria - rcc)
     199           6 :                            charge_set(imol_now) = charge_set(imol_now) + zeff
     200             :                         END IF
     201             :                      END DO
     202             :                   END IF
     203             :                   ! Charges of the wfc involved
     204          10 :                   DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
     205           8 :                      i = molecule_set(imol_now)%lmi(ispin)%states(istate)
     206           8 :                      ria = pbc(center(1:3, i), cell)
     207          32 :                      dipole = dipole - zwfc*(ria - rcc)
     208          10 :                      charge_set(imol_now) = charge_set(imol_now) - zwfc
     209             :                   END DO
     210             :                END IF
     211         161 :                dipole_set(:, imol_now) = dipole_set(:, imol_now) + dipole ! a.u.
     212             :             END DO
     213             :          END DO
     214             :       END DO
     215          16 :       CALL para_env%sum(dipole_set)
     216          16 :       CALL para_env%sum(charge_set)
     217             : 
     218             :       iounit = cp_print_key_unit_nr(logger, loc_print_key, "MOLECULAR_DIPOLES", &
     219          16 :                                     extension=".MolDip", middle_name="MOLECULAR_DIPOLES")
     220          16 :       IF (iounit > 0) THEN
     221             :          WRITE (UNIT=iounit, FMT='(A80)') &
     222           8 :             "# molecule nr,      charge,           dipole vector,           dipole[Debye]"
     223          84 :          dipole_set(:, :) = dipole_set(:, :)*debye ! Debye
     224          27 :          DO I = 1, SIZE(dipole_set, 2)
     225          19 :             WRITE (UNIT=iounit, FMT='(T8,I6,T21,5F12.6)') I, charge_set(I), dipole_set(1:3, I), &
     226         103 :                SQRT(DOT_PRODUCT(dipole_set(1:3, I), dipole_set(1:3, I)))
     227             :          END DO
     228          84 :          WRITE (UNIT=iounit, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum  =', SUM(dipole_set)
     229             :       END IF
     230             :       CALL cp_print_key_finished_output(iounit, logger, loc_print_key, &
     231          16 :                                         "MOLECULAR_DIPOLES")
     232             : 
     233          16 :       DEALLOCATE (dipole_set, charge_set)
     234             : 
     235          32 :    END SUBROUTINE calculate_molecular_dipole
     236             :    !------------------------------------------------------------------------------
     237             : 
     238             : END MODULE molecular_dipoles
     239             : 

Generated by: LCOV version 1.15