diff --git a/src/qs_linres_op.F b/src/qs_linres_op.F
index bbe91d24a1..f94188ae9b 100644
--- a/src/qs_linres_op.F
+++ b/src/qs_linres_op.F
@@ -45,7 +45,9 @@ MODULE qs_linres_op
    USE cp_para_types,                   ONLY: cp_para_env_type
    USE dbcsr_api,                       ONLY: &
         dbcsr_checksum, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
-        dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_p_type, dbcsr_set, dbcsr_type, &
+        dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_get_block_p, &
+        dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
+        dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_set, dbcsr_type, &
         dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_type
@@ -636,14 +638,23 @@ CONTAINS
       TYPE(qs_environment_type), POINTER                 :: qs_env
 
       LOGICAL                                            :: do_periodic
+      TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(polar_env_type), POINTER                      :: polar_env
 
-      CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
+      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, polar_env=polar_env)
       CALL get_polar_env(polar_env=polar_env, do_periodic=do_periodic)
-      IF (do_periodic) THEN
-         CALL polar_operators_berry(qs_env)
+      IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
+         IF (do_periodic) THEN
+            CALL polar_tb_operators_berry(qs_env)
+         ELSE
+            CALL polar_tb_operators_local(qs_env)
+         END IF
       ELSE
-         CALL polar_operators_local(qs_env)
+         IF (do_periodic) THEN
+            CALL polar_operators_berry(qs_env)
+         ELSE
+            CALL polar_operators_local(qs_env)
+         END IF
       END IF
 
    END SUBROUTINE polar_operators
@@ -823,6 +834,108 @@ CONTAINS
 
    END SUBROUTINE polar_operators_berry
 
+! **************************************************************************************************
+!> \brief Calculate the Berry phase operator in the AO basis and
+!>         then the derivative of the Berry phase operator with respect to
+!>         the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
+!>        afterwards multiply with the ground state MO coefficients
+!>
+!> \param qs_env ...
+!> \par History
+!>      01.2013 created [SL]
+!>      06.2018 polar_env integrated into qs_env (MK)
+!>      08.2020 adapt for xTB/DFTB (JHU)
+!> \author SL
+! **************************************************************************************************
+
+   SUBROUTINE polar_tb_operators_berry(qs_env)
+
+      TYPE(qs_environment_type), POINTER                 :: qs_env
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_tb_operators_berry'
+
+      COMPLEX(dp)                                        :: zdeta
+      INTEGER                                            :: blk, handle, i, icol, idir, irow, ispin, &
+                                                            nmo, nspins
+      LOGICAL                                            :: do_raman, found
+      REAL(dp)                                           :: dd, fdir
+      REAL(dp), DIMENSION(3)                             :: kvec, ria, rib
+      REAL(dp), DIMENSION(3, 3)                          :: hmat
+      REAL(dp), DIMENSION(:, :), POINTER                 :: d_block, s_block
+      TYPE(cell_type), POINTER                           :: cell
+      TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER       :: dBerry_psi0
+      TYPE(cp_fm_type), POINTER                          :: mo_coeff
+      TYPE(dbcsr_iterator_type)                          :: iter
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dipmat, matrix_s
+      TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
+      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
+      TYPE(polar_env_type), POINTER                      :: polar_env
+
+      CALL timeset(routineN, handle)
+
+      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
+                      cell=cell, particle_set=particle_set, &
+                      polar_env=polar_env, mos=mos, matrix_s=matrix_s)
+
+      nspins = dft_control%nspins
+
+      CALL get_polar_env(polar_env=polar_env, &
+                         do_raman=do_raman, &
+                         dBerry_psi0=dBerry_psi0)
+
+      IF (do_raman) THEN
+
+         ALLOCATE (dipmat(3))
+         DO i = 1, 3
+            ALLOCATE (dipmat(i)%matrix)
+            CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
+            CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
+         END DO
+
+         hmat = cell%hmat(:, :)/twopi
+
+         CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
+         DO WHILE (dbcsr_iterator_blocks_left(iter))
+            NULLIFY (s_block, d_block)
+            CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
+            ria = particle_set(irow)%r
+            rib = particle_set(icol)%r
+            DO idir = 1, 3
+               kvec(:) = twopi*cell%h_inv(idir, :)
+               dd = SUM(kvec(:)*ria(:))
+               zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
+               fdir = AIMAG(LOG(zdeta))
+               dd = SUM(kvec(:)*rib(:))
+               zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
+               fdir = fdir + AIMAG(LOG(zdeta))
+               CALL dbcsr_get_block_p(matrix=dipmat(idir)%matrix, &
+                                      row=irow, col=icol, BLOCK=d_block, found=found)
+               CPASSERT(found)
+               d_block = d_block + 0.5_dp*fdir*s_block
+            END DO
+         ENDDO
+         CALL dbcsr_iterator_stop(iter)
+
+         ! Compute the derivative and add the result to mo_derivatives
+         DO ispin = 1, dft_control%nspins ! spin
+            CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo)
+            DO i = 1, 3
+               CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
+                                            dBerry_psi0(i, ispin)%matrix, ncol=nmo)
+            ENDDO !x/y/z-direction
+         END DO
+
+         DO i = 1, 3
+            CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
+         END DO
+         DEALLOCATE (dipmat)
+
+      ENDIF ! do_raman
+
+      CALL timestop(handle)
+   END SUBROUTINE polar_tb_operators_berry
+
 ! **************************************************************************************************
 !> \brief Calculate the Berry phase operator in the AO basis and
 !>         then the derivative of the Berry phase operator with respect to
@@ -895,6 +1008,97 @@ CONTAINS
 
    END SUBROUTINE polar_operators_local
 
+! **************************************************************************************************
+!> \brief Calculate the local dipole operator in the AO basis
+!>        afterwards multiply with the ground state MO coefficients
+!>
+!> \param qs_env ...
+!> \par History
+!>      01.2013 created [SL]
+!>      06.2018 polar_env integrated into qs_env (MK)
+!>      08.2020 TB version (JHU)
+!> \author SL
+! **************************************************************************************************
+   SUBROUTINE polar_tb_operators_local(qs_env)
+
+      TYPE(qs_environment_type), POINTER                 :: qs_env
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_tb_operators_local'
+
+      INTEGER                                            :: blk, handle, i, icol, irow, ispin, nmo, &
+                                                            nspins
+      LOGICAL                                            :: do_raman, found
+      REAL(dp)                                           :: fdir
+      REAL(dp), DIMENSION(3)                             :: ria, rib
+      REAL(dp), DIMENSION(:, :), POINTER                 :: d_block, s_block
+      TYPE(cell_type), POINTER                           :: cell
+      TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER       :: dBerry_psi0
+      TYPE(cp_fm_type), POINTER                          :: mo_coeff
+      TYPE(dbcsr_iterator_type)                          :: iter
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dipmat, matrix_s
+      TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
+      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
+      TYPE(polar_env_type), POINTER                      :: polar_env
+
+      CALL timeset(routineN, handle)
+
+      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
+                      cell=cell, particle_set=particle_set, &
+                      polar_env=polar_env, mos=mos, matrix_s=matrix_s)
+
+      nspins = dft_control%nspins
+
+      CALL get_polar_env(polar_env=polar_env, &
+                         do_raman=do_raman, &
+                         dBerry_psi0=dBerry_psi0)
+
+      IF (do_raman) THEN
+
+         ALLOCATE (dipmat(3))
+         DO i = 1, 3
+            ALLOCATE (dipmat(i)%matrix)
+            CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
+         END DO
+
+         CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
+         DO WHILE (dbcsr_iterator_blocks_left(iter))
+            NULLIFY (s_block, d_block)
+            CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
+            ria = particle_set(irow)%r
+            ria = pbc(ria, cell)
+            rib = particle_set(icol)%r
+            rib = pbc(rib, cell)
+            DO i = 1, 3
+               CALL dbcsr_get_block_p(matrix=dipmat(i)%matrix, &
+                                      row=irow, col=icol, BLOCK=d_block, found=found)
+               CPASSERT(found)
+               fdir = 0.5_dp*(ria(i) + rib(i))
+               d_block = s_block*fdir
+            END DO
+         ENDDO
+         CALL dbcsr_iterator_stop(iter)
+
+         ! Compute the derivative and add the result to mo_derivatives
+         DO ispin = 1, dft_control%nspins ! spin
+            CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo)
+            DO i = 1, 3
+               CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
+                                            dBerry_psi0(i, ispin)%matrix, ncol=nmo)
+            ENDDO !x/y/z-direction
+         END DO
+
+         DO i = 1, 3
+            CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
+         END DO
+         DEALLOCATE (dipmat)
+
+      ENDIF ! do_raman
+
+      CALL timestop(handle)
+
+   END SUBROUTINE polar_tb_operators_local
+
 ! **************************************************************************************************
 !> \brief ...
 !> \param a ...
diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F
index a2d620931c..72f01cb86a 100644
--- a/src/qs_scf_post_tb.F
+++ b/src/qs_scf_post_tb.F
@@ -734,8 +734,8 @@ CONTAINS
          ggamma = ggamma*zphase
          IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
             tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
-            ci = ATAN(tmp)
-            dci = (1.0_dp/(1.0_dp + tmp**2))* &
+            ci = -ATAN(tmp)
+            dci = -(1.0_dp/(1.0_dp + tmp**2))* &
                   (AIMAG(dggamma)*REAL(ggamma, KIND=dp) - AIMAG(ggamma)*REAL(dggamma, KIND=dp))/(REAL(ggamma, KIND=dp))**2
             dipole = MATMUL(cell%hmat, ci)/twopi
             dipole_deriv = MATMUL(cell%hmat, dci)/twopi
@@ -746,8 +746,8 @@ CONTAINS
             ! no pbc(particle_set(i)%r(:),cell) so that the total dipole is the sum of the molecular dipoles
             ria = particle_set(i)%r(:)
             q = charges(i)
-            dipole = dipole - q*(ria - rcc)
-            dipole_deriv(:) = dipole_deriv(:) - q*(particle_set(i)%v(:) - drcc)
+            dipole = dipole + q*(ria - rcc)
+            dipole_deriv(:) = dipole_deriv(:) + q*(particle_set(i)%v(:) - drcc)
          END DO
       END IF
       CALL cp_results_erase(results=results, description=description)
diff --git a/src/xtb_ehess.F b/src/xtb_ehess.F
index 0289cc8aae..c176f344be 100644
--- a/src/xtb_ehess.F
+++ b/src/xtb_ehess.F
@@ -220,7 +220,7 @@ CONTAINS
 
       IF (do_ewald) THEN
          ! add self charge interaction and background charge contribution
-         gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
+         gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:)
          IF (ANY(periodic(:) == 1)) THEN
             gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
          END IF
@@ -250,10 +250,10 @@ CONTAINS
             DO j = 1, nj
                la = laoa(i) + 1
                lb = laob(j) + 1
-               gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
+               gcij(i, j) = gchrg(irow, la, 1) + gchrg(icol, lb, 1)
             END DO
          END DO
-         gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
+         gmij = gmcharge(irow, 1) + gmcharge(icol, 1)
          DO is = 1, SIZE(ks_matrix)
             NULLIFY (ksblock)
             CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
@@ -335,7 +335,7 @@ CONTAINS
             CALL dbcsr_get_block_p(matrix=ks_matrix(is)%matrix, &
                                    row=irow, col=icol, block=ksblock, found=found)
             CPASSERT(found)
-            ksblock = ksblock + 0.5_dp*gmij*sblock
+            ksblock = ksblock + gmij*sblock
          END DO
       ENDDO
       CALL dbcsr_iterator_stop(iter)
diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F
index c8924a88c2..478b16e359 100644
--- a/src/xtb_matrices.F
+++ b/src/xtb_matrices.F
@@ -959,7 +959,7 @@ CONTAINS
                                                             nimg, nkind, ns, nsgf, nspins
       INTEGER, DIMENSION(25)                             :: lao
       INTEGER, DIMENSION(5)                              :: occ
-      LOGICAL                                            :: pass_check
+      LOGICAL                                            :: do_efield, pass_check
       REAL(KIND=dp)                                      :: achrg, chmax, pc_ener, qmmm_el
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mcharge
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: aocg, charges
@@ -1009,6 +1009,7 @@ CONTAINS
 
       energy%hartree = 0.0_dp
       energy%qmmm_el = 0.0_dp
+      energy%efield = 0.0_dp
 
       scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
       nspins = dft_control%nspins
@@ -1024,7 +1025,14 @@ CONTAINS
          END DO
       END DO
 
-      IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
+      IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
+          dft_control%apply_efield_field) THEN
+         do_efield = .TRUE.
+      ELSE
+         do_efield = .FALSE.
+      END IF
+
+      IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
          ! Mulliken charges
          CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
@@ -1056,7 +1064,6 @@ CONTAINS
             END DO
          END DO
          DEALLOCATE (aocg)
-
          ! charge mixing
          IF (dft_control%qs_control%do_ls_scf) THEN
             !
@@ -1069,12 +1076,18 @@ CONTAINS
          DO iatom = 1, natom
             mcharge(iatom) = SUM(charges(iatom, :))
          END DO
+      END IF
 
+      IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
          CALL build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
                                 calculate_forces, just_energy)
+      END IF
 
+      IF (do_efield) THEN
          CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
+      END IF
 
+      IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
          IF (dft_control%qs_control%xtb_control%check_atomic_charges) THEN
             pass_check = .TRUE.
             DO ikind = 1, nkind
@@ -1100,7 +1113,9 @@ CONTAINS
                CPABORT("xTB Charges")
             END IF
          END IF
+      END IF
 
+      IF (dft_control%qs_control%xtb_control%coulomb_interaction .OR. do_efield) THEN
          DEALLOCATE (mcharge, charges)
       END IF
 
diff --git a/tests/DFTB/regtest-debug/TEST_FILES b/tests/DFTB/regtest-debug/TEST_FILES
new file mode 100644
index 0000000000..f8ea2b46b5
--- /dev/null
+++ b/tests/DFTB/regtest-debug/TEST_FILES
@@ -0,0 +1,10 @@
+# runs are executed in the same order as in this file
+# the second field tells which test should be run in order to compare with the last available output
+# e.g. 0 means do not compare anything, running is enough
+#      1 compares the last total energy in the file
+#      for details see cp2k/tools/do_regtest
+ch2o_t01.inp                              86    1.0E-12               0.712276774819E+00
+ch2o_t02.inp                              86    1.0E-12               0.712276774819E+00
+ch2o_t03.inp                              86    1.0E-12               0.755534217819E+00
+ch2o_t04.inp                              86    1.0E-12               0.718756134707E+00
+#EOF
diff --git a/tests/DFTB/regtest-debug/ch2o_t01.inp b/tests/DFTB/regtest-debug/ch2o_t01.inp
new file mode 100644
index 0000000000..cebfb5ee4d
--- /dev/null
+++ b/tests/DFTB/regtest-debug/ch2o_t01.inp
@@ -0,0 +1,59 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD DFTB
+      &DFTB
+        SELF_CONSISTENT    T
+        DISPERSION         F
+        ORTHOGONAL_BASIS   F
+        DO_EWALD           F
+        &PARAMETER
+          PARAM_FILE_PATH  DFTB/scc
+          PARAM_FILE_NAME  scc_parameter
+        &END PARAMETER
+      &END DFTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 20.0 20.0 20.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/DFTB/regtest-debug/ch2o_t02.inp b/tests/DFTB/regtest-debug/ch2o_t02.inp
new file mode 100644
index 0000000000..1317048f27
--- /dev/null
+++ b/tests/DFTB/regtest-debug/ch2o_t02.inp
@@ -0,0 +1,59 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD DFTB
+      &DFTB
+        SELF_CONSISTENT    T
+        DISPERSION         F
+        ORTHOGONAL_BASIS   F
+        DO_EWALD           F
+        &PARAMETER
+          PARAM_FILE_PATH  DFTB/scc
+          PARAM_FILE_NAME  scc_parameter
+        &END PARAMETER
+      &END DFTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 20.0 20.0 20.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/DFTB/regtest-debug/ch2o_t03.inp b/tests/DFTB/regtest-debug/ch2o_t03.inp
new file mode 100644
index 0000000000..51563226ce
--- /dev/null
+++ b/tests/DFTB/regtest-debug/ch2o_t03.inp
@@ -0,0 +1,65 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD DFTB
+      &DFTB
+        SELF_CONSISTENT    T
+        DISPERSION         F
+        ORTHOGONAL_BASIS   F
+        DO_EWALD           T
+        &PARAMETER
+          PARAM_FILE_PATH  DFTB/scc
+          PARAM_FILE_NAME  scc_parameter
+        &END PARAMETER
+      &END DFTB
+    &END QS
+    &POISSON
+      &EWALD
+       EWALD_TYPE SPME
+       ALPHA 0.35
+       GMAX 25
+      &END EWALD
+    &END POISSON
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 5.0 5.0 5.0
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/DFTB/regtest-debug/ch2o_t04.inp b/tests/DFTB/regtest-debug/ch2o_t04.inp
new file mode 100644
index 0000000000..ce60cad5ac
--- /dev/null
+++ b/tests/DFTB/regtest-debug/ch2o_t04.inp
@@ -0,0 +1,66 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD DFTB
+      &DFTB
+        SELF_CONSISTENT    T
+        DISPERSION         F
+        ORTHOGONAL_BASIS   F
+        DO_EWALD           T
+        &PARAMETER
+          PARAM_FILE_PATH  DFTB/scc
+          PARAM_FILE_NAME  scc_parameter
+        &END PARAMETER
+      &END DFTB
+    &END QS
+    &POISSON
+      &EWALD
+       EWALD_TYPE SPME
+       ALPHA 0.35
+       GMAX 25
+      &END EWALD
+    &END POISSON
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/QS/regtest-polar/TEST_FILES b/tests/QS/regtest-polar/TEST_FILES
index 318e2dbcb9..b99550dddf 100644
--- a/tests/QS/regtest-polar/TEST_FILES
+++ b/tests/QS/regtest-polar/TEST_FILES
@@ -7,7 +7,7 @@ h2o_LRraman.inp                                       60      2e-05
 h2o_LRraman_loc.inp                                   60      2e-05                        0.00915
 h2o_LRraman_noort.inp                                 60      2e-05                        0.02330
 H2O_md_polar.inp                                      60      4e-04                        0.63555
-xTB_LRraman.inp                                       60      2e-05                        0.00624
-xTB_LRraman_loc.inp                                   60      2e-05                        0.00868
+xTB_LRraman.inp                                       60      2e-05                        0.80485
+xTB_LRraman_loc.inp                                   60      2e-05                        0.80485
 h2o_LRraman_LRI.inp                                   60      2e-05                        0.00426
 #EOF
diff --git a/tests/QS/regtest-polar/xTB_LRraman.inp b/tests/QS/regtest-polar/xTB_LRraman.inp
index 475503f176..2ab60c038e 100644
--- a/tests/QS/regtest-polar/xTB_LRraman.inp
+++ b/tests/QS/regtest-polar/xTB_LRraman.inp
@@ -29,9 +29,9 @@
     &END CELL
     &COORD
       UNIT bohr
-      O           0.000000       0.000000       0.224953
-      H           0.000000       1.451310      -0.899812
-      H           0.000000      -1.451310      -0.899812
+      O              0.000000       0.224953     0.000000
+      H              1.451310      -0.899812     0.000000
+      H             -1.451310      -0.899812     0.000000
     &END COORD
   &END SUBSYS
 &END FORCE_EVAL
diff --git a/tests/QS/regtest-polar/xTB_LRraman_loc.inp b/tests/QS/regtest-polar/xTB_LRraman_loc.inp
index 3bbac6595f..441a52f068 100644
--- a/tests/QS/regtest-polar/xTB_LRraman_loc.inp
+++ b/tests/QS/regtest-polar/xTB_LRraman_loc.inp
@@ -30,9 +30,9 @@
     &END CELL
     &COORD
       UNIT bohr
-      O           0.000000       0.000000       0.224953
-      H           0.000000       1.451310      -0.899812
-      H           0.000000      -1.451310      -0.899812
+      O              0.000000   0.000000    0.224953
+      H              1.451310   0.000000   -0.899812
+      H             -1.451310   0.000000   -0.899812
     &END COORD
   &END SUBSYS
 &END FORCE_EVAL
diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS
index fc34734edf..e225eb9fdc 100644
--- a/tests/TEST_DIRS
+++ b/tests/TEST_DIRS
@@ -28,6 +28,7 @@ xTB/regtest-2
 xTB/regtest-3
 xTB/regtest-4
 xTB/regtest-5
+xTB/regtest-debug
 QS/regtest-almo-strong
 QS/regtest-elpa-2                                           elpa
 QS/regtest-cdft-5
@@ -62,6 +63,7 @@ QS/regtest-slab
 ATOM/regtest-pseudo
 DFTB/regtest-scc
 DFTB/regtest-scc-2
+DFTB/regtest-debug
 QS/regtest-gpw-3
 QS/regtest-ot
 QS/regtest-gpw-1                                            libint
diff --git a/tests/xTB/regtest-debug/TEST_FILES b/tests/xTB/regtest-debug/TEST_FILES
new file mode 100644
index 0000000000..c1978147d7
--- /dev/null
+++ b/tests/xTB/regtest-debug/TEST_FILES
@@ -0,0 +1,12 @@
+# runs are executed in the same order as in this file
+# the second field tells which test should be run in order to compare with the last available output
+# e.g. 0 means do not compare anything, running is enough
+#      1 compares the last total energy in the file
+#      for details see cp2k/tools/do_regtest
+ch2o_t01.inp                              86    1.0E-12               0.108915750509E+01
+ch2o_t02.inp                              86    1.0E-12               0.108915750509E+01
+ch2o_t03.inp                              86    1.0E-12               0.114591982142E+01
+ch2o_t04.inp                              86    1.0E-12               0.109829662200E+01
+ch2o_t05.inp                              87    1.0E-12               0.226843652140E+02
+ch2o_t06.inp                              87    1.0E-12               0.230117849452E+02
+#EOF
diff --git a/tests/xTB/regtest-debug/ch2o_t01.inp b/tests/xTB/regtest-debug/ch2o_t01.inp
new file mode 100644
index 0000000000..ad7c6f9c8a
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t01.inp
@@ -0,0 +1,51 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 20.0 20.0 20.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/xTB/regtest-debug/ch2o_t02.inp b/tests/xTB/regtest-debug/ch2o_t02.inp
new file mode 100644
index 0000000000..d42d4d4237
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t02.inp
@@ -0,0 +1,51 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 20.0 20.0 20.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/xTB/regtest-debug/ch2o_t03.inp b/tests/xTB/regtest-debug/ch2o_t03.inp
new file mode 100644
index 0000000000..045ef9b8f3
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t03.inp
@@ -0,0 +1,51 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+        DO_EWALD T
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 5.0 5.0 5.0
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/xTB/regtest-debug/ch2o_t04.inp b/tests/xTB/regtest-debug/ch2o_t04.inp
new file mode 100644
index 0000000000..37c91366ea
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t04.inp
@@ -0,0 +1,52 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+        DO_EWALD T
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .TRUE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE T
+   DEBUG_POLARIZABILITY F
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/xTB/regtest-debug/ch2o_t05.inp b/tests/xTB/regtest-debug/ch2o_t05.inp
new file mode 100644
index 0000000000..2e1351adf1
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t05.inp
@@ -0,0 +1,61 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC F
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &PROPERTIES
+    &LINRES
+       PRECONDITIONER FULL_ALL
+       EPS 1.e-10
+       &POLAR
+          DO_RAMAN T
+          PERIODIC_DIPOLE_OPERATOR F
+       &END
+    &END
+  &END
+  &SUBSYS
+    &CELL
+      ABC 20.0 20.0 20.0
+      PERIODIC NONE
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE F
+   DEBUG_POLARIZABILITY T
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END
diff --git a/tests/xTB/regtest-debug/ch2o_t06.inp b/tests/xTB/regtest-debug/ch2o_t06.inp
new file mode 100644
index 0000000000..ab0130b005
--- /dev/null
+++ b/tests/xTB/regtest-debug/ch2o_t06.inp
@@ -0,0 +1,61 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+         DO_EWALD T
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+      &MIXING
+          METHOD DIRECT_P_MIXING
+          ALPHA   0.2
+      &END
+    &END SCF
+    &EFIELD
+    &END
+    &PRINT
+      &MOMENTS ON
+         PERIODIC T
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &PROPERTIES
+    &LINRES
+       PRECONDITIONER FULL_ALL
+       EPS 1.e-10
+       &POLAR
+          DO_RAMAN T
+          PERIODIC_DIPOLE_OPERATOR T
+       &END
+    &END
+  &END
+  &SUBSYS
+    &CELL
+      ABC 8.0 8.0 8.0
+    &END CELL
+    &COORD
+     O     0.051368    0.000000    0.000000     
+     C     1.278612    0.000000    0.000000    
+     H     1.870460    0.939607    0.000000   
+     H     1.870460   -0.939607    0.000000  
+    &END COORD
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT CH2O
+  RUN_TYPE DEBUG
+  PRINT_LEVEL LOW
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES F
+   DEBUG_STRESS_TENSOR F
+   DEBUG_DIPOLE F
+   DEBUG_POLARIZABILITY T
+   STOP_ON_MISMATCH T
+   DE 0.0002
+&END