1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7MODULE qs_gapw_densities
8   USE atomic_kind_types,               ONLY: atomic_kind_type,&
9                                              get_atomic_kind
10   USE cp_control_types,                ONLY: dft_control_type,&
11                                              gapw_control_type
12   USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
13   USE cp_para_types,                   ONLY: cp_para_env_type
14   USE kinds,                           ONLY: dp
15   USE message_passing,                 ONLY: mp_sum
16   USE qs_charges_types,                ONLY: qs_charges_type
17   USE qs_environment_types,            ONLY: get_qs_env,&
18                                              qs_environment_type
19   USE qs_kind_types,                   ONLY: get_qs_kind,&
20                                              qs_kind_type
21   USE qs_local_rho_types,              ONLY: local_rho_type
22   USE qs_rho0_ggrid,                   ONLY: put_rho0_on_grid
23   USE qs_rho0_methods,                 ONLY: calculate_rho0_atom
24   USE qs_rho0_types,                   ONLY: rho0_atom_type,&
25                                              rho0_mpole_type
26   USE qs_rho_atom_methods,             ONLY: calculate_rho_atom
27   USE qs_rho_atom_types,               ONLY: rho_atom_type
28#include "./base/base_uses.f90"
29
30   IMPLICIT NONE
31
32   PRIVATE
33
34   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities'
35
36   PUBLIC :: prepare_gapw_den
37
38CONTAINS
39
40! **************************************************************************************************
41!> \brief ...
42!> \param qs_env ...
43!> \param local_rho_set ...
44!> \param do_rho0 ...
45!> \param kind_set_external can be provided to use different projectors/grids/basis than the default
46! **************************************************************************************************
47   SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
48
49      TYPE(qs_environment_type), POINTER                 :: qs_env
50      TYPE(local_rho_type), OPTIONAL, POINTER            :: local_rho_set
51      LOGICAL, INTENT(IN), OPTIONAL                      :: do_rho0
52      TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
53         POINTER                                         :: kind_set_external
54
55      CHARACTER(len=*), PARAMETER :: routineN = 'prepare_gapw_den', &
56         routineP = moduleN//':'//routineN
57
58      INTEGER                                            :: handle, ikind, ispin, natom, nspins, &
59                                                            output_unit
60      INTEGER, DIMENSION(:), POINTER                     :: atom_list
61      LOGICAL                                            :: extern, my_do_rho0, paw_atom
62      REAL(dp)                                           :: rho0_h_tot, tot_rs_int
63      REAL(dp), DIMENSION(:), POINTER                    :: rho1_h_tot, rho1_s_tot
64      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
65      TYPE(cp_para_env_type), POINTER                    :: para_env
66      TYPE(dft_control_type), POINTER                    :: dft_control
67      TYPE(gapw_control_type), POINTER                   :: gapw_control
68      TYPE(qs_charges_type), POINTER                     :: qs_charges
69      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: my_kind_set
70      TYPE(rho0_atom_type), DIMENSION(:), POINTER        :: rho0_atom_set
71      TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
72      TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
73
74      CALL timeset(routineN, handle)
75
76      NULLIFY (atomic_kind_set)
77      NULLIFY (my_kind_set)
78      NULLIFY (dft_control)
79      NULLIFY (gapw_control)
80      NULLIFY (para_env)
81      NULLIFY (atom_list)
82      NULLIFY (rho0_mpole)
83      NULLIFY (qs_charges)
84      NULLIFY (rho1_h_tot, rho1_s_tot)
85      NULLIFY (rho_atom_set)
86      NULLIFY (rho0_atom_set)
87
88      my_do_rho0 = .TRUE.
89      IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
90
91      output_unit = cp_logger_get_default_io_unit()
92
93      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
94                      para_env=para_env, &
95                      qs_charges=qs_charges, &
96                      qs_kind_set=my_kind_set, &
97                      atomic_kind_set=atomic_kind_set, &
98                      rho0_mpole=rho0_mpole, &
99                      rho_atom_set=rho_atom_set, &
100                      rho0_atom_set=rho0_atom_set)
101
102      gapw_control => dft_control%qs_control%gapw_control
103
104      IF (PRESENT(local_rho_set)) THEN
105         rho_atom_set => local_rho_set%rho_atom_set
106         IF (my_do_rho0) THEN
107            rho0_mpole => local_rho_set%rho0_mpole
108            rho0_atom_set => local_rho_set%rho0_atom_set
109         END IF
110      END IF
111
112      extern = .FALSE.
113      IF (PRESENT(kind_set_external)) THEN
114         CPASSERT(ASSOCIATED(kind_set_external))
115         my_kind_set => kind_set_external
116         extern = .TRUE.
117      END IF
118
119      nspins = dft_control%nspins
120
121      rho0_h_tot = 0.0_dp
122      ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
123      rho1_h_tot = 0.0_dp
124      rho1_s_tot = 0.0_dp
125
126      DO ikind = 1, SIZE(atomic_kind_set)
127         CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
128         CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
129
130!     Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
131         IF (paw_atom) &
132            CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
133                                    atom_list, natom, nspins, rho1_h_tot, rho1_s_tot)
134
135!     Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
136         IF (my_do_rho0) &
137            CALL calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, rho0_mpole, &
138                                     atom_list, natom, ikind, my_kind_set(ikind), rho0_h_tot)
139
140      ENDDO
141
142      !Do not mess with charges if using a non-default kind_set
143      IF (.NOT. extern) THEN
144         CALL mp_sum(rho1_h_tot, para_env%group)
145         CALL mp_sum(rho1_s_tot, para_env%group)
146         DO ispin = 1, nspins
147            qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
148            qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
149         END DO
150
151         IF (my_do_rho0) THEN
152            rho0_mpole%total_rho0_h = -rho0_h_tot
153!         Put the rho0_soft on the global grid
154            CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
155            IF (ABS(rho0_h_tot) .GE. 1.0E-5_dp) THEN
156               IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) .GT. 1.0E-3_dp) THEN
157                  IF (output_unit > 0) THEN
158                     WRITE (output_unit, '(/,72("*"))')
159                     WRITE (output_unit, '(T2,A,T66,1E20.8)') &
160                        "WARNING: rho0 calculated on the local grid is  :", -rho0_h_tot, &
161                        "         rho0 calculated on the global grid is :", tot_rs_int
162                     WRITE (output_unit, '(T2,A)') &
163                        "         bad integration"
164                     WRITE (output_unit, '(72("*"),/)')
165                  END IF
166               END IF
167            END IF
168            qs_charges%total_rho0_soft_rspace = tot_rs_int
169            qs_charges%total_rho0_hard_lebedev = rho0_h_tot
170         ELSE
171            qs_charges%total_rho0_hard_lebedev = 0.0_dp
172         END IF
173      END IF
174
175      DEALLOCATE (rho1_h_tot, rho1_s_tot)
176
177      CALL timestop(handle)
178
179   END SUBROUTINE prepare_gapw_den
180
181END MODULE qs_gapw_densities
182