19 USE dbcsr_api,
ONLY: dbcsr_add,&
21 dbcsr_iterator_blocks_left,&
22 dbcsr_iterator_next_block,&
23 dbcsr_iterator_start,&
29 ewald_environment_type
49 #include "./base/base_uses.f90"
55 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qmmm_tb_coulomb'
72 calculate_forces, just_energy)
74 TYPE(qs_environment_type),
INTENT(IN) :: qs_env
75 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
76 TYPE(qs_rho_type),
POINTER :: rho
77 REAL(
dp),
DIMENSION(:) :: mcharge
78 TYPE(qs_energy_type),
POINTER :: energy
79 LOGICAL,
INTENT(in) :: calculate_forces, just_energy
81 CHARACTER(len=*),
PARAMETER :: routinen =
'build_tb_coulomb_qmqm'
83 INTEGER :: atom_i, atom_j, blk, ewald_type, handle, &
84 i, ia, iatom, ikind, jatom, jkind, &
86 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
87 INTEGER,
DIMENSION(3) :: periodic
88 LOGICAL :: found, use_virial
89 REAL(kind=
dp) :: alpha, deth, dfr, dr, fi, fr, gmij
90 REAL(kind=
dp),
DIMENSION(3) :: fij, rij
91 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dsblock, gmcharge, ksblock, ksblock_2, &
93 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
94 TYPE(atprop_type),
POINTER :: atprop
95 TYPE(cell_type),
POINTER :: cell, mm_cell
96 TYPE(dbcsr_iterator_type) :: iter
97 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p, matrix_s
98 TYPE(dft_control_type),
POINTER :: dft_control
99 TYPE(distribution_1d_type),
POINTER :: local_particles
100 TYPE(ewald_environment_type),
POINTER :: ewald_env
101 TYPE(ewald_pw_type),
POINTER :: ewald_pw
102 TYPE(mp_para_env_type),
POINTER :: para_env
103 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
104 TYPE(qmmm_env_qm_type),
POINTER :: qmmm_env_qm
105 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
106 TYPE(virial_type),
POINTER :: virial
108 CALL timeset(routinen, handle)
110 NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
114 IF (calculate_forces)
THEN
120 natom =
SIZE(mcharge)
121 ALLOCATE (gmcharge(natom, nmat))
125 particle_set=particle_set, &
129 dft_control=dft_control)
131 IF (calculate_forces)
THEN
132 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
137 cpassert(.NOT. atprop%energy)
139 cpassert(.NOT. use_virial)
140 qmmm_env_qm => qs_env%qmmm_env_qm
141 ewald_env => qmmm_env_qm%ewald_env
142 ewald_pw => qmmm_env_qm%ewald_pw
143 CALL get_qs_env(qs_env=qs_env, super_cell=mm_cell)
144 CALL get_cell(cell=mm_cell, periodic=periodic, deth=deth)
145 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
149 atomic_kind_set=atomic_kind_set, &
150 local_particles=local_particles, &
151 force=force, para_env=para_env)
152 DO ikind = 1,
SIZE(local_particles%n_el)
153 DO ia = 1, local_particles%n_el(ikind)
154 iatom = local_particles%list(ikind)%array(ia)
155 DO jatom = 1, iatom - 1
156 rij = particle_set(iatom)%r - particle_set(jatom)%r
159 dr = sqrt(sum(rij(:)**2))
161 gmcharge(iatom, 1) = gmcharge(iatom, 1) - mcharge(jatom)/dr
162 gmcharge(jatom, 1) = gmcharge(jatom, 1) - mcharge(iatom)/dr
164 gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)/dr**3
165 gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)/dr**3
168 fr = erfc(alpha*dr)/dr
169 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
170 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
172 dfr = -2._dp*alpha*exp(-alpha*alpha*dr*dr)*
oorootpi/dr - fr/dr
175 gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
176 gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
183 SELECT CASE (ewald_type)
185 cpabort(
"Invalid Ewald type")
187 cpabort(
"Not allowed with DFTB")
189 cpabort(
"Standard Ewald not implemented in DFTB")
191 cpabort(
"PME not implemented in DFTB")
194 gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
197 CALL para_env%sum(gmcharge(:, 1))
200 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*
oorootpi*mcharge(:)
201 IF (any(periodic(:) == 1))
THEN
202 gmcharge(:, 1) = gmcharge(:, 1) -
pi/alpha**2/deth
205 energy%qmmm_el = energy%qmmm_el + 0.5_dp*sum(mcharge(:)*gmcharge(:, 1))
207 IF (calculate_forces)
THEN
210 atom_of_kind=atom_of_kind)
213 IF (.NOT. just_energy)
THEN
214 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
217 IF (calculate_forces .AND.
SIZE(matrix_p) == 2)
THEN
218 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
219 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
222 CALL dbcsr_iterator_start(iter, ks_matrix(1, 1)%matrix)
223 DO WHILE (dbcsr_iterator_blocks_left(iter))
224 CALL dbcsr_iterator_next_block(iter, iatom, jatom, ksblock, blk)
225 NULLIFY (sblock, ksblock_2)
226 IF (
SIZE(ks_matrix, 1) > 1)
THEN
227 CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, &
228 row=iatom, col=jatom, block=ksblock_2, found=found)
230 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
231 row=iatom, col=jatom, block=sblock, found=found)
232 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
233 ksblock = ksblock - gmij*sblock
234 IF (
SIZE(ks_matrix, 1) > 1) ksblock_2 = ksblock_2 - gmij*sblock
235 IF (calculate_forces)
THEN
236 ikind = kind_of(iatom)
237 atom_i = atom_of_kind(iatom)
238 jkind = kind_of(jatom)
239 atom_j = atom_of_kind(jatom)
241 CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
242 row=iatom, col=jatom, block=pblock, found=found)
245 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
246 row=iatom, col=jatom, block=dsblock, found=found)
247 fi = -2.0_dp*gmij*sum(pblock*dsblock)
248 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
249 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
254 CALL dbcsr_iterator_stop(iter)
255 IF (calculate_forces .AND.
SIZE(matrix_p) == 2)
THEN
256 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
257 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
261 DEALLOCATE (gmcharge)
263 CALL timestop(handle)
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Holds information on atomic properties.
Handles all functions related to the CELL.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Calculation of Ewald contributions in DFTB.
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Define the data structure for the particle information.
functions related to the poisson solver on regular grids
integer, parameter, public do_ewald_pme
integer, parameter, public do_ewald_ewald
integer, parameter, public do_ewald_none
integer, parameter, public do_ewald_spme
Calculation of QMMM Coulomb contributions in TB.
subroutine, public build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
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.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...