NEWTONX interface - writing CP2K output

src/input_cp2k_properties_dft.F

@@ -1559,6 +1559,17 @@ SUBROUTINE create_tddfpt2_section(section)
       CALL section_add_subsection(subsection, print_key)
       CALL section_release(print_key)
 
+      CALL cp_print_key_section_create(print_key, __LOCATION__, name="NAMD_PRINT", &
+                                       description="Controls the printout required for NAMD with NEWTONX.", &
+                                       print_level=debug_print_level + 1, filename="CP2K_NEWTONX")
+      CALL keyword_create(keyword, __LOCATION__, name="PRINT_VIRTUALS", &
+                          description="Print occupied AND virtual molecular orbital coefficients", &
+                          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)
+
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 

src/qs_tddfpt2_methods.F

@@ -74,6 +74,7 @@ MODULE qs_tddfpt2_methods
                                               tddfpt_print_nto_analysis,&
                                               tddfpt_print_summary
    USE qs_tddfpt2_restart,              ONLY: tddfpt_read_restart,&
+                                              tddfpt_write_newtonx_output,&
                                               tddfpt_write_restart
    USE qs_tddfpt2_stda_types,           ONLY: allocate_stda_env,&
                                               deallocate_stda_env,&
@@ -404,6 +405,8 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
       ! *** print summary information ***
       log_unit = cp_logger_get_default_io_unit()
 
+      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 "// &

src/qs_tddfpt2_restart.F

@@ -16,6 +16,8 @@ MODULE qs_tddfpt2_restart
    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_log_handling,                 ONLY: cp_logger_type
    USE cp_output_handling,              ONLY: cp_p_file,&
@@ -30,6 +32,8 @@ MODULE qs_tddfpt2_restart
    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 qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
    USE string_utilities,                ONLY: integer_to_string
 #include "./base/base_uses.f90"
@@ -45,7 +49,7 @@ MODULE qs_tddfpt2_restart
    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
 
-   PUBLIC :: tddfpt_write_restart, tddfpt_read_restart
+   PUBLIC :: tddfpt_write_restart, tddfpt_read_restart, tddfpt_write_newtonx_output
 
 ! **************************************************************************************************
 
@@ -287,5 +291,124 @@ FUNCTION tddfpt_read_restart(evects, evals, gs_mos, logger, tddfpt_section, tddf
       CALL timestop(handle)
 
    END FUNCTION tddfpt_read_restart
+! **************************************************************************************************
+!> \brief Write Ritz vectors to a binary restart file.
+!> \param evects               vectors to store
+!> \param evals                TDDFPT eigenvalues
+!> \param gs_mos               structure that holds ground state occupied and virtual
+!>                             molecular orbitals
+!> \param logger               a logger object
+!> \param tddfpt_print_section TDDFPT%PRINT input section
+!> \param qs_env ...
+! **************************************************************************************************
+   SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_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
+
+      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
+
+      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)
+
+         nspins = SIZE(evects, 1)
+         nstates = SIZE(evects, 2)
+
+         IF (debug_this_module) THEN
+            CPASSERT(SIZE(evals) == nstates)
+            CPASSERT(nspins > 0)
+            CPASSERT(nstates > 0)
+         END IF
+
+         CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
+
+         ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "NAMD_PRINT", &
+                                      extension=".inp", file_form="FORMATTED", file_action="WRITE", file_status="REPLACE")
+
+         ! print eigenvectors
+         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)
+               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)
+            END DO
+         END DO
+
+         ! print molecular orbitals
+
+         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"
+               CALL cp_fm_write_info(gs_mos(ispin)%mos_occ, ounit)
+            END IF
+            CALL cp_fm_write_formatted(gs_mos(ispin)%mos_occ, ounit, "OCCUPIED MO COEFFICIENTS")
+         END DO
+
+         IF (ounit > 0) THEN
+            WRITE (ounit, "(A)") "OCCUPIED MO EIGENVALUES"
+            DO ispin = 1, nspins
+               nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
+               DO iocc = 1, nmo_occ(ispin)
+                  WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_occ(iocc)
+               END DO
+            END DO
+         END IF
+!
+         IF (print_virtuals) THEN
+            DO ispin = 1, nspins
+               IF (ounit > 0) THEN
+                  WRITE (ounit, "(/,A)") "VIRTUAL MOS SIZE"
+                  CALL cp_fm_write_info(gs_mos(ispin)%mos_virt, ounit)
+               END IF
+               CALL cp_fm_write_formatted(gs_mos(ispin)%mos_virt, ounit, "VIRTUAL MO COEFFICIENTS")
+            END DO
+
+            IF (ounit > 0) THEN
+               WRITE (ounit, "(A)") "VIRTUAL MO EIGENVALUES"
+               DO ispin = 1, nspins
+                  nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
+                  DO ivirt = 1, nmo_virt(ispin)
+                     WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_virt(ivirt)
+                  END DO
+               END DO
+            END IF
+         END IF
+
+         CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "NAMD_PRINT")
+
+         CALL timestop(handle)
+      END IF
+
+   END SUBROUTINE tddfpt_write_newtonx_output
+! **************************************************************************************************
 
 END MODULE qs_tddfpt2_restart

src/qs_tddfpt2_utils.F

@@ -53,7 +53,8 @@ MODULE qs_tddfpt2_utils
                                               section_vals_release,&
                                               section_vals_retain,&
                                               section_vals_set_subs_vals,&
-                                              section_vals_type
+                                              section_vals_type,&
+                                              section_vals_val_get
    USE kinds,                           ONLY: dp,&
                                               int_8
    USE message_passing,                 ONLY: mp_min,&
@@ -109,22 +110,28 @@ SUBROUTINE tddfpt_init_mos(qs_env, gs_mos)
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_init_mos'
 
-      INTEGER                                            :: handle, ispin, nmo_avail, nmo_occ, &
-                                                            nmo_virt, nspins
+      INTEGER                                            :: handle, iounit, ispin, nmo_avail, &
+                                                            nmo_occ, nmo_virt, nspins
+      LOGICAL                                            :: print_virtuals_newtonx
       REAL(kind=dp), DIMENSION(:), POINTER               :: evals_virt_spin
       TYPE(cell_type), POINTER                           :: cell
       TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER        :: evals_virt
       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
       TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: mos_virt
       TYPE(cp_fm_type), POINTER                          :: mos_virt_spin
+      TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
       TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
       TYPE(qs_scf_env_type), POINTER                     :: scf_env
+      TYPE(section_vals_type), POINTER                   :: print_section
       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
 
       CALL timeset(routineN, handle)
 
+      logger => cp_get_default_logger()
+      iounit = cp_logger_get_default_io_unit(logger)
+
       CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
                       matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
       tddfpt_control => dft_control%tddfpt2_control
@@ -134,12 +141,17 @@ SUBROUTINE tddfpt_init_mos(qs_env, gs_mos)
       nspins = dft_control%nspins
       ALLOCATE (gs_mos(nspins))
 
+      ! check if virtuals should be constructed for NAMD interface with NEWTONX
+      print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
+      CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
+
       ! when the number of unoccupied orbitals is limited and OT has been used
       ! for the ground-state DFT calculation,
       ! compute the missing unoccupied orbitals using OT as well.
       NULLIFY (evals_virt, evals_virt_spin, mos_virt, mos_virt_spin)
       IF (ASSOCIATED(scf_env)) THEN
-         IF (scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) THEN
+         IF ((scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
+             (scf_env%method == ot_method_nr .AND. print_virtuals_newtonx)) THEN
             ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to:
             ! nmo_virt = tddfpt_control%nlumo
             ! number of already computed unoccupied orbitals (added_mos) .
@@ -150,10 +162,12 @@ SUBROUTINE tddfpt_init_mos(qs_env, gs_mos)
             END DO
             ! number of unoccupied orbitals to compute
             nmo_virt = tddfpt_control%nlumo - nmo_virt
-            IF (nmo_virt > 0) THEN
-               ALLOCATE (evals_virt(nspins), mos_virt(nspins))
-               ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
-               CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
+            IF (.NOT. print_virtuals_newtonx) THEN
+               IF (nmo_virt > 0) THEN
+                  ALLOCATE (evals_virt(nspins), mos_virt(nspins))
+                  ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
+                  CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
+               END IF
             END IF
          END IF
       END IF

tests/QS/regtest-tddfpt-force/TEST_FILES

@@ -15,4 +15,5 @@ h2o_f08.inp                                                1    1.0E-11
 h2o_f09.inp                                                1    1.0E-11             -17.21975280565217
 h2o_f10.inp                                                1    1.0E-11             -17.02954820220558
 h2o_f10b.inp                                               1    1.0E-11             -17.08503047402411
+h2o_tddfpt_newtonx.inp                                     1    1.0E-11             -17.22181320160399 
 #EOF

tests/QS/regtest-tddfpt-force/h2o_tddfpt_newtonx.inp

@@ -0,0 +1,92 @@
+&GLOBAL
+  PROJECT test
+  RUN_TYPE ENERGY_FORCE
+  PREFERRED_DIAG_LIBRARY SL
+  PRINT_LEVEL medium
+&END GLOBAL
+&FORCE_EVAL
+ &PRINT 
+ &FORCES
+ &END FORCES
+ &END PRINT 
+  METHOD Quickstep
+&PROPERTIES
+  &TDDFPT
+   KERNEL FULL
+   NSTATES 1
+   MAX_ITER   100
+   CONVERGENCE [eV] 1.0e-7
+   RKS_TRIPLETS F
+   &PRINT
+    &NAMD_PRINT
+     PRINT_VIRTUALS T
+    &END NAMD_PRINT
+   &END PRINT
+  &END TDDFPT
+&END PROPERTIES
+  &DFT
+    &QS
+      METHOD GPW
+    &END QS
+ &EXCITED_STATES T
+  STATE 1
+ &END EXCITED_STATES
+    &SCF
+      SCF_GUESS restart
+      &OT
+         PRECONDITIONER FULL_ALL
+         MINIMIZER DIIS
+      &END OT
+      &OUTER_SCF
+         MAX_SCF 900
+         EPS_SCF 1.0E-8
+      &END OUTER_SCF
+      MAX_SCF 10
+      EPS_SCF 1.0E-8
+    &END SCF
+    POTENTIAL_FILE_NAME POTENTIAL_UZH
+    BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
+    BASIS_SET_FILE_NAME BASIS_ccGRB_UZH
+    BASIS_SET_FILE_NAME BASIS_ADMM_UZH
+    &MGRID
+      CUTOFF 600
+      REL_CUTOFF  60
+    &END MGRID
+    &POISSON
+       PERIODIC NONE
+       POISSON_SOLVER WAVELET 
+    &END
+    &XC
+     &XC_FUNCTIONAL PBE
+     &END XC_FUNCTIONAL
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 6.0 6.0 6.0
+      PERIODIC NONE
+    &END CELL
+     &COORD
+       UNIT angstrom
+       O   0.000000    0.000000    0.000000
+       H   0.000000   -0.757136    0.580545
+       H   0.000000    0.757136    0.580545
+     &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES T
+     &END
+     NATOMS 3
+     CONNECTIVITY OFF
+    &END TOPOLOGY
+    &KIND H 
+      BASIS_SET AUX_FIT admm-dzp-q1
+      BASIS_SET ORB DZVP-MOLOPT-PBE0-GTH-q1
+      POTENTIAL GTH-PBE0-q1
+    &END KIND
+    &KIND O
+      BASIS_SET AUX_FIT admm-dzp-q6
+      BASIS_SET ORB DZVP-MOLOPT-PBE0-GTH-q6
+      POTENTIAL GTH-PBE0-q6
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL