1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief g tensor calculation by dfpt
8!>      Initialization of the epr_env, creation of the special neighbor lists
9!>      Perturbation Hamiltonians by application of the p and rxp oprtators to  psi0
10!>      Write output
11!>      Deallocate everything
12!> \note
13!>      The psi0 should be localized
14!>      the Sebastiani method works within the assumption that the orbitals are
15!>      completely contained in the simulation box
16!> \par History
17!>       created 07-2005 [MI]
18!> \author MI
19! **************************************************************************************************
20MODULE qs_linres_epr_utils
21   USE atomic_kind_types,               ONLY: atomic_kind_type
22   USE cell_types,                      ONLY: cell_type
23   USE cp_control_types,                ONLY: dft_control_type
24   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
25                                              cp_logger_type
26   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
27                                              cp_print_key_unit_nr
28   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
29                                              section_vals_type
30   USE kinds,                           ONLY: dp
31   USE mathconstants,                   ONLY: fourpi,&
32                                              twopi
33   USE particle_types,                  ONLY: particle_type
34   USE physcon,                         ONLY: a_fine,&
35                                              e_gfactor
36   USE pw_env_types,                    ONLY: pw_env_get,&
37                                              pw_env_type
38   USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
39                                              pw_pool_type
40   USE pw_types,                        ONLY: COMPLEXDATA1D,&
41                                              REALDATA3D,&
42                                              REALSPACE,&
43                                              RECIPROCALSPACE,&
44                                              pw_p_type
45   USE qs_environment_types,            ONLY: get_qs_env,&
46                                              qs_environment_type
47   USE qs_kind_types,                   ONLY: qs_kind_type
48   USE qs_linres_types,                 ONLY: deallocate_nablavks_atom_set,&
49                                              epr_env_create,&
50                                              epr_env_type,&
51                                              init_nablavks_atom_set,&
52                                              linres_control_type,&
53                                              nablavks_atom_type,&
54                                              set_epr_env
55   USE qs_matrix_pools,                 ONLY: qs_matrix_pools_type
56   USE qs_mo_types,                     ONLY: mo_set_p_type
57   USE qs_rho_atom_types,               ONLY: deallocate_rho_atom_set
58   USE qs_rho_types,                    ONLY: qs_rho_clear,&
59                                              qs_rho_create,&
60                                              qs_rho_set
61   USE scf_control_types,               ONLY: scf_control_type
62#include "./base/base_uses.f90"
63
64   IMPLICIT NONE
65
66   PRIVATE
67
68   PUBLIC :: epr_env_cleanup, epr_env_init
69
70   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_utils'
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief Initialize the epr environment
76!> \param epr_env ...
77!> \param qs_env ...
78!> \par History
79!>      07.2006 created [MI]
80!> \author MI
81! **************************************************************************************************
82   SUBROUTINE epr_env_init(epr_env, qs_env)
83      !
84      TYPE(epr_env_type)                                 :: epr_env
85      TYPE(qs_environment_type), POINTER                 :: qs_env
86
87      CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_env_init', routineP = moduleN//':'//routineN
88
89      INTEGER                                            :: handle, i_B, idir, ispin, n_mo(2), nao, &
90                                                            natom, nmoloc, nspins, output_unit
91      LOGICAL                                            :: gapw
92      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
93      TYPE(cell_type), POINTER                           :: cell
94      TYPE(cp_logger_type), POINTER                      :: logger
95      TYPE(dft_control_type), POINTER                    :: dft_control
96      TYPE(linres_control_type), POINTER                 :: linres_control
97      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
98      TYPE(nablavks_atom_type), DIMENSION(:), POINTER    :: nablavks_atom_set
99      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
100      TYPE(pw_env_type), POINTER                         :: pw_env
101      TYPE(pw_p_type), DIMENSION(:), POINTER             :: rho_g, rho_r
102      TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
103      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
104      TYPE(qs_matrix_pools_type), POINTER                :: mpools
105      TYPE(scf_control_type), POINTER                    :: scf_control
106      TYPE(section_vals_type), POINTER                   :: lr_section
107
108      CALL timeset(routineN, handle)
109
110      NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, linres_control, scf_control)
111      NULLIFY (logger, mos, mpools, particle_set)
112      NULLIFY (auxbas_pw_pool, pw_env)
113      NULLIFY (nablavks_atom_set)
114
115      n_mo(1:2) = 0
116      nao = 0
117      nmoloc = 0
118
119      logger => cp_get_default_logger()
120      !ionode = logger%para_env%mepos==logger%para_env%source
121      lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
122
123      output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
124                                         extension=".linresLog")
125
126      IF (epr_env%ref_count /= 0) THEN
127         CALL epr_env_cleanup(epr_env)
128      END IF
129
130      IF (output_unit > 0) THEN
131         WRITE (output_unit, "(/,T20,A,/)") "*** Start EPR g tensor calculation ***"
132         WRITE (output_unit, "(T10,A,/)") "Initialization of the EPR environment"
133      ENDIF
134
135      CALL epr_env_create(epr_env)
136
137      CALL get_qs_env(qs_env=qs_env, &
138                      atomic_kind_set=atomic_kind_set, &
139                      qs_kind_set=qs_kind_set, &
140                      cell=cell, &
141                      dft_control=dft_control, &
142                      linres_control=linres_control, &
143                      mos=mos, &
144                      mpools=mpools, &
145                      particle_set=particle_set, &
146                      pw_env=pw_env, &
147                      scf_control=scf_control)
148      !
149      ! Check if restat also psi0 should be restarted
150      !IF(epr_env%restart_epr .AND. scf_control%density_guess/=restart_guess)THEN
151      !   CPABORT("restart_epr requires density_guess=restart")
152      !ENDIF
153      !
154      ! check that the psi0 are localized and you have all the centers
155      CPASSERT(linres_control%localized_psi0)
156      IF (output_unit > 0) THEN
157         WRITE (output_unit, '(A)') &
158            ' To get EPR parameters within PBC you need localized zero order orbitals '
159      ENDIF
160      gapw = dft_control%qs_control%gapw
161      nspins = dft_control%nspins
162      natom = SIZE(particle_set, 1)
163      !
164      ! Conversion factors
165      ! Magical constant twopi/cell%deth just like in NMR shift (basically undo scale_fac in qs_linres_nmr_current.F)
166      epr_env%g_free_factor = -1.0_dp*e_gfactor
167      epr_env%g_zke_factor = e_gfactor*(a_fine)**2
168      epr_env%g_so_factor = (a_fine)**2*(-1.0_dp*e_gfactor - 1.0_dp)/2.0_dp*twopi/cell%deth
169      epr_env%g_so_factor_gapw = (a_fine)**2*(-1.0_dp*e_gfactor - 1.0_dp)/2.0_dp
170      ! * 2 because B_ind = 2 * B_beta
171      epr_env%g_soo_factor = 2.0_dp*fourpi*(a_fine)**2*twopi/cell%deth
172      ! 2 * 2 * 1/4 * e^2 / m * a_0^2 * 2/3 * mu_0 / (omega * 1e-30 )
173      epr_env%g_soo_chicorr_factor = 2.0/3.0_dp*fourpi*(a_fine)**2/cell%deth
174      !
175      ! If the current density on the grid needs to be stored
176      CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
177      !
178      ! Initialize local current density if GAPW calculation
179      IF (gapw) THEN
180         CALL init_nablavks_atom_set(nablavks_atom_set, atomic_kind_set, qs_kind_set, nspins)
181         CALL set_epr_env(epr_env=epr_env, &
182                          nablavks_atom_set=nablavks_atom_set)
183      ENDIF
184      !
185      ! Bind
186      ALLOCATE (epr_env%bind_set(3, 3))
187      DO i_B = 1, 3
188         DO idir = 1, 3
189            NULLIFY (epr_env%bind_set(idir, i_B)%rho, rho_r, rho_g)
190            CALL qs_rho_create(epr_env%bind_set(idir, i_B)%rho)
191            ALLOCATE (rho_r(1), rho_g(1))
192            CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(1)%pw, &
193                                   use_data=REALDATA3D, in_space=REALSPACE)
194            CALL pw_pool_create_pw(auxbas_pw_pool, rho_g(1)%pw, &
195                                   use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
196            CALL qs_rho_set(epr_env%bind_set(idir, i_B)%rho, rho_r=rho_r, rho_g=rho_g)
197         END DO
198      END DO
199
200      ! Nabla_V_ks
201      ALLOCATE (epr_env%nablavks_set(3, dft_control%nspins))
202      DO idir = 1, 3
203         DO ispin = 1, nspins
204            NULLIFY (epr_env%nablavks_set(idir, ispin)%rho, rho_r, rho_g)
205            CALL qs_rho_create(epr_env%nablavks_set(idir, ispin)%rho)
206            ALLOCATE (rho_r(1), rho_g(1))
207            CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(1)%pw, &
208                                   use_data=REALDATA3D, in_space=REALSPACE)
209            CALL pw_pool_create_pw(auxbas_pw_pool, rho_g(1)%pw, &
210                                   use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
211            CALL qs_rho_set(epr_env%nablavks_set(idir, ispin)%rho, &
212                            rho_r=rho_r, rho_g=rho_g)
213         END DO
214      END DO
215
216      ! Initialize the g tensor components
217      ALLOCATE (epr_env%g_total(3, 3))
218      ALLOCATE (epr_env%g_so(3, 3))
219      ALLOCATE (epr_env%g_soo(3, 3))
220      epr_env%g_total = 0.0_dp
221      epr_env%g_zke = 0.0_dp
222      epr_env%g_so = 0.0_dp
223      epr_env%g_soo = 0.0_dp
224
225      CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
226           &                            "PRINT%PROGRAM_RUN_INFO")
227
228      CALL timestop(handle)
229
230   END SUBROUTINE epr_env_init
231
232! **************************************************************************************************
233!> \brief Deallocate the epr environment
234!> \param epr_env ...
235!> \par History
236!>      07.2005 created [MI]
237!> \author MI
238! **************************************************************************************************
239   SUBROUTINE epr_env_cleanup(epr_env)
240
241      TYPE(epr_env_type)                                 :: epr_env
242
243      CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_env_cleanup', &
244         routineP = moduleN//':'//routineN
245
246      INTEGER                                            :: i_B, idir, ispin
247
248      epr_env%ref_count = epr_env%ref_count - 1
249      IF (epr_env%ref_count == 0) THEN
250         ! nablavks_set
251         IF (ASSOCIATED(epr_env%nablavks_set)) THEN
252            DO ispin = 1, SIZE(epr_env%nablavks_set, 2)
253            DO idir = 1, SIZE(epr_env%nablavks_set, 1)
254               CALL qs_rho_clear(epr_env%nablavks_set(idir, ispin)%rho)
255               DEALLOCATE (epr_env%nablavks_set(idir, ispin)%rho)
256            ENDDO
257            ENDDO
258            DEALLOCATE (epr_env%nablavks_set)
259         END IF
260         ! nablavks_atom_set
261         IF (ASSOCIATED(epr_env%nablavks_atom_set)) THEN
262            CALL deallocate_nablavks_atom_set(epr_env%nablavks_atom_set)
263         END IF
264         ! vks_atom_set
265         IF (ASSOCIATED(epr_env%vks_atom_set)) THEN
266            CALL deallocate_rho_atom_set(epr_env%vks_atom_set)
267         END IF
268         ! bind_set
269         IF (ASSOCIATED(epr_env%bind_set)) THEN
270            DO i_B = 1, SIZE(epr_env%bind_set, 2)
271            DO idir = 1, SIZE(epr_env%bind_set, 1)
272               CALL qs_rho_clear(epr_env%bind_set(idir, i_B)%rho)
273               DEALLOCATE (epr_env%bind_set(idir, i_B)%rho)
274            ENDDO
275            ENDDO
276            DEALLOCATE (epr_env%bind_set)
277         END IF
278         ! bind_atom_set
279         IF (ASSOCIATED(epr_env%bind_atom_set)) THEN
280            DEALLOCATE (epr_env%bind_atom_set)
281         END IF
282         ! g_total
283         IF (ASSOCIATED(epr_env%g_total)) THEN
284            DEALLOCATE (epr_env%g_total)
285         END IF
286         ! g_so
287         IF (ASSOCIATED(epr_env%g_so)) THEN
288            DEALLOCATE (epr_env%g_so)
289         END IF
290         ! g_soo
291         IF (ASSOCIATED(epr_env%g_soo)) THEN
292            DEALLOCATE (epr_env%g_soo)
293         END IF
294      END IF ! ref count
295
296   END SUBROUTINE epr_env_cleanup
297
298END MODULE qs_linres_epr_utils
299