diff --git a/Makefile b/Makefile index 9ae4c71..6cec241 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ -IRPF90 = irpf90 #-a -d +IRPF90 = irpf90 --codelet=factor_een:100000 FC = gfortran -FCFLAGS= -O2 -ffree-line-length-none -I . +FCFLAGS= -O2 -march=native -ffree-line-length-none -I . NINJA = ninja ARCHIVE= ar crs RANLIB = ranlib diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f index 07a1e3a..038c4fc 100644 --- a/el_nuc_el.irp.f +++ b/el_nuc_el.irp.f @@ -1,208 +1,121 @@ -BEGIN_PROVIDER [double precision, factor_een] +BEGIN_PROVIDER [ double precision, factor_een ] implicit none BEGIN_DOC - ! Electron-electron nucleus contribution to Jastrow factor + ! END_DOC - integer :: i, j, alpha, p, k, l, lmax, cindex - double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0 + integer :: i,j,a,p,k,l,lmax,m + double precision :: riam, rjam_cn, rial, rjal, rijk + double precision :: cn - PROVIDE cord_vect factor_een = 0.0d0 - do alpha = 1, nnuc - do j = 1, nelec - b = rescale_een_n(j, alpha) - do i = 1, nelec - u = rescale_een_e(i, j) - a = rescale_een_n(i, alpha) - a2 = a * a - b2 = b * b - c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha) - c_inv = 1.0d0 / c - cindex = 0 - do p = 2, ncord - x = 1.0d0 - do k = 0, p - 1 - if ( k /= 0 ) then - lmax = p - k - else - lmax = p - k - 2 - end if - t = x - do l = 1, rshift(p - k, 1) - t = t * c - end do - ! We have suppressed this from the following loop: - ! if ( iand(p - k - l, 1) == 0 ) then - ! - ! Start from l=0 when p-k is even - ! Start from l=1 when p-k is odd - if (iand(p - k, 1) == 0) then - y = 1.0d0 - z = 1.0d0 - else - y = a - z = b - endif - do l = iand(p - k, 1), lmax, 2 - ! This can be used in case of a flatten cord_vect - ! cidx = 1 + l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * (p - 1) + & - ! (ncord + 1) * (ncord + 1) * ncord * (alpha - 1) - cindex = cindex + 1 - factor_een = factor_een + cord_vect(cindex, typenuc_arr(alpha)) * (y + z) * t - t = t * c_inv - y = y * a2 - z = z * b2 - end do - x = x * u - end do - end do - end do - end do - end do + do p=2,ncord + do k=0,p-1 + if (k /= 0) then + lmax = p-k + else + lmax = p-k-2 + endif - factor_een = 0.5d0 * factor_een + do l=0,lmax + if ( iand(p-k-l,1) == 1) then + cycle + endif + m = (p-k-l)/2 + + do a=1, nnuc + cn = cord_vect_lkp(l,k,p,typenuc_arr(a)) + do j=1, nelec + rjal = rescale_een_n(j,a,l) + rjam_cn = rescale_een_n(j,a,m) * cn + do i=1, j-1 + rial = rescale_een_n(i,a,l) + riam = rescale_een_n(i,a,m) + rijk = rescale_een_e(i,j,k) + factor_een = factor_een + & + rijk * (rial+rjal) * riam * rjam_cn + enddo + enddo + enddo + + enddo + enddo + enddo END_PROVIDER -BEGIN_PROVIDER [double precision, factor_een_naive] +BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4,nelec) ] implicit none BEGIN_DOC - ! Electron-electron nucleus contribution to Jastrow factor in a naive way + ! dimensions 1-3 : dx,dy,dz + ! + ! Rdimension 4 : d2x + d2y + d2z END_DOC - integer :: i, j, alpha, p, k, l, lmax, cindex - double precision :: ria, rja, rij + integer :: i,j,a,p,k,l,lmax,m + double precision :: riam, rjam_cn, rial, rjal, rijk + double precision, dimension(4) :: driam, drjam_cn, drial, drjal, drijk + double precision :: cn, v1, v2, l1, l2, d1, d2 - PROVIDE cord_vect - factor_een_naive = 0.0d0 - - do alpha = 1, nnuc - do j = 2, nelec - rja = rescale_een_n(j, alpha) - do i = 1, j - 1 - ria = rescale_een_n(i, alpha) - rij = rescale_een_e(i, j) - cindex = 0 - do p = 2, ncord - do k = p - 1, 0, -1 - if ( k /= 0 ) then - lmax = p - k - else - lmax = p - k - 2 - end if - do l = lmax, 0, -1 - if ( iand(p - k - l, 1) == 0 ) then - cindex = cindex + 1 - factor_een_naive = factor_een_naive + & - cord_vect(cindex, typenuc_arr(alpha)) * & - rij ** k * (ria ** l + rja ** l) * (ria * rja) ** rshift(p - k - l, 1) - !factor_een_naive = factor_een_naive + & - ! cord_vect(cindex, typenuc_arr(alpha)) * & - ! rij(i, j, k) * (ria(i, alpha, l) + rja(j, alpha, l)) & - ! * (ria(i, alpha, l) * rja(j, alpha, l)) ** rshift(p - k - l, 1) - end if - end do - end do - end do - end do - end do - end do + factor_een_deriv_e = 0.0d0 + + do p=2,ncord + do k=0,p-1 + if (k /= 0) then + lmax = p-k + else + lmax = p-k-2 + endif + + do l=0,lmax + if ( iand(p-k-l,1) == 1) then + cycle + endif + m = (p-k-l)/2 + + do a=1, nnuc + cn = cord_vect_lkp(l,k,p,typenuc_arr(a)) + do j=1, nelec + rjal = rescale_een_n(j,a,l) + rjam_cn = rescale_een_n(j,a,m) * cn + do ii=1,4 + drjal(ii) = rescale_een_n_deriv_e(ii,j,a,l) + drjam_cn(ii) = rescale_een_n_deriv_e(ii,j,a,m) * cn + enddo + factor_een_deriv_e(:,j) = 0.d0 + do i=1, nelec + rial = rescale_een_n(i,a,l) + riam = rescale_een_n(i,a,m) + rijk = rescale_een_e(i,j,k) + do ii=1,4 + drijk(ii) = rescale_een_e_deriv_e(ii,i,j,k) + enddo + l1 = 0.d0 + l2 = 0.d0 + x(1:3) = 0.d0 + x(4) = 2.d0 + do ii=1,4 + v1 = rijk * (rial+rjal) + v2 = rjam_cn * riam + + d1 = drijk(ii) * (rial+rjal) + rijk * (rial+drjal(ii)) + d2 = drjam_cn(ii) * riam + + l1 = l1 + drijk(ii) * (rial+drjal(ii)) + l2 = l2 + drjam_cn(ii) * riam + + factor_een_deriv_e(ii,j) = factor_een_deriv_e(ii,j) + & + v1 * d2 + d1 * v2 + x(ii) * (l1 + l2) + enddo + factor_een_deriv_e(ii,j) = factor_een_deriv_e(ii,j) + & + v1 * d2 + d1 * v2 + enddo + enddo + enddo + + enddo + enddo + enddo + factor_een_deriv_e = 0.5d0 * factor_een_deriv_e END_PROVIDER -!BEGIN_PROVIDER [double precision, factor_een_prog] -! implicit none -! BEGIN_DOC -! ! Electron-electron nucleus contribution to Jastrow factor in a naive way -! END_DOC -! integer :: alpha, i, j, p, k, l, lmax, m, cindex -! double precision :: ria, rja, rij, rij_inv -! double precision :: c, c_inv, t, x, y, z ! Placeholders for optimization -! -! PROVIDE cord_vect -! factor_een_prog = 0.0d0 -! -! do alpha = 1, nnuc -! do j = 2, nelec -! rja = rescale_een_n(j, alpha) -! do i = 1, j - 1 -! ria = rescale_een_n(i, alpha) -! rij = rescale_een_e(i, j) -! rij_inv = 1.0d0 / (rij * rij) -! c = ria * rja -! c_inv = 1.0d0 / c -! cindex = 0 -! do p = 2, ncord -! -! x = 1.0d0 -! do l = 1, p -! x = x * rij -! end do -! -! do k = p - 1, 0, -1 -! if ( k /= 0 ) then -! lmax = p - k -! else -! lmax = p - k - 2 -! end if -! -! t = 1.0d0 -! do l = 1, rshift(p - k, 1) -! t = t * c -! end do -! -! do l = lmax, iand(p - k, 1), -2 -! cindex = cindex + 1 -! factor_een_prog = factor_een_prog + & -! cord_vect(cindex, typenuc_arr(alpha)) * & -! x * (ria ** l + rja ** l) * t -! t = t * c_inv -! x = x * rij_inv -! end do -! -! end do -! end do -! end do -! end do -! end do -! -!END_PROVIDER - -!BEGIN_PROVIDER [double precision, rij, (nelec, nelec, ncord)] -!&BEGIN_PROVIDER [double precision, ria, (nelec, nnuc, ncord)] -!&BEGIN_PROVIDER [double precision, rja, (nelec, nnuc, ncord)] -! BEGIN_DOC -! ! Tables with powers -! END_DOC -! integer :: i, j, k, alpha -! double precision :: x, y, z -! -! rij(:, :, :) = 0.0d0 -! ria(:, :, :) = 0.0d0 -! rja(:, :, :) = 0.0d0 -! -! implicit none -! do alpha = 1, nnuc -! do j = 2, nelec -! z = 1.0d0 -! do k = 1, ncord -! rja(j, alpha, k) = z -! z = z * rescale_een_n(j, alpha) -! end do -! do i = 1, j - 1 -! y = 1.0d0 -! do k = 1, ncord -! ria(i, alpha, k) = y -! y = y * rescale_een_n(i, alpha) -! end do -! x = 1.0d0 -! do k = 1, ncord -! rij(i, j, k) = x -! x = x * rescale_een_e(i, j) -! end do -! end do -! end do -! end do -! -!END_PROVIDER - diff --git a/jastrow b/jastrow deleted file mode 100755 index c6b9e8e..0000000 Binary files a/jastrow and /dev/null differ diff --git a/orders.irp.f b/orders.irp.f index 444f53f..ab13944 100644 --- a/orders.irp.f +++ b/orders.irp.f @@ -69,3 +69,31 @@ BEGIN_PROVIDER [double precision, aord_vect, (naord + 1, typenuc)] close(fu) END_PROVIDER + +BEGIN_PROVIDER [ double precision, cord_vect_lkp, (0:ncord-1, 0:ncord-1, 2:ncord, typenuc) ] + implicit none + BEGIN_DOC + ! + END_DOC + integer :: alpha, l,k,p,lmax,cindex + + cord_vect_lkp = 0.d0 + cindex = 0 + do alpha=1,typenuc + do p = 2, ncord + do k = p-1, 0, -1 + if ( k /= 0 ) then + lmax = p - k + else + lmax = p - k - 2 + end if + do l = lmax, 0, -1 + if (iand(p-k-l,1) == 1) cycle + cindex = cindex + 1 + cord_vect_lkp(l,k,p,alpha) = cord_vect(cindex, alpha) + end do + end do + end do + end do + +END_PROVIDER diff --git a/rescale.irp.f b/rescale.irp.f index 56b0d66..2a55d30 100644 --- a/rescale.irp.f +++ b/rescale.irp.f @@ -42,30 +42,79 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ] enddo END_PROVIDER -BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec)] +BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec, 0:ncord)] implicit none BEGIN_DOC ! R = exp(-kappa r) for electron-electron for $J_{een}$ END_DOC - integer :: i, j + integer :: i, j, l + double precision :: kappa_l - do j = 1, nelec - do i = 1, nelec - rescale_een_e(i, j) = dexp(-kappa * elec_dist(i, j)) + do l=0,ncord + kappa_l = -dble(l) * kappa + do j = 1, nelec + do i = 1, nelec + rescale_een_e(i, j, l) = kappa_l * elec_dist(i, j) + enddo + enddo + enddo + ! More efficient to compute the exp of array than to do it in the loops + rescale_een_e = dexp(rescale_een_e) + + ! Later we use a formula looping on i and j=1->j-1. We need to set Rjj=0 to + ! enable looping of j=1,nelec do l=0,ncord + do l=0,ncord + do j=1,nelec + rescale_een_e(j, j, l) = 0.d0 enddo enddo END_PROVIDER -BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc)] +BEGIN_PROVIDER [double precision, rescale_een_n, (4, nelec, nnuc, 0:ncord)] implicit none BEGIN_DOC ! R = exp(-kappa r) for electron-electron for $J_{een}$ END_DOC - integer :: i, j + integer :: i, j, l + double precision :: kappa_l - do j = 1, nnuc - do i = 1, nelec - rescale_een_n(i, j) = dexp(-kappa * elnuc_dist(i, j)) + do l=0,ncord + kappa_l = - dble(l) * kappa + do j = 1, nnuc + do i = 1, nelec + rescale_een_n(i, j, l) = kappa_l * elnuc_dist(i, j) + enddo + enddo + enddo + rescale_een_n = dexp(rescale_een_n) +END_PROVIDER + +BEGIN_PROVIDER [double precision, rescale_een_n_deriv_e, (4,nelec, nnuc, 0:ncord)] + implicit none + BEGIN_DOC + ! R = exp(-kappa r) for electron-electron for $J_{een}$ + END_DOC + integer :: i, j, l + double precision :: kappa_l + + do l=0,ncord + kappa_l = - dble(l) * kappa + do j = 1, nnuc + do i = 1, nelec + do ii=1,4 + rescale_een_n_deriv_e(ii, i, j, l) = & + kappa_l * elnuc_dist_deriv_e(ii,i,j) + enddo + rescale_een_n_deriv_e(4, i, j, l) = rescale_een_n_deriv_e(4, i, j, l) + & + rescale_een_n_deriv_e(1, i, j, l) * rescale_een_n_deriv_e(1, i, j, l) + & + rescale_een_n_deriv_e(2, i, j, l) * rescale_een_n_deriv_e(2, i, j, l) + & + rescale_een_n_deriv_e(3, i, j, l) * rescale_een_n_deriv_e(3, i, j, l) + do ii=1,4 + rescale_een_n_deriv_e(ii, i, j, l) = & + rescale_een_n_deriv_e(ii,i,j, l) * rescale_een_n(i, j, l) + enddo + enddo enddo enddo END_PROVIDER +