LCOV - code coverage report
Current view: top level - src/motion - vibrational_analysis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.2 % 595 584
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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 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            8 :                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            3 :                ALLOCATE (Dfull(ncoord, ncoord))
     533            4 :                ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2)))
     534            3 :                ALLOCATE (HeigvalDfull(SIZE(Dfull, 2)))
     535            3 :                ALLOCATE (MatM(ncoord, ncoord))
     536            2 :                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         3280 :                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         3280 :                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              :             nvib = ncoord - nRotTrM
     570           60 :             ALLOCATE (H_eigval1(ncoord))
     571           90 :             ALLOCATE (H_eigval2(SIZE(D, 2)))
     572           90 :             ALLOCATE (Hint1(ncoord, ncoord))
     573          120 :             ALLOCATE (Hint2(SIZE(D, 2), SIZE(D, 2)))
     574           60 :             ALLOCATE (rmass(SIZE(D, 2)))
     575           60 :             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           60 :             ALLOCATE (intensities_d(SIZE(D, 2)))
     583           60 :             ALLOCATE (intensities_p(SIZE(D, 2)))
     584           60 :             ALLOCATE (depol_p(SIZE(D, 2)))
     585           60 :             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        25188 :                Hint2(:, :) = MATMUL(TRANSPOSE(D), MATMUL(Hessian, D))
     605           30 :                IF (calc_intens) THEN
     606           64 :                   DO i = 1, 3
     607         4630 :                      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        13890 :                         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            1 :                   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            8 :          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              :       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           60 :    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 2.0-1