(git:d3d49ac)
Loading...
Searching...
No Matches
mode_selective.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Module performing a mdoe selective vibrational analysis
10!> \note
11!> Numerical accuracy for parallel runs:
12!> Each replica starts the SCF run from the one optimized
13!> in a previous run. It may happen then energies and derivatives
14!> of a serial run and a parallel run could be slightly different
15!> 'cause of a different starting density matrix.
16!> Exact results are obtained using:
17!> EXTRAPOLATION USE_GUESS in QS section (Teo 08.2006)
18!> \author Florian Schiffmann 08.2006
19! **************************************************************************************************
20MODULE mode_selective
21 USE 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 /= 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) > 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 > 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 >= 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 /= 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) > 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) THEN
334 cpabort("no frequency, range or involved atoms specified ")
335 END IF
336 ionode = para_env%is_source()
337 SELECT CASE (guess)
338 CASE (ms_guess_atomic)
339 ms_vib%initial_guess = 1
340 CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
341 jj = 0
342 DO k = 1, n_rep_val
343 CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
344 DO j = 1, SIZE(tmplist)
345 jj = jj + 1
346 END DO
347 END DO
348 IF (jj < 1) THEN
349 natoms = SIZE(particles)
350 ALLOCATE (map_atoms(natoms))
351 DO j = 1, natoms
352 map_atoms(j) = j
353 END DO
354 ELSE
355 natoms = jj
356 ALLOCATE (map_atoms(natoms))
357 jj = 0
358 DO m = 1, n_rep_val
359 CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
360 DO j = 1, SIZE(tmplist)
361 map_atoms(j) = tmplist(j)
362 END DO
363 END DO
364 END IF
365
366 ! apply random displacement along the mass weighted nuclear cartesian coordinates
367 ms_vib%b_vec = 0._dp
368 ms_vib%delta_vec = 0._dp
369 jj = 0
370
371 DO i = 1, nrep
372 DO j = 1, natoms
373 DO k = 1, 3
374 jj = (map_atoms(j) - 1)*3 + k
375 ms_vib%b_vec(jj, i) = abs(globenv%gaussian_rng_stream%next())
376 END DO
377 END DO
378 norm = norm2(ms_vib%b_vec(:, i))
379 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
380 END DO
381
382 IF (nrep > 1) THEN
383 DO k = 1, 10
384 DO j = 1, nrep
385 DO i = 1, nrep
386 IF (i /= j) THEN
387 ms_vib%b_vec(:, j) = &
388 ms_vib%b_vec(:, j) - dot_product(ms_vib%b_vec(:, j), ms_vib%b_vec(:, i))*ms_vib%b_vec(:, i)
389 ms_vib%b_vec(:, j) = &
390 ms_vib%b_vec(:, j)/norm2(ms_vib%b_vec(:, j))
391 END IF
392 END DO
393 END DO
394 END DO
395 END IF
396
397 ms_vib%mat_size = 0
398 DO i = 1, SIZE(ms_vib%b_vec, 1)
399 ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
400 END DO
401 CASE (ms_guess_bfgs)
402
403 ms_vib%initial_guess = 2
404 CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
405 ms_vib%mat_size = 0
406
407 CASE (ms_guess_restart_vec)
408
409 ms_vib%initial_guess = 3
410 ncoord = 3*SIZE(particles)
411 CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
412
413 ms_vib%mat_size = 0
414 CASE (ms_guess_restart)
415 ms_vib%initial_guess = 4
416 ncoord = 3*SIZE(particles)
417 CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
418
419 CASE (ms_guess_molden)
420 ms_vib%initial_guess = 5
421 ncoord = 3*SIZE(particles)
422 CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
423 ms_vib%mat_size = 0
424 END SELECT
425 CALL para_env%bcast(ms_vib%b_vec)
426 CALL para_env%bcast(ms_vib%delta_vec)
427 DO i = 1, nrep
428 ms_vib%step_r(i) = dx/norm2(ms_vib%delta_vec(:, i))
429 ms_vib%step_b(i) = norm2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
430 END DO
431 CALL timestop(handle)
432
433 END SUBROUTINE ms_initial_moves
434
435! **************************************************************************************************
436!> \brief ...
437!> \param ms_vib_section ...
438!> \param ms_vib ...
439!> \param particles ...
440!> \param mass ...
441!> \param para_env ...
442!> \param nrep ...
443!> \author Florian Schiffmann 11.2007
444! **************************************************************************************************
445 SUBROUTINE bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
446
447 TYPE(section_vals_type), POINTER :: ms_vib_section
448 TYPE(ms_vib_type) :: ms_vib
449 TYPE(particle_type), DIMENSION(:), POINTER :: particles
450 REAL(kind=dp), DIMENSION(:) :: mass
451 TYPE(mp_para_env_type), POINTER :: para_env
452 INTEGER :: nrep
453
454 CHARACTER(LEN=default_path_length) :: hes_filename
455 INTEGER :: hesunit, i, istat, j, jj, k, natoms, &
456 ncoord, output_unit, stat
457 INTEGER, DIMENSION(:), POINTER :: tmplist
458 REAL(kind=dp) :: my_val, norm
459 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
460 TYPE(cp_logger_type), POINTER :: logger
461
462 logger => cp_get_default_logger()
463 output_unit = cp_logger_get_default_io_unit(logger)
464
465 natoms = SIZE(particles)
466 ncoord = 3*natoms
467
468 ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
469 ALLOCATE (ms_vib%eig_bfgs(ncoord))
470
471 IF (para_env%is_source()) THEN
472 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
473 IF (hes_filename == "") hes_filename = "HESSIAN"
474 CALL open_file(file_name=hes_filename, file_status="OLD", &
475 file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
476 ALLOCATE (tmp(ncoord))
477 ALLOCATE (tmplist(ncoord))
478
479 ! should use the cp_fm_read_unformatted...
480 istat = 0
481 DO i = 1, ncoord
482 READ (unit=hesunit, iostat=stat) ms_vib%hes_bfgs(:, i)
483 istat = istat + stat
484 END DO
485 CALL close_file(hesunit)
486 IF (output_unit > 0) THEN
487 IF (istat /= 0) THEN
488 WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading HESSIAN **"
489 ELSE
490 WRITE (output_unit, fmt="(/,T2,A)") &
491 "*** Initial Hessian has been read successfully ***"
492 END IF
493 END IF
494 DO i = 1, ncoord
495 DO j = 1, ncoord
496 ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
497 END DO
498 END DO
499
500 CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
501 tmp(:) = 0._dp
502 IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
503 IF (ms_vib%select_id == 2) my_val = (((ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp)/vibfac)**2/massunit
504 IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
505 DO i = 1, ncoord
506 tmp(i) = abs(my_val - ms_vib%eig_bfgs(i))
507 END DO
508 ELSE IF (ms_vib%select_id == 3) THEN
509 DO i = 1, ncoord
510 DO j = 1, SIZE(ms_vib%inv_atoms)
511 DO k = 1, 3
512 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
513 tmp(i) = tmp(i) + sqrt(ms_vib%hes_bfgs(jj, i)**2)
514 END DO
515 END DO
516 IF ((sign(1._dp, ms_vib%eig_bfgs(i))*sqrt(abs(ms_vib%eig_bfgs(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
517 END DO
518 tmp(:) = -tmp(:)
519 END IF
520 CALL sort(tmp, ncoord, tmplist)
521 DO i = 1, nrep
522 ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
523 norm = norm2(ms_vib%b_vec(:, i))
524 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
525 END DO
526 DO i = 1, SIZE(ms_vib%b_vec, 1)
527 ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
528 END DO
529 DEALLOCATE (tmp)
530 DEALLOCATE (tmplist)
531 END IF
532
533 CALL para_env%bcast(ms_vib%b_vec)
534 CALL para_env%bcast(ms_vib%delta_vec)
535
536 DEALLOCATE (ms_vib%hes_bfgs)
537 DEALLOCATE (ms_vib%eig_bfgs)
538 ms_vib%mat_size = 0
539
540 END SUBROUTINE bfgs_guess
541
542! **************************************************************************************************
543!> \brief ...
544!> \param ms_vib_section ...
545!> \param para_env ...
546!> \param ms_vib ...
547!> \param mass ...
548!> \param ionode ...
549!> \param particles ...
550!> \param nrep ...
551!> \param calc_intens ...
552!> \author Florian Schiffmann 11.2007
553! **************************************************************************************************
554 SUBROUTINE rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
555
556 TYPE(section_vals_type), POINTER :: ms_vib_section
557 TYPE(mp_para_env_type), POINTER :: para_env
558 TYPE(ms_vib_type) :: ms_vib
559 REAL(kind=dp), DIMENSION(:) :: mass
560 LOGICAL :: ionode
561 TYPE(particle_type), DIMENSION(:), POINTER :: particles
562 INTEGER :: nrep
563 LOGICAL :: calc_intens
564
565 CHARACTER(LEN=default_path_length) :: ms_filename
566 INTEGER :: hesunit, i, j, mat, natoms, ncoord, &
567 output_unit, stat, statint
568 INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
569 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
570 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h
571 TYPE(cp_logger_type), POINTER :: logger
572
573 logger => cp_get_default_logger()
574 output_unit = cp_logger_get_default_io_unit(logger)
575
576 natoms = SIZE(particles)
577 ncoord = 3*natoms
578 IF (calc_intens) THEN
579 DEALLOCATE (ms_vib%dip_deriv)
580 END IF
581
582 IF (ionode) THEN
583
584 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
585 IF (ms_filename == "") ms_filename = "MS_RESTART"
586 CALL open_file(file_name=ms_filename, &
587 file_status="UNKNOWN", &
588 file_form="UNFORMATTED", &
589 file_action="READ", &
590 unit_number=hesunit)
591 READ (unit=hesunit, iostat=stat) mat
592 cpassert(stat == 0)
593 ms_vib%mat_size = mat
594 END IF
595 CALL para_env%bcast(ms_vib%mat_size)
596 ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
597 ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
598 IF (calc_intens) THEN
599 ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
600 END IF
601 IF (ionode) THEN
602 statint = 0
603 READ (unit=hesunit) ms_vib%b_mat
604 READ (unit=hesunit, iostat=stat) ms_vib%s_mat
605 IF (stat /= 0 .AND. output_unit > 0) THEN
606 WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART **"
607 END IF
608 IF (calc_intens) THEN
609 READ (unit=hesunit, iostat=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
610 IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MS_RESTART,", &
611 "intensities are requested but not present in restart file **"
612 END IF
613 CALL close_file(hesunit)
614 IF (stat == 0 .AND. statint == 0 .AND. output_unit > 0) THEN
615 WRITE (output_unit, fmt="(/,T2,A)") "*** MS_RESTART has been read successfully ***"
616 END IF
617 END IF
618 CALL para_env%bcast(ms_vib%b_mat)
619 CALL para_env%bcast(ms_vib%s_mat)
620 IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
621 ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
622 ALLOCATE (eigenval(ms_vib%mat_size))
623 ALLOCATE (ind(ms_vib%mat_size))
624
625 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), &
626 ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
627 CALL diamat_all(approx_h, eigenval)
628
629 CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, ms_vib%b_vec)
630 IF (ms_vib%initial_guess /= 4) THEN
631
632 ms_vib%b_vec = 0._dp
633 DO i = 1, nrep
634 DO j = 1, ms_vib%mat_size
635 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_h(j, ind(i))*ms_vib%b_mat(:, j)
636 END DO
637 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm2(ms_vib%b_vec(:, i))
638 END DO
639
640 DEALLOCATE (ms_vib%s_mat)
641 DEALLOCATE (ms_vib%b_mat)
642 IF (calc_intens) THEN
643 DEALLOCATE (ms_vib%dip_deriv)
644 ALLOCATE (ms_vib%dip_deriv(3, nrep))
645 END IF
646 END IF
647 DEALLOCATE (approx_h)
648 DEALLOCATE (eigenval)
649 DEALLOCATE (ind)
650 DO i = 1, nrep
651 ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
652 END DO
653
654 END SUBROUTINE rest_guess
655
656! **************************************************************************************************
657!> \brief ...
658!> \param ms_vib_section ...
659!> \param input ...
660!> \param para_env ...
661!> \param ms_vib ...
662!> \param mass ...
663!> \param ncoord ...
664!> \param nrep ...
665!> \param logger ...
666!> \author Florian Schiffmann 11.2007
667! **************************************************************************************************
668 SUBROUTINE molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
669 TYPE(section_vals_type), POINTER :: ms_vib_section, input
670 TYPE(mp_para_env_type), POINTER :: para_env
671 TYPE(ms_vib_type) :: ms_vib
672 REAL(kind=dp), DIMENSION(:) :: mass
673 INTEGER :: ncoord, nrep
674 TYPE(cp_logger_type), POINTER :: logger
675
676 CHARACTER(LEN=2) :: at_name
677 CHARACTER(LEN=default_path_length) :: ms_filename
678 CHARACTER(LEN=max_line_length) :: info
679 INTEGER :: i, istat, iw, j, jj, k, nvibs, &
680 output_molden, output_unit, stat
681 INTEGER, DIMENSION(:), POINTER :: tmplist
682 LOGICAL :: reading_vib
683 REAL(kind=dp) :: my_val, norm
684 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
685 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
686 REAL(kind=dp), DIMENSION(3, ncoord/3) :: pos
687
688 output_unit = cp_logger_get_default_io_unit(logger)
689
690 CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
691 IF (ms_filename == "") output_molden = &
692 cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
693 extension=".mol", file_status='UNKNOWN', &
694 file_action="READ")
695 IF (para_env%is_source()) THEN
696
697 IF (ms_filename == "") THEN
698 iw = output_molden
699 ELSE
700 CALL open_file(file_name=trim(ms_filename), &
701 file_status="UNKNOWN", &
702 file_form="FORMATTED", &
703 file_action="READ", &
704 unit_number=iw)
705 END IF
706 info = ""
707 READ (iw, *) info
708 READ (iw, *) info
709 istat = 0
710 nvibs = 0
711 reading_vib = .false.
712 DO
713 READ (iw, *, iostat=stat) info
714 istat = istat + stat
715 IF (trim(adjustl(info)) == "[FR-COORD]") EXIT
716
717 cpassert(stat == 0)
718
719 IF (reading_vib) nvibs = nvibs + 1
720 IF (trim(adjustl(info)) == "[FREQ]") reading_vib = .true.
721 END DO
722 rewind(iw)
723 istat = 0
724 READ (iw, *, iostat=stat) info
725 istat = istat + stat
726 READ (iw, *, iostat=stat) info
727 istat = istat + stat
728 ! Skip [Atoms] section
729 DO
730 READ (iw, *, iostat=stat) info
731 istat = istat + stat
732 cpassert(stat == 0)
733 IF (trim(adjustl(info)) == "[FREQ]") EXIT
734 END DO
735 ! Read frequencies and modes
736 ALLOCATE (freq(nvibs))
737 ALLOCATE (modes(ncoord, nvibs))
738
739 DO i = 1, nvibs
740 READ (iw, *, iostat=stat) freq(i)
741 istat = istat + stat
742 END DO
743 READ (iw, *) info
744 DO i = 1, ncoord/3
745 READ (iw, *, iostat=stat) at_name, pos(:, i)
746 istat = istat + stat
747 END DO
748 READ (iw, *) info
749 DO i = 1, nvibs
750 READ (iw, *) info
751 istat = istat + stat
752 DO j = 1, ncoord/3
753 k = (j - 1)*3 + 1
754 READ (iw, *, iostat=stat) modes(k:k + 2, i)
755 istat = istat + stat
756 END DO
757 END DO
758 IF (ms_filename /= "") CALL close_file(iw)
759 IF (output_unit > 0) THEN
760 IF (istat /= 0) THEN
761 WRITE (output_unit, fmt="(/,T2,A)") "** Error while reading MOLDEN file **"
762 ELSE
763 WRITE (output_unit, fmt="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
764 END IF
765 END IF
766 !!!!!!! select modes !!!!!!
767 ALLOCATE (tmp(nvibs))
768 tmp(:) = 0.0_dp
769 ALLOCATE (tmplist(nvibs))
770 IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
771 IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
772 IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
773 DO i = 1, nvibs
774 tmp(i) = abs(my_val - freq(i))
775 END DO
776 ELSE IF (ms_vib%select_id == 3) THEN
777 DO i = 1, nvibs
778 DO j = 1, SIZE(ms_vib%inv_atoms)
779 DO k = 1, 3
780 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
781 tmp(i) = tmp(i) + sqrt(modes(jj, i)**2)
782 END DO
783 END DO
784 IF (freq(i) <= 400._dp) tmp(i) = 0._dp
785 END DO
786 tmp(:) = -tmp(:)
787 END IF
788 CALL sort(tmp, nvibs, tmplist)
789 DO i = 1, nrep
790 ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
791 norm = norm2(ms_vib%b_vec(:, i))
792 ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
793 END DO
794 DO i = 1, nrep
795 ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
796 END DO
797
798 DEALLOCATE (freq)
799 DEALLOCATE (modes)
800 DEALLOCATE (tmp)
801 DEALLOCATE (tmplist)
802
803 END IF
804 CALL para_env%bcast(ms_vib%b_vec)
805 CALL para_env%bcast(ms_vib%delta_vec)
806
807 IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
808 "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
809 END SUBROUTINE molden_guess
810
811! **************************************************************************************************
812!> \brief Davidson algorithm for to generate a approximate Hessian for mode
813!> selective vibrational analysis
814!> \param rep_env ...
815!> \param ms_vib ...
816!> \param input ...
817!> \param nrep ...
818!> \param particles ...
819!> \param mass ...
820!> \param converged ...
821!> \param dx ...
822!> \param calc_intens ...
823!> \param output_unit_ms ...
824!> \param logger ...
825!> \author Florian Schiffmann 11.2007
826! **************************************************************************************************
827 SUBROUTINE evaluate_h_update_b(rep_env, ms_vib, input, nrep, &
828 particles, &
829 mass, &
830 converged, dx, &
831 calc_intens, output_unit_ms, logger)
832 TYPE(replica_env_type), POINTER :: rep_env
833 TYPE(ms_vib_type) :: ms_vib
834 TYPE(section_vals_type), POINTER :: input
835 INTEGER :: nrep
836 TYPE(particle_type), DIMENSION(:), POINTER :: particles
837 REAL(kind=dp), DIMENSION(:) :: mass
838 LOGICAL :: converged
839 REAL(kind=dp) :: dx
840 LOGICAL :: calc_intens
841 INTEGER :: output_unit_ms
842 TYPE(cp_logger_type), POINTER :: logger
843
844 INTEGER :: i, j, jj, k, natoms, ncoord
845 INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
846 LOGICAL :: dump_only_positive
847 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
848 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_h, h_save, residuum, tmp_b, tmp_s
849 REAL(kind=dp), DIMENSION(2, nrep) :: criteria
850 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
851
852 natoms = SIZE(particles)
853 ncoord = 3*natoms
854 nrep = SIZE(rep_env%f, 2)
855
856 !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
857 IF (ms_vib%mat_size /= 0) THEN
858
859 ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
860 ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
861
862 tmp_b(:, :) = ms_vib%b_mat
863 tmp_s(:, :) = ms_vib%s_mat
864
865 DEALLOCATE (ms_vib%b_mat)
866 DEALLOCATE (ms_vib%s_mat)
867 END IF
868
869 ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
870 ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
871
872 ms_vib%s_mat = 0.0_dp
873
874 DO i = 1, 3*natoms
875 IF (ms_vib%mat_size /= 0) THEN
876 DO j = 1, ms_vib%mat_size
877 ms_vib%b_mat(i, j) = tmp_b(i, j)
878 ms_vib%s_mat(i, j) = tmp_s(i, j)
879 END DO
880 END IF
881 DO j = 1, nrep
882 ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
883 END DO
884 END DO
885
886 IF (ms_vib%mat_size /= 0) THEN
887 DEALLOCATE (tmp_s)
888 DEALLOCATE (tmp_b)
889 END IF
890
891 ms_vib%mat_size = ms_vib%mat_size + nrep
892
893 ALLOCATE (approx_h(ms_vib%mat_size, ms_vib%mat_size))
894 ALLOCATE (h_save(ms_vib%mat_size, ms_vib%mat_size))
895 ALLOCATE (eigenval(ms_vib%mat_size))
896
897 !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
898
899 DO i = 1, nrep
900 DO j = 1, 3*natoms
901 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))
902 END DO
903 END DO
904
905 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), &
906 ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_h, ms_vib%mat_size)
907 h_save(:, :) = approx_h
908
909 CALL diamat_all(approx_h, eigenval)
910
911 !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
912 ALLOCATE (ind(ms_vib%mat_size))
913 ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
914
915 CALL select_vector(ms_vib, nrep, mass, ncoord, approx_h, eigenval, ind, residuum, criteria)
916
917 DO i = 1, nrep
918 DO j = 1, natoms
919 DO k = 1, 3
920 jj = (j - 1)*3 + k
921 ms_vib%delta_vec(jj, i) = ms_vib%b_vec(jj, i)/mass(jj)
922 END DO
923 END DO
924 END DO
925
926 DO i = 1, nrep
927 ms_vib%step_r(i) = dx/norm2(ms_vib%delta_vec(:, i))
928 ms_vib%step_b(i) = norm2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
929 END DO
930 converged = .false.
931 IF (maxval(criteria(1, :)) <= ms_vib%eps(1) .AND. maxval(criteria(2, :)) &
932 <= ms_vib%eps(2) .OR. ms_vib%mat_size >= ncoord) converged = .true.
933 ALLOCATE (freq(nrep))
934 DO i = 1, nrep
935 freq(i) = sqrt(abs(eigenval(ind(i)))*massunit)*vibfac
936 END DO
937
938 !!! write information and output !!!
939 IF (converged) THEN
940 eigenval(:) = sign(1._dp, eigenval(:))*sqrt(abs(eigenval(:))*massunit)*vibfac
941 ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
942 tmp_b = 0._dp
943 ALLOCATE (tmp_s(3, ms_vib%mat_size))
944 tmp_s = 0._dp
945 IF (calc_intens) THEN
946 ALLOCATE (intensities(ms_vib%mat_size))
947 intensities = 0._dp
948 END IF
949 DO i = 1, ms_vib%mat_size
950 DO j = 1, ms_vib%mat_size
951 tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
952 END DO
953 tmp_b(:, i) = tmp_b(:, i)/norm2(tmp_b(:, i))
954 END DO
955 IF (calc_intens) THEN
956 DO i = 1, ms_vib%mat_size
957 DO j = 1, ms_vib%mat_size
958 tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_h(j, i)
959 END DO
960 IF (calc_intens) intensities(i) = norm2(tmp_s(:, i))
961 END DO
962 END IF
963 IF (calc_intens) THEN
964 CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
965 input, nrep, approx_h, eigenval, calc_intens, &
966 intensities=intensities, logger=logger)
967 ELSE
968 CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
969 input, nrep, approx_h, eigenval, calc_intens, logger=logger)
970 END IF
971 dump_only_positive = ms_vib%low_freq > 0.0_dp
972 CALL write_vibrations_molden(input, particles, eigenval, tmp_b, intensities, calc_intens, &
973 dump_only_positive=dump_only_positive, logger=logger)
974 IF (calc_intens) THEN
975 DEALLOCATE (intensities)
976 END IF
977 DEALLOCATE (tmp_b)
978 DEALLOCATE (tmp_s)
979 END IF
980
981 IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
982 ms_vib, input, nrep, approx_h, eigenval, calc_intens, logger=logger)
983
984 DEALLOCATE (freq)
985 DEALLOCATE (approx_h)
986 DEALLOCATE (eigenval)
987 DEALLOCATE (residuum)
988 DEALLOCATE (ind)
989
990 END SUBROUTINE evaluate_h_update_b
991
992! **************************************************************************************************
993!> \brief writes the output for a mode tracking calculation
994!> \param ms_vib ...
995!> \param nrep ...
996!> \param mass ...
997!> \param ncoord ...
998!> \param approx_H ...
999!> \param eigenval ...
1000!> \param ind ...
1001!> \param residuum ...
1002!> \param criteria ...
1003!> \author Florian Schiffmann 11.2007
1004! **************************************************************************************************
1005 SUBROUTINE select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1006
1007 TYPE(ms_vib_type) :: ms_vib
1008 INTEGER :: nrep
1009 REAL(kind=dp), DIMENSION(:) :: mass
1010 INTEGER :: ncoord
1011 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1012 REAL(kind=dp), DIMENSION(:) :: eigenval
1013 INTEGER, DIMENSION(:) :: ind
1014 REAL(kind=dp), DIMENSION(:, :) :: residuum
1015 REAL(kind=dp), DIMENSION(2, nrep), OPTIONAL :: criteria
1016
1017 INTEGER :: i, j, jj, k
1018 REAL(kind=dp) :: my_val, norm
1019 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1020 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1021
1022 ALLOCATE (tmp(ms_vib%mat_size))
1023
1024 SELECT CASE (ms_vib%select_id)
1025 CASE (1)
1026 my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1027 DO i = 1, ms_vib%mat_size
1028 tmp(i) = abs(my_val - eigenval(i))
1029 END DO
1030 CALL sort(tmp, (ms_vib%mat_size), ind)
1031 residuum = 0._dp
1032 DO j = 1, nrep
1033 DO i = 1, ms_vib%mat_size
1034 residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1035 END DO
1036 END DO
1037 CASE (2)
1038 CALL get_vibs_in_range(ms_vib, approx_h, eigenval, residuum, nrep, ind)
1039 CASE (3)
1040
1041 ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1042 tmp_b = 0._dp
1043
1044 DO i = 1, ms_vib%mat_size
1045 DO j = 1, ms_vib%mat_size
1046 tmp_b(:, i) = tmp_b(:, i) + approx_h(j, i)*ms_vib%b_mat(:, j)/mass(:)
1047 END DO
1048 tmp_b(:, i) = tmp_b(:, i)/norm2(tmp_b(:, i))
1049 END DO
1050 tmp = 0._dp
1051 DO i = 1, ms_vib%mat_size
1052 DO j = 1, SIZE(ms_vib%inv_atoms)
1053 DO k = 1, 3
1054 jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1055 tmp(i) = tmp(i) + sqrt(tmp_b(jj, i)**2)
1056 END DO
1057 END DO
1058 IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1059 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
1060 ELSE
1061 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) <= ms_vib%inv_range(1)) tmp(i) = 0._dp
1062 IF ((sign(1._dp, eigenval(i))*sqrt(abs(eigenval(i))*massunit)*vibfac) >= ms_vib%inv_range(2)) tmp(i) = 0._dp
1063 END IF
1064 END DO
1065 tmp(:) = -tmp(:)
1066 CALL sort(tmp, (ms_vib%mat_size), ind)
1067 residuum(:, :) = 0._dp
1068
1069 DO j = 1, nrep
1070 DO i = 1, ms_vib%mat_size
1071 residuum(:, j) = residuum(:, j) + approx_h(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1072 END DO
1073 END DO
1074 DEALLOCATE (tmp_b)
1075 END SELECT
1076
1077 DO j = 1, nrep
1078 DO i = 1, ms_vib%mat_size
1079 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1080 END DO
1081 END DO
1082 IF (PRESENT(criteria)) THEN
1083 DO i = 1, nrep
1084 criteria(1, i) = maxval((residuum(:, i)))
1085 criteria(2, i) = norm2(residuum(:, i))
1086 END DO
1087 END IF
1088
1089 DO i = 1, nrep
1090 norm = norm2(residuum(:, i))
1091 residuum(:, i) = residuum(:, i)/norm
1092 END DO
1093
1094 DO k = 1, 10
1095 DO j = 1, nrep
1096 DO i = 1, ms_vib%mat_size
1097 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1098 residuum(:, j) = residuum(:, j)/norm2(residuum(:, j))
1099 END DO
1100 IF (nrep > 1) THEN
1101 DO i = 1, nrep
1102 IF (i /= j) THEN
1103 residuum(:, j) = residuum(:, j) - dot_product(residuum(:, j), residuum(:, i))*residuum(:, i)
1104 residuum(:, j) = residuum(:, j)/norm2(residuum(:, j))
1105 END IF
1106 END DO
1107 END IF
1108 END DO
1109 END DO
1110 ms_vib%b_vec = residuum
1111 DEALLOCATE (tmp)
1112 END SUBROUTINE select_vector
1113
1114! **************************************************************************************************
1115!> \brief writes the output for a mode tracking calculation
1116!> \param iw ...
1117!> \param converged ...
1118!> \param freq ...
1119!> \param criter ...
1120!> \param ms_vib ...
1121!> \param input ...
1122!> \param nrep ...
1123!> \param approx_H ...
1124!> \param eigenval ...
1125!> \param calc_intens ...
1126!> \param intensities ...
1127!> \param logger ...
1128!> \author Florian Schiffmann 11.2007
1129! **************************************************************************************************
1130 SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1131 approx_H, eigenval, calc_intens, intensities, logger)
1132
1133 INTEGER :: iw
1134 LOGICAL :: converged
1135 REAL(kind=dp), DIMENSION(:) :: freq
1136 REAL(kind=dp), DIMENSION(:, :) :: criter
1137 TYPE(ms_vib_type) :: ms_vib
1138 TYPE(section_vals_type), POINTER :: input
1139 INTEGER :: nrep
1140 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1141 REAL(kind=dp), DIMENSION(:) :: eigenval
1142 LOGICAL :: calc_intens
1143 REAL(kind=dp), DIMENSION(:), OPTIONAL :: intensities
1144 TYPE(cp_logger_type), POINTER :: logger
1145
1146 INTEGER :: i, j, msunit
1147 REAL(kind=dp) :: crit_a, crit_b, fint, gintval
1148 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: residuum
1149 TYPE(section_vals_type), POINTER :: ms_vib_section
1150
1151 ms_vib_section => section_vals_get_subs_vals(input, &
1152 "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1153
1154 fint = 42.255_dp*massunit*debye**2*bohr**2
1155
1156 IF (converged) THEN
1157 IF (iw > 0) THEN
1158 WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1159 DO i = 1, nrep
1160 WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1161 END DO
1162 ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1163 WRITE (iw, '( /, 1X, 79("-") )')
1164 WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1165 IF (PRESENT(intensities)) THEN
1166 WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1167 ELSE
1168 WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1169 END IF
1170 DO i = 1, SIZE(ms_vib%b_mat, 2)
1171 residuum = 0._dp
1172 DO j = 1, SIZE(ms_vib%b_mat, 2)
1173 residuum(:) = residuum(:) + approx_h(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1174 END DO
1175 DO j = 1, ms_vib%mat_size
1176 residuum(:) = residuum(:) - dot_product(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1177 END DO
1178 crit_a = maxval(residuum(:))
1179 crit_b = norm2(residuum)
1180 IF (PRESENT(intensities)) THEN
1181 gintval = fint*intensities(i)**2
1182 IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1183 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1184 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1185 ELSE
1186 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1187 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1188 END IF
1189 ELSE
1190 IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1191 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1192 'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1193 ELSE
1194 IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1195 'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1196 END IF
1197 END IF
1198 END DO
1199 DEALLOCATE (residuum)
1200
1201 msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1202 "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1203 file_status="REPLACE", file_form="UNFORMATTED", &
1204 file_action="WRITE")
1205
1206 IF (msunit > 0) THEN
1207 WRITE (unit=msunit) ms_vib%mat_size
1208 WRITE (unit=msunit) ms_vib%b_mat
1209 WRITE (unit=msunit) ms_vib%s_mat
1210 IF (calc_intens) WRITE (unit=msunit) ms_vib%dip_deriv
1211 END IF
1212
1213 CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1214 "PRINT%MS_RESTART")
1215 END IF
1216 ELSE
1217 IF (iw > 0) THEN
1218 msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1219 "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1220 file_status="REPLACE", file_form="UNFORMATTED", &
1221 file_action="WRITE")
1222
1223 IF (msunit > 0) THEN
1224 WRITE (unit=msunit) ms_vib%mat_size
1225 WRITE (unit=msunit) ms_vib%b_mat
1226 WRITE (unit=msunit) ms_vib%s_mat
1227 IF (calc_intens) WRITE (unit=msunit) ms_vib%dip_deriv
1228 END IF
1229
1230 CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1231 "PRINT%MS_RESTART")
1232
1233 WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1234 DO i = 1, nrep
1235 IF (criter(1, i) <= 1e-7 .AND. (criter(2, i)) <= 1e-6) THEN
1236 WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1237 ELSE
1238 WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 NOT CONVERGED"
1239 END IF
1240 END DO
1241 END IF
1242 END IF
1243
1244 END SUBROUTINE ms_out
1245
1246! **************************************************************************************************
1247!> \brief ...
1248!> \param ms_vib ...
1249!> \param approx_H ...
1250!> \param eigenval ...
1251!> \param residuum ...
1252!> \param nrep ...
1253!> \param ind ...
1254!> \author Florian Schiffmann 11.2007
1255! **************************************************************************************************
1256 SUBROUTINE get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1257
1258 TYPE(ms_vib_type) :: ms_vib
1259 REAL(kind=dp), DIMENSION(:, :) :: approx_h
1260 REAL(kind=dp), DIMENSION(:) :: eigenval
1261 REAL(kind=dp), DIMENSION(:, :) :: residuum
1262 INTEGER :: nrep
1263 INTEGER, DIMENSION(:) :: ind
1264
1265 INTEGER :: count1, count2, i, j
1266 INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1267 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1268 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1269 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1270 REAL(kind=dp), DIMENSION(2) :: myrange
1271
1272 myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1273 count1 = 0
1274 count2 = 0
1275 residuum = 0.0_dp
1276 ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1277 ALLOCATE (map1(SIZE(eigenval), 2))
1278 ALLOCATE (tmp(SIZE(eigenval)))
1279 DO i = 1, SIZE(eigenval)
1280 IF (abs(eigenval(i) - myrange(1)) + abs(eigenval(i) - myrange(2)) <= &
1281 abs(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1282 count1 = count1 + 1
1283 map1(count1, 1) = i
1284 ELSE
1285 count2 = count2 + 1
1286 map1(count2, 2) = i
1287 tmp(count2) = min(abs(eigenval(i) - myrange(1)), abs(eigenval(i) - myrange(2)))
1288 END IF
1289 END DO
1290
1291 IF (count1 == nrep) THEN
1292 DO j = 1, count1
1293 DO i = 1, ms_vib%mat_size
1294 residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))*(ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1295 ind(j) = map1(j, 1)
1296 END DO
1297 END DO
1298 ELSE IF (count1 > nrep) THEN
1299 ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1300 ALLOCATE (tmp1(count1))
1301 ALLOCATE (map2(count1))
1302 tmp_resid = 0._dp
1303 DO j = 1, count1
1304 DO i = 1, ms_vib%mat_size
1305 tmp_resid(:, j) = tmp_resid(:, j) + approx_h(i, map1(j, 1))* &
1306 (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1307 END DO
1308 END DO
1309
1310 DO j = 1, count1
1311 DO i = 1, ms_vib%mat_size
1312 tmp_resid(:, j) = tmp_resid(:, j) - dot_product(tmp_resid(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1313 END DO
1314 tmp(j) = maxval(tmp_resid(:, j))
1315 END DO
1316 CALL sort(tmp, count1, map2)
1317 DO j = 1, nrep
1318 residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1319 ind(j) = map1(map2(count1 + 1 - j), 1)
1320 END DO
1321 DEALLOCATE (tmp_resid)
1322 DEALLOCATE (tmp1)
1323 DEALLOCATE (map2)
1324 ELSE IF (count1 < nrep) THEN
1325
1326 ALLOCATE (map2(count2))
1327 IF (count1 /= 0) THEN
1328 DO j = 1, count1
1329 DO i = 1, ms_vib%mat_size
1330 residuum(:, j) = residuum(:, j) + approx_h(i, map1(j, 1))* &
1331 (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1332 END DO
1333 ind(j) = map1(j, 1)
1334 END DO
1335 END IF
1336 CALL sort(tmp, count2, map2)
1337 DO j = 1, nrep - count1
1338 DO i = 1, ms_vib%mat_size
1339 residuum(:, count1 + j) = residuum(:, count1 + j) + approx_h(i, map1(map2(j), 2)) &
1340 *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1341 END DO
1342 ind(count1 + j) = map1(map2(j), 2)
1343 END DO
1344
1345 DEALLOCATE (map2)
1346 END IF
1347
1348 DEALLOCATE (map1)
1349 DEALLOCATE (tmp)
1350
1351 END SUBROUTINE get_vibs_in_range
1352END MODULE mode_selective
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cp_result_methods::get_results
Definition cp_result_methods.F:43
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:311
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:122
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_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:1073
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ms_guess_restart
integer, parameter, public ms_guess_restart
Definition input_constants.F:872
input_constants::ms_guess_restart_vec
integer, parameter, public ms_guess_restart_vec
Definition input_constants.F:872
input_constants::ms_guess_molden
integer, parameter, public ms_guess_molden
Definition input_constants.F:872
input_constants::ms_guess_atomic
integer, parameter, public ms_guess_atomic
Definition input_constants.F:872
input_constants::ms_guess_bfgs
integer, parameter, public ms_guess_bfgs
Definition input_constants.F:872
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
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:605
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:367
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:750
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