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