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