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 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 3624 : SUBROUTINE print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, &
47 1812 : 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 1812 : CALL timeset(routineN, handle)
64 :
65 1812 : IF (PRESENT(electronic_charges)) THEN
66 34 : nspin = SIZE(electronic_charges, 2)
67 34 : natom = SIZE(electronic_charges, 1)
68 : ELSE
69 1778 : CPASSERT(PRESENT(atomic_charges))
70 1778 : natom = SIZE(atomic_charges, 1)
71 1778 : nspin = 0
72 : END IF
73 :
74 1812 : IF (scr > 0) THEN
75 493 : IF (SIZE(particle_set) /= natom) THEN
76 0 : CPABORT("Unexpected number of atoms/charges")
77 : END IF
78 493 : WRITE (scr, '(T2,A)') title
79 492 : SELECT CASE (nspin)
80 : CASE (0, 1)
81 492 : IF (title == "RESP charges:") THEN
82 7 : WRITE (scr, '(A)') " Type | Atom | Charge"
83 : ELSE
84 485 : WRITE (scr, '(A)') " Atom | Charge"
85 : END IF
86 : CASE DEFAULT
87 493 : WRITE (scr, '(A)') " Atom | Charge | Spin diff charge"
88 : END SELECT
89 493 : total_charge = 0.0_dp
90 : !WRITE (scr, '(A)') ""
91 1788 : DO iatom = 1, natom
92 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
93 1295 : element_symbol=element_symbol, kind_number=ikind)
94 1295 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
95 :
96 1788 : SELECT CASE (nspin)
97 : CASE (0)
98 1234 : 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 1189 : WRITE (scr, '(I6,A2,A2,F12.6)') iatom, " ", element_symbol, atomic_charges(iatom)
103 1189 : 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 1298 : total_charge = total_charge + zeff - (electronic_charges(iatom, 1) + electronic_charges(iatom, 2))
113 : END SELECT
114 : END DO
115 493 : IF (title == "RESP charges:") THEN
116 7 : WRITE (scr, '(A,F10.6)') " Total ", total_charge
117 : ELSE
118 486 : WRITE (scr, '(A,F10.6)') " Total ", total_charge
119 : END IF
120 493 : WRITE (scr, '(A)') ""
121 : END IF
122 :
123 1812 : CALL timestop(handle)
124 :
125 1812 : 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
|
