27 USE dbcsr_api,
ONLY: &
28 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
29 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
30 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_type, &
31 dbcsr_type_no_symmetry
48 #include "./base/base_uses.f90"
56 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rt_propagation_forces'
71 TYPE(qs_environment_type),
POINTER :: qs_env
73 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_c_mat_force'
74 REAL(kind=
dp),
PARAMETER ::
one = 1.0_dp,
zero = 0.0_dp
76 INTEGER :: handle, i, im, ispin, re
77 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
78 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
79 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: c_mat, rho_ao, rho_ao_im, rho_new, &
80 s_der, sinvb, sinvh, sinvh_imag
81 TYPE(dbcsr_type),
POINTER :: s_inv, tmp
82 TYPE(dft_control_type),
POINTER :: dft_control
83 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
84 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
85 TYPE(qs_rho_type),
POINTER :: rho
86 TYPE(rt_prop_type),
POINTER :: rtp
87 TYPE(rtp_control_type),
POINTER :: rtp_control
89 CALL timeset(routinen, handle)
90 NULLIFY (rtp, particle_set, atomic_kind_set, dft_control)
95 particle_set=particle_set, &
96 atomic_kind_set=atomic_kind_set, &
98 dft_control=dft_control)
100 rtp_control => dft_control%rtp_control
101 CALL get_rtp(rtp=rtp, c_mat=c_mat, s_der=s_der, s_inv=s_inv, &
102 sinvh=sinvh, sinvb=sinvb)
104 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
108 CALL dbcsr_create(tmp, template=sinvb(1)%matrix)
110 IF (rtp%linear_scaling)
THEN
111 CALL get_rtp(rtp=rtp, rho_new=rho_new)
113 CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
117 IF (rtp%propagate_complex_ks)
CALL get_rtp(rtp=rtp, sinvh_imag=sinvh_imag)
119 DO ispin = 1,
SIZE(sinvh)
122 IF (rtp%linear_scaling)
THEN
123 CALL dbcsr_multiply(
"N",
"N",
one, sinvh(ispin)%matrix, rho_new(re)%matrix,
zero, tmp, &
124 filter_eps=rtp%filter_eps)
125 IF (rtp%propagate_complex_ks) &
126 CALL dbcsr_multiply(
"N",
"N",
one, sinvh_imag(ispin)%matrix, rho_new(im)%matrix,
one, tmp, &
127 filter_eps=rtp%filter_eps)
128 CALL dbcsr_multiply(
"N",
"N",
one, sinvb(ispin)%matrix, rho_new(im)%matrix,
one, tmp, &
129 filter_eps=rtp%filter_eps)
130 CALL compute_forces(force, tmp, s_der, rho_new(im)%matrix, c_mat, kind_of, atom_of_kind)
132 CALL dbcsr_multiply(
"N",
"N",
one, sinvh(ispin)%matrix, rho_ao(ispin)%matrix,
zero, tmp)
133 IF (rtp%propagate_complex_ks) &
134 CALL dbcsr_multiply(
"N",
"N",
one, sinvh_imag(ispin)%matrix, rho_ao_im(ispin)%matrix,
one, tmp)
135 CALL dbcsr_multiply(
"N",
"N",
one, sinvb(ispin)%matrix, rho_ao_im(ispin)%matrix,
one, tmp)
136 CALL compute_forces(force, tmp, s_der, rho_ao_im(ispin)%matrix, c_mat, kind_of, atom_of_kind)
141 DO i = 1,
SIZE(force)
142 force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :)
145 CALL dbcsr_deallocate_matrix(tmp)
147 CALL timestop(handle)
161 SUBROUTINE compute_forces(force, tmp, S_der, rho_im, C_mat, kind_of, atom_of_kind)
162 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
163 TYPE(dbcsr_type),
POINTER :: tmp
164 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: s_der
165 TYPE(dbcsr_type),
POINTER :: rho_im
166 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: c_mat
167 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: kind_of, atom_of_kind
169 INTEGER :: col_atom, i, ikind, kind_atom, row_atom
171 REAL(
dp),
DIMENSION(:),
POINTER :: block_values, block_values2
172 TYPE(dbcsr_iterator_type) :: iter
178 CALL dbcsr_iterator_start(iter, tmp)
179 DO WHILE (dbcsr_iterator_blocks_left(iter))
180 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
181 CALL dbcsr_get_block_p(s_der(i)%matrix, row_atom, col_atom, block_values2, found=found)
183 ikind = kind_of(col_atom)
184 kind_atom = atom_of_kind(col_atom)
187 force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
188 2.0_dp*dot_product(block_values, block_values2)
191 CALL dbcsr_iterator_stop(iter)
195 CALL dbcsr_iterator_start(iter, rho_im)
196 DO WHILE (dbcsr_iterator_blocks_left(iter))
197 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
198 CALL dbcsr_get_block_p(c_mat(i)%matrix, row_atom, col_atom, block_values2, found=found)
200 ikind = kind_of(col_atom)
201 kind_atom = atom_of_kind(col_atom)
204 force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
205 2.0_dp*dot_product(block_values, block_values2)
208 CALL dbcsr_iterator_stop(iter)
211 END SUBROUTINE compute_forces
218 TYPE(qs_environment_type),
POINTER :: qs_env
220 TYPE(admm_type),
POINTER :: admm_env
221 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos, mos_admm
222 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_aux_im, ks_aux_re, matrix_s_aux_fit, &
223 matrix_s_aux_fit_vs_orb
224 TYPE(rt_prop_type),
POINTER :: rtp
229 CALL get_admm_env(admm_env, matrix_ks_aux_fit=ks_aux_re, &
230 matrix_ks_aux_fit_im=ks_aux_im, &
231 matrix_s_aux_fit=matrix_s_aux_fit, &
232 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
234 CALL get_rtp(rtp=rtp, mos_new=mos, admm_mos=mos_admm)
237 CALL rt_admm_forces_none(qs_env, admm_env, ks_aux_re, ks_aux_im, &
238 matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos)
253 SUBROUTINE rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos)
254 TYPE(qs_environment_type),
POINTER :: qs_env
255 TYPE(admm_type),
POINTER :: admm_env
256 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_aux_re, ks_aux_im, matrix_s_aux_fit, &
257 matrix_s_aux_fit_vs_orb
258 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos_admm, mos
260 INTEGER :: im, ispin, jspin, nao, natom, naux, nmo, &
262 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: admm_force
263 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
264 TYPE(cp_fm_struct_type),
POINTER :: mstruct
265 TYPE(cp_fm_type),
DIMENSION(2) :: tmp_aux_aux, tmp_aux_mo, tmp_aux_mo1, &
267 TYPE(dbcsr_type),
POINTER :: matrix_w_q, matrix_w_s
268 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
269 POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
270 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
271 TYPE(qs_ks_env_type),
POINTER :: ks_env
273 NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env)
276 CALL get_admm_env(admm_env, sab_aux_fit_asymm=sab_aux_fit_asymm, &
277 sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
279 ALLOCATE (matrix_w_s)
280 CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
281 name=
'W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry)
284 ALLOCATE (matrix_w_q)
285 CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
289 CALL cp_fm_create(tmp_aux_aux(jspin), admm_env%work_aux_aux%matrix_struct, name=
"taa")
290 CALL cp_fm_create(tmp_aux_nao(jspin), admm_env%work_aux_orb%matrix_struct, name=
"tao")
293 DO ispin = 1,
SIZE(ks_aux_re)
294 re = 2*ispin - 1; im = 2*ispin
295 naux = admm_env%nao_aux_fit; nmo = admm_env%nmo(ispin); nao = admm_env%nao_orb
297 mstruct => admm_env%work_aux_nmo(ispin)%matrix_struct
299 CALL cp_fm_create(tmp_aux_mo(jspin), mstruct, name=
"tam")
300 CALL cp_fm_create(tmp_aux_mo1(jspin), mstruct, name=
"tam")
310 CALL parallel_gemm(
'N',
'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(re), 0.0_dp, tmp_aux_mo1(re))
311 CALL parallel_gemm(
'N',
'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(im), 0.0_dp, tmp_aux_mo1(im))
315 CALL parallel_gemm(
"N",
"T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(re), mos(re), 0.0_dp, &
317 CALL parallel_gemm(
"N",
"T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(im), 1.0_dp, &
319 CALL parallel_gemm(
"N",
"T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(re), mos(im), 0.0_dp, &
321 CALL parallel_gemm(
"N",
"T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(re), 1.0_dp, &
325 CALL parallel_gemm(
'N',
'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(re), admm_env%A, 0.0_dp, tmp_aux_aux(re))
326 CALL parallel_gemm(
'N',
'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(im), admm_env%A, 0.0_dp, tmp_aux_aux(im))
335 CALL cp_fm_release(tmp_aux_mo(jspin))
336 CALL cp_fm_release(tmp_aux_mo1(jspin))
342 ALLOCATE (admm_force(3, natom))
345 basis_type_a=
"AUX_FIT", basis_type_b=
"AUX_FIT", &
346 sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
348 basis_type_a=
"AUX_FIT", basis_type_b=
"ORB", &
349 sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
351 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
353 CALL add_qs_force(admm_force, force,
"overlap_admm", atomic_kind_set)
354 DEALLOCATE (admm_force)
357 CALL dbcsr_deallocate_matrix(matrix_w_s)
358 CALL dbcsr_deallocate_matrix(matrix_w_q)
362 CALL cp_fm_release(tmp_aux_aux(jspin))
363 CALL cp_fm_release(tmp_aux_nao(jspin))
366 END SUBROUTINE rt_admm_forces_none
Types and set/get functions for auxiliary density matrix methods.
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
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.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
represent the structure of a full matrix
represent a full matrix distributed on many processors
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public one
real(kind=dp), parameter, public zero
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
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.
subroutine, public add_qs_force(force, qs_force, forcetype, atomic_kind_set)
Add force to a force_type variable.
Define the neighbor list data types and the corresponding functionality.
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_force(ks_env, force, basis_type_a, basis_type_b, sab_nl, matrix_p, matrixkp_p)
Calculation of the force contribution from an overlap matrix over Cartesian Gaussian functions.
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...
subroutine, public rt_admm_force(qs_env)
...
subroutine, public calc_c_mat_force(qs_env)
calculates the three additional force contributions needed in EMD P_imag*C , P_imag*B*S^-1*S_der ,...
Types and set_get for real time propagation depending on runtype and diagonalization method different...
subroutine, public get_rtp(rtp, exp_H_old, exp_H_new, H_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, S_inv, S_half, S_minus_half, B_mat, C_mat, propagator_matrix, mixing, mixing_factor, S_der, dt, nsteps, SinvH, SinvH_imag, SinvB, admm_mos)
...