Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 : ! **************************************************************************************************
70 20106 : SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, &
71 : use_motion_section, cp_subsys, cell, cell_ref, elkind)
72 : TYPE(qs_subsys_type), INTENT(OUT) :: subsys
73 : TYPE(mp_para_env_type), POINTER :: para_env
74 : TYPE(section_vals_type), OPTIONAL, POINTER :: root_section
75 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
76 : LOGICAL, INTENT(IN) :: use_motion_section
77 : TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys
78 : TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref
79 : LOGICAL, INTENT(IN), OPTIONAL :: elkind
80 :
81 : LOGICAL :: use_ref_cell
82 6702 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
83 : TYPE(cell_type), POINTER :: my_cell, my_cell_ref
84 : TYPE(cp_subsys_type), POINTER :: my_cp_subsys
85 6702 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
86 : TYPE(section_vals_type), POINTER :: cell_section, kind_section
87 :
88 6702 : NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys)
89 :
90 : ! create cp_subsys
91 6702 : IF (PRESENT(cp_subsys)) THEN
92 530 : my_cp_subsys => cp_subsys
93 6172 : ELSE IF (PRESENT(root_section)) THEN
94 : CALL cp_subsys_create(my_cp_subsys, para_env, root_section=root_section, &
95 : force_env_section=force_env_section, &
96 : subsys_section=subsys_section, &
97 : use_motion_section=use_motion_section, &
98 6172 : elkind=elkind)
99 : ELSE
100 0 : CPABORT("qs_subsys_create: cp_subsys or root_section needed")
101 : END IF
102 :
103 : ! create cp_subsys%cell
104 : !TODO: moved to cp_subsys_create(), needs further disentanglement of cell_ref.
105 6702 : use_ref_cell = .FALSE.
106 6702 : IF (PRESENT(cell)) THEN
107 394 : my_cell => cell
108 394 : IF (PRESENT(cell_ref)) THEN
109 0 : my_cell_ref => cell_ref
110 0 : use_ref_cell = .TRUE.
111 : ELSE
112 394 : CALL cell_create(my_cell_ref)
113 394 : CALL cell_clone(my_cell, my_cell_ref, tag="CELL_REF")
114 : END IF
115 : ELSE
116 6308 : cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
117 : CALL read_cell(my_cell, my_cell_ref, use_ref_cell=use_ref_cell, &
118 6308 : cell_section=cell_section, para_env=para_env)
119 : END IF
120 6702 : CALL cp_subsys_set(my_cp_subsys, cell=my_cell)
121 6702 : CALL write_cell(my_cell, subsys_section)
122 6702 : CALL write_cell(my_cell_ref, subsys_section)
123 :
124 : ! setup qs_kinds
125 6702 : CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set)
126 6702 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
127 : CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
128 6702 : para_env, force_env_section)
129 :
130 : CALL num_ao_el_per_molecule(my_cp_subsys%molecule_kinds%els, &
131 6702 : qs_kind_set)
132 :
133 : CALL qs_subsys_set(subsys, &
134 : cp_subsys=my_cp_subsys, &
135 : cell_ref=my_cell_ref, &
136 : use_ref_cell=use_ref_cell, &
137 6702 : qs_kind_set=qs_kind_set)
138 :
139 6702 : IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell)
140 6702 : IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref)
141 6702 : IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys)
142 6702 : END SUBROUTINE qs_subsys_create
143 :
144 : ! **************************************************************************************************
145 : !> \brief Read a molecule kind data set from the input file.
146 : !> \param molecule_kind_set ...
147 : !> \param qs_kind_set ...
148 : !> \date 22.11.2004
149 : !> \par History
150 : !> Rustam Z. Khaliullin 10.2014 - charges and electrons of molecules
151 : !> are now in agreement with atomic guess
152 : !> \author MI
153 : !> \version 1.0
154 : ! **************************************************************************************************
155 6702 : SUBROUTINE num_ao_el_per_molecule(molecule_kind_set, qs_kind_set)
156 :
157 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
158 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
159 :
160 : INTEGER :: arbitrary_spin, iatom, ikind, imol, &
161 : n_ao, natom, nmol_kind, nsgf, nspins, &
162 : z_molecule
163 : INTEGER, DIMENSION(0:lmat, 10) :: ne_core, ne_elem, ne_explicit
164 : INTEGER, DIMENSION(2) :: n_occ_alpha_and_beta
165 : REAL(KIND=dp) :: charge_molecule, zeff, zeff_correction
166 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta
167 : TYPE(all_potential_type), POINTER :: all_potential
168 : TYPE(atomic_kind_type), POINTER :: atomic_kind
169 : TYPE(gth_potential_type), POINTER :: gth_potential
170 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
171 : TYPE(molecule_kind_type), POINTER :: molecule_kind
172 : TYPE(sgp_potential_type), POINTER :: sgp_potential
173 :
174 6702 : IF (ASSOCIATED(molecule_kind_set)) THEN
175 :
176 6702 : nspins = 2
177 6702 : nmol_kind = SIZE(molecule_kind_set, 1)
178 : natom = 0
179 :
180 : ! *** Initialize the molecule kind data structure ***
181 6702 : ARBITRARY_SPIN = 1
182 38623 : DO imol = 1, nmol_kind
183 :
184 31921 : molecule_kind => molecule_kind_set(imol)
185 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
186 31921 : natom=natom)
187 : !nelectron = 0
188 31921 : n_ao = 0
189 95763 : n_occ_alpha_and_beta(1:nspins) = 0
190 31921 : z_molecule = 0
191 :
192 65666 : DO iatom = 1, natom
193 :
194 33745 : atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind
195 33745 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
196 : CALL get_qs_kind(qs_kind_set(ikind), &
197 : basis_set=orb_basis_set, &
198 : all_potential=all_potential, &
199 : gth_potential=gth_potential, &
200 33745 : sgp_potential=sgp_potential)
201 :
202 : ! Obtain the electronic state of the atom
203 : ! The same state is used to calculate the ATOMIC GUESS
204 : ! It is great that we are consistent with ATOMIC_GUESS
205 : CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
206 : qs_kind=qs_kind_set(ikind), &
207 : ncalc=ne_explicit, &
208 : ncore=ne_core, &
209 : nelem=ne_elem, &
210 33745 : edelta=edelta)
211 :
212 : ! If &BS section is used ATOMIC_GUESS is calculated twice
213 : ! for two separate wfns with their own alpha-beta combinations
214 : ! This is done to break the spin symmetry of the initial wfn
215 : ! For now, only alpha part of &BS is used to count electrons on
216 : ! molecules
217 : ! Get the number of explicit electrons (i.e. with orbitals)
218 : ! For now, only the total number of electrons can be obtained
219 : ! from init_atom_electronic_state
220 : n_occ_alpha_and_beta(ARBITRARY_SPIN) = &
221 : n_occ_alpha_and_beta(ARBITRARY_SPIN) + SUM(ne_explicit) + &
222 4758045 : SUM(NINT(2*edelta(:, :, ARBITRARY_SPIN)))
223 : ! We need a way to specify the number of alpha and beta electrons
224 : ! on each molecule (i.e. multiplicity is not enough)
225 : !n_occ(ispin) = n_occ(ispin) + SUM(ne_explicit) + SUM(NINT(2*edelta(:, :, ispin)))
226 :
227 33745 : IF (ASSOCIATED(all_potential)) THEN
228 : CALL get_potential(potential=all_potential, zeff=zeff, &
229 16724 : zeff_correction=zeff_correction)
230 17021 : ELSE IF (ASSOCIATED(gth_potential)) THEN
231 : CALL get_potential(potential=gth_potential, zeff=zeff, &
232 16713 : zeff_correction=zeff_correction)
233 308 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
234 : CALL get_potential(potential=sgp_potential, zeff=zeff, &
235 34 : zeff_correction=zeff_correction)
236 : ELSE
237 274 : zeff = 0.0_dp
238 274 : zeff_correction = 0.0_dp
239 : END IF
240 33745 : z_molecule = z_molecule + NINT(zeff - zeff_correction)
241 :
242 : ! this one does not work because nelem is not adjusted in the symmetry breaking code
243 : !CALL get_atomic_kind(atomic_kind,z=z)
244 : !z_molecule=z_molecule+z
245 :
246 33745 : IF (ASSOCIATED(orb_basis_set)) THEN
247 30765 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf)
248 : ELSE
249 2980 : nsgf = 0
250 : END IF
251 99411 : n_ao = n_ao + nsgf
252 :
253 : END DO ! iatom
254 :
255 : ! At this point we have the number of electrons (alpha+beta) on the molecule
256 : ! as they are seen by the ATOMIC GUESS routines
257 31921 : charge_molecule = REAL(z_molecule - n_occ_alpha_and_beta(ARBITRARY_SPIN), dp)
258 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
259 : nelectron=n_occ_alpha_and_beta(ARBITRARY_SPIN), &
260 : charge=charge_molecule, &
261 70544 : nsgf=n_ao)
262 :
263 : END DO ! imol
264 : END IF
265 :
266 6702 : END SUBROUTINE num_ao_el_per_molecule
267 :
268 : END MODULE qs_subsys_methods
|
