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