Skip to content

Commit 7444654

Browse files
Anna Hehnoschuett
Anna Hehn
authored and
oschuett
committed
Further changes for NewtonX interface
1 parent 1ed20a1 commit 7444654

File tree

11 files changed

+514
-100
lines changed

11 files changed

+514
-100
lines changed

‎data/EMSL_BASIS_SETS‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2858,8 +2858,8 @@ S 6-311Gxx 6-311G**
28582858
902.45600000 0.11902800 0.00000000
28592859
297.15800000 0.36843200 0.00000000
28602860
108.70200000 0.57729900 0.14318600
2861-
43.15530000 0.00000000 0.62446500
2862-
18.10790000 0.00000000 0.28336600
2861+
43.155300000 0.00000000 0.62446500
2862+
18.107900000 0.00000000 0.28336600
28632863
1 0 0 1 1
28642864
5.56009000 1.00000000
28652865
1 0 0 1 1

‎src/input_cp2k_exstate.F‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,14 @@ SUBROUTINE create_exstate_section(section)
7474
CALL section_add_keyword(section, keyword)
7575
CALL keyword_release(keyword)
7676

77+
CALL keyword_create(keyword, __LOCATION__, name="OVERLAP_DELTAT", &
78+
description="Keyword for the computation of the overlap matrix between two consecutive time steps.", &
79+
usage="OVERLAP_DELTAT", &
80+
default_l_val=.FALSE., &
81+
lone_keyword_l_val=.TRUE.)
82+
CALL section_add_keyword(section, keyword)
83+
CALL keyword_release(keyword)
84+
7785
END SUBROUTINE create_exstate_section
7886

7987
END MODULE input_cp2k_exstate

‎src/input_cp2k_properties_dft.F‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1600,6 +1600,16 @@ SUBROUTINE create_tddfpt2_section(section)
16001600
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
16011601
CALL section_add_keyword(print_key, keyword)
16021602
CALL keyword_release(keyword)
1603+
CALL keyword_create(keyword, __LOCATION__, name="PRINT_PHASES", &
1604+
description="Print phases of occupied and virtuals MOs.", &
1605+
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
1606+
CALL section_add_keyword(print_key, keyword)
1607+
CALL keyword_release(keyword)
1608+
CALL keyword_create(keyword, __LOCATION__, name="SCALE_WITH_PHASES", &
1609+
description="Scale ES eigenvectors with phases of occupied and virtuals MOs.", &
1610+
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
1611+
CALL section_add_keyword(print_key, keyword)
1612+
CALL keyword_release(keyword)
16031613
CALL section_add_subsection(subsection, print_key)
16041614
CALL section_release(print_key)
16051615

‎src/motion/input_cp2k_vib.F‎

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -134,12 +134,6 @@ SUBROUTINE create_vib_section(section)
134134
CALL section_add_keyword(section, keyword)
135135
CALL keyword_release(keyword)
136136

137-
CALL keyword_create(keyword, __LOCATION__, name="NAMD_PRINT", &
138-
description="Adjust cartesian eigenvalues / vectors to NewtonX format.", &
139-
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
140-
CALL section_add_keyword(section, keyword)
141-
CALL keyword_release(keyword)
142-
143137
CALL create_mode_selective_section(subsection)
144138
CALL section_add_subsection(section, subsection)
145139
CALL section_release(subsection)
@@ -217,6 +211,19 @@ SUBROUTINE create_print_vib_section(section)
217211
CALL section_add_subsection(section, print_key)
218212
CALL section_release(print_key)
219213

214+
CALL cp_print_key_section_create(print_key, __LOCATION__, name="NAMD_PRINT", &
215+
description="Adjust cartesian eigenvalues / vectors to NewtonX format.", &
216+
print_level=debug_print_level + 1, add_last=add_last_numeric, &
217+
filename="FullNormalizedCartesian")
218+
CALL keyword_create(keyword, __LOCATION__, name="BACKUP_COPIES", &
219+
description="Specifies the maximum index of backup copies.", &
220+
usage="BACKUP_COPIES {int}", &
221+
default_i_val=3)
222+
CALL section_add_keyword(print_key, keyword)
223+
CALL keyword_release(keyword)
224+
CALL section_add_subsection(section, print_key)
225+
CALL section_release(print_key)
226+
220227
CALL cp_print_key_section_create(print_key, __LOCATION__, "HESSIAN", &
221228
description="Write the Hessian matrix from a vibrational analysis calculation "// &
222229
"into a binary file.", &

‎src/motion/vibrational_analysis.F‎

Lines changed: 50 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -110,10 +110,10 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
110110
CHARACTER(LEN=default_string_length) :: description_d, description_p
111111
INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, iounit, ip1, ip2, iparticle1, &
112112
iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nRotTrM, nvib, &
113-
output_unit, output_unit_eig, prep, print_grrm, print_scine, proc_dist_type
113+
output_unit, output_unit_eig, prep, print_grrm, print_namd, print_scine, proc_dist_type
114114
INTEGER, DIMENSION(:), POINTER :: Clist, Mlist
115115
LOGICAL :: calc_intens, calc_thchdata, do_mode_tracking, intens_ir, intens_raman, &
116-
keep_rotations, lprint_namd, row_force, something_frozen
116+
keep_rotations, row_force, something_frozen
117117
REAL(KIND=dp) :: a1, a2, a3, conver, dummy, dx, &
118118
inertia(3), minimum_energy, norm, &
119119
tc_press, tc_temp, tmp
@@ -167,7 +167,7 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
167167
CALL section_vals_val_get(vib_section, "THERMOCHEMISTRY", l_val=calc_thchdata)
168168
CALL section_vals_val_get(vib_section, "TC_TEMPERATURE", r_val=tc_temp)
169169
CALL section_vals_val_get(vib_section, "TC_PRESSURE", r_val=tc_press)
170-
CALL section_vals_val_get(vib_section, "NAMD_PRINT", l_val=lprint_namd)
170+
171171
tc_temp = tc_temp*kelvin
172172
tc_press = tc_press*pascal
173173
@@ -520,6 +520,52 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
520520
END IF
521521
CALL cp_print_key_finished_output(print_scine, logger, force_env_section, "PRINT%SCINE")
522522
!
523+
! Print NEWTONX interface file
524+
print_namd = cp_print_key_unit_nr(logger, print_section, "NAMD_PRINT", &
525+
extension=".eig", file_status="REPLACE", &
526+
file_action="WRITE", do_backup=.TRUE., &
527+
file_form="UNFORMATTED")
528+
IF (print_namd > 0) THEN
529+
! NewtonX requires normalized Cartesian frequencies and eigenvectors
530+
! in full matrix format (ncoord x ncoord)
531+
NULLIFY (Dfull)
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)))
537+
Dfull = 0.0_dp
538+
! Dfull in dimension of degrees of freedom
539+
CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms)
540+
! TEST MatM = MATMUL(TRANSPOSE(Dfull),Dfull)= 1
541+
! Hessian in MWC -> Hessian in INT (Hint2Dfull)
542+
Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull))
543+
! Heig = L^T Hint2Dfull L
544+
CALL diamat_all(Hint2Dfull, HeigvalDfull)
545+
! TEST MatM = MATMUL(TRANSPOSE(Hint2Dfull),Hint2Dfull) = 1
546+
! TEST MatM=MATMUL(TRANSPOSE(MATMUL(Dfull,Hint2Dfull)),MATMUL(Dfull,Hint2Dfull)) = 1
547+
MatM = 0.0_dp
548+
DO i = 1, natoms
549+
DO j = 1, 3
550+
MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i) ! mass is sqrt(mass)
551+
END DO
552+
END DO
553+
! Dfull = Cartesian displacements of the normal modes
554+
Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !Dfull=D L / sqrt(m)
555+
DO i = 1, ncoord
556+
! Renormalize displacements
557+
norm = 1.0_dp/SUM(Dfull(:, i)*Dfull(:, i))
558+
rmass(i) = norm/massunit
559+
Dfull(:, i) = SQRT(norm)*(Dfull(:, i))
560+
END DO
561+
CALL write_eigs_unformatted(print_namd, ncoord, HeigvalDfull, Dfull)
562+
DEALLOCATE (HeigvalDfull)
563+
DEALLOCATE (Hint2Dfull)
564+
DEALLOCATE (Dfull)
565+
DEALLOCATE (MatM)
566+
DEALLOCATE (rmass)
567+
END IF !print_namd
568+
!
523569
nvib = ncoord - nRotTrM
524570
ALLOCATE (H_eigval1(ncoord))
525571
ALLOCATE (H_eigval2(SIZE(D, 2)))
@@ -552,35 +598,7 @@ SUBROUTINE vb_anal(input, input_declaration, para_env, globenv)
552598
END IF
553599
! write frequencies and eigenvectors to cartesian eig file
554600
IF (output_unit_eig > 0) THEN
555-
IF (lprint_namd) THEN
556-
! NewtonX requires normalized Cartesian frequencies and eigenvectors
557-
! in full matrix format (ncoord x ncoord)
558-
NULLIFY (Dfull)
559-
ALLOCATE (Dfull(ncoord, ncoord))
560-
ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2)))
561-
ALLOCATE (HeigvalDfull(SIZE(Dfull, 2)))
562-
ALLOCATE (MatM(ncoord, ncoord))
563-
Dfull = 0.0_dp
564-
CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms)
565-
Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull)) ! MWC -> INT
566-
CALL diamat_all(Hint2Dfull, HeigvalDfull)
567-
MatM = 0.0_dp
568-
DO i = 1, natoms
569-
DO j = 1, 3
570-
MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i)
571-
END DO
572-
END DO
573-
! Cartesian displacements of the normal modes
574-
Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !!! this is !D L / sqrt(m)
575-
! CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Dfull)
576-
CALL write_eigs_unformatted(output_unit_eig, ncoord, HeigvalDfull, Dfull)
577-
DEALLOCATE (HeigvalDfull)
578-
DEALLOCATE (Hint2Dfull)
579-
DEALLOCATE (Dfull)
580-
DEALLOCATE (MatM)
581-
END IF
582-
CALL write_eigs_unformatted(output_unit_eig, &
583-
ncoord, H_eigval1, Hint1)
601+
CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Hint1)
584602
END IF
585603
IF (nvib /= 0) THEN
586604
Hint2(:, :) = MATMUL(TRANSPOSE(D), MATMUL(Hessian, D))

‎src/qs_energy.F‎

Lines changed: 30 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,10 @@ MODULE qs_energy
1818
USE dm_ls_scf, ONLY: ls_scf
1919
USE energy_corrections, ONLY: energy_correction
2020
USE excited_states, ONLY: excited_state_energy
21+
USE input_section_types, ONLY: section_vals_get,&
22+
section_vals_get_subs_vals,&
23+
section_vals_type,&
24+
section_vals_val_get
2125
USE lri_environment_methods, ONLY: lri_print_stat
2226
USE mp2, ONLY: mp2_main
2327
USE qs_energy_init, ONLY: qs_energies_init
@@ -63,10 +67,12 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces)
6367
CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies'
6468

6569
INTEGER :: handle
66-
LOGICAL :: do_consistent_energies, my_calc_forces, &
67-
run_rtp
70+
LOGICAL :: do_consistent_energies, &
71+
do_excited_state, loverlap_deltat, &
72+
my_calc_forces, run_rtp
6873
TYPE(dft_control_type), POINTER :: dft_control
6974
TYPE(qs_energy_type), POINTER :: energy
75+
TYPE(section_vals_type), POINTER :: excited_state_section
7076

7177
CALL timeset(routineN, handle)
7278

@@ -85,13 +91,25 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces)
8591
CALL qs_energies_init(qs_env, my_calc_forces)
8692
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, energy=energy)
8793

94+
! *** check if only overlap matrix is needed for couplings
95+
loverlap_deltat = .FALSE.
96+
NULLIFY (excited_state_section)
97+
excited_state_section => section_vals_get_subs_vals(qs_env%input, "DFT%EXCITED_STATES")
98+
CALL section_vals_get(excited_state_section, explicit=do_excited_state)
99+
IF (do_excited_state) THEN
100+
CALL section_vals_val_get(excited_state_section, "OVERLAP_DELTAT", &
101+
l_val=loverlap_deltat)
102+
END IF
103+
88104
! *** Perform a SCF run ***
89-
IF (dft_control%qs_control%do_ls_scf) THEN
90-
CALL ls_scf(qs_env=qs_env)
91-
ELSE IF (dft_control%qs_control%do_almo_scf) THEN
92-
CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces)
93-
ELSE
94-
CALL scf(qs_env=qs_env)
105+
IF (.NOT. loverlap_deltat) THEN
106+
IF (dft_control%qs_control%do_ls_scf) THEN
107+
CALL ls_scf(qs_env=qs_env)
108+
ELSE IF (dft_control%qs_control%do_almo_scf) THEN
109+
CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces)
110+
ELSE
111+
CALL scf(qs_env=qs_env)
112+
END IF
95113
END IF
96114

97115
IF (do_consistent_energies) THEN
@@ -114,9 +132,11 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces)
114132
CALL energy_correction(qs_env, ec_init=.TRUE., calculate_forces=.FALSE.)
115133
END IF
116134

117-
CALL qs_energies_properties(qs_env, calc_forces)
135+
IF (.NOT. loverlap_deltat) THEN
136+
CALL qs_energies_properties(qs_env, calc_forces)
118137

119-
CALL excited_state_energy(qs_env, calculate_forces=.FALSE.)
138+
CALL excited_state_energy(qs_env, calculate_forces=.FALSE.)
139+
END IF
120140

121141
IF (dft_control%qs_control%lrigpw) THEN
122142
CALL lri_print_stat(qs_env)

‎src/qs_tddfpt2_methods.F‎

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
143143
old_state
144144
INTEGER, DIMENSION(maxspins) :: nactive
145145
LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, &
146-
is_restarted, state_change
146+
explicit, is_restarted, state_change
147147
REAL(kind=dp) :: conv
148148
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
149149
TYPE(admm_type), POINTER :: admm_env
@@ -159,8 +159,9 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
159159
TYPE(kernel_env_type) :: kernel_env
160160
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
161161
TYPE(qs_scf_env_type), POINTER :: scf_env
162-
TYPE(section_vals_type), POINTER :: lri_section, tddfpt_print_section, &
163-
tddfpt_section, xc_section
162+
TYPE(section_vals_type), POINTER :: lri_section, namd_print_section, &
163+
tddfpt_print_section, tddfpt_section, &
164+
xc_section
164165
TYPE(stda_env_type), TARGET :: stda_kernel
165166
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
166167
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
@@ -433,8 +434,6 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
433434
CALL cp_warn(__LOCATION__, "Skipping TDDFPT wavefunction optimization")
434435
END IF
435436

436-
CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_env)
437-
438437
IF (ASSOCIATED(matrix_ks_oep) .AND. tddfpt_control%dipole_form == tddfpt_dipole_velocity) THEN
439438
CALL cp_warn(__LOCATION__, &
440439
"Transition dipole moments and oscillator strengths are likely to be incorrect "// &
@@ -444,6 +443,13 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
444443

445444
! *** print summary information ***
446445
log_unit = cp_logger_get_default_io_unit()
446+
447+
namd_print_section => section_vals_get_subs_vals(tddfpt_print_section, "NAMD_PRINT")
448+
CALL section_vals_get(namd_print_section, explicit=explicit)
449+
IF (explicit) THEN
450+
CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, &
451+
matrix_s(1)%matrix, S_evects, sub_env) !qs_env, sub_env)
452+
END IF
447453
ALLOCATE (ostrength(nstates))
448454
ostrength = 0.0_dp
449455
CALL tddfpt_print_summary(log_unit, evects, evals, ostrength, mult, &

0 commit comments

Comments
 (0)