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 Writer for CASINO gwfn.data files.
10 : !> \par History
11 : !> 05.2026 created [Codex]
12 : ! **************************************************************************************************
13 : MODULE casino_utils
14 :
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE basis_set_types, ONLY: get_gto_basis_set,&
17 : gto_basis_set_type
18 : USE cell_types, ONLY: cell_type
19 : USE cp2k_info, ONLY: cp2k_version
20 : USE cp_blacs_env, ONLY: cp_blacs_env_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
24 : USE cp_files, ONLY: close_file,&
25 : open_file
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_release,&
28 : cp_fm_struct_type
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_get_submatrix,&
31 : cp_fm_release,&
32 : cp_fm_set_all,&
33 : cp_fm_to_fm_submat_general,&
34 : cp_fm_type
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_get_default_io_unit,&
37 : cp_logger_type
38 : USE external_potential_types, ONLY: get_potential,&
39 : gth_potential_type,&
40 : sgp_potential_type
41 : USE input_section_types, ONLY: section_vals_type,&
42 : section_vals_val_get
43 : USE kinds, ONLY: default_path_length,&
44 : default_string_length,&
45 : dp
46 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
47 : kpoint_init_cell_index,&
48 : kpoint_initialize,&
49 : kpoint_initialize_mo_set,&
50 : kpoint_initialize_mos
51 : USE kpoint_types, ONLY: get_kpoint_env,&
52 : get_kpoint_info,&
53 : kpoint_create,&
54 : kpoint_env_p_type,&
55 : kpoint_release,&
56 : kpoint_type
57 : USE mathconstants, ONLY: pi
58 : USE message_passing, ONLY: mp_para_env_type
59 : USE orbital_pointers, ONLY: nso
60 : USE particle_types, ONLY: particle_type
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_kind_types, ONLY: get_qs_kind,&
64 : get_qs_kind_set,&
65 : qs_kind_type
66 : USE qs_mo_types, ONLY: get_mo_set,&
67 : mo_set_type
68 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
69 : USE qs_scf_diagonalization, ONLY: do_general_diag_kp
70 : USE qs_scf_types, ONLY: qs_scf_env_type
71 : USE qs_wannier90, ONLY: prepare_wannier90_scf_mos
72 : USE scf_control_types, ONLY: scf_control_type
73 : USE string_utilities, ONLY: lowercase
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 :
78 : PRIVATE
79 :
80 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'casino_utils'
81 : INTEGER, PARAMETER, PRIVATE :: max_casino_l = 4
82 :
83 : PUBLIC :: write_casino
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief Write a CASINO gwfn.data file from the converged GPW/GAPW wavefunction.
89 : !> \param qs_env the QS environment
90 : !> \param casino_section the DFT%PRINT%CASINO input section
91 : ! **************************************************************************************************
92 10 : SUBROUTINE write_casino(qs_env, casino_section)
93 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
94 : TYPE(section_vals_type), INTENT(IN), POINTER :: casino_section
95 :
96 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_casino'
97 :
98 : CHARACTER(len=default_path_length) :: filename
99 : INTEGER :: ao_num, col_offset, handle, iao, iatom, ikind, ikp, ikp_loc, ikp_out, imo, ipgf, &
100 : iset, ishell, ishell_loc, ispin, iw, k, l, mo_num, nao_shell, natoms, nel_tot, &
101 : ngth_pseudo, nkp, nkp_mo, nkp_out, nmo, npseudo_atoms, nreal_k, nset, nsgf, nsgp_pseudo, &
102 : nspins, output_unit, periodicity, prim_num, shell_num, zatom
103 10 : INTEGER, ALLOCATABLE, DIMENSION(:) :: agauge, ao_to_atom, atomic_number, cp2k_to_casino_ao, &
104 10 : first_shell, kp_order, prim_per_shell, shell_ang_mom, shell_type
105 : INTEGER, DIMENSION(2) :: kp_range, nmo_spin
106 10 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
107 10 : INTEGER, DIMENSION(:, :), POINTER :: l_shell_set
108 : LOGICAL :: casino_kpoints_created, do_kpoints, &
109 : ionode, periodic, use_real_wfn, &
110 : write_pseudos
111 10 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: kp_real
112 : REAL(KIND=dp) :: cval, e_nn, eps_kpoint_real, kdotg, &
113 : pseudo_tol, sval, zeff
114 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coefficients, exponents, mo_scale, &
115 10 : valence_charge
116 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord, kvec, mo_energy, mos_sgf, &
117 10 : mos_sgf_im, shell_position
118 : REAL(KIND=dp), DIMENSION(3) :: scoord
119 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp, zetas
120 10 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
121 : TYPE(cell_type), POINTER :: cell
122 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
123 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
124 : TYPE(cp_fm_type) :: fm_dummy, fm_mo_coeff, fm_mo_coeff_im
125 : TYPE(cp_logger_type), POINTER :: logger
126 : TYPE(dft_control_type), POINTER :: dft_control
127 : TYPE(gth_potential_type), POINTER :: gth_potential
128 : TYPE(gto_basis_set_type), POINTER :: basis_set
129 10 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
130 : TYPE(kpoint_type), POINTER :: casino_kpoints, kpoints
131 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
132 10 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
133 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_inter_kp
134 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
135 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
136 : TYPE(sgp_potential_type), POINTER :: sgp_potential
137 :
138 10 : CALL timeset(routineN, handle)
139 :
140 10 : NULLIFY (basis_set, blacs_env, casino_kpoints, cell, dft_control, fm_struct, gcc, &
141 10 : gth_potential, kind_set, kp_env, kpoints, l_shell_set, logger, mos, mos_kp, npgf, &
142 10 : nshell, para_env, para_env_inter_kp, particle_set, sgp_potential, xkp, zetas)
143 :
144 10 : logger => cp_get_default_logger()
145 10 : output_unit = cp_logger_get_default_io_unit(logger)
146 :
147 10 : CPASSERT(ASSOCIATED(qs_env))
148 :
149 10 : CALL section_vals_val_get(casino_section, "FILENAME", c_val=filename)
150 10 : IF (LEN_TRIM(filename) == 0) filename = "gwfn.data"
151 10 : CALL section_vals_val_get(casino_section, "EPS_KPOINT_REAL", r_val=eps_kpoint_real)
152 10 : CALL section_vals_val_get(casino_section, "WRITE_PSEUDOPOTENTIALS", l_val=write_pseudos)
153 :
154 10 : CALL get_qs_env(qs_env, para_env=para_env)
155 10 : ionode = para_env%is_source()
156 :
157 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=kind_set, &
158 : natom=natoms, dft_control=dft_control, nelectron_total=nel_tot, &
159 10 : do_kpoints=do_kpoints, kpoints=kpoints, blacs_env=blacs_env)
160 10 : casino_kpoints => kpoints
161 : casino_kpoints_created = .FALSE.
162 : CALL prepare_casino_kpoint_grid(qs_env, casino_section, do_kpoints, kpoints, &
163 10 : casino_kpoints, casino_kpoints_created)
164 10 : nspins = dft_control%nspins
165 10 : IF (nspins > 2) CPABORT("CASINO gwfn.data supports at most two spin channels.")
166 :
167 40 : periodicity = COUNT(cell%perd /= 0)
168 10 : periodic = periodicity > 0
169 10 : pseudo_tol = 1.0E-8_dp
170 :
171 70 : ALLOCATE (coord(3, natoms), atomic_number(natoms), valence_charge(natoms))
172 10 : npseudo_atoms = 0
173 10 : ngth_pseudo = 0
174 10 : nsgp_pseudo = 0
175 28 : DO iatom = 1, natoms
176 72 : coord(:, iatom) = particle_set(iatom)%r(1:3)
177 18 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
178 : CALL get_qs_kind(kind_set(ikind), zatom=zatom, zeff=zeff, &
179 18 : gth_potential=gth_potential, sgp_potential=sgp_potential)
180 18 : IF (ABS(zeff) < pseudo_tol) zeff = REAL(zatom, KIND=dp)
181 18 : atomic_number(iatom) = zatom
182 18 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
183 8 : atomic_number(iatom) = zatom + 200
184 8 : npseudo_atoms = npseudo_atoms + 1
185 8 : IF (ASSOCIATED(gth_potential)) ngth_pseudo = ngth_pseudo + 1
186 8 : IF (ASSOCIATED(sgp_potential)) nsgp_pseudo = nsgp_pseudo + 1
187 10 : ELSE IF (ABS(zeff - REAL(zatom, KIND=dp)) > pseudo_tol) THEN
188 0 : atomic_number(iatom) = zatom + 200
189 0 : npseudo_atoms = npseudo_atoms + 1
190 : END IF
191 46 : valence_charge(iatom) = zeff
192 : END DO
193 :
194 10 : IF (ionode .AND. write_pseudos .AND. nsgp_pseudo > 0) THEN
195 1 : CALL write_casino_sgp_pseudopotentials(kind_set, particle_set, natoms, filename, output_unit)
196 : END IF
197 :
198 10 : IF (periodic) THEN
199 2 : CALL periodic_nuclear_repulsion_energy(cell, periodicity, coord, valence_charge, e_nn)
200 : ELSE
201 8 : CALL nuclear_repulsion_energy(particle_set, kind_set, e_nn)
202 : END IF
203 10 : e_nn = e_nn/REAL(natoms, KIND=dp)
204 :
205 10 : IF (do_kpoints) THEN
206 2 : CALL get_kpoint_info(casino_kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn)
207 : ELSE
208 8 : nkp = 1
209 8 : use_real_wfn = .TRUE.
210 : END IF
211 10 : nkp_mo = MERGE(nkp, 1, do_kpoints)
212 :
213 60 : ALLOCATE (kp_order(nkp_mo), kp_real(nkp_mo), kvec(3, nkp_mo))
214 : CALL build_kpoint_order(cell, periodic, do_kpoints, nkp, xkp, eps_kpoint_real, &
215 10 : kp_order, kp_real, nkp_out, nreal_k, kvec)
216 18 : IF (do_kpoints .AND. use_real_wfn .AND. ANY(.NOT. kp_real(1:nkp_out))) THEN
217 0 : CPABORT("CASINO complex k-points require CP2K complex k-point wavefunctions.")
218 : END IF
219 :
220 10 : CALL get_qs_kind_set(kind_set, nshell=shell_num, npgf_seg=prim_num, nsgf=nsgf)
221 10 : ao_num = nsgf
222 :
223 0 : ALLOCATE (shell_type(shell_num), prim_per_shell(shell_num), first_shell(natoms + 1), &
224 0 : shell_ang_mom(shell_num), shell_position(3, shell_num), &
225 0 : exponents(prim_num), coefficients(prim_num), ao_to_atom(ao_num), &
226 170 : cp2k_to_casino_ao(ao_num), mo_scale(ao_num))
227 :
228 10 : ishell = 0
229 10 : ipgf = 0
230 10 : iao = 0
231 28 : DO iatom = 1, natoms
232 18 : first_shell(iatom) = ishell + 1
233 18 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
234 18 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
235 : CALL get_gto_basis_set(basis_set, nset=nset, nshell=nshell, npgf=npgf, &
236 18 : zet=zetas, gcc=gcc, l=l_shell_set)
237 74 : DO iset = 1, nset
238 86 : DO ishell_loc = 1, nshell(iset)
239 40 : ishell = ishell + 1
240 40 : l = l_shell_set(ishell_loc, iset)
241 40 : IF (l > max_casino_l) THEN
242 0 : CPABORT("CASINO writer currently supports harmonic Gaussian shells up to g.")
243 : END IF
244 40 : shell_ang_mom(ishell) = l
245 40 : shell_type(ishell) = casino_shell_type(l)
246 40 : prim_per_shell(ishell) = npgf(iset)
247 160 : shell_position(:, ishell) = particle_set(iatom)%r(1:3)
248 : CALL casino_shell_coefficients(l, npgf(iset), zetas(1:npgf(iset), iset), &
249 : gcc(1:npgf(iset), ishell_loc, iset), &
250 : exponents(ipgf + 1:ipgf + npgf(iset)), &
251 40 : coefficients(ipgf + 1:ipgf + npgf(iset)))
252 40 : nao_shell = nso(l)
253 88 : DO k = 1, nao_shell
254 48 : cp2k_to_casino_ao(iao + k) = iao + casino_cp2k_index(l, k)
255 48 : mo_scale(iao + k) = casino_mo_scale(l, k)
256 88 : ao_to_atom(iao + k) = iatom
257 : END DO
258 40 : ipgf = ipgf + npgf(iset)
259 68 : iao = iao + nao_shell
260 : END DO
261 : END DO
262 : END DO
263 10 : first_shell(natoms + 1) = shell_num + 1
264 10 : CPASSERT(ishell == shell_num)
265 10 : CPASSERT(ipgf == prim_num)
266 10 : CPASSERT(iao == ao_num)
267 :
268 40 : ALLOCATE (mo_energy(ao_num, nkp_mo*nspins))
269 10 : mo_energy(:, :) = 0.0_dp
270 10 : nmo_spin(:) = 0
271 :
272 10 : IF (do_kpoints) THEN
273 2 : CALL get_kpoint_info(casino_kpoints, kp_env=kp_env, kp_range=kp_range, nkp=nkp)
274 2 : CALL get_kpoint_env(kp_env(1)%kpoint_env, mos=mos_kp)
275 4 : DO ispin = 1, nspins
276 2 : CALL get_mo_set(mos_kp(1, ispin), nmo=nmo)
277 2 : IF (nmo < ao_num) THEN
278 0 : CPABORT("CASINO gwfn.data requires a complete MO set. Increase ADDED_MOS.")
279 : END IF
280 6 : nmo_spin(ispin) = nmo
281 : END DO
282 6 : mo_num = nkp*SUM(nmo_spin)
283 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
284 2 : nrow_global=nsgf, ncol_global=mo_num)
285 2 : CALL cp_fm_create(fm_mo_coeff, fm_struct)
286 2 : CALL cp_fm_set_all(fm_mo_coeff, 0.0_dp)
287 2 : IF (.NOT. use_real_wfn) THEN
288 2 : CALL cp_fm_create(fm_mo_coeff_im, fm_struct)
289 2 : CALL cp_fm_set_all(fm_mo_coeff_im, 0.0_dp)
290 : END IF
291 2 : CALL cp_fm_struct_release(fm_struct)
292 :
293 4 : DO ispin = 1, nspins
294 8 : DO ikp = 1, nkp
295 4 : nmo = nmo_spin(ispin)
296 4 : col_offset = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
297 6 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
298 4 : ikp_loc = ikp - kp_range(1) + 1
299 4 : CALL get_kpoint_env(kp_env(ikp_loc)%kpoint_env, mos=mos_kp)
300 4 : IF (mos_kp(1, ispin)%use_mo_coeff_b) THEN
301 0 : CALL copy_dbcsr_to_fm(mos_kp(1, ispin)%mo_coeff_b, mos_kp(1, ispin)%mo_coeff)
302 : END IF
303 : CALL cp_fm_to_fm_submat_general(mos_kp(1, ispin)%mo_coeff, fm_mo_coeff, &
304 4 : nsgf, nmo, 1, 1, 1, col_offset + 1, blacs_env)
305 : mo_energy(1:ao_num, ikp + (ispin - 1)*nkp_mo) = &
306 20 : mos_kp(1, ispin)%eigenvalues(1:ao_num)
307 4 : IF (.NOT. use_real_wfn) THEN
308 4 : IF (mos_kp(2, ispin)%use_mo_coeff_b) THEN
309 0 : CALL copy_dbcsr_to_fm(mos_kp(2, ispin)%mo_coeff_b, mos_kp(2, ispin)%mo_coeff)
310 : END IF
311 : CALL cp_fm_to_fm_submat_general(mos_kp(2, ispin)%mo_coeff, fm_mo_coeff_im, &
312 4 : nsgf, nmo, 1, 1, 1, col_offset + 1, blacs_env)
313 : END IF
314 : ELSE
315 : CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff, &
316 0 : nsgf, nmo, 1, 1, 1, col_offset + 1, blacs_env)
317 0 : IF (.NOT. use_real_wfn) THEN
318 : CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff_im, &
319 0 : nsgf, nmo, 1, 1, 1, col_offset + 1, blacs_env)
320 : END IF
321 : END IF
322 : END DO
323 : END DO
324 2 : CALL get_kpoint_info(casino_kpoints, para_env_inter_kp=para_env_inter_kp)
325 2 : CALL para_env_inter_kp%sum(mo_energy)
326 : ELSE
327 8 : CALL get_qs_env(qs_env, mos=mos)
328 18 : DO ispin = 1, nspins
329 10 : CALL get_mo_set(mos(ispin), nmo=nmo)
330 10 : IF (nmo < ao_num) THEN
331 0 : CPABORT("CASINO gwfn.data requires a complete MO set. Increase ADDED_MOS.")
332 : END IF
333 10 : nmo_spin(ispin) = nmo
334 72 : mo_energy(1:ao_num, 1 + (ispin - 1)*nkp_mo) = mos(ispin)%eigenvalues(1:ao_num)
335 : END DO
336 : END IF
337 :
338 10 : IF (do_kpoints .AND. .NOT. use_real_wfn) THEN
339 6 : ALLOCATE (agauge(3*natoms))
340 6 : DO iatom = 1, natoms
341 52 : scoord(:) = MATMUL(cell%h_inv, particle_set(iatom)%r(1:3))
342 18 : agauge(3*(iatom - 1) + 1:3*iatom) = -NINT(scoord(:))
343 : END DO
344 : END IF
345 :
346 10 : IF (ionode) THEN
347 5 : IF (npseudo_atoms > 0) THEN
348 2 : WRITE (output_unit, "((T2,A,I0,A))") "CASINO| Marked ", npseudo_atoms, &
349 4 : " pseudopotential atoms in gwfn.data."
350 2 : IF (ngth_pseudo > 0) THEN
351 : WRITE (output_unit, "((T2,A))") &
352 1 : "CASINO| GTH pseudopotentials require matching external CASINO *_pp.data files."
353 : END IF
354 2 : IF (.NOT. write_pseudos) THEN
355 : WRITE (output_unit, "((T2,A))") &
356 1 : "CASINO| WRITE_PSEUDOPOTENTIALS is disabled; provide CASINO *_pp.data files manually."
357 : END IF
358 : END IF
359 5 : WRITE (output_unit, "((T2,A,A))") 'CASINO| Writing gwfn.data file ', TRIM(filename)
360 : CALL open_file(file_name=filename, file_status="REPLACE", file_action="WRITE", &
361 5 : file_form="FORMATTED", unit_number=iw)
362 : CALL write_casino_header(iw, periodicity, nspins, e_nn, nel_tot, natoms, coord, &
363 : atomic_number, valence_charge, cell, periodic, nkp_out, nreal_k, &
364 : kvec, shell_num, ao_num, prim_num, shell_ang_mom, shell_type, &
365 5 : prim_per_shell, first_shell, exponents, coefficients, shell_position)
366 : END IF
367 :
368 60 : ALLOCATE (mos_sgf(nsgf, ao_num), mos_sgf_im(nsgf, ao_num))
369 10 : mos_sgf(:, :) = 0.0_dp
370 10 : mos_sgf_im(:, :) = 0.0_dp
371 :
372 10 : IF (do_kpoints) THEN
373 4 : DO ispin = 1, nspins
374 6 : DO ikp_out = 1, nkp_out
375 2 : ikp = kp_order(ikp_out)
376 2 : col_offset = (ikp - 1)*nmo_spin(ispin) + (ispin - 1)*nmo_spin(1)*nkp
377 2 : CALL cp_fm_get_submatrix(fm_mo_coeff, mos_sgf, 1, col_offset + 1, nsgf, ao_num)
378 2 : IF (.NOT. use_real_wfn) THEN
379 2 : CALL cp_fm_get_submatrix(fm_mo_coeff_im, mos_sgf_im, 1, col_offset + 1, nsgf, ao_num)
380 10 : DO iao = 1, ao_num
381 8 : iatom = ao_to_atom(iao)
382 : kdotg = 2.0_dp*pi*DOT_PRODUCT(xkp(:, ikp), &
383 32 : REAL(agauge(3*(iatom - 1) + 1:3*iatom), KIND=dp))
384 8 : cval = COS(kdotg)
385 8 : sval = SIN(kdotg)
386 42 : DO imo = 1, ao_num
387 : CALL rotate_complex_pair(mos_sgf(cp2k_to_casino_ao(iao), imo), &
388 40 : mos_sgf_im(cp2k_to_casino_ao(iao), imo), cval, sval)
389 : END DO
390 : END DO
391 : ELSE
392 0 : mos_sgf_im(:, :) = 0.0_dp
393 : END IF
394 4 : IF (ionode) THEN
395 : CALL write_casino_orbitals(iw, mos_sgf, mos_sgf_im, cp2k_to_casino_ao, mo_scale, &
396 1 : ao_num,.NOT. kp_real(ikp_out))
397 : END IF
398 : END DO
399 : END DO
400 : ELSE
401 18 : DO ispin = 1, nspins
402 10 : IF (mos(ispin)%use_mo_coeff_b) THEN
403 0 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff)
404 : END IF
405 10 : CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, mos_sgf, 1, 1, nsgf, ao_num)
406 10 : mos_sgf_im(:, :) = 0.0_dp
407 18 : IF (ionode) THEN
408 : CALL write_casino_orbitals(iw, mos_sgf, mos_sgf_im, cp2k_to_casino_ao, mo_scale, &
409 5 : ao_num, .FALSE.)
410 : END IF
411 : END DO
412 : END IF
413 :
414 10 : IF (ionode) THEN
415 5 : WRITE (iw, *)
416 5 : IF (periodic) THEN
417 1 : WRITE (iw, '(A)') "EIGENVALUES"
418 1 : WRITE (iw, '(A)') "-----------"
419 2 : DO ikp_out = 1, nkp_out
420 1 : ikp = kp_order(ikp_out)
421 3 : DO ispin = 1, nspins
422 1 : IF (nspins == 1) THEN
423 1 : WRITE (iw, '(A,I6,3F14.8)') "k", ikp_out, kvec(:, ikp_out)
424 : ELSE
425 0 : WRITE (iw, '(A,I3,A,I6,3F14.8)') "spin", ispin, " k", ikp_out, kvec(:, ikp_out)
426 : END IF
427 2 : CALL write_real_vector(iw, mo_energy(1:ao_num, ikp + (ispin - 1)*nkp_mo))
428 : END DO
429 : END DO
430 : END IF
431 5 : CALL close_file(unit_number=iw)
432 : END IF
433 :
434 10 : DEALLOCATE (mos_sgf, mos_sgf_im)
435 10 : IF (do_kpoints) THEN
436 2 : CALL cp_fm_release(fm_mo_coeff)
437 2 : IF (.NOT. use_real_wfn) CALL cp_fm_release(fm_mo_coeff_im)
438 : END IF
439 10 : IF (casino_kpoints_created) CALL kpoint_release(casino_kpoints)
440 10 : IF (ALLOCATED(agauge)) DEALLOCATE (agauge)
441 0 : DEALLOCATE (ao_to_atom, atomic_number, coefficients, coord, cp2k_to_casino_ao, exponents, &
442 0 : first_shell, kp_order, kp_real, kvec, mo_energy, mo_scale, prim_per_shell, &
443 10 : shell_ang_mom, shell_position, shell_type, valence_charge)
444 :
445 10 : CALL timestop(handle)
446 60 : END SUBROUTINE write_casino
447 :
448 : ! **************************************************************************************************
449 : !> \brief Write CASINO pseudopotential files for the semilocal ECP kinds.
450 : !> \param kind_set the QS kinds
451 : !> \param particle_set the particle set
452 : !> \param natoms the number of atoms
453 : !> \param gwfn_filename the gwfn.data filename
454 : !> \param output_unit output unit for log messages
455 : ! **************************************************************************************************
456 1 : SUBROUTINE write_casino_sgp_pseudopotentials(kind_set, particle_set, natoms, gwfn_filename, output_unit)
457 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
458 : POINTER :: kind_set
459 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
460 : POINTER :: particle_set
461 : INTEGER, INTENT(IN) :: natoms
462 : CHARACTER(LEN=*), INTENT(IN) :: gwfn_filename
463 : INTEGER, INTENT(IN) :: output_unit
464 :
465 : CHARACTER(LEN=2) :: element_symbol
466 : INTEGER :: iatom, ikind, zatom
467 1 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: written
468 : REAL(KIND=dp) :: zeff
469 : TYPE(sgp_potential_type), POINTER :: sgp_potential
470 :
471 1 : NULLIFY (sgp_potential)
472 3 : ALLOCATE (written(SIZE(kind_set)))
473 1 : written(:) = .FALSE.
474 :
475 3 : DO iatom = 1, natoms
476 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, &
477 2 : kind_number=ikind, z=zatom)
478 2 : IF (written(ikind)) CYCLE
479 :
480 1 : CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential, zeff=zeff)
481 4 : IF (ASSOCIATED(sgp_potential)) THEN
482 1 : IF (ABS(zeff) < 1.0E-8_dp) zeff = REAL(zatom, KIND=dp)
483 : CALL write_casino_sgp_pseudopotential(sgp_potential, element_symbol, zatom, zeff, &
484 1 : gwfn_filename, output_unit)
485 1 : written(ikind) = .TRUE.
486 : END IF
487 : END DO
488 :
489 1 : DEALLOCATE (written)
490 1 : END SUBROUTINE write_casino_sgp_pseudopotentials
491 :
492 : ! **************************************************************************************************
493 : !> \brief Write a CASINO tabulated pseudopotential for a CP2K semilocal ECP.
494 : !> \param sgp_potential the CP2K semilocal Gaussian potential
495 : !> \param element_symbol the chemical symbol
496 : !> \param zatom the nuclear charge
497 : !> \param zeff the ECP valence charge
498 : !> \param gwfn_filename the gwfn.data filename
499 : !> \param output_unit output unit for log messages
500 : ! **************************************************************************************************
501 1 : SUBROUTINE write_casino_sgp_pseudopotential(sgp_potential, element_symbol, zatom, zeff, &
502 : gwfn_filename, output_unit)
503 : TYPE(sgp_potential_type), INTENT(IN), POINTER :: sgp_potential
504 : CHARACTER(LEN=*), INTENT(IN) :: element_symbol
505 : INTEGER, INTENT(IN) :: zatom
506 : REAL(KIND=dp), INTENT(IN) :: zeff
507 : CHARACTER(LEN=*), INTENT(IN) :: gwfn_filename
508 : INTEGER, INTENT(IN) :: output_unit
509 :
510 : CHARACTER(LEN=default_path_length) :: pp_filename
511 : INTEGER :: igrid, iw, l, local_l, ngrid, nloc, &
512 : nsemiloc, sl_lmax
513 : INTEGER, DIMENSION(0:10) :: npot
514 : LOGICAL :: ecp_local, ecp_semi_local, has_nlcc
515 : REAL(KIND=dp) :: agrid, bgrid, r, rmax, rv_local
516 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rgrid
517 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rpot
518 :
519 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, &
520 : ecp_semi_local=ecp_semi_local, nloc=nloc, sl_lmax=sl_lmax, &
521 1 : npot=npot, has_nlcc=has_nlcc)
522 1 : IF (.NOT. ecp_local .OR. nloc == 0) THEN
523 0 : WRITE (output_unit, "((T2,A,A,A))") "CASINO| Cannot write ", TRIM(element_symbol), &
524 0 : "_pp.data: only CP2K semilocal ECP potentials are supported."
525 : RETURN
526 : END IF
527 :
528 1 : local_l = MERGE(sl_lmax + 1, 0, ecp_semi_local)
529 1 : rmax = 100.0_dp
530 1 : agrid = 70.0_dp*EXP(-5.0_dp*LOG(10.0_dp))/REAL(zatom, KIND=dp)
531 1 : bgrid = 1.0_dp/70.0_dp
532 :
533 1 : ngrid = 0
534 831 : DO
535 832 : r = agrid*(EXP(bgrid*REAL(ngrid, KIND=dp)) - 1.0_dp)
536 832 : IF (r > rmax) EXIT
537 831 : ngrid = ngrid + 1
538 : END DO
539 :
540 6 : ALLOCATE (rgrid(ngrid), rpot(0:local_l, ngrid))
541 832 : DO igrid = 1, ngrid
542 831 : r = agrid*(EXP(bgrid*REAL(igrid - 1, KIND=dp)) - 1.0_dp)
543 831 : rgrid(igrid) = r
544 : rv_local = casino_sgp_r_times_v(nloc, sgp_potential%nrloc(1:nloc), &
545 : sgp_potential%bloc(1:nloc), &
546 831 : sgp_potential%aloc(1:nloc), r, zeff, .TRUE.)
547 2494 : DO l = 0, local_l
548 1662 : rpot(l, igrid) = rv_local
549 2493 : IF (l < local_l .AND. ecp_semi_local) THEN
550 831 : nsemiloc = npot(l)
551 831 : IF (nsemiloc > 0) THEN
552 : rpot(l, igrid) = rpot(l, igrid) + &
553 : casino_sgp_r_times_v(nsemiloc, sgp_potential%nrpot(1:nsemiloc, l), &
554 : sgp_potential%bpot(1:nsemiloc, l), &
555 831 : sgp_potential%apot(1:nsemiloc, l), r, zeff, .FALSE.)
556 : END IF
557 : END IF
558 : END DO
559 : END DO
560 2494 : rpot(:, :) = 2.0_dp*rpot(:, :)
561 :
562 1 : CALL casino_pp_filename(gwfn_filename, element_symbol, pp_filename)
563 : CALL open_file(file_name=pp_filename, file_status="REPLACE", file_action="WRITE", &
564 1 : file_form="FORMATTED", unit_number=iw)
565 1 : WRITE (iw, '(A)') "CP2K ECP pseudopotential in real space"
566 1 : WRITE (iw, '(A)') "Atomic number and pseudo-charge"
567 1 : WRITE (iw, '(I6,1X,F18.10)') zatom, zeff
568 1 : WRITE (iw, '(A)') "Energy units (rydberg/hartree/ev):"
569 1 : WRITE (iw, '(A)') "rydberg"
570 1 : WRITE (iw, '(A)') "Angular momentum of local component (0=s,1=p,2=d..)"
571 1 : WRITE (iw, '(I6)') local_l
572 1 : WRITE (iw, '(A)') "NLRULE override (1) VMC/DMC (2) config gen (0 ==> input/default value)"
573 1 : WRITE (iw, '(2I6)') 0, 0
574 1 : WRITE (iw, '(A)') "Number of grid points"
575 1 : WRITE (iw, '(I8)') ngrid
576 1 : WRITE (iw, '(A)') "R(i) in atomic units"
577 832 : DO igrid = 1, ngrid
578 832 : WRITE (iw, '(ES20.12)') rgrid(igrid)
579 : END DO
580 3 : DO l = 0, local_l
581 2 : WRITE (iw, '(A,I0,A)') "r*potential (L=", l, ") in Ry"
582 1665 : DO igrid = 1, ngrid
583 1664 : WRITE (iw, '(ES20.12)') rpot(l, igrid)
584 : END DO
585 : END DO
586 1 : CALL close_file(unit_number=iw)
587 :
588 1 : WRITE (output_unit, "((T2,A,A))") "CASINO| Wrote pseudopotential file ", TRIM(pp_filename)
589 1 : IF (has_nlcc) THEN
590 0 : WRITE (output_unit, "((T2,A,A,A))") "CASINO| NLCC terms for ", TRIM(element_symbol), &
591 0 : " are not represented in CASINO *_pp.data."
592 : END IF
593 :
594 1 : DEALLOCATE (rgrid, rpot)
595 1 : END SUBROUTINE write_casino_sgp_pseudopotential
596 :
597 : ! **************************************************************************************************
598 : !> \brief Return r times a CP2K semilocal Gaussian ECP channel in Hartree.
599 : !> \param nterm the number of Gaussian terms
600 : !> \param nr the CP2K r**(n-2) exponents
601 : !> \param gaussian_exponent the Gaussian exponents
602 : !> \param coefficient the Gaussian coefficients
603 : !> \param r the radial grid point
604 : !> \param zeff the ECP valence charge
605 : !> \param local_channel true for the local Coulomb-tailed channel
606 : !> \return r times the potential value
607 : ! **************************************************************************************************
608 1662 : FUNCTION casino_sgp_r_times_v(nterm, nr, gaussian_exponent, coefficient, r, zeff, local_channel) RESULT(r_times_v)
609 : INTEGER, INTENT(IN) :: nterm
610 : INTEGER, DIMENSION(:), INTENT(IN) :: nr
611 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gaussian_exponent, coefficient
612 : REAL(KIND=dp), INTENT(IN) :: r, zeff
613 : LOGICAL, INTENT(IN) :: local_channel
614 : REAL(KIND=dp) :: r_times_v
615 :
616 : INTEGER :: iterm
617 :
618 1662 : IF (r == 0.0_dp) THEN
619 1662 : r_times_v = 0.0_dp
620 : RETURN
621 : END IF
622 :
623 1660 : r_times_v = 0.0_dp
624 1660 : IF (local_channel) r_times_v = -zeff
625 4980 : DO iterm = 1, nterm
626 3320 : CPASSERT(nr(iterm) >= 1)
627 : r_times_v = r_times_v + coefficient(iterm)*r**(nr(iterm) - 1)* &
628 4980 : EXP(-gaussian_exponent(iterm)*r*r)
629 : END DO
630 : END FUNCTION casino_sgp_r_times_v
631 :
632 : ! **************************************************************************************************
633 : !> \brief Build the CASINO pseudopotential filename next to gwfn.data.
634 : !> \param gwfn_filename the gwfn.data filename
635 : !> \param element_symbol the chemical symbol
636 : !> \param pp_filename the CASINO pseudopotential filename
637 : ! **************************************************************************************************
638 1 : SUBROUTINE casino_pp_filename(gwfn_filename, element_symbol, pp_filename)
639 : CHARACTER(LEN=*), INTENT(IN) :: gwfn_filename, element_symbol
640 : CHARACTER(LEN=*), INTENT(OUT) :: pp_filename
641 :
642 : CHARACTER(LEN=2) :: symbol
643 : INTEGER :: slash
644 :
645 1 : symbol = ADJUSTL(element_symbol)
646 1 : CALL lowercase(symbol)
647 1 : slash = INDEX(TRIM(gwfn_filename), "/", BACK=.TRUE.)
648 1 : IF (slash > 0) THEN
649 0 : pp_filename = gwfn_filename(1:slash)//TRIM(symbol)//"_pp.data"
650 : ELSE
651 1 : pp_filename = TRIM(symbol)//"_pp.data"
652 : END IF
653 1 : END SUBROUTINE casino_pp_filename
654 :
655 : ! **************************************************************************************************
656 : !> \brief Prepare the k-point object used for CASINO export.
657 : !> \param qs_env the QS environment
658 : !> \param casino_section the CASINO print section
659 : !> \param do_kpoints true when the SCF used k-points
660 : !> \param kpoints_scf the converged SCF k-point object
661 : !> \param kpoints_out the k-point object to write
662 : !> \param created true if kpoints_out must be released by the caller
663 : ! **************************************************************************************************
664 14 : SUBROUTINE prepare_casino_kpoint_grid(qs_env, casino_section, do_kpoints, kpoints_scf, &
665 : kpoints_out, created)
666 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
667 : TYPE(section_vals_type), INTENT(IN), POINTER :: casino_section
668 : LOGICAL, INTENT(IN) :: do_kpoints
669 : TYPE(kpoint_type), POINTER :: kpoints_scf, kpoints_out
670 : LOGICAL, INTENT(OUT) :: created
671 :
672 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_casino_kpoint_grid'
673 :
674 : CHARACTER(LEN=default_string_length) :: kp_scheme, reuse_reason
675 : INTEGER :: aligned_blocks, aligned_max_size, &
676 : handle, nfull, output_unit
677 : INTEGER, DIMENSION(3) :: nkp_grid
678 10 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
679 : LOGICAL :: diis_step, full_grid, full_kpoint_grid, &
680 : gamma_centered, reuse_scf_mos, &
681 : reused_scf_mos, symmetry
682 : REAL(KIND=dp) :: aligned_min_svalue, eps_geo, wsum
683 : REAL(KIND=dp), DIMENSION(3) :: kp_shift
684 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp_source
685 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp_source
686 : TYPE(cell_type), POINTER :: cell
687 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
688 : TYPE(cp_logger_type), POINTER :: logger
689 10 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
690 : TYPE(dft_control_type), POINTER :: dft_control
691 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
692 : TYPE(mp_para_env_type), POINTER :: para_env
693 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
694 10 : POINTER :: sab_nl
695 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
696 : TYPE(qs_scf_env_type), POINTER :: scf_env
697 : TYPE(scf_control_type), POINTER :: scf_control
698 :
699 10 : CALL timeset(routineN, handle)
700 :
701 10 : created = .FALSE.
702 10 : kpoints_out => kpoints_scf
703 10 : NULLIFY (blacs_env, cell, cell_to_index, dft_control, logger, matrix_ks, matrix_s, mos, &
704 10 : para_env, particle_set, sab_nl, scf_control, scf_env, wkp_source, xkp_source)
705 :
706 10 : IF (.NOT. do_kpoints) THEN
707 8 : CALL timestop(handle)
708 8 : RETURN
709 : END IF
710 2 : CPASSERT(ASSOCIATED(kpoints_scf))
711 :
712 : CALL get_kpoint_info(kpoints_scf, kp_scheme=kp_scheme, symmetry=symmetry, &
713 : full_grid=full_grid, nkp_grid=nkp_grid, kp_shift=kp_shift, &
714 2 : gamma_centered=gamma_centered, eps_geo=eps_geo)
715 2 : IF (.NOT. symmetry .OR. full_grid) THEN
716 0 : CALL timestop(handle)
717 0 : RETURN
718 : END IF
719 :
720 2 : CALL section_vals_val_get(casino_section, "FULL_KPOINT_GRID", l_val=full_kpoint_grid)
721 2 : IF (.NOT. full_kpoint_grid) THEN
722 0 : CPABORT("CASINO export requires a full k-point grid. Use PRINT%CASINO%FULL_KPOINT_GRID.")
723 : END IF
724 :
725 2 : SELECT CASE (TRIM(kp_scheme))
726 : CASE ("MONKHORST-PACK", "MACDONALD", "GENERAL")
727 : ! supported below
728 : CASE DEFAULT
729 2 : CPABORT("CASINO%FULL_KPOINT_GRID supports only MONKHORST-PACK, MACDONALD, and GENERAL k-points.")
730 : END SELECT
731 :
732 2 : logger => cp_get_default_logger()
733 2 : output_unit = cp_logger_get_default_io_unit(logger)
734 2 : CALL section_vals_val_get(casino_section, "REUSE_SCF_MOS", l_val=reuse_scf_mos)
735 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, cell=cell, &
736 : particle_set=particle_set, mos=mos, dft_control=dft_control, &
737 : sab_orb=sab_nl, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
738 2 : scf_env=scf_env, scf_control=scf_control)
739 2 : CPASSERT(ASSOCIATED(para_env))
740 2 : CPASSERT(ASSOCIATED(blacs_env))
741 2 : CPASSERT(ASSOCIATED(cell))
742 2 : CPASSERT(ASSOCIATED(particle_set))
743 2 : CPASSERT(ASSOCIATED(mos))
744 2 : CPASSERT(ASSOCIATED(dft_control))
745 2 : CPASSERT(ASSOCIATED(sab_nl))
746 2 : CPASSERT(ASSOCIATED(matrix_ks))
747 2 : CPASSERT(ASSOCIATED(matrix_s))
748 2 : CPASSERT(ASSOCIATED(scf_env))
749 2 : CPASSERT(ASSOCIATED(scf_control))
750 :
751 2 : NULLIFY (kpoints_out)
752 2 : CALL kpoint_create(kpoints_out)
753 2 : kpoints_out%kp_scheme = kp_scheme
754 2 : kpoints_out%symmetry = .FALSE.
755 2 : kpoints_out%full_grid = .TRUE.
756 2 : kpoints_out%verbose = .FALSE.
757 2 : kpoints_out%use_real_wfn = .FALSE.
758 2 : kpoints_out%eps_geo = eps_geo
759 2 : kpoints_out%parallel_group_size = para_env%num_pe
760 :
761 4 : SELECT CASE (TRIM(kp_scheme))
762 : CASE ("MONKHORST-PACK", "MACDONALD")
763 8 : kpoints_out%nkp_grid(1:3) = nkp_grid(1:3)
764 8 : kpoints_out%kp_shift(1:3) = kp_shift(1:3)
765 2 : kpoints_out%gamma_centered = gamma_centered
766 2 : CALL kpoint_initialize(kpoints_out, particle_set, cell)
767 : CASE ("GENERAL")
768 0 : IF (.NOT. ASSOCIATED(kpoints_scf%xkp_input) .OR. &
769 : .NOT. ASSOCIATED(kpoints_scf%wkp_input)) THEN
770 0 : CPABORT("CASINO%FULL_KPOINT_GRID cannot recover the unreduced GENERAL k-point set.")
771 : END IF
772 0 : xkp_source => kpoints_scf%xkp_input
773 0 : wkp_source => kpoints_scf%wkp_input
774 0 : nfull = SIZE(wkp_source)
775 0 : wsum = SUM(wkp_source)
776 0 : IF (wsum <= 0.0_dp) CPABORT("CASINO%FULL_KPOINT_GRID found invalid GENERAL k-point weights.")
777 0 : kpoints_out%nkp = nfull
778 0 : ALLOCATE (kpoints_out%xkp(3, nfull), kpoints_out%wkp(nfull))
779 0 : kpoints_out%xkp(1:3, 1:nfull) = xkp_source(1:3, 1:nfull)
780 2 : kpoints_out%wkp(1:nfull) = wkp_source(1:nfull)/wsum
781 : END SELECT
782 :
783 2 : CALL kpoint_env_initialize(kpoints_out, para_env, blacs_env)
784 2 : CALL kpoint_initialize_mos(kpoints_out, mos)
785 2 : CALL kpoint_initialize_mo_set(kpoints_out)
786 2 : CALL kpoint_init_cell_index(kpoints_out, sab_nl, para_env, dft_control%nimages)
787 :
788 2 : reused_scf_mos = .FALSE.
789 2 : reuse_reason = ""
790 2 : aligned_blocks = 0
791 2 : aligned_max_size = 0
792 2 : aligned_min_svalue = 0.0_dp
793 2 : diis_step = .FALSE.
794 2 : IF (reuse_scf_mos) THEN
795 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints_scf, scf_env, scf_control, .FALSE., &
796 2 : diis_step)
797 2 : CALL get_kpoint_info(kpoints_out, cell_to_index=cell_to_index)
798 : CALL prepare_wannier90_scf_mos(kpoints_out, kpoints_scf, matrix_s, matrix_ks, &
799 : cell_to_index, sab_nl, para_env, reused_scf_mos, &
800 : reuse_reason, aligned_blocks, aligned_max_size, &
801 2 : aligned_min_svalue)
802 : END IF
803 2 : IF (reused_scf_mos) THEN
804 2 : IF (output_unit > 0) THEN
805 : WRITE (output_unit, '(T2,A)') &
806 1 : "CASINO| Reused SCF MO coefficients for the full k-point grid."
807 1 : IF (aligned_blocks > 0) THEN
808 : WRITE (output_unit, '(T2,A,I0,A,I0,A,ES10.3)') &
809 0 : "CASINO| Ritz-stabilized ", aligned_blocks, &
810 0 : " degenerate SCF MO subspace(s); largest block has ", aligned_max_size, &
811 0 : " band(s), min metric eigenvalue ", aligned_min_svalue
812 : END IF
813 : END IF
814 : ELSE
815 0 : IF (output_unit > 0) THEN
816 0 : IF (reuse_scf_mos) THEN
817 : WRITE (output_unit, '(T2,A,A)') &
818 0 : "CASINO| Could not reuse SCF MOs: ", TRIM(reuse_reason)
819 : END IF
820 : WRITE (output_unit, '(T2,A)') &
821 0 : "CASINO| Diagonalizing the full k-point grid for export."
822 : END IF
823 0 : diis_step = .FALSE.
824 : CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints_out, scf_env, scf_control, .FALSE., &
825 0 : diis_step)
826 : END IF
827 2 : created = .TRUE.
828 :
829 2 : CALL timestop(handle)
830 10 : END SUBROUTINE prepare_casino_kpoint_grid
831 :
832 : ! **************************************************************************************************
833 : !> \brief Build the CASINO k-point order with all real k-points first.
834 : !> \param cell ...
835 : !> \param periodic ...
836 : !> \param do_kpoints ...
837 : !> \param nkp_total ...
838 : !> \param xkp ...
839 : !> \param eps_kpoint_real ...
840 : !> \param kp_order ...
841 : !> \param kp_real ...
842 : !> \param nkp_out ...
843 : !> \param nreal_k ...
844 : !> \param kvec ...
845 : ! **************************************************************************************************
846 10 : SUBROUTINE build_kpoint_order(cell, periodic, do_kpoints, nkp_total, xkp, eps_kpoint_real, &
847 10 : kp_order, kp_real, nkp_out, nreal_k, kvec)
848 : TYPE(cell_type), INTENT(IN), POINTER :: cell
849 : LOGICAL, INTENT(IN) :: periodic, do_kpoints
850 : INTEGER, INTENT(IN) :: nkp_total
851 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
852 : OPTIONAL, POINTER :: xkp
853 : REAL(KIND=dp), INTENT(IN) :: eps_kpoint_real
854 : INTEGER, DIMENSION(:), INTENT(OUT) :: kp_order
855 : LOGICAL, DIMENSION(:), INTENT(OUT) :: kp_real
856 : INTEGER, INTENT(OUT) :: nkp_out, nreal_k
857 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kvec
858 :
859 : INTEGER :: ikp, jkp, ncomplex_k
860 20 : INTEGER, DIMENSION(nkp_total) :: complex_order, real_order
861 8 : LOGICAL, DIMENSION(nkp_total) :: used
862 :
863 10 : IF (.NOT. periodic) THEN
864 8 : kp_order(1) = 1
865 8 : kp_real(1) = .TRUE.
866 8 : nkp_out = 1
867 8 : nreal_k = 1
868 32 : kvec(:, 1) = 0.0_dp
869 : RETURN
870 : END IF
871 :
872 2 : IF (.NOT. do_kpoints) THEN
873 0 : kp_order(1) = 1
874 0 : kp_real(1) = .TRUE.
875 0 : nkp_out = 1
876 0 : nreal_k = 1
877 0 : kvec(:, 1) = 0.0_dp
878 : RETURN
879 : END IF
880 :
881 6 : used(:) = .FALSE.
882 2 : nreal_k = 0
883 2 : ncomplex_k = 0
884 :
885 6 : DO ikp = 1, nkp_total
886 6 : IF (is_real_kpoint(cell, xkp(:, ikp), eps_kpoint_real)) THEN
887 0 : used(ikp) = .TRUE.
888 0 : nreal_k = nreal_k + 1
889 0 : real_order(nreal_k) = ikp
890 : END IF
891 : END DO
892 :
893 6 : DO ikp = 1, nkp_total
894 6 : IF (.NOT. used(ikp)) THEN
895 2 : used(ikp) = .TRUE.
896 2 : ncomplex_k = ncomplex_k + 1
897 2 : complex_order(ncomplex_k) = ikp
898 2 : DO jkp = ikp + 1, nkp_total
899 2 : IF (.NOT. used(jkp)) THEN
900 2 : IF (is_conjugate_kpoint(cell, xkp(:, ikp), xkp(:, jkp), eps_kpoint_real)) THEN
901 2 : used(jkp) = .TRUE.
902 2 : EXIT
903 : END IF
904 : END IF
905 : END DO
906 : END IF
907 : END DO
908 :
909 2 : nkp_out = nreal_k + ncomplex_k
910 2 : DO ikp = 1, nreal_k
911 0 : kp_order(ikp) = real_order(ikp)
912 0 : kp_real(ikp) = .TRUE.
913 2 : kvec(:, ikp) = 2.0_dp*pi*MATMUL(TRANSPOSE(cell%h_inv), xkp(:, real_order(ikp)))
914 : END DO
915 4 : DO ikp = 1, ncomplex_k
916 2 : jkp = nreal_k + ikp
917 2 : kp_order(jkp) = complex_order(ikp)
918 2 : kp_real(jkp) = .FALSE.
919 10 : kvec(:, jkp) = 2.0_dp*pi*MATMUL(TRANSPOSE(cell%h_inv), xkp(:, complex_order(ikp)))
920 : END DO
921 2 : END SUBROUTINE build_kpoint_order
922 :
923 : ! **************************************************************************************************
924 : !> \brief Returns true for conjugate k-points modulo reciprocal lattice vectors.
925 : !> \param cell ...
926 : !> \param xk1 ...
927 : !> \param xk2 ...
928 : !> \param eps_kpoint_real ...
929 : !> \return ...
930 : ! **************************************************************************************************
931 2 : FUNCTION is_conjugate_kpoint(cell, xk1, xk2, eps_kpoint_real) RESULT(is_conjugate)
932 : TYPE(cell_type), INTENT(IN), POINTER :: cell
933 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xk1, xk2
934 : REAL(KIND=dp), INTENT(IN) :: eps_kpoint_real
935 : LOGICAL :: is_conjugate
936 :
937 : INTEGER :: idir
938 : REAL(KIND=dp) :: reduced
939 :
940 2 : is_conjugate = .TRUE.
941 8 : DO idir = 1, 3
942 8 : IF (cell%perd(idir) /= 0) THEN
943 6 : reduced = xk1(idir) + xk2(idir)
944 6 : reduced = reduced - REAL(NINT(reduced), KIND=dp)
945 6 : IF (ABS(reduced) > eps_kpoint_real) THEN
946 : is_conjugate = .FALSE.
947 : EXIT
948 : END IF
949 : END IF
950 : END DO
951 2 : END FUNCTION is_conjugate_kpoint
952 :
953 : ! **************************************************************************************************
954 : !> \brief Returns true for Gamma/BZ-edge k-points where real Bloch orbitals can be used.
955 : !> \param cell ...
956 : !> \param xk ...
957 : !> \param eps_kpoint_real ...
958 : !> \return ...
959 : ! **************************************************************************************************
960 4 : FUNCTION is_real_kpoint(cell, xk, eps_kpoint_real) RESULT(is_real)
961 : TYPE(cell_type), INTENT(IN), POINTER :: cell
962 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xk
963 : REAL(KIND=dp), INTENT(IN) :: eps_kpoint_real
964 : LOGICAL :: is_real
965 :
966 : INTEGER :: idir
967 : REAL(KIND=dp) :: reduced
968 :
969 4 : is_real = .TRUE.
970 16 : DO idir = 1, 3
971 16 : IF (cell%perd(idir) /= 0) THEN
972 12 : reduced = xk(idir) - REAL(NINT(xk(idir)), KIND=dp)
973 4 : IF (ABS(reduced) > eps_kpoint_real .AND. &
974 16 : ABS(ABS(reduced) - 0.5_dp) > eps_kpoint_real) is_real = .FALSE.
975 : END IF
976 : END DO
977 4 : END FUNCTION is_real_kpoint
978 :
979 : ! **************************************************************************************************
980 : !> \brief Write all non-orbital CASINO gwfn.data sections.
981 : !> \param iw ...
982 : !> \param periodicity ...
983 : !> \param nspins ...
984 : !> \param e_nn ...
985 : !> \param nel_tot ...
986 : !> \param natoms ...
987 : !> \param coord ...
988 : !> \param atomic_number ...
989 : !> \param valence_charge ...
990 : !> \param cell ...
991 : !> \param periodic ...
992 : !> \param nkp ...
993 : !> \param nreal_k ...
994 : !> \param kvec ...
995 : !> \param shell_num ...
996 : !> \param ao_num ...
997 : !> \param prim_num ...
998 : !> \param shell_ang_mom ...
999 : !> \param shell_type ...
1000 : !> \param prim_per_shell ...
1001 : !> \param first_shell ...
1002 : !> \param exponents ...
1003 : !> \param coefficients ...
1004 : !> \param shell_position ...
1005 : ! **************************************************************************************************
1006 15 : SUBROUTINE write_casino_header(iw, periodicity, nspins, e_nn, nel_tot, natoms, coord, &
1007 5 : atomic_number, valence_charge, cell, periodic, nkp, nreal_k, kvec, &
1008 10 : shell_num, ao_num, prim_num, shell_ang_mom, shell_type, prim_per_shell, &
1009 5 : first_shell, exponents, coefficients, shell_position)
1010 : INTEGER, INTENT(IN) :: iw, periodicity, nspins
1011 : REAL(KIND=dp), INTENT(IN) :: e_nn
1012 : INTEGER, INTENT(IN) :: nel_tot, natoms
1013 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coord
1014 : INTEGER, DIMENSION(:), INTENT(IN) :: atomic_number
1015 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: valence_charge
1016 : TYPE(cell_type), INTENT(IN), POINTER :: cell
1017 : LOGICAL, INTENT(IN) :: periodic
1018 : INTEGER, INTENT(IN) :: nkp, nreal_k
1019 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kvec
1020 : INTEGER, INTENT(IN) :: shell_num, ao_num, prim_num
1021 : INTEGER, DIMENSION(:), INTENT(IN) :: shell_ang_mom, shell_type, &
1022 : prim_per_shell, first_shell
1023 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exponents, coefficients
1024 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: shell_position
1025 :
1026 : INTEGER :: highest_ang_mom, i
1027 :
1028 25 : highest_ang_mom = MAXVAL(shell_ang_mom) + 1
1029 :
1030 5 : WRITE (iw, '(A)') "CP2K CASINO gwfn.data"
1031 5 : WRITE (iw, *)
1032 5 : WRITE (iw, '(A)') "BASIC INFO"
1033 5 : WRITE (iw, '(A)') "----------"
1034 5 : WRITE (iw, '(A)') "Generated by:"
1035 5 : WRITE (iw, '(1X,A)') TRIM(cp2k_version)
1036 5 : WRITE (iw, '(A)') "Method:"
1037 5 : WRITE (iw, '(A)') " DFT"
1038 5 : WRITE (iw, '(A)') "DFT functional:"
1039 5 : WRITE (iw, '(A)') " CP2K"
1040 5 : WRITE (iw, '(A)') "Periodicity:"
1041 5 : WRITE (iw, '(1X,I0)') periodicity
1042 5 : WRITE (iw, '(A)') "Spin unrestricted:"
1043 9 : WRITE (iw, '(1X,A)') MERGE(".true. ", ".false.", nspins > 1)
1044 5 : WRITE (iw, '(A)') "Nuclear repulsion energy (au/atom):"
1045 5 : WRITE (iw, '(1PE20.13)') e_nn
1046 5 : WRITE (iw, '(A)') "Number of electrons per primitive cell"
1047 5 : WRITE (iw, '(1X,I0)') nel_tot
1048 5 : WRITE (iw, *)
1049 :
1050 5 : WRITE (iw, '(A)') "GEOMETRY"
1051 5 : WRITE (iw, '(A)') "--------"
1052 5 : WRITE (iw, '(A)') "Number of atoms"
1053 5 : WRITE (iw, '(1X,I0)') natoms
1054 5 : WRITE (iw, '(A)') "Atomic positions (au)"
1055 14 : DO i = 1, natoms
1056 14 : WRITE (iw, '(3(1PE20.13))') coord(:, i)
1057 : END DO
1058 5 : WRITE (iw, '(A)') "Atomic numbers for each atom"
1059 5 : CALL write_integer_vector(iw, atomic_number)
1060 5 : WRITE (iw, '(A)') "Valence charges for each atom"
1061 5 : CALL write_real_vector(iw, valence_charge)
1062 5 : IF (.NOT. periodic) WRITE (iw, *)
1063 5 : IF (periodic) THEN
1064 1 : WRITE (iw, '(A)') "Primitive lattice vectors (au)"
1065 4 : DO i = 1, 3
1066 13 : WRITE (iw, '(3(1PE20.13))') cell%hmat(:, i)
1067 : END DO
1068 1 : WRITE (iw, *)
1069 :
1070 1 : WRITE (iw, '(A)') "K SPACE NET"
1071 1 : WRITE (iw, '(A)') "-----------"
1072 1 : WRITE (iw, '(A)') "Number of k points"
1073 1 : WRITE (iw, '(1X,I0)') nkp
1074 1 : WRITE (iw, '(A)') "Number of 'real' k points on BZ edge"
1075 1 : WRITE (iw, '(1X,I0)') nreal_k
1076 1 : WRITE (iw, '(A)') "k point coordinates (au)"
1077 2 : DO i = 1, nkp
1078 2 : WRITE (iw, '(3(1PE20.13))') kvec(:, i)
1079 : END DO
1080 1 : WRITE (iw, *)
1081 : END IF
1082 :
1083 5 : WRITE (iw, '(A)') "BASIS SET"
1084 5 : WRITE (iw, '(A)') "---------"
1085 5 : WRITE (iw, '(A)') "Number of Gaussian centres"
1086 5 : WRITE (iw, '(1X,I0)') natoms
1087 5 : WRITE (iw, '(A)') "Number of shells per primitive cell"
1088 5 : WRITE (iw, '(1X,I0)') shell_num
1089 5 : WRITE (iw, '(A)') "Number of basis functions ('AO') per primitive cell"
1090 5 : WRITE (iw, '(1X,I0)') ao_num
1091 5 : WRITE (iw, '(A)') "Number of Gaussian primitives per primitive cell"
1092 5 : WRITE (iw, '(1X,I0)') prim_num
1093 5 : WRITE (iw, '(A)') "Highest shell angular momentum (s/p/d/f/g... 1/2/3/4/5...)"
1094 5 : WRITE (iw, '(1X,I0)') highest_ang_mom
1095 5 : WRITE (iw, '(A)') "Code for shell types (s/sp/p/d/f... 1/2/3/4/5...)"
1096 5 : CALL write_integer_vector(iw, shell_type)
1097 5 : WRITE (iw, '(A)') "Number of primitive Gaussians in each shell"
1098 5 : CALL write_integer_vector(iw, prim_per_shell)
1099 5 : WRITE (iw, '(A)') "Sequence number of first shell on each centre"
1100 5 : CALL write_integer_vector(iw, first_shell)
1101 5 : WRITE (iw, '(A)') "Exponents of Gaussian primitives"
1102 5 : CALL write_real_vector(iw, exponents)
1103 5 : WRITE (iw, '(A)') "Correctly normalised contraction coefficients"
1104 5 : CALL write_real_vector(iw, coefficients)
1105 5 : WRITE (iw, '(A)') "Position of each shell (au)"
1106 25 : DO i = 1, shell_num
1107 25 : WRITE (iw, '(3(1PE20.13))') shell_position(:, i)
1108 : END DO
1109 5 : WRITE (iw, *)
1110 :
1111 5 : WRITE (iw, '(A)') "MULTIDETERMINANT INFORMATION"
1112 5 : WRITE (iw, '(A)') "----------------------------"
1113 5 : WRITE (iw, '(A)') "GS"
1114 5 : WRITE (iw, *)
1115 5 : WRITE (iw, '(A)') "ORBITAL COEFFICIENTS"
1116 5 : WRITE (iw, '(A)') "---------------------------"
1117 5 : END SUBROUTINE write_casino_header
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief Write one CASINO MO block for a spin/k-point.
1121 : !> \param iw ...
1122 : !> \param mos_sgf ...
1123 : !> \param mos_sgf_im ...
1124 : !> \param cp2k_to_casino_ao ...
1125 : !> \param mo_scale ...
1126 : !> \param ao_num ...
1127 : !> \param complex_orbitals ...
1128 : ! **************************************************************************************************
1129 6 : SUBROUTINE write_casino_orbitals(iw, mos_sgf, mos_sgf_im, cp2k_to_casino_ao, mo_scale, ao_num, &
1130 : complex_orbitals)
1131 : INTEGER, INTENT(IN) :: iw
1132 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: mos_sgf, mos_sgf_im
1133 : INTEGER, DIMENSION(:), INTENT(IN) :: cp2k_to_casino_ao
1134 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: mo_scale
1135 : INTEGER, INTENT(IN) :: ao_num
1136 : LOGICAL, INTENT(IN) :: complex_orbitals
1137 :
1138 : INTEGER :: iao, imo, nbuffer
1139 : REAL(KIND=dp), DIMENSION(4) :: buffer
1140 :
1141 6 : nbuffer = 0
1142 6 : buffer(:) = 0.0_dp
1143 32 : DO imo = 1, ao_num
1144 188 : DO iao = 1, ao_num
1145 156 : CALL push_real(iw, buffer, nbuffer, mo_scale(iao)*mos_sgf(cp2k_to_casino_ao(iao), imo))
1146 182 : IF (complex_orbitals) THEN
1147 16 : CALL push_real(iw, buffer, nbuffer, mo_scale(iao)*mos_sgf_im(cp2k_to_casino_ao(iao), imo))
1148 : END IF
1149 : END DO
1150 : END DO
1151 6 : IF (nbuffer > 0) WRITE (iw, '(4(1PE20.13))') buffer(1:nbuffer)
1152 6 : END SUBROUTINE write_casino_orbitals
1153 :
1154 : ! **************************************************************************************************
1155 : !> \brief Append one real number to a four-column output buffer.
1156 : !> \param iw ...
1157 : !> \param buffer ...
1158 : !> \param nbuffer ...
1159 : !> \param value ...
1160 : ! **************************************************************************************************
1161 172 : SUBROUTINE push_real(iw, buffer, nbuffer, value)
1162 : INTEGER, INTENT(IN) :: iw
1163 : REAL(KIND=dp), DIMENSION(4), INTENT(INOUT) :: buffer
1164 : INTEGER, INTENT(INOUT) :: nbuffer
1165 : REAL(KIND=dp), INTENT(IN) :: value
1166 :
1167 172 : nbuffer = nbuffer + 1
1168 172 : buffer(nbuffer) = value
1169 172 : IF (nbuffer == SIZE(buffer)) THEN
1170 43 : WRITE (iw, '(4(1PE20.13))') buffer
1171 43 : nbuffer = 0
1172 : END IF
1173 172 : END SUBROUTINE push_real
1174 :
1175 : ! **************************************************************************************************
1176 : !> \brief Write an integer vector in CASINO-friendly fixed-width columns.
1177 : !> \param iw ...
1178 : !> \param values ...
1179 : ! **************************************************************************************************
1180 20 : SUBROUTINE write_integer_vector(iw, values)
1181 : INTEGER, INTENT(IN) :: iw
1182 : INTEGER, DIMENSION(:), INTENT(IN) :: values
1183 :
1184 : INTEGER :: i, ilast
1185 :
1186 40 : DO i = 1, SIZE(values), 8
1187 20 : ilast = MIN(i + 7, SIZE(values))
1188 40 : WRITE (iw, '(8I10)') values(i:ilast)
1189 : END DO
1190 20 : END SUBROUTINE write_integer_vector
1191 :
1192 : ! **************************************************************************************************
1193 : !> \brief Write a real vector in CASINO-friendly fixed-width columns.
1194 : !> \param iw ...
1195 : !> \param values ...
1196 : ! **************************************************************************************************
1197 16 : SUBROUTINE write_real_vector(iw, values)
1198 : INTEGER, INTENT(IN) :: iw
1199 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
1200 :
1201 : INTEGER :: i, ilast
1202 :
1203 62 : DO i = 1, SIZE(values), 4
1204 46 : ilast = MIN(i + 3, SIZE(values))
1205 62 : WRITE (iw, '(4(1PE20.13))') values(i:ilast)
1206 : END DO
1207 16 : END SUBROUTINE write_real_vector
1208 :
1209 : ! **************************************************************************************************
1210 : !> \brief Rotate a complex value by exp(-i*k.g) using the TREXIO gauge convention.
1211 : !> \param re ...
1212 : !> \param im ...
1213 : !> \param cval ...
1214 : !> \param sval ...
1215 : ! **************************************************************************************************
1216 32 : SUBROUTINE rotate_complex_pair(re, im, cval, sval)
1217 : REAL(KIND=dp), INTENT(INOUT) :: re, im
1218 : REAL(KIND=dp), INTENT(IN) :: cval, sval
1219 :
1220 : REAL(KIND=dp) :: im_old, re_old
1221 :
1222 32 : re_old = re
1223 32 : im_old = im
1224 32 : re = cval*re_old + sval*im_old
1225 32 : im = -sval*re_old + cval*im_old
1226 32 : END SUBROUTINE rotate_complex_pair
1227 :
1228 : ! **************************************************************************************************
1229 : !> \brief Convert CP2K normalized primitive data to CASINO contraction coefficients.
1230 : !> \param l ...
1231 : !> \param nprim ...
1232 : !> \param zetas ...
1233 : !> \param gcc ...
1234 : !> \param exponents ...
1235 : !> \param coefficients ...
1236 : ! **************************************************************************************************
1237 40 : SUBROUTINE casino_shell_coefficients(l, nprim, zetas, gcc, exponents, coefficients)
1238 : INTEGER, INTENT(IN) :: l, nprim
1239 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetas, gcc
1240 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: exponents, coefficients
1241 :
1242 : INTEGER :: i
1243 : REAL(KIND=dp) :: contraction_norm, expzet, prefac, &
1244 : prim_cart_fac
1245 80 : REAL(KIND=dp), DIMENSION(nprim) :: raw_coeff
1246 :
1247 40 : expzet = 0.25_dp*REAL(2*l + 3, KIND=dp)
1248 40 : prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
1249 186 : DO i = 1, nprim
1250 146 : prim_cart_fac = prefac*zetas(i)**expzet
1251 186 : raw_coeff(i) = gcc(i)/prim_cart_fac
1252 : END DO
1253 40 : contraction_norm = casino_contraction_norm(l, nprim, zetas, raw_coeff)
1254 186 : DO i = 1, nprim
1255 146 : exponents(i) = zetas(i)
1256 186 : coefficients(i) = raw_coeff(i)*contraction_norm*casino_primitive_norm(l, zetas(i))
1257 : END DO
1258 40 : END SUBROUTINE casino_shell_coefficients
1259 :
1260 : ! **************************************************************************************************
1261 : !> \brief Whole-contraction normalization used by CASINO's molden2qmc converter.
1262 : !> \param l ...
1263 : !> \param nprim ...
1264 : !> \param zetas ...
1265 : !> \param raw_coeff ...
1266 : !> \return ...
1267 : ! **************************************************************************************************
1268 40 : FUNCTION casino_contraction_norm(l, nprim, zetas, raw_coeff) RESULT(norm)
1269 : INTEGER, INTENT(IN) :: l, nprim
1270 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetas, raw_coeff
1271 : REAL(KIND=dp) :: norm
1272 :
1273 : INTEGER :: i, j
1274 : REAL(KIND=dp) :: overlap
1275 :
1276 40 : overlap = 0.0_dp
1277 186 : DO i = 1, nprim
1278 952 : DO j = 1, nprim
1279 : overlap = overlap + raw_coeff(i)*raw_coeff(j)* &
1280 912 : (2.0_dp*SQRT(zetas(i)*zetas(j))/(zetas(i) + zetas(j)))**(l + 1.5_dp)
1281 : END DO
1282 : END DO
1283 40 : norm = 1.0_dp/SQRT(overlap)
1284 40 : END FUNCTION casino_contraction_norm
1285 :
1286 : ! **************************************************************************************************
1287 : !> \brief Primitive m-independent normalization used by CASINO's Gaussian evaluator.
1288 : !> \param l ...
1289 : !> \param alpha ...
1290 : !> \return ...
1291 : ! **************************************************************************************************
1292 146 : FUNCTION casino_primitive_norm(l, alpha) RESULT(norm)
1293 : INTEGER, INTENT(IN) :: l
1294 : REAL(KIND=dp), INTENT(IN) :: alpha
1295 : REAL(KIND=dp) :: norm
1296 :
1297 146 : norm = SQRT(2.0_dp**(l + 1.5_dp)*alpha**(l + 1.5_dp))/pi**0.75_dp
1298 146 : IF (l > 0) norm = norm*SQRT(2.0_dp**l/odd_double_factorial(2*l - 1))
1299 146 : END FUNCTION casino_primitive_norm
1300 :
1301 : ! **************************************************************************************************
1302 : !> \brief CASINO shell type code.
1303 : !> \param l ...
1304 : !> \return ...
1305 : ! **************************************************************************************************
1306 40 : FUNCTION casino_shell_type(l) RESULT(shell_type)
1307 : INTEGER, INTENT(IN) :: l
1308 : INTEGER :: shell_type
1309 :
1310 40 : IF (l == 0) THEN
1311 : shell_type = 1
1312 : ELSE
1313 4 : shell_type = l + 2
1314 : END IF
1315 40 : END FUNCTION casino_shell_type
1316 :
1317 : ! **************************************************************************************************
1318 : !> \brief CP2K AO index for a CASINO/MOLDEN ordered harmonic shell.
1319 : !> \param l ...
1320 : !> \param k ...
1321 : !> \return ...
1322 : ! **************************************************************************************************
1323 48 : FUNCTION casino_cp2k_index(l, k) RESULT(idx)
1324 : INTEGER, INTENT(IN) :: l, k
1325 : INTEGER :: idx
1326 :
1327 : INTEGER, DIMENSION(9, 0:max_casino_l), PARAMETER :: map = RESHAPE([1, 0, 0, 0, 0, 0, 0, 0, 0 &
1328 : , 3, 1, 2, 0, 0, 0, 0, 0, 0, 3, 4, 2, 5, 1, 0, 0, 0, 0, 4, 5, 3, 6, 2, 7, 1, 0, 0, 5, 6, 4&
1329 : , 7, 3, 8, 2, 9, 1], [9, max_casino_l + 1])
1330 :
1331 48 : idx = map(k, l)
1332 48 : END FUNCTION casino_cp2k_index
1333 :
1334 : ! **************************************************************************************************
1335 : !> \brief Scale factors converting MOLDEN harmonic MO coefficients to CASINO conventions.
1336 : !> \param l ...
1337 : !> \param k ...
1338 : !> \return ...
1339 : ! **************************************************************************************************
1340 48 : FUNCTION casino_mo_scale(l, k) RESULT(scale)
1341 : INTEGER, INTENT(IN) :: l, k
1342 : REAL(KIND=dp) :: scale
1343 :
1344 : REAL(KIND=dp), DIMENSION(5), PARAMETER :: d_factor = [0.5_dp, 3.0_dp, 3.0_dp, 3.0_dp, 6.0_dp]
1345 :
1346 : INTEGER :: m
1347 :
1348 48 : IF (l <= 1) THEN
1349 : scale = 1.0_dp
1350 : ELSE
1351 0 : m = casino_m_quantum_number(k)
1352 0 : scale = casino_m_dependent_factor(l, m)
1353 0 : IF (l == 2) scale = scale*d_factor(k)
1354 : END IF
1355 48 : END FUNCTION casino_mo_scale
1356 :
1357 : ! **************************************************************************************************
1358 : !> \brief m sequence in CASINO/MOLDEN harmonic order: 0,+1,-1,+2,-2,...
1359 : !> \param k ...
1360 : !> \return ...
1361 : ! **************************************************************************************************
1362 0 : FUNCTION casino_m_quantum_number(k) RESULT(m)
1363 : INTEGER, INTENT(IN) :: k
1364 : INTEGER :: m
1365 :
1366 0 : IF (k == 1) THEN
1367 : m = 0
1368 0 : ELSE IF (MOD(k, 2) == 0) THEN
1369 0 : m = k/2
1370 : ELSE
1371 0 : m = -(k/2)
1372 : END IF
1373 0 : END FUNCTION casino_m_quantum_number
1374 :
1375 : ! **************************************************************************************************
1376 : !> \brief CASINO m-dependent normalization factor.
1377 : !> \param l ...
1378 : !> \param m ...
1379 : !> \return ...
1380 : ! **************************************************************************************************
1381 0 : FUNCTION casino_m_dependent_factor(l, m) RESULT(factor)
1382 : INTEGER, INTENT(IN) :: l, m
1383 : REAL(KIND=dp) :: factor
1384 :
1385 : INTEGER :: am
1386 : REAL(KIND=dp) :: prefactor
1387 :
1388 0 : am = ABS(m)
1389 0 : prefactor = MERGE(1.0_dp, 2.0_dp, am == 0)
1390 0 : factor = SQRT(prefactor*factorial(l - am)/factorial(l + am))
1391 0 : END FUNCTION casino_m_dependent_factor
1392 :
1393 : ! **************************************************************************************************
1394 : !> \brief Real factorial for small non-negative integers.
1395 : !> \param n ...
1396 : !> \return ...
1397 : ! **************************************************************************************************
1398 0 : FUNCTION factorial(n) RESULT(value)
1399 : INTEGER, INTENT(IN) :: n
1400 : REAL(KIND=dp) :: value
1401 :
1402 : INTEGER :: i
1403 :
1404 0 : value = 1.0_dp
1405 0 : DO i = 2, n
1406 0 : value = value*REAL(i, KIND=dp)
1407 : END DO
1408 0 : END FUNCTION factorial
1409 :
1410 : ! **************************************************************************************************
1411 : !> \brief Odd double factorial.
1412 : !> \param n ...
1413 : !> \return ...
1414 : ! **************************************************************************************************
1415 28 : FUNCTION odd_double_factorial(n) RESULT(value)
1416 : INTEGER, INTENT(IN) :: n
1417 : REAL(KIND=dp) :: value
1418 :
1419 : INTEGER :: i
1420 :
1421 28 : value = 1.0_dp
1422 28 : DO i = MAX(1, n), 1, -2
1423 28 : value = value*REAL(i, KIND=dp)
1424 : END DO
1425 28 : END FUNCTION odd_double_factorial
1426 :
1427 : ! **************************************************************************************************
1428 : !> \brief Computes the nuclear repulsion energy of a molecular system.
1429 : !> \param particle_set ...
1430 : !> \param kind_set ...
1431 : !> \param e_nn ...
1432 : ! **************************************************************************************************
1433 8 : SUBROUTINE nuclear_repulsion_energy(particle_set, kind_set, e_nn)
1434 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1435 : POINTER :: particle_set
1436 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1437 : POINTER :: kind_set
1438 : REAL(KIND=dp), INTENT(OUT) :: e_nn
1439 :
1440 : INTEGER :: i, ikind, j, jkind, natoms
1441 : REAL(KIND=dp) :: r_ij, zeff_i, zeff_j
1442 :
1443 8 : natoms = SIZE(particle_set)
1444 8 : e_nn = 0.0_dp
1445 22 : DO i = 1, natoms
1446 14 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
1447 14 : CALL get_qs_kind(kind_set(ikind), zeff=zeff_i)
1448 28 : DO j = i + 1, natoms
1449 24 : r_ij = NORM2(particle_set(i)%r - particle_set(j)%r)
1450 6 : CALL get_atomic_kind(particle_set(j)%atomic_kind, kind_number=jkind)
1451 6 : CALL get_qs_kind(kind_set(jkind), zeff=zeff_j)
1452 20 : e_nn = e_nn + zeff_i*zeff_j/r_ij
1453 : END DO
1454 : END DO
1455 8 : END SUBROUTINE nuclear_repulsion_energy
1456 :
1457 : ! **************************************************************************************************
1458 : !> \brief Computes the CASINO-compatible 3D periodic nuclear repulsion energy.
1459 : !> \param cell ...
1460 : !> \param periodicity ...
1461 : !> \param coord ...
1462 : !> \param charge ...
1463 : !> \param e_nn ...
1464 : ! **************************************************************************************************
1465 2 : SUBROUTINE periodic_nuclear_repulsion_energy(cell, periodicity, coord, charge, e_nn)
1466 : TYPE(cell_type), INTENT(IN), POINTER :: cell
1467 : INTEGER, INTENT(IN) :: periodicity
1468 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coord
1469 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charge
1470 : REAL(KIND=dp), INTENT(OUT) :: e_nn
1471 :
1472 : INTEGER :: gmax, i, ig1, ig2, ig3, j, n1, n2, n3, &
1473 : natoms, nmax
1474 : REAL(KIND=dp) :: alpha, alpha2, cutoff_arg, g_cut, g_sq, min_g, min_h, neut_energy, phase, &
1475 : r, real_cut, real_energy, recip_energy, self_energy, struc_im, struc_re, volume
1476 : REAL(KIND=dp), DIMENSION(3) :: delta, g_index, gvec, lattice_shift
1477 :
1478 2 : e_nn = 0.0_dp
1479 2 : IF (periodicity /= 3) RETURN
1480 :
1481 2 : volume = ABS(cell%deth)
1482 2 : IF (volume <= 0.0_dp) CPABORT("CASINO periodic nuclear repulsion requires a non-zero cell volume.")
1483 :
1484 2 : natoms = SIZE(charge)
1485 2 : IF (natoms == 0) RETURN
1486 :
1487 : min_h = HUGE(1.0_dp)
1488 : min_g = HUGE(1.0_dp)
1489 8 : DO i = 1, 3
1490 24 : min_h = MIN(min_h, NORM2(cell%hmat(:, i)))
1491 6 : g_index = 0.0_dp
1492 6 : g_index(i) = 1.0_dp
1493 30 : gvec = 2.0_dp*pi*MATMUL(TRANSPOSE(cell%h_inv), g_index)
1494 26 : min_g = MIN(min_g, NORM2(gvec))
1495 : END DO
1496 2 : IF (min_h <= 0.0_dp .OR. min_g <= 0.0_dp) THEN
1497 0 : CPABORT("CASINO periodic nuclear repulsion requires non-zero lattice vectors.")
1498 : END IF
1499 :
1500 2 : cutoff_arg = SQRT(-LOG(1.0E-12_dp))
1501 2 : alpha = SQRT(pi)*(REAL(natoms, KIND=dp)/volume)**(1.0_dp/3.0_dp)
1502 2 : alpha2 = alpha*alpha
1503 2 : real_cut = cutoff_arg/alpha
1504 2 : g_cut = 2.0_dp*alpha*cutoff_arg
1505 2 : nmax = MAX(1, CEILING(real_cut/min_h) + 1)
1506 2 : gmax = MAX(1, CEILING(g_cut/min_g) + 1)
1507 :
1508 2 : real_energy = 0.0_dp
1509 6 : DO i = 1, natoms
1510 14 : DO j = 1, natoms
1511 84 : DO n1 = -nmax, nmax
1512 728 : DO n2 = -nmax, nmax
1513 6552 : DO n3 = -nmax, nmax
1514 5832 : IF (i == j .AND. n1 == 0 .AND. n2 == 0 .AND. n3 == 0) CYCLE
1515 : lattice_shift = REAL(n1, KIND=dp)*cell%hmat(:, 1) + &
1516 : REAL(n2, KIND=dp)*cell%hmat(:, 2) + &
1517 23312 : REAL(n3, KIND=dp)*cell%hmat(:, 3)
1518 23312 : delta = coord(:, i) - coord(:, j) + lattice_shift
1519 23312 : r = NORM2(delta)
1520 6476 : IF (r <= real_cut) real_energy = real_energy + charge(i)*charge(j)*ERFC(alpha*r)/r
1521 : END DO
1522 : END DO
1523 : END DO
1524 : END DO
1525 : END DO
1526 2 : real_energy = 0.5_dp*real_energy
1527 :
1528 2 : recip_energy = 0.0_dp
1529 24 : DO ig1 = -gmax, gmax
1530 266 : DO ig2 = -gmax, gmax
1531 2926 : DO ig3 = -gmax, gmax
1532 2662 : IF (ig1 == 0 .AND. ig2 == 0 .AND. ig3 == 0) CYCLE
1533 10640 : g_index = [REAL(ig1, KIND=dp), REAL(ig2, KIND=dp), REAL(ig3, KIND=dp)]
1534 13300 : gvec = 2.0_dp*pi*MATMUL(TRANSPOSE(cell%h_inv), g_index)
1535 10640 : g_sq = DOT_PRODUCT(gvec, gvec)
1536 2660 : IF (SQRT(g_sq) > g_cut) CYCLE
1537 : struc_re = 0.0_dp
1538 : struc_im = 0.0_dp
1539 1212 : DO i = 1, natoms
1540 3232 : phase = DOT_PRODUCT(gvec, coord(:, i))
1541 808 : struc_re = struc_re + charge(i)*COS(phase)
1542 1212 : struc_im = struc_im + charge(i)*SIN(phase)
1543 : END DO
1544 : recip_energy = recip_energy + EXP(-g_sq/(4.0_dp*alpha2))/g_sq* &
1545 2904 : (struc_re*struc_re + struc_im*struc_im)
1546 : END DO
1547 : END DO
1548 : END DO
1549 2 : recip_energy = 2.0_dp*pi*recip_energy/volume
1550 :
1551 6 : self_energy = -alpha*SUM(charge*charge)/SQRT(pi)
1552 6 : neut_energy = -pi*SUM(charge)**2/(2.0_dp*alpha2*volume)
1553 2 : e_nn = real_energy + recip_energy + self_energy + neut_energy
1554 : END SUBROUTINE periodic_nuclear_repulsion_energy
1555 :
1556 : END MODULE casino_utils
|
