74 a_bohr,
angstrom,
bohr,
boltzmann,
debye,
e_mass,
h_bar,
hertz,
kelvin,
kjmol,
massunit, &
82 #include "../base/base_uses.f90"
86 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'vibrational_analysis'
87 LOGICAL,
PARAMETER :: debug_this_module = .false.
101 SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
102 TYPE(section_vals_type),
POINTER :: input
103 TYPE(section_type),
POINTER :: input_declaration
104 TYPE(mp_para_env_type),
POINTER :: para_env
105 TYPE(global_environment_type),
POINTER :: globenv
107 CHARACTER(len=*),
PARAMETER :: routinen =
'vb_anal'
108 CHARACTER(LEN=1),
DIMENSION(3),
PARAMETER :: lab = (/
"X",
"Y",
"Z"/)
110 CHARACTER(LEN=default_string_length) :: description_d, description_p
111 INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, iounit, ip1, ip2, iparticle1, &
112 iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nrottrm, nvib, &
113 output_unit, output_unit_eig, prep, print_grrm, print_namd, print_scine, proc_dist_type
114 INTEGER,
DIMENSION(:),
POINTER :: clist, mlist
115 LOGICAL :: calc_intens, calc_thchdata, do_mode_tracking, intens_ir, intens_raman, &
116 keep_rotations, row_force, something_frozen
117 REAL(kind=
dp) :: a1, a2, a3, conver, dummy, dx, &
118 inertia(3), minimum_energy, norm, &
119 tc_press, tc_temp, tmp
120 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: h_eigval1, h_eigval2, heigvaldfull, &
121 konst, mass, pos0, rmass
122 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: hessian, hint1, hint2, hint2dfull, matm
123 REAL(kind=
dp),
DIMENSION(3) :: d_deriv, d_print
124 REAL(kind=
dp),
DIMENSION(3, 3) :: p_deriv, p_print
125 REAL(kind=
dp),
DIMENSION(:),
POINTER :: depol_p, depol_u, depp, depu, din, &
126 intensities_d, intensities_p, pin
127 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: d, dfull, dip_deriv, rottrm
128 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: polar_deriv, tmp_dip
129 REAL(kind=
dp),
DIMENSION(:, :, :, :),
POINTER :: tmp_polar
130 TYPE(cp_logger_type),
POINTER :: logger
131 TYPE(cp_subsys_type),
POINTER :: subsys
132 TYPE(f_env_type),
POINTER :: f_env
133 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
134 TYPE(replica_env_type),
POINTER :: rep_env
135 TYPE(section_vals_type),
POINTER :: force_env_section, &
136 mode_tracking_section, print_section, &
139 CALL timeset(routinen, handle)
140 NULLIFY (d, rottrm, logger, subsys, f_env, particles, rep_env, intensities_d, intensities_p, &
141 vib_section, print_section, depol_p, depol_u)
147 "PROGRAM_RUN_INFO", &
156 file_status=
"REPLACE", &
157 file_action=
"WRITE", &
159 file_form=
"UNFORMATTED")
172 tc_press = tc_press*
pascal
175 intens_raman = .false.
179 nrep = max(1, para_env%num_pe/prep)
180 prep = para_env%num_pe/nrep
188 input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=row_force)
189 IF (
ASSOCIATED(rep_env))
THEN
192 particles => subsys%particles%els
194 IF (do_mode_tracking)
THEN
195 CALL ms_vb_anal(input, rep_env, para_env, globenv, particles, &
196 nrep, calc_intens, dx, output_unit, logger)
199 CALL get_moving_atoms(force_env=f_env%force_env, ilist=mlist)
200 something_frozen =
SIZE(particles) .NE.
SIZE(mlist)
203 ALLOCATE (clist(ncoord))
204 ALLOCATE (mass(natoms))
205 ALLOCATE (pos0(ncoord))
206 ALLOCATE (hessian(ncoord, ncoord))
207 IF (calc_intens)
THEN
208 description_d =
'[DIPOLE]'
209 ALLOCATE (tmp_dip(ncoord, 3, 2))
211 description_p =
'[POLAR]'
212 ALLOCATE (tmp_polar(ncoord, 3, 3, 2))
218 clist((i - 1)*3 + 1) = (imap - 1)*3 + 1
219 clist((i - 1)*3 + 2) = (imap - 1)*3 + 2
220 clist((i - 1)*3 + 3) = (imap - 1)*3 + 3
221 mass(i) = particles(imap)%atomic_kind%mass
222 cpassert(mass(i) > 0.0_dp)
223 mass(i) = sqrt(mass(i))
224 pos0((i - 1)*3 + 1) = particles(imap)%r(1)
225 pos0((i - 1)*3 + 2) = particles(imap)%r(2)
226 pos0((i - 1)*3 + 3) = particles(imap)%r(3)
232 IF (something_frozen)
THEN
234 ALLOCATE (rottrm(natoms*3, nrottrm))
236 CALL rot_ana(particles, rottrm, nrottrm, print_section, &
237 keep_rotations, mass_weighted=.true., natoms=natoms, inertia=inertia)
240 nvib = 3*natoms - nrottrm
243 ALLOCATE (d(3*natoms, 3*natoms))
245 ALLOCATE (d(3*natoms, nvib))
247 CALL build_d_matrix(rottrm, nrottrm, d, full=.false., &
252 hessian = huge(0.0_dp)
253 IF (output_unit > 0)
WRITE (output_unit,
'(/,T2,A)')
"VIB| Vibrational Analysis Info"
254 DO icoordp = 1, ncoord, nrep
259 rep_env%r(imap, j) = pos0(i)
261 IF (icoord + j <= ncoord)
THEN
262 imap = clist(icoord + j)
263 rep_env%r(imap, j) = rep_env%r(imap, j) + dx
269 IF (calc_intens)
THEN
270 IF (icoord + j <= ncoord)
THEN
272 description=description_d))
THEN
273 CALL get_results(results=rep_env%results(j)%results, &
274 description=description_d, &
276 CALL get_results(results=rep_env%results(j)%results, &
277 description=description_d, &
278 values=tmp_dip(icoord + j, :, 1), &
281 d_print(:) = tmp_dip(icoord + j, :, 1)
284 description=description_p))
THEN
285 CALL get_results(results=rep_env%results(j)%results, &
286 description=description_p, &
288 CALL get_results(results=rep_env%results(j)%results, &
289 description=description_p, &
290 values=tmp_polar(icoord + j, :, :, 1), &
292 intens_raman = .true.
293 p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
297 IF (icoord + j <= ncoord)
THEN
300 hessian(i, icoord + j) = rep_env%f(imap, j)
302 imap = clist(icoord + j)
304 IF (output_unit > 0)
THEN
306 IF (mod(imap, 3) /= 0) iparticle1 = iparticle1 + 1
307 WRITE (output_unit,
'(T2,A,I5,A,I5,3A)') &
308 "VIB| REPLICA Nr.", j,
"- Energy and Forces for particle:", &
309 iparticle1,
" coordinate: ", lab(imap - (iparticle1 - 1)*3), &
310 " + D"//trim(lab(imap - (iparticle1 - 1)*3))
311 WRITE (output_unit,
'(T2,A,T43,A,T57,F24.12)') &
312 "VIB|",
"Total energy:", rep_env%f(rep_env%ndim + 1, j)
313 WRITE (output_unit,
'(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
316 WRITE (output_unit,
'(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
317 particles(imap)%atomic_kind%name, &
318 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
321 WRITE (output_unit,
'(T3,A)')
'Dipole moment [Debye]'
322 WRITE (output_unit,
'(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
323 'X=', d_print(1)*
debye,
'Y=', d_print(2)*
debye,
'Z=', d_print(3)*
debye, &
324 'Total=', sqrt(sum(d_print(1:3)**2))*
debye
326 IF (intens_raman)
THEN
327 WRITE (output_unit,
'(T2,A)') &
328 'POLAR| Polarizability tensor [a.u.]'
329 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12))') &
330 'POLAR| xx,yy,zz', p_print(1, 1), p_print(2, 2), p_print(3, 3)
331 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12))') &
332 'POLAR| xy,xz,yz', p_print(1, 2), p_print(1, 3), p_print(2, 3)
333 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12),/)') &
334 'POLAR| yx,zx,zy', p_print(2, 1), p_print(3, 1), p_print(3, 2)
340 DO icoordm = 1, ncoord, nrep
345 rep_env%r(imap, j) = pos0(i)
347 IF (icoord + j <= ncoord)
THEN
348 imap = clist(icoord + j)
349 rep_env%r(imap, j) = rep_env%r(imap, j) - dx
355 IF (calc_intens)
THEN
356 IF (icoord + j <= ncoord)
THEN
357 k = (icoord + j + 2)/3
359 description=description_d))
THEN
360 CALL get_results(results=rep_env%results(j)%results, &
361 description=description_d, &
363 CALL get_results(results=rep_env%results(j)%results, &
364 description=description_d, &
365 values=tmp_dip(icoord + j, :, 2), &
367 tmp_dip(icoord + j, :, 1) = (tmp_dip(icoord + j, :, 1) - &
368 tmp_dip(icoord + j, :, 2))/(2.0_dp*dx*mass(k))
369 d_print(:) = tmp_dip(icoord + j, :, 1)
372 description=description_p))
THEN
373 CALL get_results(results=rep_env%results(j)%results, &
374 description=description_p, &
376 CALL get_results(results=rep_env%results(j)%results, &
377 description=description_p, &
378 values=tmp_polar(icoord + j, :, :, 2), &
380 tmp_polar(icoord + j, :, :, 1) = (tmp_polar(icoord + j, :, :, 1) - &
381 tmp_polar(icoord + j, :, :, 2))/(2.0_dp*dx*mass(k))
382 p_print(:, :) = tmp_polar(icoord + j, :, :, 1)
386 IF (icoord + j <= ncoord)
THEN
387 imap = clist(icoord + j)
389 IF (mod(imap, 3) /= 0) iparticle1 = iparticle1 + 1
391 IF (mod(icoord + j, 3) /= 0) ip1 = ip1 + 1
393 IF (output_unit > 0)
THEN
394 WRITE (output_unit,
'(T2,A,I5,A,I5,3A)') &
395 "VIB| REPLICA Nr.", j,
"- Energy and Forces for particle:", &
396 iparticle1,
" coordinate: ", lab(imap - (iparticle1 - 1)*3), &
397 " - D"//trim(lab(imap - (iparticle1 - 1)*3))
398 WRITE (output_unit,
'(T2,A,T43,A,T57,F24.12)') &
399 "VIB|",
"Total energy:", rep_env%f(rep_env%ndim + 1, j)
400 WRITE (output_unit,
'(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
403 WRITE (output_unit,
'(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
404 particles(imap)%atomic_kind%name, &
405 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
408 WRITE (output_unit,
'(T3,A)')
'Dipole moment [Debye]'
409 WRITE (output_unit,
'(T5,3(A,F14.8,1X),T60,A,T67,F14.8)') &
410 'X=', d_print(1)*
debye,
'Y=', d_print(2)*
debye,
'Z=', d_print(3)*
debye, &
411 'Total=', sqrt(sum(d_print(1:3)**2))*
debye
413 IF (intens_raman)
THEN
414 WRITE (output_unit,
'(T2,A)') &
415 'POLAR| Polarizability tensor [a.u.]'
416 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12))') &
417 'POLAR| xx,yy,zz', p_print(1, 1), p_print(2, 2), p_print(3, 3)
418 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12))') &
419 'POLAR| xy,xz,yz', p_print(1, 2), p_print(1, 3), p_print(2, 3)
420 WRITE (output_unit,
'(T2,A,T24,3(1X,F18.12),/)') &
421 'POLAR| yx,zx,zy', p_print(2, 1), p_print(3, 1), p_print(3, 2)
427 IF (mod(imap, 3) /= 0) iparticle2 = iparticle2 + 1
429 IF (mod(iseq, 3) /= 0) ip2 = ip2 + 1
430 tmp = hessian(iseq, icoord + j) - rep_env%f(imap, j)
431 tmp = -tmp/(2.0_dp*dx*mass(ip1)*mass(ip2))*1e6_dp
433 hessian(iseq, icoord + j) = tmp
442 particles(imap)%r(1:3) = pos0((i - 1)*3 + 1:(i - 1)*3 + 3)
447 rep_env%r(imap, j) = pos0(i)
452 minimum_energy = rep_env%f(rep_env%ndim + 1, j)
453 IF (output_unit > 0)
THEN
454 WRITE (output_unit,
'(T2,A)') &
455 "VIB| ",
" Minimum Structure - Energy and Forces:"
456 WRITE (output_unit,
'(T2,A,T43,A,T57,F24.12)') &
457 "VIB|",
"Total energy:", rep_env%f(rep_env%ndim + 1, j)
458 WRITE (output_unit,
'(T2,"VIB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
461 WRITE (output_unit,
'(T2,"VIB|",T12,A,T30,3(2X,F15.9))') &
462 particles(imap)%atomic_kind%name, &
463 rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, j)
468 IF (output_unit > 0)
THEN
469 WRITE (output_unit,
'(T2,A)')
"VIB| Hessian in cartesian coordinates"
474 CALL write_va_hessian(vib_section, para_env, ncoord, globenv, hessian, logger)
480 hessian(j, i) = hessian(i, j)
486 file_position=
"REWIND", extension=
".rrm")
487 IF (print_grrm > 0)
THEN
490 particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
492 ALLOCATE (hint1(ncoord, ncoord), rmass(ncoord))
495 rmass(3*(imap - 1) + 1:3*(imap - 1) + 3) = mass(imap)
499 hint1(j, i) = hessian(j, i)*rmass(i)*rmass(j)*1.0e-6_dp
502 nfrozen =
SIZE(particles) - natoms
503 CALL write_grrm(print_grrm, f_env%force_env, particles, minimum_energy, &
504 hessian=hint1, fixed_atoms=nfrozen)
505 DEALLOCATE (hint1, rmass)
511 file_position=
"REWIND", extension=
".scine")
512 IF (print_scine > 0)
THEN
515 particles(imap)%f(1:3) = rep_env%f((imap - 1)*3 + 1:(imap - 1)*3 + 3, 1)
517 nfrozen =
SIZE(particles) - natoms
518 cpassert(nfrozen == 0)
519 CALL write_scine(print_scine, f_env%force_env, particles, minimum_energy, hessian=hessian)
525 extension=
".eig", file_status=
"REPLACE", &
526 file_action=
"WRITE", do_backup=.true., &
527 file_form=
"UNFORMATTED")
528 IF (print_namd > 0)
THEN
532 ALLOCATE (dfull(ncoord, ncoord))
533 ALLOCATE (hint2dfull(
SIZE(dfull, 2),
SIZE(dfull, 2)))
534 ALLOCATE (heigvaldfull(
SIZE(dfull, 2)))
535 ALLOCATE (matm(ncoord, ncoord))
536 ALLOCATE (rmass(
SIZE(dfull, 2)))
539 CALL build_d_matrix(rottrm, nrottrm, dfull, full=.true., natoms=natoms)
542 hint2dfull(:, :) = matmul(transpose(dfull), matmul(hessian, dfull))
550 matm((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i)
554 dfull = matmul(matm, matmul(dfull, hint2dfull))
557 norm = 1.0_dp/sum(dfull(:, i)*dfull(:, i))
559 dfull(:, i) = sqrt(norm)*(dfull(:, i))
561 CALL write_eigs_unformatted(print_namd, ncoord, heigvaldfull, dfull)
562 DEALLOCATE (heigvaldfull)
563 DEALLOCATE (hint2dfull)
569 nvib = ncoord - nrottrm
570 ALLOCATE (h_eigval1(ncoord))
571 ALLOCATE (h_eigval2(
SIZE(d, 2)))
572 ALLOCATE (hint1(ncoord, ncoord))
573 ALLOCATE (hint2(
SIZE(d, 2),
SIZE(d, 2)))
574 ALLOCATE (rmass(
SIZE(d, 2)))
575 ALLOCATE (konst(
SIZE(d, 2)))
576 IF (calc_intens)
THEN
577 ALLOCATE (dip_deriv(3,
SIZE(d, 2)))
579 ALLOCATE (polar_deriv(3, 3,
SIZE(d, 2)))
582 ALLOCATE (intensities_d(
SIZE(d, 2)))
583 ALLOCATE (intensities_p(
SIZE(d, 2)))
584 ALLOCATE (depol_p(
SIZE(d, 2)))
585 ALLOCATE (depol_u(
SIZE(d, 2)))
586 intensities_d = 0._dp
587 intensities_p = 0._dp
590 hint1(:, :) = hessian
592 IF (output_unit > 0)
THEN
593 WRITE (output_unit,
'(T2,"VIB| Cartesian Low frequencies ---",4G12.5)') &
594 (h_eigval1(i), i=1, min(9, ncoord))
595 WRITE (output_unit,
'(T2,A)')
"VIB| Eigenvectors before removal of rotations and translations"
600 IF (output_unit_eig > 0)
THEN
601 CALL write_eigs_unformatted(output_unit_eig, ncoord, h_eigval1, hint1)
604 hint2(:, :) = matmul(transpose(d), matmul(hessian, d))
605 IF (calc_intens)
THEN
607 dip_deriv(i, :) = matmul(tmp_dip(:, i, 1), d)
611 polar_deriv(i, j, :) = matmul(tmp_polar(:, i, j, 1), d)
616 IF (output_unit > 0)
THEN
617 WRITE (output_unit,
'(T2,"VIB| Frequencies after removal of the rotations and translations")')
619 WRITE (output_unit,
'(T2,"VIB| Internal Low frequencies ---",4G12.5)') h_eigval2
624 hessian((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i)
628 d = matmul(hessian, matmul(d, hint2))
630 norm = 1.0_dp/sum(d(:, i)*d(:, i))
634 d(:, i) = sqrt(norm)*d(:, i)
636 konst(i) = sign(1.0_dp, h_eigval2(i))*2.0_dp*
pi**2*(abs(h_eigval2(i))/
massunit)**2*rmass(i)
637 IF (calc_intens)
THEN
640 d_deriv(:) = d_deriv(:) + dip_deriv(:, j)*hint2(j, i)
642 intensities_d(i) = sqrt(dot_product(d_deriv, d_deriv))
646 p_deriv(:, :) = p_deriv(:, :) + polar_deriv(:, :, j)*hint2(j, i)
650 p_deriv(:, :) = p_deriv(:, :)*conver
652 a1 = (p_deriv(1, 1) + p_deriv(2, 2) + p_deriv(3, 3))/3.0_dp
653 a2 = (p_deriv(1, 1) - p_deriv(2, 2))**2 + &
654 (p_deriv(2, 2) - p_deriv(3, 3))**2 + &
655 (p_deriv(3, 3) - p_deriv(1, 1))**2
656 a3 = (p_deriv(1, 2)**2 + p_deriv(2, 3)**2 + p_deriv(3, 1)**2)
657 intensities_p(i) = 45.0_dp*a1*a1 + 7.0_dp/2.0_dp*(a2 + 6.0_dp*a3)
659 dummy = 45.0_dp*a1*a1 + 4.0_dp/2.0_dp*(a2 + 6.0_dp*a3)
660 IF (dummy > 5.e-7_dp)
THEN
662 depol_p(i) = 3.0_dp/2.0_dp*(a2 + 6.0_dp*a3)/(45.0_dp*a1*a1 + &
663 4.0_dp/2.0_dp*(a2 + 6.0_dp*a3))
665 depol_u(i) = 6.0_dp/2.0_dp*(a2 + 6.0_dp*a3)/(45.0_dp*a1*a1 + &
666 7.0_dp/2.0_dp*(a2 + 6.0_dp*a3))
673 h_eigval2(i) = sign(1.0_dp, h_eigval2(i))*sqrt(abs(h_eigval2(i))*
massunit)*
vibfac/1000.0_dp
675 IF (calc_intens)
THEN
677 IF (.NOT. intens_ir)
THEN
678 WRITE (iounit,
'(T2,"VIB| No IR intensities available. Check input")')
680 IF (.NOT. intens_raman)
THEN
681 WRITE (iounit,
'(T2,"VIB| No Raman intensities available. Check input")')
688 NULLIFY (din, pin, depp, depu)
689 IF (intens_ir) din => intensities_d
690 IF (intens_raman) pin => intensities_p
691 IF (intens_raman) depp => depol_p
692 IF (intens_raman) depu => depol_u
693 CALL vib_out(iw, nvib, d, konst, rmass, h_eigval2, particles, mlist, din, pin, depp, depu)
695 IF (.NOT. something_frozen .AND. calc_thchdata)
THEN
696 CALL get_thch_values(h_eigval2, iw, mass, nvib, inertia, 1, minimum_energy, tc_temp, tc_press)
699 dump_only_positive=.false., logger=logger,
list=mlist)
701 IF (output_unit > 0)
THEN
702 WRITE (output_unit,
'(T2,"VIB| No further vibrational info. Detected a single atom")')
709 DEALLOCATE (h_eigval1)
710 DEALLOCATE (h_eigval2)
719 IF (calc_intens)
THEN
720 DEALLOCATE (dip_deriv)
721 DEALLOCATE (polar_deriv)
723 DEALLOCATE (tmp_polar)
725 DEALLOCATE (intensities_d)
726 DEALLOCATE (intensities_p)
735 CALL timestop(handle)
744 SUBROUTINE get_moving_atoms(force_env, Ilist)
745 TYPE(force_env_type),
POINTER :: force_env
746 INTEGER,
DIMENSION(:),
POINTER :: ilist
748 CHARACTER(len=*),
PARAMETER :: routinen =
'get_moving_atoms'
750 INTEGER :: handle, i, ii, ikind, j, ndim, &
751 nfixed_atoms, nfixed_atoms_total, nkind
752 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ifixd_list, work
753 TYPE(cp_subsys_type),
POINTER :: subsys
754 TYPE(fixd_constraint_type),
DIMENSION(:),
POINTER :: fixd_list
755 TYPE(molecule_kind_list_type),
POINTER :: molecule_kinds
756 TYPE(molecule_kind_type),
DIMENSION(:),
POINTER :: molecule_kind_set
757 TYPE(molecule_kind_type),
POINTER :: molecule_kind
758 TYPE(particle_list_type),
POINTER :: particles
759 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
761 CALL timeset(routinen, handle)
765 molecule_kinds=molecule_kinds)
767 nkind = molecule_kinds%n_els
768 molecule_kind_set => molecule_kinds%els
769 particle_set => particles%els
772 nfixed_atoms_total = 0
774 molecule_kind => molecule_kind_set(ikind)
776 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
778 ndim =
SIZE(particle_set) - nfixed_atoms_total
780 ALLOCATE (ilist(ndim))
782 IF (nfixed_atoms_total /= 0)
THEN
783 ALLOCATE (ifixd_list(nfixed_atoms_total))
784 ALLOCATE (work(nfixed_atoms_total))
785 nfixed_atoms_total = 0
787 molecule_kind => molecule_kind_set(ikind)
789 IF (
ASSOCIATED(fixd_list))
THEN
790 DO ii = 1,
SIZE(fixd_list)
791 IF (.NOT. fixd_list(ii)%restraint%active)
THEN
792 nfixed_atoms_total = nfixed_atoms_total + 1
793 ifixd_list(nfixed_atoms_total) = fixd_list(ii)%fixd
798 CALL sort(ifixd_list, nfixed_atoms_total, work)
802 loop_count:
DO i = 1,
SIZE(particle_set)
803 DO WHILE (i > ifixd_list(j))
805 IF (j > nfixed_atoms_total)
EXIT loop_count
807 IF (i /= ifixd_list(j))
THEN
812 DEALLOCATE (ifixd_list)
818 DO j = i,
SIZE(particle_set)
822 CALL timestop(handle)
824 END SUBROUTINE get_moving_atoms
842 SUBROUTINE vib_out(iw, nvib, D, k, m, freq, particles, Mlist, intensities_d, intensities_p, &
844 INTEGER,
INTENT(IN) :: iw, nvib
845 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: d
846 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: k, m, freq
847 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
848 INTEGER,
DIMENSION(:),
POINTER :: mlist
849 REAL(kind=
dp),
DIMENSION(:),
POINTER :: intensities_d, intensities_p, depol_p, &
852 CHARACTER(LEN=2) :: element_symbol
853 INTEGER :: from, iatom, icol, j, jatom, katom, &
855 REAL(kind=
dp) :: fint, pint
860 WRITE (unit=iw, fmt=
"(/,T2,'VIB|',T30,'NORMAL MODES - CARTESIAN DISPLACEMENTS')")
861 WRITE (unit=iw, fmt=
"(T2,'VIB|')")
862 DO jatom = 1, nvib, 3
864 to = min(from + 2, nvib)
865 WRITE (unit=iw, fmt=
"(T2,'VIB|',13X,3(8X,I5,8X))") &
866 (icol, icol=from, to)
867 WRITE (unit=iw, fmt=
"(T2,'VIB|Frequency (cm^-1)',3(1X,F12.6,8X))") &
868 (freq(icol), icol=from, to)
869 IF (
ASSOCIATED(intensities_d))
THEN
870 WRITE (unit=iw, fmt=
"(T2,'VIB|IR int (KM/Mole) ',3(1X,F12.6,8X))") &
871 (fint*intensities_d(icol)**2, icol=from, to)
873 IF (
ASSOCIATED(intensities_p))
THEN
874 WRITE (unit=iw, fmt=
"(T2,'VIB|Raman (A^4/amu) ',3(1X,F12.6,8X))") &
875 (pint*intensities_p(icol), icol=from, to)
876 WRITE (unit=iw, fmt=
"(T2,'VIB|Depol Ratio (P) ',3(1X,F12.6,8X))") &
877 (depol_p(icol), icol=from, to)
878 WRITE (unit=iw, fmt=
"(T2,'VIB|Depol Ratio (U) ',3(1X,F12.6,8X))") &
879 (depol_u(icol), icol=from, to)
881 WRITE (unit=iw, fmt=
"(T2,'VIB|Red.Masses (a.u.)',3(1X,F12.6,8X))") &
882 (m(icol), icol=from, to)
883 WRITE (unit=iw, fmt=
"(T2,'VIB|Frc consts (a.u.)',3(1X,F12.6,8X))") &
884 (k(icol), icol=from, to)
885 WRITE (unit=iw, fmt=
"(T2,' ATOM',2X,'EL',7X,3(4X,' X ',1X,' Y ',1X,' Z '))")
886 DO iatom = 1, natom, 3
888 IF (mod(iatom, 3) /= 0) katom = katom + 1
889 CALL get_atomic_kind(atomic_kind=particles(mlist(katom))%atomic_kind, &
890 element_symbol=element_symbol)
891 WRITE (unit=iw, fmt=
"(T2,I5,2X,A2,7X,3(4X,2(F5.2,1X),F5.2))") &
892 mlist(katom), element_symbol, &
893 ((d(iatom + j, icol), j=0, 2), icol=from, to)
895 WRITE (unit=iw, fmt=
"(/)")
898 END SUBROUTINE vib_out
910 SUBROUTINE build_d_matrix(mat, dof, Dout, full, natoms)
911 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: mat
912 INTEGER,
INTENT(IN) :: dof
913 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: dout
914 LOGICAL,
OPTIONAL :: full
915 INTEGER,
INTENT(IN) :: natoms
917 CHARACTER(len=*),
PARAMETER :: routinen =
'build_D_matrix'
919 INTEGER :: handle, i, ifound, iseq, j, nvib
921 REAL(kind=
dp) :: norm
922 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: work
923 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: d
925 CALL timeset(routinen, handle)
927 IF (
PRESENT(full)) my_full = full
929 nvib = 3*natoms - dof
930 ALLOCATE (work(3*natoms))
931 ALLOCATE (d(3*natoms, 3*natoms))
936 norm = dot_product(mat(:, i), mat(:, j))
938 cpwarn(
"Orthogonality error in transformation matrix")
945 DO WHILE (ifound /= nvib)
947 cpassert(iseq <= 3*natoms)
951 DO i = 1, dof + ifound
952 norm = dot_product(work, d(:, i))
953 work(:) = work - norm*d(:, i)
956 norm = sqrt(dot_product(work, work))
960 d(:, dof + ifound) = work/norm
963 cpassert(dof + ifound == 3*natoms)
967 dout = d(:, dof + 1:)
971 CALL timestop(handle)
972 END SUBROUTINE build_d_matrix
989 SUBROUTINE get_thch_values(freqs, iw, mass, nvib, inertia, spin, totene, temp, pressure)
991 REAL(kind=
dp),
DIMENSION(:) :: freqs
992 INTEGER,
INTENT(IN) :: iw
993 REAL(kind=
dp),
DIMENSION(:) :: mass
994 INTEGER,
INTENT(IN) :: nvib
995 REAL(kind=
dp),
INTENT(IN) :: inertia(3)
996 INTEGER,
INTENT(IN) :: spin
997 REAL(kind=
dp),
INTENT(IN) :: totene, temp, pressure
999 INTEGER :: i, natoms, sym_num
1000 REAL(kind=
dp) :: el_entropy, entropy, exp_min_one, fact, fact2, freq_arg, freq_arg2, &
1001 freqsum, gibbs, heat_capacity, inertia_kg(3), mass_tot, one_min_exp, partition_function, &
1002 rot_cv, rot_energy, rot_entropy, rot_part_func, rotvibtra, tran_cv, tran_energy, &
1003 tran_enthalpy, tran_entropy, tran_part_func, vib_cv, vib_energy, vib_entropy, &
1005 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: mass_kg
1012 freqsum = freqsum + freqs(i)
1021 ALLOCATE (mass_kg(natoms))
1022 mass_kg(:) = mass(:)**2*
e_mass
1023 mass_tot = sum(mass_kg)
1029 IF (inertia_kg(1)*inertia_kg(2)*inertia_kg(3) > 1.0_dp)
THEN
1030 rot_part_func = fact*fact*fact*inertia_kg(1)*inertia_kg(2)*inertia_kg(3)*
pi
1031 rot_part_func = sqrt(rot_part_func)
1037 IF (inertia_kg(1) > 1.0_dp)
THEN
1038 rot_part_func = fact*inertia_kg(1)
1039 ELSE IF (inertia_kg(2) > 1.0_dp)
THEN
1040 rot_part_func = fact*inertia_kg(2)
1042 rot_part_func = fact*inertia_kg(3)
1050 tran_part_func = (
boltzmann*temp)**2.5_dp/(pressure*(
h_bar*2.0_dp*
pi)**3.0_dp)*(2.0_dp*
pi*mass_tot)**1.5_dp
1057 vib_part_func = 1.0_dp
1059 vib_entropy = 0.0_dp
1064 freq_arg = fact*freqs(i)
1065 freq_arg2 = fact2*freqs(i)
1066 exp_min_one = exp(freq_arg) - 1.0_dp
1067 one_min_exp = 1.0_dp - exp(-freq_arg)
1071 vib_part_func = vib_part_func*(1.0_dp/one_min_exp)
1073 vib_energy = vib_energy + freq_arg2*0.5_dp + freq_arg2/exp_min_one
1075 vib_entropy = vib_entropy + freq_arg/exp_min_one - log(one_min_exp)
1077 vib_cv = vib_cv + freq_arg*freq_arg*exp(freq_arg)/exp_min_one/exp_min_one
1086 partition_function = rot_part_func*tran_part_func*vib_part_func
1090 entropy = el_entropy + rot_entropy + tran_entropy + vib_entropy
1094 rotvibtra = rot_energy + tran_enthalpy + vib_energy
1097 heat_capacity = vib_cv + tran_cv + rot_cv
1100 gibbs = vib_energy + rot_energy + tran_enthalpy - temp*entropy
1102 DEALLOCATE (mass_kg)
1105 WRITE (unit=iw, fmt=
"(/,T2,'VIB|',T30,'NORMAL MODES - THERMOCHEMICAL DATA')")
1106 WRITE (unit=iw, fmt=
"(T2,'VIB|')")
1108 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Symmetry number:',T70,I16)") sym_num
1109 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Temperature [K]:',T70,F16.2)") temp
1110 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Pressure [Pa]:',T70,F16.2)") pressure
1112 WRITE (unit=iw, fmt=
"(/)")
1114 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Electronic energy (U) [kJ/mol]:',T60,F26.8)") totene*
kjmol
1115 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Zero-point correction [kJ/mol]:',T60,F26.8)") zpe/1000.0_dp
1116 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Entropy [kJ/(mol K)]:',T60,F26.8)") entropy/1000.0_dp
1117 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Enthalpy correction (H-U) [kJ/mol]:',T60,F26.8)") rotvibtra/1000.0_dp
1118 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Gibbs energy correction [kJ/mol]:',T60,F26.8)") gibbs/1000.0_dp
1119 WRITE (unit=iw, fmt=
"(T2,'VIB|', T20, 'Heat capacity [kJ/(mol*K)]:',T70,F16.8)") heat_capacity/1000.0_dp
1120 WRITE (unit=iw, fmt=
"(/)")
1123 END SUBROUTINE get_thch_values
1134 SUBROUTINE write_eigs_unformatted(unit, dof, eigenvalues, eigenvectors)
1135 INTEGER,
INTENT(IN) :: unit, dof
1136 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenvalues
1137 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: eigenvectors
1139 CHARACTER(len=*),
PARAMETER :: routinen =
'write_eigs_unformatted'
1141 INTEGER :: handle, jj
1143 CALL timeset(routinen, handle)
1144 IF (unit .GT. 0)
THEN
1148 WRITE (unit) eigenvalues(1:dof)
1151 WRITE (unit) eigenvectors(1:dof, jj)
1154 CALL timestop(handle)
1156 END SUBROUTINE write_eigs_unformatted
1167 SUBROUTINE write_va_hessian(vib_section, para_env, ncoord, globenv, Hessian, logger)
1169 TYPE(section_vals_type),
POINTER :: vib_section
1170 TYPE(mp_para_env_type),
POINTER :: para_env
1172 TYPE(global_environment_type),
POINTER :: globenv
1173 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: hessian
1174 TYPE(cp_logger_type),
POINTER :: logger
1176 CHARACTER(LEN=*),
PARAMETER :: routinen =
'write_va_hessian'
1178 INTEGER :: handle, hesunit, i, j, ndf
1179 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
1180 TYPE(cp_fm_struct_type),
POINTER :: fm_struct_hes
1181 TYPE(cp_fm_type) :: hess_mat
1183 CALL timeset(routinen, handle)
1186 extension=
".hess", file_form=
"UNFORMATTED", file_action=
"WRITE", &
1187 file_position=
"REWIND")
1191 globenv%blacs_repeatable)
1194 nrow_global=ndf, ncol_global=ndf)
1195 CALL cp_fm_create(hess_mat, fm_struct_hes, name=
"hess_mat")
1208 CALL cp_fm_release(hess_mat)
1211 CALL timestop(handle)
1213 END SUBROUTINE write_va_hessian
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_write_unformatted(fm, unit)
...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
set of type/routines to handle the storage of results in force_envs
logical function, public test_for_result(results, description)
test for a certain result in the result_list
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
interface to use cp2k as library
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Define type storing the global information of a run. Keep the amount of stored data small....
subroutine, public write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, fixed_atoms)
Write GRRM interface file.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Collection of simple mathematical functions and subroutines.
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...
Interface to the message passing library MPI.
Module performing a mdoe selective vibrational analysis.
subroutine, public ms_vb_anal(input, rep_env, para_env, globenv, particles, nrep, calc_intens, dx, output_unit, logger)
Module performing a vibrational analysis.
Functions handling the MOLDEN format. Split from mode_selective.
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
represent a simple array based list of the given type
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Output Utilities for MOTION_SECTION.
real(kind=dp), parameter, public thrs_motion
subroutine, public rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weighted, natoms, rot_dof, inertia)
Performs an analysis of the principal inertia axis Getting back the generators of the translating and...
represent a simple array based list of the given type
Define methods related to particle_type.
subroutine, public write_particle_matrix(matrix, particle_set, iw, el_per_part, Ilist, parts_per_line)
...
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public boltzmann
real(kind=dp), parameter, public a_bohr
real(kind=dp), parameter, public vibfac
real(kind=dp), parameter, public n_avogadro
real(kind=dp), parameter, public kelvin
real(kind=dp), parameter, public h_bar
real(kind=dp), parameter, public hertz
real(kind=dp), parameter, public angstrom
real(kind=dp), parameter, public e_mass
real(kind=dp), parameter, public wavenumbers
real(kind=dp), parameter, public massunit
real(kind=dp), parameter, public kjmol
real(kind=dp), parameter, public pascal
real(kind=dp), parameter, public bohr
real(kind=dp), parameter, public debye
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
subroutine, public rep_env_create(rep_env, para_env, input, input_declaration, nrep, prep, sync_v, keep_wf_history, row_force)
creates a replica environment together with its force environment
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
types used to handle many replica of the same system that differ only in atom positions,...
subroutine, public rep_env_release(rep_env)
releases the given replica environment
subroutine, public write_scine(iounit, force_env, particles, energy, hessian)
Write SCINE interface file.
All kind of helpful little routines.
Module performing a vibrational analysis.
subroutine, public vb_anal(input, input_declaration, para_env, globenv)
Module performing a vibrational analysis.