Line data Source code
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 Calculates the energy contribution and the mo_derivative of
10 : !> a static periodic electric field
11 : !> \par History
12 : !> none
13 : !> \author fschiff (06.2010)
14 : ! **************************************************************************************************
15 : MODULE qs_efield_berry
16 : USE ai_moments, ONLY: cossin
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind,&
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE block_p_types, ONLY: block_p_type
23 : USE cell_types, ONLY: cell_type,&
24 : pbc
25 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm,&
26 : cp_cfm_solve
27 : USE cp_cfm_types, ONLY: cp_cfm_create,&
28 : cp_cfm_release,&
29 : cp_cfm_set_all,&
30 : cp_cfm_type
31 : USE cp_control_types, ONLY: dft_control_type
32 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
33 : dbcsr_get_block_p,&
34 : dbcsr_p_type,&
35 : dbcsr_set,&
36 : dbcsr_type
37 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
38 : copy_fm_to_dbcsr,&
39 : cp_dbcsr_plus_fm_fm_t,&
40 : cp_dbcsr_sm_fm_multiply,&
41 : dbcsr_deallocate_matrix_set
42 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
43 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
44 : cp_fm_struct_release,&
45 : cp_fm_struct_type
46 : USE cp_fm_types, ONLY: cp_fm_create,&
47 : cp_fm_release,&
48 : cp_fm_set_all,&
49 : cp_fm_type
50 : USE kinds, ONLY: dp
51 : USE mathconstants, ONLY: gaussi,&
52 : pi,&
53 : twopi,&
54 : z_one,&
55 : z_zero
56 : USE message_passing, ONLY: mp_para_env_type
57 : USE orbital_pointers, ONLY: ncoset
58 : USE parallel_gemm_api, ONLY: parallel_gemm
59 : USE particle_types, ONLY: particle_type
60 : USE qs_energy_types, ONLY: qs_energy_type
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type,&
63 : set_qs_env
64 : USE qs_force_types, ONLY: qs_force_type
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : get_qs_kind_set,&
67 : qs_kind_type
68 : USE qs_mo_types, ONLY: get_mo_set,&
69 : mo_set_type
70 : USE qs_moments, ONLY: build_berry_moment_matrix
71 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
72 : neighbor_list_iterate,&
73 : neighbor_list_iterator_create,&
74 : neighbor_list_iterator_p_type,&
75 : neighbor_list_iterator_release,&
76 : neighbor_list_set_p_type
77 : USE qs_period_efield_types, ONLY: efield_berry_type,&
78 : init_efield_matrices,&
79 : set_efield_matrices
80 : USE virial_methods, ONLY: virial_pair_force
81 : USE virial_types, ONLY: virial_type
82 : #include "./base/base_uses.f90"
83 :
84 : IMPLICIT NONE
85 :
86 : PRIVATE
87 :
88 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_efield_berry'
89 :
90 : ! *** Public subroutines ***
91 :
92 : PUBLIC :: qs_efield_berry_phase
93 :
94 : ! **************************************************************************************************
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 :
100 : ! **************************************************************************************************
101 : !> \brief ...
102 : !> \param qs_env ...
103 : !> \param just_energy ...
104 : !> \param calculate_forces ...
105 : ! **************************************************************************************************
106 98089 : SUBROUTINE qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
107 :
108 : TYPE(qs_environment_type), POINTER :: qs_env
109 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
110 :
111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_berry_phase'
112 :
113 : INTEGER :: handle
114 : LOGICAL :: s_mstruct_changed
115 : TYPE(dft_control_type), POINTER :: dft_control
116 :
117 98089 : CALL timeset(routineN, handle)
118 :
119 98089 : NULLIFY (dft_control)
120 : CALL get_qs_env(qs_env, s_mstruct_changed=s_mstruct_changed, &
121 98089 : dft_control=dft_control)
122 :
123 98089 : IF (dft_control%apply_period_efield) THEN
124 2058 : IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env)
125 2058 : IF (dft_control%period_efield%displacement_field) THEN
126 898 : CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
127 : ELSE
128 1160 : CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces)
129 : END IF
130 : END IF
131 :
132 98089 : CALL timestop(handle)
133 :
134 98089 : END SUBROUTINE qs_efield_berry_phase
135 :
136 : ! **************************************************************************************************
137 : !> \brief ...
138 : !> \param qs_env ...
139 : ! **************************************************************************************************
140 144 : SUBROUTINE qs_efield_integrals(qs_env)
141 :
142 : TYPE(qs_environment_type), POINTER :: qs_env
143 :
144 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_integrals'
145 :
146 : INTEGER :: handle, i
147 : REAL(dp), DIMENSION(3) :: kvec
148 : TYPE(cell_type), POINTER :: cell
149 144 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: cosmat, matrix_s, sinmat
150 : TYPE(dft_control_type), POINTER :: dft_control
151 : TYPE(efield_berry_type), POINTER :: efield
152 :
153 144 : CALL timeset(routineN, handle)
154 144 : CPASSERT(ASSOCIATED(qs_env))
155 :
156 144 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
157 144 : NULLIFY (matrix_s)
158 144 : CALL get_qs_env(qs_env=qs_env, efield=efield, cell=cell, matrix_s=matrix_s)
159 144 : CALL init_efield_matrices(efield)
160 1008 : ALLOCATE (cosmat(3), sinmat(3))
161 576 : DO i = 1, 3
162 432 : ALLOCATE (cosmat(i)%matrix, sinmat(i)%matrix)
163 :
164 432 : CALL dbcsr_copy(cosmat(i)%matrix, matrix_s(1)%matrix, 'COS MAT')
165 432 : CALL dbcsr_copy(sinmat(i)%matrix, matrix_s(1)%matrix, 'SIN MAT')
166 432 : CALL dbcsr_set(cosmat(i)%matrix, 0.0_dp)
167 432 : CALL dbcsr_set(sinmat(i)%matrix, 0.0_dp)
168 :
169 1728 : kvec(:) = twopi*cell%h_inv(i, :)
170 576 : CALL build_berry_moment_matrix(qs_env, cosmat(i)%matrix, sinmat(i)%matrix, kvec)
171 : END DO
172 144 : CALL set_efield_matrices(efield=efield, cosmat=cosmat, sinmat=sinmat)
173 144 : CALL set_qs_env(qs_env=qs_env, efield=efield)
174 144 : CALL timestop(handle)
175 :
176 144 : END SUBROUTINE qs_efield_integrals
177 :
178 : ! **************************************************************************************************
179 : !> \brief ...
180 : !> \param qs_env ...
181 : !> \param just_energy ...
182 : !> \param calculate_forces ...
183 : ! **************************************************************************************************
184 1160 : SUBROUTINE qs_efield_derivatives(qs_env, just_energy, calculate_forces)
185 : TYPE(qs_environment_type), POINTER :: qs_env
186 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
187 :
188 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_derivatives'
189 :
190 : COMPLEX(dp) :: zdet, zdeta, zi(3)
191 : INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, j, &
192 : jatom, jkind, jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, &
193 : nseta, nsetb, sgfa, sgfb
194 1160 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
195 1160 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
196 1160 : npgfb, nsgfa, nsgfb
197 1160 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
198 : LOGICAL :: found, uniform, use_virial
199 : REAL(dp) :: charge, ci(3), cqi(3), dab, dd, &
200 : ener_field, f0, fab, fieldpol(3), &
201 : focc, fpolvec(3), hmat(3, 3), occ, &
202 : qi(3), ti(3)
203 : REAL(dp), DIMENSION(3) :: forcea, forceb, kvec, ra, rab, rb, ria
204 2320 : REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
205 2320 : REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab
206 1160 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
207 1160 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
208 1160 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
209 20880 : TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
210 : TYPE(cell_type), POINTER :: cell
211 1160 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
212 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
213 1160 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp, mo_derivs_tmp
214 1160 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, op_fm_set, opvec
215 : TYPE(cp_fm_type), POINTER :: mo_coeff
216 1160 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
217 1160 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
218 : TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
219 : TYPE(dft_control_type), POINTER :: dft_control
220 : TYPE(efield_berry_type), POINTER :: efield
221 1160 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
222 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
223 1160 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
224 : TYPE(mp_para_env_type), POINTER :: para_env
225 : TYPE(neighbor_list_iterator_p_type), &
226 1160 : DIMENSION(:), POINTER :: nl_iterator
227 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
228 1160 : POINTER :: sab_orb
229 1160 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
230 : TYPE(qs_energy_type), POINTER :: energy
231 1160 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
232 1160 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
233 : TYPE(qs_kind_type), POINTER :: qs_kind
234 : TYPE(virial_type), POINTER :: virial
235 :
236 1160 : CALL timeset(routineN, handle)
237 :
238 1160 : NULLIFY (dft_control, cell, particle_set)
239 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
240 1160 : particle_set=particle_set, virial=virial)
241 1160 : NULLIFY (qs_kind_set, efield, para_env, sab_orb)
242 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
243 1160 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
244 :
245 : ! calculate stress only if forces requested also
246 1160 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
247 0 : use_virial = use_virial .AND. calculate_forces
248 : ! disable stress calculation
249 : IF (use_virial) THEN
250 0 : CPABORT("Stress tensor for periodic E-field not implemented")
251 : END IF
252 :
253 4640 : fieldpol = dft_control%period_efield%polarisation
254 8120 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
255 4640 : fieldpol = -fieldpol*dft_control%period_efield%strength
256 15080 : hmat = cell%hmat(:, :)/twopi
257 4640 : DO idir = 1, 3
258 4640 : fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
259 : END DO
260 :
261 : ! nuclear contribution
262 1160 : natom = SIZE(particle_set)
263 1160 : IF (calculate_forces) THEN
264 28 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
265 28 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
266 : END IF
267 4640 : zi(:) = CMPLX(1._dp, 0._dp, dp)
268 3960 : DO ia = 1, natom
269 2800 : CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
270 2800 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
271 11200 : ria = particle_set(ia)%r
272 11200 : ria = pbc(ria, cell)
273 11200 : DO idir = 1, 3
274 33600 : kvec(:) = twopi*cell%h_inv(idir, :)
275 33600 : dd = SUM(kvec(:)*ria(:))
276 8400 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
277 11200 : zi(idir) = zi(idir)*zdeta
278 : END DO
279 2800 : IF (calculate_forces) THEN
280 74 : IF (para_env%mepos == 0) THEN
281 37 : iatom = atom_of_kind(ia)
282 148 : forcea(:) = fieldpol(:)*charge
283 148 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + forcea(:)
284 : END IF
285 : END IF
286 6760 : IF (use_virial) THEN
287 0 : IF (para_env%mepos == 0) &
288 0 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria)
289 : END IF
290 : END DO
291 4640 : qi = AIMAG(LOG(zi))
292 :
293 : ! check uniform occupation
294 1160 : NULLIFY (mos)
295 1160 : CALL get_qs_env(qs_env=qs_env, mos=mos)
296 2380 : DO ispin = 1, dft_control%nspins
297 1220 : CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
298 2380 : IF (.NOT. uniform) THEN
299 0 : CPABORT("Berry phase moments for non uniform MOs' occupation numbers not implemented")
300 : END IF
301 : END DO
302 :
303 1160 : NULLIFY (mo_derivs)
304 1160 : CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
305 : ! initialize all work matrices needed
306 7140 : ALLOCATE (op_fm_set(2, dft_control%nspins))
307 7140 : ALLOCATE (opvec(2, dft_control%nspins))
308 4700 : ALLOCATE (eigrmat(dft_control%nspins))
309 4700 : ALLOCATE (inv_mat(dft_control%nspins))
310 7140 : ALLOCATE (inv_work(2, dft_control%nspins))
311 4700 : ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs)))
312 3540 : ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
313 :
314 : ! Allocate temp matrices for the wavefunction derivatives
315 2380 : DO ispin = 1, dft_control%nspins
316 1220 : NULLIFY (tmp_fm_struct, mo_coeff)
317 1220 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
318 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
319 1220 : ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
320 1220 : CALL cp_fm_create(mo_derivs_tmp(ispin), mo_coeff%matrix_struct)
321 1220 : CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
322 1220 : CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_tmp(ispin))
323 3660 : DO i = 1, SIZE(op_fm_set, 1)
324 2440 : CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
325 2440 : CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
326 3660 : CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
327 : END DO
328 1220 : CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
329 1220 : CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
330 3600 : CALL cp_fm_struct_release(tmp_fm_struct)
331 : END DO
332 : ! temp matrices for force calculation
333 1160 : IF (calculate_forces) THEN
334 28 : NULLIFY (matrix_s)
335 28 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
336 192 : ALLOCATE (tempmat(2, dft_control%nspins))
337 64 : DO ispin = 1, dft_control%nspins
338 36 : ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
339 36 : CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
340 36 : CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
341 36 : CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
342 64 : CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
343 : END DO
344 : ! integration
345 28 : CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
346 224 : ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
347 196 : ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
348 28 : lsab = MAX(ldab, lsab)
349 140 : DO i = 1, 3
350 504 : ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
351 448 : ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
352 : END DO
353 : END IF
354 :
355 : !Start the MO derivative calculation
356 : !loop over all cell vectors
357 4640 : DO idir = 1, 3
358 3480 : ci(idir) = 0.0_dp
359 : zi(idir) = z_zero
360 4640 : IF (ABS(fpolvec(idir)) .GT. 1.0E-12_dp) THEN
361 2442 : cosmat => efield%cosmat(idir)%matrix
362 2442 : sinmat => efield%sinmat(idir)%matrix
363 : !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
364 : !first step S_berry * C and C^T S_berry C
365 4944 : DO ispin = 1, dft_control%nspins ! spin
366 2502 : IF (mos(ispin)%use_mo_coeff_b) THEN
367 2502 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
368 2502 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
369 : ELSE
370 0 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
371 0 : mo_coeff_tmp(ispin) = mo_coeff
372 : END IF
373 2502 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
374 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
375 2502 : op_fm_set(1, ispin))
376 2502 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
377 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
378 4944 : op_fm_set(2, ispin))
379 : END DO
380 : !second step invert C^T S_berry C
381 2442 : zdet = z_one
382 4944 : DO ispin = 1, dft_control%nspins
383 2502 : CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
384 2502 : CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
385 2502 : CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
386 2502 : CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
387 4944 : zdet = zdet*zdeta
388 : END DO
389 2442 : zi(idir) = zdet**occ
390 2442 : ci(idir) = AIMAG(LOG(zdet**occ))
391 :
392 2442 : IF (.NOT. just_energy) THEN
393 : !compute the orbital derivative
394 2390 : focc = fpolvec(idir)
395 4818 : DO ispin = 1, dft_control%nspins
396 31564 : inv_work(1, ispin)%local_data(:, :) = REAL(inv_mat(ispin)%local_data(:, :), dp)
397 31564 : inv_work(2, ispin)%local_data(:, :) = AIMAG(inv_mat(ispin)%local_data(:, :))
398 2428 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
399 : CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
400 2428 : 1.0_dp, mo_derivs_tmp(ispin))
401 : CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
402 7246 : 1.0_dp, mo_derivs_tmp(ispin))
403 : END DO
404 : END IF
405 :
406 : !compute nuclear forces
407 2442 : IF (calculate_forces) THEN
408 46 : nkind = SIZE(qs_kind_set)
409 46 : natom = SIZE(particle_set)
410 184 : kvec(:) = twopi*cell%h_inv(idir, :)
411 :
412 : ! calculate: C [C^T S_berry C]^(-1) C^T
413 : ! Store this matrix in DBCSR form (only S overlap blocks)
414 100 : DO ispin = 1, dft_control%nspins
415 54 : CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
416 54 : CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
417 54 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
418 : CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
419 54 : opvec(1, ispin))
420 : CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
421 54 : opvec(2, ispin))
422 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
423 54 : matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
424 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
425 154 : matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
426 : END DO
427 :
428 : ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
429 230 : ALLOCATE (basis_set_list(nkind))
430 138 : DO ikind = 1, nkind
431 92 : qs_kind => qs_kind_set(ikind)
432 92 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
433 138 : IF (ASSOCIATED(basis_set_a)) THEN
434 92 : basis_set_list(ikind)%gto_basis_set => basis_set_a
435 : ELSE
436 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
437 : END IF
438 : END DO
439 : !
440 46 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
441 1472 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
442 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
443 1426 : iatom=iatom, jatom=jatom, r=rab)
444 1426 : basis_set_a => basis_set_list(ikind)%gto_basis_set
445 1426 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
446 1426 : basis_set_b => basis_set_list(jkind)%gto_basis_set
447 1426 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
448 : ! basis ikind
449 1426 : first_sgfa => basis_set_a%first_sgf
450 1426 : la_max => basis_set_a%lmax
451 1426 : la_min => basis_set_a%lmin
452 1426 : npgfa => basis_set_a%npgf
453 1426 : nseta = basis_set_a%nset
454 1426 : nsgfa => basis_set_a%nsgf_set
455 1426 : rpgfa => basis_set_a%pgf_radius
456 1426 : set_radius_a => basis_set_a%set_radius
457 1426 : sphi_a => basis_set_a%sphi
458 1426 : zeta => basis_set_a%zet
459 : ! basis jkind
460 1426 : first_sgfb => basis_set_b%first_sgf
461 1426 : lb_max => basis_set_b%lmax
462 1426 : lb_min => basis_set_b%lmin
463 1426 : npgfb => basis_set_b%npgf
464 1426 : nsetb = basis_set_b%nset
465 1426 : nsgfb => basis_set_b%nsgf_set
466 1426 : rpgfb => basis_set_b%pgf_radius
467 1426 : set_radius_b => basis_set_b%set_radius
468 1426 : sphi_b => basis_set_b%sphi
469 1426 : zetb => basis_set_b%zet
470 :
471 1426 : atom_a = atom_of_kind(iatom)
472 1426 : atom_b = atom_of_kind(jatom)
473 :
474 1426 : ldsa = SIZE(sphi_a, 1)
475 1426 : ldsb = SIZE(sphi_b, 1)
476 1426 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
477 5704 : rb(:) = ra + rab
478 1426 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
479 :
480 1426 : IF (iatom <= jatom) THEN
481 960 : irow = iatom
482 960 : icol = jatom
483 : ELSE
484 466 : irow = jatom
485 466 : icol = iatom
486 : END IF
487 :
488 1426 : IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
489 : fab = 1.0_dp*occ
490 : ELSE
491 1372 : fab = 2.0_dp*occ
492 : END IF
493 :
494 5704 : DO i = 1, 3
495 1167894 : dcost(i, 1)%block = 0.0_dp
496 1167894 : dsint(i, 1)%block = 0.0_dp
497 1167894 : dcost(i, 2)%block = 0.0_dp
498 1169320 : dsint(i, 2)%block = 0.0_dp
499 : END DO
500 :
501 4278 : DO iset = 1, nseta
502 2852 : ncoa = npgfa(iset)*ncoset(la_max(iset))
503 2852 : sgfa = first_sgfa(1, iset)
504 9982 : DO jset = 1, nsetb
505 5704 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
506 2744 : ncob = npgfb(jset)*ncoset(lb_max(jset))
507 2744 : sgfb = first_sgfb(1, jset)
508 : ! Calculate the primitive integrals (da|b)
509 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
510 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
511 2744 : ra, rb, kvec, cosab, sinab, dcosab, dsinab)
512 10976 : DO i = 1, 3
513 : CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
514 : ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
515 : ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
516 10976 : dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
517 : END DO
518 : ! Calculate the primitive integrals (a|db)
519 : CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
520 : la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
521 2744 : rb, ra, kvec, cosab, sinab, dcosab, dsinab)
522 13828 : DO i = 1, 3
523 1264188 : dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i))
524 1264188 : dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i))
525 : CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
526 : ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
527 : ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
528 13936 : dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
529 : END DO
530 : END DO
531 : END DO
532 1426 : forcea = 0.0_dp
533 1426 : forceb = 0.0_dp
534 3252 : DO ispin = 1, dft_control%nspins
535 1826 : NULLIFY (rblock, iblock)
536 : CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
537 1826 : row=irow, col=icol, BLOCK=rblock, found=found)
538 1826 : CPASSERT(found)
539 : CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
540 1826 : row=irow, col=icol, BLOCK=iblock, found=found)
541 1826 : CPASSERT(found)
542 1826 : n1 = SIZE(rblock, 1)
543 1826 : n2 = SIZE(rblock, 2)
544 1826 : CPASSERT(SIZE(iblock, 1) == n1)
545 1826 : CPASSERT(SIZE(iblock, 2) == n2)
546 1826 : CPASSERT(lsab >= n1)
547 1826 : CPASSERT(lsab >= n2)
548 6904 : IF (iatom <= jatom) THEN
549 4976 : DO i = 1, 3
550 : forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
551 572148 : - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
552 : forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
553 573392 : - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
554 : END DO
555 : ELSE
556 2328 : DO i = 1, 3
557 : forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
558 211626 : - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
559 : forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
560 212208 : - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
561 : END DO
562 : END IF
563 : END DO
564 5704 : force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fab*fpolvec(idir)*forcea(1:3)
565 5704 : force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fab*fpolvec(idir)*forceb(1:3)
566 1472 : IF (use_virial) THEN
567 0 : f0 = -fab*fpolvec(idir)
568 0 : CALL virial_pair_force(virial%pv_virial, f0, forcea, ra)
569 0 : CALL virial_pair_force(virial%pv_virial, f0, forceb, rb)
570 : END IF
571 :
572 : END DO
573 46 : CALL neighbor_list_iterator_release(nl_iterator)
574 46 : DEALLOCATE (basis_set_list)
575 :
576 : END IF
577 : END IF
578 : END DO
579 :
580 : ! Energy
581 4640 : ener_field = 0.0_dp
582 : ti = 0.0_dp
583 4640 : DO idir = 1, 3
584 : ! make sure the total normalized polarization is within [-1:1]
585 3480 : cqi(idir) = qi(idir) + ci(idir)
586 3480 : IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
587 3480 : IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
588 : ! now check for log branch
589 3480 : IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
590 0 : ti(idir) = (efield%polarisation(idir) - cqi(idir))/pi
591 0 : DO i = 1, 10
592 0 : cqi(idir) = cqi(idir) + SIGN(1.0_dp, ti(idir))*twopi
593 0 : IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
594 : END DO
595 : END IF
596 4640 : ener_field = ener_field + fpolvec(idir)*cqi(idir)
597 : END DO
598 :
599 : ! update the references
600 1160 : IF (calculate_forces) THEN
601 : ! check for smoothness of energy surface
602 112 : IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(fpolvec))) THEN
603 0 : CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
604 : END IF
605 28 : efield%field_energy = ener_field
606 112 : efield%polarisation(:) = cqi(:)
607 : END IF
608 1160 : energy%efield = ener_field
609 :
610 1160 : IF (.NOT. just_energy) THEN
611 : ! Add the result to mo_derivativs
612 2254 : DO ispin = 1, dft_control%nspins
613 2254 : CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin), mo_derivs(ispin)%matrix)
614 : END DO
615 1108 : IF (use_virial) THEN
616 0 : ti = 0.0_dp
617 0 : DO i = 1, 3
618 0 : DO j = 1, 3
619 0 : ti(j) = ti(j) + hmat(j, i)*cqi(i)
620 : END DO
621 : END DO
622 0 : DO i = 1, 3
623 0 : DO j = 1, 3
624 0 : virial%pv_virial(i, j) = virial%pv_virial(i, j) - fieldpol(i)*ti(j)
625 : END DO
626 : END DO
627 : END IF
628 : END IF
629 :
630 2380 : DO ispin = 1, dft_control%nspins
631 1220 : CALL cp_cfm_release(eigrmat(ispin))
632 1220 : CALL cp_cfm_release(inv_mat(ispin))
633 1220 : CALL cp_fm_release(mo_derivs_tmp(ispin))
634 1220 : IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
635 4820 : DO i = 1, SIZE(op_fm_set, 1)
636 2440 : CALL cp_fm_release(opvec(i, ispin))
637 2440 : CALL cp_fm_release(op_fm_set(i, ispin))
638 3660 : CALL cp_fm_release(inv_work(i, ispin))
639 : END DO
640 : END DO
641 1160 : DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
642 1160 : DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
643 :
644 1160 : IF (calculate_forces) THEN
645 84 : DO ikind = 1, SIZE(atomic_kind_set)
646 676 : CALL para_env%sum(force(ikind)%efield)
647 : END DO
648 28 : DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
649 112 : DO i = 1, 3
650 84 : DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
651 112 : DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
652 : END DO
653 28 : CALL dbcsr_deallocate_matrix_set(tempmat)
654 : END IF
655 1160 : CALL timestop(handle)
656 :
657 4640 : END SUBROUTINE qs_efield_derivatives
658 :
659 : ! **************************************************************************************************
660 : !> \brief ...
661 : !> \param qs_env ...
662 : !> \param just_energy ...
663 : !> \param calculate_forces ...
664 : ! **************************************************************************************************
665 898 : SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
666 : TYPE(qs_environment_type), POINTER :: qs_env
667 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
668 :
669 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_dispfield_derivatives'
670 :
671 : COMPLEX(dp) :: zdet, zdeta, zi(3)
672 : INTEGER :: handle, i, ia, iatom, icol, idir, ikind, iodeb, irow, iset, ispin, jatom, jkind, &
673 : jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, &
674 : sgfa, sgfb
675 898 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
676 898 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
677 898 : npgfb, nsgfa, nsgfb
678 898 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
679 : LOGICAL :: found, uniform, use_virial
680 : REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), &
681 : ener_field, fab, fieldpol(3), focc, &
682 : hmat(3, 3), occ, omega, qi(3), &
683 : rlog(3), zlog(3)
684 : REAL(dp), DIMENSION(3) :: dfilter, forcea, forceb, kvec, ra, rab, &
685 : rb, ria
686 1796 : REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
687 2694 : REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab, force_tmp
688 898 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
689 898 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
690 898 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
691 16164 : TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
692 : TYPE(cell_type), POINTER :: cell
693 898 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
694 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
695 898 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp
696 898 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, mo_derivs_tmp, op_fm_set, opvec
697 : TYPE(cp_fm_type), POINTER :: mo_coeff
698 898 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
699 898 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
700 : TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
701 : TYPE(dft_control_type), POINTER :: dft_control
702 : TYPE(efield_berry_type), POINTER :: efield
703 898 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
704 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
705 898 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
706 : TYPE(mp_para_env_type), POINTER :: para_env
707 : TYPE(neighbor_list_iterator_p_type), &
708 898 : DIMENSION(:), POINTER :: nl_iterator
709 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
710 898 : POINTER :: sab_orb
711 898 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
712 : TYPE(qs_energy_type), POINTER :: energy
713 898 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
714 898 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
715 : TYPE(qs_kind_type), POINTER :: qs_kind
716 : TYPE(virial_type), POINTER :: virial
717 :
718 898 : CALL timeset(routineN, handle)
719 :
720 898 : NULLIFY (dft_control, cell, particle_set)
721 : CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
722 898 : particle_set=particle_set, virial=virial)
723 898 : NULLIFY (qs_kind_set, efield, para_env, sab_orb)
724 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
725 898 : efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
726 :
727 : ! calculate stress only if forces requested also
728 898 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
729 0 : use_virial = use_virial .AND. calculate_forces
730 : ! disable stress calculation
731 : IF (use_virial) THEN
732 0 : CPABORT("Stress tensor for periodic D-field not implemented")
733 : END IF
734 :
735 3592 : dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
736 :
737 3592 : fieldpol = dft_control%period_efield%polarisation
738 6286 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
739 3592 : fieldpol = fieldpol*dft_control%period_efield%strength
740 :
741 898 : omega = cell%deth
742 11674 : hmat = cell%hmat(:, :)/(twopi*omega)
743 :
744 : ! nuclear contribution to polarization
745 898 : natom = SIZE(particle_set)
746 898 : IF (calculate_forces) THEN
747 10 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
748 10 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
749 40 : ALLOCATE (force_tmp(natom, 3, 3))
750 310 : force_tmp = 0.0_dp
751 : END IF
752 3592 : zi(:) = CMPLX(1._dp, 0._dp, dp)
753 2694 : DO ia = 1, natom
754 1796 : CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
755 1796 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
756 7184 : ria = particle_set(ia)%r
757 7184 : ria = pbc(ria, cell)
758 7184 : DO idir = 1, 3
759 21552 : kvec(:) = twopi*cell%h_inv(idir, :)
760 21552 : dd = SUM(kvec(:)*ria(:))
761 5388 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
762 7184 : zi(idir) = zi(idir)*zdeta
763 : END DO
764 4490 : IF (calculate_forces) THEN
765 20 : IF (para_env%mepos == 0) THEN
766 40 : DO i = 1, 3
767 40 : force_tmp(ia, i, i) = force_tmp(ia, i, i) + charge/omega
768 : END DO
769 : END IF
770 : END IF
771 : END DO
772 3592 : rlog = AIMAG(LOG(zi))
773 :
774 : ! check uniform occupation
775 898 : NULLIFY (mos)
776 898 : CALL get_qs_env(qs_env=qs_env, mos=mos)
777 1796 : DO ispin = 1, dft_control%nspins
778 898 : CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
779 1796 : IF (.NOT. uniform) THEN
780 0 : CPABORT("Berry phase moments for non uniform MO occupation numbers not implemented")
781 : END IF
782 : END DO
783 :
784 : ! initialize all work matrices needed
785 898 : NULLIFY (mo_derivs)
786 898 : CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
787 5388 : ALLOCATE (op_fm_set(2, dft_control%nspins))
788 5388 : ALLOCATE (opvec(2, dft_control%nspins))
789 3592 : ALLOCATE (eigrmat(dft_control%nspins))
790 3592 : ALLOCATE (inv_mat(dft_control%nspins))
791 5388 : ALLOCATE (inv_work(2, dft_control%nspins))
792 6286 : ALLOCATE (mo_derivs_tmp(3, SIZE(mo_derivs)))
793 3592 : ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
794 :
795 : ! Allocate temp matrices for the wavefunction derivatives
796 1796 : DO ispin = 1, dft_control%nspins
797 898 : NULLIFY (tmp_fm_struct, mo_coeff)
798 898 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
799 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
800 898 : ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
801 898 : CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
802 3592 : DO i = 1, 3
803 2694 : CALL cp_fm_create(mo_derivs_tmp(i, ispin), mo_coeff%matrix_struct)
804 3592 : CALL cp_fm_set_all(matrix=mo_derivs_tmp(i, ispin), alpha=0.0_dp)
805 : END DO
806 2694 : DO i = 1, SIZE(op_fm_set, 1)
807 1796 : CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
808 1796 : CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
809 2694 : CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
810 : END DO
811 898 : CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
812 898 : CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
813 2694 : CALL cp_fm_struct_release(tmp_fm_struct)
814 : END DO
815 : ! temp matrices for force calculation
816 898 : IF (calculate_forces) THEN
817 10 : NULLIFY (matrix_s)
818 10 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
819 60 : ALLOCATE (tempmat(2, dft_control%nspins))
820 20 : DO ispin = 1, dft_control%nspins
821 10 : ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
822 10 : CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
823 10 : CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
824 10 : CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
825 20 : CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
826 : END DO
827 : ! integration
828 10 : CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
829 80 : ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
830 70 : ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
831 10 : lsab = MAX(lsab, ldab)
832 50 : DO i = 1, 3
833 180 : ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
834 160 : ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
835 : END DO
836 : END IF
837 :
838 : !Start the MO derivative calculation
839 : !loop over all cell vectors
840 3592 : DO idir = 1, 3
841 2694 : zi(idir) = z_zero
842 2694 : cosmat => efield%cosmat(idir)%matrix
843 2694 : sinmat => efield%sinmat(idir)%matrix
844 : !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
845 : !first step S_berry * C and C^T S_berry C
846 5388 : DO ispin = 1, dft_control%nspins ! spin
847 2694 : IF (mos(ispin)%use_mo_coeff_b) THEN
848 2694 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
849 2694 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
850 : ELSE
851 0 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
852 0 : mo_coeff_tmp(ispin) = mo_coeff
853 : END IF
854 2694 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
855 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
856 2694 : op_fm_set(1, ispin))
857 2694 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
858 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
859 5388 : op_fm_set(2, ispin))
860 : END DO
861 : !second step invert C^T S_berry C
862 2694 : zdet = z_one
863 5388 : DO ispin = 1, dft_control%nspins
864 2694 : CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
865 2694 : CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
866 2694 : CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
867 2694 : CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
868 5388 : zdet = zdet*zdeta
869 : END DO
870 2694 : zi(idir) = zdet**occ
871 2694 : zlog(idir) = AIMAG(LOG(zi(idir)))
872 :
873 2694 : IF (.NOT. just_energy) THEN
874 : !compute the orbital derivative
875 5388 : DO ispin = 1, dft_control%nspins
876 35022 : inv_work(1, ispin)%local_data(:, :) = REAL(inv_mat(ispin)%local_data(:, :), dp)
877 35022 : inv_work(2, ispin)%local_data(:, :) = AIMAG(inv_mat(ispin)%local_data(:, :))
878 2694 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
879 13470 : DO i = 1, 3
880 8082 : focc = hmat(idir, i)
881 : CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
882 8082 : 1.0_dp, mo_derivs_tmp(idir, ispin))
883 : CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
884 10776 : 1.0_dp, mo_derivs_tmp(idir, ispin))
885 : END DO
886 : END DO
887 : END IF
888 :
889 : !compute nuclear forces
890 3592 : IF (calculate_forces) THEN
891 30 : nkind = SIZE(qs_kind_set)
892 30 : natom = SIZE(particle_set)
893 120 : kvec(:) = twopi*cell%h_inv(idir, :)
894 :
895 : ! calculate: C [C^T S_berry C]^(-1) C^T
896 : ! Store this matrix in DBCSR form (only S overlap blocks)
897 60 : DO ispin = 1, dft_control%nspins
898 30 : CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
899 30 : CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
900 30 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
901 : CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
902 30 : opvec(1, ispin))
903 : CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
904 30 : opvec(2, ispin))
905 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
906 30 : matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
907 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
908 90 : matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
909 : END DO
910 :
911 : ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
912 150 : ALLOCATE (basis_set_list(nkind))
913 90 : DO ikind = 1, nkind
914 60 : qs_kind => qs_kind_set(ikind)
915 60 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
916 90 : IF (ASSOCIATED(basis_set_a)) THEN
917 60 : basis_set_list(ikind)%gto_basis_set => basis_set_a
918 : ELSE
919 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
920 : END IF
921 : END DO
922 : !
923 30 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
924 585 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
925 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
926 555 : iatom=iatom, jatom=jatom, r=rab)
927 555 : basis_set_a => basis_set_list(ikind)%gto_basis_set
928 555 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
929 555 : basis_set_b => basis_set_list(jkind)%gto_basis_set
930 555 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
931 : ! basis ikind
932 555 : first_sgfa => basis_set_a%first_sgf
933 555 : la_max => basis_set_a%lmax
934 555 : la_min => basis_set_a%lmin
935 555 : npgfa => basis_set_a%npgf
936 555 : nseta = basis_set_a%nset
937 555 : nsgfa => basis_set_a%nsgf_set
938 555 : rpgfa => basis_set_a%pgf_radius
939 555 : set_radius_a => basis_set_a%set_radius
940 555 : sphi_a => basis_set_a%sphi
941 555 : zeta => basis_set_a%zet
942 : ! basis jkind
943 555 : first_sgfb => basis_set_b%first_sgf
944 555 : lb_max => basis_set_b%lmax
945 555 : lb_min => basis_set_b%lmin
946 555 : npgfb => basis_set_b%npgf
947 555 : nsetb = basis_set_b%nset
948 555 : nsgfb => basis_set_b%nsgf_set
949 555 : rpgfb => basis_set_b%pgf_radius
950 555 : set_radius_b => basis_set_b%set_radius
951 555 : sphi_b => basis_set_b%sphi
952 555 : zetb => basis_set_b%zet
953 :
954 555 : ldsa = SIZE(sphi_a, 1)
955 555 : ldsb = SIZE(sphi_b, 1)
956 555 : ra(:) = pbc(particle_set(iatom)%r(:), cell)
957 2220 : rb(:) = ra + rab
958 555 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
959 :
960 555 : IF (iatom <= jatom) THEN
961 354 : irow = iatom
962 354 : icol = jatom
963 : ELSE
964 201 : irow = jatom
965 201 : icol = iatom
966 : END IF
967 :
968 555 : IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
969 : fab = 1.0_dp*occ
970 : ELSE
971 525 : fab = 2.0_dp*occ
972 : END IF
973 :
974 2220 : DO i = 1, 3
975 454545 : dcost(i, 1)%block = 0.0_dp
976 454545 : dsint(i, 1)%block = 0.0_dp
977 454545 : dcost(i, 2)%block = 0.0_dp
978 455100 : dsint(i, 2)%block = 0.0_dp
979 : END DO
980 :
981 1665 : DO iset = 1, nseta
982 1110 : ncoa = npgfa(iset)*ncoset(la_max(iset))
983 1110 : sgfa = first_sgfa(1, iset)
984 3885 : DO jset = 1, nsetb
985 2220 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
986 1140 : ncob = npgfb(jset)*ncoset(lb_max(jset))
987 1140 : sgfb = first_sgfb(1, jset)
988 : ! Calculate the primitive integrals (da|b)
989 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
990 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
991 1140 : ra, rb, kvec, cosab, sinab, dcosab, dsinab)
992 4560 : DO i = 1, 3
993 : CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
994 : ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
995 : ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
996 4560 : dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
997 : END DO
998 : ! Calculate the primitive integrals (a|db)
999 : CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1000 : la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1001 1140 : rb, ra, kvec, cosab, sinab, dcosab, dsinab)
1002 5670 : DO i = 1, 3
1003 560556 : dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i))
1004 560556 : dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i))
1005 : CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
1006 : ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1007 : ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1008 5640 : dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
1009 : END DO
1010 : END DO
1011 : END DO
1012 555 : forcea = 0.0_dp
1013 555 : forceb = 0.0_dp
1014 1110 : DO ispin = 1, dft_control%nspins
1015 555 : NULLIFY (rblock, iblock)
1016 : CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
1017 555 : row=irow, col=icol, BLOCK=rblock, found=found)
1018 555 : CPASSERT(found)
1019 : CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
1020 555 : row=irow, col=icol, BLOCK=iblock, found=found)
1021 555 : CPASSERT(found)
1022 555 : n1 = SIZE(rblock, 1)
1023 555 : n2 = SIZE(rblock, 2)
1024 555 : CPASSERT(SIZE(iblock, 1) == n1)
1025 555 : CPASSERT(SIZE(iblock, 2) == n2)
1026 555 : CPASSERT(lsab >= n1)
1027 555 : CPASSERT(lsab >= n2)
1028 2220 : IF (iatom <= jatom) THEN
1029 1416 : DO i = 1, 3
1030 : forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
1031 160542 : - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
1032 : forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
1033 160896 : - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
1034 : END DO
1035 : ELSE
1036 804 : DO i = 1, 3
1037 : forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
1038 85023 : - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
1039 : forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
1040 85224 : - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
1041 : END DO
1042 : END IF
1043 : END DO
1044 2250 : DO i = 1, 3
1045 6660 : force_tmp(iatom, :, i) = force_tmp(iatom, :, i) - fab*hmat(i, idir)*forcea(:)
1046 7215 : force_tmp(jatom, :, i) = force_tmp(jatom, :, i) - fab*hmat(i, idir)*forceb(:)
1047 : END DO
1048 : END DO
1049 30 : CALL neighbor_list_iterator_release(nl_iterator)
1050 30 : DEALLOCATE (basis_set_list)
1051 : END IF
1052 : END DO
1053 :
1054 : ! make sure the total normalized polarization is within [-1:1]
1055 3592 : DO idir = 1, 3
1056 2694 : cqi(idir) = rlog(idir) + zlog(idir)
1057 2694 : IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
1058 2694 : IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
1059 : ! now check for log branch
1060 3592 : IF (calculate_forces) THEN
1061 30 : IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN
1062 0 : di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
1063 0 : DO i = 1, 10
1064 0 : cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi
1065 0 : IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
1066 : END DO
1067 : END IF
1068 : END IF
1069 : END DO
1070 3592 : DO idir = 1, 3
1071 2694 : qi(idir) = 0.0_dp
1072 2694 : ci(idir) = 0.0_dp
1073 11674 : DO i = 1, 3
1074 10776 : ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
1075 : END DO
1076 : END DO
1077 :
1078 : ! update the references
1079 898 : IF (calculate_forces) THEN
1080 40 : ener_field = SUM(ci)
1081 : ! check for smoothness of energy surface
1082 130 : IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN
1083 0 : CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
1084 : END IF
1085 10 : efield%field_energy = ener_field
1086 40 : efield%polarisation(:) = cqi(:)
1087 : END IF
1088 :
1089 : ! Energy
1090 898 : ener_field = 0.0_dp
1091 3592 : DO i = 1, 3
1092 3592 : ener_field = ener_field + dfilter(i)*(fieldpol(i) - 2._dp*twopi*ci(i))**2
1093 : END DO
1094 898 : energy%efield = 0.25_dp*omega/twopi*ener_field
1095 :
1096 : ! debugging output
1097 : IF (para_env%is_source()) THEN
1098 898 : iodeb = -1
1099 : IF (iodeb > 0) THEN
1100 : WRITE (iodeb, '(A,T61,F20.10)') " Polarisation Quantum: ", 2._dp*twopi*twopi*hmat(3, 3)
1101 : WRITE (iodeb, '(A,T21,3F20.10)') " Polarisation: ", 2._dp*twopi*ci(1:3)
1102 : WRITE (iodeb, '(A,T21,3F20.10)') " Displacement: ", fieldpol(1:3)
1103 : WRITE (iodeb, '(A,T21,3F20.10)') " E-Field: ", ((fieldpol(i) - 2._dp*twopi*ci(i)), i=1, 3)
1104 : WRITE (iodeb, '(A,T61,F20.10)') " Disp Free Energy:", energy%efield
1105 : END IF
1106 : END IF
1107 :
1108 898 : IF (.NOT. just_energy) THEN
1109 3592 : DO i = 1, 3
1110 3592 : di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i)
1111 : END DO
1112 : ! Add the result to mo_derivativs
1113 1796 : DO ispin = 1, dft_control%nspins
1114 898 : CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_coeff_tmp(ispin))
1115 4490 : DO idir = 1, 3
1116 : CALL cp_fm_scale_and_add(1.0_dp, mo_coeff_tmp(ispin), di(idir), &
1117 3592 : mo_derivs_tmp(idir, ispin))
1118 : END DO
1119 : END DO
1120 1796 : DO ispin = 1, dft_control%nspins
1121 1796 : CALL copy_fm_to_dbcsr(mo_coeff_tmp(ispin), mo_derivs(ispin)%matrix)
1122 : END DO
1123 : END IF
1124 :
1125 898 : IF (calculate_forces) THEN
1126 40 : DO i = 1, 3
1127 100 : DO ia = 1, natom
1128 60 : CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
1129 60 : iatom = atom_of_kind(ia)
1130 450 : force(ikind)%efield(1:3, iatom) = force(ikind)%efield(1:3, iatom) + di(i)*force_tmp(ia, 1:3, i)
1131 : END DO
1132 : END DO
1133 : END IF
1134 :
1135 1796 : DO ispin = 1, dft_control%nspins
1136 898 : CALL cp_cfm_release(eigrmat(ispin))
1137 898 : CALL cp_cfm_release(inv_mat(ispin))
1138 898 : IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
1139 3592 : DO i = 1, 3
1140 3592 : CALL cp_fm_release(mo_derivs_tmp(i, ispin))
1141 : END DO
1142 3592 : DO i = 1, SIZE(op_fm_set, 1)
1143 1796 : CALL cp_fm_release(opvec(i, ispin))
1144 1796 : CALL cp_fm_release(op_fm_set(i, ispin))
1145 2694 : CALL cp_fm_release(inv_work(i, ispin))
1146 : END DO
1147 : END DO
1148 898 : DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
1149 898 : DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
1150 :
1151 898 : IF (calculate_forces) THEN
1152 30 : DO ikind = 1, SIZE(atomic_kind_set)
1153 190 : CALL para_env%sum(force(ikind)%efield)
1154 : END DO
1155 10 : DEALLOCATE (force_tmp)
1156 10 : DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
1157 40 : DO i = 1, 3
1158 30 : DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
1159 40 : DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
1160 : END DO
1161 10 : CALL dbcsr_deallocate_matrix_set(tempmat)
1162 : END IF
1163 898 : CALL timestop(handle)
1164 :
1165 3592 : END SUBROUTINE qs_dispfield_derivatives
1166 :
1167 : ! **************************************************************************************************
1168 : !> \brief ...
1169 : !> \param cos_block ...
1170 : !> \param sin_block ...
1171 : !> \param ncoa ...
1172 : !> \param nsgfa ...
1173 : !> \param sgfa ...
1174 : !> \param sphi_a ...
1175 : !> \param ldsa ...
1176 : !> \param ncob ...
1177 : !> \param nsgfb ...
1178 : !> \param sgfb ...
1179 : !> \param sphi_b ...
1180 : !> \param ldsb ...
1181 : !> \param cosab ...
1182 : !> \param sinab ...
1183 : !> \param ldab ...
1184 : !> \param work ...
1185 : !> \param ldwork ...
1186 : ! **************************************************************************************************
1187 23304 : SUBROUTINE contract_all(cos_block, sin_block, &
1188 46608 : ncoa, nsgfa, sgfa, sphi_a, ldsa, &
1189 46608 : ncob, nsgfb, sgfb, sphi_b, ldsb, &
1190 23304 : cosab, sinab, ldab, work, ldwork)
1191 :
1192 : REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
1193 : INTEGER, INTENT(IN) :: ncoa, nsgfa, sgfa
1194 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
1195 : INTEGER, INTENT(IN) :: ldsa, ncob, nsgfb, sgfb
1196 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
1197 : INTEGER, INTENT(IN) :: ldsb
1198 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
1199 : INTEGER, INTENT(IN) :: ldab
1200 : REAL(dp), DIMENSION(:, :) :: work
1201 : INTEGER, INTENT(IN) :: ldwork
1202 :
1203 : ! Calculate cosine
1204 :
1205 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, cosab(1, 1), ldab, &
1206 23304 : sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1207 :
1208 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1209 23304 : work(1, 1), ldwork, 1.0_dp, cos_block(sgfa, sgfb), SIZE(cos_block, 1))
1210 :
1211 : ! Calculate sine
1212 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sinab(1, 1), ldab, &
1213 23304 : sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1214 :
1215 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1216 23304 : work(1, 1), ldwork, 1.0_dp, sin_block(sgfa, sgfb), SIZE(sin_block, 1))
1217 :
1218 23304 : END SUBROUTINE contract_all
1219 :
1220 : END MODULE qs_efield_berry
|
