(git:b279b6b)
atomic_charges.F
Go to the documentation of this file.
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  SUBROUTINE print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, &
47  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  CALL timeset(routinen, handle)
64 
65  IF (PRESENT(electronic_charges)) THEN
66  nspin = SIZE(electronic_charges, 2)
67  natom = SIZE(electronic_charges, 1)
68  ELSE
69  cpassert(PRESENT(atomic_charges))
70  natom = SIZE(atomic_charges, 1)
71  nspin = 0
72  END IF
73 
74  IF (scr > 0) THEN
75  IF (SIZE(particle_set) /= natom) THEN
76  cpabort("Unexpected number of atoms/charges")
77  END IF
78  WRITE (scr, '(T2,A)') title
79  SELECT CASE (nspin)
80  CASE (0, 1)
81  IF (title == "RESP charges:") THEN
82  WRITE (scr, '(A)') " Type | Atom | Charge"
83  ELSE
84  WRITE (scr, '(A)') " Atom | Charge"
85  END IF
86  CASE DEFAULT
87  WRITE (scr, '(A)') " Atom | Charge | Spin diff charge"
88  END SELECT
89  total_charge = 0.0_dp
90  !WRITE (scr, '(A)') ""
91  DO iatom = 1, natom
92  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
93  element_symbol=element_symbol, kind_number=ikind)
94  CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
95 
96  SELECT CASE (nspin)
97  CASE (0)
98  IF (title == "RESP charges:") THEN
99  WRITE (scr, '(T3,A4,2X,I6,A2,A2,F12.6)') "RESP", iatom, " ", element_symbol, atomic_charges(iatom)
100  total_charge = total_charge + atomic_charges(iatom)
101  ELSE
102  WRITE (scr, '(I6,A2,A2,F12.6)') iatom, " ", element_symbol, atomic_charges(iatom)
103  total_charge = total_charge + atomic_charges(iatom)
104  END IF
105  CASE (1)
106  WRITE (scr, '(I6,A2,A2,F12.6)') iatom, " ", element_symbol, zeff - electronic_charges(iatom, 1)
107  total_charge = total_charge + zeff - electronic_charges(iatom, 1)
108  CASE DEFAULT
109  WRITE (scr, '(I6,A2,A2,2F12.6)') iatom, " ", element_symbol, &
110  zeff - (electronic_charges(iatom, 1) + electronic_charges(iatom, 2)), &
111  (electronic_charges(iatom, 1) - electronic_charges(iatom, 2))
112  total_charge = total_charge + zeff - (electronic_charges(iatom, 1) + electronic_charges(iatom, 2))
113  END SELECT
114  END DO
115  IF (title == "RESP charges:") THEN
116  WRITE (scr, '(A,F10.6)') " Total ", total_charge
117  ELSE
118  WRITE (scr, '(A,F10.6)') " Total ", total_charge
119  END IF
120  WRITE (scr, '(A)') ""
121  END IF
122 
123  CALL timestop(handle)
124 
125  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  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  CALL timeset(routinen, handle)
154 
155  natom = 0
156  IF (PRESENT(charge)) THEN
157  natom = SIZE(charge)
158  END IF
159 
160  IF (scr > 0) THEN
161 
162  WRITE (scr, '(T2,A)') 'multipoles:'
163 
164  DO iatom = 1, natom
165  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
166  element_symbol=element_symbol, kind_number=ikind)
167  CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
168 
169  WRITE (scr, '(a,i5)') ' iatom= ', iatom
170  WRITE (scr, '(a,a2)') ' element_symbol= ', element_symbol
171  WRITE (scr, '(a,f20.10)') ' zeff= ', zeff
172 
173  WRITE (scr, '(a, f20.10)') 'charge = ', charge(iatom)
174  WRITE (scr, '(a,3f20.10)') 'dipole = ', dipole(:, iatom)
175  WRITE (scr, '(a)') 'quadrupole = '
176  DO i = 1, 3
177  WRITE (scr, '(3f20.10)') quadrupole(i, :, iatom)
178  END DO
179 
180  END DO
181  WRITE (scr, '(A)') ""
182  END IF
183 
184  CALL timestop(handle)
185 
186  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  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  IF (scr > 0) THEN
206  natom = SIZE(bond_orders, 1)
207  IF (SIZE(particle_set) /= natom) THEN
208  cpabort("Unexpected number of atoms/charges")
209  END IF
210  WRITE (scr, '(/,T2,A)') "Mayer Bond Orders"
211  WRITE (scr, '(T2,A,T20,A,T40,A)') " Type Atom 1 ", " Type Atom 2 ", " Bond Order "
212  DO iatom = 1, natom
213  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=el1)
214  DO jatom = iatom + 1, natom
215  CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, element_symbol=el2)
216  IF (bond_orders(iatom, jatom) > 0.1_dp) THEN
217  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  END SUBROUTINE print_bond_orders
224 
225 END MODULE atomic_charges
atomic_charges
simple routine to print charges for all atomic charge methods (currently mulliken,...
Definition: atomic_charges.F:14
atomic_charges::print_bond_orders
subroutine, public print_bond_orders(particle_set, scr, bond_orders)
Print Mayer bond orders.
Definition: atomic_charges.F:197
atomic_charges::print_atomic_charges
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Definition: atomic_charges.F:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451