(git:ed6f26b)
Loading...
Searching...
No Matches
vibrational_analysis.F
Go to the documentation of this file.
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! **************************************************************************************************
20MODULE 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
91CONTAINS
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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: h_eigval1, h_eigval2, heigvaldfull, &
121 konst, mass, pos0, rmass
122 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 REAL(kind=dp), DIMENSION(:), POINTER :: depol_p, depol_u, depp, depu, din, &
126 intensities_d, intensities_p, pin
127 REAL(kind=dp), DIMENSION(:, :), POINTER :: d, dfull, dip_deriv, rottrm
128 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: polar_deriv, tmp_dip
129 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 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 CALL timeset(routinen, handle)
140 NULLIFY (d, rottrm, logger, subsys, f_env, particles, rep_env, intensities_d, intensities_p, &
141 vib_section, print_section, depol_p, depol_u)
142 logger => cp_get_default_logger()
143 vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS")
144 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 extension=".vibLog")
149 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 file_form="UNFORMATTED")
160
161 CALL section_vals_val_get(vib_section, "DX", r_val=dx)
162 CALL section_vals_val_get(vib_section, "NPROC_REP", i_val=prep)
163 CALL section_vals_val_get(vib_section, "PROC_DIST_TYPE", i_val=proc_dist_type)
164 row_force = (proc_dist_type == do_rep_blocked)
165 CALL section_vals_val_get(vib_section, "FULLY_PERIODIC", l_val=keep_rotations)
166 CALL section_vals_val_get(vib_section, "INTENSITIES", l_val=calc_intens)
167 CALL section_vals_val_get(vib_section, "THERMOCHEMISTRY", l_val=calc_thchdata)
168 CALL section_vals_val_get(vib_section, "TC_TEMPERATURE", r_val=tc_temp)
169 CALL section_vals_val_get(vib_section, "TC_PRESSURE", r_val=tc_press)
170
171 tc_temp = tc_temp*kelvin
172 tc_press = tc_press*pascal
173
174 intens_ir = .false.
175 intens_raman = .false.
176
177 mode_tracking_section => section_vals_get_subs_vals(vib_section, "MODE_SELECTIVE")
178 CALL section_vals_get(mode_tracking_section, explicit=do_mode_tracking)
179 nrep = max(1, para_env%num_pe/prep)
180 prep = para_env%num_pe/nrep
181 iw = cp_print_key_unit_nr(logger, print_section, "BANNER", extension=".vibLog")
182 CALL vib_header(iw, nrep, prep)
183 CALL cp_print_key_finished_output(iw, logger, print_section, "BANNER")
184 ! Just one force_env allowed
185 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 input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=row_force)
189 IF (ASSOCIATED(rep_env)) THEN
190 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env)
191 CALL force_env_get(f_env%force_env, subsys=subsys)
192 particles => subsys%particles%els
193 ! Decide which kind of Vibrational Analysis to perform
194 IF (do_mode_tracking) THEN
195 CALL ms_vb_anal(input, rep_env, para_env, globenv, particles, &
196 nrep, calc_intens, dx, output_unit, logger)
197 CALL f_env_rm_defaults(f_env, ierr)
198 ELSE
199 CALL get_moving_atoms(force_env=f_env%force_env, ilist=mlist)
200 something_frozen = SIZE(particles) .NE. SIZE(mlist)
201 natoms = SIZE(mlist)
202 ncoord = natoms*3
203 ALLOCATE (clist(ncoord))
204 ALLOCATE (mass(natoms))
205 ALLOCATE (pos0(ncoord))
206 ALLOCATE (hessian(ncoord, ncoord))
207 IF (calc_intens) THEN
208 description_d = '[DIPOLE]'
209 ALLOCATE (tmp_dip(ncoord, 3, 2))
210 tmp_dip = 0._dp
211 description_p = '[POLAR]'
212 ALLOCATE (tmp_polar(ncoord, 3, 3, 2))
213 tmp_polar = 0._dp
214 END IF
215 clist = 0
216 DO i = 1, natoms
217 imap = mlist(i)
218 clist((i - 1)*3 + 1) = (imap - 1)*3 + 1
219 clist((i - 1)*3 + 2) = (imap - 1)*3 + 2
220 clist((i - 1)*3 + 3) = (imap - 1)*3 + 3
221 mass(i) = particles(imap)%atomic_kind%mass
222 cpassert(mass(i) > 0.0_dp)
223 mass(i) = sqrt(mass(i))
224 pos0((i - 1)*3 + 1) = particles(imap)%r(1)
225 pos0((i - 1)*3 + 2) = particles(imap)%r(2)
226 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 IF (something_frozen) THEN
233 nrottrm = 0
234 ALLOCATE (rottrm(natoms*3, nrottrm))
235 ELSE
236 CALL rot_ana(particles, rottrm, nrottrm, print_section, &
237 keep_rotations, mass_weighted=.true., natoms=natoms, inertia=inertia)
238 END IF
239 ! Generate the suitable rototranslating basis set
240 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 ALLOCATE (d(3*natoms, nvib))
246 END IF
247 CALL build_d_matrix(rottrm, nrottrm, d, full=.false., &
248 natoms=natoms)
249 !
250 ! Loop on atoms and coordinates
251 !
252 hessian = huge(0.0_dp)
253 IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "VIB| Vibrational Analysis Info"
254 DO icoordp = 1, ncoord, nrep
255 icoord = icoordp - 1
256 DO j = 1, nrep
257 DO i = 1, ncoord
258 imap = clist(i)
259 rep_env%r(imap, j) = pos0(i)
260 END DO
261 IF (icoord + j <= ncoord) THEN
262 imap = clist(icoord + j)
263 rep_env%r(imap, j) = rep_env%r(imap, j) + dx
264 END IF
265 END DO
266 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
267
268 DO j = 1, nrep
269 IF (calc_intens) THEN
270 IF (icoord + j <= ncoord) THEN
271 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 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 nval=nres)
280 intens_ir = .true.
281 d_print(:) = tmp_dip(icoord + j, :, 1)
282 END IF
283 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 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 nval=nres)
292 intens_raman = .true.
293 p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
294 END IF
295 END IF
296 END IF
297 IF (icoord + j <= ncoord) THEN
298 DO i = 1, ncoord
299 imap = clist(i)
300 hessian(i, icoord + j) = rep_env%f(imap, j)
301 END DO
302 imap = clist(icoord + j)
303 ! Dump Info
304 IF (output_unit > 0) THEN
305 iparticle1 = imap/3
306 IF (mod(imap, 3) /= 0) iparticle1 = iparticle1 + 1
307 WRITE (output_unit, '(T2,A,I5,A,I5,3A)') &
308 "VIB| REPLICA Nr.", j, "- Energy and Forces for particle:", &
309 iparticle1, " coordinate: ", lab(imap - (iparticle1 - 1)*3), &
310 " + D"//trim(lab(imap - (iparticle1 - 1)*3))
311 WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
312 "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
313 WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
314 DO i = 1, natoms
315 imap = mlist(i)
316 WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
317 particles(imap)%atomic_kind%name, &
318 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
319 END DO
320 IF (intens_ir) THEN
321 WRITE (output_unit, '(T3,A)') 'Dipole moment [Debye]'
322 WRITE (output_unit, '(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
323 'X=', d_print(1)*debye, 'Y=', d_print(2)*debye, 'Z=', d_print(3)*debye, &
324 'Total=', sqrt(sum(d_print(1:3)**2))*debye
325 END IF
326 IF (intens_raman) THEN
327 WRITE (output_unit, '(T2,A)') &
328 'POLAR| Polarizability tensor [a.u.]'
329 WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
330 '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 '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 '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 DO icoordm = 1, ncoord, nrep
341 icoord = icoordm - 1
342 DO j = 1, nrep
343 DO i = 1, ncoord
344 imap = clist(i)
345 rep_env%r(imap, j) = pos0(i)
346 END DO
347 IF (icoord + j <= ncoord) THEN
348 imap = clist(icoord + j)
349 rep_env%r(imap, j) = rep_env%r(imap, j) - dx
350 END IF
351 END DO
352 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
353
354 DO j = 1, nrep
355 IF (calc_intens) THEN
356 IF (icoord + j <= ncoord) THEN
357 k = (icoord + j + 2)/3
358 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 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 nval=nres)
367 tmp_dip(icoord + j, :, 1) = (tmp_dip(icoord + j, :, 1) - &
368 tmp_dip(icoord + j, :, 2))/(2.0_dp*dx*mass(k))
369 d_print(:) = tmp_dip(icoord + j, :, 1)
370 END IF
371 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 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 nval=nres)
380 tmp_polar(icoord + j, :, :, 1) = (tmp_polar(icoord + j, :, :, 1) - &
381 tmp_polar(icoord + j, :, :, 2))/(2.0_dp*dx*mass(k))
382 p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
383 END IF
384 END IF
385 END IF
386 IF (icoord + j <= ncoord) THEN
387 imap = clist(icoord + j)
388 iparticle1 = imap/3
389 IF (mod(imap, 3) /= 0) iparticle1 = iparticle1 + 1
390 ip1 = (icoord + j)/3
391 IF (mod(icoord + j, 3) /= 0) ip1 = ip1 + 1
392 ! Dump Info
393 IF (output_unit > 0) THEN
394 WRITE (output_unit, '(T2,A,I5,A,I5,3A)') &
395 "VIB| REPLICA Nr.", j, "- Energy and Forces for particle:", &
396 iparticle1, " coordinate: ", lab(imap - (iparticle1 - 1)*3), &
397 " - D"//trim(lab(imap - (iparticle1 - 1)*3))
398 WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
399 "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
400 WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
401 DO i = 1, natoms
402 imap = mlist(i)
403 WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
404 particles(imap)%atomic_kind%name, &
405 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
406 END DO
407 IF (intens_ir) THEN
408 WRITE (output_unit, '(T3,A)') 'Dipole moment [Debye]'
409 WRITE (output_unit, '(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
410 'X=', d_print(1)*debye, 'Y=', d_print(2)*debye, 'Z=', d_print(3)*debye, &
411 'Total=', sqrt(sum(d_print(1:3)**2))*debye
412 END IF
413 IF (intens_raman) THEN
414 WRITE (output_unit, '(T2,A)') &
415 'POLAR| Polarizability tensor [a.u.]'
416 WRITE (output_unit, '(T2,A,T24,3(1X,F18.12))') &
417 '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 '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 'POLAR| yx,zx,zy', p_print(2, 1), p_print(3, 1), p_print(3, 2)
422 END IF
423 END IF
424 DO iseq = 1, ncoord
425 imap = clist(iseq)
426 iparticle2 = imap/3
427 IF (mod(imap, 3) /= 0) iparticle2 = iparticle2 + 1
428 ip2 = iseq/3
429 IF (mod(iseq, 3) /= 0) ip2 = ip2 + 1
430 tmp = hessian(iseq, icoord + j) - rep_env%f(imap, j)
431 tmp = -tmp/(2.0_dp*dx*mass(ip1)*mass(ip2))*1e6_dp
432 ! Mass weighted Hessian
433 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 DO i = 1, natoms
441 imap = mlist(i)
442 particles(imap)%r(1:3) = pos0((i - 1)*3 + 1:(i - 1)*3 + 3)
443 END DO
444 DO j = 1, nrep
445 DO i = 1, ncoord
446 imap = clist(i)
447 rep_env%r(imap, j) = pos0(i)
448 END DO
449 END DO
450 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
451 j = 1
452 minimum_energy = rep_env%f(rep_env%ndim + 1, j)
453 IF (output_unit > 0) THEN
454 WRITE (output_unit, '(T2,A)') &
455 "VIB| ", " Minimum Structure - Energy and Forces:"
456 WRITE (output_unit, '(T2,A,T43,A,T57,F24.12)') &
457 "VIB|", "Total energy:", rep_env%f(rep_env%ndim + 1, j)
458 WRITE (output_unit, '(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
459 DO i = 1, natoms
460 imap = mlist(i)
461 WRITE (output_unit, '(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
462 particles(imap)%atomic_kind%name, &
463 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
464 END DO
465 END IF
466
467 ! Dump Info
468 IF (output_unit > 0) THEN
469 WRITE (output_unit, '(T2,A)') "VIB| Hessian in cartesian coordinates"
470 CALL write_particle_matrix(hessian, particles, output_unit, el_per_part=3, &
471 ilist=mlist)
472 END IF
473
474 CALL write_va_hessian(vib_section, para_env, ncoord, globenv, hessian, logger)
475
476 ! Enforce symmetry in the Hessian
477 DO i = 1, ncoord
478 DO j = i, ncoord
479 ! Take the upper diagonal part
480 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 file_position="REWIND", extension=".rrm")
487 IF (print_grrm > 0) THEN
488 DO i = 1, natoms
489 imap = mlist(i)
490 particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
491 END DO
492 ALLOCATE (hint1(ncoord, ncoord), rmass(ncoord))
493 DO i = 1, natoms
494 imap = mlist(i)
495 rmass(3*(imap - 1) + 1:3*(imap - 1) + 3) = mass(imap)
496 END DO
497 DO i = 1, ncoord
498 DO j = 1, ncoord
499 hint1(j, i) = hessian(j, i)*rmass(i)*rmass(j)*1.0e-6_dp
500 END DO
501 END DO
502 nfrozen = SIZE(particles) - natoms
503 CALL write_grrm(print_grrm, f_env%force_env, particles, minimum_energy, &
504 hessian=hint1, fixed_atoms=nfrozen)
505 DEALLOCATE (hint1, rmass)
506 END IF
507 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 file_position="REWIND", extension=".scine")
512 IF (print_scine > 0) THEN
513 DO i = 1, natoms
514 imap = mlist(i)
515 particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
516 END DO
517 nfrozen = SIZE(particles) - natoms
518 cpassert(nfrozen == 0)
519 CALL write_scine(print_scine, f_env%force_env, particles, minimum_energy, hessian=hessian)
520 END IF
521 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 file_form="UNFORMATTED")
528 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 ALLOCATE (dfull(ncoord, ncoord))
533 ALLOCATE (hint2dfull(SIZE(dfull, 2), SIZE(dfull, 2)))
534 ALLOCATE (heigvaldfull(SIZE(dfull, 2)))
535 ALLOCATE (matm(ncoord, ncoord))
536 ALLOCATE (rmass(SIZE(dfull, 2)))
537 dfull = 0.0_dp
538 ! Dfull in dimension of degrees of freedom
539 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 hint2dfull(:, :) = matmul(transpose(dfull), matmul(hessian, dfull))
543 ! Heig = L^T Hint2Dfull L
544 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 matm = 0.0_dp
548 DO i = 1, natoms
549 DO j = 1, 3
550 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 dfull = matmul(matm, matmul(dfull, hint2dfull)) !Dfull=D L / sqrt(m)
555 DO i = 1, ncoord
556 ! Renormalize displacements
557 norm = 1.0_dp/sum(dfull(:, i)*dfull(:, i))
558 rmass(i) = norm/massunit
559 dfull(:, i) = sqrt(norm)*(dfull(:, i))
560 END DO
561 CALL write_eigs_unformatted(print_namd, ncoord, heigvaldfull, dfull)
562 DEALLOCATE (heigvaldfull)
563 DEALLOCATE (hint2dfull)
564 DEALLOCATE (dfull)
565 DEALLOCATE (matm)
566 DEALLOCATE (rmass)
567 END IF !print_namd
568 !
569 nvib = ncoord - nrottrm
570 ALLOCATE (h_eigval1(ncoord))
571 ALLOCATE (h_eigval2(SIZE(d, 2)))
572 ALLOCATE (hint1(ncoord, ncoord))
573 ALLOCATE (hint2(SIZE(d, 2), SIZE(d, 2)))
574 ALLOCATE (rmass(SIZE(d, 2)))
575 ALLOCATE (konst(SIZE(d, 2)))
576 IF (calc_intens) THEN
577 ALLOCATE (dip_deriv(3, SIZE(d, 2)))
578 dip_deriv = 0.0_dp
579 ALLOCATE (polar_deriv(3, 3, SIZE(d, 2)))
580 polar_deriv = 0.0_dp
581 END IF
582 ALLOCATE (intensities_d(SIZE(d, 2)))
583 ALLOCATE (intensities_p(SIZE(d, 2)))
584 ALLOCATE (depol_p(SIZE(d, 2)))
585 ALLOCATE (depol_u(SIZE(d, 2)))
586 intensities_d = 0._dp
587 intensities_p = 0._dp
588 depol_p = 0._dp
589 depol_u = 0._dp
590 hint1(:, :) = hessian
591 CALL diamat_all(hint1, h_eigval1)
592 IF (output_unit > 0) THEN
593 WRITE (output_unit, '(T2,"VIB| Cartesian Low frequencies ---",4G12.5)') &
594 (h_eigval1(i), i=1, min(9, ncoord))
595 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 ilist=mlist)
598 END IF
599 ! write frequencies and eigenvectors to cartesian eig file
600 IF (output_unit_eig > 0) THEN
601 CALL write_eigs_unformatted(output_unit_eig, ncoord, h_eigval1, hint1)
602 END IF
603 IF (nvib /= 0) THEN
604 hint2(:, :) = matmul(transpose(d), matmul(hessian, d))
605 IF (calc_intens) THEN
606 DO i = 1, 3
607 dip_deriv(i, :) = matmul(tmp_dip(:, i, 1), d)
608 END DO
609 DO i = 1, 3
610 DO j = 1, 3
611 polar_deriv(i, j, :) = matmul(tmp_polar(:, i, j, 1), d)
612 END DO
613 END DO
614 END IF
615 CALL diamat_all(hint2, h_eigval2)
616 IF (output_unit > 0) THEN
617 WRITE (output_unit, '(T2,"VIB| Frequencies after removal of the rotations and translations")')
618 ! Frequency at the moment are in a.u.
619 WRITE (output_unit, '(T2,"VIB| Internal Low frequencies ---",4G12.5)') h_eigval2
620 END IF
621 hessian = 0.0_dp
622 DO i = 1, natoms
623 DO j = 1, 3
624 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 d = matmul(hessian, matmul(d, hint2))
629 DO i = 1, nvib
630 norm = 1.0_dp/sum(d(:, i)*d(:, i))
631 ! Reduced Masess
632 rmass(i) = norm/massunit
633 ! Renormalize displacements and convert in Angstrom
634 d(:, i) = sqrt(norm)*d(:, i)
635 ! Force constants
636 konst(i) = sign(1.0_dp, h_eigval2(i))*2.0_dp*pi**2*(abs(h_eigval2(i))/massunit)**2*rmass(i)
637 IF (calc_intens) THEN
638 d_deriv = 0._dp
639 DO j = 1, nvib
640 d_deriv(:) = d_deriv(:) + dip_deriv(:, j)*hint2(j, i)
641 END DO
642 intensities_d(i) = sqrt(dot_product(d_deriv, d_deriv))
643 p_deriv = 0._dp
644 DO j = 1, nvib
645 ! P_deriv has units bohr^2/sqrt(a.u.)
646 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 p_deriv(:, :) = p_deriv(:, :)*conver
651 ! this is wron, just for testing
652 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 (p_deriv(3, 3) - p_deriv(1, 1))**2
656 a3 = (p_deriv(1, 2)**2 + p_deriv(2, 3)**2 + p_deriv(3, 1)**2)
657 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 dummy = 45.0_dp*a1*a1 + 4.0_dp/2.0_dp*(a2 + 6.0_dp*a3)
660 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 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 7.0_dp/2.0_dp*(a2 + 6.0_dp*a3))
667 ELSE
668 depol_p(i) = -1.0_dp
669 depol_u(i) = -1.0_dp
670 END IF
671 END IF
672 ! Convert frequencies to cm^-1
673 h_eigval2(i) = sign(1.0_dp, h_eigval2(i))*sqrt(abs(h_eigval2(i))*massunit)*vibfac/1000.0_dp
674 END DO
675 IF (calc_intens) THEN
676 IF (iounit > 0) THEN
677 IF (.NOT. intens_ir) THEN
678 WRITE (iounit, '(T2,"VIB| No IR intensities available. Check input")')
679 END IF
680 IF (.NOT. intens_raman) THEN
681 WRITE (iounit, '(T2,"VIB| No Raman intensities available. Check input")')
682 END IF
683 END IF
684 END IF
685 ! Dump Info
686 iw = cp_logger_get_default_io_unit(logger)
687 IF (iw > 0) THEN
688 NULLIFY (din, pin, depp, depu)
689 IF (intens_ir) din => intensities_d
690 IF (intens_raman) pin => intensities_p
691 IF (intens_raman) depp => depol_p
692 IF (intens_raman) depu => depol_u
693 CALL vib_out(iw, nvib, d, konst, rmass, h_eigval2, particles, mlist, din, pin, depp, depu)
694 END IF
695 IF (.NOT. something_frozen .AND. calc_thchdata) THEN
696 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 dump_only_positive=.false., logger=logger, list=mlist)
700 ELSE
701 IF (output_unit > 0) THEN
702 WRITE (output_unit, '(T2,"VIB| No further vibrational info. Detected a single atom")')
703 END IF
704 END IF
705 ! Deallocate working arrays
706 DEALLOCATE (rottrm)
707 DEALLOCATE (clist)
708 DEALLOCATE (mlist)
709 DEALLOCATE (h_eigval1)
710 DEALLOCATE (h_eigval2)
711 DEALLOCATE (hint1)
712 DEALLOCATE (hint2)
713 DEALLOCATE (rmass)
714 DEALLOCATE (konst)
715 DEALLOCATE (mass)
716 DEALLOCATE (pos0)
717 DEALLOCATE (d)
718 DEALLOCATE (hessian)
719 IF (calc_intens) THEN
720 DEALLOCATE (dip_deriv)
721 DEALLOCATE (polar_deriv)
722 DEALLOCATE (tmp_dip)
723 DEALLOCATE (tmp_polar)
724 END IF
725 DEALLOCATE (intensities_d)
726 DEALLOCATE (intensities_p)
727 DEALLOCATE (depol_p)
728 DEALLOCATE (depol_u)
729 CALL f_env_rm_defaults(f_env, ierr)
730 END IF
731 END IF
732 CALL cp_print_key_finished_output(output_unit, logger, print_section, "PROGRAM_RUN_INFO")
733 CALL cp_print_key_finished_output(output_unit_eig, logger, print_section, "CARTESIAN_EIGS")
734 CALL rep_env_release(rep_env)
735 CALL timestop(handle)
736 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: ifixd_list, work
753 TYPE(cp_subsys_type), POINTER :: subsys
754 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
755 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
756 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
760
761 CALL timeset(routinen, handle)
762 CALL force_env_get(force_env=force_env, subsys=subsys)
763
764 CALL cp_subsys_get(subsys=subsys, particles=particles, &
765 molecule_kinds=molecule_kinds)
766
767 nkind = molecule_kinds%n_els
768 molecule_kind_set => molecule_kinds%els
769 particle_set => particles%els
770
771 ! Count the number of fixed atoms
772 nfixed_atoms_total = 0
773 DO ikind = 1, nkind
774 molecule_kind => molecule_kind_set(ikind)
775 CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
776 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
777 END DO
778 ndim = SIZE(particle_set) - nfixed_atoms_total
779 cpassert(ndim >= 0)
780 ALLOCATE (ilist(ndim))
781
782 IF (nfixed_atoms_total /= 0) THEN
783 ALLOCATE (ifixd_list(nfixed_atoms_total))
784 ALLOCATE (work(nfixed_atoms_total))
785 nfixed_atoms_total = 0
786 DO ikind = 1, nkind
787 molecule_kind => molecule_kind_set(ikind)
788 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
789 IF (ASSOCIATED(fixd_list)) THEN
790 DO ii = 1, SIZE(fixd_list)
791 IF (.NOT. fixd_list(ii)%restraint%active) THEN
792 nfixed_atoms_total = nfixed_atoms_total + 1
793 ifixd_list(nfixed_atoms_total) = fixd_list(ii)%fixd
794 END IF
795 END DO
796 END IF
797 END DO
798 CALL sort(ifixd_list, nfixed_atoms_total, work)
799
800 ndim = 0
801 j = 1
802 loop_count: DO i = 1, SIZE(particle_set)
803 DO WHILE (i > ifixd_list(j))
804 j = j + 1
805 IF (j > nfixed_atoms_total) EXIT loop_count
806 END DO
807 IF (i /= ifixd_list(j)) THEN
808 ndim = ndim + 1
809 ilist(ndim) = i
810 END IF
811 END DO loop_count
812 DEALLOCATE (ifixd_list)
813 DEALLOCATE (work)
814 ELSE
815 i = 1
816 ndim = 0
817 END IF
818 DO j = i, SIZE(particle_set)
819 ndim = ndim + 1
820 ilist(ndim) = j
821 END DO
822 CALL timestop(handle)
823
824 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 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 fint = 42.255_dp*massunit*debye**2*bohr**2
858 pint = 1.0_dp
859 natom = SIZE(d, 1)
860 WRITE (unit=iw, fmt="(/,T2,'VIB|',T30,'NORMAL MODES - CARTESIAN DISPLACEMENTS')")
861 WRITE (unit=iw, fmt="(T2,'VIB|')")
862 DO jatom = 1, nvib, 3
863 from = jatom
864 to = min(from + 2, nvib)
865 WRITE (unit=iw, fmt="(T2,'VIB|',13X,3(8X,I5,8X))") &
866 (icol, icol=from, to)
867 WRITE (unit=iw, fmt="(T2,'VIB|Frequency (cm^-1)',3(1X,F12.6,8X))") &
868 (freq(icol), icol=from, to)
869 IF (ASSOCIATED(intensities_d)) THEN
870 WRITE (unit=iw, fmt="(T2,'VIB|IR int (KM/Mole) ',3(1X,F12.6,8X))") &
871 (fint*intensities_d(icol)**2, icol=from, to)
872 END IF
873 IF (ASSOCIATED(intensities_p)) THEN
874 WRITE (unit=iw, fmt="(T2,'VIB|Raman (A^4/amu) ',3(1X,F12.6,8X))") &
875 (pint*intensities_p(icol), icol=from, to)
876 WRITE (unit=iw, fmt="(T2,'VIB|Depol Ratio (P) ',3(1X,F12.6,8X))") &
877 (depol_p(icol), icol=from, to)
878 WRITE (unit=iw, fmt="(T2,'VIB|Depol Ratio (U) ',3(1X,F12.6,8X))") &
879 (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 (m(icol), icol=from, to)
883 WRITE (unit=iw, fmt="(T2,'VIB|Frc consts (a.u.)',3(1X,F12.6,8X))") &
884 (k(icol), icol=from, to)
885 WRITE (unit=iw, fmt="(T2,' ATOM',2X,'EL',7X,3(4X,' X ',1X,' Y ',1X,' Z '))")
886 DO iatom = 1, natom, 3
887 katom = iatom/3
888 IF (mod(iatom, 3) /= 0) katom = katom + 1
889 CALL get_atomic_kind(atomic_kind=particles(mlist(katom))%atomic_kind, &
890 element_symbol=element_symbol)
891 WRITE (unit=iw, fmt="(T2,I5,2X,A2,7X,3(4X,2(F5.2,1X),F5.2))") &
892 mlist(katom), element_symbol, &
893 ((d(iatom + j, icol), j=0, 2), icol=from, to)
894 END DO
895 WRITE (unit=iw, fmt="(/)")
896 END DO
897
898 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
923 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d
924
925 CALL timeset(routinen, handle)
926 my_full = .true.
927 IF (PRESENT(full)) my_full = full
928 ! Generate the missing vectors of the orthogonal basis set
929 nvib = 3*natoms - dof
930 ALLOCATE (work(3*natoms))
931 ALLOCATE (d(3*natoms, 3*natoms))
932 ! Check First orthogonality in the first element of the basis set
933 DO i = 1, dof
934 d(:, i) = mat(:, i)
935 DO j = i + 1, dof
936 norm = dot_product(mat(:, i), mat(:, j))
937 IF (abs(norm) > thrs_motion) THEN
938 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 DO WHILE (ifound /= nvib)
946 iseq = iseq + 1
947 cpassert(iseq <= 3*natoms)
948 work = 0.0_dp
949 work(iseq) = 1.0_dp
950 ! Gram Schmidt orthogonalization
951 DO i = 1, dof + ifound
952 norm = dot_product(work, d(:, i))
953 work(:) = work - norm*d(:, i)
954 END DO
955 ! Check norm of the new generated vector
956 norm = sqrt(dot_product(work, work))
957 IF (norm >= 10e4_dp*thrs_motion) THEN
958 ! Accept new vector
959 ifound = ifound + 1
960 d(:, dof + ifound) = work/norm
961 END IF
962 END DO
963 cpassert(dof + ifound == 3*natoms)
964 IF (my_full) THEN
965 dout = d
966 ELSE
967 dout = d(:, dof + 1:)
968 END IF
969 DEALLOCATE (work)
970 DEALLOCATE (d)
971 CALL timestop(handle)
972 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mass_kg
1006
1007! temp = 273.150_dp ! in Kelvin
1008! pressure = 101325.0_dp ! in Pascal
1009
1010 freqsum = 0.0_dp
1011 DO i = 1, nvib
1012 freqsum = freqsum + freqs(i)
1013 END DO
1014
1015! ZPE
1016 zpe = 0.5_dp*(h_bar*2._dp*pi)*freqsum*(hertz/wavenumbers)*n_avogadro
1017
1018 el_entropy = (n_avogadro*boltzmann)*log(real(spin, kind=dp))
1019!
1020 natoms = SIZE(mass)
1021 ALLOCATE (mass_kg(natoms))
1022 mass_kg(:) = mass(:)**2*e_mass
1023 mass_tot = sum(mass_kg)
1024 inertia_kg = inertia*e_mass*(a_bohr**2)
1025
1026! ROTATIONAL: Partition function and Entropy
1027 sym_num = 1
1028 fact = temp*2.0_dp*boltzmann/(h_bar*h_bar)
1029 IF (inertia_kg(1)*inertia_kg(2)*inertia_kg(3) > 1.0_dp) THEN
1030 rot_part_func = fact*fact*fact*inertia_kg(1)*inertia_kg(2)*inertia_kg(3)*pi
1031 rot_part_func = sqrt(rot_part_func)
1032 rot_entropy = n_avogadro*boltzmann*(log(rot_part_func) + 1.5_dp)
1033 rot_energy = 1.5_dp*n_avogadro*boltzmann*temp
1034 rot_cv = 1.5_dp*n_avogadro*boltzmann
1035 ELSE
1036 !linear molecule
1037 IF (inertia_kg(1) > 1.0_dp) THEN
1038 rot_part_func = fact*inertia_kg(1)
1039 ELSE IF (inertia_kg(2) > 1.0_dp) THEN
1040 rot_part_func = fact*inertia_kg(2)
1041 ELSE
1042 rot_part_func = fact*inertia_kg(3)
1043 END IF
1044 rot_entropy = n_avogadro*boltzmann*(log(rot_part_func) + 1.0_dp)
1045 rot_energy = n_avogadro*boltzmann*temp
1046 rot_cv = n_avogadro*boltzmann
1047 END IF
1048
1049! TRANSLATIONAL: Partition function and Entropy
1050 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 tran_entropy = n_avogadro*boltzmann*(log(tran_part_func) + 2.5_dp)
1052 tran_energy = 1.5_dp*n_avogadro*boltzmann*temp
1053 tran_enthalpy = 2.5_dp*n_avogadro*boltzmann*temp
1054 tran_cv = 2.5_dp*n_avogadro*boltzmann
1055
1056! VIBRATIONAL: Partition function and Entropy
1057 vib_part_func = 1.0_dp
1058 vib_energy = 0.0_dp
1059 vib_entropy = 0.0_dp
1060 vib_cv = 0.0_dp
1061 fact = 2.0_dp*pi*h_bar/boltzmann/temp*hertz/wavenumbers
1062 fact2 = 2.0_dp*pi*h_bar*hertz/wavenumbers
1063 DO i = 1, nvib
1064 freq_arg = fact*freqs(i)
1065 freq_arg2 = fact2*freqs(i)
1066 exp_min_one = exp(freq_arg) - 1.0_dp
1067 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 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 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 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 vib_cv = vib_cv + freq_arg*freq_arg*exp(freq_arg)/exp_min_one/exp_min_one
1078 END DO
1079 vib_energy = vib_energy*n_avogadro ! it contains already ZPE
1080 vib_entropy = vib_entropy*(n_avogadro*boltzmann)
1081 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 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 rotvibtra = rot_energy + tran_enthalpy + vib_energy
1095!dbg
1096! write(*,*) 'cv ', rot_cv, tran_cv, vib_cv
1097 heat_capacity = vib_cv + tran_cv + rot_cv
1098
1099! Free energy in J/mol: internal energy + PV - TS
1100 gibbs = vib_energy + rot_energy + tran_enthalpy - temp*entropy
1101
1102 DEALLOCATE (mass_kg)
1103
1104 IF (iw > 0) THEN
1105 WRITE (unit=iw, fmt="(/,T2,'VIB|',T30,'NORMAL MODES - THERMOCHEMICAL DATA')")
1106 WRITE (unit=iw, fmt="(T2,'VIB|')")
1107
1108 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Symmetry number:',T70,I16)") sym_num
1109 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Temperature [K]:',T70,F16.2)") temp
1110 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Pressure [Pa]:',T70,F16.2)") pressure
1111
1112 WRITE (unit=iw, fmt="(/)")
1113
1114 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Electronic energy (U) [kJ/mol]:',T60,F26.8)") totene*kjmol
1115 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Zero-point correction [kJ/mol]:',T60,F26.8)") zpe/1000.0_dp
1116 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Entropy [kJ/(mol K)]:',T60,F26.8)") entropy/1000.0_dp
1117 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Enthalpy correction (H-U) [kJ/mol]:',T60,F26.8)") rotvibtra/1000.0_dp
1118 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Gibbs energy correction [kJ/mol]:',T60,F26.8)") gibbs/1000.0_dp
1119 WRITE (unit=iw, fmt="(T2,'VIB|', T20, 'Heat capacity [kJ/(mol*K)]:',T70,F16.8)") heat_capacity/1000.0_dp
1120 WRITE (unit=iw, fmt="(/)")
1121 END IF
1122
1123 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 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 CALL timeset(routinen, handle)
1144 IF (unit .GT. 0) THEN
1145 ! degrees of freedom, i.e. the rank
1146 WRITE (unit) dof
1147 ! eigenvalues in one record
1148 WRITE (unit) eigenvalues(1:dof)
1149 ! eigenvectors: each record contains an eigenvector
1150 DO jj = 1, dof
1151 WRITE (unit) eigenvectors(1:dof, jj)
1152 END DO
1153 END IF
1154 CALL timestop(handle)
1155
1156 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 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 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 file_position="REWIND")
1188
1189 NULLIFY (blacs_env)
1190 CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
1191 globenv%blacs_repeatable)
1192 ndf = ncoord
1193 CALL cp_fm_struct_create(fm_struct_hes, para_env=para_env, context=blacs_env, &
1194 nrow_global=ndf, ncol_global=ndf)
1195 CALL cp_fm_create(hess_mat, fm_struct_hes, name="hess_mat")
1196 CALL cp_fm_set_all(hess_mat, alpha=0.0_dp, beta=0.0_dp)
1197
1198 DO i = 1, ncoord
1199 DO j = 1, ncoord
1200 CALL cp_fm_set_element(hess_mat, i, j, hessian(i, j))
1201 END DO
1202 END DO
1203 CALL cp_fm_write_unformatted(hess_mat, hesunit)
1204
1205 CALL cp_print_key_finished_output(hesunit, logger, vib_section, "PRINT%HESSIAN")
1206
1207 CALL cp_fm_struct_release(fm_struct_hes)
1208 CALL cp_fm_release(hess_mat)
1209 CALL cp_blacs_env_release(blacs_env)
1210
1211 CALL timestop(handle)
1212
1213 END SUBROUTINE write_va_hessian
1214
1215END MODULE vibrational_analysis
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_result_methods::get_results
Definition cp_result_methods.F:43
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2146
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:693
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::test_for_result
logical function, public test_for_result(results, description)
test for a certain result in the result_list
Definition cp_result_methods.F:127
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
f77_interface
interface to use cp2k as library
Definition f77_interface.F:20
f77_interface::f_env_add_defaults
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
Definition f77_interface.F:452
f77_interface::f_env_rm_defaults
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Definition f77_interface.F:496
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
grrm_utils
GRRM interface.
Definition grrm_utils.F:12
grrm_utils::write_grrm
subroutine, public write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, fixed_atoms)
Write GRRM interface file.
Definition grrm_utils.F:50
header::vib_header
subroutine, public vib_header(iw, nr, np)
...
Definition header.F:378
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_rep_blocked
integer, parameter, public do_rep_blocked
Definition input_constants.F:92
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:373
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mode_selective
Module performing a mdoe selective vibrational analysis.
Definition mode_selective.F:20
mode_selective::ms_vb_anal
subroutine, public ms_vb_anal(input, rep_env, para_env, globenv, particles, nrep, calc_intens, dx, output_unit, logger)
Module performing a vibrational analysis.
Definition mode_selective.F:103
molden_utils
Functions handling the MOLDEN format. Split from mode_selective.
Definition molden_utils.F:12
molden_utils::write_vibrations_molden
subroutine, public write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, dump_only_positive, logger, list)
writes the output for vibrational analysis in MOLDEN format
Definition molden_utils.F:399
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
motion_utils
Output Utilities for MOTION_SECTION.
Definition motion_utils.F:13
motion_utils::thrs_motion
real(kind=dp), parameter, public thrs_motion
Definition motion_utils.F:60
motion_utils::rot_ana
subroutine, public rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weighted, natoms, rot_dof, inertia)
Performs an analysis of the principal inertia axis Getting back the generators of the translating and...
Definition motion_utils.F:81
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_particle_matrix
subroutine, public write_particle_matrix(matrix, particle_set, iw, el_per_part, ilist, parts_per_line)
...
Definition particle_methods.F:716
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition physcon.F:129
physcon::a_bohr
real(kind=dp), parameter, public a_bohr
Definition physcon.F:136
physcon::vibfac
real(kind=dp), parameter, public vibfac
Definition physcon.F:189
physcon::n_avogadro
real(kind=dp), parameter, public n_avogadro
Definition physcon.F:126
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
physcon::h_bar
real(kind=dp), parameter, public h_bar
Definition physcon.F:103
physcon::hertz
real(kind=dp), parameter, public hertz
Definition physcon.F:186
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
physcon::e_mass
real(kind=dp), parameter, public e_mass
Definition physcon.F:109
physcon::wavenumbers
real(kind=dp), parameter, public wavenumbers
Definition physcon.F:192
physcon::massunit
real(kind=dp), parameter, public massunit
Definition physcon.F:141
physcon::kjmol
real(kind=dp), parameter, public kjmol
Definition physcon.F:168
physcon::pascal
real(kind=dp), parameter, public pascal
Definition physcon.F:174
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
physcon::debye
real(kind=dp), parameter, public debye
Definition physcon.F:201
replica_methods
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
Definition replica_methods.F:16
replica_methods::rep_env_create
subroutine, public rep_env_create(rep_env, para_env, input, input_declaration, nrep, prep, sync_v, keep_wf_history, row_force)
creates a replica environment together with its force environment
Definition replica_methods.F:91
replica_methods::rep_env_calc_e_f
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
Definition replica_methods.F:366
replica_types
types used to handle many replica of the same system that differ only in atom positions,...
Definition replica_types.F:21
replica_types::rep_env_release
subroutine, public rep_env_release(rep_env)
releases the given replica environment
Definition replica_types.F:138
scine_utils
SCINE interface.
Definition scine_utils.F:12
scine_utils::write_scine
subroutine, public write_scine(iounit, force_env, particles, energy, hessian)
Write SCINE interface file.
Definition scine_utils.F:46
util
All kind of helpful little routines.
Definition util.F:14
vibrational_analysis
Module performing a vibrational analysis.
Definition vibrational_analysis.F:20
vibrational_analysis::vb_anal
subroutine, public vb_anal(input, input_declaration, para_env, globenv)
Module performing a vibrational analysis.
Definition vibrational_analysis.F:102
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
f77_interface::f_env_type
Definition f77_interface.F:163
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:149
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
replica_types::replica_env_type
keeps replicated information about the replicas
Definition replica_types.F:96