(git:0de0cc2)
molecular_dipoles.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 Set of routines handling the localization for molecular properties
10 ! **************************************************************************************************
12  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  USE cell_types, ONLY: cell_type,&
15  pbc
16  USE cp_control_types, ONLY: dft_control_type
18  cp_logger_type
21  USE distribution_1d_types, ONLY: distribution_1d_type
23  section_vals_type,&
25  USE kinds, ONLY: dp
26  USE mathconstants, ONLY: twopi
27  USE message_passing, ONLY: mp_para_env_type
29  molecule_kind_type
30  USE molecule_types, ONLY: molecule_type
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 ***
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  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  REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
78  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  TYPE(particle_type), POINTER :: particle_set(:)
87  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
88 
89  logger => cp_get_default_logger()
90 
91  CALL get_qs_env(qs_env, dft_control=dft_control)
92  nspins = dft_control%nspins
93 
94  ! Setup reference point
95  reference = section_get_ival(loc_print_key, keyword_name="MOLECULAR_DIPOLES%REFERENCE")
96  CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%REF_POINT", r_vals=ref_point)
97  CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%PERIODIC", l_val=do_berry)
98 
99  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell)
100  particle_set => qs_loc_env%particle_set
101  para_env => qs_loc_env%para_env
102  local_molecules => qs_loc_env%local_molecules
103  nkind = SIZE(local_molecules%n_el)
104  zwfc = 3.0_dp - real(nspins, kind=dp)
105 
106  ALLOCATE (dipole_set(3, SIZE(molecule_set)))
107  ALLOCATE (charge_set(SIZE(molecule_set)))
108  dipole_set = 0.0_dp
109  charge_set = 0.0_dp
110 
111  DO ispin = 1, nspins
112  center => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
113  nstate = SIZE(center, 2)
114  DO ikind = 1, nkind ! loop over different molecules
115  nmol = SIZE(local_molecules%list(ikind)%array)
116  DO imol = 1, nmol ! all the molecules of the kind
117  imol_now = local_molecules%list(ikind)%array(imol) ! index in the global array
118  IF (.NOT. ASSOCIATED(molecule_set(imol_now)%lmi(ispin)%states)) cycle
119  molecule_kind => molecule_set(imol_now)%molecule_kind
120  first_atom = molecule_set(imol_now)%first_atom
121  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  ilast=first_atom + natom - 1)
127 
128  dipole = 0.0_dp
129  IF (do_berry) THEN
130  rcc = pbc(rcc, cell)
131  ! Find out the total charge of the molecule
132  DO iatom = 1, natom
133  i = first_atom + iatom - 1
134  atomic_kind => particle_set(i)%atomic_kind
135  CALL get_atomic_kind(atomic_kind, kind_number=akind)
136  CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
137  IF (.NOT. ghost .AND. .NOT. floating) THEN
138  CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
139  charge_set(imol_now) = charge_set(imol_now) + zeff
140  END IF
141  END DO
142  ! Charges of the wfc involved
143  DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
144  charge_set(imol_now) = charge_set(imol_now) - zwfc
145  END DO
146 
147  charge_tot = charge_set(imol_now)
148  ria = twopi*matmul(cell%h_inv, rcc)
149  zphase = cmplx(cos(ria), sin(ria), kind=dp)**charge_tot
150  ggamma = cmplx(1.0_dp, 0.0_dp, kind=dp)
151 
152  ! Nuclear charges
153  IF (ispin == 1) THEN
154  DO iatom = 1, natom
155  i = first_atom + iatom - 1
156  atomic_kind => particle_set(i)%atomic_kind
157  CALL get_atomic_kind(atomic_kind, kind_number=akind)
158  CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
159  IF (.NOT. ghost .AND. .NOT. floating) THEN
160  CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
161  ria = pbc(particle_set(i)%r, cell)
162  DO j = 1, 3
163  gvec = twopi*cell%h_inv(j, :)
164  theta = sum(ria(:)*gvec(:))
165  zeta = cmplx(cos(theta), sin(theta), kind=dp)**(zeff)
166  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  DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
174  i = molecule_set(imol_now)%lmi(ispin)%states(istate)
175  ria = pbc(center(1:3, i), cell)
176  DO j = 1, 3
177  gvec = twopi*cell%h_inv(j, :)
178  theta = sum(ria(:)*gvec(:))
179  zeta = cmplx(cos(theta), sin(theta), kind=dp)**(-zwfc)
180  ggamma(j) = ggamma(j)*zeta
181  END DO
182  END DO
183 
184  ggamma = ggamma*zphase
185  ci = aimag(log(ggamma))/twopi
186  dipole = matmul(cell%hmat, ci)
187  ELSE
188  IF (ispin == 1) THEN
189  ! Nuclear charges
190  DO iatom = 1, natom
191  i = first_atom + iatom - 1
192  atomic_kind => particle_set(i)%atomic_kind
193  CALL get_atomic_kind(atomic_kind, kind_number=akind)
194  CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
195  IF (.NOT. ghost .AND. .NOT. floating) THEN
196  CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
197  ria = pbc(particle_set(i)%r, cell) - rcc
198  dipole = dipole + zeff*(ria - rcc)
199  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  DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states)
205  i = molecule_set(imol_now)%lmi(ispin)%states(istate)
206  ria = pbc(center(1:3, i), cell)
207  dipole = dipole - zwfc*(ria - rcc)
208  charge_set(imol_now) = charge_set(imol_now) - zwfc
209  END DO
210  END IF
211  dipole_set(:, imol_now) = dipole_set(:, imol_now) + dipole ! a.u.
212  END DO
213  END DO
214  END DO
215  CALL para_env%sum(dipole_set)
216  CALL para_env%sum(charge_set)
217 
218  iounit = cp_print_key_unit_nr(logger, loc_print_key, "MOLECULAR_DIPOLES", &
219  extension=".MolDip", middle_name="MOLECULAR_DIPOLES")
220  IF (iounit > 0) THEN
221  WRITE (unit=iounit, fmt='(A80)') &
222  "# molecule nr, charge, dipole vector, dipole[Debye]"
223  dipole_set(:, :) = dipole_set(:, :)*debye ! Debye
224  DO i = 1, SIZE(dipole_set, 2)
225  WRITE (unit=iounit, fmt='(T8,I6,T21,5F12.6)') i, charge_set(i), dipole_set(1:3, i), &
226  sqrt(dot_product(dipole_set(1:3, i), dipole_set(1:3, i)))
227  END DO
228  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  "MOLECULAR_DIPOLES")
232 
233  DEALLOCATE (dipole_set, charge_set)
234 
235  END SUBROUTINE calculate_molecular_dipole
236  !------------------------------------------------------------------------------
237 
238 END MODULE molecular_dipoles
239 
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 ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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,...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
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)
...
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
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
Set of routines handling the localization for molecular properties.
subroutine, public calculate_molecular_dipole(qs_env, qs_loc_env, loc_print_key, molecule_set)
maps wfc's to molecules and also prints molecular dipoles
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Define the data structure for the molecule information.
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.
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