Merge branch 'dftd3+pbc' into 'develop'

Refold atoms around the origin before computing DFT-D3 energies and forces

See merge request QEF/q-e!2133

PP/src/d3hess.f90

@@ -15,7 +15,7 @@ program d3hess
   USE mp_world,         ONLY: world_comm
   USE environment,      ONLY: environment_start, environment_end
   !
-  USE cell_base,        ONLY: alat, at
+  USE cell_base,        ONLY: alat, at, bg
   USE ions_base,        ONLY: nat, tau, ityp, atm
   USE funct,            ONLY: get_dft_short
   USE dftd3_api,        ONLY: dftd3_init, dftd3_set_functional, get_atomic_number, dftd3_pbc_dispersion
@@ -103,7 +103,12 @@ program d3hess
   ALLOCATE( xyz(3,nat), atnum(nat), force_d3(3,nat) )
   force_d3(:,:) = 0.0_DP
   latvecs(:,:)=at(:,:)*alat
-  xyz(:,:)=tau(:,:)*alat
+  ! xyz are atomic positions centered around r=0 (in bohr units)
+  xyz(:,:) = tau(:,:)
+  CALL cryst_to_cart( nat, xyz, bg, -1 ) 
+  xyz(:,:) = xyz(:,:) - NINT(xyz(:,:))
+  CALL cryst_to_cart( nat, xyz, at,  1 )
+  xyz(:,:) = xyz(:,:)*alat
   DO iat = 1, nat
      atnum(iat) = get_atomic_number(TRIM(atm(ityp(iat))))
   ENDDO

PW/src/electrons.f90

@@ -539,8 +539,8 @@ SUBROUTINE electrons_scf ( printout, exxen )
   REAL(DP), EXTERNAL :: ewald, get_clock
   REAL(DP) :: etot_cmp_paw(nat,2,2)
   ! 
-  REAL(DP) :: latvecs(3,3)
-  !! auxiliary variables for grimme-d3
+  REAL(DP), ALLOCATABLE :: taupbc(:,:)
+  !! atomic positions centered around r=0 - for grimme-d3
   INTEGER:: atnum(1:nat), na
   !! auxiliary variables for grimme-d3
   LOGICAL :: lhb
@@ -586,14 +586,18 @@ SUBROUTINE electrons_scf ( printout, exxen )
   !
   IF (ldftd3) THEN
      CALL start_clock('energy_dftd3')
-     latvecs(:,:)=at(:,:)*alat
-     tau(:,:)=tau(:,:)*alat
+     ! taupbc are atomic positions in alat units, centered around r=0
+     ALLOCATE ( taupbc(3,nat) )
+     taupbc(:,:) = tau(:,:)
+     CALL cryst_to_cart( nat, taupbc, bg, -1 ) 
+     taupbc(:,:) = taupbc(:,:) - NINT(taupbc(:,:))
+     CALL cryst_to_cart( nat, taupbc, at,  1 ) 
      DO na = 1, nat
         atnum(na) = get_atomic_number(TRIM(atm(ityp(na))))
      ENDDO
-     call dftd3_pbc_dispersion(dftd3,tau,atnum,latvecs,edftd3)
+     call dftd3_pbc_dispersion(dftd3, alat*taupbc, atnum, alat*at, edftd3)
      edftd3=edftd3*2.d0
-     tau(:,:)=tau(:,:)/alat
+     DEALLOCATE( taupbc)
      CALL stop_clock('energy_dftd3')
   ELSE
      edftd3= 0.0

PW/src/forces.f90

@@ -99,7 +99,7 @@ SUBROUTINE forces()
   ! counter on polarization
   ! counter on atoms
   !
-  REAL(DP) :: latvecs(3,3)
+  REAL(DP), ALLOCATABLE :: taupbc(:,:)
   INTEGER :: atnum(1:nat)
   REAL(DP) :: stress_dftd3(3,3)
   !
@@ -176,13 +176,17 @@ SUBROUTINE forces()
     CALL start_clock('force_dftd3')
     ALLOCATE( force_d3(3, nat) )
     force_d3(:,:) = 0.0_DP
-    latvecs(:,:) = at(:,:)*alat
-    tau(:,:) = tau(:,:)*alat
+    ! taupbc are atomic positions in alat units, centered around r=0
+    ALLOCATE ( taupbc(3,nat) )
+    taupbc(:,:) = tau(:,:)
+    CALL cryst_to_cart( nat, taupbc, bg, -1 ) 
+    taupbc(:,:) = taupbc(:,:) - NINT(taupbc(:,:))
+    CALL cryst_to_cart( nat, taupbc, at,  1 ) 
     atnum(:) = get_atomic_number(atm(ityp(:)))
-    CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, &
+    CALL dftd3_pbc_gdisp( dftd3, alat*taupbc, atnum, alat*at, &
                           force_d3, stress_dftd3 )
     force_d3 = -2.d0*force_d3
-    tau(:,:) = tau(:,:)/alat
+    DEALLOCATE( taupbc)
     CALL stop_clock('force_dftd3')
   ENDIF
   !

PW/src/stress.f90

@@ -58,7 +58,7 @@ SUBROUTINE stress( sigma )
   ! ... Auxiliary variables for Grimme-D3
   !
   INTEGER  :: atnum(1:nat)
-  REAL(DP) :: latvecs(3,3)
+  REAL(DP), ALLOCATABLE :: taupbc(:,:)
   REAL(DP), ALLOCATABLE :: force_d3(:,:)
   !
   WRITE( stdout, '(//5x,"Computing stress (Cartesian axis) and pressure"/)' )
@@ -137,13 +137,17 @@ SUBROUTINE stress( sigma )
     CALL start_clock('stres_dftd3')
     ALLOCATE( force_d3(3,nat) )
     force_d3( : , : ) = 0.0_DP
-    latvecs(:,:) = at(:,:)*alat
-    tau(:,:) = tau(:,:)*alat
+    ! taupbc are atomic positions in alat units, centered around r=0
+    ALLOCATE ( taupbc(3,nat) )
+    taupbc(:,:) = tau(:,:)
+    CALL cryst_to_cart( nat, taupbc, bg, -1 ) 
+    taupbc(:,:) = taupbc(:,:) - NINT(taupbc(:,:))
+    CALL cryst_to_cart( nat, taupbc, at,  1 ) 
     atnum(:) = get_atomic_number(atm(ityp(:)))
-    CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, &
+    CALL dftd3_pbc_gdisp( dftd3, alat*taupbc, atnum, alat*at, &
                          force_d3, sigmad23 )
     sigmad23 = 2.d0*sigmad23
-    tau(:,:)=tau(:,:)/alat
+    DEALLOCATE( taupbc )
     DEALLOCATE( force_d3 )
     CALL stop_clock('stres_dftd3')
   END IF