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