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