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