Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief tblite matrix build
10 : !> \author JVP
11 : !> \history creation 09.2024
12 : ! **************************************************************************************************
13 :
14 : MODULE tblite_ks_matrix
15 :
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
18 : dbcsr_copy,&
19 : dbcsr_multiply,&
20 : dbcsr_p_type,&
21 : dbcsr_type
22 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
23 : USE kinds, ONLY: dp
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE qs_energy_types, ONLY: qs_energy_type
26 : USE qs_environment_types, ONLY: get_qs_env,&
27 : qs_environment_type
28 : USE qs_ks_types, ONLY: qs_ks_env_type
29 : USE qs_mo_types, ONLY: get_mo_set,&
30 : mo_set_type
31 : USE qs_rho_types, ONLY: qs_rho_get,&
32 : qs_rho_type
33 : USE tblite_interface, ONLY: tb_derive_dH_diag,&
34 : tb_derive_dH_off,&
35 : tb_ham_add_coulomb,&
36 : tb_update_charges
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_ks_matrix'
44 :
45 : PUBLIC :: build_tblite_ks_matrix
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief ...
51 : !> \param qs_env ...
52 : !> \param calculate_forces ...
53 : !> \param just_energy ...
54 : !> \param ext_ks_matrix ...
55 : ! **************************************************************************************************
56 25710 : SUBROUTINE build_tblite_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
57 : TYPE(qs_environment_type), POINTER :: qs_env
58 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
59 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
60 : POINTER :: ext_ks_matrix
61 :
62 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tblite_ks_matrix'
63 :
64 : INTEGER :: handle, img, ispin, nimg, ns, nspins
65 : LOGICAL :: do_efield, force_use_rho
66 : REAL(KIND=dp) :: pc_ener, qmmm_el
67 25710 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
68 25710 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h
69 : TYPE(dbcsr_type), POINTER :: mo_coeff
70 : TYPE(dft_control_type), POINTER :: dft_control
71 : TYPE(mp_para_env_type), POINTER :: para_env
72 : TYPE(qs_energy_type), POINTER :: energy
73 : TYPE(qs_ks_env_type), POINTER :: ks_env
74 : TYPE(qs_rho_type), POINTER :: rho
75 :
76 25710 : CALL timeset(routineN, handle)
77 :
78 25710 : NULLIFY (dft_control, ks_env, ks_matrix, rho, energy)
79 25710 : CPASSERT(ASSOCIATED(qs_env))
80 :
81 : CALL get_qs_env(qs_env, &
82 : dft_control=dft_control, &
83 : matrix_h_kp=matrix_h, &
84 : para_env=para_env, &
85 : ks_env=ks_env, &
86 : matrix_ks_kp=ks_matrix, &
87 : rho=rho, &
88 25710 : energy=energy)
89 :
90 25710 : IF (PRESENT(ext_ks_matrix)) THEN
91 : ! remap pointer to allow for non-kpoint external ks matrix
92 : ! ext_ks_matrix is used in linear response code
93 2 : ns = SIZE(ext_ks_matrix)
94 2 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
95 : END IF
96 :
97 25710 : energy%qmmm_el = 0.0_dp
98 :
99 25710 : nspins = dft_control%nspins
100 25710 : nimg = dft_control%nimages
101 25710 : CPASSERT(ASSOCIATED(matrix_h))
102 25710 : CPASSERT(ASSOCIATED(rho))
103 77130 : CPASSERT(SIZE(ks_matrix) > 0)
104 :
105 54028 : DO ispin = 1, nspins
106 450302 : DO img = 1, nimg
107 : ! copy the core matrix into the fock matrix
108 424592 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
109 : END DO
110 : END DO
111 :
112 25710 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
113 : dft_control%apply_efield_field) THEN
114 0 : do_efield = .TRUE.
115 0 : CPABORT("Not implemented yet. Use CP2K routines for GFN1")
116 : ELSE
117 25710 : do_efield = .FALSE.
118 : END IF
119 :
120 25710 : force_use_rho = .NOT. (calculate_forces .AND. nspins > 1)
121 25710 : CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, calculate_forces, force_use_rho)
122 :
123 25710 : CALL tb_ham_add_coulomb(qs_env, qs_env%tb_tblite, dft_control)
124 :
125 25710 : IF (qs_env%qmmm) THEN
126 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
127 0 : DO ispin = 1, nspins
128 : ! If QM/MM sumup the 1el Hamiltonian
129 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
130 0 : 1.0_dp, 1.0_dp)
131 0 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
132 : ! Compute QM/MM Energy
133 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
134 0 : matrix_p1(ispin)%matrix, qmmm_el)
135 0 : energy%qmmm_el = energy%qmmm_el + qmmm_el
136 : END DO
137 0 : pc_ener = qs_env%ks_qmmm_env%pc_ener
138 0 : energy%qmmm_el = energy%qmmm_el + pc_ener
139 : END IF
140 :
141 25710 : IF (calculate_forces) THEN
142 114 : CALL tb_derive_dH_diag(qs_env, force_use_rho, nimg)
143 114 : CALL tb_derive_dH_off(qs_env, force_use_rho, nimg)
144 : END IF
145 :
146 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
147 25710 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
148 144 : CPASSERT(SIZE(ks_matrix, 2) == 1)
149 : BLOCK
150 144 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
151 144 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
152 374 : DO ispin = 1, SIZE(mo_derivs)
153 230 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
154 230 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
155 0 : CPABORT("")
156 : END IF
157 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
158 374 : 0.0_dp, mo_derivs(ispin)%matrix)
159 : END DO
160 : END BLOCK
161 : END IF
162 :
163 25710 : CALL timestop(handle)
164 :
165 25710 : END SUBROUTINE build_tblite_ks_matrix
166 :
167 : END MODULE tblite_ks_matrix
|
