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