1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Environment storing all data that is needed in order to call the DFT
8!>        driver of the PEXSI library with data from the linear scaling quickstep
9!>        SCF environment, mainly parameters and intermediate data for the matrix
10!>        conversion between DBCSR and CSR format.
11!> \par History
12!>       2014.11 created [Patrick Seewald]
13!> \author Patrick Seewald
14! **************************************************************************************************
15
16MODULE pexsi_types
17   USE bibliography,                    ONLY: Lin2009,&
18                                              Lin2013,&
19                                              cite_reference
20   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
21                                              cp_logger_get_default_unit_nr,&
22                                              cp_logger_type
23   USE dbcsr_api,                       ONLY: dbcsr_csr_type,&
24                                              dbcsr_p_type,&
25                                              dbcsr_scale,&
26                                              dbcsr_type
27   USE kinds,                           ONLY: dp,&
28                                              int_8
29   USE message_passing,                 ONLY: mp_dims_create,&
30                                              mp_environ
31   USE pexsi_interface,                 ONLY: cp_pexsi_get_options,&
32                                              cp_pexsi_options,&
33                                              cp_pexsi_plan_finalize,&
34                                              cp_pexsi_plan_initialize,&
35                                              cp_pexsi_set_options
36#include "./base/base_uses.f90"
37
38   IMPLICIT NONE
39
40   PRIVATE
41
42   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pexsi_types'
43
44   LOGICAL, PARAMETER, PRIVATE          :: careful_mod = .FALSE.
45
46   INTEGER, PARAMETER, PUBLIC           :: cp2k_to_pexsi = 1, pexsi_to_cp2k = 2
47
48   PUBLIC :: lib_pexsi_env, lib_pexsi_init, lib_pexsi_finalize, &
49             convert_nspin_cp2k_pexsi
50
51! **************************************************************************************************
52!> \brief All PEXSI related data
53!> \param options                       PEXSI options
54!> \param plan                          PEXSI plan
55!> \param mp_group                      message-passing group ID
56!> \param mp_dims                       dimensions of the MPI cartesian grid used
57!>                                      for PEXSI
58!> \param num_ranks_per_pole            number of MPI ranks per pole in PEXSI
59!> \param kTS                           entropic energy contribution
60!> \param matrix_w                      energy-weighted density matrix as needed
61!>                                      for the forces
62!> \param csr_mat                       intermediate matrices in CSR format
63!> \param dbcsr_template_matrix_sym     Symmetric template matrix fixing DBCSR
64!>                                      sparsity pattern
65!> \param dbcsr_template_matrix_nonsym  Nonsymmetric template matrix fixing
66!>                                      DBCSR sparsity pattern
67!> \param csr_sparsity                  DBCSR matrix defining CSR sparsity
68!> \param csr_screening                 whether distance screening should be
69!>                                      applied to CSR matrices
70!> \param max_ev_vector                 eigenvector corresponding to the largest
71!>                                      energy eigenvalue,
72!>                                      returned by the Arnoldi method used to
73!>                                      determine the spectral radius deltaE
74!> \param nspin                         number of spins
75!> \param do_adaptive_tol_nel           Whether or not to use adaptive threshold
76!>                                      for PEXSI convergence
77!> \param adaptive_nel_alpha            constants for adaptive thresholding
78!> \param adaptive_nel_beta             ...
79!> \param tol_nel_initial               Initial convergence threshold (in number of electrons)
80!> \param tol_nel_target                Target convergence threshold (in number of electrons)
81!> \par History
82!>      11.2014 created [Patrick Seewald]
83!> \author Patrick Seewald
84! **************************************************************************************************
85   TYPE lib_pexsi_env
86      TYPE(dbcsr_type)                      :: dbcsr_template_matrix_sym, &
87                                               dbcsr_template_matrix_nonsym
88      TYPE(dbcsr_csr_type)                           :: csr_mat_p, csr_mat_ks, csr_mat_s, &
89                                                        csr_mat_E, csr_mat_F
90      TYPE(cp_pexsi_options)                   :: options
91      REAL(KIND=dp), DIMENSION(:), POINTER     :: kTS => NULL()
92
93      TYPE(dbcsr_p_type), DIMENSION(:), &
94         POINTER                                :: matrix_w => NULL()
95
96      INTEGER(kind=int_8)                      :: plan
97      INTEGER                                  :: nspin, num_ranks_per_pole, mp_group
98
99      TYPE(dbcsr_type), DIMENSION(:), &
100         POINTER                                :: max_ev_vector
101      TYPE(dbcsr_type)                      :: csr_sparsity
102      INTEGER, DIMENSION(2)                    :: mp_dims
103
104      LOGICAL                                  :: csr_screening, do_adaptive_tol_nel
105      REAL(KIND=dp)                            :: adaptive_nel_alpha, adaptive_nel_beta, &
106                                                  tol_nel_initial, tol_nel_target
107   END TYPE lib_pexsi_env
108
109CONTAINS
110
111! **************************************************************************************************
112!> \brief Initialize PEXSI
113!> \param pexsi_env All data needed by PEXSI
114!> \param mp_group message-passing group ID
115!> \param nspin number of spins
116!> \par History
117!>      11.2014 created [Patrick Seewald]
118!> \author Patrick Seewald
119! **************************************************************************************************
120   SUBROUTINE lib_pexsi_init(pexsi_env, mp_group, nspin)
121      TYPE(lib_pexsi_env), INTENT(INOUT)                 :: pexsi_env
122      INTEGER, INTENT(IN)                                :: mp_group, nspin
123
124      CHARACTER(len=*), PARAMETER :: routineN = 'lib_pexsi_init', routineP = moduleN//':'//routineN
125
126      INTEGER                                            :: handle, ispin, mynode, npSymbFact, &
127                                                            numnodes, unit_nr
128      TYPE(cp_logger_type), POINTER                      :: logger
129
130      logger => cp_get_default_logger()
131      IF (logger%para_env%ionode) THEN
132         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
133      ELSE
134         unit_nr = -1
135      ENDIF
136
137      CALL timeset(routineN, handle)
138
139      CALL cite_reference(Lin2009)
140      CALL cite_reference(Lin2013)
141
142      pexsi_env%mp_group = mp_group
143      CALL mp_environ(numnodes, mynode, pexsi_env%mp_group)
144
145      ! Use all nodes available per pole by default or if the user tries to use
146      ! more nodes than MPI size
147      IF ((pexsi_env%num_ranks_per_pole .GT. numnodes) &
148          .OR. (pexsi_env%num_ranks_per_pole .EQ. 0)) THEN
149         pexsi_env%num_ranks_per_pole = numnodes
150      ENDIF
151
152      ! Find num_ranks_per_pole from user input MIN_RANKS_PER_POLE s.t. num_ranks_per_pole
153      ! is the smallest number greater or equal to min_ranks_per_pole that divides
154      ! MPI size without remainder.
155      DO WHILE (MOD(numnodes, pexsi_env%num_ranks_per_pole) .NE. 0)
156         pexsi_env%num_ranks_per_pole = pexsi_env%num_ranks_per_pole + 1
157      ENDDO
158
159      CALL cp_pexsi_get_options(pexsi_env%options, npSymbFact=npSymbFact)
160      IF ((npSymbFact .GT. pexsi_env%num_ranks_per_pole) .OR. (npSymbFact .EQ. 0)) THEN
161         ! Use maximum possible number of ranks for symbolic factorization
162         CALL cp_pexsi_set_options(pexsi_env%options, npSymbFact=pexsi_env%num_ranks_per_pole)
163      ENDIF
164
165      ! Create dimensions for MPI cartesian grid for PEXSI
166      pexsi_env%mp_dims = 0
167      CALL mp_dims_create(pexsi_env%num_ranks_per_pole, pexsi_env%mp_dims)
168
169      ! allocations with nspin
170      pexsi_env%nspin = nspin
171      ALLOCATE (pexsi_env%kTS(nspin))
172      pexsi_env%kTS(:) = 0.0_dp
173      ALLOCATE (pexsi_env%max_ev_vector(nspin))
174      ALLOCATE (pexsi_env%matrix_w(nspin))
175      DO ispin = 1, pexsi_env%nspin
176         ALLOCATE (pexsi_env%matrix_w(ispin)%matrix)
177      ENDDO
178
179      ! Initialize PEXSI
180      pexsi_env%plan = cp_pexsi_plan_initialize(pexsi_env%mp_group, pexsi_env%mp_dims(1), &
181                                                pexsi_env%mp_dims(2), mynode)
182
183      pexsi_env%do_adaptive_tol_nel = .FALSE.
184
185      ! Print PEXSI infos
186      IF (unit_nr > 0) CALL print_pexsi_info(pexsi_env, unit_nr)
187
188      CALL timestop(handle)
189   END SUBROUTINE lib_pexsi_init
190
191! **************************************************************************************************
192!> \brief Release all PEXSI data
193!> \param pexsi_env ...
194!> \par History
195!>      11.2014 created [Patrick Seewald]
196!> \author Patrick Seewald
197! **************************************************************************************************
198   SUBROUTINE lib_pexsi_finalize(pexsi_env)
199      TYPE(lib_pexsi_env), INTENT(INOUT)                 :: pexsi_env
200
201      CHARACTER(len=*), PARAMETER :: routineN = 'lib_pexsi_finalize', &
202         routineP = moduleN//':'//routineN
203
204      INTEGER                                            :: handle, ispin
205
206      CALL timeset(routineN, handle)
207      CALL cp_pexsi_plan_finalize(pexsi_env%plan)
208      DEALLOCATE (pexsi_env%kTS)
209      DEALLOCATE (pexsi_env%max_ev_vector)
210      DO ispin = 1, pexsi_env%nspin
211         DEALLOCATE (pexsi_env%matrix_w(ispin)%matrix)
212      ENDDO
213      DEALLOCATE (pexsi_env%matrix_w)
214      CALL timestop(handle)
215   END SUBROUTINE lib_pexsi_finalize
216
217! **************************************************************************************************
218!> \brief Scale various quantities with factors of 2. This converts spin restricted
219!>        DFT calculations (PEXSI) to the unrestricted case (as is the case where
220!>        the density matrix method is called in the linear scaling code).
221!> \param[in] direction ...
222!> \param[in,out] numElectron ...
223!> \param matrix_p ...
224!> \param matrix_w ...
225!> \param kTS ...
226!> \par History
227!>      01.2015 created [Patrick Seewald]
228!> \author Patrick Seewald
229! **************************************************************************************************
230   SUBROUTINE convert_nspin_cp2k_pexsi(direction, numElectron, matrix_p, matrix_w, kTS)
231      INTEGER, INTENT(IN)                                :: direction
232      REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: numElectron
233      TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL          :: matrix_p
234      TYPE(dbcsr_p_type), INTENT(INOUT), OPTIONAL        :: matrix_w
235      REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: kTS
236
237      CHARACTER(len=*), PARAMETER :: routineN = 'convert_nspin_cp2k_pexsi', &
238         routineP = moduleN//':'//routineN
239
240      INTEGER                                            :: handle
241      REAL(KIND=dp)                                      :: scaling
242
243      CALL timeset(routineN, handle)
244
245      SELECT CASE (direction)
246      CASE (cp2k_to_pexsi)
247         scaling = 2.0_dp
248      CASE (pexsi_to_cp2k)
249         scaling = 0.5_dp
250      END SELECT
251
252      IF (PRESENT(numElectron)) numElectron = scaling*numElectron
253      IF (PRESENT(matrix_p)) CALL dbcsr_scale(matrix_p, scaling)
254      IF (PRESENT(matrix_w)) CALL dbcsr_scale(matrix_w%matrix, scaling)
255      IF (PRESENT(kTS)) kTS = scaling*kTS
256
257      CALL timestop(handle)
258   END SUBROUTINE convert_nspin_cp2k_pexsi
259
260! **************************************************************************************************
261!> \brief Print relevant options of PEXSI
262!> \param pexsi_env ...
263!> \param unit_nr ...
264! **************************************************************************************************
265   SUBROUTINE print_pexsi_info(pexsi_env, unit_nr)
266      TYPE(lib_pexsi_env), INTENT(IN)                    :: pexsi_env
267      INTEGER, INTENT(IN)                                :: unit_nr
268
269      INTEGER                                            :: mynode, npSymbFact, numnodes, numPole, &
270                                                            ordering, rowOrdering
271      REAL(KIND=dp)                                      :: gap, muInertiaExpansion, &
272                                                            muInertiaTolerance, muMax0, muMin0, &
273                                                            muPEXSISafeGuard, temperature
274
275      CALL mp_environ(numnodes, mynode, pexsi_env%mp_group)
276
277      CALL cp_pexsi_get_options(pexsi_env%options, temperature=temperature, gap=gap, &
278                                numPole=numPole, muMin0=muMin0, muMax0=muMax0, muInertiaTolerance= &
279                                muInertiaTolerance, muInertiaExpansion=muInertiaExpansion, &
280                                muPEXSISafeGuard=muPEXSISafeGuard, ordering=ordering, rowOrdering=rowOrdering, &
281                                npSymbFact=npSymbFact)
282
283      WRITE (unit_nr, '(/A)') " PEXSI| Initialized with parameters"
284      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Electronic temperature", temperature
285      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Spectral gap", gap
286      WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of poles", numPole
287      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Target tolerance in number of electrons", &
288         pexsi_env%tol_nel_target
289      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Convergence tolerance for inertia counting in mu", &
290         muInertiaTolerance
291      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Restart tolerance for inertia counting in mu", &
292         muPEXSISafeGuard
293      WRITE (unit_nr, '(A,T61,E20.3)') " PEXSI|   Expansion of mu interval for inertia counting", &
294         muInertiaExpansion
295
296      WRITE (unit_nr, '(/A)') " PEXSI| Parallelization scheme"
297      WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of poles processed in parallel", &
298         numnodes/pexsi_env%num_ranks_per_pole
299      WRITE (unit_nr, '(A,T61,I20)') " PEXSI|   Number of processors per pole", &
300         pexsi_env%num_ranks_per_pole
301      WRITE (unit_nr, '(A,T71,I5,T76,I5)') " PEXSI|   Process grid dimensions", &
302         pexsi_env%mp_dims(1), pexsi_env%mp_dims(2)
303      SELECT CASE (ordering)
304      CASE (0)
305         WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "parallel"
306      CASE (1)
307         WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "sequential"
308      CASE (2)
309         WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Ordering strategy", "multiple minimum degree"
310      END SELECT
311      SELECT CASE (rowOrdering)
312      CASE (0)
313         WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Row permutation strategy", "no row permutation"
314      CASE (1)
315         WRITE (unit_nr, '(A,T61,A20)') " PEXSI|   Row permutation strategy", "make diagonal entry larger than off diagonal"
316      END SELECT
317      IF (ordering .EQ. 0) WRITE (unit_nr, '(A,T61,I20/)') &
318         " PEXSI|   Number of processors for symbolic factorization", npSymbFact
319
320   END SUBROUTINE print_pexsi_info
321END MODULE pexsi_types
322