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