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_kind_types, ONLY: create_qs_kind_set,&
41 : get_qs_kind,&
42 : init_atom_electronic_state,&
43 : qs_kind_type
44 : USE qs_subsys_types, ONLY: qs_subsys_set,&
45 : qs_subsys_type
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_subsys_methods'
52 :
53 : PUBLIC :: qs_subsys_create
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Creates a qs_subsys. Optionally an existsing cp_subsys is used.
59 : !> \param subsys ...
60 : !> \param para_env ...
61 : !> \param root_section ...
62 : !> \param force_env_section ...
63 : !> \param subsys_section ...
64 : !> \param use_motion_section ...
65 : !> \param cp_subsys ...
66 : !> \param cell ...
67 : !> \param cell_ref ...
68 : !> \param elkind ...
69 : !> \param silent ...
70 : ! **************************************************************************************************
71 22260 : SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, &
72 : use_motion_section, cp_subsys, cell, cell_ref, elkind, silent)
73 : TYPE(qs_subsys_type), INTENT(OUT) :: subsys
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : TYPE(section_vals_type), OPTIONAL, POINTER :: root_section
76 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
77 : LOGICAL, INTENT(IN) :: use_motion_section
78 : TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys
79 : TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref
80 : LOGICAL, INTENT(IN), OPTIONAL :: elkind, silent
81 :
82 : LOGICAL :: be_silent, use_ref_cell
83 7420 : 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 7420 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
87 : TYPE(section_vals_type), POINTER :: cell_section, kind_section
88 :
89 7420 : NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys)
90 :
91 7420 : be_silent = .FALSE.
92 7420 : IF (PRESENT(silent)) be_silent = silent
93 : ! create cp_subsys
94 7420 : IF (PRESENT(cp_subsys)) THEN
95 544 : my_cp_subsys => cp_subsys
96 6876 : 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 6876 : elkind=elkind)
102 : ELSE
103 0 : 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 7420 : use_ref_cell = .FALSE.
109 7420 : IF (PRESENT(cell)) THEN
110 394 : my_cell => cell
111 394 : IF (PRESENT(cell_ref)) THEN
112 0 : my_cell_ref => cell_ref
113 0 : use_ref_cell = .TRUE.
114 : ELSE
115 394 : CALL cell_create(my_cell_ref)
116 394 : CALL cell_clone(my_cell, my_cell_ref, tag="CELL_REF")
117 : END IF
118 : ELSE
119 7026 : 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 7026 : cell_section=cell_section, para_env=para_env)
122 : END IF
123 7420 : CALL cp_subsys_set(my_cp_subsys, cell=my_cell)
124 7420 : CALL write_cell(my_cell, subsys_section)
125 7420 : CALL write_cell(my_cell_ref, subsys_section)
126 :
127 : ! setup qs_kinds
128 7420 : CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set)
129 7420 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
130 : CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
131 7420 : para_env, force_env_section, be_silent)
132 :
133 : CALL num_ao_el_per_molecule(my_cp_subsys%molecule_kinds%els, &
134 7420 : qs_kind_set)
135 :
136 : CALL qs_subsys_set(subsys, &
137 : cp_subsys=my_cp_subsys, &
138 : cell_ref=my_cell_ref, &
139 : use_ref_cell=use_ref_cell, &
140 7420 : qs_kind_set=qs_kind_set)
141 :
142 7420 : IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell)
143 7420 : IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref)
144 7420 : IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys)
145 :
146 7420 : 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 7420 : 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 : INTEGER :: arbitrary_spin, iatom, ikind, imol, &
165 : n_ao, natom, nmol_kind, nsgf, nspins, &
166 : z_molecule
167 : INTEGER, DIMENSION(0:lmat, 10) :: ne_core, ne_elem, ne_explicit
168 : INTEGER, DIMENSION(2) :: n_occ_alpha_and_beta
169 : REAL(KIND=dp) :: charge_molecule, zeff, zeff_correction
170 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta
171 : TYPE(all_potential_type), POINTER :: all_potential
172 : TYPE(atomic_kind_type), POINTER :: atomic_kind
173 : TYPE(gth_potential_type), POINTER :: gth_potential
174 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
175 : TYPE(molecule_kind_type), POINTER :: molecule_kind
176 : TYPE(sgp_potential_type), POINTER :: sgp_potential
177 :
178 7420 : IF (ASSOCIATED(molecule_kind_set)) THEN
179 :
180 7420 : nspins = 2
181 7420 : nmol_kind = SIZE(molecule_kind_set, 1)
182 : natom = 0
183 :
184 : ! *** Initialize the molecule kind data structure ***
185 7420 : ARBITRARY_SPIN = 1
186 42389 : DO imol = 1, nmol_kind
187 :
188 34969 : molecule_kind => molecule_kind_set(imol)
189 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
190 34969 : natom=natom)
191 : !nelectron = 0
192 34969 : n_ao = 0
193 104907 : n_occ_alpha_and_beta(1:nspins) = 0
194 34969 : z_molecule = 0
195 :
196 71762 : DO iatom = 1, natom
197 :
198 36793 : atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind
199 36793 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
200 : CALL get_qs_kind(qs_kind_set(ikind), &
201 : basis_set=orb_basis_set, &
202 : all_potential=all_potential, &
203 : gth_potential=gth_potential, &
204 36793 : sgp_potential=sgp_potential)
205 :
206 : ! Obtain the electronic state of the atom
207 : ! The same state is used to calculate the ATOMIC GUESS
208 : ! It is great that we are consistent with ATOMIC_GUESS
209 : CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
210 : qs_kind=qs_kind_set(ikind), &
211 : ncalc=ne_explicit, &
212 : ncore=ne_core, &
213 : nelem=ne_elem, &
214 36793 : edelta=edelta)
215 :
216 : ! If &BS section is used ATOMIC_GUESS is calculated twice
217 : ! for two separate wfns with their own alpha-beta combinations
218 : ! This is done to break the spin symmetry of the initial wfn
219 : ! For now, only alpha part of &BS is used to count electrons on
220 : ! molecules
221 : ! Get the number of explicit electrons (i.e. with orbitals)
222 : ! For now, only the total number of electrons can be obtained
223 : ! from init_atom_electronic_state
224 : n_occ_alpha_and_beta(ARBITRARY_SPIN) = &
225 : n_occ_alpha_and_beta(ARBITRARY_SPIN) + SUM(ne_explicit) + &
226 5187813 : SUM(NINT(2*edelta(:, :, ARBITRARY_SPIN)))
227 : ! We need a way to specify the number of alpha and beta electrons
228 : ! on each molecule (i.e. multiplicity is not enough)
229 : !n_occ(ispin) = n_occ(ispin) + SUM(ne_explicit) + SUM(NINT(2*edelta(:, :, ispin)))
230 :
231 36793 : IF (ASSOCIATED(all_potential)) THEN
232 : CALL get_potential(potential=all_potential, zeff=zeff, &
233 19630 : zeff_correction=zeff_correction)
234 17163 : ELSE IF (ASSOCIATED(gth_potential)) THEN
235 : CALL get_potential(potential=gth_potential, zeff=zeff, &
236 16827 : zeff_correction=zeff_correction)
237 336 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
238 : CALL get_potential(potential=sgp_potential, zeff=zeff, &
239 50 : zeff_correction=zeff_correction)
240 : ELSE
241 286 : zeff = 0.0_dp
242 286 : zeff_correction = 0.0_dp
243 : END IF
244 36793 : z_molecule = z_molecule + NINT(zeff - zeff_correction)
245 :
246 : ! this one does not work because nelem is not adjusted in the symmetry breaking code
247 : !CALL get_atomic_kind(atomic_kind,z=z)
248 : !z_molecule=z_molecule+z
249 :
250 36793 : IF (ASSOCIATED(orb_basis_set)) THEN
251 30955 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf)
252 : ELSE
253 5838 : nsgf = 0
254 : END IF
255 108555 : n_ao = n_ao + nsgf
256 :
257 : END DO ! iatom
258 :
259 : ! At this point we have the number of electrons (alpha+beta) on the molecule
260 : ! as they are seen by the ATOMIC GUESS routines
261 34969 : charge_molecule = REAL(z_molecule - n_occ_alpha_and_beta(ARBITRARY_SPIN), dp)
262 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
263 : nelectron=n_occ_alpha_and_beta(ARBITRARY_SPIN), &
264 : charge=charge_molecule, &
265 77358 : nsgf=n_ao)
266 :
267 : END DO ! imol
268 : END IF
269 :
270 7420 : END SUBROUTINE num_ao_el_per_molecule
271 :
272 : END MODULE qs_subsys_methods
|
