comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Further changes for NewtonX interface

Browse source
This commit is contained in:
Anna Hehn 2022-06-21 13:30:01 +02:00 • committed by Ole Schütt
parent 1ed20a10a6
commit 7444654cf1
11 changed files with 515 additions and 101 deletions
Download patch file Download diff file Expand all files Collapse all files

4
data/EMSL_BASIS_SETS
View file

@ -2858,8 +2858,8 @@ S 6-311Gxx 6-311G**
902.45600000 0.11902800 0.00000000
297.15800000 0.36843200 0.00000000
108.70200000 0.57729900 0.14318600
43.15530000 0.00000000 0.62446500
18.10790000 0.00000000 0.28336600
43.155300000 0.00000000 0.62446500
18.107900000 0.00000000 0.28336600
1 0 0 1 1
5.56009000 1.00000000
1 0 0 1 1

8
src/input_cp2k_exstate.F
View file

@ -74,6 +74,14 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OVERLAP_DELTAT", &
description="Keyword for the computation of the overlap matrix between two consecutive time steps.", &
usage="OVERLAP_DELTAT", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_exstate_section
END MODULE input_cp2k_exstate

10
src/input_cp2k_properties_dft.F
View file

@ -1600,6 +1600,16 @@ CONTAINS
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PRINT_PHASES", &
description="Print phases of occupied and virtuals MOs.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SCALE_WITH_PHASES", &
description="Scale ES eigenvectors with phases of occupied and virtuals MOs.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection, print_key)
CALL section_release(print_key)

19
src/motion/input_cp2k_vib.F
View file

@ -134,12 +134,6 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NAMD_PRINT", &
description="Adjust cartesian eigenvalues / vectors to NewtonX format.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL create_mode_selective_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -217,6 +211,19 @@ CONTAINS
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key, __LOCATION__, name="NAMD_PRINT", &
description="Adjust cartesian eigenvalues / vectors to NewtonX format.", &
print_level=debug_print_level + 1, add_last=add_last_numeric, &
filename="FullNormalizedCartesian")
CALL keyword_create(keyword, __LOCATION__, name="BACKUP_COPIES", &
description="Specifies the maximum index of backup copies.", &
usage="BACKUP_COPIES {int}", &
default_i_val=3)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key, __LOCATION__, "HESSIAN", &
description="Write the Hessian matrix from a vibrational analysis calculation "// &
"into a binary file.", &

82
src/motion/vibrational_analysis.F
View file

@ -110,10 +110,10 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: description_d, description_p
INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, iounit, ip1, ip2, iparticle1, &
iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nRotTrM, nvib, &
output_unit, output_unit_eig, prep, print_grrm, print_scine, proc_dist_type
output_unit, output_unit_eig, prep, print_grrm, print_namd, print_scine, proc_dist_type
INTEGER, DIMENSION(:), POINTER :: Clist, Mlist
LOGICAL :: calc_intens, calc_thchdata, do_mode_tracking, intens_ir, intens_raman, &
keep_rotations, lprint_namd, row_force, something_frozen
keep_rotations, row_force, something_frozen
REAL(KIND=dp) :: a1, a2, a3, conver, dummy, dx, &
inertia(3), minimum_energy, norm, &
tc_press, tc_temp, tmp
@ -167,7 +167,7 @@ CONTAINS
CALL section_vals_val_get(vib_section, "THERMOCHEMISTRY", l_val=calc_thchdata)
CALL section_vals_val_get(vib_section, "TC_TEMPERATURE", r_val=tc_temp)
CALL section_vals_val_get(vib_section, "TC_PRESSURE", r_val=tc_press)
CALL section_vals_val_get(vib_section, "NAMD_PRINT", l_val=lprint_namd)
tc_temp = tc_temp*kelvin
tc_press = tc_press*pascal
@ -520,6 +520,52 @@ CONTAINS
END IF
CALL cp_print_key_finished_output(print_scine, logger, force_env_section, "PRINT%SCINE")
!
! Print NEWTONX interface file
print_namd = cp_print_key_unit_nr(logger, print_section, "NAMD_PRINT", &
extension=".eig", file_status="REPLACE", &
file_action="WRITE", do_backup=.TRUE., &
file_form="UNFORMATTED")
IF (print_namd > 0) THEN
! NewtonX requires normalized Cartesian frequencies and eigenvectors
! in full matrix format (ncoord x ncoord)
NULLIFY (Dfull)
ALLOCATE (Dfull(ncoord, ncoord))
ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2)))
ALLOCATE (HeigvalDfull(SIZE(Dfull, 2)))
ALLOCATE (MatM(ncoord, ncoord))
ALLOCATE (rmass(SIZE(Dfull, 2)))
Dfull = 0.0_dp
! Dfull in dimension of degrees of freedom
CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms)
! TEST MatM = MATMUL(TRANSPOSE(Dfull),Dfull)= 1
! Hessian in MWC -> Hessian in INT (Hint2Dfull)
Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull))
! Heig = L^T Hint2Dfull L
CALL diamat_all(Hint2Dfull, HeigvalDfull)
! TEST MatM = MATMUL(TRANSPOSE(Hint2Dfull),Hint2Dfull) = 1
! TEST MatM=MATMUL(TRANSPOSE(MATMUL(Dfull,Hint2Dfull)),MATMUL(Dfull,Hint2Dfull)) = 1
MatM = 0.0_dp
DO i = 1, natoms
DO j = 1, 3
MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i) ! mass is sqrt(mass)
END DO
END DO
! Dfull = Cartesian displacements of the normal modes
Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !Dfull=D L / sqrt(m)
DO i = 1, ncoord
! Renormalize displacements
norm = 1.0_dp/SUM(Dfull(:, i)*Dfull(:, i))
rmass(i) = norm/massunit
Dfull(:, i) = SQRT(norm)*(Dfull(:, i))
END DO
CALL write_eigs_unformatted(print_namd, ncoord, HeigvalDfull, Dfull)
DEALLOCATE (HeigvalDfull)
DEALLOCATE (Hint2Dfull)
DEALLOCATE (Dfull)
DEALLOCATE (MatM)
DEALLOCATE (rmass)
END IF !print_namd
!
nvib = ncoord - nRotTrM
ALLOCATE (H_eigval1(ncoord))
ALLOCATE (H_eigval2(SIZE(D, 2)))
@ -552,35 +598,7 @@ CONTAINS
END IF
! write frequencies and eigenvectors to cartesian eig file
IF (output_unit_eig > 0) THEN
IF (lprint_namd) THEN
! NewtonX requires normalized Cartesian frequencies and eigenvectors
! in full matrix format (ncoord x ncoord)
NULLIFY (Dfull)
ALLOCATE (Dfull(ncoord, ncoord))
ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2)))
ALLOCATE (HeigvalDfull(SIZE(Dfull, 2)))
ALLOCATE (MatM(ncoord, ncoord))
Dfull = 0.0_dp
CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms)
Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull)) ! MWC -> INT
CALL diamat_all(Hint2Dfull, HeigvalDfull)
MatM = 0.0_dp
DO i = 1, natoms
DO j = 1, 3
MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i)
END DO
END DO
! Cartesian displacements of the normal modes
Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !!! this is !D L / sqrt(m)
! CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Dfull)
CALL write_eigs_unformatted(output_unit_eig, ncoord, HeigvalDfull, Dfull)
DEALLOCATE (HeigvalDfull)
DEALLOCATE (Hint2Dfull)
DEALLOCATE (Dfull)
DEALLOCATE (MatM)
END IF
CALL write_eigs_unformatted(output_unit_eig, &
ncoord, H_eigval1, Hint1)
CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Hint1)
END IF
IF (nvib /= 0) THEN
Hint2(:, :) = MATMUL(TRANSPOSE(D), MATMUL(Hessian, D))

40
src/qs_energy.F
View file

@ -18,6 +18,10 @@ MODULE qs_energy
USE dm_ls_scf, ONLY: ls_scf
USE energy_corrections, ONLY: energy_correction
USE excited_states, ONLY: excited_state_energy
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE lri_environment_methods, ONLY: lri_print_stat
USE mp2, ONLY: mp2_main
USE qs_energy_init, ONLY: qs_energies_init
@ -63,10 +67,12 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies'
INTEGER :: handle
LOGICAL :: do_consistent_energies, my_calc_forces, &
run_rtp
LOGICAL :: do_consistent_energies, &
do_excited_state, loverlap_deltat, &
my_calc_forces, run_rtp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_energy_type), POINTER :: energy
TYPE(section_vals_type), POINTER :: excited_state_section
CALL timeset(routineN, handle)
@ -85,13 +91,25 @@ CONTAINS
CALL qs_energies_init(qs_env, my_calc_forces)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, energy=energy)
! *** check if only overlap matrix is needed for couplings
loverlap_deltat = .FALSE.
NULLIFY (excited_state_section)
excited_state_section => section_vals_get_subs_vals(qs_env%input, "DFT%EXCITED_STATES")
CALL section_vals_get(excited_state_section, explicit=do_excited_state)
IF (do_excited_state) THEN
CALL section_vals_val_get(excited_state_section, "OVERLAP_DELTAT", &
l_val=loverlap_deltat)
END IF
! *** Perform a SCF run ***
IF (dft_control%qs_control%do_ls_scf) THEN
CALL ls_scf(qs_env=qs_env)
ELSE IF (dft_control%qs_control%do_almo_scf) THEN
CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces)
ELSE
CALL scf(qs_env=qs_env)
IF (.NOT. loverlap_deltat) THEN
IF (dft_control%qs_control%do_ls_scf) THEN
CALL ls_scf(qs_env=qs_env)
ELSE IF (dft_control%qs_control%do_almo_scf) THEN
CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces)
ELSE
CALL scf(qs_env=qs_env)
END IF
END IF
IF (do_consistent_energies) THEN
@ -114,9 +132,11 @@ CONTAINS
CALL energy_correction(qs_env, ec_init=.TRUE., calculate_forces=.FALSE.)
END IF
CALL qs_energies_properties(qs_env, calc_forces)
IF (.NOT. loverlap_deltat) THEN
CALL qs_energies_properties(qs_env, calc_forces)
CALL excited_state_energy(qs_env, calculate_forces=.FALSE.)
CALL excited_state_energy(qs_env, calculate_forces=.FALSE.)
END IF
IF (dft_control%qs_control%lrigpw) THEN
CALL lri_print_stat(qs_env)

16
src/qs_tddfpt2_methods.F
View file

@ -143,7 +143,7 @@ CONTAINS
old_state
INTEGER, DIMENSION(maxspins) :: nactive
LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, &
is_restarted, state_change
explicit, is_restarted, state_change
REAL(kind=dp) :: conv
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
TYPE(admm_type), POINTER :: admm_env
@ -159,8 +159,9 @@ CONTAINS
TYPE(kernel_env_type) :: kernel_env
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(section_vals_type), POINTER :: lri_section, tddfpt_print_section, &
tddfpt_section, xc_section
TYPE(section_vals_type), POINTER :: lri_section, namd_print_section, &
tddfpt_print_section, tddfpt_section, &
xc_section
TYPE(stda_env_type), TARGET :: stda_kernel
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
@ -433,8 +434,6 @@ CONTAINS
CALL cp_warn(__LOCATION__, "Skipping TDDFPT wavefunction optimization")
END IF
CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_env)
IF (ASSOCIATED(matrix_ks_oep) .AND. tddfpt_control%dipole_form == tddfpt_dipole_velocity) THEN
CALL cp_warn(__LOCATION__, &
"Transition dipole moments and oscillator strengths are likely to be incorrect "// &
@ -444,6 +443,13 @@ CONTAINS
! *** print summary information ***
log_unit = cp_logger_get_default_io_unit()
namd_print_section => section_vals_get_subs_vals(tddfpt_print_section, "NAMD_PRINT")
CALL section_vals_get(namd_print_section, explicit=explicit)
IF (explicit) THEN
CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, &
matrix_s(1)%matrix, S_evects, sub_env) !qs_env, sub_env)
END IF
ALLOCATE (ostrength(nstates))
ostrength = 0.0_dp
CALL tddfpt_print_summary(log_unit, evects, evals, ostrength, mult, &

334
src/qs_tddfpt2_restart.F
View file

@ -8,17 +8,20 @@
MODULE qs_tddfpt2_restart
USE cp_blacs_env, ONLY: cp_blacs_env_type,&
get_blacs_info
USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
USE cp_files, ONLY: close_file,&
open_file
USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
cp_fm_scale_and_add,&
cp_fm_trace
USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
fm_pool_create_fm
USE cp_fm_types, ONLY: cp_fm_get_info,&
cp_fm_p_type,&
cp_fm_read_unformatted,&
cp_fm_write_formatted,&
cp_fm_write_info,&
cp_fm_write_unformatted
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: &
cp_fm_create, cp_fm_get_info, cp_fm_p_type, cp_fm_read_unformatted, cp_fm_release, &
cp_fm_type, cp_fm_write_formatted, cp_fm_write_info, cp_fm_write_unformatted
USE cp_log_handling, ONLY: cp_logger_type
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
@ -26,14 +29,17 @@ MODULE qs_tddfpt2_restart
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
dp
USE message_passing, ONLY: mp_bcast
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE message_passing, ONLY: mp_bcast,&
mp_min,&
mp_sum
USE parallel_gemm_api, ONLY: parallel_gemm
USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
USE string_utilities, ONLY: integer_to_string
#include "./base/base_uses.f90"
@ -300,30 +306,39 @@ CONTAINS
!> molecular orbitals
!> \param logger a logger object
!> \param tddfpt_print_section TDDFPT%PRINT input section
!> \param qs_env ...
!> \param matrix_s ...
!> \param S_evects ...
!> \param sub_env ...
! **************************************************************************************************
SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_env)
SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, &
matrix_s, S_evects, sub_env)
TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects
REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: tddfpt_print_section
TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
TYPE(cp_logger_type), INTENT(in), POINTER :: logger
TYPE(section_vals_type), INTENT(in), POINTER :: tddfpt_print_section
TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects
TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_newtonx_output'
INTEGER :: handle, iocc, ispin, istate, ivirt, nao, &
nspins, nstates, ounit
INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
LOGICAL :: print_virtuals
TYPE(cp_para_env_type), POINTER :: para_env
LOGICAL :: print_phases, print_virtuals, &
scale_with_phases
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase_evects
TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: evects_mo
TYPE(cp_fm_struct_type), POINTER :: fmstruct
IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "NAMD_PRINT"), cp_p_file)) THEN
CALL timeset(routineN, handle)
NULLIFY (para_env)
CALL get_qs_env(qs_env, para_env=para_env)
CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals)
CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%SCALE_WITH_PHASES", l_val=scale_with_phases)
nspins = SIZE(evects, 1)
nstates = SIZE(evects, 2)
@ -336,36 +351,101 @@ CONTAINS
CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
IF (sub_env%is_split) THEN
CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when states"// &
" are distributed to different CPU pools.")
END IF
ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "NAMD_PRINT", &
extension=".inp", file_form="FORMATTED", file_action="WRITE", file_status="REPLACE")
IF (debug_this_module) CALL tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
! print eigenvectors
IF (print_virtuals) THEN
ALLOCATE (evects_mo(nspins, nstates))
DO istate = 1, nstates
DO ispin = 1, nspins
! transform eigenvectors
NULLIFY (fmstruct)
nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
CALL cp_fm_struct_create(fmstruct, para_env=sub_env%para_env, &
context=sub_env%blacs_env, &
nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin))
NULLIFY (evects_mo(ispin, istate)%matrix)
ALLOCATE (evects_mo(ispin, istate)%matrix)
CALL cp_fm_create(evects_mo(ispin, istate)%matrix, fmstruct)
CALL cp_fm_struct_release(fmstruct)
CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, istate)%matrix, S_evects(ispin, istate)%matrix, &
ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
END DO
END DO
DO istate = 1, nstates
DO ispin = 1, nspins
CALL parallel_gemm("T", "N", &
nmo_virt(ispin), &
nmo_occ(ispin), &
nao, &
1.0_dp, &
gs_mos(ispin)%mos_virt, &
S_evects(ispin, istate)%matrix, & !this also needs to be orthogonalized
0.0_dp, &
evects_mo(ispin, istate)%matrix)
END DO
END DO
END IF
DO istate = 1, nstates
DO ispin = 1, nspins
! TDDFPT wave function is actually stored as a linear combination of virtual MOs
! that replaces the corresponding deoccupied MO. Unfortunately, the phase
! of the occupied MOs varies depending on the eigensolver used as well as
! how eigenvectors are distributed across computational cores. The phase is important
! because TDDFPT wave functions are used to compute a response electron density
! \rho^{-} = 1/2 * [C_{0} * evect^T + evect * C_{0}^{-}], where C_{0} is the expansion
! coefficients of the reference ground-state wave function. To make the restart file
! transferable, TDDFPT wave functions are stored in assumption that all ground state
! MOs have a positive phase.
CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ)
IF (ounit > 0) THEN
WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
CALL cp_fm_write_info(evects(ispin, istate)%matrix, ounit)
IF (.NOT. print_virtuals) THEN
CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ)
IF (ounit > 0) THEN
WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
CALL cp_fm_write_info(evects(ispin, istate)%matrix, ounit)
END IF
CALL cp_fm_write_formatted(evects(ispin, istate)%matrix, ounit, "ES EIGENVECTORS")
ELSE
CALL cp_fm_column_scale(evects_mo(ispin, istate)%matrix, gs_mos(ispin)%phases_occ)
IF (ounit > 0) THEN
WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
CALL cp_fm_write_info(evects_mo(ispin, istate)%matrix, ounit)
END IF
CALL cp_fm_write_formatted(evects_mo(ispin, istate)%matrix, ounit, "ES EIGENVECTORS")
END IF
CALL cp_fm_write_formatted(evects(ispin, istate)%matrix, ounit, "ES EIGENVECTORS")
CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ)
! compute and print phase of eigenvectors
nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
ALLOCATE (phase_evects(nmo_occ(ispin)))
IF (print_virtuals) THEN
CALL compute_phase_eigenvectors(evects_mo(ispin, istate)%matrix, phase_evects, sub_env)
ELSE
CALL compute_phase_eigenvectors(evects(ispin, istate)%matrix, phase_evects, sub_env)
END IF
IF (ounit > 0) THEN
WRITE (ounit, "(/,A,/)") "PHASES ES EIGENVECTORS"
DO iocc = 1, nmo_occ(ispin)
WRITE (ounit, "(F20.14)") phase_evects(iocc)
END DO
END IF
DEALLOCATE (phase_evects)
END DO
END DO
! print molecular orbitals
IF (print_virtuals) THEN
IF (ASSOCIATED(evects_mo)) THEN
DO istate = 1, nstates
DO ispin = 1, nspins
CALL cp_fm_release(matrix=evects_mo(ispin, istate)%matrix)
DEALLOCATE (evects_mo(ispin, istate)%matrix)
END DO
END DO
END IF
DEALLOCATE (evects_mo)
END IF
CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals)
DO ispin = 1, nspins
IF (ounit > 0) THEN
WRITE (ounit, "(/,A)") "OCCUPIED MOS SIZE"
@ -404,6 +484,27 @@ CONTAINS
END IF
END IF
! print phases of molecular orbitals
IF (print_phases) THEN
IF (ounit > 0) THEN
WRITE (ounit, "(A)") "PHASES OCCUPIED ORBITALS"
DO ispin = 1, nspins
DO iocc = 1, nmo_occ(ispin)
WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_occ(iocc)
END DO
END DO
IF (print_virtuals) THEN
WRITE (ounit, "(A)") "PHASES VIRTUAL ORBITALS"
DO ispin = 1, nspins
DO ivirt = 1, nmo_virt(ispin)
WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_virt(ivirt)
END DO
END DO
END IF
END IF
END IF
CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "NAMD_PRINT")
CALL timestop(handle)
@ -411,5 +512,166 @@ CONTAINS
END SUBROUTINE tddfpt_write_newtonx_output
! **************************************************************************************************
!> \brief ...
!> \param evects ...
!> \param ounit ...
!> \param S_evects ...
!> \param matrix_s ...
! **************************************************************************************************
SUBROUTINE tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects
INTEGER, INTENT(in) :: ounit
TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects
TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_check_orthonormality'
INTEGER :: handle, ispin, ivect, jvect, nspins, &
nvects_total
INTEGER, DIMENSION(maxspins) :: nactive
REAL(kind=dp) :: norm
REAL(kind=dp), DIMENSION(maxspins) :: weights
CALL timeset(routineN, handle)
nspins = SIZE(evects, 1)
nvects_total = SIZE(evects, 2)
IF (debug_this_module) THEN
CPASSERT(SIZE(S_evects, 1) == nspins)
CPASSERT(SIZE(S_evects, 2) == nvects_total)
END IF
DO ispin = 1, nspins
CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, ncol_global=nactive(ispin))
END DO
DO jvect = 1, nvects_total
! <psi1_i | psi1_j>
DO ivect = 1, jvect - 1
CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
norm = SUM(weights(1:nspins))
DO ispin = 1, nspins
CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect)%matrix, -norm, evects(ispin, ivect)%matrix)
END DO
END DO
! <psi1_j | psi1_j>
DO ispin = 1, nspins
CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect)%matrix, S_evects(ispin, jvect)%matrix, &
ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
END DO
CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
norm = SUM(weights(1:nspins))
norm = 1.0_dp/SQRT(norm)
IF ((ounit > 0) .AND. debug_this_module) WRITE (ounit, '(A,F10.8)') "norm", norm
END DO
CALL timestop(handle)
END SUBROUTINE tddfpt_check_orthonormality
! **************************************************************************************************
!> \brief ...
!> \param evects ...
!> \param phase_evects ...
!> \param sub_env ...
! **************************************************************************************************
SUBROUTINE compute_phase_eigenvectors(evects, phase_evects, sub_env)
! copied from parts of tddgpt_init_ground_state_mos by S. Chulkov
TYPE(cp_fm_type), INTENT(in) :: evects
REAL(kind=dp), DIMENSION(:), INTENT(out) :: phase_evects
TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
CHARACTER(len=*), PARAMETER :: routineN = 'compute_phase_eigenvectors'
REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp)
INTEGER :: handle, icol_global, icol_local, irow_global, irow_local, ncol_global, &
ncol_local, nrow_global, nrow_local, sign_int
INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
sum_sign_array
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
REAL(kind=dp) :: element
REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
POINTER :: my_block
CALL timeset(routineN, handle)
! compute and print the phase of excited-state eigenvectors:
CALL cp_fm_get_info(evects, nrow_global=nrow_global, ncol_global=ncol_global, &
nrow_local=nrow_local, ncol_local=ncol_local, local_data=my_block, &
row_indices=row_indices, col_indices=col_indices) ! nrow_global either nao or nocc
ALLOCATE (minrow_neg_array(ncol_global), minrow_pos_array(ncol_global), sum_sign_array(ncol_global))
minrow_neg_array(:) = nrow_global
minrow_pos_array(:) = nrow_global
sum_sign_array(:) = 0
DO icol_local = 1, ncol_local
icol_global = col_indices(icol_local)
DO irow_local = 1, nrow_local
irow_global = row_indices(irow_local)
element = my_block(irow_local, icol_local)
sign_int = 0
IF (element >= eps_dp) THEN
sign_int = 1
ELSE IF (element <= -eps_dp) THEN
sign_int = -1
END IF
sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
IF (sign_int > 0) THEN
IF (minrow_pos_array(icol_global) > irow_global) &
minrow_pos_array(icol_global) = irow_global
ELSE IF (sign_int < 0) THEN
IF (minrow_neg_array(icol_global) > irow_global) &
minrow_neg_array(icol_global) = irow_global
END IF
END DO
END DO
CALL mp_sum(sum_sign_array, sub_env%para_env%group)
CALL mp_min(minrow_neg_array, sub_env%para_env%group)
CALL mp_min(minrow_pos_array, sub_env%para_env%group)
DO icol_global = 1, ncol_global
IF (sum_sign_array(icol_global) > 0) THEN
! most of the expansion coefficients are positive => MO's phase = +1
phase_evects(icol_global) = 1.0_dp
ELSE IF (sum_sign_array(icol_global) < 0) THEN
! most of the expansion coefficients are negative => MO's phase = -1
phase_evects(icol_global) = -1.0_dp
ELSE
! equal number of positive and negative expansion coefficients
IF (minrow_pos_array(icol_global) <= minrow_neg_array(icol_global)) THEN
! the first positive expansion coefficient has a lower index then
! the first negative expansion coefficient; MO's phase = +1
phase_evects(icol_global) = 1.0_dp
ELSE
! MO's phase = -1
phase_evects(icol_global) = -1.0_dp
END IF
END IF
END DO
DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
CALL timestop(handle)
END SUBROUTINE compute_phase_eigenvectors
END MODULE qs_tddfpt2_restart

1
src/qs_tddfpt2_types.F
View file

@ -106,6 +106,7 @@ MODULE qs_tddfpt2_types
!> phase of occupied MOs; +1.0 -- positive, -1.0 -- negative;
!> it is mainly needed to make the restart file transferable
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phases_occ
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phases_virt
END TYPE tddfpt_ground_state_mos
! **************************************************************************************************

97
src/qs_tddfpt2_utils.F
View file

@ -182,7 +182,8 @@ CONTAINS
nlumo=tddfpt_control%nlumo, &
blacs_env=blacs_env, cholesky_method=cholesky_restore, &
matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
mos_virt=mos_virt_spin, evals_virt=evals_virt_spin)
mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
qs_env=qs_env)
END DO
moc = 0
@ -236,13 +237,14 @@ CONTAINS
!> (in addition to the ADDED_MOS, if present). NULL when no OT is in use.
!> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals.
!> NULL when no OT is in use.
!> \param qs_env ...
!> \par History
!> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov]
!> * 06.2016 renamed, altered prototype [Sergey Chulkov]
!> * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov]
! **************************************************************************************************
SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, &
mos_virt, evals_virt)
mos_virt, evals_virt, qs_env)
TYPE(tddfpt_ground_state_mos) :: gs_mos
TYPE(mo_set_type), INTENT(IN) :: mo_set
INTEGER, INTENT(in) :: nlumo
@ -251,16 +253,17 @@ CONTAINS
TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s
TYPE(cp_fm_type), INTENT(IN), POINTER :: mos_virt
REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt
TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_ground_state_mos'
REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp)
INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, irow_global, &
INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
sum_sign_array
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
LOGICAL :: do_eigen
LOGICAL :: do_eigen, print_phases
REAL(kind=dp) :: element, maxocc
REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
POINTER :: my_block
@ -268,18 +271,28 @@ CONTAINS
mo_occ_scf
TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, &
ao_mo_virt_fm_struct, wfn_fm_struct
TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm
TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
work_fm_virt
TYPE(cp_fm_type), POINTER :: mo_coeff_extended
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(mo_set_type), POINTER :: mos_extended
TYPE(section_vals_type), POINTER :: print_section
CALL timeset(routineN, handle)
NULLIFY (logger)
logger => cp_get_default_logger()
iounit = cp_logger_get_default_io_unit(logger)
CALL get_blacs_info(blacs_env, para_env=para_env)
CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
nelectron=nelectrons, occupation_numbers=mo_occ_scf)
print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
nmo_virt = nao - nmo_occ
IF (nlumo >= 0) &
nmo_virt = MIN(nmo_virt, nlumo)
@ -291,6 +304,7 @@ CONTAINS
do_eigen = .FALSE.
! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
IF (ASSOCIATED(evals_virt)) THEN
CPASSERT(ASSOCIATED(mos_virt))
IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE.
ELSE
IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE.
@ -309,8 +323,7 @@ CONTAINS
ALLOCATE (gs_mos%evals_occ(nmo_occ))
ALLOCATE (gs_mos%evals_virt(nmo_virt))
ALLOCATE (gs_mos%phases_occ(nmo_occ))
CALL cp_fm_struct_release(ao_mo_virt_fm_struct)
ALLOCATE (gs_mos%phases_virt(nmo_virt))
! ++ nullify pointers
NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
@ -457,6 +470,73 @@ CONTAINS
gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
END IF
IF (print_phases) THEN
! compute the phase of molecular orbitals;
! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors
!CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
row_indices=row_indices, col_indices=col_indices, local_data=my_block)
ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
minrow_neg_array(:) = nao
minrow_pos_array(:) = nao
sum_sign_array(:) = 0
DO icol_local = 1, ncol_local
icol_global = col_indices(icol_local)
DO irow_local = 1, nrow_local
element = my_block(irow_local, icol_local)
sign_int = 0
IF (element >= eps_dp) THEN
sign_int = 1
ELSE IF (element <= -eps_dp) THEN
sign_int = -1
END IF
sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
irow_global = row_indices(irow_local)
IF (sign_int > 0) THEN
IF (minrow_pos_array(icol_global) > irow_global) &
minrow_pos_array(icol_global) = irow_global
ELSE IF (sign_int < 0) THEN
IF (minrow_neg_array(icol_global) > irow_global) &
minrow_neg_array(icol_global) = irow_global
END IF
END DO
END DO
CALL mp_sum(sum_sign_array, para_env%group)
CALL mp_min(minrow_neg_array, para_env%group)
CALL mp_min(minrow_pos_array, para_env%group)
DO icol_local = 1, nmo_virt
IF (sum_sign_array(icol_local) > 0) THEN
! most of the expansion coefficients are positive => MO's phase = +1
gs_mos%phases_virt(icol_local) = 1.0_dp
ELSE IF (sum_sign_array(icol_local) < 0) THEN
! most of the expansion coefficients are negative => MO's phase = -1
gs_mos%phases_virt(icol_local) = -1.0_dp
ELSE
! equal number of positive and negative expansion coefficients
IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
! the first positive expansion coefficient has a lower index then
! the first negative expansion coefficient; MO's phase = +1
gs_mos%phases_virt(icol_local) = 1.0_dp
ELSE
! MO's phase = -1
gs_mos%phases_virt(icol_local) = -1.0_dp
END IF
END IF
END DO
DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
CALL cp_fm_release(work_fm_virt)
END IF !print_phases
CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
CALL cp_fm_release(work_fm)
IF (do_eigen) THEN
@ -492,6 +572,9 @@ CONTAINS
IF (ALLOCATED(gs_mos%evals_occ)) &
DEALLOCATE (gs_mos%evals_occ)
IF (ALLOCATED(gs_mos%phases_virt)) &
DEALLOCATE (gs_mos%phases_virt)
IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
CALL cp_fm_release(gs_mos%evals_occ_matrix)
DEALLOCATE (gs_mos%evals_occ_matrix)

5
tests/QS/regtest-gpw-9/h2o_NEWTONX.inp
View file

@ -43,9 +43,8 @@
&END FORCE_EVAL
&VIBRATIONAL_ANALYSIS
&PRINT
&CARTESIAN_EIGS
&END CARTESIAN_EIGS
&NAMD_PRINT
&END NAMD_PRINT
&END PRINT
NAMD_PRINT
DX 0.001
&END VIBRATIONAL_ANALYSIS