xTB vibrational analysis: enable intensity calculation, 1 new regtest (…

…#313)
cp2k · Apr 18, 2019 · c9da12c · c9da12c
1 parent 5fed092
commit c9da12c
Show file tree

Hide file tree

Showing 3 changed files with 63 additions and 2 deletions.
diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F
@@ -41,6 +41,9 @@ MODULE qs_scf_post_tb
                                               cp_print_key_unit_nr
    USE cp_para_types,                   ONLY: cp_para_env_type
    USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
+   USE cp_result_methods,               ONLY: cp_results_erase,&
+                                              put_results
+   USE cp_result_types,                 ONLY: cp_result_type
    USE dbcsr_api,                       ONLY: dbcsr_p_type,&
                                               dbcsr_type
    USE input_constants,                 ONLY: ot_precond_full_all
@@ -615,11 +618,12 @@ SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
       TYPE(cell_type), POINTER                           :: cell
+      TYPE(cp_result_type), POINTER                      :: results
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
-      NULLIFY (atomic_kind_set, cell)
+      NULLIFY (atomic_kind_set, cell, results)
       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
-                      particle_set=particle_set, cell=cell)
+                      particle_set=particle_set, cell=cell, results=results)
 
       ! Reference point
       reference = section_get_ival(input, keyword_name="REFERENCE")
@@ -685,6 +689,9 @@ SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
             dipole_deriv(:) = dipole_deriv(:)-q*(particle_set(i)%v(:)-drcc)
          END DO
       END IF
+      CALL cp_results_erase(results=results, description=description)
+      CALL put_results(results=results, description=description, &
+                       values=dipole(1:3))
       IF (unit_nr > 0) THEN
          WRITE (unit_nr, *)
          WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//"(A.U.)|", dipole

diff --git a/tests/xTB/regtest-4/TEST_FILES b/tests/xTB/regtest-4/TEST_FILES
@@ -15,4 +15,5 @@ H2O-xtb-tc2-1.inp                         11    1.0E-12             -92.31598419
 H2O-xtb-tc2-2.inp                         11    1.0E-12             -92.315984190681831
 H2O-32-xtb-trs4.inp                       11    1.0E-12             -46.150675543227536
 H2O-32-xtb-trs4-dyn.inp                   11    1.0E-12             -46.150675543227536
+h2o_vib.inp                                8    1.0E-09                     1467.812892          
 #EOF
diff --git a/tests/xTB/regtest-4/h2o_vib.inp b/tests/xTB/regtest-4/h2o_vib.inp
@@ -0,0 +1,53 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &DFT
+    &QS
+      METHOD xTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      &OT
+         PRECONDITIONER FULL_SINGLE_INVERSE
+         MINIMIZER DIIS
+      &END
+      &OUTER_SCF
+         MAX_SCF  10
+      &END
+      MAX_SCF  10
+    &END SCF
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 6 6 6
+      PERIODIC NONE
+    &END
+    &COORD
+    O   0.000000    0.000000   -0.065587
+    H   0.000000   -0.757136    0.520545
+    H   0.000000    0.757136    0.520545
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT h2oVib
+  RUN_TYPE VIBRATIONAL_ANALYSIS
+&END GLOBAL
+
+&VIBRATIONAL_ANALYSIS
+ INTENSITIES .TRUE.
+ NPROC_REP 1
+ DX 0.001
+&END
+
+