Line data Source code
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 = -1
64 : INTEGER :: select_id = -1
65 : INTEGER, DIMENSION(:), POINTER :: inv_atoms => NULL()
66 : REAL(KIND=dp) :: eps(2) = 0.0_dp
67 : REAL(KIND=dp) :: sel_freq = 0.0_dp
68 : REAL(KIND=dp) :: low_freq = 0.0_dp
69 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: b_vec => NULL()
70 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: delta_vec => NULL()
71 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ms_force => NULL()
72 : REAL(KIND=dp), DIMENSION(:), POINTER :: eig_bfgs => NULL()
73 : REAL(KIND=dp), DIMENSION(:), POINTER :: f_range => NULL()
74 : REAL(KIND=dp), DIMENSION(:), POINTER :: inv_range => NULL()
75 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_b => NULL()
76 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_r => NULL()
77 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: b_mat => NULL()
78 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_deriv => NULL()
79 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hes_bfgs => NULL()
80 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_mat => NULL()
81 : INTEGER :: initial_guess = -1
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 24 : 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 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mass, pos0
120 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: tmp_deriv
121 24 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: tmp_dip
122 : TYPE(ms_vib_type) :: ms_vib
123 :
124 24 : CALL timeset(routineN, handle)
125 24 : converged = .FALSE.
126 24 : natoms = SIZE(particles)
127 24 : ncoord = 3*natoms
128 96 : ALLOCATE (mass(3*natoms))
129 886 : DO i = 1, natoms
130 3472 : DO j = 1, 3
131 2586 : mass((i - 1)*3 + j) = particles(i)%atomic_kind%mass
132 3448 : mass((i - 1)*3 + j) = SQRT(mass((i - 1)*3 + j))
133 : END DO
134 : END DO
135 : ! Allocate working arrays
136 96 : ALLOCATE (ms_vib%delta_vec(ncoord, nrep))
137 72 : ALLOCATE (ms_vib%b_vec(ncoord, nrep))
138 72 : ALLOCATE (ms_vib%step_r(nrep))
139 48 : ALLOCATE (ms_vib%step_b(nrep))
140 24 : IF (calc_intens) THEN
141 20 : description = '[DIPOLE]'
142 80 : ALLOCATE (tmp_dip(nrep, 3, 2))
143 60 : 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 24 : calc_intens, logger)
150 24 : ncoord = 3*natoms
151 72 : ALLOCATE (pos0(ncoord))
152 96 : ALLOCATE (ms_vib%ms_force(ncoord, nrep))
153 886 : DO i = 1, natoms
154 3472 : DO j = 1, 3
155 3448 : pos0((i - 1)*3 + j) = particles((i))%r(j)
156 : END DO
157 : END DO
158 162 : ncoord = 3*natoms
159 : DO
160 23178 : ms_vib%ms_force = HUGE(0.0_dp)
161 336 : DO i = 1, nrep
162 23178 : DO j = 1, ncoord
163 23016 : 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 162 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
167 :
168 336 : DO i = 1, nrep
169 174 : IF (calc_intens) THEN
170 : CALL get_results(results=rep_env%results(i)%results, &
171 : description=description, &
172 150 : n_rep=ip1)
173 : CALL get_results(results=rep_env%results(i)%results, &
174 : description=description, &
175 : values=tmp_dip(i, :, 1), &
176 150 : nval=ip1)
177 : END IF
178 23178 : DO j = 1, ncoord
179 23016 : ms_vib%ms_force(j, i) = rep_env%f(j, i)
180 : END DO
181 : END DO
182 336 : DO i = 1, nrep
183 23178 : DO j = 1, ncoord
184 23016 : 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 162 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
188 162 : IF (calc_intens) THEN
189 300 : DO i = 1, nrep
190 : CALL get_results(results=rep_env%results(i)%results, &
191 : description=description, &
192 150 : n_rep=ip1)
193 : CALL get_results(results=rep_env%results(i)%results, &
194 : description=description, &
195 : values=tmp_dip(i, :, 2), &
196 150 : nval=ip1)
197 900 : 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 162 : output_unit, logger)
207 162 : IF (converged) EXIT
208 162 : IF (calc_intens) THEN
209 390 : ALLOCATE (tmp_deriv(3, ms_vib%mat_size))
210 3730 : tmp_deriv = ms_vib%dip_deriv
211 130 : DEALLOCATE (ms_vib%dip_deriv)
212 390 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
213 3730 : ms_vib%dip_deriv(:, 1:ms_vib%mat_size) = tmp_deriv(:, 1:ms_vib%mat_size)
214 130 : DEALLOCATE (tmp_deriv)
215 : END IF
216 : END DO
217 24 : DEALLOCATE (ms_vib%ms_force)
218 24 : DEALLOCATE (pos0)
219 24 : DEALLOCATE (ms_vib%step_r)
220 24 : DEALLOCATE (ms_vib%step_b)
221 24 : DEALLOCATE (ms_vib%b_vec)
222 24 : DEALLOCATE (ms_vib%delta_vec)
223 24 : DEALLOCATE (mass)
224 24 : DEALLOCATE (ms_vib%b_mat)
225 24 : DEALLOCATE (ms_vib%s_mat)
226 24 : IF (ms_vib%select_id == 3) THEN
227 8 : DEALLOCATE (ms_vib%inv_atoms)
228 : END IF
229 24 : IF (ASSOCIATED(ms_vib%eig_bfgs)) THEN
230 0 : DEALLOCATE (ms_vib%eig_bfgs)
231 : END IF
232 24 : IF (ASSOCIATED(ms_vib%hes_bfgs)) THEN
233 0 : DEALLOCATE (ms_vib%hes_bfgs)
234 : END IF
235 24 : IF (calc_intens) THEN
236 20 : DEALLOCATE (ms_vib%dip_deriv)
237 20 : DEALLOCATE (tmp_dip)
238 : END IF
239 24 : CALL timestop(handle)
240 72 : 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 24 : SUBROUTINE MS_initial_moves(para_env, nrep, input, globenv, ms_vib, particles, &
257 24 : 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 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map_atoms
275 24 : 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 24 : CALL timeset(routineN, handle)
281 24 : NULLIFY (ms_vib%eig_bfgs, ms_vib%f_range, ms_vib%hes_bfgs, ms_vib%inv_range)
282 24 : ms_vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
283 24 : CALL section_vals_val_get(ms_vib_section, "INITIAL_GUESS", i_val=guess)
284 24 : CALL section_vals_val_get(ms_vib_section, "EPS_MAX_VAL", r_val=ms_vib%eps(1))
285 24 : CALL section_vals_val_get(ms_vib_section, "EPS_NORM", r_val=ms_vib%eps(2))
286 24 : CALL section_vals_val_get(ms_vib_section, "RANGE", n_rep_val=n_rep_val)
287 24 : ms_vib%select_id = 0
288 24 : IF (n_rep_val .NE. 0) THEN
289 2 : CALL section_vals_val_get(ms_vib_section, "RANGE", r_vals=ms_vib%f_range)
290 2 : IF (ms_vib%f_range(1) .GT. ms_vib%f_range(2)) THEN
291 0 : my_val = ms_vib%f_range(2)
292 0 : ms_vib%f_range(2) = ms_vib%f_range(1)
293 0 : ms_vib%f_range(1) = my_val
294 : END IF
295 2 : ms_vib%select_id = 2
296 : END IF
297 24 : CALL section_vals_val_get(ms_vib_section, "FREQUENCY", r_val=ms_vib%sel_freq)
298 24 : CALL section_vals_val_get(ms_vib_section, "LOWEST_FREQUENCY", r_val=ms_vib%low_freq)
299 24 : IF (ms_vib%sel_freq .GT. 0._dp) ms_vib%select_id = 1
300 24 : involved_at_section => section_vals_get_subs_vals(ms_vib_section, "INVOLVED_ATOMS")
301 24 : CALL section_vals_get(involved_at_section, explicit=do_involved_atoms)
302 24 : IF (do_involved_atoms) THEN
303 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", n_rep_val=n_rep_val)
304 8 : jj = 0
305 16 : DO k = 1, n_rep_val
306 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=k, i_vals=tmplist)
307 32 : DO j = 1, SIZE(tmplist)
308 24 : jj = jj + 1
309 : END DO
310 : END DO
311 8 : IF (jj .GE. 1) THEN
312 8 : natoms = jj
313 24 : ALLOCATE (ms_vib%inv_atoms(natoms))
314 8 : jj = 0
315 16 : DO m = 1, n_rep_val
316 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=m, i_vals=tmplist)
317 32 : DO j = 1, SIZE(tmplist)
318 24 : ms_vib%inv_atoms(j) = tmplist(j)
319 : END DO
320 : END DO
321 8 : ms_vib%select_id = 3
322 : END IF
323 8 : CALL section_vals_val_get(involved_at_section, "RANGE", n_rep_val=n_rep_val)
324 8 : IF (n_rep_val .NE. 0) THEN
325 0 : CALL section_vals_val_get(involved_at_section, "RANGE", r_vals=ms_vib%inv_range)
326 0 : IF (ms_vib%inv_range(1) .GT. ms_vib%inv_range(2)) THEN
327 0 : ms_vib%inv_range(2) = my_val
328 0 : ms_vib%inv_range(2) = ms_vib%inv_range(1)
329 0 : ms_vib%inv_range(1) = my_val
330 : END IF
331 : END IF
332 : END IF
333 24 : IF (ms_vib%select_id == 0) &
334 0 : CPABORT("no frequency, range or involved atoms specified ")
335 24 : ionode = para_env%is_source()
336 12 : SELECT CASE (guess)
337 : CASE (ms_guess_atomic)
338 12 : ms_vib%initial_guess = 1
339 12 : CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
340 12 : jj = 0
341 22 : DO k = 1, n_rep_val
342 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
343 42 : DO j = 1, SIZE(tmplist)
344 30 : jj = jj + 1
345 : END DO
346 : END DO
347 12 : IF (jj < 1) THEN
348 2 : natoms = SIZE(particles)
349 6 : ALLOCATE (map_atoms(natoms))
350 14 : DO j = 1, natoms
351 14 : map_atoms(j) = j
352 : END DO
353 : ELSE
354 10 : natoms = jj
355 30 : ALLOCATE (map_atoms(natoms))
356 10 : jj = 0
357 20 : DO m = 1, n_rep_val
358 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
359 40 : DO j = 1, SIZE(tmplist)
360 30 : 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 526 : ms_vib%b_vec = 0._dp
367 526 : ms_vib%delta_vec = 0._dp
368 12 : jj = 0
369 :
370 28 : DO i = 1, nrep
371 56 : DO j = 1, natoms
372 176 : DO k = 1, 3
373 120 : jj = (map_atoms(j) - 1)*3 + k
374 160 : ms_vib%b_vec(jj, i) = ABS(globenv%gaussian_rng_stream%next())
375 : END DO
376 : END DO
377 514 : norm = SQRT(DOT_PRODUCT(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
378 526 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
379 : END DO
380 :
381 12 : IF (nrep .GT. 1) THEN
382 44 : DO k = 1, 10
383 124 : DO j = 1, nrep
384 280 : DO i = 1, nrep
385 240 : IF (i .NE. j) THEN
386 : ms_vib%b_vec(:, j) = &
387 1520 : 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 1520 : 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 12 : ms_vib%mat_size = 0
397 474 : DO i = 1, SIZE(ms_vib%b_vec, 1)
398 972 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
399 : END DO
400 : CASE (ms_guess_bfgs)
401 :
402 4 : ms_vib%initial_guess = 2
403 4 : CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
404 4 : ms_vib%mat_size = 0
405 :
406 : CASE (ms_guess_restart_vec)
407 :
408 4 : ms_vib%initial_guess = 3
409 4 : ncoord = 3*SIZE(particles)
410 4 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
411 :
412 4 : ms_vib%mat_size = 0
413 : CASE (ms_guess_restart)
414 0 : ms_vib%initial_guess = 4
415 0 : ncoord = 3*SIZE(particles)
416 0 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
417 :
418 : CASE (ms_guess_molden)
419 4 : ms_vib%initial_guess = 5
420 4 : ncoord = 3*SIZE(particles)
421 4 : CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
422 28 : ms_vib%mat_size = 0
423 : END SELECT
424 5324 : CALL para_env%bcast(ms_vib%b_vec)
425 5324 : CALL para_env%bcast(ms_vib%delta_vec)
426 52 : DO i = 1, nrep
427 2650 : ms_vib%step_r(i) = dx/SQRT(DOT_PRODUCT(ms_vib%delta_vec(:, i), ms_vib%delta_vec(:, i)))
428 2674 : 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 24 : CALL timestop(handle)
431 :
432 48 : 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 4 : 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 4 : INTEGER, DIMENSION(:), POINTER :: tmplist
457 : REAL(KIND=dp) :: my_val, norm
458 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
459 : TYPE(cp_logger_type), POINTER :: logger
460 :
461 8 : logger => cp_get_default_logger()
462 4 : output_unit = cp_logger_get_default_io_unit(logger)
463 :
464 4 : natoms = SIZE(particles)
465 4 : ncoord = 3*natoms
466 :
467 16 : ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
468 12 : ALLOCATE (ms_vib%eig_bfgs(ncoord))
469 :
470 4 : IF (para_env%is_source()) THEN
471 2 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
472 2 : IF (hes_filename == "") hes_filename = "HESSIAN"
473 : CALL open_file(file_name=hes_filename, file_status="OLD", &
474 2 : file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
475 6 : ALLOCATE (tmp(ncoord))
476 6 : ALLOCATE (tmplist(ncoord))
477 :
478 : ! should use the cp_fm_read_unformatted...
479 356 : DO i = 1, ncoord
480 356 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%hes_bfgs(:, i)
481 : END DO
482 2 : CALL close_file(hesunit)
483 2 : IF (output_unit > 0) THEN
484 2 : IF (stat /= 0) THEN
485 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading HESSIAN **"
486 : ELSE
487 : WRITE (output_unit, FMT="(/,T2,A)") &
488 2 : "*** Initial Hessian has been read successfully ***"
489 : END IF
490 : END IF
491 356 : DO i = 1, ncoord
492 63014 : DO j = 1, ncoord
493 63012 : ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
494 : END DO
495 : END DO
496 :
497 2 : CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
498 356 : tmp(:) = 0._dp
499 2 : IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
500 2 : 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 2 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
502 178 : DO i = 1, ncoord
503 178 : tmp(i) = ABS(my_val - ms_vib%eig_bfgs(i))
504 : END DO
505 1 : ELSE IF (ms_vib%select_id == 3) THEN
506 178 : DO i = 1, ncoord
507 531 : DO j = 1, SIZE(ms_vib%inv_atoms)
508 1593 : DO k = 1, 3
509 1062 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
510 1416 : tmp(i) = tmp(i) + SQRT(ms_vib%hes_bfgs(jj, i)**2)
511 : END DO
512 : END DO
513 178 : 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 178 : tmp(:) = -tmp(:)
516 : END IF
517 2 : CALL sort(tmp, ncoord, tmplist)
518 4 : DO i = 1, nrep
519 356 : ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
520 356 : norm = SQRT(DOT_PRODUCT(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
521 358 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
522 : END DO
523 356 : DO i = 1, SIZE(ms_vib%b_vec, 1)
524 710 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
525 : END DO
526 2 : DEALLOCATE (tmp)
527 2 : DEALLOCATE (tmplist)
528 : END IF
529 :
530 1428 : CALL para_env%bcast(ms_vib%b_vec)
531 1428 : CALL para_env%bcast(ms_vib%delta_vec)
532 :
533 4 : DEALLOCATE (ms_vib%hes_bfgs)
534 4 : DEALLOCATE (ms_vib%eig_bfgs)
535 4 : ms_vib%mat_size = 0
536 :
537 4 : 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 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
566 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
567 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H
568 : TYPE(cp_logger_type), POINTER :: logger
569 :
570 4 : logger => cp_get_default_logger()
571 4 : output_unit = cp_logger_get_default_io_unit(logger)
572 :
573 4 : natoms = SIZE(particles)
574 4 : ncoord = 3*natoms
575 4 : IF (calc_intens) THEN
576 4 : DEALLOCATE (ms_vib%dip_deriv)
577 : END IF
578 :
579 4 : IF (ionode) THEN
580 :
581 2 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
582 2 : 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 2 : unit_number=hesunit)
588 2 : READ (UNIT=hesunit, IOSTAT=stat) mat
589 2 : CPASSERT(stat == 0)
590 2 : ms_vib%mat_size = mat
591 : END IF
592 4 : CALL para_env%bcast(ms_vib%mat_size)
593 16 : ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
594 12 : ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
595 4 : IF (calc_intens) THEN
596 12 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
597 : END IF
598 4 : IF (ionode) THEN
599 2 : statint = 0
600 3384 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%b_mat
601 3384 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%s_mat
602 2 : IF (calc_intens) READ (UNIT=hesunit, IOSTAT=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
603 2 : IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART,", &
604 0 : "intensities are requested but not present in restart file **"
605 2 : CALL close_file(hesunit)
606 2 : IF (output_unit > 0) THEN
607 2 : IF (stat /= 0) THEN
608 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART **"
609 : ELSE
610 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** RESTART has been read successfully ***"
611 : END IF
612 : END IF
613 : END IF
614 13532 : CALL para_env%bcast(ms_vib%b_mat)
615 13532 : CALL para_env%bcast(ms_vib%s_mat)
616 340 : IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
617 16 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
618 12 : ALLOCATE (eigenval(ms_vib%mat_size))
619 12 : 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 4 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
623 4 : CALL diamat_all(approx_H, eigenval)
624 :
625 4 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, ms_vib%b_vec)
626 4 : IF (ms_vib%initial_guess .NE. 4) THEN
627 :
628 716 : ms_vib%b_vec = 0._dp
629 8 : DO i = 1, nrep
630 42 : DO j = 1, ms_vib%mat_size
631 6768 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_H(j, ind(i))*ms_vib%b_mat(:, j)
632 : END DO
633 1424 : 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 4 : DEALLOCATE (ms_vib%s_mat)
637 4 : DEALLOCATE (ms_vib%b_mat)
638 4 : IF (calc_intens) THEN
639 4 : DEALLOCATE (ms_vib%dip_deriv)
640 12 : ALLOCATE (ms_vib%dip_deriv(3, nrep))
641 : END IF
642 : END IF
643 4 : DEALLOCATE (approx_H)
644 4 : DEALLOCATE (eigenval)
645 4 : DEALLOCATE (ind)
646 8 : DO i = 1, nrep
647 716 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
648 : END DO
649 :
650 4 : 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 4 : 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 4 : INTEGER, DIMENSION(:), POINTER :: tmplist
678 : LOGICAL :: reading_vib
679 : REAL(KIND=dp) :: my_val, norm
680 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
681 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
682 8 : REAL(KIND=dp), DIMENSION(3, ncoord/3) :: pos
683 :
684 8 : output_unit = cp_logger_get_default_io_unit(logger)
685 :
686 4 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
687 4 : 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 0 : file_action="READ")
691 4 : IF (para_env%is_source()) THEN
692 :
693 2 : IF (ms_filename == "") THEN
694 0 : 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 2 : unit_number=iw)
701 : END IF
702 2 : info = ""
703 2 : READ (iw, *, IOSTAT=stat) info
704 2 : READ (iw, *, IOSTAT=stat) info
705 2 : istat = 0
706 2 : nvibs = 0
707 2 : reading_vib = .FALSE.
708 140 : DO
709 142 : READ (iw, *, IOSTAT=stat) info
710 142 : istat = istat + stat
711 142 : IF (TRIM(ADJUSTL(info)) == "[FR-COORD]") EXIT
712 :
713 140 : CPASSERT(stat == 0)
714 :
715 140 : IF (reading_vib) nvibs = nvibs + 1
716 140 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") reading_vib = .TRUE.
717 : END DO
718 2 : REWIND (iw)
719 2 : istat = 0
720 2 : READ (iw, *, IOSTAT=stat) info
721 2 : istat = istat + stat
722 2 : READ (iw, *, IOSTAT=stat) info
723 2 : istat = istat + stat
724 : ! Skip [Atoms] section
725 118 : DO
726 120 : READ (iw, *, IOSTAT=stat) info
727 120 : istat = istat + stat
728 120 : CPASSERT(stat == 0)
729 120 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") EXIT
730 : END DO
731 : ! Read frequencies and modes
732 6 : ALLOCATE (freq(nvibs))
733 8 : ALLOCATE (modes(ncoord, nvibs))
734 :
735 22 : DO i = 1, nvibs
736 20 : READ (iw, *, IOSTAT=stat) freq(i)
737 22 : istat = istat + stat
738 : END DO
739 2 : READ (iw, *) info
740 120 : DO i = 1, ncoord/3
741 118 : READ (iw, *, IOSTAT=stat) at_name, pos(:, i)
742 120 : istat = istat + stat
743 : END DO
744 2 : READ (iw, *) info
745 22 : DO i = 1, nvibs
746 20 : READ (iw, *) info
747 20 : istat = istat + stat
748 1202 : DO j = 1, ncoord/3
749 1180 : k = (j - 1)*3 + 1
750 1180 : READ (iw, *, IOSTAT=stat) modes(k:k + 2, i)
751 1200 : istat = istat + stat
752 : END DO
753 : END DO
754 2 : IF (ms_filename .NE. "") CALL close_file(iw)
755 2 : IF (output_unit > 0) THEN
756 2 : IF (istat /= 0) THEN
757 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MOLDEN file **"
758 : ELSE
759 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
760 : END IF
761 : END IF
762 : !!!!!!! select modes !!!!!!
763 4 : ALLOCATE (tmp(nvibs))
764 22 : tmp(:) = 0.0_dp
765 6 : ALLOCATE (tmplist(nvibs))
766 2 : IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
767 2 : IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
768 2 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
769 11 : DO i = 1, nvibs
770 11 : tmp(i) = ABS(my_val - freq(i))
771 : END DO
772 1 : ELSE IF (ms_vib%select_id == 3) THEN
773 11 : DO i = 1, nvibs
774 30 : DO j = 1, SIZE(ms_vib%inv_atoms)
775 90 : DO k = 1, 3
776 60 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
777 80 : tmp(i) = tmp(i) + SQRT(modes(jj, i)**2)
778 : END DO
779 : END DO
780 11 : IF (freq(i) .LE. 400._dp) tmp(i) = 0._dp
781 : END DO
782 11 : tmp(:) = -tmp(:)
783 : END IF
784 2 : CALL sort(tmp, nvibs, tmplist)
785 4 : DO i = 1, nrep
786 356 : ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
787 356 : norm = SQRT(DOT_PRODUCT(ms_vib%b_vec(:, i), ms_vib%b_vec(:, i)))
788 358 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
789 : END DO
790 4 : DO i = 1, nrep
791 358 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
792 : END DO
793 :
794 2 : DEALLOCATE (freq)
795 2 : DEALLOCATE (modes)
796 2 : DEALLOCATE (tmp)
797 2 : DEALLOCATE (tmplist)
798 :
799 : END IF
800 1428 : CALL para_env%bcast(ms_vib%b_vec)
801 1428 : CALL para_env%bcast(ms_vib%delta_vec)
802 :
803 4 : IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
804 0 : "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
805 4 : 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 162 : SUBROUTINE evaluate_H_update_b(rep_env, ms_vib, input, nrep, &
824 : particles, &
825 162 : 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 162 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
842 : LOGICAL :: dump_only_positive
843 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
844 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H, H_save, residuum, tmp_b, tmp_s
845 324 : REAL(KIND=dp), DIMENSION(2, nrep) :: criteria
846 162 : REAL(Kind=dp), DIMENSION(:), POINTER :: intensities
847 :
848 162 : natoms = SIZE(particles)
849 162 : ncoord = 3*natoms
850 162 : nrep = SIZE(rep_env%f, 2)
851 :
852 : !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
853 162 : IF (ms_vib%mat_size .NE. 0) THEN
854 :
855 690 : ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
856 414 : ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
857 :
858 153792 : tmp_b(:, :) = ms_vib%b_mat
859 153792 : tmp_s(:, :) = ms_vib%s_mat
860 :
861 138 : DEALLOCATE (ms_vib%b_mat)
862 138 : DEALLOCATE (ms_vib%s_mat)
863 : END IF
864 :
865 810 : ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
866 486 : ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
867 :
868 176832 : ms_vib%s_mat = 0.0_dp
869 :
870 22896 : DO i = 1, 3*natoms
871 22734 : IF (ms_vib%mat_size .NE. 0) THEN
872 172878 : DO j = 1, ms_vib%mat_size
873 152730 : ms_vib%b_mat(i, j) = tmp_b(i, j)
874 172878 : ms_vib%s_mat(i, j) = tmp_s(i, j)
875 : END DO
876 : END IF
877 45738 : DO j = 1, nrep
878 45576 : ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
879 : END DO
880 : END DO
881 :
882 162 : IF (ms_vib%mat_size .NE. 0) THEN
883 138 : DEALLOCATE (tmp_s)
884 138 : DEALLOCATE (tmp_b)
885 : END IF
886 :
887 162 : ms_vib%mat_size = ms_vib%mat_size + nrep
888 :
889 648 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
890 486 : ALLOCATE (H_save(ms_vib%mat_size, ms_vib%mat_size))
891 486 : ALLOCATE (eigenval(ms_vib%mat_size))
892 :
893 : !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
894 :
895 336 : DO i = 1, nrep
896 23178 : DO j = 1, 3*natoms
897 23016 : 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 162 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
903 14006 : H_save(:, :) = approx_H
904 :
905 162 : CALL diamat_all(approx_H, eigenval)
906 :
907 : !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
908 486 : ALLOCATE (ind(ms_vib%mat_size))
909 648 : ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
910 :
911 162 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
912 :
913 336 : DO i = 1, nrep
914 7950 : DO j = 1, natoms
915 30630 : DO k = 1, 3
916 22842 : jj = (j - 1)*3 + k
917 30456 : 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 336 : DO i = 1, nrep
923 23016 : ms_vib%step_r(i) = dx/SQRT(DOT_PRODUCT(ms_vib%delta_vec(:, i), ms_vib%delta_vec(:, i)))
924 23178 : 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 162 : converged = .FALSE.
927 : IF (MAXVAL(criteria(1, :)) .LE. ms_vib%eps(1) .AND. MAXVAL(criteria(2, :)) &
928 834 : .LE. ms_vib%eps(2) .OR. ms_vib%mat_size .GE. ncoord) converged = .TRUE.
929 486 : ALLOCATE (freq(nrep))
930 336 : DO i = 1, nrep
931 336 : freq(i) = SQRT(ABS(eigenval(ind(i)))*massunit)*vibfac
932 : END DO
933 :
934 : !!! write information and output !!!
935 162 : IF (converged) THEN
936 198 : eigenval(:) = SIGN(1._dp, eigenval(:))*SQRT(ABS(eigenval(:))*massunit)*vibfac
937 96 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
938 23040 : tmp_b = 0._dp
939 72 : ALLOCATE (tmp_s(3, ms_vib%mat_size))
940 720 : tmp_s = 0._dp
941 24 : IF (calc_intens) THEN
942 60 : ALLOCATE (intensities(ms_vib%mat_size))
943 170 : intensities = 0._dp
944 : END IF
945 198 : DO i = 1, ms_vib%mat_size
946 2268 : DO j = 1, ms_vib%mat_size
947 331218 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
948 : END DO
949 45882 : tmp_b(:, i) = tmp_b(:, i)/SQRT(DOT_PRODUCT(tmp_b(:, i), tmp_b(:, i)))
950 : END DO
951 24 : IF (calc_intens) THEN
952 170 : DO i = 1, ms_vib%mat_size
953 2100 : DO j = 1, ms_vib%mat_size
954 7950 : tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_H(j, i)
955 : END DO
956 620 : IF (calc_intens) intensities(i) = SQRT(DOT_PRODUCT(tmp_s(:, i), tmp_s(:, i)))
957 : END DO
958 : END IF
959 24 : 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 20 : intensities=intensities, logger=logger)
963 : ELSE
964 : CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
965 4 : input, nrep, approx_H, eigenval, calc_intens, logger=logger)
966 : END IF
967 24 : 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 24 : dump_only_positive=dump_only_positive, logger=logger)
970 24 : IF (calc_intens) THEN
971 20 : DEALLOCATE (intensities)
972 : END IF
973 24 : DEALLOCATE (tmp_b)
974 24 : DEALLOCATE (tmp_s)
975 : END IF
976 :
977 162 : IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
978 138 : ms_vib, input, nrep, approx_H, eigenval, calc_intens, logger=logger)
979 :
980 162 : DEALLOCATE (freq)
981 162 : DEALLOCATE (approx_H)
982 162 : DEALLOCATE (eigenval)
983 162 : DEALLOCATE (residuum)
984 162 : DEALLOCATE (ind)
985 :
986 324 : 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 166 : 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 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1016 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1017 :
1018 498 : ALLOCATE (tmp(ms_vib%mat_size))
1019 :
1020 282 : SELECT CASE (ms_vib%select_id)
1021 : CASE (1)
1022 116 : my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1023 1006 : DO i = 1, ms_vib%mat_size
1024 1006 : tmp(i) = ABS(my_val - eigenval(i))
1025 : END DO
1026 116 : CALL sort(tmp, (ms_vib%mat_size), ind)
1027 15892 : residuum = 0._dp
1028 238 : DO j = 1, nrep
1029 1152 : DO i = 1, ms_vib%mat_size
1030 144040 : 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 6 : CALL get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1035 : CASE (3)
1036 :
1037 176 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1038 39560 : tmp_b = 0._dp
1039 :
1040 266 : DO i = 1, ms_vib%mat_size
1041 1728 : DO j = 1, ms_vib%mat_size
1042 268290 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
1043 : END DO
1044 78854 : tmp_b(:, i) = tmp_b(:, i)/SQRT(DOT_PRODUCT(tmp_b(:, i), tmp_b(:, i)))
1045 : END DO
1046 266 : tmp = 0._dp
1047 266 : DO i = 1, ms_vib%mat_size
1048 666 : DO j = 1, SIZE(ms_vib%inv_atoms)
1049 1998 : DO k = 1, 3
1050 1332 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1051 1776 : tmp(i) = tmp(i) + SQRT(tmp_b(jj, i)**2)
1052 : END DO
1053 : END DO
1054 266 : IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1055 222 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) .LE. 400._dp) tmp(i) = 0._dp
1056 : ELSE
1057 0 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) .LE. ms_vib%inv_range(1)) tmp(i) = 0._dp
1058 0 : 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 266 : tmp(:) = -tmp(:)
1062 44 : CALL sort(tmp, (ms_vib%mat_size), ind)
1063 7876 : residuum(:, :) = 0._dp
1064 :
1065 88 : DO j = 1, nrep
1066 310 : DO i = 1, ms_vib%mat_size
1067 39560 : 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 210 : DEALLOCATE (tmp_b)
1071 : END SELECT
1072 :
1073 344 : DO j = 1, nrep
1074 1528 : DO i = 1, ms_vib%mat_size
1075 366822 : 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 166 : IF (PRESENT(criteria)) THEN
1079 336 : DO i = 1, nrep
1080 23190 : criteria(1, i) = MAXVAL((residuum(:, i)))
1081 23178 : criteria(2, i) = SQRT(DOT_PRODUCT(residuum(:, i), residuum(:, i)))
1082 : END DO
1083 : END IF
1084 :
1085 344 : DO i = 1, nrep
1086 23728 : norm = SQRT(DOT_PRODUCT(residuum(:, i), residuum(:, i)))
1087 23894 : residuum(:, i) = residuum(:, i)/norm
1088 : END DO
1089 :
1090 1826 : DO k = 1, 10
1091 3606 : DO j = 1, nrep
1092 13620 : DO i = 1, ms_vib%mat_size
1093 3666440 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1094 3668220 : residuum(:, j) = residuum(:, j)/SQRT(DOT_PRODUCT(residuum(:, j), residuum(:, j)))
1095 : END DO
1096 3440 : IF (nrep .GT. 1) THEN
1097 720 : DO i = 1, nrep
1098 720 : IF (i .NE. j) THEN
1099 4560 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), residuum(:, i))*residuum(:, i)
1100 4560 : 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 23894 : ms_vib%b_vec = residuum
1107 166 : DEALLOCATE (tmp)
1108 166 : 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 162 : SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1127 162 : 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 162 : 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 162 : "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1149 :
1150 162 : fint = 42.255_dp*massunit*debye**2*bohr**2
1151 :
1152 162 : IF (converged) THEN
1153 24 : IF (iw .GT. 0) THEN
1154 12 : WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1155 26 : DO i = 1, nrep
1156 26 : WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1157 : END DO
1158 36 : ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1159 12 : WRITE (iw, '( /, 1X, 79("-") )')
1160 12 : WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1161 12 : IF (PRESENT(intensities)) THEN
1162 10 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1163 : ELSE
1164 2 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1165 : END IF
1166 99 : DO i = 1, SIZE(ms_vib%b_mat, 2)
1167 11508 : residuum = 0._dp
1168 1134 : DO j = 1, SIZE(ms_vib%b_mat, 2)
1169 165609 : residuum(:) = residuum(:) + approx_H(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1170 : END DO
1171 1134 : DO j = 1, ms_vib%mat_size
1172 330084 : residuum(:) = residuum(:) - DOT_PRODUCT(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1173 : END DO
1174 11595 : crit_a = MAXVAL(residuum(:))
1175 11508 : crit_b = SQRT(DOT_PRODUCT(residuum, residuum))
1176 99 : IF (PRESENT(intensities)) THEN
1177 75 : gintval = fint*intensities(i)**2
1178 75 : IF (crit_a .LE. ms_vib%eps(1) .AND. crit_b .LE. ms_vib%eps(2)) THEN
1179 28 : 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 26 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1181 : ELSE
1182 47 : 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 47 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1184 : END IF
1185 : ELSE
1186 12 : IF (crit_a .LE. ms_vib%eps(1) .AND. crit_b .LE. ms_vib%eps(2)) THEN
1187 12 : IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1188 6 : 'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1189 : ELSE
1190 0 : IF (eigenval(i) .GT. ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1191 0 : 'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1192 : END IF
1193 : END IF
1194 : END DO
1195 12 : 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 12 : file_action="WRITE")
1201 :
1202 12 : IF (msunit > 0) THEN
1203 12 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%mat_size
1204 11520 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%b_mat
1205 11520 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%s_mat
1206 312 : 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 12 : "PRINT%MS_RESTART")
1211 : END IF
1212 : ELSE
1213 138 : 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 69 : file_action="WRITE")
1218 :
1219 69 : IF (msunit > 0) THEN
1220 69 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%mat_size
1221 76896 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%b_mat
1222 76896 : WRITE (UNIT=msunit, IOSTAT=stat) ms_vib%s_mat
1223 1869 : 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 69 : "PRINT%MS_RESTART")
1228 :
1229 69 : WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1230 142 : DO i = 1, nrep
1231 142 : IF (criter(1, i) .LE. 1E-7 .AND. (criter(2, i)) .LE. 1E-6) THEN
1232 1 : WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1233 : ELSE
1234 72 : 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 162 : 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 6 : 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 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1263 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1264 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1265 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1266 : REAL(KIND=dp), DIMENSION(2) :: myrange
1267 :
1268 18 : myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1269 6 : count1 = 0
1270 6 : count2 = 0
1271 126 : residuum = 0.0_dp
1272 6 : ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1273 18 : ALLOCATE (map1(SIZE(eigenval), 2))
1274 18 : ALLOCATE (tmp(SIZE(eigenval)))
1275 30 : DO i = 1, SIZE(eigenval)
1276 24 : IF (ABS(eigenval(i) - myrange(1)) + ABS(eigenval(i) - myrange(2)) .LE. &
1277 6 : ABS(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1278 0 : count1 = count1 + 1
1279 0 : map1(count1, 1) = i
1280 : ELSE
1281 24 : count2 = count2 + 1
1282 24 : map1(count2, 2) = i
1283 24 : tmp(count2) = MIN(ABS(eigenval(i) - myrange(1)), ABS(eigenval(i) - myrange(2)))
1284 : END IF
1285 : END DO
1286 :
1287 6 : IF (count1 .EQ. nrep) THEN
1288 0 : DO j = 1, count1
1289 0 : DO i = 1, ms_vib%mat_size
1290 0 : 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 0 : ind(j) = map1(j, 1)
1292 : END DO
1293 : END DO
1294 6 : ELSE IF (count1 .GT. nrep) THEN
1295 0 : ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1296 0 : ALLOCATE (tmp1(count1))
1297 0 : ALLOCATE (map2(count1))
1298 0 : tmp_resid = 0._dp
1299 0 : DO j = 1, count1
1300 0 : DO i = 1, ms_vib%mat_size
1301 : tmp_resid(:, j) = tmp_resid(:, j) + approx_H(i, map1(j, 1))* &
1302 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1303 : END DO
1304 : END DO
1305 :
1306 0 : DO j = 1, count1
1307 0 : DO i = 1, ms_vib%mat_size
1308 0 : 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 0 : tmp(j) = MAXVAL(tmp_resid(:, j))
1311 : END DO
1312 0 : CALL sort(tmp, count1, map2)
1313 0 : DO j = 1, nrep
1314 0 : residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1315 0 : ind(j) = map1(map2(count1 + 1 - j), 1)
1316 : END DO
1317 0 : DEALLOCATE (tmp_resid)
1318 0 : DEALLOCATE (tmp1)
1319 0 : DEALLOCATE (map2)
1320 6 : ELSE IF (count1 .LT. nrep) THEN
1321 :
1322 18 : ALLOCATE (map2(count2))
1323 6 : IF (count1 .NE. 0) THEN
1324 0 : DO j = 1, count1
1325 0 : DO i = 1, ms_vib%mat_size
1326 : residuum(:, j) = residuum(:, j) + approx_H(i, map1(j, 1))* &
1327 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1328 : END DO
1329 0 : ind(j) = map1(j, 1)
1330 : END DO
1331 : END IF
1332 6 : CALL sort(tmp, count2, map2)
1333 18 : DO j = 1, nrep - count1
1334 60 : DO i = 1, ms_vib%mat_size
1335 : residuum(:, count1 + j) = residuum(:, count1 + j) + approx_H(i, map1(map2(j), 2)) &
1336 492 : *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1337 : END DO
1338 18 : ind(count1 + j) = map1(map2(j), 2)
1339 : END DO
1340 :
1341 6 : DEALLOCATE (map2)
1342 : END IF
1343 :
1344 6 : DEALLOCATE (map1)
1345 6 : DEALLOCATE (tmp)
1346 :
1347 6 : END SUBROUTINE get_vibs_in_range
1348 0 : END MODULE mode_selective
|
