(git:97501a3)
Loading...
Searching...
No Matches
qs_chargemol.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Write wfx file, works as interface to chargemol and multiwfn
10!> \note Works only for HF and KS-type wavefunctions, execution is aborted otherwise
11!> \par History
12!> 06.2020 created [Hossam Elgabarty]
13!> 01.2021 modified [Hossam]
14!> \author Hossam Elgabarty
15! **************************************************************************************************
22 USE cell_types, ONLY: cell_type,&
25 USE cp_dbcsr_api, ONLY: dbcsr_p_type
27 USE cp_files, ONLY: open_file
28 USE cp_fm_types, ONLY: cp_fm_get_info,&
37 USE kinds, ONLY: default_path_length,&
38 dp
39 USE kpoint_types, ONLY: kpoint_type
40 USE machine, ONLY: m_timestamp,&
42 USE orbital_pointers, ONLY: nco,&
43 nso
44 USE orbital_symbols, ONLY: cgf_symbol,&
53 USE qs_kind_types, ONLY: get_qs_kind,&
56 USE qs_mo_types, ONLY: mo_set_type
60#include "./base/base_uses.f90"
61
62 IMPLICIT NONE
63 PRIVATE
64
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_chargemol'
66 LOGICAL, PARAMETER :: debug_this_module = .false.
67
68 INTEGER, PARAMETER :: chargemol_lmax = 5
69
70 PUBLIC :: write_wfx
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief ...
76!> \param qs_env ...
77!> \param dft_section ...
78! **************************************************************************************************
79 SUBROUTINE write_wfx(qs_env, dft_section)
80 TYPE(qs_environment_type), POINTER :: qs_env
81 TYPE(section_vals_type), POINTER :: dft_section
82
83 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_wfx'
84
85 CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
86 CHARACTER(len=2) :: asym
87 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
88 CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
89 CHARACTER(LEN=default_path_length) :: filename
90 CHARACTER(LEN=timestamp_length) :: timestamp
91 INTEGER :: handle, i, iatom, icgf, ico, ikind, iloop, imo, ipgf, irow, iset, isgf, ishell, &
92 ispin, iwfx, lshell, max_l, ncgf, ncol_global, nelec_alpha, nelec_beta, nkpoint, noccup, &
93 nrow_global, nset, nsgf, nspins, num_atoms, output_unit, roks_high, roks_low, tot_npgf
94 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomic_numbers, kind_of
95 INTEGER, DIMENSION(:), POINTER :: lmax, npgf, nshell
96 INTEGER, DIMENSION(:, :), POINTER :: l
97 LOGICAL :: do_kpoints, newl, periodic, unrestricted
98 REAL(kind=dp) :: zeta
99 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, cmatrix, smatrix
100 REAL(kind=dp), DIMENSION(:), POINTER :: core_charges
101 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
102 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
103 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
104 TYPE(cell_type), POINTER :: cell
105 TYPE(cp_logger_type), POINTER :: logger
106 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
107 TYPE(dft_control_type), POINTER :: dft_control
108 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
109 TYPE(kpoint_type), POINTER :: kpoint_env
110 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
111 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
112 TYPE(qs_energy_type), POINTER :: energy
113 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
114 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
115 TYPE(qs_scf_env_type), POINTER :: scf_env
116 TYPE(qs_subsys_type), POINTER :: subsys
117 TYPE(scf_control_type), POINTER :: scf_control
118 TYPE(section_vals_type), POINTER :: chargemol_section
119
120 ! !--------------------------------------------------------------------------------------------!
121
122 CALL timeset(routinen, handle)
123
124 logger => cp_get_default_logger()
125 output_unit = cp_logger_get_default_io_unit(logger)
126
127 IF (output_unit > 0) THEN
128 WRITE (output_unit, '(/,T2,A)') &
129 '!-----------------------------------------------------------------------------!'
130 WRITE (output_unit, '(T32,A)') "Writing Chargemol wfx file..."
131 WRITE (output_unit, '(T2,A)') &
132 '!-----------------------------------------------------------------------------!'
133 END IF
134
135 ! Collect environment info
136 cpassert(ASSOCIATED(qs_env))
137 CALL get_qs_env(qs_env, cell=cell, mos=mos, particle_set=particle_set, &
138 qs_kind_set=qs_kind_set, scf_env=scf_env, scf_control=scf_control, &
139 dft_control=dft_control, energy=energy, force=force, subsys=subsys, &
140 atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
141 kpoints=kpoint_env, matrix_s=matrix_s)
142
143 nkpoint = kpoint_env%nkp
144
145 IF (scf_control%use_ot) THEN
146 cpabort("OT is incompatible with .wfx output. Switch off OT.")
147 END IF
148
149 IF (dft_control%roks) THEN
150 CALL cp_abort(__location__, "The output of chargemol .wfx files is currently "// &
151 "implemented only for (un)restricted HF and Kohn-Sham-like methods.")
152 END IF
153
154 IF (dft_control%lsd) THEN
155 unrestricted = .true.
156 ELSE
157 unrestricted = .false.
158 END IF
159
160 chargemol_section => section_vals_get_subs_vals(dft_section, "PRINT%CHARGEMOL")
161
162 CALL section_vals_val_get(chargemol_section, "PERIODIC", l_val=periodic)
163
164 !!---------------------------------------------------------------------------------!
165 ! output unit
166 !!---------------------------------------------------------------------------------!
167
168 filename = cp_print_key_generate_filename(logger, chargemol_section, &
169 extension=".wfx", my_local=.false.)
170
171 CALL open_file(filename, file_status="UNKNOWN", &
172 file_action="WRITE", &
173 unit_number=iwfx)
174
175 !!---------------------------------------------------------------------------------!
176
177 IF (iwfx > 0) THEN
178
179 nspins = SIZE(mos)
180 IF (nspins == 1) THEN
181 nelec_alpha = mos(1)%nelectron/2
182 nelec_beta = nelec_alpha
183 ELSE
184 nelec_alpha = mos(1)%nelectron
185 nelec_beta = mos(2)%nelectron
186 END IF
187
188 IF (dft_control%roks) THEN
189 IF (mos(1)%homo > mos(2)%homo) THEN
190 roks_high = 1
191 roks_low = 2
192 ELSE
193 roks_high = 2
194 roks_low = 1
195 END IF
196 END IF
197
198 !!---------------------------------------------------------------------------------!
199 ! Initial parsing of topology and basis set, check maximum l .le. chargemol_lmax
200 !!---------------------------------------------------------------------------------!
201 num_atoms = SIZE(particle_set)
202 ALLOCATE (atoms_cart(3, num_atoms))
203 ALLOCATE (atom_symbols(num_atoms))
204 ALLOCATE (atomic_numbers(num_atoms))
205 ALLOCATE (core_charges(num_atoms))
206
207 max_l = 0
208 DO iatom = 1, num_atoms
209 NULLIFY (orb_basis_set)
210 atoms_cart(1:3, iatom) = particle_set(iatom)%r(1:3)
211 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
212 element_symbol=asym, kind_number=ikind)
213 atom_symbols(iatom) = asym
214 CALL get_ptable_info(asym, number=atomic_numbers(iatom))
215 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
216 core_charge=core_charges(iatom))
217 IF (ASSOCIATED(orb_basis_set)) THEN
218 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
219 nset=nset, lmax=lmax)
220
221 DO iset = 1, nset
222 max_l = max(max_l, lmax(iset))
223 END DO
224 ELSE
225 cpabort("Unknown basis set type. ")
226 END IF
227 IF (max_l > chargemol_lmax) THEN
228 CALL cp_abort(__location__, "Chargemol supports basis sets with l"// &
229 " up to 5 only (h angular functions).")
230 END IF
231 END DO
232
233 ! !===========================================================================!
234 ! Header comment and Title
235 ! !===========================================================================!
236 WRITE (iwfx, "(A)") "# Chargemol input file generated by CP2K "
237 CALL m_timestamp(timestamp)
238 WRITE (iwfx, "(A,/)") "# Creation date "//timestamp(:19)
239
240 WRITE (iwfx, "(A)") "<Title>"
241 WRITE (iwfx, *) trim(logger%iter_info%project_name)
242 WRITE (iwfx, "(A,/)") "</Title>"
243
244 ! !===========================================================================!
245 ! Keywords
246 ! TODO: check: always GTO??
247 ! !===========================================================================!
248 WRITE (iwfx, "(A)") "<Keywords>"
249 WRITE (iwfx, "(A)") "GTO"
250 WRITE (iwfx, "(A,/)") "</Keywords>"
251
252 ! !===========================================================================!
253 ! Model
254 ! !===========================================================================!
255 WRITE (iwfx, "(A)") "#<Model>"
256 IF (unrestricted) THEN
257 WRITE (iwfx, "(A)") "#Unrestricted Kohn-Sham"
258 ELSE
259 WRITE (iwfx, "(A)") "#Restricted Kohn-Sham"
260 END IF
261 WRITE (iwfx, "(A,/)") "#</Model>"
262
263 ! !===========================================================================!
264 ! Number of nuclei
265 ! !===========================================================================!
266 WRITE (iwfx, "(A)") "<Number of Nuclei>"
267 WRITE (iwfx, "(I6)") num_atoms
268 WRITE (iwfx, "(A,/)") "</Number of Nuclei>"
269
270 ! !===========================================================================!
271 ! Number of Occupied MOs
272 ! !===========================================================================!
273 WRITE (iwfx, "(A)") "<Number of Occupied Molecular Orbitals>"
274 noccup = 0
275 IF (dft_control%roks .AND. nspins == 2) THEN
276 noccup = max(mos(1)%homo, mos(2)%homo)
277 ELSE
278 DO ispin = 1, dft_control%nspins
279 noccup = noccup + mos(ispin)%homo
280 END DO
281 END IF
282 WRITE (iwfx, "(2I6)") noccup
283 WRITE (iwfx, "(A,/)") "</Number of Occupied Molecular Orbitals>"
284
285 ! !===========================================================================!
286 ! Number of Perturbations
287 ! !===========================================================================!
288 WRITE (iwfx, "(A)") "<Number of Perturbations>"
289 WRITE (iwfx, "(I6)") 0
290 WRITE (iwfx, "(A,/)") "</Number of Perturbations>"
291
292 ! !===========================================================================!
293 ! Total charge
294 ! !===========================================================================!
295 WRITE (iwfx, "(A)") "<Net Charge>"
296 WRITE (iwfx, "(F8.4)") real(dft_control%charge, kind=dp)
297 WRITE (iwfx, "(A,/)") "</Net Charge>"
298
299 ! !===========================================================================!
300 ! Number of electrons
301 ! !===========================================================================!
302 WRITE (iwfx, "(A)") "<Number of Electrons>"
303 WRITE (iwfx, "(I6)") scf_env%nelectron
304 WRITE (iwfx, "(A,/)") "</Number of Electrons>"
305
306 !===========================================================================!
307 ! Number of Alpha electrons
308 !===========================================================================!
309 WRITE (iwfx, "(A)") "<Number of Alpha Electrons>"
310 WRITE (iwfx, "(I6)") nelec_alpha
311 WRITE (iwfx, "(A,/)") "</Number of Alpha Electrons>"
312
313 !===========================================================================!
314 ! Number of Beta electrons
315 !===========================================================================!
316 WRITE (iwfx, "(A)") "<Number of Beta Electrons>"
317 WRITE (iwfx, "(I6)") nelec_beta
318 WRITE (iwfx, "(A,/)") "</Number of Beta Electrons>"
319
320 !===========================================================================!
321 ! Spin multiplicity
322 !===========================================================================!
323 WRITE (iwfx, "(A)") "<Electron Spin Multiplicity>"
324 WRITE (iwfx, "(I4)") dft_control%multiplicity
325 WRITE (iwfx, "(A,/)") "</Electron Spin Multiplicity>"
326
327 ! !===========================================================================!
328 ! Number of Core Electrons
329 ! !===========================================================================!
330 WRITE (iwfx, "(A)") "<Number of Core Electrons>"
331 WRITE (iwfx, "(F16.10)") sum(atomic_numbers) - sum(core_charges)
332 WRITE (iwfx, "(A,/)") "</Number of Core Electrons>"
333
334 ! !===========================================================================!
335 ! Nuclear Names
336 ! !===========================================================================!
337 WRITE (iwfx, "(A)") "<Nuclear Names>"
338 DO iatom = 1, num_atoms
339 WRITE (iwfx, "(A4)") atom_symbols(iatom)
340 END DO
341 WRITE (iwfx, "(A,/)") "</Nuclear Names>"
342
343 ! !===========================================================================!
344 ! Atomic Numbers
345 ! !===========================================================================!
346 WRITE (iwfx, "(A)") "<Atomic Numbers>"
347 DO iatom = 1, num_atoms
348 WRITE (iwfx, "(I4)") atomic_numbers(iatom)
349 END DO
350 WRITE (iwfx, "(A,/)") "</Atomic Numbers>"
351
352 ! !===========================================================================!
353 ! Nuclear charges
354 ! !===========================================================================!
355 WRITE (iwfx, "(A)") "<Nuclear Charges>"
356 DO iatom = 1, num_atoms
357 WRITE (iwfx, "(F12.8)") core_charges(iatom)
358 END DO
359 WRITE (iwfx, "(A,/)") "</Nuclear Charges>"
360
361 ! !===========================================================================!
362 ! Nuclear coordinates
363 ! !===========================================================================!
364 WRITE (iwfx, "(A)") "<Nuclear Cartesian Coordinates>"
365 DO iatom = 1, num_atoms
366 WRITE (iwfx, "(3ES26.16E3)") atoms_cart(1:3, iatom)
367 END DO
368 WRITE (iwfx, "(A,/)") "</Nuclear Cartesian Coordinates>"
369
370 ! !===========================================================================!
371 ! Periodic boundary conditions
372 ! !===========================================================================!
373 IF (periodic) THEN
374 CALL get_cell(cell)
375 IF (sum(cell%perd) == 0) THEN
376 CALL cp_warn(__location__, "Writing of periodic cell information"// &
377 " requested in input, but system is not periodic, ignoring request.")
378 ELSE
379 WRITE (iwfx, "(A)") "<Number of Translation Vectors>"
380 WRITE (iwfx, "(I3)") sum(cell%perd)
381 WRITE (iwfx, "(A,/)") "</Number of Translation Vectors>"
382
383 WRITE (iwfx, "(A)") "<Translation Vectors>"
384 DO iatom = 1, 3
385 IF (cell%perd(iatom) == 1) THEN
386 WRITE (iwfx, "(3F12.6)") cell%hmat(1:3, iatom)
387 END IF
388 END DO
389 WRITE (iwfx, "(A,/)") "</Translation Vectors>"
390 END IF
391 WRITE (iwfx, "(A)") "<Number of Kpoints>"
392 WRITE (iwfx, "(I3)") 1
393 WRITE (iwfx, "(A,/)") "</Number of Kpoints>"
394 WRITE (iwfx, "(A)") "<Kpoint Weights>"
395 WRITE (iwfx, "(I3)") 1
396 WRITE (iwfx, "(A,/)") "</Kpoint Weights>"
397 WRITE (iwfx, "(A)") "<Kpoint Fractional Coordinates>"
398 WRITE (iwfx, "(F8.6)") 0.0
399 WRITE (iwfx, "(A,/)") "</Kpoint Fractional Coordinates>"
400 END IF
401
402 ! !===========================================================================!
403 ! Primitive centers
404 ! !===========================================================================!
405 WRITE (iwfx, "(A)") "<Primitive Centers>"
406 tot_npgf = 0
407 DO iatom = 1, num_atoms
408 iloop = 1
409 NULLIFY (orb_basis_set)
410 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
411 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
412 IF (ASSOCIATED(orb_basis_set)) THEN
413 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
414 nset=nset, &
415 npgf=npgf, &
416 nshell=nshell, &
417 l=l, &
418 gcc=gcc)
419
420 DO iset = 1, nset
421 DO ishell = 1, nshell(iset)
422 lshell = l(ishell, iset)
423 DO ico = 1, nco(lshell)
424 DO ipgf = 1, npgf(iset)
425 tot_npgf = tot_npgf + 1
426 IF (mod(iloop + 10, 10) /= 0) THEN
427 WRITE (iwfx, "(I4)", advance="no") iatom
428 newl = .true.
429 ELSE
430 WRITE (iwfx, "(I4)") iatom
431 newl = .false.
432 END IF
433 iloop = iloop + 1
434 !END IF
435 END DO
436 END DO
437 END DO
438 END DO
439 IF (newl) THEN
440 WRITE (iwfx, *)
441 END IF
442 ELSE
443 cpabort("Unknown basis set type. ")
444 END IF
445 END DO
446 WRITE (iwfx, "(A,/)") "</Primitive Centers>"
447
448 ! !===========================================================================!
449 ! Number of primitives
450 ! !===========================================================================!
451 WRITE (iwfx, "(A)") "<Number of Primitives>"
452 WRITE (iwfx, "(I8)") tot_npgf
453 WRITE (iwfx, "(A,/)") "</Number of Primitives>"
454
455 ! !===========================================================================!
456 ! Primitive Types
457 ! !===========================================================================!
458 WRITE (iwfx, "(A)") "<Primitive Types>"
459 DO iatom = 1, num_atoms
460 iloop = 1
461 NULLIFY (orb_basis_set)
462 NULLIFY (bcgf_symbol)
463 NULLIFY (gcc)
464 NULLIFY (zet)
465 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
466 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
467 IF (ASSOCIATED(orb_basis_set)) THEN
468 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
469 nset=nset, &
470 nshell=nshell, &
471 l=l, &
472 npgf=npgf, &
473 ncgf=ncgf, &
474 zet=zet, &
475 cgf_symbol=bcgf_symbol, &
476 sgf_symbol=bsgf_symbol, &
477 gcc=gcc)
478
479 icgf = 1
480 DO iset = 1, nset
481 DO ishell = 1, nshell(iset)
482 lshell = l(ishell, iset)
483 DO ico = 1, nco(lshell)
484 DO ipgf = 1, orb_basis_set%npgf(iset)
485 IF (mod(iloop + 10, 10) /= 0) THEN
486 WRITE (iwfx, '(I6)', advance="no") &
487 pgf_type_chargemol(bcgf_symbol(icgf) (3:))
488 newl = .true.
489 ELSE
490 WRITE (iwfx, '(I6)') &
491 pgf_type_chargemol(bcgf_symbol(icgf) (3:))
492 newl = .false.
493 END IF
494 iloop = iloop + 1
495 !END IF
496 END DO
497 icgf = icgf + 1
498 END DO !ico
499 END DO ! ishell
500 END DO ! iset
501
502 ELSE
503 cpabort("Unknown basis set type.")
504 END IF
505 IF (newl) THEN
506 WRITE (iwfx, *)
507 END IF
508 END DO ! iatom
509 WRITE (iwfx, "(A,/)") "</Primitive Types>"
510
511 ! !===========================================================================!
512 ! Primitive Exponents
513 ! !===========================================================================!
514 !! NOTE: CP2K orders the atomic orbitals as follows:
515 !! Cartesian orbitals: x (slowest loop), y, z (fastest loop)
516 !! Spherical orbitals: m = -l, ..., 0, ..., +l
517
518 WRITE (iwfx, "(A)") "<Primitive Exponents>"
519 DO iatom = 1, num_atoms
520 NULLIFY (orb_basis_set)
521 NULLIFY (gcc)
522 NULLIFY (zet)
523 iloop = 1
524
525 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
526 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
527 IF (ASSOCIATED(orb_basis_set)) THEN
528 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
529 nset=nset, &
530 npgf=npgf, &
531 nshell=nshell, &
532 l=l, &
533 zet=zet, &
534 gcc=gcc)
535
536 DO iset = 1, nset
537 DO ishell = 1, nshell(iset)
538 lshell = l(ishell, iset)
539 DO ico = 1, nco(lshell)
540 DO ipgf = 1, npgf(iset)
541 IF (mod(iloop + 5, 5) /= 0) THEN
542 WRITE (iwfx, "(ES20.10)", advance="no") zet(ipgf, iset)
543 newl = .true.
544 ELSE
545 WRITE (iwfx, "(ES20.10)") zet(ipgf, iset)
546 newl = .false.
547 END IF
548 iloop = iloop + 1
549 !END IF
550 END DO
551 END DO
552 END DO
553 END DO
554 IF (newl) THEN
555 WRITE (iwfx, *)
556 END IF
557 ELSE
558 cpabort("Unknown basis set type. ")
559 END IF
560 END DO
561 WRITE (iwfx, "(A,/)") "</Primitive Exponents>"
562
563 ! !===========================================================================!
564 ! MO occupation numbers
565 ! nonesense for roks at the moment
566 ! !===========================================================================!
567 WRITE (iwfx, "(A)") "<Molecular Orbital Occupation Numbers>"
568 IF (.NOT. dft_control%roks) THEN
569 DO ispin = 1, nspins
570 DO imo = 1, mos(ispin)%homo
571 WRITE (iwfx, "(f8.6)") &
572 mos(ispin)%occupation_numbers(imo)
573 END DO
574 END DO
575 ELSE
576 DO imo = 1, mos(roks_low)%homo
577 WRITE (iwfx, "(*(f8.6,1x))") &
578 mos(roks_low)%occupation_numbers(imo) &
579 + mos(roks_low)%occupation_numbers(imo)
580 END DO
581 DO imo = mos(roks_low)%homo + 1, mos(roks_high)%homo
582 WRITE (iwfx, "(f8.6)") &
583 mos(roks_high)%occupation_numbers(imo)
584 END DO
585 END IF
586 WRITE (iwfx, "(A,/)") "</Molecular Orbital Occupation Numbers>"
587
588 ! !===========================================================================!
589 ! MO energies
590 ! !===========================================================================!
591 WRITE (iwfx, "(A)") "<Molecular Orbital Energies>"
592 DO ispin = 1, nspins
593 DO imo = 1, mos(ispin)%homo
594 WRITE (iwfx, "(ES20.10)") mos(ispin)%eigenvalues(imo)
595 END DO
596 END DO
597 WRITE (iwfx, "(A,/)") "</Molecular Orbital Energies>"
598
599 ! !===========================================================================!
600 ! MO Spin types
601 ! nonesense for ROKS
602 ! !===========================================================================!
603 WRITE (iwfx, "(A)") "<Molecular Orbital Spin Types>"
604 DO imo = 1, mos(1)%homo
605 IF (nspins == 1) THEN
606 WRITE (iwfx, '(A15)') "Alpha and Beta"
607 ELSEIF (dft_control%uks) THEN
608 WRITE (iwfx, '(A6)') "Alpha"
609 ELSE
610 CALL cp_abort(__location__, "This wavefunction type is currently"// &
611 " not supported by the chargemol interface.")
612 END IF
613 END DO
614 IF (nspins == 2) THEN
615 DO imo = 1, mos(2)%homo
616 WRITE (iwfx, '(A5)') "Beta"
617 END DO
618 END IF
619 WRITE (iwfx, "(A,/)") "</Molecular Orbital Spin Types>"
620
621 ! !===========================================================================!
622 ! Kpoint index of orbitals
623 ! !===========================================================================!
624 IF (periodic) THEN
625 WRITE (iwfx, "(A)") "<Kpoint Index for Orbitals>"
626 DO ispin = 1, nspins
627 DO imo = 1, mos(ispin)%homo
628 WRITE (iwfx, "(I6)") 1
629 END DO
630 END DO
631 WRITE (iwfx, "(A,/)") "</Kpoint Index for Orbitals>"
632 END IF
633
634 ! !===========================================================================!
635 ! MO Primitive Coefficients
636 ! !===========================================================================!
637 WRITE (iwfx, "(A)") "<Molecular Orbital Primitive Coefficients>"
638 NULLIFY (orb_basis_set)
639 NULLIFY (gcc)
640
641 IF (mos(1)%use_mo_coeff_b) THEN
642 DO ispin = 1, SIZE(mos)
643 IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
644 cpassert(.false.)
645 END IF
646 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
647 mos(ispin)%mo_coeff)
648 END DO
649 END IF
650
651 DO ispin = 1, nspins
652 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
653 nrow_global=nrow_global, &
654 ncol_global=ncol_global)
655
656 ALLOCATE (smatrix(nrow_global, ncol_global))
657 smatrix = 0.0_dp
658
659 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
660
661 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
662
663 ALLOCATE (cmatrix(ncgf, ncol_global))
664
665 cmatrix = 0.0_dp
666
667 ! get MO coefficients in terms of contracted cartesian functions
668 icgf = 1
669 isgf = 1
670 DO iatom = 1, num_atoms
671 NULLIFY (orb_basis_set)
672 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
673 CALL get_qs_kind(qs_kind_set(ikind), &
674 basis_set=orb_basis_set)
675 IF (ASSOCIATED(orb_basis_set)) THEN
676 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
677 nset=nset, &
678 nshell=nshell, &
679 l=l)
680 DO iset = 1, nset
681 DO ishell = 1, nshell(iset)
682 lshell = l(ishell, iset)
683 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, &
684 nso(lshell), 1.0_dp, &
685 orbtramat(lshell)%c2s, nso(lshell), &
686 smatrix(isgf, 1), nsgf, 0.0_dp, &
687 cmatrix(icgf, 1), ncgf)
688 !IF (iatom==1 .and. debug_this_module) THEN
689 ! print *, lshell
690 ! print *, "size ", size(orbtramat(lshell)%s2c,1),",",size(orbtramat(lshell)%s2c,2)
691 ! DO itest = 1, size(orbtramat(lshell)%s2c,1)
692 ! print *, orbtramat(lshell)%s2c(itest,:)
693 ! END DO
694 ! print *, " "
695 ! DO itest = 1, 5
696 ! print *, my_d_orbtramat(itest,:)
697 ! END DO
698 !END IF
699 icgf = icgf + nco(lshell)
700 isgf = isgf + nso(lshell)
701 END DO
702 END DO
703 ELSE
704 cpabort("Unknown basis set type. ")
705 END IF
706 END DO ! iatom
707
708 ! Now write MO coefficients in terms of cartesian primitives
709 DO imo = 1, mos(ispin)%homo
710
711 WRITE (iwfx, "(A)") "<MO Number>"
712 IF (nspins == 1 .OR. ispin == 1) THEN
713 WRITE (iwfx, "(I6)") imo
714 ELSE
715 WRITE (iwfx, "(I6)") imo + mos(1)%homo
716 END IF
717 WRITE (iwfx, "(A)") "</MO Number>"
718
719 irow = 1 ! row of cmatrix, each row corresponds to
720 ! a contracted Cartesian function
721 DO iatom = 1, num_atoms
722 iloop = 1
723 NULLIFY (orb_basis_set)
724 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
725 kind_number=ikind)
726 CALL get_qs_kind(qs_kind_set(ikind), &
727 basis_set=orb_basis_set)
728 IF (ASSOCIATED(orb_basis_set)) THEN
729 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
730 nset=nset, &
731 nshell=nshell, &
732 gcc=gcc, &
733 l=l)
734
735 icgf = 1
736 DO iset = 1, nset
737 DO ishell = 1, nshell(iset)
738 lshell = orb_basis_set%l(ishell, iset)
739
740 DO ico = orb_basis_set%first_cgf(ishell, iset), &
741 orb_basis_set%last_cgf(ishell, iset)
742
743 !loop over primitives
744 DO ipgf = 1, orb_basis_set%npgf(iset)
745
746 zeta = orb_basis_set%zet(ipgf, iset)
747
748 IF (mod(iloop + 5, 5) /= 0) THEN
749 WRITE (iwfx, '(ES20.10)', advance="no") &
750 cmatrix(irow, imo)* &
751 orb_basis_set%norm_cgf(ico)* &
752 orb_basis_set%gcc(ipgf, ishell, iset)
753 newl = .true.
754 ELSE
755 WRITE (iwfx, '(ES20.10)') &
756 cmatrix(irow, imo)* &
757 orb_basis_set%norm_cgf(ico)* &
758 orb_basis_set%gcc(ipgf, ishell, iset)
759 newl = .false.
760 END IF
761 iloop = iloop + 1
762 END DO ! end loop over primitives
763 irow = irow + 1
764 END DO !ico
765 END DO ! ishell
766 END DO ! iset
767
768 ELSE
769 cpabort("Unknown basis set type.")
770 END IF
771 IF (newl) THEN
772 WRITE (iwfx, *)
773 END IF
774 END DO ! iatom
775 !Write (iwfx,*)
776 END DO ! imo
777
778 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
779 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
780
781 END DO ! ispin
782 WRITE (iwfx, "(A,/)") "</Molecular Orbital Primitive Coefficients>"
783
784 ! !===========================================================================!
785 ! Energy
786 ! !===========================================================================!
787 WRITE (iwfx, "(A)") "<Energy>"
788 WRITE (iwfx, "(ES26.16E3)") energy%total
789 WRITE (iwfx, "(A,/)") "</Energy>"
790
791 ! !===========================================================================!
792 ! Virial ratio
793 ! !===========================================================================!
794 WRITE (iwfx, "(A)") "<Virial Ratio (-V/T)>"
795 WRITE (iwfx, "(ES20.12)") - 1*(energy%total - energy%kinetic)/energy%kinetic
796 WRITE (iwfx, "(A,/)") "</Virial Ratio (-V/T)>"
797
798 ! !===========================================================================!
799 ! Force-related quantities
800 ! inactivated for now
801 ! !===========================================================================!
802 IF (ASSOCIATED(force) .AND. debug_this_module) THEN
803 WRITE (iwfx, "(A)") "<Nuclear Cartesian Energy Gradients>"
804
805 CALL get_qs_env(qs_env, force=force, atomic_kind_set=atomic_kind_set)
806 ALLOCATE (atom_of_kind(num_atoms))
807 ALLOCATE (kind_of(num_atoms))
808 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
809 atom_of_kind=atom_of_kind, &
810 kind_of=kind_of)
811
812 DO iatom = 1, num_atoms
813 ikind = kind_of(iatom)
814 i = atom_of_kind(iatom)
815 WRITE (iwfx, "(3ES20.12E3)") force(ikind)%total(1:3, i)
816 END DO
817
818 WRITE (iwfx, "(A,/)") "<Nuclear Cartesian Energy Gradients>"
819
820 ! W
821 WRITE (iwfx, "(A)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
822 WRITE (iwfx, "(A)") ""
823 WRITE (iwfx, "(A,/)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
824
825 ! Full Virial ratio
826 WRITE (iwfx, "(A)") "<Full Virial Ratio, -(V - W)/T)>"
827 WRITE (iwfx, "(A)") ""
828 WRITE (iwfx, "(A,/)") "</Full Virial Ratio, -(V - W)/T)>"
829
830 DEALLOCATE (atom_of_kind)
831 DEALLOCATE (kind_of)
832
833 END IF
834
835 ! !===========================================================================!
836 ! Core-density
837 ! !===========================================================================!
838 ! Additional Electron Density Function (EDF)
839 ! WRITE (iwfx, "(A)") "#<Additional Electron Density Function (EDF)>"
840 ! WRITE (iwfx, "(A,/)") "#</Additional Electron Density Function (EDF)>"
841
842 ! !-----------------------------------!
843 ! Done
844 ! !-----------------------------------!
845 DEALLOCATE (atoms_cart)
846 DEALLOCATE (atom_symbols)
847 DEALLOCATE (atomic_numbers)
848 DEALLOCATE (core_charges)
849
850 ELSE ! no output unit
851 iwfx = -1
852 END IF
853
854 IF (output_unit > 0) THEN
855 WRITE (output_unit, '(/,T2,A)') &
856 '!--------------------------- Chargemol wfx written ---------------------------!'
857 END IF
858
859 CALL timestop(handle)
860
861 END SUBROUTINE write_wfx
862
863! **************************************************************************************************
864!> \brief ...
865!> \param l_symb ...
866!> \return ...
867! **************************************************************************************************
868 INTEGER FUNCTION pgf_type_chargemol(l_symb) RESULT(label)
869 !INTEGER, INTENT(OUT) :: label
870 CHARACTER(len=*), INTENT(IN) :: l_symb
871
872 SELECT CASE (l_symb)
873 CASE ("s")
874 label = 1
875 CASE ("px")
876 label = 2
877 CASE ("py")
878 label = 3
879 CASE ("pz")
880 label = 4
881 CASE ("dx2")
882 label = 5
883 CASE ("dy2")
884 label = 6
885 CASE ("dz2")
886 label = 7
887 CASE ("dxy")
888 label = 8
889 CASE ("dxz")
890 label = 9
891 CASE ("dyz")
892 label = 10
893 CASE ("fx3")
894 label = 11
895 CASE ("fy3")
896 label = 12
897 CASE ("fz3")
898 label = 13
899 CASE ("fx2y")
900 label = 14
901 CASE ("fx2z")
902 label = 15
903 CASE ("fy2z")
904 label = 16
905 CASE ("fxy2")
906 label = 17
907 CASE ("fxz2")
908 label = 18
909 CASE ("fyz2")
910 label = 19
911 CASE ("fxyz")
912 label = 20
913 CASE ("gx4")
914 label = 21
915 CASE ("gy4")
916 label = 22
917 CASE ("gz4")
918 label = 23
919 CASE ("gx3y")
920 label = 24
921 CASE ("gx3z")
922 label = 25
923 CASE ("gxy3")
924 label = 26
925 CASE ("gy3z")
926 label = 27
927 CASE ("gxz3")
928 label = 28
929 CASE ("gyz3")
930 label = 29
931 CASE ("gx2y2")
932 label = 30
933 CASE ("gx2z2")
934 label = 31
935 CASE ("gy2z2")
936 label = 32
937 CASE ("gx2yz")
938 label = 33
939 CASE ("gxy2z")
940 label = 34
941 CASE ("gxyz2")
942 label = 35
943 CASE ("hz5")
944 label = 36
945 CASE ("hyz4")
946 label = 37
947 CASE ("hy2z3")
948 label = 38
949 CASE ("hy3z2")
950 label = 39
951 CASE ("hy4z")
952 label = 40
953 CASE ("hy5")
954 label = 41
955 CASE ("hxz4")
956 label = 42
957 CASE ("hxyz3")
958 label = 43
959 CASE ("hxy2z2")
960 label = 44
961 CASE ("hxy3z")
962 label = 45
963 CASE ("hxy4")
964 label = 46
965 CASE ("hx2z3")
966 label = 47
967 CASE ("hx2yz2")
968 label = 48
969 CASE ("hx2y2z")
970 label = 49
971 CASE ("hx2y3")
972 label = 50
973 CASE ("hx3z2")
974 label = 51
975 CASE ("hx3yz")
976 label = 52
977 CASE ("hx3y2")
978 label = 53
979 CASE ("hx4z")
980 label = 54
981 CASE ("hx4y")
982 label = 55
983 CASE ("hx5")
984 label = 56
985 CASE default
986 cpabort("The chargemol interface supports basis functions up to l=5 only")
987 !label = 99
988 END SELECT
989
990 END FUNCTION pgf_type_chargemol
991
992END MODULE qs_chargemol
993
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_path_length
Definition kinds.F:58
Types and basic routines needed for a kpoint calculation.
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
integer, parameter, public timestamp_length
Definition machine.F:53
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:380
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public nso
orbital_symbols
character(len=12) function, public cgf_symbol(n, lxyz)
Build a Cartesian orbital symbol (orbital labels for printing).
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
type(orbtramat_type), dimension(:), pointer, public orbtramat
Define the data structure for the particle information.
Periodic Table related data definitions.
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Write wfx file, works as interface to chargemol and multiwfn.
subroutine, public write_wfx(qs_env, dft_section)
...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
module that contains the definitions of the scf types
types that represent a quickstep subsys
parameters that control an scf iteration
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
Provides all information about a quickstep kind.