1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines that work on qs_subsys_type
8!> \author Ole Schuett
9! **************************************************************************************************
10MODULE qs_subsys_methods
11   USE atom_types,                      ONLY: lmat
12   USE atomic_kind_types,               ONLY: atomic_kind_type,&
13                                              get_atomic_kind
14   USE basis_set_types,                 ONLY: get_gto_basis_set,&
15                                              gto_basis_set_type
16   USE cell_methods,                    ONLY: read_cell,&
17                                              write_cell
18   USE cell_types,                      ONLY: cell_clone,&
19                                              cell_create,&
20                                              cell_release,&
21                                              cell_type
22   USE cp_para_types,                   ONLY: cp_para_env_type
23   USE cp_subsys_methods,               ONLY: cp_subsys_create
24   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
25                                              cp_subsys_release,&
26                                              cp_subsys_set,&
27                                              cp_subsys_type
28   USE external_potential_types,        ONLY: all_potential_type,&
29                                              get_potential,&
30                                              gth_potential_type,&
31                                              sgp_potential_type
32   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
33                                              section_vals_type
34   USE kinds,                           ONLY: dp
35   USE molecule_kind_types,             ONLY: get_molecule_kind,&
36                                              molecule_kind_type,&
37                                              set_molecule_kind
38   USE qs_kind_types,                   ONLY: create_qs_kind_set,&
39                                              get_qs_kind,&
40                                              init_atom_electronic_state,&
41                                              qs_kind_type
42   USE qs_subsys_types,                 ONLY: qs_subsys_set,&
43                                              qs_subsys_type
44#include "./base/base_uses.f90"
45
46   IMPLICIT NONE
47   PRIVATE
48
49   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_subsys_methods'
50
51   PUBLIC :: qs_subsys_create
52
53CONTAINS
54
55! **************************************************************************************************
56!> \brief Creates a qs_subsys. Optionally an existsing cp_subsys is used.
57!> \param subsys ...
58!> \param para_env ...
59!> \param root_section ...
60!> \param force_env_section ...
61!> \param subsys_section ...
62!> \param use_motion_section ...
63!> \param cp_subsys ...
64!> \param cell ...
65!> \param cell_ref ...
66!> \param elkind ...
67! **************************************************************************************************
68   SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, &
69                               use_motion_section, cp_subsys, cell, cell_ref, elkind)
70      TYPE(qs_subsys_type), POINTER                      :: subsys
71      TYPE(cp_para_env_type), POINTER                    :: para_env
72      TYPE(section_vals_type), OPTIONAL, POINTER         :: root_section
73      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
74      LOGICAL, INTENT(IN)                                :: use_motion_section
75      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: cp_subsys
76      TYPE(cell_type), OPTIONAL, POINTER                 :: cell, cell_ref
77      LOGICAL, INTENT(IN), OPTIONAL                      :: elkind
78
79      CHARACTER(len=*), PARAMETER :: routineN = 'qs_subsys_create', &
80         routineP = moduleN//':'//routineN
81
82      LOGICAL                                            :: use_ref_cell
83      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
84      TYPE(cell_type), POINTER                           :: my_cell, my_cell_ref
85      TYPE(cp_subsys_type), POINTER                      :: my_cp_subsys
86      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
87      TYPE(section_vals_type), POINTER                   :: cell_section, kind_section
88
89      NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys)
90
91      IF (ASSOCIATED(subsys)) CPABORT("qs_subsys_create: subsys already associated")
92
93      ! create cp_subsys
94      IF (PRESENT(cp_subsys)) THEN
95         my_cp_subsys => cp_subsys
96      ELSE IF (PRESENT(root_section)) THEN
97         CALL cp_subsys_create(my_cp_subsys, para_env, root_section=root_section, &
98                               force_env_section=force_env_section, &
99                               subsys_section=subsys_section, &
100                               use_motion_section=use_motion_section, &
101                               elkind=elkind)
102      ELSE
103         CPABORT("qs_subsys_create: cp_subsys or root_section needed")
104      END IF
105
106      ! create cp_subsys%cell
107      !TODO: moved to cp_subsys_create(), needs further disentanglement of cell_ref.
108      IF (PRESENT(cell)) THEN
109         my_cell => cell
110         IF (PRESENT(cell_ref)) THEN
111            my_cell_ref => cell_ref
112            use_ref_cell = .TRUE.
113         ELSE
114            CALL cell_create(my_cell_ref)
115            CALL cell_clone(my_cell, my_cell_ref)
116            use_ref_cell = .FALSE.
117         END IF
118      ELSE
119         cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
120         CALL read_cell(my_cell, my_cell_ref, use_ref_cell=use_ref_cell, &
121                        cell_section=cell_section, para_env=para_env)
122      END IF
123      CALL cp_subsys_set(my_cp_subsys, cell=my_cell)
124      CALL write_cell(my_cell, subsys_section, cell_ref=my_cell_ref)
125
126      ! setup qs_kinds
127      CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set)
128      kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
129      CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
130                              para_env, force_env_section)
131
132      CALL num_ao_el_per_molecule(my_cp_subsys%molecule_kinds%els, &
133                                  qs_kind_set)
134
135      ! finally create qs_subsys
136      ALLOCATE (subsys)
137      CALL qs_subsys_set(subsys, &
138                         cp_subsys=my_cp_subsys, &
139                         cell_ref=my_cell_ref, &
140                         use_ref_cell=use_ref_cell, &
141                         qs_kind_set=qs_kind_set)
142
143      IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell)
144      IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref)
145      IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys)
146   END SUBROUTINE qs_subsys_create
147
148! **************************************************************************************************
149!> \brief   Read a molecule kind data set from the input file.
150!> \param molecule_kind_set ...
151!> \param qs_kind_set ...
152!> \date    22.11.2004
153!> \par History
154!>      Rustam Z. Khaliullin 10.2014 - charges and electrons of molecules
155!>                                     are now in agreement with atomic guess
156!> \author  MI
157!> \version 1.0
158! **************************************************************************************************
159   SUBROUTINE num_ao_el_per_molecule(molecule_kind_set, qs_kind_set)
160
161      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
162      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
163
164      CHARACTER(len=*), PARAMETER :: routineN = 'num_ao_el_per_molecule', &
165         routineP = moduleN//':'//routineN
166
167      INTEGER                                            :: arbitrary_spin, iatom, ikind, imol, &
168                                                            n_ao, natom, nmol_kind, nsgf, nspins, &
169                                                            z_molecule
170      INTEGER, DIMENSION(0:lmat, 10)                     :: ne_core, ne_elem, ne_explicit
171      INTEGER, DIMENSION(2)                              :: n_occ_alpha_and_beta
172      REAL(KIND=dp)                                      :: charge_molecule, zeff, zeff_correction
173      REAL(KIND=dp), DIMENSION(0:lmat, 10, 2)            :: edelta
174      TYPE(all_potential_type), POINTER                  :: all_potential
175      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
176      TYPE(gth_potential_type), POINTER                  :: gth_potential
177      TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
178      TYPE(molecule_kind_type), POINTER                  :: molecule_kind
179      TYPE(sgp_potential_type), POINTER                  :: sgp_potential
180
181      IF (ASSOCIATED(molecule_kind_set)) THEN
182
183         nspins = 2
184         nmol_kind = SIZE(molecule_kind_set, 1)
185         natom = 0
186
187         !   *** Initialize the molecule kind data structure ***
188         ARBITRARY_SPIN = 1
189         DO imol = 1, nmol_kind
190
191            molecule_kind => molecule_kind_set(imol)
192            CALL get_molecule_kind(molecule_kind=molecule_kind, &
193                                   natom=natom)
194            !nelectron = 0
195            n_ao = 0
196            n_occ_alpha_and_beta(1:nspins) = 0
197            z_molecule = 0
198
199            DO iatom = 1, natom
200
201               atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind
202               CALL get_atomic_kind(atomic_kind, kind_number=ikind)
203               CALL get_qs_kind(qs_kind_set(ikind), &
204                                basis_set=orb_basis_set, &
205                                all_potential=all_potential, &
206                                gth_potential=gth_potential, &
207                                sgp_potential=sgp_potential)
208
209               ! Obtain the electronic state of the atom
210               ! The same state is used to calculate the ATOMIC GUESS
211               ! It is great that we are consistent with ATOMIC_GUESS
212               CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
213                                               qs_kind=qs_kind_set(ikind), &
214                                               ncalc=ne_explicit, &
215                                               ncore=ne_core, &
216                                               nelem=ne_elem, &
217                                               edelta=edelta)
218
219               ! If &BS section is used ATOMIC_GUESS is calculated twice
220               ! for two separate wfns with their own alpha-beta combinations
221               ! This is done to break the spin symmetry of the initial wfn
222               ! For now, only alpha part of &BS is used to count electrons on
223               ! molecules
224               ! Get the number of explicit electrons (i.e. with orbitals)
225               ! For now, only the total number of electrons can be obtained
226               ! from init_atom_electronic_state
227               n_occ_alpha_and_beta(ARBITRARY_SPIN) = &
228                  n_occ_alpha_and_beta(ARBITRARY_SPIN) + SUM(ne_explicit) + &
229                  SUM(NINT(2*edelta(:, :, ARBITRARY_SPIN)))
230               ! We need a way to specify the number of alpha and beta electrons
231               ! on each molecule (i.e. multiplicity is not enough)
232               !n_occ(ispin) = n_occ(ispin) + SUM(ne_explicit) + SUM(NINT(2*edelta(:, :, ispin)))
233
234               IF (ASSOCIATED(all_potential)) THEN
235                  CALL get_potential(potential=all_potential, zeff=zeff, &
236                                     zeff_correction=zeff_correction)
237               ELSE IF (ASSOCIATED(gth_potential)) THEN
238                  CALL get_potential(potential=gth_potential, zeff=zeff, &
239                                     zeff_correction=zeff_correction)
240               ELSE IF (ASSOCIATED(sgp_potential)) THEN
241                  CALL get_potential(potential=sgp_potential, zeff=zeff, &
242                                     zeff_correction=zeff_correction)
243               ELSE
244                  zeff = 0.0_dp
245                  zeff_correction = 0.0_dp
246               END IF
247               z_molecule = z_molecule + NINT(zeff - zeff_correction)
248
249               ! this one does not work because nelem is not adjusted in the symmetry breaking code
250               !CALL get_atomic_kind(atomic_kind,z=z)
251               !z_molecule=z_molecule+z
252
253               IF (ASSOCIATED(orb_basis_set)) THEN
254                  CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf)
255               ELSE
256                  nsgf = 0
257               ENDIF
258               n_ao = n_ao + nsgf
259
260            END DO ! iatom
261
262            ! At this point we have the number of electrons (alpha+beta) on the molecule
263            !  as they are seen by the ATOMIC GUESS routines
264            charge_molecule = REAL(z_molecule - n_occ_alpha_and_beta(ARBITRARY_SPIN), dp)
265            CALL set_molecule_kind(molecule_kind=molecule_kind, &
266                                   nelectron=n_occ_alpha_and_beta(ARBITRARY_SPIN), &
267                                   charge=charge_molecule, &
268                                   nsgf=n_ao)
269
270         END DO ! imol
271      END IF
272
273   END SUBROUTINE num_ao_el_per_molecule
274
275END MODULE qs_subsys_methods
276