Add testing of direct gradient for CN dependent function

src/mctc/ncoord/type.f90

@@ -290,7 +290,7 @@ subroutine add_coordination_number_derivs(self, mol, trans, dEdcn, gradient, sig
                countd = den * self%ncoord_dcount(izp, jzp, r1) * rij/r1
 
                gradient(:, iat) = gradient(:, iat) + countd &
-                  & * (dEdcn(iat) * self%directed_factor + dEdcn(jat))
+                  & * (dEdcn(iat) + dEdcn(jat) * self%directed_factor)
                gradient(:, jat) = gradient(:, jat) - countd &
                   & * (dEdcn(iat) + dEdcn(jat) * self%directed_factor)
 

test/test_ncoord.f90

@@ -87,6 +87,7 @@ subroutine collect_ncoord(testsuite)
       & new_unittest("dcndL-mb07_erf_en", test_dcndL_mb07_erf_en), &
       & new_unittest("dcndL-antracene_erf_en", test_dcndL_anthracene_erf_en), &
       & new_unittest("cn-mb01_erf_dftd4", test_cn_mb01_erf_dftd4), &
+      & new_unittest("dfdcn-mb01_erf_dftd4", test_dfdcn_mb01_erf_dftd4), &
       & new_unittest("cn-mb02_erf_dftd4", test_cn_mb02_erf_dftd4), &
       & new_unittest("cn-mb03_erf_dftd4", test_cn_mb03_erf_dftd4), &
       & new_unittest("cn-acetic_erf_dftd4", test_cn_acetic_erf_dftd4), &
@@ -1354,6 +1355,72 @@ subroutine test_cn_mb01_erf_dftd4(error)
 
    end subroutine test_cn_mb01_erf_dftd4
 
+   subroutine test_dfdcn_mb01_erf_dftd4(error)
+
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+
+      type(structure_type) :: mol
+      class(ncoord_type), allocatable :: erf_dftd4_ncoord
+      real(wp), allocatable :: rcov(:)
+      real(wp), allocatable :: en(:)
+      real(wp) :: gradient(3, 16), sigma(3,3), dEdcn(16)
+      real(wp), allocatable :: cn(:)
+      real(wp), allocatable :: lattr(:, :)
+
+      real(wp), parameter :: cutoff = 30.0_wp
+      real(wp), parameter :: ref_gradient(3, 16) = reshape([ &
+      &  2.3515337861584880_wp, -3.0453827536691351_wp,  0.5392536733971742_wp, & 
+      & -0.4086227335489838_wp,  0.0501074408829545_wp,  0.2137313744331271_wp, & 
+      &  0.0249303991608369_wp, -0.0373339882416463_wp, -0.0756057413224991_wp, & 
+      & -0.5436172455452489_wp,  0.4933244472902760_wp,  1.2499429094951524_wp, & 
+      &  0.0187062301843386_wp,  0.0079863446319139_wp, -0.0342050060044000_wp, & 
+      &  0.0541178064781164_wp,  0.0314514945594538_wp,  0.0256701683856457_wp, & 
+      & -1.1131746078453495_wp, -1.6288054931512879_wp,  0.4146928754871603_wp, & 
+      &  0.0887505649437715_wp,  0.0350320577127468_wp, -0.0407677557149725_wp, & 
+      & -0.0499761354250785_wp,  0.1548143128637227_wp, -0.2724294261044626_wp, & 
+      & -0.2140027394072616_wp, -0.0253209254791824_wp, -0.0023845891460968_wp, & 
+      &  0.0091192286107324_wp, -0.0054889861418698_wp, -0.0708601381531837_wp, & 
+      &  0.1710608463902656_wp,  0.5486016273258228_wp,  0.0456750463319099_wp, & 
+      & -1.4496220571725840_wp,  0.2848532829040934_wp,  2.1631929338961609_wp, & 
+      &  0.2437819446177886_wp,  0.0864993947414499_wp, -0.8074962768633803_wp, & 
+      &  0.0760143426254694_wp,  0.2932453158431059_wp, -0.1384011976765554_wp, & 
+      &  0.7410003697746987_wp,  2.7564164279275825_wp, -3.2100088504407793_wp], &
+      & shape(ref_gradient))
+      real(wp), parameter :: ref_sigma(3, 3) = reshape([ &
+      &-14.102412014163017_wp,  5.058593018200457_wp,  7.2112237247416280_wp, &      
+      &  5.058593018200459_wp,-18.111874478203490_wp,  6.8763096498221916_wp, &        
+      &  7.211223724741628_wp,  6.876309649822192_wp, -22.636906090720050_wp], &
+      & shape(ref_sigma))
+
+      dEdcn = 1.0_wp
+      gradient = 0.0_wp
+      sigma = 0.0_wp
+
+      call get_structure(mol, "mindless01")
+
+      allocate(rcov(mol%nid), en(mol%nid), cn(mol%nat))
+      rcov(:) = get_covalent_rad(mol%num)
+      en(:) = get_pauling_en(mol%num)
+
+      ! Test also the external interface
+      call new_ncoord(erf_dftd4_ncoord, mol, "dftd4", &
+         & cutoff=cutoff, rcov=rcov, en=en)
+      call get_lattice_points(mol%periodic, mol%lattice, cutoff, lattr)
+      call erf_dftd4_ncoord%add_coordination_number_derivs(mol, lattr, &
+         & dEdcn, gradient, sigma)
+
+      if (any(abs(gradient - ref_gradient) > thr)) then
+         call test_failed(error, "Coordination number gradient does not match")
+         print'(3es21.14)', gradient - ref_gradient
+      end if
+
+      if (any(abs(sigma - ref_sigma) > thr)) then
+         call test_failed(error, "Coordination numbers sigma does not match")
+         print'(3es21.14)', sigma - ref_sigma
+      end if
+
+   end subroutine test_dfdcn_mb01_erf_dftd4
 
    subroutine test_cn_mb02_erf_dftd4(error)