1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Set of routines to apply restraints to the KS hamiltonian
8! **************************************************************************************************
9MODULE qs_ks_apply_restraints
10   USE cp_control_types,                ONLY: dft_control_type
11   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
12                                              copy_fm_to_dbcsr
13   USE cp_fm_types,                     ONLY: cp_fm_create,&
14                                              cp_fm_p_type,&
15                                              cp_fm_type
16   USE cp_para_types,                   ONLY: cp_para_env_type
17   USE dbcsr_api,                       ONLY: dbcsr_p_type
18   USE input_constants,                 ONLY: cdft_charge_constraint,&
19                                              outer_scf_becke_constraint,&
20                                              outer_scf_hirshfeld_constraint
21   USE kinds,                           ONLY: dp
22   USE mulliken,                        ONLY: mulliken_restraint
23   USE pw_methods,                      ONLY: pw_scale
24   USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
25                                              pw_pool_type
26   USE pw_types,                        ONLY: REALDATA3D,&
27                                              REALSPACE
28   USE qs_cdft_methods,                 ONLY: becke_constraint,&
29                                              hirshfeld_constraint
30   USE qs_cdft_types,                   ONLY: cdft_control_type
31   USE qs_energy_types,                 ONLY: qs_energy_type
32   USE qs_environment_types,            ONLY: get_qs_env,&
33                                              qs_environment_type
34   USE qs_mo_types,                     ONLY: get_mo_set,&
35                                              mo_set_p_type
36   USE qs_rho_types,                    ONLY: qs_rho_get,&
37                                              qs_rho_type
38   USE s_square_methods,                ONLY: s2_restraint
39#include "./base/base_uses.f90"
40
41   IMPLICIT NONE
42
43   PRIVATE
44
45   LOGICAL, PARAMETER :: debug_this_module = .TRUE.
46   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_apply_restraints'
47
48   PUBLIC :: qs_ks_mulliken_restraint, qs_ks_s2_restraint
49   PUBLIC :: qs_ks_cdft_constraint
50
51CONTAINS
52
53! **************************************************************************************************
54!> \brief Apply a CDFT constraint
55!> \param qs_env the qs_env where to apply the constraint
56!> \param auxbas_pw_pool the pool that owns the real space grid where the CDFT potential is defined
57!> \param calculate_forces if forces should be calculated
58!> \param cdft_control the CDFT control type
59! **************************************************************************************************
60   SUBROUTINE qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
61      TYPE(qs_environment_type), POINTER                 :: qs_env
62      TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
63      LOGICAL, INTENT(in)                                :: calculate_forces
64      TYPE(cdft_control_type), POINTER                   :: cdft_control
65
66      CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_cdft_constraint', &
67         routineP = moduleN//':'//routineN
68
69      INTEGER                                            :: iatom, igroup, natom
70      LOGICAL                                            :: do_kpoints
71      REAL(KIND=dp)                                      :: inv_vol
72      TYPE(dft_control_type), POINTER                    :: dft_control
73
74      NULLIFY (dft_control)
75      CALL get_qs_env(qs_env, dft_control=dft_control)
76      IF (dft_control%qs_control%cdft) THEN
77         cdft_control => dft_control%qs_control%cdft_control
78         ! Test no k-points
79         CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
80         IF (do_kpoints) CPABORT("CDFT constraints with k-points not supported.")
81
82         SELECT CASE (cdft_control%type)
83         CASE (outer_scf_becke_constraint, outer_scf_hirshfeld_constraint)
84            IF (cdft_control%need_pot) THEN
85               ! First SCF iteraration => allocate storage
86               DO igroup = 1, SIZE(cdft_control%group)
87                  CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%weight%pw, &
88                                         use_data=REALDATA3D, in_space=REALSPACE)
89                  ! Sanity check
90                  IF (cdft_control%group(igroup)%constraint_type /= cdft_charge_constraint &
91                      .AND. dft_control%nspins == 1) &
92                     CALL cp_abort(__LOCATION__, &
93                                   "Spin constraints require a spin polarized calculation.")
94               END DO
95               IF (cdft_control%atomic_charges) THEN
96                  IF (.NOT. ASSOCIATED(cdft_control%charge)) &
97                     ALLOCATE (cdft_control%charge(cdft_control%natoms))
98                  DO iatom = 1, cdft_control%natoms
99                     CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%charge(iatom)%pw, &
100                                            use_data=REALDATA3D, in_space=REALSPACE)
101                  END DO
102               END IF
103               ! Another sanity check
104               CALL get_qs_env(qs_env, natom=natom)
105               IF (natom < cdft_control%natoms) &
106                  CALL cp_abort(__LOCATION__, &
107                                "The number of constraint atoms exceeds the total number of atoms.")
108            ELSE
109               DO igroup = 1, SIZE(cdft_control%group)
110                  inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw%pw_grid%dvol
111                  CALL pw_scale(cdft_control%group(igroup)%weight%pw, inv_vol)
112               END DO
113            END IF
114            ! Build/Integrate CDFT constraints with selected population analysis method
115            IF (cdft_control%type == outer_scf_becke_constraint) THEN
116               CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
117            ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
118               CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
119            END IF
120            DO igroup = 1, SIZE(cdft_control%group)
121               CALL pw_scale(cdft_control%group(igroup)%weight%pw, cdft_control%group(igroup)%weight%pw%pw_grid%dvol)
122            END DO
123            IF (cdft_control%need_pot) cdft_control%need_pot = .FALSE.
124         CASE DEFAULT
125            CPABORT("Unknown constraint type.")
126         END SELECT
127      END IF
128
129   END SUBROUTINE qs_ks_cdft_constraint
130
131! **************************************************************************************************
132!> \brief ...
133!> \param energy ...
134!> \param dft_control ...
135!> \param just_energy ...
136!> \param para_env ...
137!> \param ks_matrix ...
138!> \param matrix_s ...
139!> \param rho ...
140!> \param mulliken_order_p ...
141! **************************************************************************************************
142   SUBROUTINE qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
143                                       ks_matrix, matrix_s, rho, mulliken_order_p)
144
145      TYPE(qs_energy_type), POINTER                      :: energy
146      TYPE(dft_control_type), POINTER                    :: dft_control
147      LOGICAL, INTENT(in)                                :: just_energy
148      TYPE(cp_para_env_type), POINTER                    :: para_env
149      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix, matrix_s
150      TYPE(qs_rho_type), POINTER                         :: rho
151      REAL(KIND=dp)                                      :: mulliken_order_p
152
153      CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_mulliken_restraint', &
154         routineP = moduleN//':'//routineN
155
156      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ksmat, rho_ao
157
158      energy%mulliken = 0.0_dp
159
160      IF (dft_control%qs_control%mulliken_restraint) THEN
161
162         ! Test no k-points
163         CPASSERT(SIZE(matrix_s, 2) == 1)
164
165         CALL qs_rho_get(rho, rho_ao=rho_ao)
166
167         IF (just_energy) THEN
168            CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
169                                    para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
170                                    order_p=mulliken_order_p)
171         ELSE
172            ksmat => ks_matrix(:, 1)
173            CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
174                                    para_env, matrix_s(1, 1)%matrix, rho_ao, energy=energy%mulliken, &
175                                    ks_matrix=ksmat, order_p=mulliken_order_p)
176         ENDIF
177
178      ENDIF
179
180   END SUBROUTINE qs_ks_mulliken_restraint
181
182! **************************************************************************************************
183!> \brief ...
184!> \param dft_control ...
185!> \param qs_env ...
186!> \param matrix_s ...
187!> \param energy ...
188!> \param calculate_forces ...
189!> \param just_energy ...
190! **************************************************************************************************
191   SUBROUTINE qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
192                                 energy, calculate_forces, just_energy)
193
194      TYPE(dft_control_type), POINTER                    :: dft_control
195      TYPE(qs_environment_type), POINTER                 :: qs_env
196      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
197      TYPE(qs_energy_type), POINTER                      :: energy
198      LOGICAL, INTENT(in)                                :: calculate_forces, just_energy
199
200      CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_s2_restraint', &
201         routineP = moduleN//':'//routineN
202
203      INTEGER                                            :: i
204      TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: fm_mo_derivs
205      TYPE(cp_fm_type), POINTER                          :: mo_coeff
206      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs, smat
207      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mo_array
208
209      NULLIFY (fm_mo_derivs, mo_array, mo_coeff, mo_derivs)
210
211      IF (dft_control%qs_control%s2_restraint) THEN
212         ! Test no k-points
213         CPASSERT(SIZE(matrix_s, 2) == 1)
214         ! adds s2_restraint energy and orbital derivatives
215         CPASSERT(dft_control%nspins == 2)
216         CPASSERT(qs_env%requires_mo_derivs)
217         ! forces are not implemented (not difficult, but ... )
218         CPASSERT(.NOT. calculate_forces)
219         CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
220
221         ALLOCATE (fm_mo_derivs(SIZE(mo_derivs, 1))) !fm->dbcsr
222         DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
223            CALL get_mo_set(mo_set=mo_array(i)%mo_set, mo_coeff=mo_coeff) !fm->dbcsr
224            CALL cp_fm_create(fm_mo_derivs(i)%matrix, mo_coeff%matrix_struct) !fm->dbcsr
225            CALL copy_dbcsr_to_fm(mo_derivs(i)%matrix, fm_mo_derivs(i)%matrix) !fm->dbcsr
226         ENDDO !fm->dbcsr
227
228         smat => matrix_s(:, 1)
229         CALL s2_restraint(mo_array, smat, fm_mo_derivs, energy%s2_restraint, &
230                           dft_control%qs_control%s2_restraint_control, just_energy)
231
232         DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
233            CALL copy_fm_to_dbcsr(fm_mo_derivs(i)%matrix, mo_derivs(i)%matrix) !fm->dbcsr
234         ENDDO !fm->dbcsr
235         DEALLOCATE (fm_mo_derivs) !fm->dbcsr
236
237      ELSE
238         energy%s2_restraint = 0.0_dp
239      ENDIF
240   END SUBROUTINE qs_ks_s2_restraint
241
242END MODULE qs_ks_apply_restraints
243