(git:98357aa)
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-2026 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 ELSE IF (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 cpassert(ASSOCIATED(mos(ispin)%mo_coeff_b))
644 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
645 mos(ispin)%mo_coeff)
646 END DO
647 END IF
648
649 DO ispin = 1, nspins
650 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
651 nrow_global=nrow_global, &
652 ncol_global=ncol_global)
653
654 ALLOCATE (smatrix(nrow_global, ncol_global))
655 smatrix = 0.0_dp
656
657 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
658
659 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
660
661 ALLOCATE (cmatrix(ncgf, ncol_global))
662
663 cmatrix = 0.0_dp
664
665 ! get MO coefficients in terms of contracted cartesian functions
666 icgf = 1
667 isgf = 1
668 DO iatom = 1, num_atoms
669 NULLIFY (orb_basis_set)
670 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
671 CALL get_qs_kind(qs_kind_set(ikind), &
672 basis_set=orb_basis_set)
673 IF (ASSOCIATED(orb_basis_set)) THEN
674 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
675 nset=nset, &
676 nshell=nshell, &
677 l=l)
678 DO iset = 1, nset
679 DO ishell = 1, nshell(iset)
680 lshell = l(ishell, iset)
681 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, &
682 nso(lshell), 1.0_dp, &
683 orbtramat(lshell)%c2s, nso(lshell), &
684 smatrix(isgf, 1), nsgf, 0.0_dp, &
685 cmatrix(icgf, 1), ncgf)
686 !IF (iatom==1 .and. debug_this_module) THEN
687 ! print *, lshell
688 ! print *, "size ", size(orbtramat(lshell)%s2c,1),",",size(orbtramat(lshell)%s2c,2)
689 ! DO itest = 1, size(orbtramat(lshell)%s2c,1)
690 ! print *, orbtramat(lshell)%s2c(itest,:)
691 ! END DO
692 ! print *, " "
693 ! DO itest = 1, 5
694 ! print *, my_d_orbtramat(itest,:)
695 ! END DO
696 !END IF
697 icgf = icgf + nco(lshell)
698 isgf = isgf + nso(lshell)
699 END DO
700 END DO
701 ELSE
702 cpabort("Unknown basis set type. ")
703 END IF
704 END DO ! iatom
705
706 ! Now write MO coefficients in terms of cartesian primitives
707 DO imo = 1, mos(ispin)%homo
708
709 WRITE (iwfx, "(A)") "<MO Number>"
710 IF (nspins == 1 .OR. ispin == 1) THEN
711 WRITE (iwfx, "(I6)") imo
712 ELSE
713 WRITE (iwfx, "(I6)") imo + mos(1)%homo
714 END IF
715 WRITE (iwfx, "(A)") "</MO Number>"
716
717 irow = 1 ! row of cmatrix, each row corresponds to
718 ! a contracted Cartesian function
719 DO iatom = 1, num_atoms
720 iloop = 1
721 NULLIFY (orb_basis_set)
722 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
723 kind_number=ikind)
724 CALL get_qs_kind(qs_kind_set(ikind), &
725 basis_set=orb_basis_set)
726 IF (ASSOCIATED(orb_basis_set)) THEN
727 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
728 nset=nset, &
729 nshell=nshell, &
730 gcc=gcc, &
731 l=l)
732
733 icgf = 1
734 DO iset = 1, nset
735 DO ishell = 1, nshell(iset)
736 lshell = orb_basis_set%l(ishell, iset)
737
738 DO ico = orb_basis_set%first_cgf(ishell, iset), &
739 orb_basis_set%last_cgf(ishell, iset)
740
741 !loop over primitives
742 DO ipgf = 1, orb_basis_set%npgf(iset)
743
744 zeta = orb_basis_set%zet(ipgf, iset)
745
746 IF (mod(iloop + 5, 5) /= 0) THEN
747 WRITE (iwfx, '(ES20.10)', advance="no") &
748 cmatrix(irow, imo)* &
749 orb_basis_set%norm_cgf(ico)* &
750 orb_basis_set%gcc(ipgf, ishell, iset)
751 newl = .true.
752 ELSE
753 WRITE (iwfx, '(ES20.10)') &
754 cmatrix(irow, imo)* &
755 orb_basis_set%norm_cgf(ico)* &
756 orb_basis_set%gcc(ipgf, ishell, iset)
757 newl = .false.
758 END IF
759 iloop = iloop + 1
760 END DO ! end loop over primitives
761 irow = irow + 1
762 END DO !ico
763 END DO ! ishell
764 END DO ! iset
765
766 ELSE
767 cpabort("Unknown basis set type.")
768 END IF
769 IF (newl) THEN
770 WRITE (iwfx, *)
771 END IF
772 END DO ! iatom
773 !Write (iwfx,*)
774 END DO ! imo
775
776 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
777 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
778
779 END DO ! ispin
780 WRITE (iwfx, "(A,/)") "</Molecular Orbital Primitive Coefficients>"
781
782 ! !===========================================================================!
783 ! Energy
784 ! !===========================================================================!
785 WRITE (iwfx, "(A)") "<Energy>"
786 WRITE (iwfx, "(ES26.16E3)") energy%total
787 WRITE (iwfx, "(A,/)") "</Energy>"
788
789 ! !===========================================================================!
790 ! Virial ratio
791 ! !===========================================================================!
792 WRITE (iwfx, "(A)") "<Virial Ratio (-V/T)>"
793 WRITE (iwfx, "(ES20.12)") - 1*(energy%total - energy%kinetic)/energy%kinetic
794 WRITE (iwfx, "(A,/)") "</Virial Ratio (-V/T)>"
795
796 ! !===========================================================================!
797 ! Force-related quantities
798 ! inactivated for now
799 ! !===========================================================================!
800 IF (ASSOCIATED(force) .AND. debug_this_module) THEN
801 WRITE (iwfx, "(A)") "<Nuclear Cartesian Energy Gradients>"
802
803 CALL get_qs_env(qs_env, force=force, atomic_kind_set=atomic_kind_set)
804 ALLOCATE (atom_of_kind(num_atoms))
805 ALLOCATE (kind_of(num_atoms))
806 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
807 atom_of_kind=atom_of_kind, &
808 kind_of=kind_of)
809
810 DO iatom = 1, num_atoms
811 ikind = kind_of(iatom)
812 i = atom_of_kind(iatom)
813 WRITE (iwfx, "(3ES20.12E3)") force(ikind)%total(1:3, i)
814 END DO
815
816 WRITE (iwfx, "(A,/)") "<Nuclear Cartesian Energy Gradients>"
817
818 ! W
819 WRITE (iwfx, "(A)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
820 WRITE (iwfx, "(A)") ""
821 WRITE (iwfx, "(A,/)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
822
823 ! Full Virial ratio
824 WRITE (iwfx, "(A)") "<Full Virial Ratio, -(V - W)/T)>"
825 WRITE (iwfx, "(A)") ""
826 WRITE (iwfx, "(A,/)") "</Full Virial Ratio, -(V - W)/T)>"
827
828 DEALLOCATE (atom_of_kind)
829 DEALLOCATE (kind_of)
830
831 END IF
832
833 ! !===========================================================================!
834 ! Core-density
835 ! !===========================================================================!
836 ! Additional Electron Density Function (EDF)
837 ! WRITE (iwfx, "(A)") "#<Additional Electron Density Function (EDF)>"
838 ! WRITE (iwfx, "(A,/)") "#</Additional Electron Density Function (EDF)>"
839
840 ! !-----------------------------------!
841 ! Done
842 ! !-----------------------------------!
843 DEALLOCATE (atoms_cart)
844 DEALLOCATE (atom_symbols)
845 DEALLOCATE (atomic_numbers)
846 DEALLOCATE (core_charges)
847
848 ELSE ! no output unit
849 iwfx = -1
850 END IF
851
852 IF (output_unit > 0) THEN
853 WRITE (output_unit, '(/,T2,A)') &
854 '!--------------------------- Chargemol wfx written ---------------------------!'
855 END IF
856
857 CALL timestop(handle)
858
859 END SUBROUTINE write_wfx
860
861! **************************************************************************************************
862!> \brief ...
863!> \param l_symb ...
864!> \return ...
865! **************************************************************************************************
866 INTEGER FUNCTION pgf_type_chargemol(l_symb) RESULT(label)
867 !INTEGER, INTENT(OUT) :: label
868 CHARACTER(len=*), INTENT(IN) :: l_symb
869
870 SELECT CASE (l_symb)
871 CASE ("s")
872 label = 1
873 CASE ("px")
874 label = 2
875 CASE ("py")
876 label = 3
877 CASE ("pz")
878 label = 4
879 CASE ("dx2")
880 label = 5
881 CASE ("dy2")
882 label = 6
883 CASE ("dz2")
884 label = 7
885 CASE ("dxy")
886 label = 8
887 CASE ("dxz")
888 label = 9
889 CASE ("dyz")
890 label = 10
891 CASE ("fx3")
892 label = 11
893 CASE ("fy3")
894 label = 12
895 CASE ("fz3")
896 label = 13
897 CASE ("fx2y")
898 label = 14
899 CASE ("fx2z")
900 label = 15
901 CASE ("fy2z")
902 label = 16
903 CASE ("fxy2")
904 label = 17
905 CASE ("fxz2")
906 label = 18
907 CASE ("fyz2")
908 label = 19
909 CASE ("fxyz")
910 label = 20
911 CASE ("gx4")
912 label = 21
913 CASE ("gy4")
914 label = 22
915 CASE ("gz4")
916 label = 23
917 CASE ("gx3y")
918 label = 24
919 CASE ("gx3z")
920 label = 25
921 CASE ("gxy3")
922 label = 26
923 CASE ("gy3z")
924 label = 27
925 CASE ("gxz3")
926 label = 28
927 CASE ("gyz3")
928 label = 29
929 CASE ("gx2y2")
930 label = 30
931 CASE ("gx2z2")
932 label = 31
933 CASE ("gy2z2")
934 label = 32
935 CASE ("gx2yz")
936 label = 33
937 CASE ("gxy2z")
938 label = 34
939 CASE ("gxyz2")
940 label = 35
941 CASE ("hz5")
942 label = 36
943 CASE ("hyz4")
944 label = 37
945 CASE ("hy2z3")
946 label = 38
947 CASE ("hy3z2")
948 label = 39
949 CASE ("hy4z")
950 label = 40
951 CASE ("hy5")
952 label = 41
953 CASE ("hxz4")
954 label = 42
955 CASE ("hxyz3")
956 label = 43
957 CASE ("hxy2z2")
958 label = 44
959 CASE ("hxy3z")
960 label = 45
961 CASE ("hxy4")
962 label = 46
963 CASE ("hx2z3")
964 label = 47
965 CASE ("hx2yz2")
966 label = 48
967 CASE ("hx2y2z")
968 label = 49
969 CASE ("hx2y3")
970 label = 50
971 CASE ("hx3z2")
972 label = 51
973 CASE ("hx3yz")
974 label = 52
975 CASE ("hx3y2")
976 label = 53
977 CASE ("hx4z")
978 label = 54
979 CASE ("hx4y")
980 label = 55
981 CASE ("hx5")
982 label = 56
983 CASE default
984 cpabort("The chargemol interface supports basis functions up to l=5 only")
985 !label = 99
986 END SELECT
987
988 END FUNCTION pgf_type_chargemol
989
990END MODULE qs_chargemol
991
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, ccon)
...
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:233
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:311
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:46
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:381
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, mimic, 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, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, cneo_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, 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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
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:60
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.