diff --git a/deriv_num b/deriv_num deleted file mode 100755 index 0f56468..0000000 Binary files a/deriv_num and /dev/null differ diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f index 688ee9d..393ef20 100644 --- a/el_nuc_el.irp.f +++ b/el_nuc_el.irp.f @@ -1,111 +1,121 @@ BEGIN_PROVIDER [ double precision, factor_een ] - implicit none - BEGIN_DOC - ! ElectronE-electron-nuclei contribution to Jastrow factor - ! - ! 5436.20340250000 - END_DOC - integer :: i, j, a, p, k, l, lmax, m, n - double precision :: cn, accu2, accu + implicit none + BEGIN_DOC + ! Electron -electron-nuclei contribution to Jastrow factor + ! + ! 5436.20340250000 + END_DOC + integer :: i, j, a, p, k, l, lmax, m, n + double precision :: cn, accu2, accu -! factor_een = factor_een_blas -! return +! double precision :: ria_tmp(nelec,dim_cord_vect,nnuc) +! double precision :: rja_tmp(nelec,dim_cord_vect,nnuc) +! +! do a = 1, nnuc +! do n = 1, dim_cord_vect +! +! l = lkpm_of_cindex(1,n) +! k = lkpm_of_cindex(2,n) +! p = lkpm_of_cindex(3,n) +! m = lkpm_of_cindex(4,n) +! +! do i = 1, nelec +! ria_tmp(i,n,a) = rescale_een_n(i,a,m) +! rja_tmp(i,n,a) = rescale_een_n(i,a,m+l) +! enddo +! enddo +! +! enddo - factor_een = 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) cycle + ! factor_een = factor_een_blas + ! return - m = (p - k - l) / 2 + factor_een = 0.0d0 - do a = 1, nnuc - accu2 = 0.d0 - cn = cord_vect_lkp(l, k, p, typenuc_arr(a)) - do j = 1, nelec - accu = 0.d0 - do i = 1, nelec - accu = accu + & - rescale_een_e(i,j,k) * & - rescale_een_n(i,a,m) - enddo - accu2 = accu2 + accu*rescale_een_n(j,a,m+l) + do n = 1, dim_cord_vect + + l = lkpm_of_cindex(1,n) + k = lkpm_of_cindex(2,n) + p = lkpm_of_cindex(3,n) + m = lkpm_of_cindex(4,n) + + do a = 1, nnuc + accu2 = 0.d0 + cn = cord_vect_full(n, a) + do j = 1, nelec + accu = 0.d0 + do i = 1, nelec + accu = accu + & + rescale_een_e(i,j,k) * & + rescale_een_n(i,a,m) enddo - factor_een = factor_een + accu2 * cn + accu2 = accu2 + accu*rescale_een_n(j,a,m+l) enddo - + factor_een = factor_een + accu2 * cn enddo + enddo - enddo END_PROVIDER BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ] implicit none + BEGIN_DOC +! Derivative of the Jeen +! 35533.115255 + END_DOC integer :: i, j, a, p, k, l, lmax, m, n double precision :: cn, accu, accu2, daccu(1:4), daccu2(1:4) - - factor_een_deriv_e(1:4,1:nelec) = factor_een_deriv_e_blas(1:4,1:nelec) - return - - factor_een_deriv_e(1:4,1:nelec) = 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) cycle - - m = (p - k - l) / 2 - - do a = 1, nnuc - cn = cord_vect_lkp(l, k, p, typenuc_arr(a)) - do j = 1, nelec - accu=0.d0 - accu2 = 0.d0 - daccu (1:4) = 0.d0 - daccu2(1:4) = 0.d0 - do i = 1, nelec - accu = accu + & - rescale_een_e(i,j,k) * & - rescale_een_n(i,a,m) - accu2 = accu2 + & - rescale_een_e(i,j,k) * & - rescale_een_n(i,a,m+l) - daccu(1:4) = daccu(1:4) + & - rescale_een_e_deriv_e_t(1:4,i,j,k) * & - rescale_een_n(i,a,m) - daccu2(1:4) = daccu2(1:4) + & - rescale_een_e_deriv_e_t(1:4,i,j,k) * & - rescale_een_n(i,a,m+l) - - enddo - factor_een_deriv_e(1:4,j) = factor_een_deriv_e(1:4,j) + & - (accu * rescale_een_n_deriv_e(1:4,j,a,m+l) + daccu(1:4) * rescale_een_n(j,a,m+l) + & - daccu2(1:4)* rescale_een_n(j,a,m) + accu2*rescale_een_n_deriv_e(1:4,j,a,m)) * cn - - factor_een_deriv_e(4,j) = factor_een_deriv_e(4,j) + 2.d0*( & - daccu (1) * rescale_een_n_deriv_e(1,j,a,m+l) + & - daccu (2) * rescale_een_n_deriv_e(2,j,a,m+l) + & - daccu (3) * rescale_een_n_deriv_e(3,j,a,m+l) + & - daccu2(1) * rescale_een_n_deriv_e(1,j,a,m ) + & - daccu2(2) * rescale_een_n_deriv_e(2,j,a,m ) + & - daccu2(3) * rescale_een_n_deriv_e(3,j,a,m ) )*cn - enddo - enddo - enddo - enddo - enddo + +! factor_een_deriv_e(1:4,1:nelec) = factor_een_deriv_e_blas(1:4,1:nelec) +! return + + factor_een_deriv_e(1:4,1:nelec) = 0.0d0 + + do n = 1, dim_cord_vect + + l = lkpm_of_cindex(1,n) + k = lkpm_of_cindex(2,n) + p = lkpm_of_cindex(3,n) + m = lkpm_of_cindex(4,n) + + do a = 1, nnuc + cn = cord_vect_full(n, a) + do j = 1, nelec + accu=0.d0 + accu2 = 0.d0 + daccu (1:4) = 0.d0 + daccu2(1:4) = 0.d0 + do i = 1, nelec + accu = accu + & + rescale_een_e(i,j,k) * & + rescale_een_n(i,a,m) + accu2 = accu2 + & + rescale_een_e(i,j,k) * & + rescale_een_n(i,a,m+l) + daccu(1:4) = daccu(1:4) + & + rescale_een_e_deriv_e_t(1:4,i,j,k) * & + rescale_een_n(i,a,m) + daccu2(1:4) = daccu2(1:4) + & + rescale_een_e_deriv_e_t(1:4,i,j,k) * & + rescale_een_n(i,a,m+l) + + enddo + factor_een_deriv_e(1:4,j) = factor_een_deriv_e(1:4,j) + & + (accu * rescale_een_n_deriv_e(1:4,j,a,m+l) + daccu(1:4) * rescale_een_n(j,a,m+l) +& + daccu2(1:4)* rescale_een_n(j,a,m) + accu2*rescale_een_n_deriv_e(1:4,j,a,m)) * cn + + factor_een_deriv_e(4,j) = factor_een_deriv_e(4,j) + 2.d0*( & + daccu (1) * rescale_een_n_deriv_e(1,j,a,m+l) + & + daccu (2) * rescale_een_n_deriv_e(2,j,a,m+l) + & + daccu (3) * rescale_een_n_deriv_e(3,j,a,m+l) + & + daccu2(1) * rescale_een_n_deriv_e(1,j,a,m ) + & + daccu2(2) * rescale_een_n_deriv_e(2,j,a,m ) + & + daccu2(3) * rescale_een_n_deriv_e(3,j,a,m ) )*cn + enddo + enddo + enddo END_PROVIDER diff --git a/el_nuc_el_blas.irp.f b/el_nuc_el_blas.irp.f index 9330ff8..8814fcf 100644 --- a/el_nuc_el_blas.irp.f +++ b/el_nuc_el_blas.irp.f @@ -32,51 +32,43 @@ enddo - do p = 2, ncord - do k = 0, p - 1 - m = p-k - if (k > 0) then - lmax = m - else - lmax = m - 2 - endif + do n = 1, dim_cord_vect - n = shiftr(m,1) - do l = iand(m, 1), lmax, 2 + l = lkpm_of_cindex(1,n) + k = lkpm_of_cindex(2,n) + p = lkpm_of_cindex(3,n) + m = lkpm_of_cindex(4,n) - do a = 1, nnuc - cn(a) = cord_vect_lkp(l, k, p, typenuc_arr(a)) - enddo + do a = 1, nnuc + cn(a) = cord_vect_full(n, a) + enddo - do a = 1, nnuc - accu = 0.d0 + do a = 1, nnuc + accu = 0.d0 - do j=1,nelec - accu = accu + rescale_een_n(j,a,n) * tmp_c(j,a,n+l,k) + do j=1,nelec + accu = accu + rescale_een_n(j,a,m) * tmp_c(j,a,m+l,k) - factor_een_deriv_e_blas(1:4,j) = factor_een_deriv_e_blas(1:4,j) + (& - tmp_c(j,a,n,k) * rescale_een_n_deriv_e(1:4,j,a,n+l) + & - dtmp_c(1:4,j,a,n,k) * rescale_een_n(j,a,n+l) + & - dtmp_c(1:4,j,a,n+l,k) * rescale_een_n(j,a,n) + & - tmp_c(j,a,n+l,k)*rescale_een_n_deriv_e(1:4,j,a,n) & - ) * cn(a) + factor_een_deriv_e_blas(1:4,j) = factor_een_deriv_e_blas(1:4,j) + (& + tmp_c(j,a,m,k) * rescale_een_n_deriv_e(1:4,j,a,m+l) + & + dtmp_c(1:4,j,a,m,k) * rescale_een_n(j,a,m+l) + & + dtmp_c(1:4,j,a,m+l,k) * rescale_een_n(j,a,m) + & + tmp_c(j,a,m+l,k)*rescale_een_n_deriv_e(1:4,j,a,m) & + ) * cn(a) - factor_een_deriv_e_blas(4,j) = factor_een_deriv_e_blas(4,j) + (& - dtmp_c(1,j,a,n ,k) * rescale_een_n_deriv_e(1,j,a,n+l) +& - dtmp_c(2,j,a,n ,k) * rescale_een_n_deriv_e(2,j,a,n+l) +& - dtmp_c(3,j,a,n ,k) * rescale_een_n_deriv_e(3,j,a,n+l) +& - dtmp_c(1,j,a,n+l,k) * rescale_een_n_deriv_e(1,j,a,n ) +& - dtmp_c(2,j,a,n+l,k) * rescale_een_n_deriv_e(2,j,a,n ) +& - dtmp_c(3,j,a,n+l,k) * rescale_een_n_deriv_e(3,j,a,n )& - )*cn(a)*2.d0 - - enddo - factor_een_blas = factor_een_blas + accu * cn(a) - - enddo - n = n-1 + factor_een_deriv_e_blas(4,j) = factor_een_deriv_e_blas(4,j) + (& + dtmp_c(1,j,a,m ,k) * rescale_een_n_deriv_e(1,j,a,m+l) + & + dtmp_c(2,j,a,m ,k) * rescale_een_n_deriv_e(2,j,a,m+l) + & + dtmp_c(3,j,a,m ,k) * rescale_een_n_deriv_e(3,j,a,m+l) + & + dtmp_c(1,j,a,m+l,k) * rescale_een_n_deriv_e(1,j,a,m ) + & + dtmp_c(2,j,a,m+l,k) * rescale_een_n_deriv_e(2,j,a,m ) + & + dtmp_c(3,j,a,m+l,k) * rescale_een_n_deriv_e(3,j,a,m ) & + )*cn(a)*2.d0 enddo + factor_een_blas = factor_een_blas + accu * cn(a) + enddo enddo + END_PROVIDER diff --git a/jastrow b/jastrow deleted file mode 100755 index b1b7d3b..0000000 Binary files a/jastrow and /dev/null differ diff --git a/orders.irp.f b/orders.irp.f index 922f3be..c0bbe68 100644 --- a/orders.irp.f +++ b/orders.irp.f @@ -21,7 +21,7 @@ BEGIN_PROVIDER [integer, ncord] END_DOC ncord = 5 END_PROVIDER - + BEGIN_PROVIDER [integer, dim_cord_vect] implicit none BEGIN_DOC @@ -32,20 +32,50 @@ BEGIN_PROVIDER [integer, dim_cord_vect] dim_cord_vect = 0 do p = 2, ncord - do k = 0, p - 1 - if ( k /= 0 ) then - lmax = p - k - else - lmax = p - k - 2 - end if - do l = iand(p - k, 1), lmax, 2 - dim_cord_vect = dim_cord_vect + 1 + 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 + dim_cord_vect = dim_cord_vect + 1 end do end do end do END_PROVIDER - + + +BEGIN_PROVIDER [ integer, lkpm_of_cindex, (4,dim_cord_vect) ] + implicit none + BEGIN_DOC +! Transform l,k,p into a consecutive index + END_DOC + integer :: p,k,l,lmax,m + integer :: kk + kk=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) == 1) cycle + m = (p - k - l) / 2 + kk = kk+1 + lkpm_of_cindex(1,kk) = l + lkpm_of_cindex(2,kk) = k + lkpm_of_cindex(3,kk) = p + lkpm_of_cindex(4,kk) = m + enddo + enddo + enddo + +END_PROVIDER BEGIN_PROVIDER [double precision, aord_vect, (naord + 1, typenuc)] &BEGIN_PROVIDER [double precision, bord_vect, (nbord + 1)] @@ -59,17 +89,29 @@ BEGIN_PROVIDER [double precision, aord_vect, (naord + 1, typenuc)] PROVIDE ncord character(len=*), parameter :: FILE_NAME = "jast_coeffs.txt" integer :: i, fu, rc - + open(action='read', file=FILE_NAME, iostat=rc, newunit=fu) - + read(fu, *) aord_vect read(fu, *) bord_vect read(fu, *) cord_vect - + close(fu) END_PROVIDER +BEGIN_PROVIDER [ double precision, cord_vect_full, (dim_cord_vect, nnuc) ] + implicit none + BEGIN_DOC + ! cord_vect for all atoms + END_DOC + integer :: a + do a=1,nnuc + cord_vect_full(1:dim_cord_vect,a) = cord_vect(1:dim_cord_vect,typenuc_arr(a)) + enddo + +END_PROVIDER + BEGIN_PROVIDER [ double precision, cord_vect_lkp, (0:ncord-1, 0:ncord-1, 2:ncord, typenuc) ] implicit none BEGIN_DOC @@ -78,8 +120,8 @@ BEGIN_PROVIDER [ double precision, cord_vect_lkp, (0:ncord-1, 0:ncord-1, 2:ncord integer :: alpha, l, k, p, lmax, cindex cord_vect_lkp = 0.0d0 - cindex = 0 do alpha = 1, typenuc + cindex = 0 do p = 2, ncord do k = p - 1, 0, -1 if ( k /= 0 ) then