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