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