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