39 #include "./base/base_uses.f90"
48 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'molecular_dipoles'
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(:)
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, &
71 LOGICAL :: do_berry, floating, ghost
72 REAL(kind=
dp) :: charge_tot, dipole(3), ria(3), theta, &
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
91 CALL get_qs_env(qs_env, dft_control=dft_control)
92 nspins = dft_control%nspins
95 reference =
section_get_ival(loc_print_key, keyword_name=
"MOLECULAR_DIPOLES%REFERENCE")
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)
106 ALLOCATE (dipole_set(3,
SIZE(molecule_set)))
107 ALLOCATE (charge_set(
SIZE(molecule_set)))
112 center => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
113 nstate =
SIZE(center, 2)
115 nmol =
SIZE(local_molecules%list(ikind)%array)
117 imol_now = local_molecules%list(ikind)%array(imol)
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
125 ref_point=ref_point, ifirst=first_atom, &
126 ilast=first_atom + natom - 1)
133 i = first_atom + iatom - 1
134 atomic_kind => particle_set(i)%atomic_kind
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
143 DO istate = 1,
SIZE(molecule_set(imol_now)%lmi(ispin)%states)
144 charge_set(imol_now) = charge_set(imol_now) - zwfc
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)
155 i = first_atom + iatom - 1
156 atomic_kind => particle_set(i)%atomic_kind
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)
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
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)
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
184 ggamma = ggamma*zphase
185 ci = aimag(log(ggamma))/
twopi
186 dipole = matmul(cell%hmat, ci)
191 i = first_atom + iatom - 1
192 atomic_kind => particle_set(i)%atomic_kind
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
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
211 dipole_set(:, imol_now) = dipole_set(:, imol_now) + dipole
215 CALL para_env%sum(dipole_set)
216 CALL para_env%sum(charge_set)
219 extension=
".MolDip", middle_name=
"MOLECULAR_DIPOLES")
221 WRITE (unit=iounit, fmt=
'(A80)') &
222 "# molecule nr, charge, dipole vector, dipole[Debye]"
223 dipole_set(:, :) = dipole_set(:, :)*
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)))
228 WRITE (unit=iounit, fmt=
"(T2,A,T61,E20.12)")
' DIPOLE : CheckSum =', sum(dipole_set)
233 DEALLOCATE (dipole_set, charge_set)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
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.
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...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
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>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public debye
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.
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...