diff --git a/src/ddx_constants.f90 b/src/ddx_constants.f90 index 66577adc..d4503ab7 100644 --- a/src/ddx_constants.f90 +++ b/src/ddx_constants.f90 @@ -138,6 +138,9 @@ module ddx_constants !> Derivative of the characteristic function U at all grid points of all !! spheres. Dimension is (3, ngrid, nsph). real(dp), allocatable :: zi(:, :, :) + !> Derivative of the characteristic function U at all grid points of all + !! spheres wrt to radii. Dimension is (ngrid, nsph). + real(dp), allocatable :: zi_dr(:, :) !> Number of external Lebedev grid points on a molecular surface. integer :: ncav !> Number of external Lebedev grid points on each sphere. @@ -950,6 +953,13 @@ subroutine constants_geometry_init(params, constants, ddx_error) return endif constants % zi = zero + + allocate(constants % zi_dr(params % ngrid, params % nsph), stat=info) + if (info .ne. 0) then + call update_error(ddx_error, "`zi_dr` allocation failed") + return + endif + constants % zi_dr = zero end if ! Build arrays fi, ui, zi !$omp parallel do default(none) shared(params,constants,swthr) & @@ -983,6 +993,9 @@ subroutine constants_geometry_init(params, constants, ddx_error) & params % rsph(jsph) / vv constants % zi(:, igrid, isph) = & & constants % zi(:, igrid, isph) + vv*v + constants % zi_dr(igrid, isph) = & + & constants % zi_dr(igrid, isph) & + & + vv*dot_product(v, constants % cgrid(:, igrid)) end if end if enddo @@ -2010,6 +2023,12 @@ subroutine constants_free(constants, ddx_error) call update_error(ddx_error, "`zi` deallocation failed!") end if end if + if (allocated(constants % zi_dr)) then + deallocate(constants % zi_dr, stat=istat) + if (istat .ne. 0) then + call update_error(ddx_error, "`zi_dr` deallocation failed!") + end if + end if if (allocated(constants % ncav_sph)) then deallocate(constants % ncav_sph, stat=istat) if (istat .ne. 0) then diff --git a/src/ddx_cosmo.f90 b/src/ddx_cosmo.f90 index 1bb3b19a..5a597131 100644 --- a/src/ddx_cosmo.f90 +++ b/src/ddx_cosmo.f90 @@ -239,6 +239,68 @@ subroutine cosmo_solvation_force_terms(params, constants, workspace, & end subroutine cosmo_solvation_force_terms +!> Compute the solvation term of the forces and the radial derivatives +!> (solute aspecific). This must be summed to the solute specific term +!! to get the full forces. +!! +!> @ingroup Fortran_interface_ddcosmo +!! @param[in] params: User specified parameters +!! @param[in] constants: Precomputed constants +!! @param[inout] workspace: Preallocated workspaces +!! @param[inout] state: ddx state (contains solutions and RHSs) +!! @param[in] e_cav: electric field, size (3, ncav) +!! @param[inout] force: force term +!! @param[inout] ddx_error: ddX error +!! +subroutine cosmo_solvation_force_dr_terms(params, constants, workspace, & + & state, e_cav, force, dr, ddx_error) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_workspace_type), intent(inout) :: workspace + type(ddx_state_type), intent(inout) :: state + type(ddx_error_type), intent(inout) :: ddx_error + real(dp), intent(in) :: e_cav(3, constants % ncav) + real(dp), intent(inout) :: force(3, params % nsph) + real(dp), intent(inout) :: dr(params % nsph) + ! local variables + real(dp) :: start_time, finish_time + integer :: isph + + ! dummy operation on unused interface arguments + if (ddx_error % flag .eq. 0) continue + + start_time = omp_get_wtime() + + force = zero + do isph = 1, params % nsph + call contract_grad_l(params, constants, isph, state % xs, & + & state % sgrid, workspace % tmp_vylm(:, 1), & + & workspace % tmp_vdylm(:, :, 1), workspace % tmp_vplm(:, 1), & + & workspace % tmp_vcos(:, 1), workspace % tmp_vsin(:, 1), & + & force(:, isph), dr(isph)) + call contract_grad_u(params, constants, isph, state % sgrid, & + & state % phi_grid, force(:, isph), dr(isph)) + end do + + force = -pt5*force + + call zeta_grad(params, constants, state, e_cav, force) + + finish_time = omp_get_wtime() + state % force_time = finish_time - start_time + +end subroutine cosmo_solvation_force_dr_terms + +!> This routines precomputes the intermediates to be used in the evaluation +!! of ddCOSMO analytical derivatives. +!! +!! @param[in] params: ddx parameters +!! @param[in] constant: ddx constants +!! @param[inout] workspace: ddx workspaces +!! @param[inout] state: ddx state +!! + !> This routines precomputes the intermediates to be used in the evaluation !! of ddCOSMO analytical derivatives. !! diff --git a/src/ddx_gradients.f90 b/src/ddx_gradients.f90 index cd98a4f2..ff0f0f07 100644 --- a/src/ddx_gradients.f90 +++ b/src/ddx_gradients.f90 @@ -15,24 +15,38 @@ module ddx_gradients contains !> Compute the gradients of the ddCOSMO matrix -subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx) +subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr) type(ddx_params_type), intent(in) :: params type(ddx_constants_type), intent(in) :: constants integer, intent(in) :: isph - real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma - real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi - real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm - real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc - real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin - real(dp), dimension(3), intent(inout) :: fx + real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma + real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi + real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm + real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc + real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin + real(dp), dimension(3), intent(inout) :: fx + real(dp), optional, intent(inout) :: dr + + logical :: do_dr + real(dp) :: dr_local + + if (present(dr)) then + do_dr = .true. + else + do_dr = .false. + end if + + dr_local = 0.0d0 - call contract_gradi_Lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx ) - call contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx ) + call contract_gradi_Lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx, dr_local) + call contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr_local) + + if (do_dr) dr = dr_local end subroutine contract_grad_L !> Contribution to the gradients of the ddCOSMO matrix -subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx ) +subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr) type(ddx_params_type), intent(in) :: params type(ddx_constants_type), intent(in) :: constants integer, intent(in) :: isph @@ -43,15 +57,19 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin real(dp), dimension(3), intent(inout) :: fx integer :: ig, ij, jsph, l, ind, m - real(dp) :: vvij, tij, xij, oij, t, fac, fl, f1, f2, f3, beta, tlow, thigh + real(dp) :: vvij, tij, xij, oij, t, fac, fl, f1, fr1, f2, fr2, f3, fr3, beta, tlow, thigh real(dp) :: vij(3), sij(3), alp(3), va(3) + real(dp) :: dsij, dtij(3), alp1(3), alp2(3), alp_rad, dsij_rad(3), dtij_rad, alp1_rad, alp2_rad, va_rad real(dp), external :: dnrm2 + real(dp), intent(inout) :: dr + real(dp) :: sjac(3,3), qij tlow = one - pt5*(one - params % se)*params % eta thigh = one + pt5*(one + params % se)*params % eta do ig = 1, params % ngrid va = zero - do ij = constants % inl(isph), constants % inl(isph+1) - 1 + va_rad = zero + do ij = constants % inl(isph), constants % inl(isph+1) - 1 jsph = constants % nl(ij) vij = params % csph(:,isph) + & & params % rsph(isph)*constants % cgrid(:,ig) - & @@ -64,8 +82,37 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc else sij = one end if + + dsij = 1.0_dp / (tij*params%rsph(jsph)) !scalar + dtij(:) = sij(:)/params%rsph(jsph) !vector + + ! build the jacobian of sik + sjac = zero + sjac(1,1) = one + sjac(2,2) = one + sjac(3,3) = one + qij = one/vvij + do icomp = 1, 3 + do jcomp = 1, 3 + sjac(icomp,jcomp) = qij*(sjac(icomp,jcomp) & + & - sij(icomp)*sij(jcomp)) + end do + end do + + dsij_rad = constants%cgrid(:,ig)/vvij & + & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3) + dtij_rad = dot_product(vij, constants%cgrid(:,ig)) & + & / (vvij*params%rsph(jsph)) + call dbasis(params, constants, sij, basloc, dbsloc, vplm, vcos, vsin) alp = zero + alp1 = zero + alp2 = zero + + alp_rad = zero + alp1_rad = zero + alp2_rad = zero + t = one do l = 1, params % lmax ind = l*l + l + 1 @@ -74,7 +121,23 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc do m = -l, l f2 = fac*sigma(ind+m,jsph) f1 = f2*fl*basloc(ind+m) - alp(:) = alp(:) + f1*sij(:) + f2*dbsloc(:,ind+m) + + alp1(:) = f1*dtij + alp1_rad = f1*dtij_rad + + alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + & + & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m) + alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + & + & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m) + alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + & + & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m) + alp2(:) = alp2(:)*tij*f2 + + alp2_rad = f2*tij*dot_product(dsij_rad, dbsloc(:,ind+m)) + + alp(:) = alp(:) + alp1(:) + alp2(:) + alp_rad = alp_rad + alp1_rad + alp2_rad + end do t = t*tij end do @@ -87,20 +150,23 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc oij = xij f2 = zero end if - f1 = oij/params % rsph(jsph) + f1 = oij va(:) = va(:) + f1*alp(:) + beta*f2*constants % zi(:,ig,isph) + va_rad = va_rad + f1*alp_rad + beta*f2*constants % zi_dr(ig,isph) if (tij .gt. tlow) then - f3 = beta*dfsw(tij,params % se,params % eta)/params % rsph(jsph) + f3 = beta*dfsw(tij,params % se,params % eta) if (constants % fi(ig,isph).gt.one) f3 = f3/constants % fi(ig,isph) - va(:) = va(:) + f3*sij(:) + va(:) = va(:) + f3*dtij + va_rad = va_rad + f3*dtij_rad end if end do fx = fx - constants % wgrid(ig)*xi(ig)*va(:) + dr = dr - constants % wgrid(ig)*xi(ig)*va_rad end do end subroutine contract_gradi_Lik !> Contribution to the gradients of the ddCOSMO matrix -subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx ) +subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr) type(ddx_params_type), intent(in) :: params type(ddx_constants_type), intent(in) :: constants integer, intent(in) :: isph @@ -116,8 +182,11 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc real(dp) :: vvji, tji, xji, oji, t, fac, fl, f1, f2, beta, di, tlow, thigh real(dp) :: b, g1, g2, vvjk, tjk, xjk real(dp) :: vji(3), sji(3), alp(3), vb(3), vjk(3), sjk(3), vc(3) + real(dp) :: dsji, dtji(3), alp1(3), alp2(3), alp_rad, dsji_rad(3), dtji_rad, alp1_rad, alp2_rad, vb_rad, vc_rad real(dp) :: rho, ctheta, stheta, cphi, sphi real(dp), external :: dnrm2 + real(dp), intent(inout) :: dr + real(dp) :: sjac(3,3), qji tlow = one - pt5*(one - params % se)*params % eta thigh = one + pt5*(one + params % se)*params % eta @@ -125,6 +194,10 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc do ig = 1, params % ngrid vb = zero vc = zero + + vb_rad = zero + vc_rad = zero + do ji = constants % inl(isph), constants % inl(isph+1) - 1 jsph = constants % nl(ji) vji = params % csph(:,jsph) + & @@ -138,9 +211,36 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc else sji = one end if + + ! build the jacobian of sji + sjac = zero + sjac(1,1) = - one + sjac(2,2) = - one + sjac(3,3) = - one + qji = one/vvji + do icomp = 1, 3 + do jcomp = 1, 3 + sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) & + & + sji(icomp)*sji(jcomp)) + end do + end do + + dsji = 1.0_dp/(tji*params % rsph(isph)) + dtji(:) = -sji/params % rsph(isph) + + dsji_rad = zero + dtji_rad = -tji/params%rsph(isph) + call dbasis(params, constants, sji, basloc, dbsloc, vplm, vcos, vsin) - alp = zero - t = one + alp = zero + alp1 = zero + alp2 = zero + + alp_rad = zero + alp1_rad = zero + alp2_rad = zero + + t = one do l = 1, params % lmax ind = l*l + l + 1 fl = dble(l) @@ -148,18 +248,35 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc do m = -l, l f2 = fac*sigma(ind+m,isph) f1 = f2*fl*basloc(ind+m) - alp = alp + f1*sji + f2*dbsloc(:,ind+m) + + alp1(:) = f1*dtji + alp1_rad = f1*dtji_rad + + alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + & + & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m) + alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + & + & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m) + alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + & + & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m) + alp2(:) = alp2(:)*tji*f2 + + alp2_rad = f2*tji*dot_product(dsji_rad,dbsloc(:,ind+m)) + + alp(:) = alp(:) + alp1(:) + alp2(:) + alp_rad = alp_rad + alp1_rad + alp2_rad + end do t = t*tji end do + xji = fsw(tji, params % se, params % eta) if (constants % fi(ig,jsph).gt.one) then oji = xji/constants % fi(ig,jsph) else oji = xji end if - f1 = oji/params % rsph(isph) - vb = vb + f1*alp*xi(ig,jsph) + vb = vb + oji*alp*xi(ig,jsph) + vb_rad = vb_rad + oji*alp_rad*xi(ig,jsph) if (tji .gt. tlow) then beta = intmlp(params, constants, tji, sigma(:,isph), basloc) if (constants % fi(ig,jsph) .gt. one) then @@ -190,37 +307,55 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc end if end do if (proc) then - g1 = di*di*dfsw(tji, params % se, params % eta)/params % rsph(isph) + g1 = di*di*dfsw(tji, params % se, params % eta) g2 = g1*xi(ig,jsph)*b - vc = vc + g2*sji + vc = vc - g2*dtji + vc_rad = vc_rad - g2*dtji_rad end if else di = one fac = zero end if - f2 = (one-fac)*di*dfsw(tji, params % se, params % eta)/params % rsph(isph) - vb = vb + f2*xi(ig,jsph)*beta*sji + f2 = (one-fac)*di*dfsw(tji, params % se, params % eta) + vb = vb + f2*xi(ig,jsph)*beta*dtji + vb_rad = vb_rad + f2*xi(ig,jsph)*beta*dtji_rad end if end do - fx = fx + constants % wgrid(ig)*(vb - vc) + fx = fx - constants % wgrid(ig)*(vb + vc) + dr = dr - constants % wgrid(ig)*(vb_rad + vc_rad) end do end subroutine contract_gradi_Lji + + !> Gradient of the characteristic function U -subroutine contract_grad_U(params, constants, isph, xi, phi, fx ) +subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr) type(ddx_params_type), intent(in) :: params type(ddx_constants_type), intent(in) :: constants integer, intent(in) :: isph real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi, phi real(dp), dimension(3), intent(inout) :: fx + real(dp), optional, intent(inout) :: dr integer :: ig, ji, jsph real(dp) :: vvji, tji, fac, swthr - real(dp) :: alp(3), vji(3), sji(3) + real(dp) :: alp(3), vji(3), sji(3), dtji(3) + real(dp) :: dtji_rad, alp_rad real(dp), external :: dnrm2 + logical :: do_dr + real(dp) :: dr_local + + if (present(dr)) then + do_dr = .true. + else + do_dr = .false. + end if + do ig = 1, params % ngrid alp = zero + alp_rad = zero if (constants % ui(ig,isph) .gt. zero .and. constants % ui(ig,isph).lt.one) then alp = alp + phi(ig,isph)*xi(ig,isph)*constants % zi(:,ig,isph) + alp_rad = alp_rad + phi(ig,isph)*xi(ig,isph)*constants % zi_dr(ig,isph) end if do ji = constants % inl(isph), constants % inl(isph+1) - 1 jsph = constants % nl(ji) @@ -233,11 +368,21 @@ subroutine contract_grad_U(params, constants, isph, xi, phi, fx ) swthr = one + (params % se + 1.d0)*params % eta / 2.d0 if (tji.lt.swthr .and. tji.gt.swthr-params % eta .and. constants % ui(ig,jsph).gt.zero) then sji = vji/vvji - fac = - dfsw(tji, params % se, params % eta)/params % rsph(isph) - alp = alp + fac*phi(ig,jsph)*xi(ig,jsph)*sji + dtji = - sji / params % rsph(isph) + dtji_rad = - vvji/(params%rsph(isph)**2) + ! dtji_rad = dot_product(vji, constants%cgrid(:,ig)) / (vvji* params%rsph(isph)) + fac = dfsw(tji, params % se, params % eta) + alp = alp + fac*phi(ig,jsph)*xi(ig,jsph) * dtji + alp_rad = alp_rad + fac*phi(ig,jsph)*xi(ig,jsph) * dtji_rad end if end do - fx = fx - constants % wgrid(ig)*alp + ! fx = fx - constants % wgrid(ig)*alp + if (present(dr)) then + dr = dr - constants % wgrid(ig)*alp_rad + else + fx = fx - constants % wgrid(ig)*alp + end if + end do end subroutine contract_grad_U @@ -2237,4 +2382,28 @@ subroutine zeta_grad(params, constants, state, e_cav, forces) end do end subroutine zeta_grad +!> Force term: interaction of the external electric field with the zeta +!! intermediate. This routine is called by the gradients of ddCOSMO +!! ddPCM and ddLPB. +subroutine zeta_grad_dr(params, constants, state, e_cav, dr) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_state_type), intent(inout) :: state + real(dp), intent(inout) :: dr(params % nsph) + real(dp), intent(in) :: e_cav(3, constants % ncav) + ! local variables + integer :: icav, isph, igrid + + icav = 0 + do isph = 1, params % nsph + do igrid = 1, params % ngrid + if (constants % ui(igrid, isph) .eq. zero) cycle + icav = icav + 1 + dr(isph) = dr(isph) + pt5 & + & * dot_product(constants % cgrid(:,igrid), e_cav(:, icav))*state % zeta(icav) + end do + end do +end subroutine zeta_grad_dr + end module ddx_gradients diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 49209055..ac0a7c4a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -92,6 +92,7 @@ ddx_add_test("ddlpb_esolv.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt ddx_add_test("matrix_derivatives.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt") ddx_add_test("matrix_adjoint.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt") ddx_add_test("matrix_solvers.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt") +ddx_add_test("test_gradients.f90" "${CMAKE_CURRENT_SOURCE_DIR}/Input_cosmo_small.txt") ddx_add_test("m2l.f90") ddx_add_test("multipolar_solutes.f90") ddx_add_test("error.f90") diff --git a/tests/Input_cosmo_small.txt b/tests/Input_cosmo_small.txt new file mode 100644 index 00000000..284a3d05 --- /dev/null +++ b/tests/Input_cosmo_small.txt @@ -0,0 +1,21 @@ +'' ! Name of log-file for additional information +1 ! Number of OpenMP cores to be used +1 ! Model. 1 for COSMO, 2 for PCM and 3 for LPB +7 ! Maximal degree of modeling spherical harmonics +302 ! Approximate number of Lebedev grid points +78.3553 ! Dielectric permittivity constant +0.0 ! Shift of the regularized characteristic function +0.1 ! Regularization parameter +0.0 ! Debye H\"{u}ckel parameter +0 ! Whether to compute and store sparse matrices (1) or not (0) +1d-14 ! Relative threshold for the iterative solver +200 ! Maximum number of iterations +25 ! Number of DIIS extrapolation points +1 ! Whether to compute (1) or not (0) forces +0 ! Whether to use (1) or not (0) the FMM +20 ! Max degree of multipole spherical harmonics for the FMM +20 ! Max degree of local spherical harmonics for the FMM +3 # nsph + 1.09581 0.0 0.0 0.0 2.0 + 1.17352 0.99 0.0 0.0 3.0 + 0.33039 -3.37146 -0.33148 0.00013 4.00253 diff --git a/tests/test_gradients.f90 b/tests/test_gradients.f90 new file mode 100644 index 00000000..95c99b3e --- /dev/null +++ b/tests/test_gradients.f90 @@ -0,0 +1,375 @@ +!> @copyright (c) 2020-2021 RWTH Aachen. All rights reserved. +!! +!! ddX software +!! +!! @file src/ddx_driver.f90 +!! Main ddx driver. +!! +!! @version 1.0.0 +!! @author Aleksandr Mikhalev +!! @date 2021-02-25 + +!> Standalone application of ddX: this program can compute the solvation +!! energy and forces for a solute made of point charges. All the +!! relevant steps are here outlined. +program test_gradients +use ddx +use ddx_multipolar_solutes +use omp_lib +implicit none + +real(dp), parameter :: step = 1d-6 +character(len=255) :: dummy_file_name = '' +character(len=255) :: fname +character(len=2047) :: banner +type(ddx_type) :: ddx_data +type(ddx_state_type) :: state +type(ddx_error_type) :: ddx_error +type(ddx_electrostatics_type) :: electrostatics +real(dp), allocatable :: psi(:, :), force(:, :), charges(:), & + & multipoles(:, :), dr(:), numforce(:,:), numdr(:) +real(dp) :: tol, esolv, start_time, finish_time +integer :: i, isph, info, j +real(dp), allocatable :: x(:,:), s(:,:), tmp_lx(:,:) +real(dp) :: xpsi, slx, sphi, diff +real(dp), allocatable :: grad_xpsi(:,:), grad_slx(:,:), & + & grad_sphi(:,:), grad_xpsi_num(:,:), grad_slx_num(:,:), & + & grad_sphi_num(:,:), dr_xpsi(:), dr_slx(:), dr_sphi(:), & + & dr_xpsi_num(:), dr_slx_num(:), dr_sphi_num(:) +real(dp), external :: ddot +real(dp), parameter :: threshold = 1e-8 + + +call get_command_argument(1, fname) +call ddfromfile(fname, ddx_data, tol, charges, ddx_error) +call check_error(ddx_error) + +call allocate_state(ddx_data % params, ddx_data % constants, state, & + & ddx_error) +call check_error(ddx_error) + +allocate( & + & grad_xpsi(3,ddx_data%params%nsph), & + & grad_slx(3,ddx_data%params%nsph), & + & grad_sphi(3,ddx_data%params%nsph), & + & grad_xpsi_num(3,ddx_data%params%nsph), & + & grad_slx_num(3,ddx_data%params%nsph), & + & grad_sphi_num(3,ddx_data%params%nsph), & + & dr_xpsi(ddx_data%params%nsph), & + & dr_slx(ddx_data%params%nsph), & + & dr_sphi(ddx_data%params%nsph), & + & dr_xpsi_num(ddx_data%params%nsph), & + & dr_slx_num(ddx_data%params%nsph), & + & dr_sphi_num(ddx_data%params%nsph), & + & x(ddx_data%constants%nbasis,ddx_data%params%nsph), & + & s(ddx_data%constants%nbasis,ddx_data%params%nsph), & + & tmp_lx(ddx_data%constants%nbasis,ddx_data%params%nsph), & + & force(3, ddx_data % params % nsph), & + & dr(ddx_data % params % nsph), & + & numforce(3, ddx_data % params % nsph), & + & numdr(ddx_data % params % nsph), & + & multipoles(1, ddx_data % params % nsph), stat=info) +if (info .ne. 0) then + write(6, *) "Allocation failed in ddx_driver" + stop 1 +end if + +multipoles(1, :) = charges/sqrt4pi +call multipole_electrostatics(ddx_data % params, ddx_data % constants, & + & ddx_data % workspace, multipoles, 0, electrostatics, ddx_error) + +allocate(psi(ddx_data % constants % nbasis, ddx_data % params % nsph), & + & stat=info) +if (info .ne. 0) then + write(6, *) "Allocation failed in ddx_driver" + stop 1 +end if +call multipole_psi(ddx_data % params, multipoles, 0, psi) + + +! Reference calculation +! call draco(ddx_data%params%csph, ddx_data%params%nsph, & + ! & ddx_data%params%rsph) +call ddrun(ddx_data,state,electrostatics,psi,tol,esolv,ddx_error,force) + + +x(:,:) = state%xs(:,:) +s(:,:) = state%s(:,:) +! Compute the first Lagrangian term 0.5 X^T Psi +xpsi = pt5*ddot(ddx_data%constants%n,x,1,psi,1) + +! Compute the second Lagrangian term -0.5 S^T L X +ddx_data%constants%dodiag = .true. +call lx(ddx_data%params,ddx_data%constants,ddx_data%workspace, & + & x,tmp_lx,ddx_error) +slx = -pt5*ddot(ddx_data%constants%n,s,1,tmp_lx,1) + +! Compute the third Lagrangian term - 0.5 S^T Phi +! Note that Phi is saved as -Phi, so no minus required +sphi = pt5*ddot(ddx_data%constants%n,s,1,state%phi,1) + +write(6,*) xpsi, slx, sphi +! Check if the Lagrangian formulation returns the energy +diff = abs(esolv - (xpsi + slx + sphi)) +if (diff.gt.threshold) then + write(6, *) "Inconsistency" + stop 1 +end if + + +call sgradlx(ddx_data%params,ddx_data%constants,ddx_data%workspace, & + & state,grad_slx,ddx_error) + +call sgradphi(ddx_data%params,ddx_data%constants,ddx_data%workspace, & + & state,grad_sphi,ddx_error,electrostatics%e_cav,multipoles) + +call sdrlx(ddx_data%params,ddx_data%constants,ddx_data%workspace, & + & state,grad_slx,dr_slx,ddx_error) + +call sdrphi(ddx_data%params,ddx_data%constants,ddx_data%workspace, & + & state,grad_sphi,dr_sphi,ddx_error,electrostatics%e_cav,multipoles) + +do isph = 1, ddx_data%params%nsph + do j = 1, 3 + ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) - step + call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, & + & esolv,force,x,s) + + grad_xpsi_num(j, isph) = - xpsi + grad_slx_num(j, isph) = - slx + grad_sphi_num(j, isph) = - sphi + + ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) + 2.0d0*step + call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, & + & esolv,force,x,s) + + grad_xpsi_num(j, isph) = grad_xpsi_num(j, isph) + xpsi + grad_slx_num(j, isph) = grad_slx_num(j, isph) + slx + grad_sphi_num(j, isph) = grad_sphi_num(j, isph) + sphi + grad_xpsi_num(j, isph) = grad_xpsi_num(j, isph)/(2.0d0*step) + grad_slx_num(j, isph) = grad_slx_num(j, isph)/(2.0d0*step) + grad_sphi_num(j, isph) = grad_sphi_num(j, isph)/(2.0d0*step) + + ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) - step + end do +end do + +diff = maxval(abs(grad_slx_num - grad_slx)) +write(6,*) "Difference S^T grad L X", diff +if (diff.gt.threshold) then + write(6, *) "Inconsistency" + stop 1 +end if + +diff = maxval(abs(grad_sphi_num - grad_sphi)) +write(6,*) "Difference S^T grad Phi", diff +if (diff.gt.threshold) then + write(6, *) "Inconsistency" + stop 1 +end if + + +do isph = 1, ddx_data%params%nsph + ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) - step + call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, & + & esolv,force,x,s) + + dr_xpsi_num( isph) = - xpsi + dr_slx_num(isph) = - slx + dr_sphi_num( isph) = - sphi + + ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) + 2.0d0*step + call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, & + & esolv,force,x,s) + + dr_xpsi_num( isph) = dr_xpsi_num( isph) + xpsi + dr_slx_num(isph) = dr_slx_num(isph) + slx + dr_sphi_num( isph) = dr_sphi_num(isph) + sphi + dr_xpsi_num( isph) = dr_xpsi_num(isph)/(2.0d0*step) + dr_slx_num(isph) = dr_slx_num(isph)/(2.0d0*step) + dr_sphi_num( isph) = dr_sphi_num(isph)/(2.0d0*step) + + ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) - step +end do + +write(6,*) "-----" +diff = maxval(abs(dr_slx_num - dr_slx)) +write(6,*) "Difference S^T dr L X", diff +if (diff.gt.threshold) then + write(6, *) "Inconsistency" + stop 1 +end if +diff = maxval(abs(dr_sphi_num - dr_sphi)) +write(6,*) "Difference S^T dr Phi", diff +if (diff.gt.threshold) then + write(6, *) "Inconsistency" + stop 1 +end if + +deallocate(psi, multipoles, charges, force, numforce, dr, numdr,& + & x, s, tmp_lx, & + & grad_xpsi, grad_slx, grad_sphi, & + & grad_xpsi_num, grad_slx_num, grad_sphi_num, & + & dr_xpsi, dr_slx, dr_sphi, & + & dr_xpsi_num, dr_slx_num, dr_sphi_num, & + & stat=info) +if (info .ne. 0) then + write(6, *) "Deallocation failed in ddx_driver" + stop 1 +end if +call deallocate_electrostatics(electrostatics, ddx_error) +call deallocate_state(state, ddx_error) +call deallocate_model(ddx_data, ddx_error) + +contains + +subroutine displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi, & + & psi,esolv,force,x,s) + type(ddx_type), intent(inout) :: ddx_data + real(dp), intent(in) :: tol + real(dp), intent(out) :: xpsi, slx, sphi + type(ddx_type) :: ddx_data2 + type(ddx_error_type) :: error2 + type(ddx_state_type) :: state2 + type(ddx_electrostatics_type) :: electrostatics2 + real(dp), intent(in) :: multipoles(1,ddx_data%params%nsph) + real(dp), intent(out) :: esolv + real(dp), intent(out) :: force(3,ddx_data%params%nsph) + real(dp), intent(out) :: psi(ddx_data%constants%nbasis, & + & ddx_data%params%nsph) + real(dp), intent(in) :: & + & x(ddx_data%constants%nbasis,ddx_data%params%nsph), & + & s(ddx_data%constants%nbasis,ddx_data%params%nsph) + + ! Make a copy of the model + call allocate_model(ddx_data%params%nsph,ddx_data%params%csph(1,:), & + & ddx_data%params%csph(2,:),ddx_data%params%csph(3,:), & + & ddx_data%params%rsph,ddx_data%params%model, & + & ddx_data%params%lmax,ddx_data%params%ngrid,ddx_data%params%force, & + & ddx_data%params%fmm,ddx_data%params%pm,ddx_data%params%pl, & + & ddx_data%params%se,ddx_data%params%eta,ddx_data%params%eps, & + & ddx_data%params%kappa,ddx_data%params%matvecmem, & + & ddx_data%params%maxiter,ddx_data%params%jacobi_ndiis, & + & ddx_data%params%nproc,dummy_file_name,ddx_data2,error2) + call check_error(error2) + + call allocate_state(ddx_data2%params,ddx_data2%constants,state2, & + & error2) + call check_error(error2) + + + call multipole_electrostatics(ddx_data2%params,ddx_data2%constants, & + & ddx_data2%workspace,multipoles,0,electrostatics2,error2) + call check_error(error2) + + call multipole_psi(ddx_data2%params,multipoles,0,psi) + call ddrun(ddx_data2,state2,electrostatics2,psi,tol,esolv,ddx_error,force) + + ! Compute the first Lagrangian term 0.5 X^T Psi + xpsi = pt5*ddot(ddx_data2%constants%n,x,1,psi,1) + ! Compute the second Lagrangian term -0.5 S^T L X + ddx_data2%constants%dodiag = .true. + call lx(ddx_data2%params,ddx_data2%constants,ddx_data2%workspace, & + & x,tmp_lx,error2) + slx = -pt5*ddot(ddx_data2%constants%n,s,1,tmp_lx,1) + ! Compute the third Lagrangian term - 0.5 S^T Phi + ! Note that Phi is saved as -Phi, so no minus required + sphi = pt5*ddot(ddx_data2%constants%n,s,1,state2%phi,1) + + call deallocate_electrostatics(electrostatics2,error2) + call deallocate_state(state2,error2) + call deallocate_model(ddx_data2,error2) +end subroutine displaced_run + +subroutine sgradlx(params, constants, workspace, & + & state, force, ddx_error) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_workspace_type), intent(inout) :: workspace + type(ddx_state_type), intent(inout) :: state + type(ddx_error_type), intent(inout) :: ddx_error + real(dp), intent(inout) :: force(3, params % nsph) + integer :: isph + force = zero + do isph = 1, params % nsph + call contract_grad_l(params,constants,isph,state%xs, & + & state%sgrid,workspace%tmp_vylm(:,1), & + & workspace%tmp_vdylm(:,:,1),workspace%tmp_vplm(:,1), & + & workspace%tmp_vcos(:,1),workspace%tmp_vsin(:,1), & + & force(:,isph)) + end do + force = -pt5*force +end subroutine sgradlx + +subroutine sdrlx(params, constants, workspace, & + & state, force, dr, ddx_error) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_workspace_type), intent(inout) :: workspace + type(ddx_state_type), intent(inout) :: state + type(ddx_error_type), intent(inout) :: ddx_error + real(dp), intent(inout) :: force(3, params % nsph) + real(dp), intent(inout) :: dr(params % nsph) + integer :: isph + force = zero + dr = zero + do isph = 1, params % nsph + call contract_grad_l(params,constants,isph,state%xs, & + & state%sgrid,workspace%tmp_vylm(:,1), & + & workspace%tmp_vdylm(:,:,1),workspace%tmp_vplm(:,1), & + & workspace%tmp_vcos(:,1),workspace%tmp_vsin(:,1), & + & force(:,isph), dr=dr(isph)) + end do + dr = -pt5*dr + force = -pt5*force +end subroutine sdrlx + +subroutine sgradphi(params, constants, workspace, & + & state,force,ddx_error,e_cav,multipoles) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_workspace_type), intent(inout) :: workspace + type(ddx_state_type), intent(inout) :: state + type(ddx_error_type), intent(inout) :: ddx_error + real(dp), intent(inout) :: force(3, params % nsph) + real(dp), intent(in) :: e_cav(3, constants % ncav) + real(dp), intent(in) :: multipoles(1, ddx_data % params % nsph) + integer :: isph + force = zero + do isph = 1, params % nsph + call contract_grad_u(params, constants, isph, state % sgrid, & + & state % phi_grid, force(:, isph)) + end do + force = -pt5*force + call zeta_grad(params, constants, state, e_cav, force) + call multipole_force_terms(params,constants, & + & workspace,state,0,multipoles,force,ddx_error) +end subroutine sgradphi + + +subroutine sdrphi(params, constants, workspace, & + & state,force,dr,ddx_error,e_cav,multipoles) + implicit none + type(ddx_params_type), intent(in) :: params + type(ddx_constants_type), intent(in) :: constants + type(ddx_workspace_type), intent(inout) :: workspace + type(ddx_state_type), intent(inout) :: state + type(ddx_error_type), intent(inout) :: ddx_error + real(dp), intent(inout) :: force(3, params % nsph) + real(dp), intent(inout) :: dr(params % nsph) + real(dp), intent(in) :: e_cav(3, constants % ncav) + real(dp), intent(in) :: multipoles(1, ddx_data % params % nsph) + integer :: isph + dr = zero + do isph = 1, params % nsph + call contract_grad_u(params, constants, isph, state % sgrid, & + & state % phi_grid, force(:, isph), dr=dr(isph)) + end do + dr = -pt5*dr + call zeta_grad_dr(params, constants, state, e_cav, dr) +end subroutine sdrphi + +end program test_gradients