(git:e7e05ae)
mode_selective.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 mdoe selective 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 Florian Schiffmann 08.2006
19 ! **************************************************************************************************
20 MODULE mode_selective
21  USE cp_files, ONLY: close_file,&
22  open_file
23  USE cp_log_handling, ONLY: cp_get_default_logger,&
24  cp_logger_get_default_io_unit,&
25  cp_logger_type
26  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
27  cp_print_key_unit_nr
28  USE cp_result_methods, ONLY: get_results
29  USE global_types, ONLY: global_environment_type
30  USE input_constants, ONLY: ms_guess_atomic,&
31  ms_guess_bfgs,&
32  ms_guess_molden,&
33  ms_guess_restart,&
34  ms_guess_restart_vec
35  USE input_section_types, ONLY: section_vals_get,&
36  section_vals_get_subs_vals,&
37  section_vals_type,&
38  section_vals_val_get
39  USE kinds, ONLY: default_path_length,&
40  default_string_length,&
41  dp,&
42  max_line_length
43  USE mathlib, ONLY: diamat_all
44  USE message_passing, ONLY: mp_para_env_type
45  USE molden_utils, ONLY: write_vibrations_molden
46  USE particle_types, ONLY: particle_type
47  USE physcon, ONLY: bohr,&
48  debye,&
49  massunit,&
50  vibfac
51  USE replica_methods, ONLY: rep_env_calc_e_f
52  USE replica_types, ONLY: replica_env_type
53  USE util, ONLY: sort
54 #include "./base/base_uses.f90"
55 
56  IMPLICIT NONE
57 
58  PRIVATE
59  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mode_selective'
60  LOGICAL, PARAMETER :: debug_this_module = .false.
61 
62  TYPE ms_vib_type
63  INTEGER :: mat_size
64  INTEGER :: select_id
65  INTEGER, DIMENSION(:), POINTER :: inv_atoms
66  REAL(KIND=dp) :: eps(2)
67  REAL(KIND=dp) :: sel_freq
68  REAL(KIND=dp) :: low_freq
69  REAL(KIND=dp), POINTER, DIMENSION(:, :) :: b_vec
70  REAL(KIND=dp), POINTER, DIMENSION(:, :) :: delta_vec
71  REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ms_force
72  REAL(KIND=dp), DIMENSION(:), POINTER :: eig_bfgs
73  REAL(KIND=dp), DIMENSION(:), POINTER :: f_range
74  REAL(KIND=dp), DIMENSION(:), POINTER :: inv_range
75  REAL(KIND=dp), POINTER, DIMENSION(:) :: step_b
76  REAL(KIND=dp), POINTER, DIMENSION(:) :: step_r
77  REAL(KIND=dp), DIMENSION(:, :), POINTER :: b_mat
78  REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_deriv
79  REAL(KIND=dp), DIMENSION(:, :), POINTER :: hes_bfgs
80  REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_mat
81  INTEGER :: initial_guess
82  END TYPE ms_vib_type
83 
84  PUBLIC :: ms_vb_anal
85 
86 CONTAINS
87 ! **************************************************************************************************
88 !> \brief Module performing a vibrational analysis
89 !> \param input ...
90 !> \param rep_env ...
91 !> \param para_env ...
92 !> \param globenv ...
93 !> \param particles ...
94 !> \param nrep ...
95 !> \param calc_intens ...
96 !> \param dx ...
97 !> \param output_unit ...
98 !> \param logger ...
99 !> \author Teodoro Laino 08.2006
100 ! **************************************************************************************************
101  SUBROUTINE ms_vb_anal(input, rep_env, para_env, globenv, particles, &
102  nrep, calc_intens, dx, output_unit, logger)
103  TYPE(section_vals_type), POINTER :: input
104  TYPE(replica_env_type), POINTER :: rep_env
105  TYPE(mp_para_env_type), POINTER :: para_env
106  TYPE(global_environment_type), POINTER :: globenv
107  TYPE(particle_type), DIMENSION(:), POINTER :: particles
108  INTEGER :: nrep
109  LOGICAL :: calc_intens
110  REAL(kind=dp) :: dx
111  INTEGER :: output_unit
112  TYPE(cp_logger_type), POINTER :: logger
113 
114  CHARACTER(len=*), PARAMETER :: routinen = 'ms_vb_anal'
115 
116  CHARACTER(LEN=default_string_length) :: description
117  INTEGER :: handle, i, ip1, j, natoms, ncoord
118  LOGICAL :: converged
119  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mass, pos0
120  REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp_deriv
121  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tmp_dip
122  TYPE(ms_vib_type) :: ms_vib
123 
124  CALL timeset(routinen, handle)
125  converged = .false.
126  natoms = SIZE(particles)
127  ncoord = 3*natoms
128  ALLOCATE (mass(3*natoms))
129  DO i = 1, natoms
130  DO j = 1, 3
131  mass((i - 1)*3 + j) = particles(i)%atomic_kind%mass
132  mass((i - 1)*3 + j) = sqrt(mass((i - 1)*3 + j))
133  END DO
134  END DO
135  ! Allocate working arrays
136  ALLOCATE (ms_vib%delta_vec(ncoord, nrep))
137  ALLOCATE (ms_vib%b_vec(ncoord, nrep))
138  ALLOCATE (ms_vib%step_r(nrep))
139  ALLOCATE (ms_vib%step_b(nrep))
140  IF (calc_intens) THEN
141  description = '[DIPOLE]'
142  ALLOCATE (tmp_dip(nrep, 3, 2))
143  ALLOCATE (ms_vib%dip_deriv(3, nrep))
144  END IF
145  CALL ms_initial_moves(para_env, nrep, input, globenv, ms_vib, &
146  particles, &
147  mass, &
148  dx, &
149  calc_intens, logger)
150  ncoord = 3*natoms
151  ALLOCATE (pos0(ncoord))
152  ALLOCATE (ms_vib%ms_force(ncoord, nrep))
153  DO i = 1, natoms
154  DO j = 1, 3
155  pos0((i - 1)*3 + j) = particles((i))%r(j)
156  END DO
157  END DO
158  ncoord = 3*natoms
159  DO
160  ms_vib%ms_force = huge(0.0_dp)
161  DO i = 1, nrep
162  DO j = 1, ncoord
163  rep_env%r(j, i) = pos0(j) + ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
164  END DO
165  END DO
166  CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
167 
168  DO i = 1, nrep
169  IF (calc_intens) THEN
170  CALL get_results(results=rep_env%results(i)%results, &
171  description=description, &
172  n_rep=ip1)
173  CALL get_results(results=rep_env%results(i)%results, &
174  description=description, &
175  values=tmp_dip(i, :, 1), &
176  nval=ip1)
177  END IF
178  DO j = 1, ncoord
179  ms_vib%ms_force(j, i) = rep_env%f(j, i)
180  END DO
181  END DO
182  DO i = 1, nrep
183  DO j = 1, ncoord
184  rep_env%r(j, i) = pos0(j) - ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
185  END DO
186  END DO
187  CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
188  IF (calc_intens) THEN
189  DO i = 1, nrep
190  CALL get_results(results=rep_env%results(i)%results, &
191  description=description, &
192  n_rep=ip1)
193  CALL get_results(results=rep_env%results(i)%results, &
194  description=description, &
195  values=tmp_dip(i, :, 2), &
196  nval=ip1)
197  ms_vib%dip_deriv(:, ms_vib%mat_size + i) = (tmp_dip(i, :, 1) - tmp_dip(i, :, 2))/(2*ms_vib%step_b(i))
198  END DO
199  END IF
200 
201  CALL evaluate_h_update_b(rep_env, ms_vib, input, nrep, &
202  particles, &
203  mass, &
204  converged, &
205  dx, calc_intens, &
206  output_unit, logger)
207  IF (converged) EXIT
208  IF (calc_intens) THEN
209  ALLOCATE (tmp_deriv(3, ms_vib%mat_size))
210  tmp_deriv = ms_vib%dip_deriv
211  DEALLOCATE (ms_vib%dip_deriv)
212  ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
213  ms_vib%dip_deriv(:, 1:ms_vib%mat_size) = tmp_deriv(:, 1:ms_vib%mat_size)
214  DEALLOCATE (tmp_deriv)
215  END IF
216  END DO
217  DEALLOCATE (ms_vib%ms_force)
218  DEALLOCATE (pos0)
219  DEALLOCATE (ms_vib%step_r)
220  DEALLOCATE (ms_vib%step_b)
221  DEALLOCATE (ms_vib%b_vec)
222  DEALLOCATE (ms_vib%delta_vec)
223  DEALLOCATE (mass)
224  DEALLOCATE (ms_vib%b_mat)
225  DEALLOCATE (ms_vib%s_mat)
226  IF (ms_vib%select_id == 3) THEN
227  DEALLOCATE (ms_vib%inv_atoms)
228  END IF
229  IF (ASSOCIATED(ms_vib%eig_bfgs)) THEN
230  DEALLOCATE (ms_vib%eig_bfgs)
231  END IF
232  IF (ASSOCIATED(ms_vib%hes_bfgs)) THEN
233  DEALLOCATE (ms_vib%hes_bfgs)
234  END IF
235  IF (calc_intens) THEN
236  DEALLOCATE (ms_vib%dip_deriv)
237  DEALLOCATE (tmp_dip)
238  END IF
239  CALL timestop(handle)
240  END SUBROUTINE ms_vb_anal
241 ! **************************************************************************************************
242 !> \brief Generates the first displacement vector for a mode selctive vibrational
243 !> analysis. At the moment this is a random number for selected atoms
244 !> \param para_env ...
245 !> \param nrep ...
246 !> \param input ...
247 !> \param globenv ...
248 !> \param ms_vib ...
249 !> \param particles ...
250 !> \param mass ...
251 !> \param dx ...
252 !> \param calc_intens ...
253 !> \param logger ...
254 !> \author Florian Schiffmann 11.2007
255 ! **************************************************************************************************
256  SUBROUTINE ms_initial_moves(para_env, nrep, input, globenv, ms_vib, particles, &
257  mass, dx, &
258  calc_intens, logger)
259  TYPE(mp_para_env_type), POINTER :: para_env
260  INTEGER :: nrep
261  TYPE(section_vals_type), POINTER :: input
262  TYPE(global_environment_type), POINTER :: globenv
263  TYPE(ms_vib_type) :: ms_vib
264  TYPE(particle_type), DIMENSION(:), POINTER :: particles
265  REAL(kind=dp), DIMENSION(:) :: mass
266  REAL(kind=dp) :: dx
267  LOGICAL :: calc_intens
268  TYPE(cp_logger_type), POINTER :: logger
269 
270  CHARACTER(len=*), PARAMETER :: routinen = 'MS_initial_moves'
271 
272  INTEGER :: guess, handle, i, j, jj, k, m, &
273  n_rep_val, natoms, ncoord
274  INTEGER, ALLOCATABLE, DIMENSION(:) :: map_atoms
275  INTEGER, DIMENSION(:), POINTER :: tmplist
276  LOGICAL :: do_involved_atoms, ionode
277  REAL(kind=dp) :: my_val, norm
278  TYPE(section_vals_type), POINTER :: involved_at_section, ms_vib_section
279 
280  CALL timeset(routinen, handle)
281  NULLIFY (ms_vib%eig_bfgs, ms_vib%f_range, ms_vib%hes_bfgs, ms_vib%inv_range)
282  ms_vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
283  CALL section_vals_val_get(ms_vib_section, "INITIAL_GUESS", i_val=guess)
284  CALL section_vals_val_get(ms_vib_section, "EPS_MAX_VAL", r_val=ms_vib%eps(1))
285  CALL section_vals_val_get(ms_vib_section, "EPS_NORM", r_val=ms_vib%eps(2))
286  CALL section_vals_val_get(ms_vib_section, "RANGE", n_rep_val=n_rep_val)
287  ms_vib%select_id = 0
288  IF (n_rep_val .NE. 0) THEN
289  CALL section_vals_val_get(ms_vib_section, "RANGE", r_vals=ms_vib%f_range)
290  IF (ms_vib%f_range(1) .GT. ms_vib%f_range(2)) THEN
291  my_val = ms_vib%f_range(2)
292  ms_vib%f_range(2) = ms_vib%f_range(1)
293  ms_vib%f_range(1) = my_val
294  END IF
295  ms_vib%select_id = 2
296  END IF
297  CALL section_vals_val_get(ms_vib_section, "FREQUENCY", r_val=ms_vib%sel_freq)
298  CALL section_vals_val_get(ms_vib_section, "LOWEST_FREQUENCY", r_val=ms_vib%low_freq)
299  IF (ms_vib%sel_freq .GT. 0._dp) ms_vib%select_id = 1
300  involved_at_section => section_vals_get_subs_vals(ms_vib_section, "INVOLVED_ATOMS")
301  CALL section_vals_get(involved_at_section, explicit=do_involved_atoms)
302  IF (do_involved_atoms) THEN
303  CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", n_rep_val=n_rep_val)
304  jj = 0
305  DO k = 1, n_rep_val
306  CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=k, i_vals=tmplist)
307  DO j = 1, SIZE(tmplist)
308  jj = jj + 1
309  END DO
310  END DO
311  IF (jj .GE. 1) THEN
312  natoms = jj
313  ALLOCATE (ms_vib%inv_atoms(natoms))
314  jj = 0
315  DO m = 1, n_rep_val
316  CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=m, i_vals=tmplist)
317  DO j = 1, SIZE(tmplist)
318  ms_vib%inv_atoms(j) = tmplist(j)
319  END DO
320  END DO
321  ms_vib%select_id = 3
322  END IF
323  CALL section_vals_val_get(involved_at_section, "RANGE", n_rep_val=n_rep_val)
324  IF (n_rep_val .NE. 0) THEN
325  CALL section_vals_val_get(involved_at_section, "RANGE", r_vals=ms_vib%inv_range)
326  IF (ms_vib%inv_range(1) .GT. ms_vib%inv_range(2)) THEN
327  ms_vib%inv_range(2) = my_val
328  ms_vib%inv_range(2) = ms_vib%inv_range(1)
329  ms_vib%inv_range(1) = my_val
330  END IF
331  END IF
332  END IF
333  IF (ms_vib%select_id == 0) &
334  cpabort("no frequency, range or involved atoms specified ")
335  ionode = para_env%is_source()
336  SELECT CASE (guess)
337  CASE (ms_guess_atomic)
338  ms_vib%initial_guess = 1
339  CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
340  jj = 0
341  DO k = 1, n_rep_val
342  CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
343  DO j = 1, SIZE(tmplist)
344  jj = jj + 1
345  END DO
346  END DO
347  IF (jj < 1) THEN
348  natoms = SIZE(particles)
349  ALLOCATE (map_atoms(natoms))
350  DO j = 1, natoms
351  map_atoms(j) = j
352  END DO
353  ELSE
354  natoms = jj
355  ALLOCATE (map_atoms(natoms))
356  jj = 0
357  DO m = 1, n_rep_val
358  CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
359  DO j = 1, SIZE(tmplist)
360  map_atoms(j) = tmplist(j)
361  END DO
362  END DO
363  END IF
364 
365  ! apply random displacement along the mass weighted nuclear cartesian coordinates
366  ms_vib%b_vec = 0._dp
367  ms_vib%delta_vec = 0._dp
368  jj = 0
369 
370  DO i = 1, nrep
371  DO j = 1, natoms
372  DO k = 1, 3
373  jj = (map_atoms(j) - 1)*3 + k
374  ms_vib%b_vec(jj, i) = abs(globenv%gaussian_rng_stream%next())
375  END DO
376  END DO
377  norm = sqrt(dot_product(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
378  ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
379  END DO
380 
381  IF (nrep .GT. 1) THEN
382  DO k = 1, 10
383  DO j = 1, nrep
384  DO i = 1, nrep
385  IF (i .NE. j) THEN
386  ms_vib%b_vec(:, j) = &
387  ms_vib%b_vec(:, j) - dot_product(ms_vib%b_vec(:, j), ms_vib%b_vec(:, i))*ms_vib%b_vec(:, i)
388  ms_vib%b_vec(:, j) = &
389  ms_vib%b_vec(:, j)/sqrt(dot_product(ms_vib%b_vec(:, j), ms_vib%b_vec(:, j)))
390  END IF
391  END DO
392  END DO
393  END DO
394  END IF
395 
396  ms_vib%mat_size = 0
397  DO i = 1, SIZE(ms_vib%b_vec, 1)
398  ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
399  END DO
400  CASE (ms_guess_bfgs)
401 
402  ms_vib%initial_guess = 2
403  CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
404  ms_vib%mat_size = 0
405 
406  CASE (ms_guess_restart_vec)
407 
408  ms_vib%initial_guess = 3
409  ncoord = 3*SIZE(particles)
410  CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
411 
412  ms_vib%mat_size = 0
413  CASE (ms_guess_restart)
414  ms_vib%initial_guess = 4
415  ncoord = 3*SIZE(particles)
416  CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
417 
418  CASE (ms_guess_molden)
419  ms_vib%initial_guess = 5
420  ncoord = 3*SIZE(particles)
421  CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
422  ms_vib%mat_size = 0
423  END SELECT
424  CALL para_env%bcast(ms_vib%b_vec)
425  CALL para_env%bcast(ms_vib%delta_vec)
426  DO i = 1, nrep
427  ms_vib%step_r(i) = dx/sqrt(dot_product(ms_vib%delta_vec(:, i), ms_vib%delta_vec(:, i)))
428  ms_vib%step_b(i) = sqrt(dot_product(ms_vib%step_r(i)*ms_vib%b_vec(:, i), ms_vib%step_r(i)*ms_vib%b_vec(:, i)))
429  END DO
430  CALL timestop(handle)
431 
432  END SUBROUTINE ms_initial_moves
433 
434 ! **************************************************************************************************
435 !> \brief ...
436 !> \param ms_vib_section ...
437 !> \param ms_vib ...
438 !> \param particles ...
439 !> \param mass ...
440 !> \param para_env ...
441 !> \param nrep ...
442 !> \author Florian Schiffmann 11.2007
443 ! **************************************************************************************************
444  SUBROUTINE bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
445 
446  TYPE(section_vals_type), POINTER :: ms_vib_section
447  TYPE(ms_vib_type) :: ms_vib
448  TYPE(particle_type), DIMENSION(:), POINTER :: particles
449  REAL(kind=dp), DIMENSION(:) :: mass
450  TYPE(mp_para_env_type), POINTER :: para_env
451  INTEGER :: nrep
452 
453  CHARACTER(LEN=default_path_length) :: hes_filename
454  INTEGER :: hesunit, i, j, jj, k, natoms, ncoord, &
455  output_unit, stat
456  INTEGER, DIMENSION(:), POINTER :: tmplist
457  REAL(kind=dp) :: my_val, norm
458  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
459  TYPE(cp_logger_type), POINTER :: logger
460 
461  logger => cp_get_default_logger()
462  output_unit = cp_logger_get_default_io_unit(logger)
463 
464  natoms = SIZE(particles)
465  ncoord = 3*natoms
466 
467  ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
468  ALLOCATE (ms_vib%eig_bfgs(ncoord))
469 
470  IF (para_env%is_source()) THEN
471  CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
472  IF (hes_filename == "") hes_filename = "HESSIAN"
473  CALL open_file(file_name=hes_filename, file_status="OLD", &
474  file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
475  ALLOCATE (tmp(ncoord))
476  ALLOCATE (tmplist(ncoord))
477 
478  ! should use the cp_fm_read_unformatted...
479  DO i = 1, ncoord
480  READ (unit=hesunit, iostat=stat) ms_vib%hes_bfgs(:, i)
481  END DO
482  CALL close_file(hesunit)
483  IF (output_unit > 0) THEN
484  IF (stat /= 0) THEN
485  WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading HESSIAN **"
486  ELSE
487  WRITE (output_unit, fmt="(/,T2,A)") &
488  "*** Initial Hessian has been read successfully ***"
489  END IF
490  END IF
491  DO i = 1, ncoord
492  DO j = 1, ncoord
493  ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
494  END DO
495  END DO
496 
497  CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
498  tmp(:) = 0._dp
499  IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
500  IF (ms_vib%select_id == 2) my_val = (((ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp)/vibfac)**2/massunit
501  IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
502  DO i = 1, ncoord
503  tmp(i) = abs(my_val - ms_vib%eig_bfgs(i))
504  END DO
505  ELSE IF (ms_vib%select_id == 3) THEN
506  DO i = 1, ncoord
507  DO j = 1, SIZE(ms_vib%inv_atoms)
508  DO k = 1, 3
509  jj = (ms_vib%inv_atoms(j) - 1)*3 + k
510  tmp(i) = tmp(i) + sqrt(ms_vib%hes_bfgs(jj, i)**2)
511  END DO
512  END DO
513  IF ((sign(1._dp, ms_vib%eig_bfgs(i))*sqrt(abs(ms_vib%eig_bfgs(i))*massunit)*vibfac) .LE. 400._dp) tmp(i) = 0._dp
514  END DO
515  tmp(:) = -tmp(:)
516  END IF
517  CALL sort(tmp, ncoord, tmplist)
518  DO i = 1, nrep
519  ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
520  norm = sqrt(dot_product(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
521  ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
522  END DO
523  DO i = 1, SIZE(ms_vib%b_vec, 1)
524  ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
525  END DO
526  DEALLOCATE (tmp)
527  DEALLOCATE (tmplist)
528  END IF
529 
530  CALL para_env%bcast(ms_vib%b_vec)
531  CALL para_env%bcast(ms_vib%delta_vec)
532 
533  DEALLOCATE (ms_vib%hes_bfgs)
534  DEALLOCATE (ms_vib%eig_bfgs)
535  ms_vib%mat_size = 0
536 
537  END SUBROUTINE bfgs_guess
538 
539 ! **************************************************************************************************
540 !> \brief ...
541 !> \param ms_vib_section ...
542 !> \param para_env ...
543 !> \param ms_vib ...
544 !> \param mass ...
545 !> \param ionode ...
546 !> \param particles ...
547 !> \param nrep ...
548 !> \param calc_intens ...
549 !> \author Florian Schiffmann 11.2007
550 ! **************************************************************************************************
551  SUBROUTINE rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
552 
553  TYPE(section_vals_type), POINTER :: ms_vib_section
554  TYPE(mp_para_env_type), POINTER :: para_env
555  TYPE(ms_vib_type) :: ms_vib
556  REAL(kind=dp), DIMENSION(:) :: mass
557  LOGICAL :: ionode
558  TYPE(particle_type), DIMENSION(:), POINTER :: particles
559  INTEGER :: nrep
560  LOGICAL :: calc_intens
561 
562  CHARACTER(LEN=default_path_length) :: ms_filename
563  INTEGER :: hesunit, i, j, mat, natoms, ncoord, &
564  output_unit, stat, statint
565  INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
566  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
567  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h
568  TYPE(cp_logger_type), POINTER :: logger
569 
570  logger => cp_get_default_logger()
571  output_unit = cp_logger_get_default_io_unit(logger)
572 
573  natoms = SIZE(particles)
574  ncoord = 3*natoms
575  IF (calc_intens) THEN
576  DEALLOCATE (ms_vib%dip_deriv)
577  END IF
578 
579  IF (ionode) THEN
580 
581  CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
582  IF (ms_filename == "") ms_filename = "MS_RESTART"
583  CALL open_file(file_name=ms_filename, &
584  file_status="UNKNOWN", &
585  file_form="UNFORMATTED", &
586  file_action="READ", &
587  unit_number=hesunit)
588  READ (unit=hesunit, iostat=stat) mat
589  cpassert(stat == 0)
590  ms_vib%mat_size = mat
591  END IF
592  CALL para_env%bcast(ms_vib%mat_size)
593  ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
594  ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
595  IF (calc_intens) THEN
596  ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
597  END IF
598  IF (ionode) THEN
599  statint = 0
600  READ (unit=hesunit, iostat=stat) ms_vib%b_mat
601  READ (unit=hesunit, iostat=stat) ms_vib%s_mat
602  IF (calc_intens) READ (unit=hesunit, iostat=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
603  IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART,", &
604  "intensities are requested but not present in restart file **"
605  CALL close_file(hesunit)
606  IF (output_unit > 0) THEN
607  IF (stat /= 0) THEN
608  WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART **"
609  ELSE
610  WRITE (output_unit, fmt="(/,T2,A)") "*** RESTART has been read successfully ***"
611  END IF
612  END IF
613  END IF
614  CALL para_env%bcast(ms_vib%b_mat)
615  CALL para_env%bcast(ms_vib%s_mat)
616  IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
617  ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
618  ALLOCATE (eigenval(ms_vib%mat_size))
619  ALLOCATE (ind(ms_vib%mat_size))
620 
621  CALL dgemm('T', 'N', ms_vib%mat_size, ms_vib%mat_size, SIZE(ms_vib%s_mat, 1), 1._dp, ms_vib%b_mat, SIZE(ms_vib%b_mat, 1), &
622  ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
623  CALL diamat_all(approx_h, eigenval)
624 
625  CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, ms_vib%b_vec)
626  IF (ms_vib%initial_guess .NE. 4) THEN
627 
628  ms_vib%b_vec = 0._dp
629  DO i = 1, nrep
630  DO j = 1, ms_vib%mat_size
631  ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_h(j, ind(i))*ms_vib%b_mat(:, j)
632  END DO
633  ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/sqrt(dot_product(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
634  END DO
635 
636  DEALLOCATE (ms_vib%s_mat)
637  DEALLOCATE (ms_vib%b_mat)
638  IF (calc_intens) THEN
639  DEALLOCATE (ms_vib%dip_deriv)
640  ALLOCATE (ms_vib%dip_deriv(3, nrep))
641  END IF
642  END IF
643  DEALLOCATE (approx_h)
644  DEALLOCATE (eigenval)
645  DEALLOCATE (ind)
646  DO i = 1, nrep
647  ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
648  END DO
649 
650  END SUBROUTINE rest_guess
651 
652 ! **************************************************************************************************
653 !> \brief ...
654 !> \param ms_vib_section ...
655 !> \param input ...
656 !> \param para_env ...
657 !> \param ms_vib ...
658 !> \param mass ...
659 !> \param ncoord ...
660 !> \param nrep ...
661 !> \param logger ...
662 !> \author Florian Schiffmann 11.2007
663 ! **************************************************************************************************
664  SUBROUTINE molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
665  TYPE(section_vals_type), POINTER :: ms_vib_section, input
666  TYPE(mp_para_env_type), POINTER :: para_env
667  TYPE(ms_vib_type) :: ms_vib
668  REAL(kind=dp), DIMENSION(:) :: mass
669  INTEGER :: ncoord, nrep
670  TYPE(cp_logger_type), POINTER :: logger
671 
672  CHARACTER(LEN=2) :: at_name
673  CHARACTER(LEN=default_path_length) :: ms_filename
674  CHARACTER(LEN=max_line_length) :: info
675  INTEGER :: i, istat, iw, j, jj, k, nvibs, &
676  output_molden, output_unit, stat
677  INTEGER, DIMENSION(:), POINTER :: tmplist
678  LOGICAL :: reading_vib
679  REAL(kind=dp) :: my_val, norm
680  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
681  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
682  REAL(kind=dp), DIMENSION(3, ncoord/3) :: pos
683 
684  output_unit = cp_logger_get_default_io_unit(logger)
685 
686  CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
687  IF (ms_filename == "") output_molden = &
688  cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
689  extension=".mol", file_status='UNKNOWN', &
690  file_action="READ")
691  IF (para_env%is_source()) THEN
692 
693  IF (ms_filename == "") THEN
694  iw = output_molden
695  ELSE
696  CALL open_file(file_name=trim(ms_filename), &
697  file_status="UNKNOWN", &
698  file_form="FORMATTED", &
699  file_action="READ", &
700  unit_number=iw)
701  END IF
702  info = ""
703  READ (iw, *, iostat=stat) info
704  READ (iw, *, iostat=stat) info
705  istat = 0
706  nvibs = 0
707  reading_vib = .false.
708  DO
709  READ (iw, *, iostat=stat) info
710  istat = istat + stat
711  IF (trim(adjustl(info)) == "[FR-COORD]") EXIT
712 
713  cpassert(stat == 0)
714 
715  IF (reading_vib) nvibs = nvibs + 1
716  IF (trim(adjustl(info)) == "[FREQ]") reading_vib = .true.
717  END DO
718  rewind(iw)
719  istat = 0
720  READ (iw, *, iostat=stat) info
721  istat = istat + stat
722  READ (iw, *, iostat=stat) info
723  istat = istat + stat
724  ! Skip [Atoms] section
725  DO
726  READ (iw, *, iostat=stat) info
727  istat = istat + stat
728  cpassert(stat == 0)
729  IF (trim(adjustl(info)) == "[FREQ]") EXIT
730  END DO
731  ! Read frequencies and modes
732  ALLOCATE (freq(nvibs))
733  ALLOCATE (modes(ncoord, nvibs))
734 
735  DO i = 1, nvibs
736  READ (iw, *, iostat=stat) freq(i)
737  istat = istat + stat
738  END DO
739  READ (iw, *) info
740  DO i = 1, ncoord/3
741  READ (iw, *, iostat=stat) at_name, pos(:, i)
742  istat = istat + stat
743  END DO
744  READ (iw, *) info
745  DO i = 1, nvibs
746  READ (iw, *) info
747  istat = istat + stat
748  DO j = 1, ncoord/3
749  k = (j - 1)*3 + 1
750  READ (iw, *, iostat=stat) modes(k:k + 2, i)
751  istat = istat + stat
752  END DO
753  END DO
754  IF (ms_filename .NE. "") CALL close_file(iw)
755  IF (output_unit > 0) THEN
756  IF (istat /= 0) THEN
757  WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MOLDEN file **"
758  ELSE
759  WRITE (output_unit, fmt="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
760  END IF
761  END IF
762  !!!!!!! select modes !!!!!!
763  ALLOCATE (tmp(nvibs))
764  tmp(:) = 0.0_dp
765  ALLOCATE (tmplist(nvibs))
766  IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
767  IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
768  IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
769  DO i = 1, nvibs
770  tmp(i) = abs(my_val - freq(i))
771  END DO
772  ELSE IF (ms_vib%select_id == 3) THEN
773  DO i = 1, nvibs
774  DO j = 1, SIZE(ms_vib%inv_atoms)
775  DO k = 1, 3
776  jj = (ms_vib%inv_atoms(j) - 1)*3 + k
777  tmp(i) = tmp(i) + sqrt(modes(jj, i)**2)
778  END DO
779  END DO
780  IF (freq(i) .LE. 400._dp) tmp(i) = 0._dp
781  END DO
782  tmp(:) = -tmp(:)
783  END IF
784  CALL sort(tmp, nvibs, tmplist)
785  DO i = 1, nrep
786  ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
787  norm = sqrt(dot_product(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
788  ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
789  END DO
790  DO i = 1, nrep
791  ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
792  END DO
793 
794  DEALLOCATE (freq)
795  DEALLOCATE (modes)
796  DEALLOCATE (tmp)
797  DEALLOCATE (tmplist)
798 
799  END IF
800  CALL para_env%bcast(ms_vib%b_vec)
801  CALL para_env%bcast(ms_vib%delta_vec)
802 
803  IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
804  "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
805  END SUBROUTINE molden_guess
806 
807 ! **************************************************************************************************
808 !> \brief Davidson algorithm for to generate a approximate Hessian for mode
809 !> selective vibrational analysis
810 !> \param rep_env ...
811 !> \param ms_vib ...
812 !> \param input ...
813 !> \param nrep ...
814 !> \param particles ...
815 !> \param mass ...
816 !> \param converged ...
817 !> \param dx ...
818 !> \param calc_intens ...
819 !> \param output_unit_ms ...
820 !> \param logger ...
821 !> \author Florian Schiffmann 11.2007
822 ! **************************************************************************************************
823  SUBROUTINE evaluate_h_update_b(rep_env, ms_vib, input, nrep, &
824  particles, &
825  mass, &
826  converged, dx, &
827  calc_intens, output_unit_ms, logger)
828  TYPE(replica_env_type), POINTER :: rep_env
829  TYPE(ms_vib_type) :: ms_vib
830  TYPE(section_vals_type), POINTER :: input
831  INTEGER :: nrep
832  TYPE(particle_type), DIMENSION(:), POINTER :: particles
833  REAL(kind=dp), DIMENSION(:) :: mass
834  LOGICAL :: converged
835  REAL(kind=dp) :: dx
836  LOGICAL :: calc_intens
837  INTEGER :: output_unit_ms
838  TYPE(cp_logger_type), POINTER :: logger
839 
840  INTEGER :: i, j, jj, k, natoms, ncoord
841  INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
842  LOGICAL :: dump_only_positive
843  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
844  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h, h_save, residuum, tmp_b, tmp_s
845  REAL(kind=dp), DIMENSION(2, nrep) :: criteria
846  REAL(kind=dp), DIMENSION(:), POINTER :: intensities
847 
848  natoms = SIZE(particles)
849  ncoord = 3*natoms
850  nrep = SIZE(rep_env%f, 2)
851 
852  !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
853  IF (ms_vib%mat_size .NE. 0) THEN
854 
855  ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
856  ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
857 
858  tmp_b(:, :) = ms_vib%b_mat
859  tmp_s(:, :) = ms_vib%s_mat
860 
861  DEALLOCATE (ms_vib%b_mat)
862  DEALLOCATE (ms_vib%s_mat)
863  END IF
864 
865  ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
866  ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
867 
868  ms_vib%s_mat = 0.0_dp
869 
870  DO i = 1, 3*natoms
871  IF (ms_vib%mat_size .NE. 0) THEN
872  DO j = 1, ms_vib%mat_size
873  ms_vib%b_mat(i, j) = tmp_b(i, j)
874  ms_vib%s_mat(i, j) = tmp_s(i, j)
875  END DO
876  END IF
877  DO j = 1, nrep
878  ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
879  END DO
880  END DO
881 
882  IF (ms_vib%mat_size .NE. 0) THEN
883  DEALLOCATE (tmp_s)
884  DEALLOCATE (tmp_b)
885  END IF
886 
887  ms_vib%mat_size = ms_vib%mat_size + nrep
888 
889  ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
890  ALLOCATE (h_save(ms_vib%mat_size, ms_vib%mat_size))
891  ALLOCATE (eigenval(ms_vib%mat_size))
892 
893  !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
894 
895  DO i = 1, nrep
896  DO j = 1, 3*natoms
897  ms_vib%s_mat(j, ms_vib%mat_size - nrep + i) = -(ms_vib%ms_force(j, i) - rep_env%f(j, i))/(2*ms_vib%step_b(i)*mass(j))
898  END DO
899  END DO
900 
901  CALL dgemm('T', 'N', ms_vib%mat_size, ms_vib%mat_size, SIZE(ms_vib%s_mat, 1), 1._dp, ms_vib%b_mat, SIZE(ms_vib%b_mat, 1), &
902  ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
903  h_save(:, :) = approx_h
904 
905  CALL diamat_all(approx_h, eigenval)
906 
907  !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
908  ALLOCATE (ind(ms_vib%mat_size))
909  ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
910 
911  CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, residuum, criteria)
912 
913  DO i = 1, nrep
914  DO j = 1, natoms
915  DO k = 1, 3
916  jj = (j - 1)*3 + k
917  ms_vib%delta_vec(jj, i) = ms_vib%b_vec(jj, i)/mass(jj)
918  END DO
919  END DO
920  END DO
921 
922  DO i = 1, nrep
923  ms_vib%step_r(i) = dx/sqrt(dot_product(ms_vib%delta_vec(:, i), ms_vib%delta_vec(:, i)))
924  ms_vib%step_b(i) = sqrt(dot_product(ms_vib%step_r(i)*ms_vib%b_vec(:, i), ms_vib%step_r(i)*ms_vib%b_vec(:, i)))
925  END DO
926  converged = .false.
927  IF (maxval(criteria(1, :)) .LE. ms_vib%eps(1) .AND. maxval(criteria(2, :)) &
928  .LE. ms_vib%eps(2) .OR. ms_vib%mat_size .GE. ncoord) converged = .true.
929  ALLOCATE (freq(nrep))
930  DO i = 1, nrep
931  freq(i) = sqrt(abs(eigenval(ind(i)))*massunit)*vibfac
932  END DO
933 
934  !!! write information and output !!!
935  IF (converged) THEN
936  eigenval(:) = sign(1._dp, eigenval(:))*sqrt(abs(eigenval(:))*massunit)*vibfac
937  ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
938  tmp_b = 0._dp
939  ALLOCATE (tmp_s(3, ms_vib%mat_size))
940  tmp_s = 0._dp
941  IF (calc_intens) THEN
942  ALLOCATE (intensities(ms_vib%mat_size))
943  intensities = 0._dp
944  END IF
945  DO i = 1, ms_vib%mat_size
946  DO j = 1, ms_vib%mat_size
947  tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
948  END DO
949  tmp_b(:, i) = tmp_b(:, i)/sqrt(dot_product(tmp_b(:, i), tmp_b(:, i)))
950  END DO
951  IF (calc_intens) THEN
952  DO i = 1, ms_vib%mat_size
953  DO j = 1, ms_vib%mat_size
954  tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_h(j, i)
955  END DO
956  IF (calc_intens) intensities(i) = sqrt(dot_product(tmp_s(:, i), tmp_s(:, i)))
957  END DO
958  END IF
959  IF (calc_intens) THEN
960  CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
961  input, nrep, approx_h, eigenval, calc_intens, &
962  intensities=intensities, logger=logger)
963  ELSE
964  CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
965  input, nrep, approx_h, eigenval, calc_intens, logger=logger)
966  END IF
967  dump_only_positive = ms_vib%low_freq .GT. 0.0_dp
968  CALL write_vibrations_molden(input, particles, eigenval, tmp_b, intensities, calc_intens, &
969  dump_only_positive=dump_only_positive, logger=logger)
970  IF (calc_intens) THEN
971  DEALLOCATE (intensities)
972  END IF
973  DEALLOCATE (tmp_b)
974  DEALLOCATE (tmp_s)
975  END IF
976 
977  IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
978  ms_vib, input, nrep, approx_h, eigenval, calc_intens, logger=logger)
979 
980  DEALLOCATE (freq)
981  DEALLOCATE (approx_h)
982  DEALLOCATE (eigenval)
983  DEALLOCATE (residuum)
984  DEALLOCATE (ind)
985 
986  END SUBROUTINE evaluate_h_update_b
987 
988 ! **************************************************************************************************
989 !> \brief writes the output for a mode tracking calculation
990 !> \param ms_vib ...
991 !> \param nrep ...
992 !> \param mass ...
993 !> \param ncoord ...
994 !> \param approx_H ...
995 !> \param eigenval ...
996 !> \param ind ...
997 !> \param residuum ...
998 !> \param criteria ...
999 !> \author Florian Schiffmann 11.2007
1000 ! **************************************************************************************************
1001  SUBROUTINE select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1002 
1003  TYPE(ms_vib_type) :: ms_vib
1004  INTEGER :: nrep
1005  REAL(kind=dp), DIMENSION(:) :: mass
1006  INTEGER :: ncoord
1007  REAL(kind=dp), DIMENSION(:, :) :: approx_h
1008  REAL(kind=dp), DIMENSION(:) :: eigenval
1009  INTEGER, DIMENSION(:) :: ind
1010  REAL(kind=dp), DIMENSION(:, :) :: residuum
1011  REAL(kind=dp), DIMENSION(2, nrep), OPTIONAL :: criteria
1012 
1013  INTEGER :: i, j, jj, k
1014  REAL(kind=dp) :: my_val, norm
1015  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1016  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1017 
1018  ALLOCATE (tmp(ms_vib%mat_size))
1019 
1020  SELECT CASE (ms_vib%select_id)
1021  CASE (1)
1022  my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1023  DO i = 1, ms_vib%mat_size
1024  tmp(i) = abs(my_val - eigenval(i))
1025  END DO
1026  CALL sort(tmp, (ms_vib%mat_size), ind)
1027  residuum = 0._dp
1028  DO j = 1, nrep
1029  DO i = 1, ms_vib%mat_size
1030  residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1031  END DO
1032  END DO
1033  CASE (2)
1034  CALL get_vibs_in_range(ms_vib, approx_h, eigenval, residuum, nrep, ind)
1035  CASE (3)
1036 
1037  ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1038  tmp_b = 0._dp
1039 
1040  DO i = 1, ms_vib%mat_size
1041  DO j = 1, ms_vib%mat_size
1042  tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
1043  END DO
1044  tmp_b(:, i) = tmp_b(:, i)/sqrt(dot_product(tmp_b(:, i), tmp_b(:, i)))
1045  END DO
1046  tmp = 0._dp
1047  DO i = 1, ms_vib%mat_size
1048  DO j = 1, SIZE(ms_vib%inv_atoms)
1049  DO k = 1, 3
1050  jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1051  tmp(i) = tmp(i) + sqrt(tmp_b(jj, i)**2)
1052  END DO
1053  END DO
1054  IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1055  IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) .LE. 400._dp) tmp(i) = 0._dp
1056  ELSE
1057  IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) .LE. ms_vib%inv_range(1)) tmp(i) = 0._dp
1058  IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) .GE. ms_vib%inv_range(2)) tmp(i) = 0._dp
1059  END IF
1060  END DO
1061  tmp(:) = -tmp(:)
1062  CALL sort(tmp, (ms_vib%mat_size), ind)
1063  residuum(:, :) = 0._dp
1064 
1065  DO j = 1, nrep
1066  DO i = 1, ms_vib%mat_size
1067  residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1068  END DO
1069  END DO
1070  DEALLOCATE (tmp_b)
1071  END SELECT
1072 
1073  DO j = 1, nrep
1074  DO i = 1, ms_vib%mat_size
1075  residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1076  END DO
1077  END DO
1078  IF (PRESENT(criteria)) THEN
1079  DO i = 1, nrep
1080  criteria(1, i) = maxval((residuum(:, i)))
1081  criteria(2, i) = sqrt(dot_product(residuum(:, i), residuum(:, i)))
1082  END DO
1083  END IF
1084 
1085  DO i = 1, nrep
1086  norm = sqrt(dot_product(residuum(:, i), residuum(:, i)))
1087  residuum(:, i) = residuum(:, i)/norm
1088  END DO
1089 
1090  DO k = 1, 10
1091  DO j = 1, nrep
1092  DO i = 1, ms_vib%mat_size
1093  residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1094  residuum(:, j) = residuum(:, j)/sqrt(dot_product(residuum(:, j), residuum(:, j)))
1095  END DO
1096  IF (nrep .GT. 1) THEN
1097  DO i = 1, nrep
1098  IF (i .NE. j) THEN
1099  residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), residuum(:, i))*residuum(:, i)
1100  residuum(:, j) = residuum(:, j)/sqrt(dot_product(residuum(:, j), residuum(:, j)))
1101  END IF
1102  END DO
1103  END IF
1104  END DO
1105  END DO
1106  ms_vib%b_vec = residuum
1107  DEALLOCATE (tmp)
1108  END SUBROUTINE select_vector
1109 
1110 ! **************************************************************************************************
1111 !> \brief writes the output for a mode tracking calculation
1112 !> \param iw ...
1113 !> \param converged ...
1114 !> \param freq ...
1115 !> \param criter ...
1116 !> \param ms_vib ...
1117 !> \param input ...
1118 !> \param nrep ...
1119 !> \param approx_H ...
1120 !> \param eigenval ...
1121 !> \param calc_intens ...
1122 !> \param intensities ...
1123 !> \param logger ...
1124 !> \author Florian Schiffmann 11.2007
1125 ! **************************************************************************************************
1126  SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1127  approx_H, eigenval, calc_intens, intensities, logger)
1128 
1129  INTEGER :: iw
1130  LOGICAL :: converged
1131  REAL(kind=dp), DIMENSION(:) :: freq
1132  REAL(kind=dp), DIMENSION(:, :) :: criter
1133  TYPE(ms_vib_type) :: ms_vib
1134  TYPE(section_vals_type), POINTER :: input
1135  INTEGER :: nrep
1136  REAL(kind=dp), DIMENSION(:, :) :: approx_h
1137  REAL(kind=dp), DIMENSION(:) :: eigenval
1138  LOGICAL :: calc_intens
1139  REAL(kind=dp), DIMENSION(:), OPTIONAL :: intensities
1140  TYPE(cp_logger_type), POINTER :: logger
1141 
1142  INTEGER :: i, j, msunit, stat
1143  REAL(kind=dp) :: crit_a, crit_b, fint, gintval
1144  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: residuum
1145  TYPE(section_vals_type), POINTER :: ms_vib_section
1146 
1147  ms_vib_section => section_vals_get_subs_vals(input, &
1148  "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1149 
1150  fint = 42.255_dp*massunit*debye**2*bohr**2
1151 
1152  IF (converged) THEN
1153  IF (iw .GT. 0) THEN
1154  WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1155  DO i = 1, nrep
1156  WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1157  END DO
1158  ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1159  WRITE (iw, '( /, 1X, 79("-") )')
1160  WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1161  IF (PRESENT(intensities)) THEN
1162  WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1163  ELSE
1164  WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1165  END IF
1166  DO i = 1, SIZE(ms_vib%b_mat, 2)
1167  residuum = 0._dp
1168  DO j = 1, SIZE(ms_vib%b_mat, 2)
1169  residuum(:) = residuum(:) + approx_h(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1170  END DO
1171  DO j = 1, ms_vib%mat_size
1172  residuum(:) = residuum(:) - dot_product(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1173  END DO
1174  crit_a = maxval(residuum(:))
1175  crit_b = sqrt(dot_product(residuum, residuum))
1176  IF (PRESENT(intensities)) THEN
1177  gintval = fint*intensities(i)**2
1178  IF (crit_a .LE. ms_vib%eps(1) .AND. crit_b .LE. ms_vib%eps(2)) THEN
1179  IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1180  'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1181  ELSE
1182  IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1183  'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1184  END IF
1185  ELSE
1186  IF (crit_a .LE. ms_vib%eps(1) .AND. crit_b .LE. ms_vib%eps(2)) THEN
1187  IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1188  'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1189  ELSE
1190  IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1191  'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1192  END IF
1193  END IF
1194  END DO
1195  DEALLOCATE (residuum)
1196 
1197  msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1198  "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1199  file_status="REPLACE", file_form="UNFORMATTED", &
1200  file_action="WRITE")
1201 
1202  IF (msunit > 0) THEN
1203  WRITE (unit=msunit, iostat=stat) ms_vib%mat_size
1204  WRITE (unit=msunit, iostat=stat) ms_vib%b_mat
1205  WRITE (unit=msunit, iostat=stat) ms_vib%s_mat
1206  IF (calc_intens) WRITE (unit=msunit, iostat=stat) ms_vib%dip_deriv
1207  END IF
1208 
1209  CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1210  "PRINT%MS_RESTART")
1211  END IF
1212  ELSE
1213  IF (iw .GT. 0) THEN
1214  msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1215  "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1216  file_status="REPLACE", file_form="UNFORMATTED", &
1217  file_action="WRITE")
1218 
1219  IF (msunit > 0) THEN
1220  WRITE (unit=msunit, iostat=stat) ms_vib%mat_size
1221  WRITE (unit=msunit, iostat=stat) ms_vib%b_mat
1222  WRITE (unit=msunit, iostat=stat) ms_vib%s_mat
1223  IF (calc_intens) WRITE (unit=msunit, iostat=stat) ms_vib%dip_deriv
1224  END IF
1225 
1226  CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1227  "PRINT%MS_RESTART")
1228 
1229  WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1230  DO i = 1, nrep
1231  IF (criter(1, i) .LE. 1e-7 .AND. (criter(2, i)) .LE. 1e-6) THEN
1232  WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1233  ELSE
1234  WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 NOT CONVERGED"
1235  END IF
1236  END DO
1237  END IF
1238  END IF
1239 
1240  END SUBROUTINE ms_out
1241 
1242 ! **************************************************************************************************
1243 !> \brief ...
1244 !> \param ms_vib ...
1245 !> \param approx_H ...
1246 !> \param eigenval ...
1247 !> \param residuum ...
1248 !> \param nrep ...
1249 !> \param ind ...
1250 !> \author Florian Schiffmann 11.2007
1251 ! **************************************************************************************************
1252  SUBROUTINE get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1253 
1254  TYPE(ms_vib_type) :: ms_vib
1255  REAL(kind=dp), DIMENSION(:, :) :: approx_h
1256  REAL(kind=dp), DIMENSION(:) :: eigenval
1257  REAL(kind=dp), DIMENSION(:, :) :: residuum
1258  INTEGER :: nrep
1259  INTEGER, DIMENSION(:) :: ind
1260 
1261  INTEGER :: count1, count2, i, j
1262  INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1263  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1264  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1265  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1266  REAL(kind=dp), DIMENSION(2) :: myrange
1267 
1268  myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1269  count1 = 0
1270  count2 = 0
1271  residuum = 0.0_dp
1272  ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1273  ALLOCATE (map1(SIZE(eigenval), 2))
1274  ALLOCATE (tmp(SIZE(eigenval)))
1275  DO i = 1, SIZE(eigenval)
1276  IF (abs(eigenval(i) - myrange(1)) + abs(eigenval(i) - myrange(2)) .LE. &
1277  abs(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1278  count1 = count1 + 1
1279  map1(count1, 1) = i
1280  ELSE
1281  count2 = count2 + 1
1282  map1(count2, 2) = i
1283  tmp(count2) = min(abs(eigenval(i) - myrange(1)), abs(eigenval(i) - myrange(2)))
1284  END IF
1285  END DO
1286 
1287  IF (count1 .EQ. nrep) THEN
1288  DO j = 1, count1
1289  DO i = 1, ms_vib%mat_size
1290  residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))*(ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1291  ind(j) = map1(j, 1)
1292  END DO
1293  END DO
1294  ELSE IF (count1 .GT. nrep) THEN
1295  ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1296  ALLOCATE (tmp1(count1))
1297  ALLOCATE (map2(count1))
1298  tmp_resid = 0._dp
1299  DO j = 1, count1
1300  DO i = 1, ms_vib%mat_size
1301  tmp_resid(:, j) = tmp_resid(:, j) + approx_h(i, map1(j, 1))* &
1302  (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1303  END DO
1304  END DO
1305 
1306  DO j = 1, count1
1307  DO i = 1, ms_vib%mat_size
1308  tmp_resid(:, j) = tmp_resid(:, j) - dot_product(tmp_resid(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1309  END DO
1310  tmp(j) = maxval(tmp_resid(:, j))
1311  END DO
1312  CALL sort(tmp, count1, map2)
1313  DO j = 1, nrep
1314  residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1315  ind(j) = map1(map2(count1 + 1 - j), 1)
1316  END DO
1317  DEALLOCATE (tmp_resid)
1318  DEALLOCATE (tmp1)
1319  DEALLOCATE (map2)
1320  ELSE IF (count1 .LT. nrep) THEN
1321 
1322  ALLOCATE (map2(count2))
1323  IF (count1 .NE. 0) THEN
1324  DO j = 1, count1
1325  DO i = 1, ms_vib%mat_size
1326  residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))* &
1327  (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1328  END DO
1329  ind(j) = map1(j, 1)
1330  END DO
1331  END IF
1332  CALL sort(tmp, count2, map2)
1333  DO j = 1, nrep - count1
1334  DO i = 1, ms_vib%mat_size
1335  residuum(:, count1 + j) = residuum(:, count1 + j) + approx_h(i, map1(map2(j), 2)) &
1336  *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1337  END DO
1338  ind(count1 + j) = map1(map2(j), 2)
1339  END DO
1340 
1341  DEALLOCATE (map2)
1342  END IF
1343 
1344  DEALLOCATE (map1)
1345  DEALLOCATE (tmp)
1346 
1347  END SUBROUTINE get_vibs_in_range
1348 END MODULE mode_selective
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
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
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::ms_guess_restart
integer, parameter, public ms_guess_restart
Definition: input_constants.F:832
input_constants::ms_guess_restart_vec
integer, parameter, public ms_guess_restart_vec
Definition: input_constants.F:832
input_constants::ms_guess_molden
integer, parameter, public ms_guess_molden
Definition: input_constants.F:832
input_constants::ms_guess_atomic
integer, parameter, public ms_guess_atomic
Definition: input_constants.F:832
input_constants::ms_guess_bfgs
integer, parameter, public ms_guess_bfgs
Definition: input_constants.F:832
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::max_line_length
integer, parameter, public max_line_length
Definition: kinds.F:59
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
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
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::vibfac
real(kind=dp), parameter, public vibfac
Definition: physcon.F:189
physcon::massunit
real(kind=dp), parameter, public massunit
Definition: physcon.F:141
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_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
util
All kind of helpful little routines.
Definition: util.F:14