xTB: Displacement field energies and gradient

cp2k · Mar 19, 2019 · e91fdb5 · e91fdb5
1 parent f2247c1
commit e91fdb5
Show file tree

Hide file tree

Showing 4 changed files with 40 additions and 48 deletions.
diff --git a/src/efield_tb_methods.F b/src/efield_tb_methods.F
@@ -574,7 +574,7 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
       fieldpol = fieldpol*dft_control%period_efield%strength
 
       omega = cell%deth
-      hmat = cell%hmat(:, :)/(twopi*omega)
+      hmat = cell%hmat(:, :)/twopi
 
       natom = SIZE(particle_set)
       nspin = SIZE(ks_matrix, 1)
@@ -583,7 +583,6 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
       DO ia = 1, natom
          charge = mcharge(ia)
          ria = particle_set(ia)%r
-         ria = pbc(ria, cell)
          DO idir = 1, 3
             kvec(:) = twopi*cell%h_inv(idir, :)
             dd = SUM(kvec(:)*ria(:))
@@ -595,19 +594,20 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
 
       ! make sure the total normalized polarization is within [-1:1]
       DO idir = 1, 3
-         cqi(idir) = qi(idir)
-!        IF (cqi(idir) > pi) cqi(idir) = cqi(idir)-twopi
-!        IF (cqi(idir) < -pi) cqi(idir) = cqi(idir)+twopi
-!        ! now check for log branch
-!        IF (calculate_forces) THEN
-!           IF (ABS(efield%polarisation(idir)-cqi(idir)) > pi) THEN
-!              di(idir) = (efield%polarisation(idir)-cqi(idir))/pi
-!              DO i = 1, 10
-!                 cqi(idir) = cqi(idir)+SIGN(1.0_dp, di(idir))*twopi
-!                 IF (ABS(efield%polarisation(idir)-cqi(idir)) < pi) EXIT
-!              END DO
-!           END IF
-!        END IF
+         cqi(idir) = qi(idir)/omega
+         IF (cqi(idir) > pi) cqi(idir) = cqi(idir)-twopi
+         IF (cqi(idir) < -pi) cqi(idir) = cqi(idir)+twopi
+         ! now check for log branch
+         IF (calculate_forces) THEN
+            IF (ABS(efield%polarisation(idir)-cqi(idir)) > pi) THEN
+               di(idir) = (efield%polarisation(idir)-cqi(idir))/pi
+               DO i = 1, 10
+                  cqi(idir) = cqi(idir)+SIGN(1.0_dp, di(idir))*twopi
+                  IF (ABS(efield%polarisation(idir)-cqi(idir)) < pi) EXIT
+               END DO
+            END IF
+         END IF
+         cqi(idir) = cqi(idir)*omega
       END DO
       DO idir = 1, 3
          ci(idir) = 0.0_dp
@@ -623,17 +623,19 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
             CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface")
          END IF
          efield%field_energy = ener_field
-         efield%polarisation(:) = cqi(:)
+         efield%polarisation(:) = cqi(:)/omega
       END IF
 
       ! Energy
       ener_field = 0.0_dp
       DO idir = 1, 3
-         ener_field = ener_field+dfilter(idir)*(fieldpol(idir)-2._dp*twopi*ci(idir))**2
+         ener_field = ener_field+dfilter(idir)*(fieldpol(idir)-2._dp*twopi/omega*ci(idir))**2
       END DO
-      energy%efield = 0.25_dp*omega/twopi*ener_field
+      energy%efield = 0.25_dp/twopi*ener_field
 
       IF (.NOT. just_energy) THEN
+         di(:) = -(fieldpol(:)-2._dp*twopi/omega*ci(:))*dfilter(:)/omega
+
          CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
 
@@ -650,15 +652,11 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
             ALLOCATE (atom_of_kind(natom), kind_of(natom))
             CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
             IF (para_env%mepos == 0) THEN
-               DO idir = 1, 3
-                  di(idir) = -omega*(fieldpol(idir)-2._dp*twopi*ci(idir))*dfilter(idir)
-               END DO
-               forcea(1:3) = di(1:3)/omega
                DO ia = 1, natom
                   charge = mcharge(ia)
                   iatom = atom_of_kind(ia)
                   ikind = kind_of(ia)
-                  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom)+forcea(:)*charge
+                  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom)+di(:)*charge
                   IF (use_virial) THEN
                      ria = particle_set(ia)%r
                      ria = pbc(ria, cell)
@@ -675,10 +673,7 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
                CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
 
                DO idir = 1, 3
-                  di(idir) = -omega*(fieldpol(idir)-2._dp*twopi*ci(idir))*dfilter(idir)
-               END DO
-               DO idir = 1, 3
-                  hdi(idir) = SUM(di(1:3)*hmat(1:3, idir))
+                  hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
                END DO
                fdir = 0.0_dp
                ria = particle_set(irow)%r
@@ -736,10 +731,7 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
                CPASSERT(ic > 0)
 
                DO idir = 1, 3
-                  di(idir) = -omega*(fieldpol(idir)-2._dp*twopi*ci(idir))*dfilter(idir)
-               END DO
-               DO idir = 1, 3
-                  hdi(idir) = SUM(di(1:3)*hmat(1:3, idir))
+                  hdi(idir) = -SUM(di(1:3)*hmat(1:3, idir))
                END DO
                fdir = 0.0_dp
                ria = particle_set(irow)%r

diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F
@@ -416,16 +416,6 @@ SUBROUTINE scf_post_calculation_tb(qs_env, tb_type, no_mos)
          END IF
       ENDIF
 
-      ! EFIELD CUBE FILE
-      print_key => section_vals_get_subs_vals(print_section, "EFIELD_CUBE")
-      IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
-         IF (do_cube) THEN
-            CPWARN("Efield cube file not implemented for TB methods.")
-         ELSE
-            CPWARN("Efield cube file not implemented for TB methods.")
-         END IF
-      END IF
-
       ! DENSITY CUBE FILE
       print_key => section_vals_get_subs_vals(print_section, "E_DENSITY_CUBE")
       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
@@ -464,6 +454,16 @@ SUBROUTINE scf_post_calculation_tb(qs_env, tb_type, no_mos)
          END IF
       END IF
 
+      ! EFIELD CUBE FILE
+      print_key => section_vals_get_subs_vals(print_section, "EFIELD_CUBE")
+      IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
+         IF (do_cube) THEN
+            CPWARN("Efield cube file not implemented for TB methods.")
+         ELSE
+            CPWARN("Efield cube file not implemented for TB methods.")
+         END IF
+      END IF
+
       ! ELF
       print_key => section_vals_get_subs_vals(print_section, "ELF_CUBE")
       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN

diff --git a/tests/xTB/regtest-2/HF-dfield-debug.inp b/tests/xTB/regtest-2/HF-dfield-debug.inp
@@ -8,8 +8,8 @@
     &END
     &PERIODIC_EFIELD
        POLARISATION 1 2 5 
-       INTENSITY 0.005
-       DISPLACEMENT_FIELD F
+       INTENSITY 0.01
+       DISPLACEMENT_FIELD T
     &END    
     &QS
       METHOD xTB
@@ -25,7 +25,7 @@
       &END EWALD
     &END POISSON
     &SCF
-      EPS_SCF 1.0E-8
+      EPS_SCF 1.0E-7
       SCF_GUESS ATOMIC 
       MAX_SCF 100
     &END SCF
@@ -50,6 +50,6 @@
 &DEBUG
    DEBUG_FORCES T
    DEBUG_STRESS_TENSOR F
-   STOP_ON_MISMATCH F
-   DX 0.0001
+   STOP_ON_MISMATCH T
+   DX 0.001
 &END DEBUG
diff --git a/tests/xTB/regtest-2/TEST_FILES b/tests/xTB/regtest-2/TEST_FILES
@@ -12,8 +12,8 @@ HF-field.inp                                           1      1e-12
 HF-field-gopt.inp                                      1      1e-12              -5.66101691610414
 HF-field-debug.inp                                     0      1e-12
 HF-dfilter-debug.inp                                   0      1e-12
-HF-dfield-gopt.inp                                     1      1e-12              -5.65726523262359
-HF-dfield.inp                                          1      1e-12              -5.65669037633279
+HF-dfield-gopt.inp                                     1      1e-12              -5.66151709288324
+HF-dfield.inp                                          1      1e-12              -5.66132795820476
 HF-dfield-debug.inp                                    0      1e-12
 HF-loc-field.inp                                       1      1e-12              -5.65286303322644
 HF-loc-field-gopt.inp                                  1      1e-12              -5.66995866072124