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