(git:374b731)
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-2024 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,&
26 USE cp_files, ONLY: open_file
27 USE cp_fm_types, ONLY: cp_fm_get_info,&
33 USE dbcsr_api, ONLY: dbcsr_p_type
37 USE kinds, ONLY: default_string_length,&
38 dp
39 USE kpoint_types, ONLY: kpoint_type
40 USE machine, ONLY: m_datum
41 USE orbital_pointers, ONLY: nco,&
42 nso
43 USE orbital_symbols, ONLY: cgf_symbol,&
52 USE qs_kind_types, ONLY: get_qs_kind,&
55 USE qs_mo_types, ONLY: mo_set_type
59#include "./base/base_uses.f90"
60
61 IMPLICIT NONE
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_chargemol'
65 LOGICAL, PARAMETER :: debug_this_module = .false.
66
67 INTEGER, PARAMETER :: chargemol_lmax = 5
68
69 PUBLIC :: write_wfx
70
71CONTAINS
72
73! **************************************************************************************************
74!> \brief ...
75!> \param qs_env ...
76!> \param dft_section ...
77! **************************************************************************************************
78 SUBROUTINE write_wfx(qs_env, dft_section)
79 TYPE(qs_environment_type), POINTER :: qs_env
80 TYPE(section_vals_type), POINTER :: dft_section
81
82 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_wfx'
83
84 CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
85 CHARACTER(len=2) :: asym
86 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
87 CHARACTER(LEN=256) :: datx
88 CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
89 CHARACTER(len=default_string_length) :: filename
90 INTEGER :: handle, i, iatom, icgf, ico, ikind, iloop, imo, ipgf, irow, iset, isgf, ishell, &
91 ispin, iwfx, lshell, max_l, ncgf, ncol_global, nelec_alpha, nelec_beta, nkpoint, noccup, &
92 nrow_global, nset, nsgf, nspins, num_atoms, output_unit, roks_high, roks_low, tot_npgf
93 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomic_numbers, kind_of
94 INTEGER, DIMENSION(:), POINTER :: lmax, npgf, nshell
95 INTEGER, DIMENSION(:, :), POINTER :: l
96 LOGICAL :: do_kpoints, newl, periodic, unrestricted
97 REAL(kind=dp) :: zeta
98 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, cmatrix, smatrix
99 REAL(kind=dp), DIMENSION(:), POINTER :: core_charges
100 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
101 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
102 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
103 TYPE(cell_type), POINTER :: cell
104 TYPE(cp_logger_type), POINTER :: logger
105 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
106 TYPE(dft_control_type), POINTER :: dft_control
107 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
108 TYPE(kpoint_type), POINTER :: kpoint_env
109 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
110 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
111 TYPE(qs_energy_type), POINTER :: energy
112 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
113 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114 TYPE(qs_scf_env_type), POINTER :: scf_env
115 TYPE(qs_subsys_type), POINTER :: subsys
116 TYPE(scf_control_type), POINTER :: scf_control
117 TYPE(section_vals_type), POINTER :: chargemol_section
118
119 ! !--------------------------------------------------------------------------------------------!
120
121 CALL timeset(routinen, handle)
122
123 logger => cp_get_default_logger()
124 output_unit = cp_logger_get_default_io_unit(logger)
125
126 IF (output_unit > 0) THEN
127 WRITE (output_unit, '(/,T2,A)') &
128 '!-----------------------------------------------------------------------------!'
129 WRITE (output_unit, '(T32,A)') "Writing Chargemol wfx file..."
130 WRITE (output_unit, '(T2,A)') &
131 '!-----------------------------------------------------------------------------!'
132 END IF
133
134 ! Collect environment info
135 cpassert(ASSOCIATED(qs_env))
136 CALL get_qs_env(qs_env, cell=cell, mos=mos, particle_set=particle_set, &
137 qs_kind_set=qs_kind_set, scf_env=scf_env, scf_control=scf_control, &
138 dft_control=dft_control, energy=energy, force=force, subsys=subsys, &
139 atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
140 kpoints=kpoint_env, matrix_s=matrix_s)
141
142 nkpoint = kpoint_env%nkp
143
144 IF (scf_control%use_ot) THEN
145 cpabort("OT is incompatible with .wfx output. Switch off OT.")
146 END IF
147
148 IF (dft_control%do_tddfpt_calculation .OR. dft_control%roks) THEN
149 CALL cp_abort(__location__, "The output of chargemol .wfx files is currently "// &
150 "implemented only for (un)restricted HF and Kohn-Sham-like methods.")
151 END IF
152
153 IF (dft_control%lsd) THEN
154 unrestricted = .true.
155 ELSE
156 unrestricted = .false.
157 END IF
158
159 chargemol_section => section_vals_get_subs_vals(dft_section, "PRINT%CHARGEMOL")
160
161 CALL section_vals_val_get(chargemol_section, "PERIODIC", l_val=periodic)
162
163 !!---------------------------------------------------------------------------------!
164 ! output unit
165 !!---------------------------------------------------------------------------------!
166
167 filename = cp_print_key_generate_filename(logger, chargemol_section, &
168 extension=".wfx", my_local=.false.)
169
170 CALL open_file(filename, file_status="UNKNOWN", &
171 file_action="WRITE", &
172 unit_number=iwfx)
173
174 !!---------------------------------------------------------------------------------!
175
176 IF (iwfx > 0) THEN
177
178 nspins = SIZE(mos)
179 IF (nspins == 1) THEN
180 nelec_alpha = mos(1)%nelectron/2
181 nelec_beta = nelec_alpha
182 ELSE
183 nelec_alpha = mos(1)%nelectron
184 nelec_beta = mos(2)%nelectron
185 END IF
186
187 IF (dft_control%roks) THEN
188 IF (mos(1)%homo > mos(2)%homo) THEN
189 roks_high = 1
190 roks_low = 2
191 ELSE
192 roks_high = 2
193 roks_low = 1
194 END IF
195 END IF
196
197 !!---------------------------------------------------------------------------------!
198 ! Initial parsing of topology and basis set, check maximum l .le. chargemol_lmax
199 !!---------------------------------------------------------------------------------!
200 num_atoms = SIZE(particle_set)
201 ALLOCATE (atoms_cart(3, num_atoms))
202 ALLOCATE (atom_symbols(num_atoms))
203 ALLOCATE (atomic_numbers(num_atoms))
204 ALLOCATE (core_charges(num_atoms))
205
206 max_l = 0
207 DO iatom = 1, num_atoms
208 NULLIFY (orb_basis_set)
209 atoms_cart(1:3, iatom) = particle_set(iatom)%r(1:3)
210 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
211 element_symbol=asym, kind_number=ikind)
212 atom_symbols(iatom) = asym
213 CALL get_ptable_info(asym, number=atomic_numbers(iatom))
214 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
215 core_charge=core_charges(iatom))
216 IF (ASSOCIATED(orb_basis_set)) THEN
217 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
218 nset=nset, lmax=lmax)
219
220 DO iset = 1, nset
221 max_l = max(max_l, lmax(iset))
222 END DO
223 ELSE
224 cpabort("Unknown basis set type. ")
225 END IF
226 IF (max_l > chargemol_lmax) THEN
227 CALL cp_abort(__location__, "Chargemol supports basis sets with l"// &
228 " up to 5 only (h angular functions).")
229 END IF
230 END DO
231
232 ! !===========================================================================!
233 ! Header comment and Title
234 ! !===========================================================================!
235 WRITE (iwfx, "(A)") "# Chargemol input file generated by CP2K "
236 CALL m_datum(datx)
237 WRITE (iwfx, "(A,/)") "# Creation date "//trim(datx)
238
239 WRITE (iwfx, "(A)") "<Title>"
240 WRITE (iwfx, *) trim(logger%iter_info%project_name)
241 WRITE (iwfx, "(A,/)") "</Title>"
242
243 ! !===========================================================================!
244 ! Keywords
245 ! TODO: check: always GTO??
246 ! !===========================================================================!
247 WRITE (iwfx, "(A)") "<Keywords>"
248 WRITE (iwfx, "(A)") "GTO"
249 WRITE (iwfx, "(A,/)") "</Keywords>"
250
251 ! !===========================================================================!
252 ! Model
253 ! !===========================================================================!
254 WRITE (iwfx, "(A)") "#<Model>"
255 IF (unrestricted) THEN
256 WRITE (iwfx, "(A)") "#Unrestricted Kohn-Sham"
257 ELSE
258 WRITE (iwfx, "(A)") "#Restricted Kohn-Sham"
259 END IF
260 WRITE (iwfx, "(A,/)") "#</Model>"
261
262 ! !===========================================================================!
263 ! Number of nuclei
264 ! !===========================================================================!
265 WRITE (iwfx, "(A)") "<Number of Nuclei>"
266 WRITE (iwfx, "(I6)") num_atoms
267 WRITE (iwfx, "(A,/)") "</Number of Nuclei>"
268
269 ! !===========================================================================!
270 ! Number of Occupied MOs
271 ! !===========================================================================!
272 WRITE (iwfx, "(A)") "<Number of Occupied Molecular Orbitals>"
273 noccup = 0
274 IF (dft_control%roks .AND. nspins == 2) THEN
275 noccup = max(mos(1)%homo, mos(2)%homo)
276 ELSE
277 DO ispin = 1, dft_control%nspins
278 noccup = noccup + mos(ispin)%homo
279 END DO
280 END IF
281 WRITE (iwfx, "(2I6)") noccup
282 WRITE (iwfx, "(A,/)") "</Number of Occupied Molecular Orbitals>"
283
284 ! !===========================================================================!
285 ! Number of Perturbations
286 ! !===========================================================================!
287 WRITE (iwfx, "(A)") "<Number of Perturbations>"
288 WRITE (iwfx, "(I6)") 0
289 WRITE (iwfx, "(A,/)") "</Number of Perturbations>"
290
291 ! !===========================================================================!
292 ! Total charge
293 ! !===========================================================================!
294 WRITE (iwfx, "(A)") "<Net Charge>"
295 WRITE (iwfx, "(F8.4)") real(dft_control%charge, kind=dp)
296 WRITE (iwfx, "(A,/)") "</Net Charge>"
297
298 ! !===========================================================================!
299 ! Number of electrons
300 ! !===========================================================================!
301 WRITE (iwfx, "(A)") "<Number of Electrons>"
302 WRITE (iwfx, "(I6)") scf_env%nelectron
303 WRITE (iwfx, "(A,/)") "</Number of Electrons>"
304
305 !===========================================================================!
306 ! Number of Alpha electrons
307 !===========================================================================!
308 WRITE (iwfx, "(A)") "<Number of Alpha Electrons>"
309 WRITE (iwfx, "(I6)") nelec_alpha
310 WRITE (iwfx, "(A,/)") "</Number of Alpha Electrons>"
311
312 !===========================================================================!
313 ! Number of Beta electrons
314 !===========================================================================!
315 WRITE (iwfx, "(A)") "<Number of Beta Electrons>"
316 WRITE (iwfx, "(I6)") nelec_beta
317 WRITE (iwfx, "(A,/)") "</Number of Beta Electrons>"
318
319 !===========================================================================!
320 ! Spin multiplicity
321 !===========================================================================!
322 WRITE (iwfx, "(A)") "<Electron Spin Multiplicity>"
323 WRITE (iwfx, "(I4)") dft_control%multiplicity
324 WRITE (iwfx, "(A,/)") "</Electron Spin Multiplicity>"
325
326 ! !===========================================================================!
327 ! Number of Core Electrons
328 ! !===========================================================================!
329 WRITE (iwfx, "(A)") "<Number of Core Electrons>"
330 WRITE (iwfx, "(F16.10)") sum(atomic_numbers) - sum(core_charges)
331 WRITE (iwfx, "(A,/)") "</Number of Core Electrons>"
332
333 ! !===========================================================================!
334 ! Nuclear Names
335 ! !===========================================================================!
336 WRITE (iwfx, "(A)") "<Nuclear Names>"
337 DO iatom = 1, num_atoms
338 WRITE (iwfx, "(A4)") atom_symbols(iatom)
339 END DO
340 WRITE (iwfx, "(A,/)") "</Nuclear Names>"
341
342 ! !===========================================================================!
343 ! Atomic Numbers
344 ! !===========================================================================!
345 WRITE (iwfx, "(A)") "<Atomic Numbers>"
346 DO iatom = 1, num_atoms
347 WRITE (iwfx, "(I4)") atomic_numbers(iatom)
348 END DO
349 WRITE (iwfx, "(A,/)") "</Atomic Numbers>"
350
351 ! !===========================================================================!
352 ! Nuclear charges
353 ! !===========================================================================!
354 WRITE (iwfx, "(A)") "<Nuclear Charges>"
355 DO iatom = 1, num_atoms
356 WRITE (iwfx, "(F12.8)") core_charges(iatom)
357 END DO
358 WRITE (iwfx, "(A,/)") "</Nuclear Charges>"
359
360 ! !===========================================================================!
361 ! Nuclear coordinates
362 ! !===========================================================================!
363 WRITE (iwfx, "(A)") "<Nuclear Cartesian Coordinates>"
364 DO iatom = 1, num_atoms
365 WRITE (iwfx, "(3ES26.16E3)") atoms_cart(1:3, iatom)
366 END DO
367 WRITE (iwfx, "(A,/)") "</Nuclear Cartesian Coordinates>"
368
369 ! !===========================================================================!
370 ! Periodic boundary conditions
371 ! !===========================================================================!
372 IF (periodic) THEN
373 CALL get_cell(cell)
374 IF (sum(cell%perd) == 0) THEN
375 CALL cp_warn(__location__, "Writing of periodic cell information"// &
376 " requested in input, but system is not periodic, ignoring request.")
377 ELSE
378 WRITE (iwfx, "(A)") "<Number of Translation Vectors>"
379 WRITE (iwfx, "(I3)") sum(cell%perd)
380 WRITE (iwfx, "(A,/)") "</Number of Translation Vectors>"
381
382 WRITE (iwfx, "(A)") "<Translation Vectors>"
383 DO iatom = 1, 3
384 IF (cell%perd(iatom) == 1) THEN
385 WRITE (iwfx, "(3F12.6)") cell%hmat(1:3, iatom)
386 END IF
387 END DO
388 WRITE (iwfx, "(A,/)") "</Translation Vectors>"
389 END IF
390 WRITE (iwfx, "(A)") "<Number of Kpoints>"
391 WRITE (iwfx, "(I3)") 1
392 WRITE (iwfx, "(A,/)") "</Number of Kpoints>"
393 WRITE (iwfx, "(A)") "<Kpoint Weights>"
394 WRITE (iwfx, "(I3)") 1
395 WRITE (iwfx, "(A,/)") "</Kpoint Weights>"
396 WRITE (iwfx, "(A)") "<Kpoint Fractional Coordinates>"
397 WRITE (iwfx, "(F8.6)") 0.0
398 WRITE (iwfx, "(A,/)") "</Kpoint Fractional Coordinates>"
399 END IF
400
401 ! !===========================================================================!
402 ! Primitive centers
403 ! !===========================================================================!
404 WRITE (iwfx, "(A)") "<Primitive Centers>"
405 tot_npgf = 0
406 DO iatom = 1, num_atoms
407 iloop = 1
408 NULLIFY (orb_basis_set)
409 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
410 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
411 IF (ASSOCIATED(orb_basis_set)) THEN
412 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
413 nset=nset, &
414 npgf=npgf, &
415 nshell=nshell, &
416 l=l, &
417 gcc=gcc)
418
419 DO iset = 1, nset
420 DO ishell = 1, nshell(iset)
421 lshell = l(ishell, iset)
422 DO ico = 1, nco(lshell)
423 DO ipgf = 1, npgf(iset)
424 tot_npgf = tot_npgf + 1
425 IF (mod(iloop + 10, 10) /= 0) THEN
426 WRITE (iwfx, "(I4)", advance="no") iatom
427 newl = .true.
428 ELSE
429 WRITE (iwfx, "(I4)") iatom
430 newl = .false.
431 END IF
432 iloop = iloop + 1
433 !END IF
434 END DO
435 END DO
436 END DO
437 END DO
438 IF (newl) THEN
439 WRITE (iwfx, *)
440 END IF
441 ELSE
442 cpabort("Unknown basis set type. ")
443 END IF
444 END DO
445 WRITE (iwfx, "(A,/)") "</Primitive Centers>"
446
447 ! !===========================================================================!
448 ! Number of primitives
449 ! !===========================================================================!
450 WRITE (iwfx, "(A)") "<Number of Primitives>"
451 WRITE (iwfx, "(I8)") tot_npgf
452 WRITE (iwfx, "(A,/)") "</Number of Primitives>"
453
454 ! !===========================================================================!
455 ! Primitive Types
456 ! !===========================================================================!
457 WRITE (iwfx, "(A)") "<Primitive Types>"
458 DO iatom = 1, num_atoms
459 iloop = 1
460 NULLIFY (orb_basis_set)
461 NULLIFY (bcgf_symbol)
462 NULLIFY (gcc)
463 NULLIFY (zet)
464 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
465 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
466 IF (ASSOCIATED(orb_basis_set)) THEN
467 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
468 nset=nset, &
469 nshell=nshell, &
470 l=l, &
471 npgf=npgf, &
472 ncgf=ncgf, &
473 zet=zet, &
474 cgf_symbol=bcgf_symbol, &
475 sgf_symbol=bsgf_symbol, &
476 gcc=gcc)
477
478 icgf = 1
479 DO iset = 1, nset
480 DO ishell = 1, nshell(iset)
481 lshell = l(ishell, iset)
482 DO ico = 1, nco(lshell)
483 DO ipgf = 1, orb_basis_set%npgf(iset)
484 IF (mod(iloop + 10, 10) /= 0) THEN
485 WRITE (iwfx, '(I6)', advance="no") &
486 pgf_type_chargemol(bcgf_symbol(icgf) (3:))
487 newl = .true.
488 ELSE
489 WRITE (iwfx, '(I6)') &
490 pgf_type_chargemol(bcgf_symbol(icgf) (3:))
491 newl = .false.
492 END IF
493 iloop = iloop + 1
494 !END IF
495 END DO
496 icgf = icgf + 1
497 END DO !ico
498 END DO ! ishell
499 END DO ! iset
500
501 ELSE
502 cpabort("Unknown basis set type.")
503 END IF
504 IF (newl) THEN
505 WRITE (iwfx, *)
506 END IF
507 END DO ! iatom
508 WRITE (iwfx, "(A,/)") "</Primitive Types>"
509
510 ! !===========================================================================!
511 ! Primitive Exponents
512 ! !===========================================================================!
513 !! NOTE: CP2K orders the atomic orbitals as follows:
514 !! Cartesian orbitals: x (slowest loop), y, z (fastest loop)
515 !! Spherical orbitals: m = -l, ..., 0, ..., +l
516
517 WRITE (iwfx, "(A)") "<Primitive Exponents>"
518 DO iatom = 1, num_atoms
519 NULLIFY (orb_basis_set)
520 NULLIFY (gcc)
521 NULLIFY (zet)
522 iloop = 1
523
524 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
525 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
526 IF (ASSOCIATED(orb_basis_set)) THEN
527 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
528 nset=nset, &
529 npgf=npgf, &
530 nshell=nshell, &
531 l=l, &
532 zet=zet, &
533 gcc=gcc)
534
535 DO iset = 1, nset
536 DO ishell = 1, nshell(iset)
537 lshell = l(ishell, iset)
538 DO ico = 1, nco(lshell)
539 DO ipgf = 1, npgf(iset)
540 IF (mod(iloop + 5, 5) /= 0) THEN
541 WRITE (iwfx, "(ES20.10)", advance="no") zet(ipgf, iset)
542 newl = .true.
543 ELSE
544 WRITE (iwfx, "(ES20.10)") zet(ipgf, iset)
545 newl = .false.
546 END IF
547 iloop = iloop + 1
548 !END IF
549 END DO
550 END DO
551 END DO
552 END DO
553 IF (newl) THEN
554 WRITE (iwfx, *)
555 END IF
556 ELSE
557 cpabort("Unknown basis set type. ")
558 END IF
559 END DO
560 WRITE (iwfx, "(A,/)") "</Primitive Exponents>"
561
562 ! !===========================================================================!
563 ! MO occupation numbers
564 ! nonesense for roks at the moment
565 ! !===========================================================================!
566 WRITE (iwfx, "(A)") "<Molecular Orbital Occupation Numbers>"
567 IF (.NOT. dft_control%roks) THEN
568 DO ispin = 1, nspins
569 DO imo = 1, mos(ispin)%homo
570 WRITE (iwfx, "(f8.6)") &
571 mos(ispin)%occupation_numbers(imo)
572 END DO
573 END DO
574 ELSE
575 DO imo = 1, mos(roks_low)%homo
576 WRITE (iwfx, "(*(f8.6,1x))") &
577 mos(roks_low)%occupation_numbers(imo) &
578 + mos(roks_low)%occupation_numbers(imo)
579 END DO
580 DO imo = mos(roks_low)%homo + 1, mos(roks_high)%homo
581 WRITE (iwfx, "(f8.6)") &
582 mos(roks_high)%occupation_numbers(imo)
583 END DO
584 END IF
585 WRITE (iwfx, "(A,/)") "</Molecular Orbital Occupation Numbers>"
586
587 ! !===========================================================================!
588 ! MO energies
589 ! !===========================================================================!
590 WRITE (iwfx, "(A)") "<Molecular Orbital Energies>"
591 DO ispin = 1, nspins
592 DO imo = 1, mos(ispin)%homo
593 WRITE (iwfx, "(ES20.10)") mos(ispin)%eigenvalues(imo)
594 END DO
595 END DO
596 WRITE (iwfx, "(A,/)") "</Molecular Orbital Energies>"
597
598 ! !===========================================================================!
599 ! MO Spin types
600 ! nonesense for ROKS
601 ! !===========================================================================!
602 WRITE (iwfx, "(A)") "<Molecular Orbital Spin Types>"
603 DO imo = 1, mos(1)%homo
604 IF (nspins == 1) THEN
605 WRITE (iwfx, '(A15)') "Alpha and Beta"
606 ELSEIF (dft_control%uks) THEN
607 WRITE (iwfx, '(A6)') "Alpha"
608 ELSE
609 CALL cp_abort(__location__, "This wavefunction type is currently"// &
610 " not supported by the chargemol interface.")
611 END IF
612 END DO
613 IF (nspins == 2) THEN
614 DO imo = 1, mos(2)%homo
615 WRITE (iwfx, '(A5)') "Beta"
616 END DO
617 END IF
618 WRITE (iwfx, "(A,/)") "</Molecular Orbital Spin Types>"
619
620 ! !===========================================================================!
621 ! Kpoint index of orbitals
622 ! !===========================================================================!
623 IF (periodic) THEN
624 WRITE (iwfx, "(A)") "<Kpoint Index for Orbitals>"
625 DO ispin = 1, nspins
626 DO imo = 1, mos(ispin)%homo
627 WRITE (iwfx, "(I6)") 1
628 END DO
629 END DO
630 WRITE (iwfx, "(A,/)") "</Kpoint Index for Orbitals>"
631 END IF
632
633 ! !===========================================================================!
634 ! MO Primitive Coefficients
635 ! !===========================================================================!
636 WRITE (iwfx, "(A)") "<Molecular Orbital Primitive Coefficients>"
637 NULLIFY (orb_basis_set)
638 NULLIFY (gcc)
639
640 IF (mos(1)%use_mo_coeff_b) THEN
641 DO ispin = 1, SIZE(mos)
642 IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
643 cpassert(.false.)
644 END IF
645 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
646 mos(ispin)%mo_coeff)
647 END DO
648 END IF
649
650 DO ispin = 1, nspins
651 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
652 nrow_global=nrow_global, &
653 ncol_global=ncol_global)
654
655 ALLOCATE (smatrix(nrow_global, ncol_global))
656 smatrix = 0.0_dp
657
658 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
659
660 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
661
662 ALLOCATE (cmatrix(ncgf, ncol_global))
663
664 cmatrix = 0.0_dp
665
666 ! get MO coefficients in terms of contracted cartesian functions
667 icgf = 1
668 isgf = 1
669 DO iatom = 1, num_atoms
670 NULLIFY (orb_basis_set)
671 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
672 CALL get_qs_kind(qs_kind_set(ikind), &
673 basis_set=orb_basis_set)
674 IF (ASSOCIATED(orb_basis_set)) THEN
675 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
676 nset=nset, &
677 nshell=nshell, &
678 l=l)
679 DO iset = 1, nset
680 DO ishell = 1, nshell(iset)
681 lshell = l(ishell, iset)
682 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, &
683 nso(lshell), 1.0_dp, &
684 orbtramat(lshell)%c2s, nso(lshell), &
685 smatrix(isgf, 1), nsgf, 0.0_dp, &
686 cmatrix(icgf, 1), ncgf)
687 !IF (iatom==1 .and. debug_this_module) THEN
688 ! print *, lshell
689 ! print *, "size ", size(orbtramat(lshell)%s2c,1),",",size(orbtramat(lshell)%s2c,2)
690 ! DO itest = 1, size(orbtramat(lshell)%s2c,1)
691 ! print *, orbtramat(lshell)%s2c(itest,:)
692 ! END DO
693 ! print *, " "
694 ! DO itest = 1, 5
695 ! print *, my_d_orbtramat(itest,:)
696 ! END DO
697 !END IF
698 icgf = icgf + nco(lshell)
699 isgf = isgf + nso(lshell)
700 END DO
701 END DO
702 ELSE
703 cpabort("Unknown basis set type. ")
704 END IF
705 END DO ! iatom
706
707 ! Now write MO coefficients in terms of cartesian primitives
708 DO imo = 1, mos(ispin)%homo
709
710 WRITE (iwfx, "(A)") "<MO Number>"
711 IF (nspins == 1 .OR. ispin == 1) THEN
712 WRITE (iwfx, "(I6)") imo
713 ELSE
714 WRITE (iwfx, "(I6)") imo + mos(1)%homo
715 END IF
716 WRITE (iwfx, "(A)") "</MO Number>"
717
718 irow = 1 ! row of cmatrix, each row corresponds to
719 ! a contracted Cartesian function
720 DO iatom = 1, num_atoms
721 iloop = 1
722 NULLIFY (orb_basis_set)
723 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
724 kind_number=ikind)
725 CALL get_qs_kind(qs_kind_set(ikind), &
726 basis_set=orb_basis_set)
727 IF (ASSOCIATED(orb_basis_set)) THEN
728 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
729 nset=nset, &
730 nshell=nshell, &
731 gcc=gcc, &
732 l=l)
733
734 icgf = 1
735 DO iset = 1, nset
736 DO ishell = 1, nshell(iset)
737 lshell = orb_basis_set%l(ishell, iset)
738
739 DO ico = orb_basis_set%first_cgf(ishell, iset), &
740 orb_basis_set%last_cgf(ishell, iset)
741
742 !loop over primitives
743 DO ipgf = 1, orb_basis_set%npgf(iset)
744
745 zeta = orb_basis_set%zet(ipgf, iset)
746
747 IF (mod(iloop + 5, 5) /= 0) THEN
748 WRITE (iwfx, '(ES20.10)', advance="no") &
749 cmatrix(irow, imo)* &
750 orb_basis_set%norm_cgf(ico)* &
751 orb_basis_set%gcc(ipgf, ishell, iset)
752 newl = .true.
753 ELSE
754 WRITE (iwfx, '(ES20.10)') &
755 cmatrix(irow, imo)* &
756 orb_basis_set%norm_cgf(ico)* &
757 orb_basis_set%gcc(ipgf, ishell, iset)
758 newl = .false.
759 END IF
760 iloop = iloop + 1
761 END DO ! end loop over primitives
762 irow = irow + 1
763 END DO !ico
764 END DO ! ishell
765 END DO ! iset
766
767 ELSE
768 cpabort("Unknown basis set type.")
769 END IF
770 IF (newl) THEN
771 WRITE (iwfx, *)
772 END IF
773 END DO ! iatom
774 !Write (iwfx,*)
775 END DO ! imo
776
777 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
778 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
779
780 END DO ! ispin
781 WRITE (iwfx, "(A,/)") "</Molecular Orbital Primitive Coefficients>"
782
783 ! !===========================================================================!
784 ! Energy
785 ! !===========================================================================!
786 WRITE (iwfx, "(A)") "<Energy>"
787 WRITE (iwfx, "(ES26.16E3)") energy%total
788 WRITE (iwfx, "(A,/)") "</Energy>"
789
790 ! !===========================================================================!
791 ! Virial ratio
792 ! !===========================================================================!
793 WRITE (iwfx, "(A)") "<Virial Ratio (-V/T)>"
794 WRITE (iwfx, "(ES20.12)") - 1*(energy%total - energy%kinetic)/energy%kinetic
795 WRITE (iwfx, "(A,/)") "</Virial Ratio (-V/T)>"
796
797 ! !===========================================================================!
798 ! Force-related quantities
799 ! inactivated for now
800 ! !===========================================================================!
801 IF (ASSOCIATED(force) .AND. debug_this_module) THEN
802 WRITE (iwfx, "(A)") "<Nuclear Cartesian Energy Gradients>"
803
804 CALL get_qs_env(qs_env, force=force, atomic_kind_set=atomic_kind_set)
805 ALLOCATE (atom_of_kind(num_atoms))
806 ALLOCATE (kind_of(num_atoms))
807 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
808 atom_of_kind=atom_of_kind, &
809 kind_of=kind_of)
810
811 DO iatom = 1, num_atoms
812 ikind = kind_of(iatom)
813 i = atom_of_kind(iatom)
814 WRITE (iwfx, "(3ES20.12E3)") force(ikind)%total(1:3, i)
815 END DO
816
817 WRITE (iwfx, "(A,/)") "<Nuclear Cartesian Energy Gradients>"
818
819 ! W
820 WRITE (iwfx, "(A)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
821 WRITE (iwfx, "(A)") ""
822 WRITE (iwfx, "(A,/)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
823
824 ! Full Virial ratio
825 WRITE (iwfx, "(A)") "<Full Virial Ratio, -(V - W)/T)>"
826 WRITE (iwfx, "(A)") ""
827 WRITE (iwfx, "(A,/)") "</Full Virial Ratio, -(V - W)/T)>"
828
829 DEALLOCATE (atom_of_kind)
830 DEALLOCATE (kind_of)
831
832 END IF
833
834 ! !===========================================================================!
835 ! Core-density
836 ! !===========================================================================!
837 ! Additional Electron Density Function (EDF)
838 ! WRITE (iwfx, "(A)") "#<Additional Electron Density Function (EDF)>"
839 ! WRITE (iwfx, "(A,/)") "#</Additional Electron Density Function (EDF)>"
840
841 ! !-----------------------------------!
842 ! Done
843 ! !-----------------------------------!
844 DEALLOCATE (atoms_cart)
845 DEALLOCATE (atom_symbols)
846 DEALLOCATE (atomic_numbers)
847 DEALLOCATE (core_charges)
848
849 ELSE ! no output unit
850 iwfx = -1
851 END IF
852
853 IF (output_unit > 0) THEN
854 WRITE (output_unit, '(/,T2,A)') &
855 '!--------------------------- Chargemol wfx written ---------------------------!'
856 END IF
857
858 CALL timestop(handle)
859
860 END SUBROUTINE write_wfx
861
862! **************************************************************************************************
863!> \brief ...
864!> \param l_symb ...
865!> \return ...
866! **************************************************************************************************
867 INTEGER FUNCTION pgf_type_chargemol(l_symb) RESULT(label)
868 !INTEGER, INTENT(OUT) :: label
869 CHARACTER(len=*), INTENT(IN) :: l_symb
870
871 SELECT CASE (l_symb)
872 CASE ("s")
873 label = 1
874 CASE ("px")
875 label = 2
876 CASE ("py")
877 label = 3
878 CASE ("pz")
879 label = 4
880 CASE ("dx2")
881 label = 5
882 CASE ("dy2")
883 label = 6
884 CASE ("dz2")
885 label = 7
886 CASE ("dxy")
887 label = 8
888 CASE ("dxz")
889 label = 9
890 CASE ("dyz")
891 label = 10
892 CASE ("fx3")
893 label = 11
894 CASE ("fy3")
895 label = 12
896 CASE ("fz3")
897 label = 13
898 CASE ("fx2y")
899 label = 14
900 CASE ("fx2z")
901 label = 15
902 CASE ("fy2z")
903 label = 16
904 CASE ("fxy2")
905 label = 17
906 CASE ("fxz2")
907 label = 18
908 CASE ("fyz2")
909 label = 19
910 CASE ("fxyz")
911 label = 20
912 CASE ("gx4")
913 label = 21
914 CASE ("gy4")
915 label = 22
916 CASE ("gz4")
917 label = 23
918 CASE ("gx3y")
919 label = 24
920 CASE ("gx3z")
921 label = 25
922 CASE ("gxy3")
923 label = 26
924 CASE ("gy3z")
925 label = 27
926 CASE ("gxz3")
927 label = 28
928 CASE ("gyz3")
929 label = 29
930 CASE ("gx2y2")
931 label = 30
932 CASE ("gx2z2")
933 label = 31
934 CASE ("gy2z2")
935 label = 32
936 CASE ("gx2yz")
937 label = 33
938 CASE ("gxy2z")
939 label = 34
940 CASE ("gxyz2")
941 label = 35
942 CASE ("hz5")
943 label = 36
944 CASE ("hyz4")
945 label = 37
946 CASE ("hy2z3")
947 label = 38
948 CASE ("hy3z2")
949 label = 39
950 CASE ("hy4z")
951 label = 40
952 CASE ("hy5")
953 label = 41
954 CASE ("hxz4")
955 label = 42
956 CASE ("hxyz3")
957 label = 43
958 CASE ("hxy2z2")
959 label = 44
960 CASE ("hxy3z")
961 label = 45
962 CASE ("hxy4")
963 label = 46
964 CASE ("hx2z3")
965 label = 47
966 CASE ("hx2yz2")
967 label = 48
968 CASE ("hx2y2z")
969 label = 49
970 CASE ("hx2y3")
971 label = 50
972 CASE ("hx3z2")
973 label = 51
974 CASE ("hx3yz")
975 label = 52
976 CASE ("hx3y2")
977 label = 53
978 CASE ("hx4z")
979 label = 54
980 CASE ("hx4y")
981 label = 55
982 CASE ("hx5")
983 label = 56
984 CASE default
985 cpabort("The chargemol interface supports basis functions up to l=5 only")
986 !label = 99
987 END SELECT
988
989 END FUNCTION pgf_type_chargemol
990
991END MODULE qs_chargemol
992
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)
...
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_string_length
Definition kinds.F:57
Types and basic routines needed for a kpoint calculation.
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_datum(cal_date)
returns a datum in human readable format using a standard Fortran routine
Definition machine.F:270
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_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, 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, rhs)
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, 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_r3d_rs_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_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)
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.