LCOV - code coverage report
Current view: top level - src/motion - vibrational_analysis.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 586 597 98.2 %
Date: 2024-04-20 06:29:22 Functions: 7 7 100.0 %

          Line data    Source code
       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 Module performing a vibrational analysis
      10             : !> \note
      11             : !>      Numerical accuracy for parallel runs:
      12             : !>       Each replica starts the SCF run from the one optimized
      13             : !>       in a previous run. It may happen then energies and derivatives
      14             : !>       of a serial run and a parallel run could be slightly different
      15             : !>       'cause of a different starting density matrix.
      16             : !>       Exact results are obtained using:
      17             : !>          EXTRAPOLATION USE_GUESS in QS section (Teo 08.2006)
      18             : !> \author Teodoro Laino 08.2006
      19             : ! **************************************************************************************************
      20             : MODULE vibrational_analysis
      21             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      22             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
      23             :                                               cp_blacs_env_release,&
      24             :                                               cp_blacs_env_type
      25             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      26             :                                               cp_fm_struct_release,&
      27             :                                               cp_fm_struct_type
      28             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      29             :                                               cp_fm_release,&
      30             :                                               cp_fm_set_all,&
      31             :                                               cp_fm_set_element,&
      32             :                                               cp_fm_type,&
      33             :                                               cp_fm_write_unformatted
      34             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      35             :                                               cp_logger_get_default_io_unit,&
      36             :                                               cp_logger_type
      37             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      38             :                                               cp_print_key_unit_nr
      39             :    USE cp_result_methods,               ONLY: get_results,&
      40             :                                               test_for_result
      41             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      42             :                                               cp_subsys_type
      43             :    USE f77_interface,                   ONLY: f_env_add_defaults,&
      44             :                                               f_env_rm_defaults,&
      45             :                                               f_env_type
      46             :    USE force_env_types,                 ONLY: force_env_get,&
      47             :                                               force_env_type
      48             :    USE global_types,                    ONLY: global_environment_type
      49             :    USE grrm_utils,                      ONLY: write_grrm
      50             :    USE header,                          ONLY: vib_header
      51             :    USE input_constants,                 ONLY: do_rep_blocked
      52             :    USE input_section_types,             ONLY: section_type,&
      53             :                                               section_vals_get,&
      54             :                                               section_vals_get_subs_vals,&
      55             :                                               section_vals_type,&
      56             :                                               section_vals_val_get
      57             :    USE kinds,                           ONLY: default_string_length,&
      58             :                                               dp
      59             :    USE mathconstants,                   ONLY: pi
      60             :    USE mathlib,                         ONLY: diamat_all
      61             :    USE message_passing,                 ONLY: mp_para_env_type
      62             :    USE mode_selective,                  ONLY: ms_vb_anal
      63             :    USE molden_utils,                    ONLY: write_vibrations_molden
      64             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      65             :    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
      66             :                                               get_molecule_kind,&
      67             :                                               molecule_kind_type
      68             :    USE motion_utils,                    ONLY: rot_ana,&
      69             :                                               thrs_motion
      70             :    USE particle_list_types,             ONLY: particle_list_type
      71             :    USE particle_methods,                ONLY: write_particle_matrix
      72             :    USE particle_types,                  ONLY: particle_type
      73             :    USE physcon,                         ONLY: &
      74             :         a_bohr, angstrom, bohr, boltzmann, debye, e_mass, h_bar, hertz, kelvin, kjmol, massunit, &
      75             :         n_avogadro, pascal, vibfac, wavenumbers
      76             :    USE replica_methods,                 ONLY: rep_env_calc_e_f,&
      77             :                                               rep_env_create
      78             :    USE replica_types,                   ONLY: rep_env_release,&
      79             :                                               replica_env_type
      80             :    USE scine_utils,                     ONLY: write_scine
      81             :    USE util,                            ONLY: sort
      82             : #include "../base/base_uses.f90"
      83             : 
      84             :    IMPLICIT NONE
      85             :    PRIVATE
      86             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'vibrational_analysis'
      87             :    LOGICAL, PARAMETER                   :: debug_this_module = .FALSE.
      88             : 
      89             :    PUBLIC :: vb_anal
      90             : 
      91             : CONTAINS
      92             : 
      93             : ! **************************************************************************************************
      94             : !> \brief Module performing a vibrational analysis
      95             : !> \param input ...
      96             : !> \param input_declaration ...
      97             : !> \param para_env ...
      98             : !> \param globenv ...
      99             : !> \author Teodoro Laino 08.2006
     100             : ! **************************************************************************************************
     101          54 :    SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
     102             :       TYPE(section_vals_type), POINTER                   :: input
     103             :       TYPE(section_type), POINTER                        :: input_declaration
     104             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     105             :       TYPE(global_environment_type), POINTER             :: globenv
     106             : 
     107             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'vb_anal'
     108             :       CHARACTER(LEN=1), DIMENSION(3), PARAMETER          :: lab = (/"X", "Y", "Z"/)
     109             : 
     110             :       CHARACTER(LEN=default_string_length)               :: description_d, description_p
     111             :       INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, iounit, ip1, ip2, iparticle1, &
     112             :          iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nRotTrM, nvib, &
     113             :          output_unit, output_unit_eig, prep, print_grrm, print_namd, print_scine, proc_dist_type
     114          54 :       INTEGER, DIMENSION(:), POINTER                     :: Clist, Mlist
     115             :       LOGICAL :: calc_intens, calc_thchdata, do_mode_tracking, intens_ir, intens_raman, &
     116             :          keep_rotations, row_force, something_frozen
     117             :       REAL(KIND=dp)                                      :: a1, a2, a3, conver, dummy, dx, &
     118             :                                                             inertia(3), minimum_energy, norm, &
     119             :                                                             tc_press, tc_temp, tmp
     120          54 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: H_eigval1, H_eigval2, HeigvalDfull, &
     121          54 :                                                             konst, mass, pos0, rmass
     122          54 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Hessian, Hint1, Hint2, Hint2Dfull, MatM
     123             :       REAL(KIND=dp), DIMENSION(3)                        :: D_deriv, d_print
     124             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: P_deriv, p_print
     125          54 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: depol_p, depol_u, depp, depu, din, &
     126          54 :                                                             intensities_d, intensities_p, pin
     127          54 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: D, Dfull, dip_deriv, RotTrM
     128          54 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: polar_deriv, tmp_dip
     129          54 :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: tmp_polar
     130             :       TYPE(cp_logger_type), POINTER                      :: logger
     131             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     132             :       TYPE(f_env_type), POINTER                          :: f_env
     133          54 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     134             :       TYPE(replica_env_type), POINTER                    :: rep_env
     135             :       TYPE(section_vals_type), POINTER                   :: force_env_section, &
     136             :                                                             mode_tracking_section, print_section, &
     137             :                                                             vib_section
     138             : 
     139          54 :       CALL timeset(routineN, handle)
     140          54 :       NULLIFY (D, RotTrM, logger, subsys, f_env, particles, rep_env, intensities_d, intensities_p, &
     141          54 :                vib_section, print_section, depol_p, depol_u)
     142          54 :       logger => cp_get_default_logger()
     143          54 :       vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS")
     144          54 :       print_section => section_vals_get_subs_vals(vib_section, "PRINT")
     145             :       output_unit = cp_print_key_unit_nr(logger, &
     146             :                                          print_section, &
     147             :                                          "PROGRAM_RUN_INFO", &
     148          54 :                                          extension=".vibLog")
     149          54 :       iounit = cp_logger_get_default_io_unit(logger)
     150             :       ! for output of cartesian frequencies and eigenvectors of the
     151             :       ! Hessian that can be used for initialisation of MD calculations
     152             :       output_unit_eig = cp_print_key_unit_nr(logger, &
     153             :                                              print_section, &
     154             :                                              "CARTESIAN_EIGS", &
     155             :                                              extension=".eig", &
     156             :                                              file_status="REPLACE", &
     157             :                                              file_action="WRITE", &
     158             :                                              do_backup=.TRUE., &
     159          54 :                                              file_form="UNFORMATTED")
     160             : 
     161          54 :       CALL section_vals_val_get(vib_section, "DX", r_val=dx)
     162          54 :       CALL section_vals_val_get(vib_section, "NPROC_REP", i_val=prep)
     163          54 :       CALL section_vals_val_get(vib_section, "PROC_DIST_TYPE", i_val=proc_dist_type)
     164          54 :       row_force = (proc_dist_type == do_rep_blocked)
     165          54 :       CALL section_vals_val_get(vib_section, "FULLY_PERIODIC", l_val=keep_rotations)
     166          54 :       CALL section_vals_val_get(vib_section, "INTENSITIES", l_val=calc_intens)
     167          54 :       CALL section_vals_val_get(vib_section, "THERMOCHEMISTRY", l_val=calc_thchdata)
     168          54 :       CALL section_vals_val_get(vib_section, "TC_TEMPERATURE", r_val=tc_temp)
     169          54 :       CALL section_vals_val_get(vib_section, "TC_PRESSURE", r_val=tc_press)
     170             : 
     171          54 :       tc_temp = tc_temp*kelvin
     172          54 :       tc_press = tc_press*pascal
     173             : 
     174          54 :       intens_ir = .FALSE.
     175          54 :       intens_raman = .FALSE.
     176             : 
     177          54 :       mode_tracking_section => section_vals_get_subs_vals(vib_section, "MODE_SELECTIVE")
     178          54 :       CALL section_vals_get(mode_tracking_section, explicit=do_mode_tracking)
     179          54 :       nrep = MAX(1, para_env%num_pe/prep)
     180          54 :       prep = para_env%num_pe/nrep
     181          54 :       iw = cp_print_key_unit_nr(logger, print_section, "BANNER", extension=".vibLog")
     182          54 :       CALL vib_header(iw, nrep, prep)
     183          54 :       CALL cp_print_key_finished_output(iw, logger, print_section, "BANNER")
     184             :       ! Just one force_env allowed
     185          54 :       force_env_section => section_vals_get_subs_vals(input, "FORCE_EVAL")
     186             :       ! Create Replica Environments
     187             :       CALL rep_env_create(rep_env, para_env=para_env, input=input, &
     188          54 :                           input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=row_force)
     189          54 :       IF (ASSOCIATED(rep_env)) THEN
     190          54 :          CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env)
     191          54 :          CALL force_env_get(f_env%force_env, subsys=subsys)
     192          54 :          particles => subsys%particles%els
     193             :          ! Decide which kind of Vibrational Analysis to perform
     194          54 :          IF (do_mode_tracking) THEN
     195             :             CALL ms_vb_anal(input, rep_env, para_env, globenv, particles, &
     196          24 :                             nrep, calc_intens, dx, output_unit, logger)
     197          24 :             CALL f_env_rm_defaults(f_env, ierr)
     198             :          ELSE
     199          30 :             CALL get_moving_atoms(force_env=f_env%force_env, Ilist=Mlist)
     200          30 :             something_frozen = SIZE(particles) .NE. SIZE(Mlist)
     201          30 :             natoms = SIZE(Mlist)
     202          30 :             ncoord = natoms*3
     203          90 :             ALLOCATE (Clist(ncoord))
     204          90 :             ALLOCATE (mass(natoms))
     205          90 :             ALLOCATE (pos0(ncoord))
     206         120 :             ALLOCATE (Hessian(ncoord, ncoord))
     207          30 :             IF (calc_intens) THEN
     208          16 :                description_d = '[DIPOLE]'
     209          64 :                ALLOCATE (tmp_dip(ncoord, 3, 2))
     210         936 :                tmp_dip = 0._dp
     211          16 :                description_p = '[POLAR]'
     212          80 :                ALLOCATE (tmp_polar(ncoord, 3, 3, 2))
     213        2808 :                tmp_polar = 0._dp
     214             :             END IF
     215         264 :             Clist = 0
     216         108 :             DO i = 1, natoms
     217          78 :                imap = Mlist(i)
     218          78 :                Clist((i - 1)*3 + 1) = (imap - 1)*3 + 1
     219          78 :                Clist((i - 1)*3 + 2) = (imap - 1)*3 + 2
     220          78 :                Clist((i - 1)*3 + 3) = (imap - 1)*3 + 3
     221          78 :                mass(i) = particles(imap)%atomic_kind%mass
     222          78 :                CPASSERT(mass(i) > 0.0_dp)
     223          78 :                mass(i) = SQRT(mass(i))
     224          78 :                pos0((i - 1)*3 + 1) = particles(imap)%r(1)
     225          78 :                pos0((i - 1)*3 + 2) = particles(imap)%r(2)
     226         108 :                pos0((i - 1)*3 + 3) = particles(imap)%r(3)
     227             :             END DO
     228             :             !
     229             :             ! Determine the principal axes of inertia.
     230             :             ! Generation of coordinates in the rotating and translating frame
     231             :             !
     232          30 :             IF (something_frozen) THEN
     233           4 :                nRotTrM = 0
     234          12 :                ALLOCATE (RotTrM(natoms*3, nRotTrM))
     235             :             ELSE
     236             :                CALL rot_ana(particles, RotTrM, nRotTrM, print_section, &
     237          26 :                             keep_rotations, mass_weighted=.TRUE., natoms=natoms, inertia=inertia)
     238             :             END IF
     239             :             ! Generate the suitable rototranslating basis set
     240          30 :             nvib = 3*natoms - nRotTrM
     241             :             IF (.FALSE.) THEN !option full in build_D_matrix, at the moment not enabled
     242             :                !but dimensions of D must be adjusted in this case
     243             :                ALLOCATE (D(3*natoms, 3*natoms))
     244             :             ELSE
     245         150 :                ALLOCATE (D(3*natoms, nvib))
     246             :             END IF
     247             :             CALL build_D_matrix(RotTrM, nRotTrM, D, full=.FALSE., &
     248          30 :                                 natoms=natoms)
     249             :             !
     250             :             ! Loop on atoms and coordinates
     251             :             !
     252        2190 :             Hessian = HUGE(0.0_dp)
     253          30 :             IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "VIB| Vibrational Analysis Info"
     254         188 :             DO icoordp = 1, ncoord, nrep
     255         158 :                icoord = icoordp - 1
     256         408 :                DO j = 1, nrep
     257        2308 :                   DO i = 1, ncoord
     258        2058 :                      imap = Clist(i)
     259        2308 :                      rep_env%r(imap, j) = pos0(i)
     260             :                   END DO
     261         408 :                   IF (icoord + j <= ncoord) THEN
     262         234 :                      imap = Clist(icoord + j)
     263         234 :                      rep_env%r(imap, j) = rep_env%r(imap, j) + Dx
     264             :                   END IF
     265             :                END DO
     266         158 :                CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
     267             : 
     268         438 :                DO j = 1, nrep
     269         250 :                   IF (calc_intens) THEN
     270         140 :                      IF (icoord + j <= ncoord) THEN
     271         132 :                         IF (test_for_result(results=rep_env%results(j)%results, &
     272             :                                             description=description_d)) THEN
     273             :                            CALL get_results(results=rep_env%results(j)%results, &
     274             :                                             description=description_d, &
     275         132 :                                             n_rep=nres)
     276             :                            CALL get_results(results=rep_env%results(j)%results, &
     277             :                                             description=description_d, &
     278             :                                             values=tmp_dip(icoord + j, :, 1), &
     279         132 :                                             nval=nres)
     280         132 :                            intens_ir = .TRUE.
     281         528 :                            d_print(:) = tmp_dip(icoord + j, :, 1)
     282             :                         END IF
     283         132 :                         IF (test_for_result(results=rep_env%results(j)%results, &
     284             :                                             description=description_p)) THEN
     285             :                            CALL get_results(results=rep_env%results(j)%results, &
     286             :                                             description=description_p, &
     287          12 :                                             n_rep=nres)
     288             :                            CALL get_results(results=rep_env%results(j)%results, &
     289             :                                             description=description_p, &
     290             :                                             values=tmp_polar(icoord + j, :, :, 1), &
     291          12 :                                             nval=nres)
     292          12 :                            intens_raman = .TRUE.
     293         156 :                            p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
     294             :                         END IF
     295             :                      END IF
     296             :                   END IF
     297         408 :                   IF (icoord + j <= ncoord) THEN
     298        2160 :                      DO i = 1, ncoord
     299        1926 :                         imap = Clist(i)
     300        2160 :                         Hessian(i, icoord + j) = rep_env%f(imap, j)
     301             :                      END DO
     302         234 :                      imap = Clist(icoord + j)
     303             :                      ! Dump Info
     304         234 :                      IF (output_unit > 0) THEN
     305          75 :                         iparticle1 = imap/3
     306          75 :                         IF (MOD(imap, 3) /= 0) iparticle1 = iparticle1 + 1
     307             :                         WRITE (output_unit, '(T2,A,I5,A,I5,3A)') &
     308          75 :                            "VIB| REPLICA Nr.", j, "- Energy and Forces for particle:", &
     309          75 :                            iparticle1, "  coordinate: ", lab(imap - (iparticle1 - 1)*3), &
     310         150 :                            " + D"//TRIM(lab(imap - (iparticle1 - 1)*3))
     311             :                         WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
     312          75 :                            "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
     313          75 :                         WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
     314         276 :                         DO i = 1, natoms
     315         201 :                            imap = Mlist(i)
     316             :                            WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
     317         201 :                               particles(imap)%atomic_kind%name, &
     318         477 :                               rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
     319             :                         END DO
     320          75 :                         IF (intens_ir) THEN
     321          33 :                            WRITE (output_unit, '(T3,A)') 'Dipole moment [Debye]'
     322             :                            WRITE (output_unit, '(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
     323          33 :                               'X=', d_print(1)*debye, 'Y=', d_print(2)*debye, 'Z=', d_print(3)*debye, &
     324         165 :                               'Total=', SQRT(SUM(d_print(1:3)**2))*debye
     325             :                         END IF
     326          75 :                         IF (intens_raman) THEN
     327             :                            WRITE (output_unit, '(T2,A)') &
     328           6 :                               'POLAR| Polarizability tensor [a.u.]'
     329             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
     330           6 :                               'POLAR| xx,yy,zz', p_print(1, 1), p_print(2, 2), p_print(3, 3)
     331             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
     332           6 :                               'POLAR| xy,xz,yz', p_print(1, 2), p_print(1, 3), p_print(2, 3)
     333             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12),/)') &
     334           6 :                               'POLAR| yx,zx,zy', p_print(2, 1), p_print(3, 1), p_print(3, 2)
     335             :                         END IF
     336             :                      END IF
     337             :                   END IF
     338             :                END DO
     339             :             END DO
     340         188 :             DO icoordm = 1, ncoord, nrep
     341         158 :                icoord = icoordm - 1
     342         408 :                DO j = 1, nrep
     343        2308 :                   DO i = 1, ncoord
     344        2058 :                      imap = Clist(i)
     345        2308 :                      rep_env%r(imap, j) = pos0(i)
     346             :                   END DO
     347         408 :                   IF (icoord + j <= ncoord) THEN
     348         234 :                      imap = Clist(icoord + j)
     349         234 :                      rep_env%r(imap, j) = rep_env%r(imap, j) - Dx
     350             :                   END IF
     351             :                END DO
     352         158 :                CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
     353             : 
     354         438 :                DO j = 1, nrep
     355         250 :                   IF (calc_intens) THEN
     356         140 :                      IF (icoord + j <= ncoord) THEN
     357         132 :                         k = (icoord + j + 2)/3
     358         132 :                         IF (test_for_result(results=rep_env%results(j)%results, &
     359             :                                             description=description_d)) THEN
     360             :                            CALL get_results(results=rep_env%results(j)%results, &
     361             :                                             description=description_d, &
     362         132 :                                             n_rep=nres)
     363             :                            CALL get_results(results=rep_env%results(j)%results, &
     364             :                                             description=description_d, &
     365             :                                             values=tmp_dip(icoord + j, :, 2), &
     366         132 :                                             nval=nres)
     367             :                            tmp_dip(icoord + j, :, 1) = (tmp_dip(icoord + j, :, 1) - &
     368         528 :                                                         tmp_dip(icoord + j, :, 2))/(2.0_dp*Dx*mass(k))
     369         528 :                            d_print(:) = tmp_dip(icoord + j, :, 1)
     370             :                         END IF
     371         132 :                         IF (test_for_result(results=rep_env%results(j)%results, &
     372             :                                             description=description_p)) THEN
     373             :                            CALL get_results(results=rep_env%results(j)%results, &
     374             :                                             description=description_p, &
     375          12 :                                             n_rep=nres)
     376             :                            CALL get_results(results=rep_env%results(j)%results, &
     377             :                                             description=description_p, &
     378             :                                             values=tmp_polar(icoord + j, :, :, 2), &
     379          12 :                                             nval=nres)
     380             :                            tmp_polar(icoord + j, :, :, 1) = (tmp_polar(icoord + j, :, :, 1) - &
     381         156 :                                                              tmp_polar(icoord + j, :, :, 2))/(2.0_dp*Dx*mass(k))
     382         156 :                            p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
     383             :                         END IF
     384             :                      END IF
     385             :                   END IF
     386         408 :                   IF (icoord + j <= ncoord) THEN
     387         234 :                      imap = Clist(icoord + j)
     388         234 :                      iparticle1 = imap/3
     389         234 :                      IF (MOD(imap, 3) /= 0) iparticle1 = iparticle1 + 1
     390         234 :                      ip1 = (icoord + j)/3
     391         234 :                      IF (MOD(icoord + j, 3) /= 0) ip1 = ip1 + 1
     392             :                      ! Dump Info
     393         234 :                      IF (output_unit > 0) THEN
     394             :                         WRITE (output_unit, '(T2,A,I5,A,I5,3A)') &
     395          75 :                            "VIB| REPLICA Nr.", j, "- Energy and Forces for particle:", &
     396          75 :                            iparticle1, "  coordinate: ", lab(imap - (iparticle1 - 1)*3), &
     397         150 :                            " - D"//TRIM(lab(imap - (iparticle1 - 1)*3))
     398             :                         WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
     399          75 :                            "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
     400          75 :                         WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
     401         276 :                         DO i = 1, natoms
     402         201 :                            imap = Mlist(i)
     403             :                            WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
     404         201 :                               particles(imap)%atomic_kind%name, &
     405         477 :                               rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
     406             :                         END DO
     407          75 :                         IF (intens_ir) THEN
     408          33 :                            WRITE (output_unit, '(T3,A)') 'Dipole moment [Debye]'
     409             :                            WRITE (output_unit, '(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
     410          33 :                               'X=', d_print(1)*debye, 'Y=', d_print(2)*debye, 'Z=', d_print(3)*debye, &
     411         165 :                               'Total=', SQRT(SUM(d_print(1:3)**2))*debye
     412             :                         END IF
     413          75 :                         IF (intens_raman) THEN
     414             :                            WRITE (output_unit, '(T2,A)') &
     415           6 :                               'POLAR| Polarizability tensor [a.u.]'
     416             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
     417           6 :                               'POLAR| xx,yy,zz', p_print(1, 1), p_print(2, 2), p_print(3, 3)
     418             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
     419           6 :                               'POLAR| xy,xz,yz', p_print(1, 2), p_print(1, 3), p_print(2, 3)
     420             :                            WRITE (output_unit, '(T2,A,T24,3(1X,F18.12),/)') &
     421           6 :                               'POLAR| yx,zx,zy', p_print(2, 1), p_print(3, 1), p_print(3, 2)
     422             :                         END IF
     423             :                      END IF
     424        2160 :                      DO iseq = 1, ncoord
     425        1926 :                         imap = Clist(iseq)
     426        1926 :                         iparticle2 = imap/3
     427        1926 :                         IF (MOD(imap, 3) /= 0) iparticle2 = iparticle2 + 1
     428        1926 :                         ip2 = iseq/3
     429        1926 :                         IF (MOD(iseq, 3) /= 0) ip2 = ip2 + 1
     430        1926 :                         tmp = Hessian(iseq, icoord + j) - rep_env%f(imap, j)
     431        1926 :                         tmp = -tmp/(2.0_dp*Dx*mass(ip1)*mass(ip2))*1E6_dp
     432             :                         ! Mass weighted Hessian
     433        2160 :                         Hessian(iseq, icoord + j) = tmp
     434             :                      END DO
     435             :                   END IF
     436             :                END DO
     437             :             END DO
     438             : 
     439             :             ! restore original particle positions for output
     440         108 :             DO i = 1, natoms
     441          78 :                imap = Mlist(i)
     442         342 :                particles(imap)%r(1:3) = pos0((i - 1)*3 + 1:(i - 1)*3 + 3)
     443             :             END DO
     444          82 :             DO j = 1, nrep
     445         484 :                DO i = 1, ncoord
     446         402 :                   imap = Clist(i)
     447         454 :                   rep_env%r(imap, j) = pos0(i)
     448             :                END DO
     449             :             END DO
     450          30 :             CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
     451          30 :             j = 1
     452          30 :             minimum_energy = rep_env%f(rep_env%ndim + 1, j)
     453          30 :             IF (output_unit > 0) THEN
     454             :                WRITE (output_unit, '(T2,A)') &
     455          10 :                   "VIB| ", " Minimum Structure - Energy and Forces:"
     456             :                WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
     457          10 :                   "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
     458          10 :                WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
     459          35 :                DO i = 1, natoms
     460          25 :                   imap = Mlist(i)
     461             :                   WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
     462          25 :                      particles(imap)%atomic_kind%name, &
     463          60 :                      rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
     464             :                END DO
     465             :             END IF
     466             : 
     467             :             ! Dump Info
     468          30 :             IF (output_unit > 0) THEN
     469          10 :                WRITE (output_unit, '(T2,A)') "VIB| Hessian in cartesian coordinates"
     470             :                CALL write_particle_matrix(Hessian, particles, output_unit, el_per_part=3, &
     471          10 :                                           Ilist=Mlist)
     472             :             END IF
     473             : 
     474          30 :             CALL write_va_hessian(vib_section, para_env, ncoord, globenv, Hessian, logger)
     475             : 
     476             :             ! Enforce symmetry in the Hessian
     477         264 :             DO i = 1, ncoord
     478        1344 :                DO j = i, ncoord
     479             :                   ! Take the upper diagonal part
     480        1314 :                   Hessian(j, i) = Hessian(i, j)
     481             :                END DO
     482             :             END DO
     483             :             !
     484             :             ! Print GRMM interface file
     485             :             print_grrm = cp_print_key_unit_nr(logger, force_env_section, "PRINT%GRRM", &
     486          30 :                                               file_position="REWIND", extension=".rrm")
     487          30 :             IF (print_grrm > 0) THEN
     488           7 :                DO i = 1, natoms
     489           5 :                   imap = Mlist(i)
     490          37 :                   particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
     491             :                END DO
     492          12 :                ALLOCATE (Hint1(ncoord, ncoord), rmass(ncoord))
     493           7 :                DO i = 1, natoms
     494           5 :                   imap = Mlist(i)
     495          22 :                   rmass(3*(imap - 1) + 1:3*(imap - 1) + 3) = mass(imap)
     496             :                END DO
     497          17 :                DO i = 1, ncoord
     498         134 :                   DO j = 1, ncoord
     499         132 :                      Hint1(j, i) = Hessian(j, i)*rmass(i)*rmass(j)*1.0E-6_dp
     500             :                   END DO
     501             :                END DO
     502           2 :                nfrozen = SIZE(particles) - natoms
     503             :                CALL write_grrm(print_grrm, f_env%force_env, particles, minimum_energy, &
     504           2 :                                hessian=Hint1, fixed_atoms=nfrozen)
     505           2 :                DEALLOCATE (Hint1, rmass)
     506             :             END IF
     507          30 :             CALL cp_print_key_finished_output(print_grrm, logger, force_env_section, "PRINT%GRRM")
     508             :             !
     509             :             ! Print SCINE interface file
     510             :             print_scine = cp_print_key_unit_nr(logger, force_env_section, "PRINT%SCINE", &
     511          30 :                                                file_position="REWIND", extension=".scine")
     512          30 :             IF (print_scine > 0) THEN
     513           4 :                DO i = 1, natoms
     514           3 :                   imap = Mlist(i)
     515          22 :                   particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
     516             :                END DO
     517           1 :                nfrozen = SIZE(particles) - natoms
     518           1 :                CPASSERT(nfrozen == 0)
     519           1 :                CALL write_scine(print_scine, f_env%force_env, particles, minimum_energy, hessian=Hessian)
     520             :             END IF
     521          30 :             CALL cp_print_key_finished_output(print_scine, logger, force_env_section, "PRINT%SCINE")
     522             :             !
     523             :             ! Print NEWTONX interface file
     524             :             print_namd = cp_print_key_unit_nr(logger, print_section, "NAMD_PRINT", &
     525             :                                               extension=".eig", file_status="REPLACE", &
     526             :                                               file_action="WRITE", do_backup=.TRUE., &
     527          30 :                                               file_form="UNFORMATTED")
     528          30 :             IF (print_namd > 0) THEN
     529             :                ! NewtonX requires normalized Cartesian frequencies and eigenvectors
     530             :                ! in full matrix format (ncoord x ncoord)
     531             :                NULLIFY (Dfull)
     532           4 :                ALLOCATE (Dfull(ncoord, ncoord))
     533           4 :                ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2)))
     534           3 :                ALLOCATE (HeigvalDfull(SIZE(Dfull, 2)))
     535           4 :                ALLOCATE (MatM(ncoord, ncoord))
     536           3 :                ALLOCATE (rmass(SIZE(Dfull, 2)))
     537          91 :                Dfull = 0.0_dp
     538             :                ! Dfull in dimension of degrees of freedom
     539           1 :                CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms)
     540             :                ! TEST MatM = MATMUL(TRANSPOSE(Dfull),Dfull)= 1
     541             :                ! Hessian in MWC -> Hessian in INT (Hint2Dfull)
     542        6196 :                Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull))
     543             :                ! Heig = L^T Hint2Dfull L
     544           1 :                CALL diamat_all(Hint2Dfull, HeigvalDfull)
     545             :                ! TEST  MatM = MATMUL(TRANSPOSE(Hint2Dfull),Hint2Dfull) = 1
     546             :                ! TEST MatM=MATMUL(TRANSPOSE(MATMUL(Dfull,Hint2Dfull)),MATMUL(Dfull,Hint2Dfull)) = 1
     547          91 :                MatM = 0.0_dp
     548           4 :                DO i = 1, natoms
     549          13 :                   DO j = 1, 3
     550          12 :                      MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i) ! mass is sqrt(mass)
     551             :                   END DO
     552             :                END DO
     553             :                ! Dfull = Cartesian displacements of the normal modes
     554        6196 :                Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull))  !Dfull=D L / sqrt(m)
     555          10 :                DO i = 1, ncoord
     556             :                   ! Renormalize displacements
     557          90 :                   norm = 1.0_dp/SUM(Dfull(:, i)*Dfull(:, i))
     558           9 :                   rmass(i) = norm/massunit
     559          91 :                   Dfull(:, i) = SQRT(norm)*(Dfull(:, i))
     560             :                END DO
     561           1 :                CALL write_eigs_unformatted(print_namd, ncoord, HeigvalDfull, Dfull)
     562           1 :                DEALLOCATE (HeigvalDfull)
     563           1 :                DEALLOCATE (Hint2Dfull)
     564           1 :                DEALLOCATE (Dfull)
     565           1 :                DEALLOCATE (MatM)
     566           1 :                DEALLOCATE (rmass)
     567             :             END IF !print_namd
     568             :             !
     569          30 :             nvib = ncoord - nRotTrM
     570          90 :             ALLOCATE (H_eigval1(ncoord))
     571          90 :             ALLOCATE (H_eigval2(SIZE(D, 2)))
     572         120 :             ALLOCATE (Hint1(ncoord, ncoord))
     573         120 :             ALLOCATE (Hint2(SIZE(D, 2), SIZE(D, 2)))
     574          90 :             ALLOCATE (rmass(SIZE(D, 2)))
     575          90 :             ALLOCATE (konst(SIZE(D, 2)))
     576          30 :             IF (calc_intens) THEN
     577          48 :                ALLOCATE (dip_deriv(3, SIZE(D, 2)))
     578         216 :                dip_deriv = 0.0_dp
     579          48 :                ALLOCATE (polar_deriv(3, 3, SIZE(D, 2)))
     580         666 :                polar_deriv = 0.0_dp
     581             :             END IF
     582          90 :             ALLOCATE (intensities_d(SIZE(D, 2)))
     583          90 :             ALLOCATE (intensities_p(SIZE(D, 2)))
     584          90 :             ALLOCATE (depol_p(SIZE(D, 2)))
     585          90 :             ALLOCATE (depol_u(SIZE(D, 2)))
     586         120 :             intensities_d = 0._dp
     587         120 :             intensities_p = 0._dp
     588         120 :             depol_p = 0._dp
     589         120 :             depol_u = 0._dp
     590        2190 :             Hint1(:, :) = Hessian
     591          30 :             CALL diamat_all(Hint1, H_eigval1)
     592          30 :             IF (output_unit > 0) THEN
     593             :                WRITE (output_unit, '(T2,"VIB| Cartesian Low frequencies ---",4G12.5)') &
     594          10 :                   (H_eigval1(i), i=1, MIN(9, ncoord))
     595          10 :                WRITE (output_unit, '(T2,A)') "VIB| Eigenvectors before removal of rotations and translations"
     596             :                CALL write_particle_matrix(Hint1, particles, output_unit, el_per_part=3, &
     597          10 :                                           Ilist=Mlist)
     598             :             END IF
     599             :             ! write frequencies and eigenvectors to cartesian eig file
     600          30 :             IF (output_unit_eig > 0) THEN
     601          15 :                CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Hint1)
     602             :             END IF
     603          30 :             IF (nvib /= 0) THEN
     604       37284 :                Hint2(:, :) = MATMUL(TRANSPOSE(D), MATMUL(Hessian, D))
     605          30 :                IF (calc_intens) THEN
     606          64 :                   DO i = 1, 3
     607        5854 :                      dip_deriv(i, :) = MATMUL(tmp_dip(:, i, 1), D)
     608             :                   END DO
     609          64 :                   DO i = 1, 3
     610         208 :                      DO j = 1, 3
     611       17562 :                         polar_deriv(i, j, :) = MATMUL(tmp_polar(:, i, j, 1), D)
     612             :                      END DO
     613             :                   END DO
     614             :                END IF
     615          30 :                CALL diamat_all(Hint2, H_eigval2)
     616          30 :                IF (output_unit > 0) THEN
     617          10 :                   WRITE (output_unit, '(T2,"VIB| Frequencies after removal of the rotations and translations")')
     618             :                   ! Frequency at the moment are in a.u.
     619          10 :                   WRITE (output_unit, '(T2,"VIB| Internal  Low frequencies ---",4G12.5)') H_eigval2
     620             :                END IF
     621        2190 :                Hessian = 0.0_dp
     622         108 :                DO i = 1, natoms
     623         342 :                   DO j = 1, 3
     624         312 :                      Hessian((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i)
     625             :                   END DO
     626             :                END DO
     627             :                ! Cartesian displacements of the normal modes
     628       49770 :                D = MATMUL(Hessian, MATMUL(D, Hint2))
     629         120 :                DO i = 1, nvib
     630         810 :                   norm = 1.0_dp/SUM(D(:, i)*D(:, i))
     631             :                   ! Reduced Masess
     632          90 :                   rmass(i) = norm/massunit
     633             :                   ! Renormalize displacements and convert in Angstrom
     634         810 :                   D(:, i) = SQRT(norm)*D(:, i)
     635             :                   ! Force constants
     636          90 :                   konst(i) = SIGN(1.0_dp, H_eigval2(i))*2.0_dp*pi**2*(ABS(H_eigval2(i))/massunit)**2*rmass(i)
     637          90 :                   IF (calc_intens) THEN
     638          50 :                      D_deriv = 0._dp
     639         232 :                      DO j = 1, nvib
     640         778 :                         D_deriv(:) = D_deriv(:) + dip_deriv(:, j)*Hint2(j, i)
     641             :                      END DO
     642         200 :                      intensities_d(i) = SQRT(DOT_PRODUCT(D_deriv, D_deriv))
     643          50 :                      P_deriv = 0._dp
     644         232 :                      DO j = 1, nvib
     645             :                         ! P_deriv has units bohr^2/sqrt(a.u.)
     646        2416 :                         P_deriv(:, :) = P_deriv(:, :) + polar_deriv(:, :, j)*Hint2(j, i)
     647             :                      END DO
     648             :                      ! P_deriv now has units A^2/sqrt(amu)
     649             :                      conver = angstrom**2*SQRT(massunit)
     650         650 :                      P_deriv(:, :) = P_deriv(:, :)*conver
     651             :                      ! this is wron, just for testing
     652          50 :                      a1 = (P_deriv(1, 1) + P_deriv(2, 2) + P_deriv(3, 3))/3.0_dp
     653             :                      a2 = (P_deriv(1, 1) - P_deriv(2, 2))**2 + &
     654             :                           (P_deriv(2, 2) - P_deriv(3, 3))**2 + &
     655          50 :                           (P_deriv(3, 3) - P_deriv(1, 1))**2
     656          50 :                      a3 = (P_deriv(1, 2)**2 + P_deriv(2, 3)**2 + P_deriv(3, 1)**2)
     657          50 :                      intensities_p(i) = 45.0_dp*a1*a1 + 7.0_dp/2.0_dp*(a2 + 6.0_dp*a3)
     658             :                      ! to avoid division by zero:
     659          50 :                      dummy = 45.0_dp*a1*a1 + 4.0_dp/2.0_dp*(a2 + 6.0_dp*a3)
     660          50 :                      IF (dummy > 5.E-7_dp) THEN
     661             :                         ! depolarization of plane polarized incident light
     662             :                         depol_p(i) = 3.0_dp/2.0_dp*(a2 + 6.0_dp*a3)/(45.0_dp*a1*a1 + &
     663           2 :                                                                      4.0_dp/2.0_dp*(a2 + 6.0_dp*a3))
     664             :                         ! depolarization of unpolarized (natural) incident light
     665             :                         depol_u(i) = 6.0_dp/2.0_dp*(a2 + 6.0_dp*a3)/(45.0_dp*a1*a1 + &
     666           2 :                                                                      7.0_dp/2.0_dp*(a2 + 6.0_dp*a3))
     667             :                      ELSE
     668          48 :                         depol_p(i) = -1.0_dp
     669          48 :                         depol_u(i) = -1.0_dp
     670             :                      END IF
     671             :                   END IF
     672             :                   ! Convert frequencies to cm^-1
     673         120 :                   H_eigval2(i) = SIGN(1.0_dp, H_eigval2(i))*SQRT(ABS(H_eigval2(i))*massunit)*vibfac/1000.0_dp
     674             :                END DO
     675          30 :                IF (calc_intens) THEN
     676          16 :                   IF (iounit > 0) THEN
     677           8 :                      IF (.NOT. intens_ir) THEN
     678           0 :                         WRITE (iounit, '(T2,"VIB| No IR intensities available. Check input")')
     679             :                      END IF
     680           8 :                      IF (.NOT. intens_raman) THEN
     681           7 :                         WRITE (iounit, '(T2,"VIB| No Raman intensities available. Check input")')
     682             :                      END IF
     683             :                   END IF
     684             :                END IF
     685             :                ! Dump Info
     686          30 :                iw = cp_logger_get_default_io_unit(logger)
     687          30 :                IF (iw > 0) THEN
     688          15 :                   NULLIFY (din, pin, depp, depu)
     689          15 :                   IF (intens_ir) din => intensities_d
     690          15 :                   IF (intens_raman) pin => intensities_p
     691          15 :                   IF (intens_raman) depp => depol_p
     692          15 :                   IF (intens_raman) depu => depol_u
     693          15 :                   CALL vib_out(iw, nvib, D, konst, rmass, H_eigval2, particles, Mlist, din, pin, depp, depu)
     694             :                END IF
     695          30 :                IF (.NOT. something_frozen .AND. calc_thchdata) THEN
     696           2 :                   CALL get_thch_values(H_eigval2, iw, mass, nvib, inertia, 1, minimum_energy, tc_temp, tc_press)
     697             :                END IF
     698             :                CALL write_vibrations_molden(input, particles, H_eigval2, D, intensities_d, calc_intens, &
     699          30 :                                             dump_only_positive=.FALSE., logger=logger, list=Mlist)
     700             :             ELSE
     701           0 :                IF (output_unit > 0) THEN
     702           0 :                   WRITE (output_unit, '(T2,"VIB| No further vibrational info. Detected a single atom")')
     703             :                END IF
     704             :             END IF
     705             :             ! Deallocate working arrays
     706          30 :             DEALLOCATE (RotTrM)
     707          30 :             DEALLOCATE (Clist)
     708          30 :             DEALLOCATE (Mlist)
     709          30 :             DEALLOCATE (H_eigval1)
     710          30 :             DEALLOCATE (H_eigval2)
     711          30 :             DEALLOCATE (Hint1)
     712          30 :             DEALLOCATE (Hint2)
     713          30 :             DEALLOCATE (rmass)
     714          30 :             DEALLOCATE (konst)
     715          30 :             DEALLOCATE (mass)
     716          30 :             DEALLOCATE (pos0)
     717          30 :             DEALLOCATE (D)
     718          30 :             DEALLOCATE (Hessian)
     719          30 :             IF (calc_intens) THEN
     720          16 :                DEALLOCATE (dip_deriv)
     721          16 :                DEALLOCATE (polar_deriv)
     722          16 :                DEALLOCATE (tmp_dip)
     723          16 :                DEALLOCATE (tmp_polar)
     724             :             END IF
     725          30 :             DEALLOCATE (intensities_d)
     726          30 :             DEALLOCATE (intensities_p)
     727          30 :             DEALLOCATE (depol_p)
     728          30 :             DEALLOCATE (depol_u)
     729          30 :             CALL f_env_rm_defaults(f_env, ierr)
     730             :          END IF
     731             :       END IF
     732          54 :       CALL cp_print_key_finished_output(output_unit, logger, print_section, "PROGRAM_RUN_INFO")
     733          54 :       CALL cp_print_key_finished_output(output_unit_eig, logger, print_section, "CARTESIAN_EIGS")
     734          54 :       CALL rep_env_release(rep_env)
     735          54 :       CALL timestop(handle)
     736         108 :    END SUBROUTINE vb_anal
     737             : 
     738             : ! **************************************************************************************************
     739             : !> \brief give back a list of moving atoms
     740             : !> \param force_env ...
     741             : !> \param Ilist ...
     742             : !> \author Teodoro Laino 08.2006
     743             : ! **************************************************************************************************
     744          30 :    SUBROUTINE get_moving_atoms(force_env, Ilist)
     745             :       TYPE(force_env_type), POINTER                      :: force_env
     746             :       INTEGER, DIMENSION(:), POINTER                     :: Ilist
     747             : 
     748             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_moving_atoms'
     749             : 
     750             :       INTEGER                                            :: handle, i, ii, ikind, j, ndim, &
     751             :                                                             nfixed_atoms, nfixed_atoms_total, nkind
     752          30 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ifixd_list, work
     753             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     754          30 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     755             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     756          30 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     757             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     758             :       TYPE(particle_list_type), POINTER                  :: particles
     759          30 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     760             : 
     761          30 :       CALL timeset(routineN, handle)
     762          30 :       CALL force_env_get(force_env=force_env, subsys=subsys)
     763             : 
     764             :       CALL cp_subsys_get(subsys=subsys, particles=particles, &
     765          30 :                          molecule_kinds=molecule_kinds)
     766             : 
     767          30 :       nkind = molecule_kinds%n_els
     768          30 :       molecule_kind_set => molecule_kinds%els
     769          30 :       particle_set => particles%els
     770             : 
     771             :       ! Count the number of fixed atoms
     772          30 :       nfixed_atoms_total = 0
     773         102 :       DO ikind = 1, nkind
     774          72 :          molecule_kind => molecule_kind_set(ikind)
     775          72 :          CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
     776         102 :          nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
     777             :       END DO
     778          30 :       ndim = SIZE(particle_set) - nfixed_atoms_total
     779          30 :       CPASSERT(ndim >= 0)
     780          90 :       ALLOCATE (Ilist(ndim))
     781             : 
     782          30 :       IF (nfixed_atoms_total /= 0) THEN
     783          12 :          ALLOCATE (ifixd_list(nfixed_atoms_total))
     784          12 :          ALLOCATE (work(nfixed_atoms_total))
     785           4 :          nfixed_atoms_total = 0
     786          12 :          DO ikind = 1, nkind
     787           8 :             molecule_kind => molecule_kind_set(ikind)
     788           8 :             CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
     789          12 :             IF (ASSOCIATED(fixd_list)) THEN
     790          14 :                DO ii = 1, SIZE(fixd_list)
     791          14 :                   IF (.NOT. fixd_list(ii)%restraint%active) THEN
     792           6 :                      nfixed_atoms_total = nfixed_atoms_total + 1
     793           6 :                      ifixd_list(nfixed_atoms_total) = fixd_list(ii)%fixd
     794             :                   END IF
     795             :                END DO
     796             :             END IF
     797             :          END DO
     798           4 :          CALL sort(ifixd_list, nfixed_atoms_total, work)
     799             : 
     800           4 :          ndim = 0
     801           4 :          j = 1
     802          14 :          Loop_count: DO i = 1, SIZE(particle_set)
     803          14 :             DO WHILE (i > ifixd_list(j))
     804           4 :                j = j + 1
     805          14 :                IF (j > nfixed_atoms_total) EXIT Loop_count
     806             :             END DO
     807          14 :             IF (i /= ifixd_list(j)) THEN
     808           4 :                ndim = ndim + 1
     809           4 :                Ilist(ndim) = i
     810             :             END IF
     811             :          END DO Loop_count
     812           4 :          DEALLOCATE (ifixd_list)
     813           4 :          DEALLOCATE (work)
     814             :       ELSE
     815             :          i = 1
     816             :          ndim = 0
     817             :       END IF
     818         104 :       DO j = i, SIZE(particle_set)
     819          74 :          ndim = ndim + 1
     820         104 :          Ilist(ndim) = j
     821             :       END DO
     822          30 :       CALL timestop(handle)
     823             : 
     824          30 :    END SUBROUTINE get_moving_atoms
     825             : 
     826             : ! **************************************************************************************************
     827             : !> \brief Dumps results of the vibrational analysis
     828             : !> \param iw ...
     829             : !> \param nvib ...
     830             : !> \param D ...
     831             : !> \param k ...
     832             : !> \param m ...
     833             : !> \param freq ...
     834             : !> \param particles ...
     835             : !> \param Mlist ...
     836             : !> \param intensities_d ...
     837             : !> \param intensities_p ...
     838             : !> \param depol_p ...
     839             : !> \param depol_u ...
     840             : !> \author Teodoro Laino 08.2006
     841             : ! **************************************************************************************************
     842          15 :    SUBROUTINE vib_out(iw, nvib, D, k, m, freq, particles, Mlist, intensities_d, intensities_p, &
     843             :                       depol_p, depol_u)
     844             :       INTEGER, INTENT(IN)                                :: iw, nvib
     845             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: D
     846             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: k, m, freq
     847             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     848             :       INTEGER, DIMENSION(:), POINTER                     :: Mlist
     849             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: intensities_d, intensities_p, depol_p, &
     850             :                                                             depol_u
     851             : 
     852             :       CHARACTER(LEN=2)                                   :: element_symbol
     853             :       INTEGER                                            :: from, iatom, icol, j, jatom, katom, &
     854             :                                                             natom, to
     855             :       REAL(KIND=dp)                                      :: fint, pint
     856             : 
     857          15 :       fint = 42.255_dp*massunit*debye**2*bohr**2
     858          15 :       pint = 1.0_dp
     859          15 :       natom = SIZE(D, 1)
     860          15 :       WRITE (UNIT=iw, FMT="(/,T2,'VIB|',T30,'NORMAL MODES - CARTESIAN DISPLACEMENTS')")
     861          15 :       WRITE (UNIT=iw, FMT="(T2,'VIB|')")
     862          32 :       DO jatom = 1, nvib, 3
     863          17 :          from = jatom
     864          17 :          to = MIN(from + 2, nvib)
     865             :          WRITE (UNIT=iw, FMT="(T2,'VIB|',13X,3(8X,I5,8X))") &
     866          62 :             (icol, icol=from, to)
     867             :          WRITE (UNIT=iw, FMT="(T2,'VIB|Frequency (cm^-1)',3(1X,F12.6,8X))") &
     868          17 :             (freq(icol), icol=from, to)
     869          17 :          IF (ASSOCIATED(intensities_d)) THEN
     870             :             WRITE (UNIT=iw, FMT="(T2,'VIB|IR int (KM/Mole) ',3(1X,F12.6,8X))") &
     871          34 :                (fint*intensities_d(icol)**2, icol=from, to)
     872             :          END IF
     873          17 :          IF (ASSOCIATED(intensities_p)) THEN
     874             :             WRITE (UNIT=iw, FMT="(T2,'VIB|Raman (A^4/amu)  ',3(1X,F12.6,8X))") &
     875           2 :                (pint*intensities_p(icol), icol=from, to)
     876             :             WRITE (UNIT=iw, FMT="(T2,'VIB|Depol Ratio (P)  ',3(1X,F12.6,8X))") &
     877           1 :                (depol_p(icol), icol=from, to)
     878             :             WRITE (UNIT=iw, FMT="(T2,'VIB|Depol Ratio (U)  ',3(1X,F12.6,8X))") &
     879           1 :                (depol_u(icol), icol=from, to)
     880             :          END IF
     881             :          WRITE (UNIT=iw, FMT="(T2,'VIB|Red.Masses (a.u.)',3(1X,F12.6,8X))") &
     882          17 :             (m(icol), icol=from, to)
     883             :          WRITE (UNIT=iw, FMT="(T2,'VIB|Frc consts (a.u.)',3(1X,F12.6,8X))") &
     884          17 :             (k(icol), icol=from, to)
     885          17 :          WRITE (UNIT=iw, FMT="(T2,' ATOM',2X,'EL',7X,3(4X,'  X  ',1X,'  Y  ',1X,'  Z  '))")
     886          61 :          DO iatom = 1, natom, 3
     887          44 :             katom = iatom/3
     888          44 :             IF (MOD(iatom, 3) /= 0) katom = katom + 1
     889             :             CALL get_atomic_kind(atomic_kind=particles(Mlist(katom))%atomic_kind, &
     890          44 :                                  element_symbol=element_symbol)
     891             :             WRITE (UNIT=iw, FMT="(T2,I5,2X,A2,7X,3(4X,2(F5.2,1X),F5.2))") &
     892          44 :                Mlist(katom), element_symbol, &
     893         585 :                ((D(iatom + j, icol), j=0, 2), icol=from, to)
     894             :          END DO
     895          32 :          WRITE (UNIT=iw, FMT="(/)")
     896             :       END DO
     897             : 
     898          15 :    END SUBROUTINE vib_out
     899             : 
     900             : ! **************************************************************************************************
     901             : !> \brief Generates the transformation matrix from hessian in cartesian into
     902             : !>      internal coordinates (based on Gram-Schmidt orthogonalization)
     903             : !> \param mat ...
     904             : !> \param dof ...
     905             : !> \param Dout ...
     906             : !> \param full ...
     907             : !> \param natoms ...
     908             : !> \author Teodoro Laino 08.2006
     909             : ! **************************************************************************************************
     910          31 :    SUBROUTINE build_D_matrix(mat, dof, Dout, full, natoms)
     911             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: mat
     912             :       INTEGER, INTENT(IN)                                :: dof
     913             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: Dout
     914             :       LOGICAL, OPTIONAL                                  :: full
     915             :       INTEGER, INTENT(IN)                                :: natoms
     916             : 
     917             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'build_D_matrix'
     918             : 
     919             :       INTEGER                                            :: handle, i, ifound, iseq, j, nvib
     920             :       LOGICAL                                            :: my_full
     921             :       REAL(KIND=dp)                                      :: norm
     922          31 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work
     923          31 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: D
     924             : 
     925          31 :       CALL timeset(routineN, handle)
     926          31 :       my_full = .TRUE.
     927          31 :       IF (PRESENT(full)) my_full = full
     928             :       ! Generate the missing vectors of the orthogonal basis set
     929          31 :       nvib = 3*natoms - dof
     930         124 :       ALLOCATE (work(3*natoms))
     931         124 :       ALLOCATE (D(3*natoms, 3*natoms))
     932             :       ! Check First orthogonality in the first element of the basis set
     933         181 :       DO i = 1, dof
     934        1410 :          D(:, i) = mat(:, i)
     935         532 :          DO j = i + 1, dof
     936        3330 :             norm = DOT_PRODUCT(mat(:, i), mat(:, j))
     937         501 :             IF (ABS(norm) > thrs_motion) THEN
     938           0 :                CPWARN("Orthogonality error in transformation matrix")
     939             :             END IF
     940             :          END DO
     941             :       END DO
     942             :       ! Generate the nvib orthogonal vectors
     943             :       iseq = 0
     944             :       ifound = 0
     945         158 :       DO WHILE (ifound /= nvib)
     946         127 :          iseq = iseq + 1
     947         127 :          CPASSERT(iseq <= 3*natoms)
     948        1168 :          work = 0.0_dp
     949         127 :          work(iseq) = 1.0_dp
     950             :          ! Gram Schmidt orthogonalization
     951         896 :          DO i = 1, dof + ifound
     952        7414 :             norm = DOT_PRODUCT(work, D(:, i))
     953        7541 :             work(:) = work - norm*D(:, i)
     954             :          END DO
     955             :          ! Check norm of the new generated vector
     956        1168 :          norm = SQRT(DOT_PRODUCT(work, work))
     957         158 :          IF (norm >= 10E4_dp*thrs_motion) THEN
     958             :             ! Accept new vector
     959          93 :             ifound = ifound + 1
     960         840 :             D(:, dof + ifound) = work/norm
     961             :          END IF
     962             :       END DO
     963          31 :       CPASSERT(dof + ifound == 3*natoms)
     964          31 :       IF (my_full) THEN
     965          91 :          Dout = D
     966             :       ELSE
     967         840 :          Dout = D(:, dof + 1:)
     968             :       END IF
     969          31 :       DEALLOCATE (work)
     970          31 :       DEALLOCATE (D)
     971          31 :       CALL timestop(handle)
     972          31 :    END SUBROUTINE build_D_matrix
     973             : 
     974             : ! **************************************************************************************************
     975             : !> \brief Calculate a few thermochemical  properties from vibrational analysis
     976             : !>         It is supposed to work for molecules in the gas phase and without constraints
     977             : !> \param freqs ...
     978             : !> \param iw ...
     979             : !> \param mass ...
     980             : !> \param nvib ...
     981             : !> \param inertia ...
     982             : !> \param spin ...
     983             : !> \param totene ...
     984             : !> \param temp ...
     985             : !> \param pressure ...
     986             : !> \author MI 10:2015
     987             : ! **************************************************************************************************
     988             : 
     989           2 :    SUBROUTINE get_thch_values(freqs, iw, mass, nvib, inertia, spin, totene, temp, pressure)
     990             : 
     991             :       REAL(KIND=dp), DIMENSION(:)                        :: freqs
     992             :       INTEGER, INTENT(IN)                                :: iw
     993             :       REAL(KIND=dp), DIMENSION(:)                        :: mass
     994             :       INTEGER, INTENT(IN)                                :: nvib
     995             :       REAL(KIND=dp), INTENT(IN)                          :: inertia(3)
     996             :       INTEGER, INTENT(IN)                                :: spin
     997             :       REAL(KIND=dp), INTENT(IN)                          :: totene, temp, pressure
     998             : 
     999             :       INTEGER                                            :: i, natoms, sym_num
    1000             :       REAL(KIND=dp) :: el_entropy, entropy, exp_min_one, fact, fact2, freq_arg, freq_arg2, &
    1001             :          freqsum, Gibbs, heat_capacity, inertia_kg(3), mass_tot, one_min_exp, partition_function, &
    1002             :          rot_cv, rot_energy, rot_entropy, rot_part_func, rotvibtra, tran_cv, tran_energy, &
    1003             :          tran_enthalpy, tran_entropy, tran_part_func, vib_cv, vib_energy, vib_entropy, &
    1004             :          vib_part_func, zpe
    1005           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mass_kg
    1006             : 
    1007             : !    temp = 273.150_dp ! in Kelvin
    1008             : !    pressure = 101325.0_dp ! in Pascal
    1009             : 
    1010           2 :       freqsum = 0.0_dp
    1011           4 :       DO i = 1, nvib
    1012           4 :          freqsum = freqsum + freqs(i)
    1013             :       END DO
    1014             : 
    1015             : !   ZPE
    1016           2 :       zpe = 0.5_dp*(h_bar*2._dp*pi)*freqsum*(hertz/wavenumbers)*n_avogadro
    1017             : 
    1018           2 :       el_entropy = (n_avogadro*boltzmann)*LOG(REAL(spin, KIND=dp))
    1019             : !
    1020           2 :       natoms = SIZE(mass)
    1021           6 :       ALLOCATE (mass_kg(natoms))
    1022           6 :       mass_kg(:) = mass(:)**2*e_mass
    1023           6 :       mass_tot = SUM(mass_kg)
    1024           8 :       inertia_kg = inertia*e_mass*(a_bohr**2)
    1025             : 
    1026             : !   ROTATIONAL: Partition function and Entropy
    1027           2 :       sym_num = 1
    1028           2 :       fact = temp*2.0_dp*boltzmann/(h_bar*h_bar)
    1029           2 :       IF (inertia_kg(1)*inertia_kg(2)*inertia_kg(3) > 1.0_dp) THEN
    1030           0 :          rot_part_func = fact*fact*fact*inertia_kg(1)*inertia_kg(2)*inertia_kg(3)*pi
    1031           0 :          rot_part_func = SQRT(rot_part_func)
    1032           0 :          rot_entropy = n_avogadro*boltzmann*(LOG(rot_part_func) + 1.5_dp)
    1033           0 :          rot_energy = 1.5_dp*n_avogadro*boltzmann*temp
    1034           0 :          rot_cv = 1.5_dp*n_avogadro*boltzmann
    1035             :       ELSE
    1036             :          !linear molecule
    1037           2 :          IF (inertia_kg(1) > 1.0_dp) THEN
    1038           0 :             rot_part_func = fact*inertia_kg(1)
    1039           2 :          ELSE IF (inertia_kg(2) > 1.0_dp) THEN
    1040           0 :             rot_part_func = fact*inertia_kg(2)
    1041             :          ELSE
    1042           2 :             rot_part_func = fact*inertia_kg(3)
    1043             :          END IF
    1044           2 :          rot_entropy = n_avogadro*boltzmann*(LOG(rot_part_func) + 1.0_dp)
    1045           2 :          rot_energy = n_avogadro*boltzmann*temp
    1046           2 :          rot_cv = n_avogadro*boltzmann
    1047             :       END IF
    1048             : 
    1049             : !   TRANSLATIONAL: Partition function and Entropy
    1050           2 :       tran_part_func = (boltzmann*temp)**2.5_dp/(pressure*(h_bar*2.0_dp*pi)**3.0_dp)*(2.0_dp*pi*mass_tot)**1.5_dp
    1051           2 :       tran_entropy = n_avogadro*boltzmann*(LOG(tran_part_func) + 2.5_dp)
    1052           2 :       tran_energy = 1.5_dp*n_avogadro*boltzmann*temp
    1053           2 :       tran_enthalpy = 2.5_dp*n_avogadro*boltzmann*temp
    1054           2 :       tran_cv = 2.5_dp*n_avogadro*boltzmann
    1055             : 
    1056             : !   VIBRATIONAL:  Partition function and Entropy
    1057           2 :       vib_part_func = 1.0_dp
    1058           2 :       vib_energy = 0.0_dp
    1059           2 :       vib_entropy = 0.0_dp
    1060           2 :       vib_cv = 0.0_dp
    1061           2 :       fact = 2.0_dp*pi*h_bar/boltzmann/temp*hertz/wavenumbers
    1062           2 :       fact2 = 2.0_dp*pi*h_bar*hertz/wavenumbers
    1063           4 :       DO i = 1, nvib
    1064           2 :          freq_arg = fact*freqs(i)
    1065           2 :          freq_arg2 = fact2*freqs(i)
    1066           2 :          exp_min_one = EXP(freq_arg) - 1.0_dp
    1067           2 :          one_min_exp = 1.0_dp - EXP(-freq_arg)
    1068             : !dbg
    1069             : !  write(*,*) 'freq ', i, freqs(i), exp_min_one , one_min_exp
    1070             : !      vib_part_func = vib_part_func*(1.0_dp/(1.0_dp - exp(-fact*freqs(i))))
    1071           2 :          vib_part_func = vib_part_func*(1.0_dp/one_min_exp)
    1072             : !      vib_energy = vib_energy + fact2*freqs(i)*0.5_dp+fact2*freqs(i)/(exp(fact*freqs(i))-1.0_dp)
    1073           2 :          vib_energy = vib_energy + freq_arg2*0.5_dp + freq_arg2/exp_min_one
    1074             : !      vib_entropy = vib_entropy +fact*freqs(i)/(exp(fact*freqs(i))-1.0_dp)-log(1.0_dp - exp(-fact*freqs(i)))
    1075           2 :          vib_entropy = vib_entropy + freq_arg/exp_min_one - LOG(one_min_exp)
    1076             : !      vib_cv = vib_cv + fact*fact*freqs(i)*freqs(i)*exp(fact*freqs(i))/(exp(fact*freqs(i))-1.0_dp)/(exp(fact*freqs(i))-1.0_dp)
    1077           4 :          vib_cv = vib_cv + freq_arg*freq_arg*EXP(freq_arg)/exp_min_one/exp_min_one
    1078             :       END DO
    1079           2 :       vib_energy = vib_energy*n_avogadro ! it contains already ZPE
    1080           2 :       vib_entropy = vib_entropy*(n_avogadro*boltzmann)
    1081           2 :       vib_cv = vib_cv*(n_avogadro*boltzmann)
    1082             : 
    1083             : !   SUMMARY
    1084             : !dbg
    1085             : !    write(*,*) 'part ', rot_part_func,tran_part_func,vib_part_func
    1086           2 :       partition_function = rot_part_func*tran_part_func*vib_part_func
    1087             : !dbg
    1088             : !    write(*,*) 'entropy ', el_entropy,rot_entropy,tran_entropy,vib_entropy
    1089             : 
    1090           2 :       entropy = el_entropy + rot_entropy + tran_entropy + vib_entropy
    1091             : !dbg
    1092             : !    write(*,*) 'energy ', rot_energy , tran_enthalpy , vib_energy, totene*kjmol*1000.0_dp
    1093             : 
    1094           2 :       rotvibtra = rot_energy + tran_enthalpy + vib_energy
    1095             : !dbg
    1096             : !    write(*,*) 'cv ', rot_cv, tran_cv, vib_cv
    1097           2 :       heat_capacity = vib_cv + tran_cv + rot_cv
    1098             : 
    1099             : !   Free energy in J/mol: internal energy + PV - TS
    1100           2 :       Gibbs = vib_energy + rot_energy + tran_enthalpy - temp*entropy
    1101             : 
    1102           2 :       DEALLOCATE (mass_kg)
    1103             : 
    1104           2 :       IF (iw > 0) THEN
    1105           1 :          WRITE (UNIT=iw, FMT="(/,T2,'VIB|',T30,'NORMAL MODES - THERMOCHEMICAL DATA')")
    1106           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|')")
    1107             : 
    1108           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Symmetry number:',T70,I16)") sym_num
    1109           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Temperature [K]:',T70,F16.2)") temp
    1110           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Pressure [Pa]:',T70,F16.2)") pressure
    1111             : 
    1112           1 :          WRITE (UNIT=iw, FMT="(/)")
    1113             : 
    1114           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Electronic energy (U) [kJ/mol]:',T60,F26.8)") totene*kjmol
    1115           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Zero-point correction [kJ/mol]:',T60,F26.8)") zpe/1000.0_dp
    1116           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Entropy [kJ/(mol K)]:',T60,F26.8)") entropy/1000.0_dp
    1117           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Enthalpy correction (H-U) [kJ/mol]:',T60,F26.8)") rotvibtra/1000.0_dp
    1118           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Gibbs energy correction [kJ/mol]:',T60,F26.8)") Gibbs/1000.0_dp
    1119           1 :          WRITE (UNIT=iw, FMT="(T2,'VIB|', T20, 'Heat capacity [kJ/(mol*K)]:',T70,F16.8)") heat_capacity/1000.0_dp
    1120           1 :          WRITE (UNIT=iw, FMT="(/)")
    1121             :       END IF
    1122             : 
    1123           2 :    END SUBROUTINE get_thch_values
    1124             : 
    1125             : ! **************************************************************************************************
    1126             : !> \brief write out the non-orthogalized, i.e. without rotation and translational symmetry removed,
    1127             : !>        eigenvalues and eigenvectors of the Cartesian Hessian in unformatted binary file
    1128             : !> \param unit : the output unit to write to
    1129             : !> \param dof  : total degrees of freedom, i.e. the rank of the Hessian matrix
    1130             : !> \param eigenvalues  : eigenvalues of the Hessian matrix
    1131             : !> \param eigenvectors : matrix with each column being the eigenvectors of the Hessian matrix
    1132             : !> \author Lianheng Tong - 2016/04/20
    1133             : ! **************************************************************************************************
    1134          16 :    SUBROUTINE write_eigs_unformatted(unit, dof, eigenvalues, eigenvectors)
    1135             :       INTEGER, INTENT(IN)                                :: unit, dof
    1136             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvalues
    1137             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: eigenvectors
    1138             : 
    1139             :       CHARACTER(len=*), PARAMETER :: routineN = 'write_eigs_unformatted'
    1140             : 
    1141             :       INTEGER                                            :: handle, jj
    1142             : 
    1143          16 :       CALL timeset(routineN, handle)
    1144          16 :       IF (unit .GT. 0) THEN
    1145             :          ! degrees of freedom, i.e. the rank
    1146          16 :          WRITE (unit) dof
    1147             :          ! eigenvalues in one record
    1148          16 :          WRITE (unit) eigenvalues(1:dof)
    1149             :          ! eigenvectors: each record contains an eigenvector
    1150         142 :          DO jj = 1, dof
    1151         142 :             WRITE (unit) eigenvectors(1:dof, jj)
    1152             :          END DO
    1153             :       END IF
    1154          16 :       CALL timestop(handle)
    1155             : 
    1156          16 :    END SUBROUTINE write_eigs_unformatted
    1157             : 
    1158             : !**************************************************************************************************
    1159             : !> \brief Write the Hessian matrix into a (unformatted) binary file
    1160             : !> \param vib_section vibrational analysis section
    1161             : !> \param para_env mpi environment
    1162             : !> \param ncoord 3 times the number of atoms
    1163             : !> \param globenv global environment
    1164             : !> \param Hessian the Hessian matrix
    1165             : !> \param logger the logger
    1166             : ! **************************************************************************************************
    1167          30 :    SUBROUTINE write_va_hessian(vib_section, para_env, ncoord, globenv, Hessian, logger)
    1168             : 
    1169             :       TYPE(section_vals_type), POINTER                   :: vib_section
    1170             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1171             :       INTEGER                                            :: ncoord
    1172             :       TYPE(global_environment_type), POINTER             :: globenv
    1173             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Hessian
    1174             :       TYPE(cp_logger_type), POINTER                      :: logger
    1175             : 
    1176             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_va_hessian'
    1177             : 
    1178             :       INTEGER                                            :: handle, hesunit, i, j, ndf
    1179             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1180             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_hes
    1181             :       TYPE(cp_fm_type)                                   :: hess_mat
    1182             : 
    1183          30 :       CALL timeset(routineN, handle)
    1184             : 
    1185             :       hesunit = cp_print_key_unit_nr(logger, vib_section, "PRINT%HESSIAN", &
    1186             :                                      extension=".hess", file_form="UNFORMATTED", file_action="WRITE", &
    1187          30 :                                      file_position="REWIND")
    1188             : 
    1189          30 :       NULLIFY (blacs_env)
    1190             :       CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
    1191          30 :                                globenv%blacs_repeatable)
    1192          30 :       ndf = ncoord
    1193             :       CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
    1194          30 :                                nrow_global=ndf, ncol_global=ndf)
    1195          30 :       CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
    1196          30 :       CALL cp_fm_set_all(hess_mat, alpha=0.0_dp, beta=0.0_dp)
    1197             : 
    1198         264 :       DO i = 1, ncoord
    1199        2190 :          DO j = 1, ncoord
    1200        2160 :             CALL cp_fm_set_element(hess_mat, i, j, Hessian(i, j))
    1201             :          END DO
    1202             :       END DO
    1203          30 :       CALL cp_fm_write_unformatted(hess_mat, hesunit)
    1204             : 
    1205          30 :       CALL cp_print_key_finished_output(hesunit, logger, vib_section, "PRINT%HESSIAN")
    1206             : 
    1207          30 :       CALL cp_fm_struct_release(fm_struct_hes)
    1208          30 :       CALL cp_fm_release(hess_mat)
    1209          30 :       CALL cp_blacs_env_release(blacs_env)
    1210             : 
    1211          30 :       CALL timestop(handle)
    1212             : 
    1213          30 :    END SUBROUTINE write_va_hessian
    1214             : 
    1215         254 : END MODULE vibrational_analysis

Generated by: LCOV version 1.15