add_coordination_number_hessian Subroutine

public subroutine add_coordination_number_hessian(mol, trans, cutoff, rcov, dEdcn, hessian)

Add the second derivative of the exponential coordination number contracted with the derivative of the energy w.r.t. the coordination number.

The counting function is reproduced here rather than taken from mctc-lib, which only provides derivatives up to first order.

Arguments

Type IntentOptional Attributes Name
type(structure_type), intent(in) :: mol

Molecular structure data
real(kind=wp), intent(in) :: trans(:,:)

Lattice points
real(kind=wp), intent(in) :: cutoff

Real space cutoff
real(kind=wp), intent(in) :: rcov(:)

Covalent radius
real(kind=wp), intent(in) :: dEdcn(:)

Derivative of expression with respect to the coordination number
real(kind=wp), intent(inout) :: hessian(:,:)

Second derivative of the energy w.r.t. the Cartesian coordinates

Source Code

subroutine add_coordination_number_hessian(mol, trans, cutoff, rcov, dEdcn, hessian)

   !> Molecular structure data
   type(structure_type), intent(in) :: mol

   !> Lattice points
   real(wp), intent(in) :: trans(:, :)

   !> Real space cutoff
   real(wp), intent(in) :: cutoff

   !> Covalent radius
   real(wp), intent(in) :: rcov(:)

   !> Derivative of expression with respect to the coordination number
   real(wp), intent(in) :: dEdcn(:)

   !> Second derivative of the energy w.r.t. the Cartesian coordinates
   real(wp), intent(inout) :: hessian(:, :)

   integer :: iat, jat, izp, jzp, itr, ic, jc, ii, jj
   real(wp) :: vec(3), r2, r1, cutoff2, rc, expterm, cf, dcf, d2cf
   real(wp) :: dEdcnij, block(3, 3), tmp

   cutoff2 = cutoff*cutoff

   do iat = 1, mol%nat
      izp = mol%id(iat)
      do jat = 1, iat - 1
         jzp = mol%id(jat)
         rc = rcov(izp) + rcov(jzp)
         dEdcnij = dEdcn(iat) + dEdcn(jat)
         do itr = 1, size(trans, 2)
            vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr))
            r2 = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
            if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle
            r1 = sqrt(r2)

            expterm = exp(-default_kcn*(rc/r1 - 1.0_wp))
            cf = 1.0_wp/(1.0_wp + expterm)
            ! first and second derivative of the counting function w.r.t. r
            tmp = cf * (1.0_wp - cf)
            dcf = -tmp * default_kcn * rc / r2
            d2cf = tmp * (1.0_wp - 2.0_wp*cf) * (default_kcn*rc)**2 / (r2*r2) &
               & + 2.0_wp * tmp * default_kcn * rc / (r2*r1)

            do ic = 1, 3
               do jc = 1, 3
                  block(ic, jc) = dEdcnij * ( &
                     & d2cf * vec(ic) * vec(jc) / r2 &
                     & - dcf * vec(ic) * vec(jc) / (r2*r1))
               end do
               block(ic, ic) = block(ic, ic) + dEdcnij * dcf / r1
            end do

            do ic = 1, 3
               ii = 3*(iat - 1) + ic
               do jc = 1, 3
                  jj = 3*(jat - 1) + jc
                  hessian(ii, 3*(iat - 1) + jc) = hessian(ii, 3*(iat - 1) + jc) + block(ic, jc)
                  hessian(3*(jat - 1) + ic, jj) = hessian(3*(jat - 1) + ic, jj) + block(ic, jc)
                  hessian(ii, jj) = hessian(ii, jj) - block(ic, jc)
                  hessian(3*(jat - 1) + ic, 3*(iat - 1) + jc) = &
                     & hessian(3*(jat - 1) + ic, 3*(iat - 1) + jc) - block(ic, jc)
               end do
            end do
         end do
      end do
   end do

end subroutine add_coordination_number_hessian