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 Set disperson types for DFT calculations
10 : !> \author JGH (04.2014)
11 : ! **************************************************************************************************
12 : MODULE qs_gcp_utils
13 :
14 : USE basis_set_types, ONLY: get_gto_basis_set,&
15 : gto_basis_set_type
16 : USE cp_parser_methods, ONLY: parser_get_next_line,&
17 : parser_get_object
18 : USE cp_parser_types, ONLY: cp_parser_type,&
19 : parser_create,&
20 : parser_release
21 : USE input_section_types, ONLY: section_vals_get,&
22 : section_vals_get_subs_vals,&
23 : section_vals_type,&
24 : section_vals_val_get
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE mathconstants, ONLY: pi
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE periodic_table, ONLY: get_ptable_info
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_gcp_types, ONLY: qs_gcp_type
33 : USE qs_kind_types, ONLY: get_qs_kind,&
34 : qs_kind_type
35 : USE sto_ng, ONLY: get_sto_ng
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gcp_utils'
43 :
44 : INTEGER, DIMENSION(106) :: nshell = [ &
45 : 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 1-30
46 : 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, & ! 31-60
47 : 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, & ! 61-90
48 : 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7] ! 91-106
49 :
50 : INTEGER, DIMENSION(106) :: nll = [ &
51 : 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, & ! 1-30
52 : 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, & ! 31-60
53 : 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, & ! 61-90
54 : 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 0, 0] ! 91-106
55 :
56 : !
57 : ! Slater exponents for valence states
58 : ! Aleksander Herman
59 : ! Empirically adjusted and consistent set of EHT valence orbital parameters for all elements of the periodic table
60 : ! Modelling and Simulation in Materials Science and Engineering, 12, 21-32 (2004)
61 : !
62 : ! Hydrogen uses 1.2000, not the original 1.0000
63 : !
64 : REAL(KIND=dp), DIMENSION(4, 106) :: sexp = RESHAPE([ &
65 : 1.2000, 0.0000, 0.0000, 0.0000, 1.6469, 0.0000, 0.0000, 0.0000, 0.6534, 0.5305, 0.0000, 0.0000, & ! 1 2 3
66 : 1.0365, 0.8994, 0.0000, 0.0000, 1.3990, 1.2685, 0.0000, 0.0000, 1.7210, 1.6105, 0.0000, 0.0000, & ! 4 5 6
67 : 2.0348, 1.9398, 0.0000, 0.0000, 2.2399, 2.0477, 0.0000, 0.0000, 2.5644, 2.4022, 0.0000, 0.0000, & ! 7 8 9
68 : 2.8812, 2.7421, 0.0000, 0.0000, 0.8675, 0.6148, 0.0000, 0.0000, 1.1935, 0.8809, 0.0000, 0.0000, & ! 10 11 12
69 : 1.5143, 1.1660, 0.0000, 0.0000, 1.7580, 1.4337, 0.0000, 0.0000, 1.9860, 1.6755, 0.0000, 0.0000, & ! 13 14 15
70 : 2.1362, 1.7721, 0.0000, 0.0000, 2.3617, 2.0176, 0.0000, 0.0000, 2.5796, 2.2501, 0.0000, 0.0000, & ! 16 17 18
71 : 0.9362, 0.6914, 0.0000, 0.0000, 1.2112, 0.9329, 0.0000, 0.0000, 1.2870, 0.9828, 2.4341, 0.0000, & ! 19 20 21
72 : 1.3416, 1.0104, 2.6439, 0.0000, 1.3570, 0.9947, 2.7809, 0.0000, 1.3804, 0.9784, 2.9775, 0.0000, & ! 22 23 24
73 : 1.4761, 1.0641, 3.2208, 0.0000, 1.5465, 1.1114, 3.4537, 0.0000, 1.5650, 1.1001, 3.6023, 0.0000, & ! 25 26 27
74 : 1.5532, 1.0594, 3.7017, 0.0000, 1.5791, 1.0527, 3.8962, 0.0000, 1.7778, 1.2448, 0.0000, 0.0000, & ! 28 29 30
75 : 2.0675, 1.5073, 0.0000, 0.0000, 2.2702, 1.7680, 0.0000, 0.0000, 2.4546, 1.9819, 0.0000, 0.0000, & ! 31 32 33
76 : 2.5680, 2.0548, 0.0000, 0.0000, 2.7523, 2.2652, 0.0000, 0.0000, 2.9299, 2.4617, 0.0000, 0.0000, & ! 34 35 36
77 : 1.0963, 0.7990, 0.0000, 0.0000, 1.3664, 1.0415, 0.0000, 0.0000, 1.4613, 1.1100, 2.1576, 0.0000, & ! 37 38 39
78 : 1.5393, 1.1647, 2.3831, 0.0000, 1.5926, 1.1738, 2.6256, 0.0000, 1.6579, 1.2186, 2.8241, 0.0000, & ! 40 41 42
79 : 1.6930, 1.2490, 2.9340, 0.0000, 1.7347, 1.2514, 3.1524, 0.0000, 1.7671, 1.2623, 3.3113, 0.0000, & ! 43 44 45
80 : 1.6261, 1.1221, 3.0858, 0.0000, 1.8184, 1.2719, 3.6171, 0.0000, 1.9900, 1.4596, 0.0000, 0.0000, & ! 46 47 48
81 : 2.4649, 1.6848, 0.0000, 0.0000, 2.4041, 1.9128, 0.0000, 0.0000, 2.5492, 2.0781, 0.0000, 0.0000, & ! 49 50 51
82 : 2.6576, 2.1718, 0.0000, 0.0000, 2.8080, 2.3390, 0.0000, 0.0000, 2.9595, 2.5074, 0.0000, 0.0000, & ! 52 53 54
83 : 1.1993, 0.8918, 0.0000, 0.0000, 1.4519, 1.1397, 0.0000, 0.0000, 1.5331, 1.1979, 2.2743, 4.4161, & ! 55 56 57
84 : 1.5379, 1.1930, 2.2912, 4.9478, 1.5162, 1.1834, 2.0558, 4.8982, 1.5322, 1.1923, 2.0718, 5.0744, & ! 58 59 60
85 : 1.5486, 1.2018, 2.0863, 5.2466, 1.5653, 1.2118, 2.0999, 5.4145, 1.5762, 1.2152, 2.0980, 5.5679, & ! 61 62 63
86 : 1.6703, 1.2874, 2.4862, 5.9888, 1.6186, 1.2460, 2.1383, 5.9040, 1.6358, 1.2570, 2.1472, 6.0598, & ! 64 65 66
87 : 1.6536, 1.2687, 2.1566, 6.2155, 1.6723, 1.2813, 2.1668, 6.3703, 1.6898, 1.2928, 2.1731, 6.5208, & ! 67 68 69
88 : 1.7063, 1.3030, 2.1754, 6.6686, 1.6647, 1.2167, 2.3795, 0.0000, 1.8411, 1.3822, 2.7702, 0.0000, & ! 70 71 72
89 : 1.9554, 1.4857, 3.0193, 0.0000, 2.0190, 1.5296, 3.1936, 0.0000, 2.0447, 1.5276, 3.3237, 0.0000, & ! 73 74 75
90 : 2.1361, 1.6102, 3.5241, 0.0000, 2.2167, 1.6814, 3.7077, 0.0000, 2.2646, 1.6759, 3.8996, 0.0000, & ! 76 77 78
91 : 2.3185, 1.7126, 4.0525, 0.0000, 2.4306, 1.8672, 0.0000, 0.0000, 2.5779, 1.9899, 0.0000, 0.0000, & ! 79 80 81
92 : 2.7241, 2.1837, 0.0000, 0.0000, 2.7869, 2.2146, 0.0000, 0.0000, 2.9312, 2.3830, 0.0000, 0.0000, & ! 82 83 84
93 : 3.1160, 2.6200, 0.0000, 0.0000, 3.2053, 2.6866, 0.0000, 0.0000, 1.4160, 1.0598, 0.0000, 0.0000, & ! 85 86 87
94 : 1.6336, 1.3011, 0.0000, 0.0000, 1.6540, 1.2890, 2.3740, 3.7960, 1.8381, 1.4726, 2.6584, 4.3613, & ! 88 89 90
95 : 1.7770, 1.4120, 2.5710, 4.5540, 1.8246, 1.4588, 2.6496, 4.7702, 1.8451, 1.4739, 2.6940, 4.9412, & ! 91 92 93
96 : 1.7983, 1.4366, 2.5123, 4.9882, 1.8011, 1.4317, 2.5170, 5.1301, 1.8408, 1.4418, 2.7349, 5.3476, & ! 94 95 96
97 : 1.8464, 1.4697, 2.5922, 5.4596, 1.8647, 1.4838, 2.6205, 5.6140, 1.8890, 1.5050, 2.6590, 5.7740, & ! 97 98 99
98 : 1.9070, 1.5190, 2.6850, 5.9220, 1.9240, 1.5320, 2.7090, 6.0690, 1.9400, 1.5440, 2.7300, 6.2130, & ! 100 101 102
99 : 2.1300, 1.7200, 2.9900, 0.0000, 1.9200, 1.4500, 2.9700, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, & ! 103 104 105
100 : 0.0000, 0.0000, 0.0000, 0.0000], & ! 106
101 : [4, 106])
102 :
103 : PUBLIC :: qs_gcp_env_set, qs_gcp_init
104 :
105 : ! **************************************************************************************************
106 : CONTAINS
107 : ! **************************************************************************************************
108 : !> \brief ...
109 : !> \param gcp_env ...
110 : !> \param xc_section ...
111 : ! **************************************************************************************************
112 10560 : SUBROUTINE qs_gcp_env_set(gcp_env, xc_section)
113 : TYPE(qs_gcp_type), POINTER :: gcp_env
114 : TYPE(section_vals_type), POINTER :: xc_section
115 :
116 : CHARACTER(LEN=default_string_length), &
117 5280 : DIMENSION(:), POINTER :: tmpstringlist
118 : INTEGER :: i_rep, n_rep
119 : LOGICAL :: explicit
120 5280 : REAL(dp), POINTER :: params(:)
121 : TYPE(section_vals_type), POINTER :: gcp_section
122 :
123 0 : CPASSERT(ASSOCIATED(gcp_env))
124 :
125 5280 : gcp_section => section_vals_get_subs_vals(xc_section, "GCP_POTENTIAL")
126 5280 : CALL section_vals_get(gcp_section, explicit=explicit)
127 5280 : IF (explicit) THEN
128 4 : CALL section_vals_val_get(gcp_section, "VERBOSE", l_val=gcp_env%verbose)
129 4 : gcp_env%do_gcp = .TRUE.
130 : CALL section_vals_val_get(gcp_section, "PARAMETER_FILE_NAME", &
131 4 : c_val=gcp_env%parameter_file_name)
132 4 : CALL section_vals_val_get(gcp_section, "GLOBAL_PARAMETERS", r_vals=params)
133 4 : gcp_env%sigma = params(1)
134 4 : gcp_env%alpha = params(2)
135 4 : gcp_env%beta = params(3)
136 4 : gcp_env%eta = params(4)
137 : ! eamiss definitions
138 4 : CALL section_vals_val_get(gcp_section, "DELTA_ENERGY", n_rep_val=n_rep)
139 4 : IF (n_rep > 0) THEN
140 12 : ALLOCATE (gcp_env%kind_type(n_rep))
141 12 : ALLOCATE (gcp_env%ea(n_rep))
142 16 : DO i_rep = 1, n_rep
143 : CALL section_vals_val_get(gcp_section, "DELTA_ENERGY", i_rep_val=i_rep, &
144 12 : c_vals=tmpstringlist)
145 12 : READ (tmpstringlist(1), *) gcp_env%kind_type(i_rep)
146 16 : READ (tmpstringlist(2), *) gcp_env%ea(i_rep)
147 : END DO
148 : END IF
149 : ELSE
150 5276 : gcp_env%do_gcp = .FALSE.
151 : END IF
152 :
153 5280 : END SUBROUTINE qs_gcp_env_set
154 :
155 : ! **************************************************************************************************
156 : !> \brief ...
157 : !> \param qs_env ...
158 : !> \param gcp_env ...
159 : ! **************************************************************************************************
160 5280 : SUBROUTINE qs_gcp_init(qs_env, gcp_env)
161 : TYPE(qs_environment_type), POINTER :: qs_env
162 : TYPE(qs_gcp_type), POINTER :: gcp_env
163 :
164 : REAL(KIND=dp), PARAMETER :: epsc = 1.e-6_dp
165 :
166 : CHARACTER(LEN=10) :: aname
167 : CHARACTER(LEN=2) :: element_symbol
168 : INTEGER :: i, ikind, nbas, nel, nkind, nsto, za
169 : LOGICAL :: at_end
170 : REAL(KIND=dp) :: ea
171 : REAL(KIND=dp), DIMENSION(10) :: al, cl
172 : TYPE(gto_basis_set_type), POINTER :: orb_basis
173 : TYPE(mp_para_env_type), POINTER :: para_env
174 5280 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
175 : TYPE(qs_kind_type), POINTER :: qs_kind
176 :
177 5280 : IF (gcp_env%do_gcp) THEN
178 4 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
179 66 : ALLOCATE (gcp_env%gcp_kind(nkind))
180 10 : DO ikind = 1, nkind
181 6 : qs_kind => qs_kind_set(ikind)
182 6 : gcp_env%gcp_kind(ikind)%rcsto = 0.0_dp
183 6 : CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
184 6 : CALL get_ptable_info(element_symbol, number=za)
185 6 : gcp_env%gcp_kind(ikind)%za = za
186 30 : gcp_env%gcp_kind(ikind)%asto = gcp_env%eta*SUM(sexp(1:4, za))/REAL(nll(za), KIND=dp)
187 6 : gcp_env%gcp_kind(ikind)%nq = nshell(za)
188 6 : gcp_env%gcp_kind(ikind)%rcsto = ((nshell(za) - 1)*2.5_dp - LOG(epsc))/gcp_env%gcp_kind(ikind)%asto
189 : ! basis
190 6 : NULLIFY (orb_basis)
191 6 : CALL get_qs_kind(qs_kind, basis_set=orb_basis, basis_type="ORB")
192 6 : CALL get_gto_basis_set(gto_basis_set=orb_basis, nsgf=nbas)
193 24 : nel = SUM(qs_kind%elec_conf)
194 6 : gcp_env%gcp_kind(ikind)%nbvirt = REAL(nbas, KIND=dp) - 0.5_dp*REAL(nel, KIND=dp)
195 : ! STO-nG
196 6 : nsto = SIZE(gcp_env%gcp_kind(ikind)%al)
197 6 : CALL get_sto_ng(gcp_env%gcp_kind(ikind)%asto, nsto, nshell(za), 0, al, cl)
198 58 : DO i = 1, nsto
199 36 : gcp_env%gcp_kind(ikind)%al(i) = al(i)
200 42 : gcp_env%gcp_kind(ikind)%cl(i) = cl(i)*(2._dp*al(i)/pi)**0.75_dp
201 : END DO
202 : END DO
203 : ! eamiss from data file
204 4 : IF (gcp_env%parameter_file_name /= "---") THEN
205 : BLOCK
206 : TYPE(cp_parser_type) :: parser
207 0 : CALL get_qs_env(qs_env, para_env=para_env)
208 0 : DO ikind = 1, nkind
209 0 : qs_kind => qs_kind_set(ikind)
210 0 : CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
211 0 : CALL get_ptable_info(element_symbol, number=za)
212 : !
213 0 : CALL parser_create(parser, gcp_env%parameter_file_name, para_env=para_env)
214 0 : ea = 0.0_dp
215 : DO
216 : at_end = .FALSE.
217 0 : CALL parser_get_next_line(parser, 1, at_end)
218 0 : IF (at_end) EXIT
219 0 : CALL parser_get_object(parser, aname)
220 0 : IF (TRIM(aname) == element_symbol) THEN
221 0 : CALL parser_get_object(parser, ea)
222 0 : EXIT
223 : END IF
224 : END DO
225 0 : CALL parser_release(parser)
226 0 : gcp_env%gcp_kind(ikind)%eamiss = ea
227 : END DO
228 : END BLOCK
229 : END IF
230 : !
231 : ! eamiss from input
232 4 : IF (ASSOCIATED(gcp_env%kind_type)) THEN
233 16 : DO i = 1, SIZE(gcp_env%kind_type)
234 12 : IF (TRIM(gcp_env%kind_type(i)) == "XX") CYCLE
235 12 : element_symbol = TRIM(gcp_env%kind_type(i))
236 12 : CALL get_ptable_info(element_symbol, number=za)
237 12 : ea = gcp_env%ea(i)
238 32 : DO ikind = 1, nkind
239 28 : IF (za == gcp_env%gcp_kind(ikind)%za) THEN
240 6 : gcp_env%gcp_kind(ikind)%eamiss = ea
241 : END IF
242 : END DO
243 : END DO
244 : END IF
245 : END IF
246 :
247 5280 : END SUBROUTINE qs_gcp_init
248 : ! **************************************************************************************************
249 :
250 : END MODULE qs_gcp_utils
251 :
|
