Fix some recent issues (#1469)

* dipole update TB (minor, still problems with PBC) * Issue #1462 vdW LONG_RANGE_CORRECTION parallel bug fix and regtest Issue #1450 not associated pointer in qs_vxc_atom.F Issue #1449 undefined variable in qs_kind_types.F
cp2k · Apr 14, 2021 · a4ef93f · a4ef93f
1 parent 3bec98f
commit a4ef93f
Show file tree

Hide file tree

Showing 8 changed files with 136 additions and 16 deletions.
diff --git a/src/efield_tb_methods.F b/src/efield_tb_methods.F
@@ -290,8 +290,8 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
       LOGICAL                                            :: found, use_virial
       REAL(KIND=dp)                                      :: charge, dd, dr, fdir, fi
-      REAL(KIND=dp), DIMENSION(3)                        :: fieldpol, fij, forcea, fpolvec, kvec, &
-                                                            qi, rab, ria, rib, rij
+      REAL(KIND=dp), DIMENSION(3)                        :: aqi, fieldpol, fij, forcea, fpolvec, &
+                                                            kvec, qi, rab, ria, rib, rij
       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ds_block, ks_block, p_block, s_block
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dsint
@@ -349,6 +349,27 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
       END DO
       qi = AIMAG(LOG(zi))
       energy%efield = -SUM(fpolvec(:)*qi(:))
+!deb
+      aqi = ATAN2(AIMAG(zi), REAL(zi, KIND=dp))
+      WRITE (60, *) "AIMAG-LOG ", qi
+      WRITE (60, *) "ATAN2     ", aqi
+!     ener_field = 0.0_dp
+!     ti = 0.0_dp
+!     DO idir = 1, 3
+!        ! make sure the total normalized polarization is within [-1:1]
+!        IF (qi(idir) > pi) qi(idir) = qi(idir) - twopi
+!        IF (qi(idir) < -pi) qi(idir) = qi(idir) + twopi
+!        ! now check for log branch
+!        IF (ABS(efield%polarisation(idir) - qi(idir)) > pi) THEN
+!           ti(idir) = (efield%polarisation(idir) - qi(idir))/pi
+!           DO i = 1, 10
+!              qi(idir) = cqi(idir) + SIGN(1.0_dp, ti(idir))*twopi
+!              IF (ABS(efield%polarisation(idir) - qi(idir)) < pi) EXIT
+!           END DO
+!        END IF
+!        ener_field = ener_field + fpolvec(idir)*qi(idir)
+!     END DO
+!deb
 
       IF (.NOT. just_energy) THEN
          CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)

diff --git a/src/input_cp2k_xc.F b/src/input_cp2k_xc.F
@@ -989,7 +989,8 @@ SUBROUTINE create_vdw_potential_section(section)
       CALL section_add_keyword(subsection, keyword)
       CALL keyword_release(keyword)
       CALL keyword_create(keyword, __LOCATION__, name="LONG_RANGE_CORRECTION", &
-                          description="Calculate a long range correction to the DFT-D3 model", &
+                          description="Calculate a long range correction to the DFT-D3 model."// &
+                          " WARNING: Use with care! Only for isotropic dense systems.", &
                           usage="LONG_RANGE_CORRECTION", default_l_val=.FALSE., &
                           lone_keyword_l_val=.TRUE.)
       CALL section_add_keyword(subsection, keyword)

diff --git a/src/qs_dispersion_pairpot.F b/src/qs_dispersion_pairpot.F
@@ -1393,10 +1393,10 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
             WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
             IF (dispersion_env%lrc) THEN
                WRITE (unit_nr, *) " Apply a long range correction"
-               WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
             END IF
             IF (dispersion_env%srb) THEN
                WRITE (unit_nr, *) " Apply a short range bond correction"
+               WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
             END IF
             IF (domol) THEN
                WRITE (unit_nr, *) " Inter-molecule scaling parameter (s8) ", kgc8
@@ -2003,9 +2003,11 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
                END DO
             END DO
             IF (use_virial) THEN
-               DO i = 1, 3
-                  virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
-               END DO
+               IF (para_env%ionode) THEN
+                  DO i = 1, 3
+                     virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
+                  END DO
+               END IF
             END IF
             DEALLOCATE (cnkind)
          END IF

diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F
@@ -1978,7 +1978,7 @@ SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_f
                          "Information provided in the input file regarding POTENTIAL for KIND <"// &
                          TRIM(qs_kind%name)//"> will be ignored!")
 
-         check = ((nb_rep > 0) .OR. explicit_basis)
+         check = ((k_rep > 0) .OR. explicit_basis)
          IF (check) &
             CALL cp_warn(__LOCATION__, &
                          "Information provided in the input file regarding BASIS for KIND <"// &
@@ -3275,6 +3275,7 @@ SUBROUTINE set_pseudo_state(econf, z, ncalc, ncore, nelem)
             END IF
          END IF
       ELSE
+         lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
          ncore = 0
          ncalc = 0
          DO l = 0, lmin

diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F
@@ -721,6 +721,9 @@ SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
                ria = pbc(ria, cell)
                via = particle_set(iat)%v(:)
                q = charges(iat)
+!deb
+! WRITE(6,*) "DEBUG ",q,particle_set(iat)%r(2),ria(2)
+!deb
                DO j = 1, 3
                   gvec = twopi*cell%h_inv(j, :)
                   theta = SUM(ria(:)*gvec(:))
@@ -729,11 +732,22 @@ SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
                   dzeta = -q*CMPLX(-SIN(theta), COS(theta), KIND=dp)**(-q - 1.0_dp)*dtheta
                   dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
                   ggamma(j) = ggamma(j)*zeta
+!deb
+! if(j==2) then
+!     WRITE(6,*) "DEBUG ",iat,gvec(:)
+!     WRITE(6,*) "DEBUG ",iat,theta,zeta,ggamma(j)
+! endif
+!deb
                END DO
             ENDDO
          END DO
          dggamma = dggamma*zphase + ggamma*dzphase
          ggamma = ggamma*zphase
+!deb
+! WRITE(6,*) "DEBUG 1  ",ggamma(1)
+! WRITE(6,*) "DEBUG 2  ",ggamma(2)
+! WRITE(6,*) "DEBUG 3  ",ggamma(3)
+!deb
          IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
             tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
             ci = -ATAN(tmp)
@@ -742,6 +756,24 @@ SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
             dipole = MATMUL(cell%hmat, ci)/twopi
             dipole_deriv = MATMUL(cell%hmat, dci)/twopi
          END IF
+!deb
+! WRITE(6,*) "DIPOLE   ",dipole
+!deb
+
+!           zi(:) = CMPLX(1._dp, 0._dp, dp)
+!           DO ia = 1, SIZE(particle_set)
+!              ria = particle_set(ia)%r
+!              ria = pbc(ria, cell)
+!              DO i = 1, 3
+!                 gvec(:) = twopi*cell%h_inv(i, :)
+!                 dd = SUM(gvec(:)*ria(:))
+!                 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
+!                 zi(i) = zi(i)*zdeta
+!              END DO
+!           END DO
+!           zi = zi*zphase
+!           ci = AIMAG(LOG(zi))/twopi
+
       ELSE
          dipole_type = "non-periodic"
          DO i = 1, SIZE(particle_set)

diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F
@@ -636,14 +636,15 @@ SUBROUTINE calculate_xc_2nd_deriv_atom(local_rho_set, qs_env, xc_section, para_e
          IF (gradient_functional) THEN
             ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho1_h(1:4, 1:na, 1:nr, 1:nspins), &
                       drho_s(1:4, 1:na, 1:nr, 1:nspins), drho1_s(1:4, 1:na, 1:nr, 1:nspins))
-
             ALLOCATE (vxg_h(1:3, 1:na, 1:nr, 1:nspins), vxg_s(1:3, 1:na, 1:nr, 1:nspins))
-            vxg_h = 0.0_dp
-            vxg_s = 0.0_dp
          ELSE
-            NULLIFY (vxg_h, vxg_s)
+            ALLOCATE (drho_h(1, 1, 1, 1), drho1_h(1, 1, 1, 1), &
+                      drho_s(1, 1, 1, 1), drho1_s(1, 1, 1, 1))
+            ALLOCATE (vxg_h(1, 1, 1, 1), vxg_s(1, 1, 1, 1))
             rrho = 0.0_dp
          END IF
+         vxg_h = 0.0_dp
+         vxg_s = 0.0_dp
 
          ! parallelization
          local_loop_limit = get_limit(natom, para_env%num_pe, para_env%mepos)
@@ -732,10 +733,8 @@ SUBROUTINE calculate_xc_2nd_deriv_atom(local_rho_set, qs_env, xc_section, para_e
          ! some cleanup
          NULLIFY (weight)
          DEALLOCATE (rho_h, rho_s, rho1_h, rho1_s, vxc_h, vxc_s)
-         IF (gradient_functional) THEN
-            DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
-            DEALLOCATE (drho1_h, drho1_s)
-         END IF
+         DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
+         DEALLOCATE (drho1_h, drho1_s)
 
          CALL xc_dset_release(deriv_set)
          CALL xc_rho_set_release(rho_set_h)

diff --git a/tests/QS/regtest-dft-vdw-corr-2/TEST_FILES b/tests/QS/regtest-dft-vdw-corr-2/TEST_FILES
@@ -11,6 +11,7 @@ dftd3_t12.inp                                         33    1.0E-14
 dftd3_t13.inp                                         33    1.0E-14              -0.00045037056834
 dftd3_t14.inp                                         33    1.0E-14              -0.00022349123009
 dftd3_t15.inp                                         33    1.0E-14              -0.00022349123009
+dftd3_t16.inp                                         31    1.0E-10              -2.42360772711E-02
 #
 argon-vdW-DF1.inp                                      1      3e-13             -85.04054534692069
 argon-vdW-DF2.inp                                      1      2e-13             -85.20254545254821

diff --git a/tests/QS/regtest-dft-vdw-corr-2/dftd3_t16.inp b/tests/QS/regtest-dft-vdw-corr-2/dftd3_t16.inp
@@ -0,0 +1,63 @@
+&GLOBAL
+  PROJECT dftd3_t16
+  RUN_TYPE ENERGY_FORCE
+  PRINT_LEVEL LOW
+&END GLOBAL
+&FORCE_EVAL
+  METHOD QS
+  STRESS_TENSOR ANALYTICAL
+  &PRINT
+    &STRESS_TENSOR ON
+    &END
+  &END
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_MOLOPT
+    POTENTIAL_FILE_NAME POTENTIAL
+    &MGRID
+      CUTOFF 200
+    &END MGRID
+    &QS
+      METHOD GPW
+    &END QS
+    &SCF
+      SCF_GUESS ATOMIC
+      MAX_SCF 10
+      EPS_SCF 1.0e-5
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE
+      &END XC_FUNCTIONAL
+      &vdW_POTENTIAL
+         DISPERSION_FUNCTIONAL PAIR_POTENTIAL
+         &PAIR_POTENTIAL
+            TYPE DFTD3
+            CALCULATE_C9_TERM .TRUE.
+            REFERENCE_C9_TERM .TRUE.
+            LONG_RANGE_CORRECTION .TRUE.
+            PARAMETER_FILE_NAME dftd3.dat
+            VERBOSE_OUTPUT
+            REFERENCE_FUNCTIONAL PBE
+            R_CUTOFF 8.
+            EPS_CN 0.01
+            KIND_COORDINATION_NUMBERS   0.00   1
+         &END PAIR_POTENTIAL
+      &END vdW_POTENTIAL
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 8.30 8.30 8.30
+    &END
+    &COORD
+ SCALED
+ Ar    0.0 0.0 0.0
+ Ar    0.5 0.5 0.0
+ Ar    0.5 0.0 0.5
+ Ar    0.0 0.5 0.5
+    &END COORD
+    &KIND Ar
+      BASIS_SET SZV-MOLOPT-SR-GTH
+      POTENTIAL GTH-PBE-q8
+    &END
+  &END SUBSYS
+&END