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