1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6MODULE optimize_dmfet_potential
7   USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
8                                              cp_blacs_env_release,&
9                                              cp_blacs_env_type
10   USE cp_control_types,                ONLY: dft_control_type
11   USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr
12   USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
13                                              cp_fm_struct_release,&
14                                              cp_fm_struct_type
15   USE cp_fm_types,                     ONLY: cp_fm_copy_general,&
16                                              cp_fm_create,&
17                                              cp_fm_maxabsval,&
18                                              cp_fm_release,&
19                                              cp_fm_set_all,&
20                                              cp_fm_type
21   USE cp_gemm_interface,               ONLY: cp_gemm
22   USE cp_para_types,                   ONLY: cp_para_env_type
23   USE dbcsr_api,                       ONLY: dbcsr_create,&
24                                              dbcsr_init_p,&
25                                              dbcsr_multiply,&
26                                              dbcsr_p_type,&
27                                              dbcsr_release,&
28                                              dbcsr_trace,&
29                                              dbcsr_type,&
30                                              dbcsr_type_no_symmetry
31   USE embed_types,                     ONLY: opt_dmfet_pot_type
32   USE input_section_types,             ONLY: section_vals_type,&
33                                              section_vals_val_get
34   USE kinds,                           ONLY: dp
35   USE qs_environment_types,            ONLY: get_qs_env,&
36                                              qs_environment_type
37   USE qs_mo_types,                     ONLY: get_mo_set,&
38                                              mo_set_p_type
39#include "./base/base_uses.f90"
40
41   IMPLICIT NONE
42
43   PRIVATE
44
45   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_dmfet_potential'
46
47   PUBLIC :: build_full_dm, prepare_dmfet_opt, release_dmfet_opt, subsys_spin, check_dmfet
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief ...
53!> \param opt_dmfet ...
54!> \param opt_dmfet_section ...
55! **************************************************************************************************
56   SUBROUTINE read_opt_dmfet_section(opt_dmfet, opt_dmfet_section)
57      TYPE(opt_dmfet_pot_type)                           :: opt_dmfet
58      TYPE(section_vals_type), POINTER                   :: opt_dmfet_section
59
60      ! Read keywords
61
62      CALL section_vals_val_get(opt_dmfet_section, "N_ITER", &
63                                i_val=opt_dmfet%n_iter)
64
65      CALL section_vals_val_get(opt_dmfet_section, "TRUST_RAD", &
66                                r_val=opt_dmfet%trust_rad)
67
68      CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_MAX", &
69                                r_val=opt_dmfet%conv_max)
70
71      CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_INT", &
72                                r_val=opt_dmfet%conv_int)
73
74      CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_MAX", &
75                                r_val=opt_dmfet%conv_max_beta)
76
77      CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_INT", &
78                                r_val=opt_dmfet%conv_int_beta)
79
80      CALL section_vals_val_get(opt_dmfet_section, "READ_DMFET_POT", &
81                                l_val=opt_dmfet%read_dmfet_pot)
82
83   END SUBROUTINE read_opt_dmfet_section
84
85! **************************************************************************************************
86!> \brief ...
87!> \param qs_env ...
88!> \return ...
89! **************************************************************************************************
90   FUNCTION subsys_spin(qs_env) RESULT(subsys_open_shell)
91      TYPE(qs_environment_type), POINTER                 :: qs_env
92      LOGICAL                                            :: subsys_open_shell
93
94      TYPE(dft_control_type), POINTER                    :: dft_control
95
96      NULLIFY (dft_control)
97      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
98      subsys_open_shell = .FALSE.
99      IF (dft_control%nspins == 2) subsys_open_shell = .TRUE.
100
101   END FUNCTION subsys_spin
102
103! **************************************************************************************************
104!> \brief ...
105!> \param qs_env ...
106!> \param opt_dmfet ...
107!> \param opt_dmfet_section ...
108! **************************************************************************************************
109   SUBROUTINE prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section)
110      TYPE(qs_environment_type), POINTER                 :: qs_env
111      TYPE(opt_dmfet_pot_type)                           :: opt_dmfet
112      TYPE(section_vals_type), POINTER                   :: opt_dmfet_section
113
114      INTEGER                                            :: diff_size, nao
115      TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
116      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
117      TYPE(cp_fm_type), POINTER                          :: mo_coeff
118      TYPE(cp_para_env_type), POINTER                    :: para_env
119      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
120
121      ! Read the input
122      CALL read_opt_dmfet_section(opt_dmfet, opt_dmfet_section)
123
124      ! Get the orbital coefficients
125      CALL get_qs_env(qs_env=qs_env, mos=mos, para_env=para_env)
126      CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao)
127
128      ! Make cp_fm matrices
129      NULLIFY (blacs_env)
130      CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
131
132      NULLIFY (opt_dmfet%dmfet_pot, opt_dmfet%dm_1, opt_dmfet%dm_2, opt_dmfet%dm_total, opt_dmfet%dm_diff)
133      NULLIFY (fm_struct)
134
135      CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
136                               nrow_global=nao, ncol_global=nao)
137      CALL cp_fm_create(opt_dmfet%dmfet_pot, fm_struct, name="dmfet_pot")
138      CALL cp_fm_create(opt_dmfet%dm_subsys, fm_struct, name="dm_subsys")
139      CALL cp_fm_create(opt_dmfet%dm_total, fm_struct, name="dm_total")
140      CALL cp_fm_create(opt_dmfet%dm_diff, fm_struct, name="dm_diff")
141
142      CALL cp_fm_set_all(opt_dmfet%dmfet_pot, 0.0_dp)
143      CALL cp_fm_set_all(opt_dmfet%dm_subsys, 0.0_dp)
144      CALL cp_fm_set_all(opt_dmfet%dm_total, 0.0_dp)
145      CALL cp_fm_set_all(opt_dmfet%dm_diff, 0.0_dp)
146
147      ! Beta spin counterparts
148      IF (opt_dmfet%open_shell_embed) THEN
149         CALL cp_fm_create(opt_dmfet%dmfet_pot_beta, fm_struct, name="dmfet_pot_beta")
150         CALL cp_fm_create(opt_dmfet%dm_subsys_beta, fm_struct, name="dm_subsys_beta")
151         CALL cp_fm_create(opt_dmfet%dm_total_beta, fm_struct, name="dm_total_beta")
152         CALL cp_fm_create(opt_dmfet%dm_diff_beta, fm_struct, name="dm_diff_beta")
153
154         CALL cp_fm_set_all(opt_dmfet%dmfet_pot_beta, 0.0_dp)
155         CALL cp_fm_set_all(opt_dmfet%dm_subsys_beta, 0.0_dp)
156         CALL cp_fm_set_all(opt_dmfet%dm_total_beta, 0.0_dp)
157         CALL cp_fm_set_all(opt_dmfet%dm_diff_beta, 0.0_dp)
158      ENDIF
159
160      ! Release structure and context
161      CALL cp_fm_struct_release(fm_struct)
162      CALL cp_blacs_env_release(blacs_env)
163
164      ! Array to store functional values
165      ALLOCATE (opt_dmfet%w_func(opt_dmfet%n_iter))
166      opt_dmfet%w_func = 0.0_dp
167
168      ! Allocate max_diff and int_diff
169      diff_size = 1
170      IF (opt_dmfet%open_shell_embed) diff_size = 2
171      ALLOCATE (opt_dmfet%max_diff(diff_size))
172      ALLOCATE (opt_dmfet%int_diff(diff_size))
173
174   END SUBROUTINE prepare_dmfet_opt
175
176! **************************************************************************************************
177!> \brief ...
178!> \param opt_dmfet ...
179! **************************************************************************************************
180   SUBROUTINE release_dmfet_opt(opt_dmfet)
181      TYPE(opt_dmfet_pot_type)                           :: opt_dmfet
182
183      CALL cp_fm_release(opt_dmfet%dmfet_pot)
184      CALL cp_fm_release(opt_dmfet%dm_subsys)
185      CALL cp_fm_release(opt_dmfet%dm_total)
186      CALL cp_fm_release(opt_dmfet%dm_diff)
187
188      IF (opt_dmfet%open_shell_embed) THEN
189         CALL cp_fm_release(opt_dmfet%dmfet_pot_beta)
190         CALL cp_fm_release(opt_dmfet%dm_subsys_beta)
191         CALL cp_fm_release(opt_dmfet%dm_total_beta)
192         CALL cp_fm_release(opt_dmfet%dm_diff_beta)
193      ENDIF
194
195      DEALLOCATE (opt_dmfet%w_func)
196      DEALLOCATE (opt_dmfet%max_diff)
197      DEALLOCATE (opt_dmfet%int_diff)
198      DEALLOCATE (opt_dmfet%all_nspins)
199
200   END SUBROUTINE release_dmfet_opt
201
202! **************************************************************************************************
203!> \brief Builds density matrices from MO coefficients in full matrix format
204!> \param qs_env ...
205!> \param dm ...
206!> \param open_shell  Subsystem is open shell
207!> \param open_shell_embed  Embedding is open shell
208!> \param dm_beta ...
209! **************************************************************************************************
210   SUBROUTINE build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta)
211      TYPE(qs_environment_type), POINTER                 :: qs_env
212      TYPE(cp_fm_type), POINTER                          :: dm
213      LOGICAL                                            :: open_shell, open_shell_embed
214      TYPE(cp_fm_type), POINTER                          :: dm_beta
215
216      CHARACTER(LEN=*), PARAMETER :: routineN = 'build_full_dm', routineP = moduleN//':'//routineN
217
218      INTEGER                                            :: handle, homo, nao
219      REAL(KIND=dp)                                      :: coeff
220      TYPE(cp_fm_type), POINTER                          :: mo_coeff
221      TYPE(cp_para_env_type), POINTER                    :: para_env
222      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
223
224      CALL timeset(routineN, handle)
225
226      coeff = 2.0_dp
227      IF (open_shell_embed) coeff = 1.0_dp
228
229      ! Get the orbital coefficients
230      CALL get_qs_env(qs_env=qs_env, mos=mos)
231      CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao, homo=homo)
232
233      ! Build the density matrix
234      CALL cp_gemm(transa="N", transb="T", m=nao, n=nao, &
235                   k=homo, alpha=coeff, &
236                   matrix_a=mo_coeff, matrix_b=mo_coeff, &
237                   beta=0.0_dp, matrix_c=dm)
238
239      ! Open shell
240      IF (open_shell) THEN
241         CALL get_mo_set(mo_set=mos(2)%mo_set, mo_coeff=mo_coeff, nao=nao, homo=homo)
242
243         ! Build the density matrix
244         CALL cp_gemm(transa="N", transb="T", m=nao, n=nao, &
245                      k=homo, alpha=coeff, &
246                      matrix_a=mo_coeff, matrix_b=mo_coeff, &
247                      beta=0.0_dp, matrix_c=dm_beta)
248      ENDIF
249
250      ! If embedding is open shell, but subsystem is not, copy alpha-spin DM intou beta-spin DM
251      IF (open_shell_embed .AND. (.NOT. open_shell)) THEN
252         CALL get_qs_env(qs_env=qs_env, para_env=para_env)
253         CALL cp_fm_copy_general(dm, dm_beta, para_env)
254
255      ENDIF
256
257      CALL timestop(handle)
258
259   END SUBROUTINE build_full_dm
260
261! **************************************************************************************************
262!> \brief ...
263!> \param opt_dmfet ...
264!> \param qs_env ...
265! **************************************************************************************************
266   SUBROUTINE check_dmfet(opt_dmfet, qs_env)
267      TYPE(opt_dmfet_pot_type)                           :: opt_dmfet
268      TYPE(qs_environment_type), POINTER                 :: qs_env
269
270      REAL(KIND=dp)                                      :: max_diff, max_diff_beta, trace
271      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
272      TYPE(dbcsr_type), POINTER                          :: diff_dbcsr, dm_s
273
274      CALL get_qs_env(qs_env, matrix_s=matrix_s)
275
276      NULLIFY (diff_dbcsr)
277      CALL dbcsr_init_p(diff_dbcsr)
278      CALL dbcsr_create(matrix=diff_dbcsr, &
279                        template=matrix_s(1)%matrix, &
280                        matrix_type=dbcsr_type_no_symmetry)
281
282      NULLIFY (dm_s)
283      CALL dbcsr_init_p(dm_s)
284      CALL dbcsr_create(matrix=dm_s, &
285                        template=matrix_s(1)%matrix, &
286                        matrix_type=dbcsr_type_no_symmetry)
287
288      CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff, diff_dbcsr, keep_sparsity=.FALSE.)
289
290      CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, &
291                          0.0_dp, dm_s)
292
293      CALL dbcsr_trace(dm_s, trace)
294
295      IF (opt_dmfet%open_shell_embed) THEN
296         CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff_beta, diff_dbcsr, keep_sparsity=.FALSE.)
297
298         CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, &
299                             0.0_dp, dm_s)
300
301         CALL dbcsr_trace(dm_s, trace)
302
303      ENDIF
304
305      ! Release dbcsr objects
306      CALL dbcsr_release(diff_dbcsr)
307      CALL dbcsr_release(dm_s)
308
309      ! Find the absolute maximum value of the DM difference
310      CALL cp_fm_maxabsval(opt_dmfet%dm_diff, max_diff)
311      IF (opt_dmfet%open_shell_embed) THEN
312         CALL cp_fm_maxabsval(opt_dmfet%dm_diff_beta, max_diff_beta)
313      ENDIF
314
315   END SUBROUTINE check_dmfet
316
317END MODULE optimize_dmfet_potential
318