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