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 simple routine to print charges for all atomic charge methods
      10             : !>      (currently mulliken, lowdin and ddapc)
      11             : !> \par History
      12             : !>      Joost VandeVondele [2006.03]
      13             : ! **************************************************************************************************
      14             : MODULE atomic_charges
      15             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16             :    USE kinds,                           ONLY: dp
      17             :    USE particle_types,                  ONLY: particle_type
      18             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      19             :                                               qs_kind_type
      20             : #include "./base/base_uses.f90"
      21             : 
      22             :    IMPLICIT NONE
      23             : 
      24             :    PRIVATE
      25             : 
      26             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atomic_charges'
      27             : 
      28             :    PUBLIC :: print_atomic_charges, print_bond_orders
      29             : 
      30             : CONTAINS
      31             : 
      32             : ! **************************************************************************************************
      33             : !> \brief generates a unified output format for atomic charges
      34             : !> \param particle_set ...
      35             : !> \param qs_kind_set ...
      36             : !> \param scr ...
      37             : !> \param title ...
      38             : !> \param electronic_charges (natom,nspin), the number of electrons of (so positive) per spin
      39             : !>                            if (nspin==1) it is the sum of alpha and beta electrons
      40             : !> \param atomic_charges truly the atomic charge (taking Z into account, atoms negative, no spin)
      41             : !> \par History
      42             : !>      03.2006 created [Joost VandeVondele]
      43             : !> \note
      44             : !>      charges are computed per spin in the LSD case
      45             : ! **************************************************************************************************
      46        2780 :    SUBROUTINE print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, &
      47        1390 :                                    atomic_charges)
      48             : 
      49             :       TYPE(particle_type), DIMENSION(:), INTENT(IN)      :: particle_set
      50             :       TYPE(qs_kind_type), DIMENSION(:), INTENT(IN)       :: qs_kind_set
      51             :       INTEGER, INTENT(IN)                                :: scr
      52             :       CHARACTER(LEN=*), INTENT(IN)                       :: title
      53             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
      54             :          OPTIONAL                                        :: electronic_charges
      55             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: atomic_charges
      56             : 
      57             :       CHARACTER(len=*), PARAMETER :: routineN = 'print_atomic_charges'
      58             : 
      59             :       CHARACTER(LEN=2)                                   :: element_symbol
      60             :       INTEGER                                            :: handle, iatom, ikind, natom, nspin
      61             :       REAL(KIND=dp)                                      :: total_charge, zeff
      62             : 
      63        1390 :       CALL timeset(routineN, handle)
      64             : 
      65        1390 :       IF (PRESENT(electronic_charges)) THEN
      66          34 :          nspin = SIZE(electronic_charges, 2)
      67          34 :          natom = SIZE(electronic_charges, 1)
      68             :       ELSE
      69        1356 :          CPASSERT(PRESENT(atomic_charges))
      70        1356 :          natom = SIZE(atomic_charges, 1)
      71        1356 :          nspin = 0
      72             :       END IF
      73             : 
      74        1390 :       IF (scr > 0) THEN
      75         325 :          IF (SIZE(particle_set) /= natom) THEN
      76           0 :             CPABORT("Unexpected number of atoms/charges")
      77             :          END IF
      78         325 :          WRITE (scr, '(T2,A)') title
      79         324 :          SELECT CASE (nspin)
      80             :          CASE (0, 1)
      81         324 :             IF (title == "RESP charges:") THEN
      82           7 :                WRITE (scr, '(A)') "  Type |   Atom   |    Charge"
      83             :             ELSE
      84         317 :                WRITE (scr, '(A)') "  Atom     |    Charge"
      85             :             END IF
      86             :          CASE DEFAULT
      87         325 :             WRITE (scr, '(A)') "  Atom     |    Charge | Spin diff charge"
      88             :          END SELECT
      89         325 :          total_charge = 0.0_dp
      90             :          !WRITE (scr, '(A)') ""
      91        1221 :          DO iatom = 1, natom
      92             :             CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
      93         896 :                                  element_symbol=element_symbol, kind_number=ikind)
      94         896 :             CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
      95             : 
      96        1221 :             SELECT CASE (nspin)
      97             :             CASE (0)
      98         835 :                IF (title == "RESP charges:") THEN
      99          45 :                   WRITE (scr, '(T3,A4,2X,I6,A2,A2,F12.6)') "RESP", iatom, "  ", element_symbol, atomic_charges(iatom)
     100          45 :                   total_charge = total_charge + atomic_charges(iatom)
     101             :                ELSE
     102         790 :                   WRITE (scr, '(I6,A2,A2,F12.6)') iatom, "  ", element_symbol, atomic_charges(iatom)
     103         790 :                   total_charge = total_charge + atomic_charges(iatom)
     104             :                END IF
     105             :             CASE (1)
     106          58 :                WRITE (scr, '(I6,A2,A2,F12.6)') iatom, "  ", element_symbol, zeff - electronic_charges(iatom, 1)
     107          58 :                total_charge = total_charge + zeff - electronic_charges(iatom, 1)
     108             :             CASE DEFAULT
     109           3 :                WRITE (scr, '(I6,A2,A2,2F12.6)') iatom, "  ", element_symbol, &
     110           3 :                   zeff - (electronic_charges(iatom, 1) + electronic_charges(iatom, 2)), &
     111           6 :                   (electronic_charges(iatom, 1) - electronic_charges(iatom, 2))
     112         899 :                total_charge = total_charge + zeff - (electronic_charges(iatom, 1) + electronic_charges(iatom, 2))
     113             :             END SELECT
     114             :          END DO
     115         325 :          IF (title == "RESP charges:") THEN
     116           7 :             WRITE (scr, '(A,F10.6)') "  Total             ", total_charge
     117             :          ELSE
     118         318 :             WRITE (scr, '(A,F10.6)') "  Total     ", total_charge
     119             :          END IF
     120         325 :          WRITE (scr, '(A)') ""
     121             :       END IF
     122             : 
     123        1390 :       CALL timestop(handle)
     124             : 
     125        1390 :    END SUBROUTINE print_atomic_charges
     126             : 
     127             : ! **************************************************************************************************
     128             : !> \brief ...
     129             : !> \param particle_set ...
     130             : !> \param qs_kind_set ...
     131             : !> \param scr ...
     132             : !> \param charge ...
     133             : !> \param dipole ...
     134             : !> \param quadrupole ...
     135             : ! **************************************************************************************************
     136           0 :    SUBROUTINE print_multipoles(particle_set, qs_kind_set, scr, charge, dipole, quadrupole)
     137             : 
     138             :       TYPE(particle_type), DIMENSION(:), INTENT(IN)      :: particle_set
     139             :       TYPE(qs_kind_type), DIMENSION(:), INTENT(IN)       :: qs_kind_set
     140             :       INTEGER, INTENT(IN)                                :: scr
     141             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: charge
     142             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     143             :          OPTIONAL                                        :: dipole
     144             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
     145             :          OPTIONAL                                        :: quadrupole
     146             : 
     147             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'print_multipoles'
     148             : 
     149             :       CHARACTER(LEN=2)                                   :: element_symbol
     150             :       INTEGER                                            :: handle, i, iatom, ikind, natom
     151             :       REAL(KIND=dp)                                      :: zeff
     152             : 
     153           0 :       CALL timeset(routineN, handle)
     154             : 
     155           0 :       natom = 0
     156           0 :       IF (PRESENT(charge)) THEN
     157           0 :          natom = SIZE(charge)
     158             :       END IF
     159             : 
     160           0 :       IF (scr > 0) THEN
     161             : 
     162           0 :          WRITE (scr, '(T2,A)') 'multipoles:'
     163             : 
     164           0 :          DO iatom = 1, natom
     165             :             CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     166           0 :                                  element_symbol=element_symbol, kind_number=ikind)
     167           0 :             CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     168             : 
     169           0 :             WRITE (scr, '(a,i5)') ' iatom= ', iatom
     170           0 :             WRITE (scr, '(a,a2)') ' element_symbol= ', element_symbol
     171           0 :             WRITE (scr, '(a,f20.10)') ' zeff= ', zeff
     172             : 
     173           0 :             WRITE (scr, '(a, f20.10)') 'charge =     ', charge(iatom)
     174           0 :             WRITE (scr, '(a,3f20.10)') 'dipole =     ', dipole(:, iatom)
     175           0 :             WRITE (scr, '(a)') 'quadrupole = '
     176           0 :             DO i = 1, 3
     177           0 :                WRITE (scr, '(3f20.10)') quadrupole(i, :, iatom)
     178             :             END DO
     179             : 
     180             :          END DO
     181           0 :          WRITE (scr, '(A)') ""
     182             :       END IF
     183             : 
     184           0 :       CALL timestop(handle)
     185             : 
     186           0 :    END SUBROUTINE print_multipoles
     187             : 
     188             : ! **************************************************************************************************
     189             : !> \brief Print Mayer bond orders
     190             : !> \param particle_set ...
     191             : !> \param scr ...
     192             : !> \param bond_orders (natom,natom)
     193             : !> \par History
     194             : !>      12.2016 created [JGH]
     195             : ! **************************************************************************************************
     196           0 :    SUBROUTINE print_bond_orders(particle_set, scr, bond_orders)
     197             : 
     198             :       TYPE(particle_type), DIMENSION(:), INTENT(IN)      :: particle_set
     199             :       INTEGER, INTENT(IN)                                :: scr
     200             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: bond_orders
     201             : 
     202             :       CHARACTER(LEN=2)                                   :: el1, el2
     203             :       INTEGER                                            :: iatom, jatom, natom
     204             : 
     205           0 :       IF (scr > 0) THEN
     206           0 :          natom = SIZE(bond_orders, 1)
     207           0 :          IF (SIZE(particle_set) /= natom) THEN
     208           0 :             CPABORT("Unexpected number of atoms/charges")
     209             :          END IF
     210           0 :          WRITE (scr, '(/,T2,A)') "Mayer Bond Orders"
     211           0 :          WRITE (scr, '(T2,A,T20,A,T40,A)') "  Type  Atom 1  ", "  Type  Atom 2  ", " Bond Order "
     212           0 :          DO iatom = 1, natom
     213           0 :             CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=el1)
     214           0 :             DO jatom = iatom + 1, natom
     215           0 :                CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, element_symbol=el2)
     216           0 :                IF (bond_orders(iatom, jatom) > 0.1_dp) THEN
     217           0 :                   WRITE (scr, '(T6,A2,I6,T24,A2,I6,T40,F12.6)') el1, iatom, el2, jatom, bond_orders(iatom, jatom)
     218             :                END IF
     219             :             END DO
     220             :          END DO
     221             :       END IF
     222             : 
     223           0 :    END SUBROUTINE print_bond_orders
     224             : 
     225             : END MODULE atomic_charges