(git:0de0cc2)
qs_loc_dipole.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
10 ! **************************************************************************************************
13  USE cell_types, ONLY: cell_type,&
14  pbc
15  USE cp_control_types, ONLY: dft_control_type
16  USE cp_log_handling, ONLY: cp_logger_type
18  cp_p_file,&
23  get_results,&
24  put_results
25  USE cp_result_types, ONLY: cp_result_type
28  section_vals_type,&
30  USE kinds, ONLY: default_string_length,&
31  dp
32  USE mathconstants, ONLY: twopi
34  USE particle_types, ONLY: particle_type
35  USE physcon, ONLY: debye
36  USE qs_environment_types, ONLY: get_qs_env,&
37  qs_environment_type
38  USE qs_kind_types, ONLY: get_qs_kind,&
39  qs_kind_type
40  USE qs_loc_types, ONLY: qs_loc_env_type
41 #include "./base/base_uses.f90"
42 
43  IMPLICIT NONE
44  PRIVATE
45 
46  ! Global parameters
47  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_dipole'
48  PUBLIC :: loc_dipole
49 
50 ! **************************************************************************************************
51 
52 CONTAINS
53 
54 ! **************************************************************************************************
55 !> \brief Computes and prints the Dipole (using localized charges)
56 !> \param input ...
57 !> \param dft_control ...
58 !> \param qs_loc_env ...
59 !> \param logger ...
60 !> \param qs_env the qs_env in which the qs_env lives
61 ! **************************************************************************************************
62  SUBROUTINE loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
63  TYPE(section_vals_type), POINTER :: input
64  TYPE(dft_control_type), POINTER :: dft_control
65  TYPE(qs_loc_env_type), POINTER :: qs_loc_env
66  TYPE(cp_logger_type), POINTER :: logger
67  TYPE(qs_environment_type), POINTER :: qs_env
68 
69  CHARACTER(len=*), PARAMETER :: routinen = 'loc_dipole'
70 
71  CHARACTER(LEN=default_string_length) :: description, descriptionthisdip, iter
72  COMPLEX(KIND=dp) :: zeta
73  COMPLEX(KIND=dp), DIMENSION(3) :: ggamma, zphase
74  INTEGER :: handle, i, ikind, ispins, j, n_rep, &
75  reference, unit_nr
76  LOGICAL :: do_berry, first_time, floating, ghost
77  REAL(kind=dp) :: charge_tot, theta, zeff, zwfc
78  REAL(kind=dp), DIMENSION(3) :: ci, dipole, dipole_old, gvec, rcc, ria
79  REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
80  TYPE(cell_type), POINTER :: cell
81  TYPE(cp_result_type), POINTER :: results
82  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
83  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
84  TYPE(section_vals_type), POINTER :: print_key
85 
86  CALL timeset(routinen, handle)
87 
88  print_key => section_vals_get_subs_vals(input, "DFT%LOCALIZE%PRINT%TOTAL_DIPOLE")
89  IF (btest(cp_print_key_should_output(logger%iter_info, print_key, first_time=first_time) &
90  , cp_p_file)) THEN
91  NULLIFY (cell, particle_set, qs_kind_set, ref_point, results)
92  CALL get_qs_env(qs_env=qs_env, &
93  cell=cell, &
94  particle_set=particle_set, &
95  qs_kind_set=qs_kind_set, &
96  results=results)
97 
98  reference = section_get_ival(print_key, keyword_name="REFERENCE")
99  CALL section_vals_val_get(print_key, "REF_POINT", r_vals=ref_point)
100  CALL section_vals_val_get(print_key, "PERIODIC", l_val=do_berry)
101  description = '[DIPOLE]'
102  descriptionthisdip = '[TOTAL_DIPOLE]'
103  CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
104 
105  dipole = 0.0_dp
106  IF (do_berry) THEN
107  rcc = pbc(rcc, cell)
108  charge_tot = real(dft_control%charge, kind=dp)
109  ria = twopi*matmul(cell%h_inv, rcc)
110  zphase = cmplx(cos(ria), sin(ria), kind=dp)**charge_tot
111  ggamma = cmplx(1.0_dp, 0.0_dp, kind=dp)
112 
113  ! Nuclear charges
114  DO i = 1, SIZE(particle_set)
115  CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
116  CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
117  IF (.NOT. ghost .AND. .NOT. floating) THEN
118  CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
119  ria = pbc(particle_set(i)%r, cell)
120  DO j = 1, 3
121  gvec = twopi*cell%h_inv(j, :)
122  theta = sum(ria(:)*gvec(:))
123  zeta = cmplx(cos(theta), sin(theta), kind=dp)**(zeff)
124  ggamma(j) = ggamma(j)*zeta
125  END DO
126  END IF
127  END DO
128 
129  ! Charges of the wfc involved
130  ! Warning, this assumes the same occupation for all states
131  zwfc = 3.0_dp - real(dft_control%nspins, dp)
132 
133  DO ispins = 1, dft_control%nspins
134  DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
135  ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
136  DO j = 1, 3
137  gvec = twopi*cell%h_inv(j, :)
138  theta = sum(ria(:)*gvec(:))
139  zeta = cmplx(cos(theta), sin(theta), kind=dp)**(-zwfc)
140  ggamma(j) = ggamma(j)*zeta
141  END DO
142  END DO
143  END DO
144  ggamma = ggamma*zphase
145  ci = aimag(log(ggamma))/twopi
146  dipole = matmul(cell%hmat, ci)
147  ELSE
148  ! Charges of the atoms involved
149  DO i = 1, SIZE(particle_set)
150  CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
151  CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
152  IF (.NOT. ghost .AND. .NOT. floating) THEN
153  CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
154  ria = pbc(particle_set(i)%r, cell)
155  dipole = dipole + zeff*(ria - rcc)
156  END IF
157  END DO
158 
159  ! Charges of the wfc involved
160  ! Warning, this assumes the same occupation for all states
161  zwfc = 3.0_dp - real(dft_control%nspins, dp)
162 
163  DO ispins = 1, dft_control%nspins
164  DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
165  ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
166  dipole = dipole - zwfc*(ria - rcc)
167  END DO
168  END DO
169  END IF
170 
171  ! Print and possibly store results
172  unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".Dipole", &
173  middle_name="TOTAL_DIPOLE")
174  IF (unit_nr > 0) THEN
175  IF (first_time) THEN
176  WRITE (unit=unit_nr, fmt="(A,T31,A,T88,A,T136,A)") &
177  "# iter_level", "dipole(x,y,z)[atomic units]", &
178  "dipole(x,y,z)[debye]", &
179  "delta_dipole(x,y,z)[atomic units]"
180  END IF
181  iter = cp_iter_string(logger%iter_info)
182  CALL get_results(results, descriptionthisdip, n_rep=n_rep)
183  IF (n_rep == 0) THEN
184  dipole_old = 0._dp
185  ELSE
186  CALL get_results(results, descriptionthisdip, dipole_old, nval=n_rep)
187  END IF
188  IF (do_berry) THEN
189  WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
190  iter(1:15), dipole, dipole*debye, pbc(dipole - dipole_old, cell)
191  ELSE
192  WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
193  iter(1:15), dipole, dipole*debye, (dipole - dipole_old)
194  END IF
195  END IF
196  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
197  CALL cp_results_erase(results, description)
198  CALL put_results(results, description, dipole)
199  CALL cp_results_erase(results, descriptionthisdip)
200  CALL put_results(results, descriptionthisdip, dipole)
201  END IF
202 
203  CALL timestop(handle)
204 
205  END SUBROUTINE loc_dipole
206 
207 END MODULE qs_loc_dipole
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
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.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
character(len=default_string_length) function, public cp_iter_string(iter_info, print_key, for_file)
returns the iteration string, a string that is useful to create unique filenames (once you trim it)
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
set of type/routines to handle the storage of results in force_envs
objects that represent the structure of input sections and the data contained in an input section
integer function, public section_get_ival(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public twopi
Calculates the moment integrals <a|r^m|b>
Definition: moments_utils.F:15
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition: moments_utils.F:61
Define the data structure for the particle information.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public debye
Definition: physcon.F:201
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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.
subroutine, public loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
Computes and prints the Dipole (using localized charges)
Definition: qs_loc_dipole.F:63
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition: qs_loc_types.F:25