(git:0761fd1)
Loading...
Searching...
No Matches
mode_selective.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Module performing a 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! **************************************************************************************************
20MODULE mode_selective
21 USE bfgs_optimizer, ONLY: bfgs_read_restart_header
22 USE cell_types, ONLY: cell_type
23 USE cp_files, ONLY: close_file,&
24 open_file
25 USE cp_log_handling, ONLY: cp_get_default_logger,&
26 cp_logger_get_default_io_unit,&
27 cp_logger_type
28 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 cp_print_key_unit_nr
30 USE cp_result_methods, ONLY: get_results
31 USE global_types, ONLY: global_environment_type
32 USE ieee_arithmetic, ONLY: ieee_is_finite
33 USE input_constants, ONLY: bfgs_direct,&
34 bfgs_inverse,&
35 ms_guess_atomic,&
36 ms_guess_bfgs,&
37 ms_guess_molden,&
38 ms_guess_restart,&
39 ms_guess_restart_vec
40 USE input_section_types, ONLY: section_vals_get,&
41 section_vals_get_subs_vals,&
42 section_vals_type,&
43 section_vals_val_get
44 USE iso_fortran_env, ONLY: iostat_end
45 USE kinds, ONLY: default_path_length,&
46 default_string_length,&
47 dp,&
48 max_line_length
49 USE mathlib, ONLY: diamat_all
50 USE message_passing, ONLY: mp_para_env_type
51 USE molden_utils, ONLY: write_vibrations_molden
52 USE particle_types, ONLY: particle_type
53 USE physcon, ONLY: bohr,&
54 debye,&
55 massunit,&
56 vibfac
57 USE replica_methods, ONLY: rep_env_calc_e_f
58 USE replica_types, ONLY: replica_env_type
59 USE util, ONLY: sort
60#include "./base/base_uses.f90"
61
62 IMPLICIT NONE
63
64 PRIVATE
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mode_selective'
66 LOGICAL, PARAMETER :: debug_this_module = .false.
67 REAL(KIND=dp), PARAMETER :: bfgs_guess_degeneracy_tol = 1.0e-8_dp
68
69 TYPE ms_vib_type
70 INTEGER :: mat_size = -1
71 INTEGER :: select_id = -1
72 INTEGER, DIMENSION(:), POINTER :: inv_atoms => null()
73 REAL(KIND=dp) :: eps(2) = 0.0_dp
74 REAL(KIND=dp) :: sel_freq = 0.0_dp
75 REAL(KIND=dp) :: low_freq = 0.0_dp
76 REAL(KIND=dp), POINTER, DIMENSION(:, :) :: b_vec => null()
77 REAL(KIND=dp), POINTER, DIMENSION(:, :) :: delta_vec => null()
78 REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ms_force => null()
79 REAL(KIND=dp), DIMENSION(:), POINTER :: eig_bfgs => null()
80 REAL(KIND=dp), DIMENSION(:), POINTER :: f_range => null()
81 REAL(KIND=dp), DIMENSION(:), POINTER :: inv_range => null()
82 REAL(KIND=dp), POINTER, DIMENSION(:) :: step_b => null()
83 REAL(KIND=dp), POINTER, DIMENSION(:) :: step_r => null()
84 REAL(KIND=dp), DIMENSION(:, :), POINTER :: b_mat => null()
85 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_deriv => null()
86 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hes_bfgs => null()
87 REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_mat => null()
88 INTEGER :: initial_guess = -1
89 END TYPE ms_vib_type
90
91 PUBLIC :: ms_vb_anal
92
93CONTAINS
94! **************************************************************************************************
95!> \brief Module performing a vibrational analysis
96!> \param input ...
97!> \param rep_env ...
98!> \param para_env ...
99!> \param globenv ...
100!> \param particles ...
101!> \param nrep ...
102!> \param calc_intens ...
103!> \param dx ...
104!> \param output_unit ...
105!> \param logger ...
106!> \param cell simulation cell
107!> \author Teodoro Laino 08.2006
108! **************************************************************************************************
109 SUBROUTINE ms_vb_anal(input, rep_env, para_env, globenv, particles, &
110 nrep, calc_intens, dx, output_unit, logger, cell)
111 TYPE(section_vals_type), POINTER :: input
112 TYPE(replica_env_type), POINTER :: rep_env
113 TYPE(mp_para_env_type), POINTER :: para_env
114 TYPE(global_environment_type), POINTER :: globenv
115 TYPE(particle_type), DIMENSION(:), POINTER :: particles
116 INTEGER :: nrep
117 LOGICAL :: calc_intens
118 REAL(kind=dp) :: dx
119 INTEGER :: output_unit
120 TYPE(cp_logger_type), POINTER :: logger
121 TYPE(cell_type), POINTER :: cell
122
123 CHARACTER(len=*), PARAMETER :: routinen = 'ms_vb_anal'
124
125 CHARACTER(LEN=default_string_length) :: description
126 INTEGER :: handle, i, ip1, j, natoms, ncoord
127 LOGICAL :: converged
128 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mass, pos0
129 REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp_deriv
130 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tmp_dip
131 TYPE(ms_vib_type) :: ms_vib
132
133 CALL timeset(routinen, handle)
134 converged = .false.
135 natoms = SIZE(particles)
136 ncoord = 3*natoms
137 ALLOCATE (mass(3*natoms))
138 DO i = 1, natoms
139 DO j = 1, 3
140 mass((i - 1)*3 + j) = particles(i)%atomic_kind%mass
141 mass((i - 1)*3 + j) = sqrt(mass((i - 1)*3 + j))
142 END DO
143 END DO
144 ! Allocate working arrays
145 ALLOCATE (ms_vib%delta_vec(ncoord, nrep))
146 ALLOCATE (ms_vib%b_vec(ncoord, nrep))
147 ALLOCATE (ms_vib%step_r(nrep))
148 ALLOCATE (ms_vib%step_b(nrep))
149 IF (calc_intens) THEN
150 description = '[DIPOLE]'
151 ALLOCATE (tmp_dip(nrep, 3, 2))
152 ALLOCATE (ms_vib%dip_deriv(3, nrep))
153 END IF
154 CALL ms_initial_moves(para_env, nrep, input, globenv, ms_vib, &
155 particles, &
156 mass, &
157 dx, &
158 calc_intens, logger)
159 ncoord = 3*natoms
160 ALLOCATE (pos0(ncoord))
161 ALLOCATE (ms_vib%ms_force(ncoord, nrep))
162 DO i = 1, natoms
163 DO j = 1, 3
164 pos0((i - 1)*3 + j) = particles((i))%r(j)
165 END DO
166 END DO
167 ncoord = 3*natoms
168 DO
169 ms_vib%ms_force = huge(0.0_dp)
170 DO i = 1, nrep
171 DO j = 1, ncoord
172 rep_env%r(j, i) = pos0(j) + ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
173 END DO
174 END DO
175 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
176
177 DO i = 1, nrep
178 IF (calc_intens) THEN
179 CALL get_results(results=rep_env%results(i)%results, &
180 description=description, &
181 n_rep=ip1)
182 CALL get_results(results=rep_env%results(i)%results, &
183 description=description, &
184 values=tmp_dip(i, :, 1), &
185 nval=ip1)
186 END IF
187 DO j = 1, ncoord
188 ms_vib%ms_force(j, i) = rep_env%f(j, i)
189 END DO
190 END DO
191 DO i = 1, nrep
192 DO j = 1, ncoord
193 rep_env%r(j, i) = pos0(j) - ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
194 END DO
195 END DO
196 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
197 IF (calc_intens) THEN
198 DO i = 1, nrep
199 CALL get_results(results=rep_env%results(i)%results, &
200 description=description, &
201 n_rep=ip1)
202 CALL get_results(results=rep_env%results(i)%results, &
203 description=description, &
204 values=tmp_dip(i, :, 2), &
205 nval=ip1)
206 ms_vib%dip_deriv(:, ms_vib%mat_size + i) = (tmp_dip(i, :, 1) - tmp_dip(i, :, 2))/(2*ms_vib%step_b(i))
207 END DO
208 END IF
209
210 CALL evaluate_h_update_b(rep_env, ms_vib, input, nrep, &
211 particles, &
212 mass, &
213 converged, &
214 dx, calc_intens, &
215 output_unit, logger, cell)
216 IF (converged) EXIT
217 IF (calc_intens) THEN
218 ALLOCATE (tmp_deriv(3, ms_vib%mat_size))
219 tmp_deriv = ms_vib%dip_deriv
220 DEALLOCATE (ms_vib%dip_deriv)
221 ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
222 ms_vib%dip_deriv(:, 1:ms_vib%mat_size) = tmp_deriv(:, 1:ms_vib%mat_size)
223 DEALLOCATE (tmp_deriv)
224 END IF
225 END DO
226 DEALLOCATE (ms_vib%ms_force)
227 DEALLOCATE (pos0)
228 DEALLOCATE (ms_vib%step_r)
229 DEALLOCATE (ms_vib%step_b)
230 DEALLOCATE (ms_vib%b_vec)
231 DEALLOCATE (ms_vib%delta_vec)
232 DEALLOCATE (mass)
233 DEALLOCATE (ms_vib%b_mat)
234 DEALLOCATE (ms_vib%s_mat)
235 IF (ms_vib%select_id == 3) THEN
236 DEALLOCATE (ms_vib%inv_atoms)
237 END IF
238 IF (ASSOCIATED(ms_vib%eig_bfgs)) THEN
239 DEALLOCATE (ms_vib%eig_bfgs)
240 END IF
241 IF (ASSOCIATED(ms_vib%hes_bfgs)) THEN
242 DEALLOCATE (ms_vib%hes_bfgs)
243 END IF
244 IF (calc_intens) THEN
245 DEALLOCATE (ms_vib%dip_deriv)
246 DEALLOCATE (tmp_dip)
247 END IF
248 CALL timestop(handle)
249 END SUBROUTINE ms_vb_anal
250! **************************************************************************************************
251!> \brief Generates the first displacement vector for a mode selctive vibrational
252!> analysis. At the moment this is a random number for selected atoms
253!> \param para_env ...
254!> \param nrep ...
255!> \param input ...
256!> \param globenv ...
257!> \param ms_vib ...
258!> \param particles ...
259!> \param mass ...
260!> \param dx ...
261!> \param calc_intens ...
262!> \param logger ...
263!> \author Florian Schiffmann 11.2007
264! **************************************************************************************************
265 SUBROUTINE ms_initial_moves(para_env, nrep, input, globenv, ms_vib, particles, &
266 mass, dx, &
267 calc_intens, logger)
268 TYPE(mp_para_env_type), POINTER :: para_env
269 INTEGER :: nrep
270 TYPE(section_vals_type), POINTER :: input
271 TYPE(global_environment_type), POINTER :: globenv
272 TYPE(ms_vib_type) :: ms_vib
273 TYPE(particle_type), DIMENSION(:), POINTER :: particles
274 REAL(kind=dp), DIMENSION(:) :: mass
275 REAL(kind=dp) :: dx
276 LOGICAL :: calc_intens
277 TYPE(cp_logger_type), POINTER :: logger
278
279 CHARACTER(len=*), PARAMETER :: routinen = 'MS_initial_moves'
280
281 INTEGER :: guess, handle, i, j, jj, k, m, &
282 n_rep_val, natoms, ncoord
283 INTEGER, ALLOCATABLE, DIMENSION(:) :: map_atoms
284 INTEGER, DIMENSION(:), POINTER :: tmplist
285 LOGICAL :: do_involved_atoms, ionode
286 REAL(kind=dp) :: my_val, norm
287 TYPE(section_vals_type), POINTER :: involved_at_section, ms_vib_section
288
289 CALL timeset(routinen, handle)
290 NULLIFY (ms_vib%eig_bfgs, ms_vib%f_range, ms_vib%hes_bfgs, ms_vib%inv_range)
291 ms_vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
292 CALL section_vals_val_get(ms_vib_section, "INITIAL_GUESS", i_val=guess)
293 CALL section_vals_val_get(ms_vib_section, "EPS_MAX_VAL", r_val=ms_vib%eps(1))
294 CALL section_vals_val_get(ms_vib_section, "EPS_NORM", r_val=ms_vib%eps(2))
295 CALL section_vals_val_get(ms_vib_section, "RANGE", n_rep_val=n_rep_val)
296 ms_vib%select_id = 0
297 IF (n_rep_val /= 0) THEN
298 CALL section_vals_val_get(ms_vib_section, "RANGE", r_vals=ms_vib%f_range)
299 IF (ms_vib%f_range(1) > ms_vib%f_range(2)) THEN
300 my_val = ms_vib%f_range(2)
301 ms_vib%f_range(2) = ms_vib%f_range(1)
302 ms_vib%f_range(1) = my_val
303 END IF
304 ms_vib%select_id = 2
305 END IF
306 CALL section_vals_val_get(ms_vib_section, "FREQUENCY", r_val=ms_vib%sel_freq)
307 CALL section_vals_val_get(ms_vib_section, "LOWEST_FREQUENCY", r_val=ms_vib%low_freq)
308 IF (ms_vib%sel_freq > 0._dp) ms_vib%select_id = 1
309 involved_at_section => section_vals_get_subs_vals(ms_vib_section, "INVOLVED_ATOMS")
310 CALL section_vals_get(involved_at_section, explicit=do_involved_atoms)
311 IF (do_involved_atoms) THEN
312 CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", n_rep_val=n_rep_val)
313 jj = 0
314 DO k = 1, n_rep_val
315 CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=k, i_vals=tmplist)
316 DO j = 1, SIZE(tmplist)
317 jj = jj + 1
318 END DO
319 END DO
320 IF (jj >= 1) THEN
321 natoms = jj
322 ALLOCATE (ms_vib%inv_atoms(natoms))
323 jj = 0
324 DO m = 1, n_rep_val
325 CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=m, i_vals=tmplist)
326 DO j = 1, SIZE(tmplist)
327 ms_vib%inv_atoms(j) = tmplist(j)
328 END DO
329 END DO
330 ms_vib%select_id = 3
331 END IF
332 CALL section_vals_val_get(involved_at_section, "RANGE", n_rep_val=n_rep_val)
333 IF (n_rep_val /= 0) THEN
334 CALL section_vals_val_get(involved_at_section, "RANGE", r_vals=ms_vib%inv_range)
335 IF (ms_vib%inv_range(1) > ms_vib%inv_range(2)) THEN
336 ms_vib%inv_range(2) = my_val
337 ms_vib%inv_range(2) = ms_vib%inv_range(1)
338 ms_vib%inv_range(1) = my_val
339 END IF
340 END IF
341 END IF
342 IF (ms_vib%select_id == 0) THEN
343 cpabort("no frequency, range or involved atoms specified ")
344 END IF
345 ionode = para_env%is_source()
346 SELECT CASE (guess)
347 CASE (ms_guess_atomic)
348 ms_vib%initial_guess = 1
349 CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
350 jj = 0
351 DO k = 1, n_rep_val
352 CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
353 DO j = 1, SIZE(tmplist)
354 jj = jj + 1
355 END DO
356 END DO
357 IF (jj < 1) THEN
358 natoms = SIZE(particles)
359 ALLOCATE (map_atoms(natoms))
360 DO j = 1, natoms
361 map_atoms(j) = j
362 END DO
363 ELSE
364 natoms = jj
365 ALLOCATE (map_atoms(natoms))
366 jj = 0
367 DO m = 1, n_rep_val
368 CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
369 DO j = 1, SIZE(tmplist)
370 map_atoms(j) = tmplist(j)
371 END DO
372 END DO
373 END IF
374
375 ! apply random displacement along the mass weighted nuclear cartesian coordinates
376 ms_vib%b_vec = 0._dp
377 ms_vib%delta_vec = 0._dp
378 jj = 0
379
380 DO i = 1, nrep
381 DO j = 1, natoms
382 DO k = 1, 3
383 jj = (map_atoms(j) - 1)*3 + k
384 ms_vib%b_vec(jj, i) = abs(globenv%gaussian_rng_stream%next())
385 END DO
386 END DO
387 norm = norm2(ms_vib%b_vec(:, i))
388 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
389 END DO
390
391 IF (nrep > 1) THEN
392 DO k = 1, 10
393 DO j = 1, nrep
394 DO i = 1, nrep
395 IF (i /= j) THEN
396 ms_vib%b_vec(:, j) = &
397 ms_vib%b_vec(:, j) - dot_product(ms_vib%b_vec(:, j), ms_vib%b_vec(:, i))*ms_vib%b_vec(:, i)
398 ms_vib%b_vec(:, j) = &
399 ms_vib%b_vec(:, j)/norm2(ms_vib%b_vec(:, j))
400 END IF
401 END DO
402 END DO
403 END DO
404 END IF
405
406 ms_vib%mat_size = 0
407 DO i = 1, SIZE(ms_vib%b_vec, 1)
408 ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
409 END DO
410 CASE (ms_guess_bfgs)
411
412 ms_vib%initial_guess = 2
413 CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
414 ms_vib%mat_size = 0
415
416 CASE (ms_guess_restart_vec)
417
418 ms_vib%initial_guess = 3
419 ncoord = 3*SIZE(particles)
420 CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
421
422 ms_vib%mat_size = 0
423 CASE (ms_guess_restart)
424 ms_vib%initial_guess = 4
425 ncoord = 3*SIZE(particles)
426 CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
427
428 CASE (ms_guess_molden)
429 ms_vib%initial_guess = 5
430 ncoord = 3*SIZE(particles)
431 CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
432 ms_vib%mat_size = 0
433 END SELECT
434 CALL para_env%bcast(ms_vib%b_vec)
435 CALL para_env%bcast(ms_vib%delta_vec)
436 DO i = 1, nrep
437 ms_vib%step_r(i) = dx/norm2(ms_vib%delta_vec(:, i))
438 ms_vib%step_b(i) = norm2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
439 END DO
440 CALL timestop(handle)
441
442 END SUBROUTINE ms_initial_moves
443
444! **************************************************************************************************
445!> \brief Read a legacy or tagged BFGS restart as a direct ionic Hessian.
446!> \param unit File unit
447!> \param ncoord Number of ionic coordinates needed by mode-selective analysis
448!> \param hessian Direct ionic Hessian
449! **************************************************************************************************
450 SUBROUTINE read_bfgs_guess_hessian(unit, ncoord, hessian)
451 INTEGER, INTENT(IN) :: unit, ncoord
452 REAL(kind=dp), DIMENSION(ncoord, ncoord), &
453 INTENT(OUT) :: hessian
454
455 INTEGER :: i, ios, stored_method, stored_ndf, &
456 stored_type
457 LOGICAL :: tagged
458 REAL(kind=dp) :: extra
459 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: column, eigenvalues
460 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: stored_matrix
461 REAL(kind=dp), DIMENSION(3, 3) :: cell_frame, reference_cell
462
463 CALL bfgs_read_restart_header(unit, stored_method, stored_ndf, stored_type, tagged, &
464 reference_cell, cell_frame)
465
466 IF (.NOT. tagged) THEN
467 ! Legacy restart: preserve the historical behaviour. Reading ncoord values from
468 ! the first ncoord records also extracts the ionic block from a CELL_OPT Hessian.
469 DO i = 1, ncoord
470 READ (unit, iostat=ios) hessian(:, i)
471 IF (ios /= 0) cpabort("Error while reading legacy BFGS Hessian")
472 END DO
473 RETURN
474 END IF
475
476 IF (stored_ndf < ncoord) THEN
477 cpabort("BFGS restart matrix is smaller than the mode-selective ionic Hessian")
478 END IF
479
480 SELECT CASE (stored_method)
481 CASE (bfgs_direct)
482 ALLOCATE (column(stored_ndf))
483 DO i = 1, stored_ndf
484 READ (unit, iostat=ios) column(:)
485 IF (ios /= 0) cpabort("Truncated BFGS restart matrix")
486 IF (.NOT. all(ieee_is_finite(column(:)))) cpabort("Non-finite BFGS restart matrix")
487 IF (i <= ncoord) hessian(:, i) = column(1:ncoord)
488 END DO
489 DEALLOCATE (column)
490
491 CASE (bfgs_inverse)
492 ALLOCATE (stored_matrix(stored_ndf, stored_ndf), eigenvalues(stored_ndf))
493 DO i = 1, stored_ndf
494 READ (unit, iostat=ios) stored_matrix(:, i)
495 IF (ios /= 0) cpabort("Truncated inverse BFGS restart matrix")
496 END DO
497 IF (.NOT. all(ieee_is_finite(stored_matrix(:, :)))) THEN
498 cpabort("Non-finite inverse BFGS restart matrix")
499 END IF
500
501 ! Invert the complete matrix before taking the ionic block. For CELL_OPT,
502 ! inverse(H_RR) is not in general equal to the RR block of inverse(H).
503 CALL diamat_all(stored_matrix, eigenvalues)
504 IF (.NOT. all(ieee_is_finite(eigenvalues(:)))) THEN
505 cpabort("Non-finite inverse BFGS restart spectrum")
506 END IF
507 IF (minval(eigenvalues(:)) <= 0.0_dp) THEN
508 cpabort("Inverse BFGS restart matrix is not positive definite")
509 END IF
510 eigenvalues(:) = 1.0_dp/eigenvalues(:)
511 IF (.NOT. all(ieee_is_finite(eigenvalues(:)))) THEN
512 cpabort("Unrepresentable inverse BFGS restart matrix")
513 END IF
514 DO i = 1, stored_ndf
515 stored_matrix(1:ncoord, i) = &
516 sqrt(eigenvalues(i))*stored_matrix(1:ncoord, i)
517 END DO
518 hessian(:, :) = matmul(stored_matrix(1:ncoord, :), &
519 transpose(stored_matrix(1:ncoord, :)))
520 DEALLOCATE (stored_matrix, eigenvalues)
521
522 CASE DEFAULT
523 cpabort("Unknown matrix representation in BFGS restart")
524 END SELECT
525
526 ! Tagged files have a known matrix dimension, so reject trailing/corrupt data.
527 READ (unit, iostat=ios) extra
528 IF (ios /= iostat_end) cpabort("Extra data in BFGS restart matrix")
529 mark_used(stored_type)
530 mark_used(reference_cell)
531 mark_used(cell_frame)
532 END SUBROUTINE read_bfgs_guess_hessian
533
534! **************************************************************************************************
535!> \brief ...
536!> \param ms_vib_section ...
537!> \param ms_vib ...
538!> \param particles ...
539!> \param mass ...
540!> \param para_env ...
541!> \param nrep ...
542!> \author Florian Schiffmann 11.2007
543! **************************************************************************************************
544 SUBROUTINE bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
545
546 TYPE(section_vals_type), POINTER :: ms_vib_section
547 TYPE(ms_vib_type) :: ms_vib
548 TYPE(particle_type), DIMENSION(:), POINTER :: particles
549 REAL(kind=dp), DIMENSION(:) :: mass
550 TYPE(mp_para_env_type), POINTER :: para_env
551 INTEGER :: nrep
552
553 CHARACTER(LEN=default_path_length) :: hes_filename
554 INTEGER :: hesunit, i, j, jj, k, natoms, ncoord, &
555 output_unit
556 INTEGER, DIMENSION(:), POINTER :: tmplist
557 REAL(kind=dp) :: my_val, norm
558 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
559 TYPE(cp_logger_type), POINTER :: logger
560
561 logger => cp_get_default_logger()
562 output_unit = cp_logger_get_default_io_unit(logger)
563
564 natoms = SIZE(particles)
565 ncoord = 3*natoms
566
567 ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
568 ALLOCATE (ms_vib%eig_bfgs(ncoord))
569
570 IF (para_env%is_source()) THEN
571 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
572 IF (hes_filename == "") hes_filename = "HESSIAN"
573 CALL open_file(file_name=hes_filename, file_status="OLD", &
574 file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
575 ALLOCATE (tmp(ncoord))
576 ALLOCATE (tmplist(ncoord))
577
578 CALL read_bfgs_guess_hessian(hesunit, ncoord, ms_vib%hes_bfgs)
579 CALL close_file(hesunit)
580 IF (output_unit > 0) THEN
581 WRITE (output_unit, fmt="(/,T2,A)") &
582 "*** Initial Hessian has been read successfully ***"
583 END IF
584 DO i = 1, ncoord
585 DO j = 1, ncoord
586 ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
587 END DO
588 END DO
589
590 CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
591 tmp(:) = 0._dp
592 IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
593 IF (ms_vib%select_id == 2) my_val = (((ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp)/vibfac)**2/massunit
594 IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
595 DO i = 1, ncoord
596 tmp(i) = abs(my_val - ms_vib%eig_bfgs(i))
597 END DO
598 ELSE IF (ms_vib%select_id == 3) THEN
599 CALL bfgs_canonicalize_degenerate_modes(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
600 DO i = 1, ncoord
601 DO j = 1, SIZE(ms_vib%inv_atoms)
602 DO k = 1, 3
603 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
604 tmp(i) = tmp(i) + sqrt(ms_vib%hes_bfgs(jj, i)**2)
605 END DO
606 END DO
607 IF ((sign(1._dp, ms_vib%eig_bfgs(i))*sqrt(abs(ms_vib%eig_bfgs(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
608 END DO
609 tmp(:) = -tmp(:)
610 END IF
611 CALL sort(tmp, ncoord, tmplist)
612 DO i = 1, nrep
613 ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
614 norm = norm2(ms_vib%b_vec(:, i))
615 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
616 END DO
617 DO i = 1, SIZE(ms_vib%b_vec, 1)
618 ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
619 END DO
620 DEALLOCATE (tmp)
621 DEALLOCATE (tmplist)
622 END IF
623
624 CALL para_env%bcast(ms_vib%b_vec)
625 CALL para_env%bcast(ms_vib%delta_vec)
626
627 DEALLOCATE (ms_vib%hes_bfgs)
628 DEALLOCATE (ms_vib%eig_bfgs)
629 ms_vib%mat_size = 0
630
631 END SUBROUTINE bfgs_guess
632
633! **************************************************************************************************
634!> \brief Replace arbitrary bases of near-degenerate BFGS eigenspaces by deterministic ones.
635!> This leaves the existing involved-atom scoring and mode ordering unchanged.
636!> \param modes Orthonormal Hessian eigenvectors
637!> \param eigenvalues Ascending Hessian eigenvalues
638! **************************************************************************************************
639 SUBROUTINE bfgs_canonicalize_degenerate_modes(modes, eigenvalues)
640 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: modes
641 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
642
643 INTEGER :: first, last, n
644 REAL(kind=dp) :: cutoff, scale
645
646 n = SIZE(eigenvalues)
647 cpassert(SIZE(modes, 1) == n .AND. SIZE(modes, 2) == n)
648 cutoff = (400.0_dp/vibfac)**2/massunit
649
650 first = 1
651 DO WHILE (first <= n)
652 IF (eigenvalues(first) <= cutoff) THEN
653 first = first + 1
654 cycle
655 END IF
656 last = first
657 DO WHILE (last < n)
658 scale = max(abs(eigenvalues(first)), abs(eigenvalues(last + 1)))
659 IF (eigenvalues(last + 1) - eigenvalues(first) > bfgs_guess_degeneracy_tol*scale) EXIT
660 last = last + 1
661 END DO
662 IF (last > first) CALL bfgs_canonical_modes(modes(:, first:last))
663 first = last + 1
664 END DO
665 END SUBROUTINE bfgs_canonicalize_degenerate_modes
666
667! **************************************************************************************************
668!> \brief Fix a subspace basis using projections of Cartesian coordinate axes.
669!> \param modes Orthonormal basis, replaced by a deterministic orthonormal basis
670! **************************************************************************************************
671 SUBROUTINE bfgs_canonical_modes(modes)
672 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: modes
673
674 INTEGER :: i, j, k, m, n, pass, pivot
675 REAL(kind=dp) :: maximum, norm
676 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: norms, projection, vector
677 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: residual, rotation
678
679 n = SIZE(modes, 1)
680 m = SIZE(modes, 2)
681 ALLOCATE (residual(n, m), rotation(m, m), norms(n), projection(n), vector(m))
682 residual(:, :) = modes(:, :)
683 rotation(:, :) = 0.0_dp
684 DO j = 1, m
685 norms(:) = sum(residual(:, :)**2, dim=2)
686 maximum = maxval(norms)
687 pivot = 1
688 DO i = 1, n
689 IF (norms(i) >= maximum*(1.0_dp - bfgs_guess_degeneracy_tol)) THEN
690 pivot = i
691 EXIT
692 END IF
693 END DO
694 vector(:) = residual(pivot, :)
695 DO pass = 1, 2
696 DO k = 1, j - 1
697 vector(:) = vector(:) - dot_product(rotation(:, k), vector)*rotation(:, k)
698 END DO
699 END DO
700 norm = norm2(vector)
701 cpassert(norm > sqrt(epsilon(1.0_dp)))
702 rotation(:, j) = vector(:)/norm
703 projection(:) = matmul(residual, rotation(:, j))
704 DO k = 1, m
705 residual(:, k) = residual(:, k) - projection(:)*rotation(k, j)
706 END DO
707 END DO
708 residual(:, :) = matmul(modes, rotation)
709 modes(:, :) = residual(:, :)
710 DEALLOCATE (residual, rotation, norms, projection, vector)
711 END SUBROUTINE bfgs_canonical_modes
712
713! **************************************************************************************************
714!> \brief ...
715!> \param ms_vib_section ...
716!> \param para_env ...
717!> \param ms_vib ...
718!> \param mass ...
719!> \param ionode ...
720!> \param particles ...
721!> \param nrep ...
722!> \param calc_intens ...
723!> \author Florian Schiffmann 11.2007
724! **************************************************************************************************
725 SUBROUTINE rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
726
727 TYPE(section_vals_type), POINTER :: ms_vib_section
728 TYPE(mp_para_env_type), POINTER :: para_env
729 TYPE(ms_vib_type) :: ms_vib
730 REAL(kind=dp), DIMENSION(:) :: mass
731 LOGICAL :: ionode
732 TYPE(particle_type), DIMENSION(:), POINTER :: particles
733 INTEGER :: nrep
734 LOGICAL :: calc_intens
735
736 CHARACTER(LEN=default_path_length) :: ms_filename
737 INTEGER :: hesunit, i, j, mat, natoms, ncoord, &
738 output_unit, stat, statint
739 INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
740 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
741 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h
742 TYPE(cp_logger_type), POINTER :: logger
743
744 logger => cp_get_default_logger()
745 output_unit = cp_logger_get_default_io_unit(logger)
746
747 natoms = SIZE(particles)
748 ncoord = 3*natoms
749 IF (calc_intens) THEN
750 DEALLOCATE (ms_vib%dip_deriv)
751 END IF
752
753 IF (ionode) THEN
754
755 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
756 IF (ms_filename == "") ms_filename = "MS_RESTART"
757 CALL open_file(file_name=ms_filename, &
758 file_status="UNKNOWN", &
759 file_form="UNFORMATTED", &
760 file_action="READ", &
761 unit_number=hesunit)
762 READ (unit=hesunit, iostat=stat) mat
763 cpassert(stat == 0)
764 ms_vib%mat_size = mat
765 END IF
766 CALL para_env%bcast(ms_vib%mat_size)
767 ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
768 ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
769 IF (calc_intens) THEN
770 ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
771 END IF
772 IF (ionode) THEN
773 statint = 0
774 READ (unit=hesunit) ms_vib%b_mat
775 READ (unit=hesunit, iostat=stat) ms_vib%s_mat
776 IF (stat /= 0 .AND. output_unit > 0) THEN
777 WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART **"
778 END IF
779 IF (calc_intens) THEN
780 READ (unit=hesunit, iostat=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
781 IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART,", &
782 "intensities are requested but not present in restart file **"
783 END IF
784 CALL close_file(hesunit)
785 IF (stat == 0 .AND. statint == 0 .AND. output_unit > 0) THEN
786 WRITE (output_unit, fmt="(/,T2,A)") "*** MS_RESTART has been read successfully ***"
787 END IF
788 END IF
789 CALL para_env%bcast(ms_vib%b_mat)
790 CALL para_env%bcast(ms_vib%s_mat)
791 IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
792 ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
793 ALLOCATE (eigenval(ms_vib%mat_size))
794 ALLOCATE (ind(ms_vib%mat_size))
795
796 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), &
797 ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
798 CALL diamat_all(approx_h, eigenval)
799
800 CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, ms_vib%b_vec)
801 IF (ms_vib%initial_guess /= 4) THEN
802
803 ms_vib%b_vec = 0._dp
804 DO i = 1, nrep
805 DO j = 1, ms_vib%mat_size
806 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_h(j, ind(i))*ms_vib%b_mat(:, j)
807 END DO
808 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm2(ms_vib%b_vec(:, i))
809 END DO
810
811 DEALLOCATE (ms_vib%s_mat)
812 DEALLOCATE (ms_vib%b_mat)
813 IF (calc_intens) THEN
814 DEALLOCATE (ms_vib%dip_deriv)
815 ALLOCATE (ms_vib%dip_deriv(3, nrep))
816 END IF
817 END IF
818 DEALLOCATE (approx_h)
819 DEALLOCATE (eigenval)
820 DEALLOCATE (ind)
821 DO i = 1, nrep
822 ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
823 END DO
824
825 END SUBROUTINE rest_guess
826
827! **************************************************************************************************
828!> \brief ...
829!> \param ms_vib_section ...
830!> \param input ...
831!> \param para_env ...
832!> \param ms_vib ...
833!> \param mass ...
834!> \param ncoord ...
835!> \param nrep ...
836!> \param logger ...
837!> \author Florian Schiffmann 11.2007
838! **************************************************************************************************
839 SUBROUTINE molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
840 TYPE(section_vals_type), POINTER :: ms_vib_section, input
841 TYPE(mp_para_env_type), POINTER :: para_env
842 TYPE(ms_vib_type) :: ms_vib
843 REAL(kind=dp), DIMENSION(:) :: mass
844 INTEGER :: ncoord, nrep
845 TYPE(cp_logger_type), POINTER :: logger
846
847 CHARACTER(LEN=2) :: at_name
848 CHARACTER(LEN=default_path_length) :: ms_filename
849 CHARACTER(LEN=max_line_length) :: info
850 INTEGER :: i, istat, iw, j, jj, k, nvibs, &
851 output_molden, output_unit, stat
852 INTEGER, DIMENSION(:), POINTER :: tmplist
853 LOGICAL :: reading_vib
854 REAL(kind=dp) :: my_val, norm
855 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
856 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
857 REAL(kind=dp), DIMENSION(3, ncoord/3) :: pos
858
859 output_unit = cp_logger_get_default_io_unit(logger)
860
861 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
862 IF (ms_filename == "") output_molden = &
863 cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
864 extension=".mol", file_status='UNKNOWN', &
865 file_action="READ")
866 IF (para_env%is_source()) THEN
867
868 IF (ms_filename == "") THEN
869 iw = output_molden
870 ELSE
871 CALL open_file(file_name=trim(ms_filename), &
872 file_status="UNKNOWN", &
873 file_form="FORMATTED", &
874 file_action="READ", &
875 unit_number=iw)
876 END IF
877 info = ""
878 READ (iw, *) info
879 READ (iw, *) info
880 istat = 0
881 nvibs = 0
882 reading_vib = .false.
883 DO
884 READ (iw, *, iostat=stat) info
885 istat = istat + stat
886 IF (trim(adjustl(info)) == "[FR-COORD]") EXIT
887
888 cpassert(stat == 0)
889
890 IF (reading_vib) nvibs = nvibs + 1
891 IF (trim(adjustl(info)) == "[FREQ]") reading_vib = .true.
892 END DO
893 rewind(iw)
894 istat = 0
895 READ (iw, *, iostat=stat) info
896 istat = istat + stat
897 READ (iw, *, iostat=stat) info
898 istat = istat + stat
899 ! Skip [Atoms] section
900 DO
901 READ (iw, *, iostat=stat) info
902 istat = istat + stat
903 cpassert(stat == 0)
904 IF (trim(adjustl(info)) == "[FREQ]") EXIT
905 END DO
906 ! Read frequencies and modes
907 ALLOCATE (freq(nvibs))
908 ALLOCATE (modes(ncoord, nvibs))
909
910 DO i = 1, nvibs
911 READ (iw, *, iostat=stat) freq(i)
912 istat = istat + stat
913 END DO
914 READ (iw, *) info
915 DO i = 1, ncoord/3
916 READ (iw, *, iostat=stat) at_name, pos(:, i)
917 istat = istat + stat
918 END DO
919 READ (iw, *) info
920 DO i = 1, nvibs
921 READ (iw, *) info
922 istat = istat + stat
923 DO j = 1, ncoord/3
924 k = (j - 1)*3 + 1
925 READ (iw, *, iostat=stat) modes(k:k + 2, i)
926 istat = istat + stat
927 END DO
928 END DO
929 IF (ms_filename /= "") CALL close_file(iw)
930 IF (output_unit > 0) THEN
931 IF (istat /= 0) THEN
932 WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MOLDEN file **"
933 ELSE
934 WRITE (output_unit, fmt="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
935 END IF
936 END IF
937 !!!!!!! select modes !!!!!!
938 ALLOCATE (tmp(nvibs))
939 tmp(:) = 0.0_dp
940 ALLOCATE (tmplist(nvibs))
941 IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
942 IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
943 IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
944 DO i = 1, nvibs
945 tmp(i) = abs(my_val - freq(i))
946 END DO
947 ELSE IF (ms_vib%select_id == 3) THEN
948 DO i = 1, nvibs
949 DO j = 1, SIZE(ms_vib%inv_atoms)
950 DO k = 1, 3
951 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
952 tmp(i) = tmp(i) + sqrt(modes(jj, i)**2)
953 END DO
954 END DO
955 IF (freq(i) <= 400._dp) tmp(i) = 0._dp
956 END DO
957 tmp(:) = -tmp(:)
958 END IF
959 CALL sort(tmp, nvibs, tmplist)
960 DO i = 1, nrep
961 ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
962 norm = norm2(ms_vib%b_vec(:, i))
963 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
964 END DO
965 DO i = 1, nrep
966 ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
967 END DO
968
969 DEALLOCATE (freq)
970 DEALLOCATE (modes)
971 DEALLOCATE (tmp)
972 DEALLOCATE (tmplist)
973
974 END IF
975 CALL para_env%bcast(ms_vib%b_vec)
976 CALL para_env%bcast(ms_vib%delta_vec)
977
978 IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
979 "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
980 END SUBROUTINE molden_guess
981
982! **************************************************************************************************
983!> \brief Davidson algorithm for to generate a approximate Hessian for mode
984!> selective vibrational analysis
985!> \param rep_env ...
986!> \param ms_vib ...
987!> \param input ...
988!> \param nrep ...
989!> \param particles ...
990!> \param mass ...
991!> \param converged ...
992!> \param dx ...
993!> \param calc_intens ...
994!> \param output_unit_ms ...
995!> \param logger ...
996!> \param cell simulation cell
997!> \author Florian Schiffmann 11.2007
998! **************************************************************************************************
999 SUBROUTINE evaluate_h_update_b(rep_env, ms_vib, input, nrep, &
1000 particles, &
1001 mass, &
1002 converged, dx, &
1003 calc_intens, output_unit_ms, logger, cell)
1004 TYPE(replica_env_type), POINTER :: rep_env
1005 TYPE(ms_vib_type) :: ms_vib
1006 TYPE(section_vals_type), POINTER :: input
1007 INTEGER :: nrep
1008 TYPE(particle_type), DIMENSION(:), POINTER :: particles
1009 REAL(kind=dp), DIMENSION(:) :: mass
1010 LOGICAL :: converged
1011 REAL(kind=dp) :: dx
1012 LOGICAL :: calc_intens
1013 INTEGER :: output_unit_ms
1014 TYPE(cp_logger_type), POINTER :: logger
1015 TYPE(cell_type), POINTER :: cell
1016
1017 INTEGER :: i, j, jj, k, natoms, ncoord
1018 INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
1019 LOGICAL :: dump_only_positive
1020 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
1021 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h, h_save, residuum, tmp_b, tmp_s
1022 REAL(kind=dp), DIMENSION(2, nrep) :: criteria
1023 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
1024
1025 natoms = SIZE(particles)
1026 ncoord = 3*natoms
1027 nrep = SIZE(rep_env%f, 2)
1028
1029 !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
1030 IF (ms_vib%mat_size /= 0) THEN
1031
1032 ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
1033 ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
1034
1035 tmp_b(:, :) = ms_vib%b_mat
1036 tmp_s(:, :) = ms_vib%s_mat
1037
1038 DEALLOCATE (ms_vib%b_mat)
1039 DEALLOCATE (ms_vib%s_mat)
1040 END IF
1041
1042 ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
1043 ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
1044
1045 ms_vib%s_mat = 0.0_dp
1046
1047 DO i = 1, 3*natoms
1048 IF (ms_vib%mat_size /= 0) THEN
1049 DO j = 1, ms_vib%mat_size
1050 ms_vib%b_mat(i, j) = tmp_b(i, j)
1051 ms_vib%s_mat(i, j) = tmp_s(i, j)
1052 END DO
1053 END IF
1054 DO j = 1, nrep
1055 ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
1056 END DO
1057 END DO
1058
1059 IF (ms_vib%mat_size /= 0) THEN
1060 DEALLOCATE (tmp_s)
1061 DEALLOCATE (tmp_b)
1062 END IF
1063
1064 ms_vib%mat_size = ms_vib%mat_size + nrep
1065
1066 ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
1067 ALLOCATE (h_save(ms_vib%mat_size, ms_vib%mat_size))
1068 ALLOCATE (eigenval(ms_vib%mat_size))
1069
1070 !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
1071
1072 DO i = 1, nrep
1073 DO j = 1, 3*natoms
1074 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))
1075 END DO
1076 END DO
1077
1078 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), &
1079 ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
1080 h_save(:, :) = approx_h
1081
1082 CALL diamat_all(approx_h, eigenval)
1083
1084 !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
1085 ALLOCATE (ind(ms_vib%mat_size))
1086 ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
1087
1088 CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, residuum, criteria)
1089
1090 DO i = 1, nrep
1091 DO j = 1, natoms
1092 DO k = 1, 3
1093 jj = (j - 1)*3 + k
1094 ms_vib%delta_vec(jj, i) = ms_vib%b_vec(jj, i)/mass(jj)
1095 END DO
1096 END DO
1097 END DO
1098
1099 DO i = 1, nrep
1100 ms_vib%step_r(i) = dx/norm2(ms_vib%delta_vec(:, i))
1101 ms_vib%step_b(i) = norm2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
1102 END DO
1103 converged = .false.
1104 IF (maxval(criteria(1, :)) <= ms_vib%eps(1) .AND. maxval(criteria(2, :)) &
1105 <= ms_vib%eps(2) .OR. ms_vib%mat_size >= ncoord) converged = .true.
1106 ALLOCATE (freq(nrep))
1107 DO i = 1, nrep
1108 freq(i) = sqrt(abs(eigenval(ind(i)))*massunit)*vibfac
1109 END DO
1110
1111 !!! write information and output !!!
1112 IF (converged) THEN
1113 eigenval(:) = sign(1._dp, eigenval(:))*sqrt(abs(eigenval(:))*massunit)*vibfac
1114 ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1115 tmp_b = 0._dp
1116 ALLOCATE (tmp_s(3, ms_vib%mat_size))
1117 tmp_s = 0._dp
1118 IF (calc_intens) THEN
1119 ALLOCATE (intensities(ms_vib%mat_size))
1120 intensities = 0._dp
1121 END IF
1122 DO i = 1, ms_vib%mat_size
1123 DO j = 1, ms_vib%mat_size
1124 tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
1125 END DO
1126 tmp_b(:, i) = tmp_b(:, i)/norm2(tmp_b(:, i))
1127 END DO
1128 IF (calc_intens) THEN
1129 DO i = 1, ms_vib%mat_size
1130 DO j = 1, ms_vib%mat_size
1131 tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_h(j, i)
1132 END DO
1133 IF (calc_intens) intensities(i) = norm2(tmp_s(:, i))
1134 END DO
1135 END IF
1136 IF (calc_intens) THEN
1137 CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
1138 input, nrep, approx_h, eigenval, calc_intens, &
1139 intensities=intensities, logger=logger)
1140 ELSE
1141 CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
1142 input, nrep, approx_h, eigenval, calc_intens, logger=logger)
1143 END IF
1144 dump_only_positive = ms_vib%low_freq > 0.0_dp
1145 CALL write_vibrations_molden(input, particles, eigenval, tmp_b, intensities, calc_intens, &
1146 dump_only_positive=dump_only_positive, logger=logger, cell=cell)
1147 IF (calc_intens) THEN
1148 DEALLOCATE (intensities)
1149 END IF
1150 DEALLOCATE (tmp_b)
1151 DEALLOCATE (tmp_s)
1152 END IF
1153
1154 IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
1155 ms_vib, input, nrep, approx_h, eigenval, calc_intens, logger=logger)
1156
1157 DEALLOCATE (freq)
1158 DEALLOCATE (approx_h)
1159 DEALLOCATE (eigenval)
1160 DEALLOCATE (residuum)
1161 DEALLOCATE (ind)
1162
1163 END SUBROUTINE evaluate_h_update_b
1164
1165! **************************************************************************************************
1166!> \brief writes the output for a mode tracking calculation
1167!> \param ms_vib ...
1168!> \param nrep ...
1169!> \param mass ...
1170!> \param ncoord ...
1171!> \param approx_H ...
1172!> \param eigenval ...
1173!> \param ind ...
1174!> \param residuum ...
1175!> \param criteria ...
1176!> \author Florian Schiffmann 11.2007
1177! **************************************************************************************************
1178 SUBROUTINE select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1179
1180 TYPE(ms_vib_type) :: ms_vib
1181 INTEGER :: nrep
1182 REAL(kind=dp), DIMENSION(:) :: mass
1183 INTEGER :: ncoord
1184 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1185 REAL(kind=dp), DIMENSION(:) :: eigenval
1186 INTEGER, DIMENSION(:) :: ind
1187 REAL(kind=dp), DIMENSION(:, :) :: residuum
1188 REAL(kind=dp), DIMENSION(2, nrep), OPTIONAL :: criteria
1189
1190 INTEGER :: i, j, jj, k
1191 REAL(kind=dp) :: my_val, norm
1192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1193 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1194
1195 ALLOCATE (tmp(ms_vib%mat_size))
1196
1197 SELECT CASE (ms_vib%select_id)
1198 CASE (1)
1199 my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1200 DO i = 1, ms_vib%mat_size
1201 tmp(i) = abs(my_val - eigenval(i))
1202 END DO
1203 CALL sort(tmp, (ms_vib%mat_size), ind)
1204 residuum = 0._dp
1205 DO j = 1, nrep
1206 DO i = 1, ms_vib%mat_size
1207 residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1208 END DO
1209 END DO
1210 CASE (2)
1211 CALL get_vibs_in_range(ms_vib, approx_h, eigenval, residuum, nrep, ind)
1212 CASE (3)
1213
1214 ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1215 tmp_b = 0._dp
1216
1217 DO i = 1, ms_vib%mat_size
1218 DO j = 1, ms_vib%mat_size
1219 tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
1220 END DO
1221 tmp_b(:, i) = tmp_b(:, i)/norm2(tmp_b(:, i))
1222 END DO
1223 tmp = 0._dp
1224 DO i = 1, ms_vib%mat_size
1225 DO j = 1, SIZE(ms_vib%inv_atoms)
1226 DO k = 1, 3
1227 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1228 tmp(i) = tmp(i) + sqrt(tmp_b(jj, i)**2)
1229 END DO
1230 END DO
1231 IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1232 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
1233 ELSE
1234 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) <= ms_vib%inv_range(1)) tmp(i) = 0._dp
1235 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) >= ms_vib%inv_range(2)) tmp(i) = 0._dp
1236 END IF
1237 END DO
1238 tmp(:) = -tmp(:)
1239 CALL sort(tmp, (ms_vib%mat_size), ind)
1240 residuum(:, :) = 0._dp
1241
1242 DO j = 1, nrep
1243 DO i = 1, ms_vib%mat_size
1244 residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1245 END DO
1246 END DO
1247 DEALLOCATE (tmp_b)
1248 END SELECT
1249
1250 DO j = 1, nrep
1251 DO i = 1, ms_vib%mat_size
1252 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1253 END DO
1254 END DO
1255 IF (PRESENT(criteria)) THEN
1256 DO i = 1, nrep
1257 criteria(1, i) = maxval((residuum(:, i)))
1258 criteria(2, i) = norm2(residuum(:, i))
1259 END DO
1260 END IF
1261
1262 DO i = 1, nrep
1263 norm = norm2(residuum(:, i))
1264 residuum(:, i) = residuum(:, i)/norm
1265 END DO
1266
1267 DO k = 1, 10
1268 DO j = 1, nrep
1269 DO i = 1, ms_vib%mat_size
1270 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1271 residuum(:, j) = residuum(:, j)/norm2(residuum(:, j))
1272 END DO
1273 IF (nrep > 1) THEN
1274 DO i = 1, nrep
1275 IF (i /= j) THEN
1276 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), residuum(:, i))*residuum(:, i)
1277 residuum(:, j) = residuum(:, j)/norm2(residuum(:, j))
1278 END IF
1279 END DO
1280 END IF
1281 END DO
1282 END DO
1283 ms_vib%b_vec = residuum
1284 DEALLOCATE (tmp)
1285 END SUBROUTINE select_vector
1286
1287! **************************************************************************************************
1288!> \brief writes the output for a mode tracking calculation
1289!> \param iw ...
1290!> \param converged ...
1291!> \param freq ...
1292!> \param criter ...
1293!> \param ms_vib ...
1294!> \param input ...
1295!> \param nrep ...
1296!> \param approx_H ...
1297!> \param eigenval ...
1298!> \param calc_intens ...
1299!> \param intensities ...
1300!> \param logger ...
1301!> \author Florian Schiffmann 11.2007
1302! **************************************************************************************************
1303 SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1304 approx_H, eigenval, calc_intens, intensities, logger)
1305
1306 INTEGER :: iw
1307 LOGICAL :: converged
1308 REAL(kind=dp), DIMENSION(:) :: freq
1309 REAL(kind=dp), DIMENSION(:, :) :: criter
1310 TYPE(ms_vib_type) :: ms_vib
1311 TYPE(section_vals_type), POINTER :: input
1312 INTEGER :: nrep
1313 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1314 REAL(kind=dp), DIMENSION(:) :: eigenval
1315 LOGICAL :: calc_intens
1316 REAL(kind=dp), DIMENSION(:), OPTIONAL :: intensities
1317 TYPE(cp_logger_type), POINTER :: logger
1318
1319 REAL(kind=dp), PARAMETER :: eps_criter1 = 1.0e-07_dp, &
1320 eps_criter2 = 1.0e-06_dp
1321
1322 INTEGER :: i, j, msunit
1323 REAL(kind=dp) :: crit_a, crit_b, fint, gintval
1324 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: residuum
1325 TYPE(section_vals_type), POINTER :: ms_vib_section
1326
1327 ms_vib_section => section_vals_get_subs_vals(input, &
1328 "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1329
1330 fint = 42.255_dp*massunit*debye**2*bohr**2
1331
1332 IF (converged) THEN
1333 IF (iw > 0) THEN
1334 WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1335 DO i = 1, nrep
1336 WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1337 END DO
1338 ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1339 WRITE (iw, '( /, 1X, 79("-") )')
1340 WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1341 IF (PRESENT(intensities)) THEN
1342 WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1343 ELSE
1344 WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1345 END IF
1346 DO i = 1, SIZE(ms_vib%b_mat, 2)
1347 residuum = 0._dp
1348 DO j = 1, SIZE(ms_vib%b_mat, 2)
1349 residuum(:) = residuum(:) + approx_h(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1350 END DO
1351 DO j = 1, ms_vib%mat_size
1352 residuum(:) = residuum(:) - dot_product(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1353 END DO
1354 crit_a = maxval(residuum(:))
1355 crit_b = norm2(residuum)
1356 IF (PRESENT(intensities)) THEN
1357 gintval = fint*intensities(i)**2
1358 IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1359 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1360 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1361 ELSE
1362 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1363 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1364 END IF
1365 ELSE
1366 IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1367 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1368 'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1369 ELSE
1370 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1371 'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1372 END IF
1373 END IF
1374 END DO
1375 DEALLOCATE (residuum)
1376
1377 msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1378 "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1379 file_status="REPLACE", file_form="UNFORMATTED", &
1380 file_action="WRITE")
1381
1382 IF (msunit > 0) THEN
1383 WRITE (unit=msunit) ms_vib%mat_size
1384 WRITE (unit=msunit) ms_vib%b_mat
1385 WRITE (unit=msunit) ms_vib%s_mat
1386 IF (calc_intens) WRITE (unit=msunit) ms_vib%dip_deriv
1387 END IF
1388
1389 CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1390 "PRINT%MS_RESTART")
1391 END IF
1392 ELSE
1393 IF (iw > 0) THEN
1394 msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1395 "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1396 file_status="REPLACE", file_form="UNFORMATTED", &
1397 file_action="WRITE")
1398
1399 IF (msunit > 0) THEN
1400 WRITE (unit=msunit) ms_vib%mat_size
1401 WRITE (unit=msunit) ms_vib%b_mat
1402 WRITE (unit=msunit) ms_vib%s_mat
1403 IF (calc_intens) WRITE (unit=msunit) ms_vib%dip_deriv
1404 END IF
1405
1406 CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1407 "PRINT%MS_RESTART")
1408
1409 WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1410 DO i = 1, nrep
1411 IF (criter(1, i) <= eps_criter1 .AND. (criter(2, i)) <= eps_criter2) THEN
1412 WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1413 ELSE
1414 WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 NOT CONVERGED"
1415 END IF
1416 END DO
1417 END IF
1418 END IF
1419
1420 END SUBROUTINE ms_out
1421
1422! **************************************************************************************************
1423!> \brief ...
1424!> \param ms_vib ...
1425!> \param approx_H ...
1426!> \param eigenval ...
1427!> \param residuum ...
1428!> \param nrep ...
1429!> \param ind ...
1430!> \author Florian Schiffmann 11.2007
1431! **************************************************************************************************
1432 SUBROUTINE get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1433
1434 TYPE(ms_vib_type) :: ms_vib
1435 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1436 REAL(kind=dp), DIMENSION(:) :: eigenval
1437 REAL(kind=dp), DIMENSION(:, :) :: residuum
1438 INTEGER :: nrep
1439 INTEGER, DIMENSION(:) :: ind
1440
1441 INTEGER :: count1, count2, i, j
1442 INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1443 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1444 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1445 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1446 REAL(kind=dp), DIMENSION(2) :: myrange
1447
1448 myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1449 count1 = 0
1450 count2 = 0
1451 residuum = 0.0_dp
1452 ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1453 ALLOCATE (map1(SIZE(eigenval), 2))
1454 ALLOCATE (tmp(SIZE(eigenval)))
1455 DO i = 1, SIZE(eigenval)
1456 IF (abs(eigenval(i) - myrange(1)) + abs(eigenval(i) - myrange(2)) <= &
1457 abs(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1458 count1 = count1 + 1
1459 map1(count1, 1) = i
1460 ELSE
1461 count2 = count2 + 1
1462 map1(count2, 2) = i
1463 tmp(count2) = min(abs(eigenval(i) - myrange(1)), abs(eigenval(i) - myrange(2)))
1464 END IF
1465 END DO
1466
1467 IF (count1 == nrep) THEN
1468 DO j = 1, count1
1469 DO i = 1, ms_vib%mat_size
1470 residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))*(ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1471 ind(j) = map1(j, 1)
1472 END DO
1473 END DO
1474 ELSE IF (count1 > nrep) THEN
1475 ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1476 ALLOCATE (tmp1(count1))
1477 ALLOCATE (map2(count1))
1478 tmp_resid = 0._dp
1479 DO j = 1, count1
1480 DO i = 1, ms_vib%mat_size
1481 tmp_resid(:, j) = tmp_resid(:, j) + approx_h(i, map1(j, 1))* &
1482 (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1483 END DO
1484 END DO
1485
1486 DO j = 1, count1
1487 DO i = 1, ms_vib%mat_size
1488 tmp_resid(:, j) = tmp_resid(:, j) - dot_product(tmp_resid(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1489 END DO
1490 tmp(j) = maxval(tmp_resid(:, j))
1491 END DO
1492 CALL sort(tmp, count1, map2)
1493 DO j = 1, nrep
1494 residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1495 ind(j) = map1(map2(count1 + 1 - j), 1)
1496 END DO
1497 DEALLOCATE (tmp_resid)
1498 DEALLOCATE (tmp1)
1499 DEALLOCATE (map2)
1500 ELSE IF (count1 < nrep) THEN
1501
1502 ALLOCATE (map2(count2))
1503 IF (count1 /= 0) THEN
1504 DO j = 1, count1
1505 DO i = 1, ms_vib%mat_size
1506 residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))* &
1507 (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1508 END DO
1509 ind(j) = map1(j, 1)
1510 END DO
1511 END IF
1512 CALL sort(tmp, count2, map2)
1513 DO j = 1, nrep - count1
1514 DO i = 1, ms_vib%mat_size
1515 residuum(:, count1 + j) = residuum(:, count1 + j) + approx_h(i, map1(map2(j), 2)) &
1516 *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1517 END DO
1518 ind(count1 + j) = map1(map2(j), 2)
1519 END DO
1520
1521 DEALLOCATE (map2)
1522 END IF
1523
1524 DEALLOCATE (map1)
1525 DEALLOCATE (tmp)
1526
1527 END SUBROUTINE get_vibs_in_range
1528END 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:225
cp_result_methods::get_results
Definition cp_result_methods.F:43
bfgs_optimizer
Routines for Geometry optimization using BFGS algorithm.
Definition bfgs_optimizer.F:14
bfgs_optimizer::bfgs_read_restart_header
subroutine, public bfgs_read_restart_header(unit, stored_method, stored_ndf, stored_type, tagged, reference_cell, cell_frame)
Read the metadata of a tagged BFGS restart, or identify a legacy raw Hessian.
Definition bfgs_optimizer.F:677
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_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:1156
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, file_access)
...
Definition cp_output_handling.F:935
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:18
input_constants::bfgs_direct
integer, parameter, public bfgs_direct
Definition input_constants.F:569
input_constants::ms_guess_restart
integer, parameter, public ms_guess_restart
Definition input_constants.F:891
input_constants::ms_guess_restart_vec
integer, parameter, public ms_guess_restart_vec
Definition input_constants.F:891
input_constants::bfgs_inverse
integer, parameter, public bfgs_inverse
Definition input_constants.F:569
input_constants::ms_guess_molden
integer, parameter, public ms_guess_molden
Definition input_constants.F:891
input_constants::ms_guess_atomic
integer, parameter, public ms_guess_atomic
Definition input_constants.F:891
input_constants::ms_guess_bfgs
integer, parameter, public ms_guess_bfgs
Definition input_constants.F:891
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::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:383
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mode_selective
Module performing a mdoe selective vibrational analysis.
Definition mode_selective.F:20
mode_selective::ms_vb_anal
subroutine, public ms_vb_anal(input, rep_env, para_env, globenv, particles, nrep, calc_intens, dx, output_unit, logger, cell)
Module performing a vibrational analysis.
Definition mode_selective.F:111
molden_utils
Functions handling the MOLDEN format. Split from mode_selective.
Definition molden_utils.F:12
molden_utils::write_vibrations_molden
subroutine, public write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, dump_only_positive, logger, list, cell)
writes the output for vibrational analysis in MOLDEN format
Definition molden_utils.F:626
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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
particle_types::particle_type
Definition particle_types.F:35
replica_types::replica_env_type
keeps replicated information about the replicas
Definition replica_types.F:96