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