LCOV - code coverage report
Current view: top level - src - molden_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 80.6 % 170 137
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            Line data    Source code
       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  Functions handling the MOLDEN format. Split from mode_selective.
      10              : !> \author Teodoro Laino, 03.2009
      11              : ! **************************************************************************************************
      12              : MODULE molden_utils
      13              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      14              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      15              :                                               gto_basis_set_type
      16              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm
      17              :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      18              :                                               cp_fm_get_submatrix
      19              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      20              :                                               cp_logger_type
      21              :    USE cp_output_handling,              ONLY: cp_p_file,&
      22              :                                               cp_print_key_finished_output,&
      23              :                                               cp_print_key_should_output,&
      24              :                                               cp_print_key_unit_nr
      25              :    USE input_constants,                 ONLY: gto_cartesian,&
      26              :                                               gto_spherical
      27              :    USE input_section_types,             ONLY: section_vals_type,&
      28              :                                               section_vals_val_get
      29              :    USE kinds,                           ONLY: dp
      30              :    USE mathconstants,                   ONLY: pi
      31              :    USE orbital_pointers,                ONLY: nco,&
      32              :                                               nso
      33              :    USE orbital_transformation_matrices, ONLY: orbtramat
      34              :    USE particle_types,                  ONLY: particle_type
      35              :    USE periodic_table,                  ONLY: get_ptable_info
      36              :    USE physcon,                         ONLY: massunit
      37              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      38              :                                               get_qs_kind_set,&
      39              :                                               qs_kind_type
      40              :    USE qs_mo_types,                     ONLY: mo_set_type
      41              : #include "./base/base_uses.f90"
      42              : 
      43              :    IMPLICIT NONE
      44              : 
      45              :    PRIVATE
      46              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molden_utils'
      47              :    LOGICAL, PARAMETER                   :: debug_this_module = .FALSE.
      48              : 
      49              :    INTEGER, PARAMETER                   :: molden_lmax = 4
      50              :    INTEGER, PARAMETER                   :: molden_ncomax = (molden_lmax + 1)*(molden_lmax + 2)/2 ! 15
      51              : 
      52              :    PUBLIC :: write_vibrations_molden, write_mos_molden
      53              : 
      54              : CONTAINS
      55              : 
      56              : ! **************************************************************************************************
      57              : !> \brief Write out the MOs in molden format for visualisation
      58              : !> \param mos the set of MOs (both spins, if UKS)
      59              : !> \param qs_kind_set for basis set info
      60              : !> \param particle_set particles data structure, for positions and kinds
      61              : !> \param print_section input section containing relevant print key
      62              : !> \author MattW, IainB
      63              : ! **************************************************************************************************
      64         9943 :    SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section)
      65              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
      66              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      67              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      68              :       TYPE(section_vals_type), POINTER                   :: print_section
      69              : 
      70              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_mos_molden'
      71              :       CHARACTER(LEN=molden_lmax+1), PARAMETER            :: angmom = "spdfg"
      72              : 
      73              :       CHARACTER(LEN=15)                                  :: fmtstr1, fmtstr2
      74              :       CHARACTER(LEN=2)                                   :: element_symbol
      75              :       INTEGER :: gto_kind, handle, i, iatom, icgf, icol, ikind, ipgf, irow, irow_in, iset, isgf, &
      76              :          ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, z
      77         9943 :       INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
      78         9943 :       INTEGER, DIMENSION(:, :), POINTER                  :: l
      79              :       INTEGER, DIMENSION(molden_ncomax, 0:molden_lmax)   :: orbmap
      80              :       LOGICAL                                            :: print_warn
      81              :       REAL(KIND=dp)                                      :: expzet, prefac
      82         9943 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cmatrix, smatrix
      83         9943 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
      84         9943 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
      85              :       TYPE(cp_logger_type), POINTER                      :: logger
      86              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
      87              : 
      88         9943 :       CALL timeset(routineN, handle)
      89              : 
      90         9943 :       logger => cp_get_default_logger()
      91         9943 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, ""), cp_p_file)) THEN
      92              : 
      93              :          iw = cp_print_key_unit_nr(logger, print_section, "", &
      94           28 :                                    extension=".molden", file_status='REPLACE')
      95              : 
      96           28 :          print_warn = .TRUE.
      97              : 
      98           28 :          CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
      99           28 :          ndigits = MIN(MAX(3, ndigits), 30)
     100           28 :          WRITE (UNIT=fmtstr1, FMT='("(I6,1X,ES",I0,".",I0,")")') ndigits + 7, ndigits
     101           28 :          WRITE (UNIT=fmtstr2, FMT='("((T51,2F",I0,".",I0,"))")') ndigits + 10, ndigits
     102              : 
     103           28 :          CALL section_vals_val_get(print_section, "GTO_KIND", i_val=gto_kind)
     104              : 
     105           28 :          IF (mos(1)%use_mo_coeff_b) THEN
     106              :             ! we are using the dbcsr mo_coeff
     107              :             ! we copy it to the fm anyway
     108            0 :             DO ispin = 1, SIZE(mos)
     109            0 :                IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
     110            0 :                   CPASSERT(.FALSE.)
     111              :                END IF
     112              :                CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
     113            0 :                                      mos(ispin)%mo_coeff) !fm->dbcsr
     114              :             END DO
     115              :          END IF
     116              : 
     117           28 :          IF (iw > 0) THEN
     118           14 :             WRITE (iw, '(T2,A)') "[Molden Format]"
     119           14 :             WRITE (iw, '(T2,A)') "[Atoms] AU"
     120          152 :             DO i = 1, SIZE(particle_set)
     121              :                CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
     122          138 :                                     element_symbol=element_symbol)
     123          138 :                CALL get_ptable_info(element_symbol, number=z)
     124              : 
     125              :                WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
     126          566 :                   element_symbol, i, z, particle_set(i)%r(:)
     127              :             END DO
     128              : 
     129           14 :             WRITE (iw, '(T2,A)') "[GTO]"
     130              : 
     131          152 :             DO i = 1, SIZE(particle_set)
     132              :                CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
     133          138 :                                     element_symbol=element_symbol)
     134          138 :                CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     135          290 :                IF (ASSOCIATED(orb_basis_set)) THEN
     136          138 :                   WRITE (iw, '(T2,I8,I8)') i, 0
     137              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     138              :                                          nset=nset, &
     139              :                                          npgf=npgf, &
     140              :                                          nshell=nshell, &
     141              :                                          l=l, &
     142              :                                          zet=zet, &
     143          138 :                                          gcc=gcc)
     144              : 
     145          434 :                   DO iset = 1, nset
     146          784 :                      DO ishell = 1, nshell(iset)
     147          350 :                         lshell = l(ishell, iset)
     148          646 :                         IF (lshell <= molden_lmax) THEN
     149              :                            WRITE (UNIT=iw, FMT='(T25,A2,4X,I4,4X,F4.2)') &
     150          350 :                               angmom(lshell + 1:lshell + 1), npgf(iset), 1.0_dp
     151              :                            ! MOLDEN expects the contraction coefficient of spherical NOT CARTESIAN NORMALISED
     152              :                            ! functions. So we undo the normalisation factors included in the gccs
     153              :                            ! Reverse engineered from basis_set_types, normalise_gcc_orb
     154          350 :                            prefac = 2_dp**lshell*(2/pi)**0.75_dp
     155          350 :                            expzet = 0.25_dp*(2*lshell + 3.0_dp)
     156              :                            WRITE (UNIT=iw, FMT=fmtstr2) &
     157         2264 :                               (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
     158         2614 :                                ipgf=1, npgf(iset))
     159              :                         ELSE
     160            0 :                            IF (print_warn) THEN
     161              :                               CALL cp_warn(__LOCATION__, &
     162            0 :                                            "MOLDEN format does not support Gaussian orbitals with l > 4.")
     163            0 :                               print_warn = .FALSE.
     164              :                            END IF
     165              :                         END IF
     166              :                      END DO
     167              :                   END DO
     168              : 
     169          138 :                   WRITE (iw, '(A4)') "    "
     170              : 
     171              :                END IF
     172              : 
     173              :             END DO
     174              : 
     175           14 :             IF (gto_kind == gto_spherical) THEN
     176           14 :                WRITE (iw, '(T2,A)') "[5D7F]"
     177           14 :                WRITE (iw, '(T2,A)') "[9G]"
     178              :             END IF
     179              : 
     180           14 :             WRITE (iw, '(T2,A)') "[MO]"
     181              :          END IF
     182              : 
     183              :          !------------------------------------------------------------------------
     184              :          ! convert from CP2K to MOLDEN format ordering
     185              :          ! http://www.cmbi.ru.nl/molden/molden_format.html
     186              :          !"The following order of D, F and G functions is expected:
     187              :          !
     188              :          !   5D: D 0, D+1, D-1, D+2, D-2
     189              :          !   6D: xx, yy, zz, xy, xz, yz
     190              :          !
     191              :          !   7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
     192              :          !  10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
     193              :          !
     194              :          !   9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
     195              :          !  15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
     196              :          !       xxyy xxzz yyzz xxyz yyxz zzxy
     197              :          !"
     198              :          ! CP2K has x in the outer (slower loop), so
     199              :          ! xx, xy, xz, yy, yz,zz for l=2, for instance
     200              :          !
     201              :          ! iorb_cp2k = orbmap(iorb_molden, l), l = 0 .. 4
     202              :          ! -----------------------------------------------------------------------
     203           28 :          IF (iw > 0) THEN
     204           14 :             IF (gto_kind == gto_cartesian) THEN
     205              :                ! -----------------------------------------------------------------
     206              :                ! Use cartesian (6D, 10F, 15G) representation.
     207              :                ! This is only format VMD can process.
     208              :                ! -----------------------------------------------------------------
     209              :                orbmap = RESHAPE((/1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     210              :                                   1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     211              :                                   1, 4, 6, 2, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     212              :                                   1, 7, 10, 4, 2, 3, 6, 9, 8, 5, 0, 0, 0, 0, 0, &
     213              :                                   1, 11, 15, 2, 3, 7, 12, 10, 14, 4, 6, 13, 5, 8, 9/), &
     214            0 :                                 (/molden_ncomax, molden_lmax + 1/))
     215           14 :             ELSE IF (gto_kind == gto_spherical) THEN
     216              :                ! -----------------------------------------------------------------
     217              :                ! Use spherical (5D, 7F, 9G) representation.
     218              :                ! -----------------------------------------------------------------
     219              :                orbmap = RESHAPE((/1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     220              :                                   3, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     221              :                                   3, 4, 2, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
     222              :                                   4, 5, 3, 6, 2, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
     223              :                                   5, 6, 4, 7, 3, 8, 2, 9, 1, 0, 0, 0, 0, 0, 0/), &
     224           14 :                                 (/molden_ncomax, molden_lmax + 1/))
     225              :             END IF
     226              :          END IF
     227              : 
     228           72 :          DO ispin = 1, SIZE(mos)
     229              :             CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
     230              :                                 nrow_global=nrow_global, &
     231           44 :                                 ncol_global=ncol_global)
     232          176 :             ALLOCATE (smatrix(nrow_global, ncol_global))
     233           44 :             CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
     234              : 
     235           44 :             IF (iw > 0) THEN
     236           22 :                IF (gto_kind == gto_cartesian) THEN
     237            0 :                   CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
     238              : 
     239            0 :                   ALLOCATE (cmatrix(ncgf, ncgf))
     240              : 
     241            0 :                   cmatrix = 0.0_dp
     242              : 
     243              :                   ! Transform spherical MOs to Cartesian MOs
     244              : 
     245            0 :                   icgf = 1
     246            0 :                   isgf = 1
     247            0 :                   DO iatom = 1, SIZE(particle_set)
     248            0 :                      NULLIFY (orb_basis_set)
     249            0 :                      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     250              :                      CALL get_qs_kind(qs_kind_set(ikind), &
     251            0 :                                       basis_set=orb_basis_set)
     252            0 :                      IF (ASSOCIATED(orb_basis_set)) THEN
     253              :                         CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     254              :                                                nset=nset, &
     255              :                                                nshell=nshell, &
     256            0 :                                                l=l)
     257            0 :                         DO iset = 1, nset
     258            0 :                            DO ishell = 1, nshell(iset)
     259            0 :                               lshell = l(ishell, iset)
     260              :                               CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, nso(lshell), 1.0_dp, &
     261              :                                          orbtramat(lshell)%c2s, nso(lshell), &
     262              :                                          smatrix(isgf, 1), nsgf, 0.0_dp, &
     263            0 :                                          cmatrix(icgf, 1), ncgf)
     264            0 :                               icgf = icgf + nco(lshell)
     265            0 :                               isgf = isgf + nso(lshell)
     266              :                            END DO
     267              :                         END DO
     268              :                      END IF
     269              :                   END DO ! iatom
     270              :                END IF
     271              : 
     272           96 :                DO icol = 1, mos(ispin)%nmo
     273              :                   ! index of the first basis function for the given atom, set, and shell
     274           74 :                   irow = 1
     275              : 
     276              :                   ! index of the first basis function in MOLDEN file.
     277              :                   ! Due to limitation of the MOLDEN format, basis functions with l > molden_lmax
     278              :                   ! cannot be exported, so we need to renumber atomic orbitals
     279           74 :                   irow_in = 1
     280              : 
     281           74 :                   WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
     282           74 :                   IF (ispin < 2) THEN
     283           63 :                      WRITE (iw, '(A)') 'Spin= Alpha'
     284              :                   ELSE
     285           11 :                      WRITE (iw, '(A)') 'Spin= Beta'
     286              :                   END IF
     287           74 :                   WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
     288              : 
     289          544 :                   DO iatom = 1, SIZE(particle_set)
     290          448 :                      NULLIFY (orb_basis_set)
     291              :                      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
     292          448 :                                           element_symbol=element_symbol, kind_number=ikind)
     293              :                      CALL get_qs_kind(qs_kind_set(ikind), &
     294          448 :                                       basis_set=orb_basis_set)
     295          970 :                      IF (ASSOCIATED(orb_basis_set)) THEN
     296              :                         CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     297              :                                                nset=nset, &
     298              :                                                nshell=nshell, &
     299          448 :                                                l=l)
     300              : 
     301          448 :                         IF (gto_kind == gto_cartesian) THEN
     302              :                            ! ----------------------------------------------
     303              :                            ! Use cartesian (6D, 10F, 15G) representation.
     304              :                            ! ----------------------------------------------
     305            0 :                            icgf = 1
     306            0 :                            DO iset = 1, nset
     307            0 :                               DO ishell = 1, nshell(iset)
     308            0 :                                  lshell = l(ishell, iset)
     309              : 
     310            0 :                                  IF (lshell <= molden_lmax) THEN
     311              :                                     CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
     312            0 :                                                       cmatrix(irow:irow + nco(lshell) - 1, icol))
     313            0 :                                     irow_in = irow_in + nco(lshell)
     314              :                                  END IF
     315              : 
     316            0 :                                  irow = irow + nco(lshell)
     317              :                               END DO ! ishell
     318              :                            END DO
     319              : 
     320          448 :                         ELSE IF (gto_kind == gto_spherical) THEN
     321              :                            ! ----------------------------------------------
     322              :                            ! Use spherical (5D, 7F, 9G) representation.
     323              :                            ! ----------------------------------------------
     324         1544 :                            DO iset = 1, nset
     325         2880 :                               DO ishell = 1, nshell(iset)
     326         1336 :                                  lshell = l(ishell, iset)
     327              : 
     328         1336 :                                  IF (lshell <= molden_lmax) THEN
     329              :                                     CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
     330         1336 :                                                       smatrix(irow:irow + nso(lshell) - 1, icol))
     331         1336 :                                     irow_in = irow_in + nso(lshell)
     332              :                                  END IF
     333              : 
     334         2432 :                                  irow = irow + nso(lshell)
     335              :                               END DO
     336              :                            END DO
     337              :                         END IF
     338              : 
     339              :                      END IF
     340              :                   END DO ! iatom
     341              :                END DO
     342              :             END IF
     343              : 
     344           44 :             IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
     345          116 :             IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
     346              :          END DO
     347              : 
     348           28 :          CALL cp_print_key_finished_output(iw, logger, print_section, "")
     349              : 
     350              :       END IF
     351              : 
     352         9943 :       CALL timestop(handle)
     353              : 
     354        19886 :    END SUBROUTINE write_mos_molden
     355              : 
     356              : ! **************************************************************************************************
     357              : !> \brief Output MO coefficients formatted correctly for MOLDEN, omitting those <= 1E(-digits)
     358              : !> \param iw       output file unit
     359              : !> \param fmtstr1  format string
     360              : !> \param ndigits  number of significant digits in MO coefficients
     361              : !> \param irow_in  index of the first atomic orbital: mo_coeff(orbmap(1))
     362              : !> \param orbmap   array to map Gaussian functions from MOLDEN to CP2K ordering
     363              : !> \param mo_coeff MO coefficients
     364              : ! **************************************************************************************************
     365         1336 :    SUBROUTINE print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap, mo_coeff)
     366              :       INTEGER, INTENT(in)                                :: iw
     367              :       CHARACTER(LEN=*), INTENT(in)                       :: fmtstr1
     368              :       INTEGER, INTENT(in)                                :: ndigits, irow_in
     369              :       INTEGER, DIMENSION(molden_ncomax), INTENT(in)      :: orbmap
     370              :       REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: mo_coeff
     371              : 
     372              :       INTEGER                                            :: orbital
     373              : 
     374        21376 :       DO orbital = 1, molden_ncomax
     375        21376 :          IF (orbmap(orbital) /= 0) THEN
     376         2504 :             IF (ABS(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
     377         1573 :                WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
     378              :             END IF
     379              :          END IF
     380              :       END DO
     381              : 
     382         1336 :    END SUBROUTINE print_coeffs
     383              : 
     384              : ! **************************************************************************************************
     385              : !> \brief writes the output for vibrational analysis in MOLDEN format
     386              : !> \param input ...
     387              : !> \param particles ...
     388              : !> \param freq ...
     389              : !> \param eigen_vec ...
     390              : !> \param intensities ...
     391              : !> \param calc_intens ...
     392              : !> \param dump_only_positive ...
     393              : !> \param logger ...
     394              : !> \param list array of mobile atom indices
     395              : !> \author Florian Schiffmann 11.2007
     396              : ! **************************************************************************************************
     397           54 :    SUBROUTINE write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, &
     398              :                                       dump_only_positive, logger, list)
     399              : 
     400              :       TYPE(section_vals_type), POINTER                   :: input
     401              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     402              :       REAL(KIND=dp), DIMENSION(:)                        :: freq
     403              :       REAL(KIND=dp), DIMENSION(:, :)                     :: eigen_vec
     404              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: intensities
     405              :       LOGICAL, INTENT(in)                                :: calc_intens, dump_only_positive
     406              :       TYPE(cp_logger_type), POINTER                      :: logger
     407              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: list
     408              : 
     409              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_vibrations_molden'
     410              : 
     411              :       CHARACTER(LEN=2)                                   :: element_symbol
     412              :       INTEGER                                            :: handle, i, iw, j, k, l, z
     413           54 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: my_list
     414              :       REAL(KIND=dp)                                      :: fint
     415              : 
     416           54 :       CALL timeset(routineN, handle)
     417              : 
     418              :       iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
     419           54 :                                 extension=".mol", file_status='REPLACE')
     420              : 
     421           54 :       IF (iw .GT. 0) THEN
     422           27 :          CPASSERT(MOD(SIZE(eigen_vec, 1), 3) == 0)
     423           27 :          CPASSERT(SIZE(freq, 1) == SIZE(eigen_vec, 2))
     424           81 :          ALLOCATE (my_list(SIZE(particles)))
     425              :          ! Either we have a list of the subset of mobile atoms,
     426              :          ! Or the eigenvectors must span the full space (all atoms)
     427           27 :          IF (PRESENT(list)) THEN
     428           57 :             my_list(:) = 0
     429           54 :             DO i = 1, SIZE(list)
     430           54 :                my_list(list(i)) = i
     431              :             END DO
     432              :          ELSE
     433           12 :             CPASSERT(SIZE(particles) == SIZE(eigen_vec, 1)/3)
     434          443 :             DO i = 1, SIZE(my_list)
     435          443 :                my_list(i) = i
     436              :             END DO
     437              :          END IF
     438           27 :          WRITE (iw, '(T2,A)') "[Molden Format]"
     439           27 :          WRITE (iw, '(T2,A)') "[Atoms] AU"
     440          500 :          DO i = 1, SIZE(particles)
     441              :             CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
     442          473 :                                  element_symbol=element_symbol)
     443          473 :             CALL get_ptable_info(element_symbol, number=z)
     444              : 
     445              :             WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
     446         1919 :                element_symbol, i, z, particles(i)%r(:)
     447              : 
     448              :          END DO
     449           27 :          WRITE (iw, '(T2,A)') "[FREQ]"
     450          159 :          DO i = 1, SIZE(freq, 1)
     451          159 :             IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
     452              :          END DO
     453           27 :          WRITE (iw, '(T2,A)') "[FR-COORD]"
     454          500 :          DO i = 1, SIZE(particles)
     455              :             CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
     456          473 :                                  element_symbol=element_symbol)
     457              :             WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
     458         1919 :                element_symbol, particles(i)%r(:)
     459              :          END DO
     460           27 :          WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
     461           27 :          l = 0
     462          159 :          DO i = 1, SIZE(eigen_vec, 2)
     463          159 :             IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
     464          132 :                l = l + 1
     465          132 :                WRITE (iw, '(T2,A,1X,I6)') "vibration", l
     466         4071 :                DO j = 1, SIZE(particles)
     467         4071 :                   IF (my_list(j) .NE. 0) THEN
     468         3927 :                      k = (my_list(j) - 1)*3
     469         3927 :                      WRITE (iw, '(T2,3(F12.6,3X))') eigen_vec(k + 1, i), eigen_vec(k + 2, i), eigen_vec(k + 3, i)
     470              :                   ELSE
     471           12 :                      WRITE (iw, '(T2,3(F12.6,3X))') 0.0_dp, 0.0_dp, 0.0_dp
     472              :                   END IF
     473              :                END DO
     474              :             END IF
     475              :          END DO
     476           27 :          IF (calc_intens) THEN
     477           18 :             fint = massunit
     478              :             ! intensity units are a.u./amu
     479           18 :             WRITE (iw, '(T2,A)') "[INT]"
     480          118 :             DO i = 1, SIZE(intensities)
     481          118 :                IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(3X,F18.6)') fint*intensities(i)**2
     482              :             END DO
     483              :          END IF
     484           27 :          DEALLOCATE (my_list)
     485              :       END IF
     486           54 :       CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
     487              : 
     488           54 :       CALL timestop(handle)
     489              : 
     490           54 :    END SUBROUTINE write_vibrations_molden
     491              : 
     492              : END MODULE molden_utils
        

Generated by: LCOV version 2.0-1