From 8ac9bcb963a2c3b7c81d8e0b376d6609408fa2d9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 14 May 2021 14:03:34 +0200 Subject: [PATCH 1/4] New scheme for pseudo explosion --- src/PROPERTIES/properties_energy.irp.f | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 56a91e6..9c63d8c 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -252,12 +252,10 @@ BEGIN_PROVIDER [ double precision, E_loc ] enddo ! Avoid divergence of E_loc and population explosion - if (do_pseudo) then + if (do_pseudo .and. qmc_method == 'DMC' ) then double precision :: delta_e -! delta_e = E_loc-E_ref -! E_loc = E_ref + erf(delta_e*time_step_sq)/time_step_sq - E_loc = max(2.d0*E_ref, E_loc) -! continue + delta_e = E_loc-E_ref + E_loc = E_ref + delta_e * dexp(-dabs(delta_e)*time_step_sq) endif E_loc_min = min(E_loc,E_loc_min) E_loc_max = max(E_loc,E_loc_max) From 216ddc4351f16ab88992aef8c10667acf36539eb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 May 2021 16:05:38 +0200 Subject: [PATCH 2/4] Fixed t_DMC --- src/PROPERTIES/properties_energy.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 9c63d8c..84554a2 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -252,7 +252,7 @@ BEGIN_PROVIDER [ double precision, E_loc ] enddo ! Avoid divergence of E_loc and population explosion - if (do_pseudo .and. qmc_method == 'DMC' ) then + if (do_pseudo .and. (qmc_method == t_DMC) ) then double precision :: delta_e delta_e = E_loc-E_ref E_loc = E_ref + delta_e * dexp(-dabs(delta_e)*time_step_sq) From df2559cf1f88bac42aeac1f5c29f967fde16cc5c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 13:57:11 +0200 Subject: [PATCH 3/4] Fix DGEMM --- configure.sh | 2 +- src/det.irp.f | 437 ++++++++++++++++++++++++-------------------------- 2 files changed, 213 insertions(+), 226 deletions(-) diff --git a/configure.sh b/configure.sh index 535eea0..ebfd996 100755 --- a/configure.sh +++ b/configure.sh @@ -44,6 +44,6 @@ export C_INCLUDE_PATH="\${QMCCHEM_PATH}/include:\${C_INCLUDE_PATH}" #export QMCCHEM_NIC=ib0 source \${QMCCHEM_PATH}/irpf90/bin/irpman #source \${QMCCHEM_PATH}/EZFIO/Bash/ezfio.sh -eval $(opam env) +eval \$(opam env) EOF diff --git a/src/det.irp.f b/src/det.irp.f index 5be8df9..e3dd661 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -1,32 +1,32 @@ BEGIN_PROVIDER [ integer, det_i ] - + BEGIN_DOC ! Current running alpha determinant END_DOC det_i=det_alpha_order(1) - + END_PROVIDER BEGIN_PROVIDER [ integer, det_j ] - + BEGIN_DOC ! Current running beta determinant END_DOC det_j=det_beta_order(1) - + END_PROVIDER subroutine det_update(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(LDS) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + if (d == 0.d0) then return endif @@ -193,15 +193,15 @@ end subroutine det_update2(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(2) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,2) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,2) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + S(1,l) = m(1) S(2,l) = m(2) S_inv(1,1) = S(1,1) @@ -209,37 +209,37 @@ subroutine det_update2(n,LDS,m,l,S,S_inv,d) S_inv(2,1) = S(1,2) S_inv(2,2) = S(2,2) call invert2(S_inv,LDS,n,d) - + end subroutine det_update1(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(1) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,1) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,1) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + S(1,l) = m(1) S_inv(1,1) = S(1,1) call invert1(S_inv,LDS,n,d) - + end subroutine det_update3(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(3) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,3) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,3) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + integer :: i do i=1,3 S(i,l) = m(i) @@ -249,22 +249,22 @@ subroutine det_update3(n,LDS,m,l,S,S_inv,d) S_inv(2,i) = S(i,2) S_inv(3,i) = S(i,3) enddo - + call invert3(S_inv,LDS,n,d) - + end subroutine det_update4(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(4) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,4) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,4) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u(4), z(4), w(4), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u integer :: i,j @@ -276,7 +276,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) z(2) = S_inv(1,2)*u(1) + S_inv(2,2)*u(2) + S_inv(3,2)*u(3) + S_inv(4,2)*u(4) z(3) = S_inv(1,3)*u(1) + S_inv(2,3)*u(2) + S_inv(3,3)*u(3) + S_inv(4,3)*u(4) z(4) = S_inv(1,4)*u(1) + S_inv(2,4)*u(2) + S_inv(3,4)*u(3) + S_inv(4,4)*u(4) - + d_inv = 1.d0/d d = d+z(l) lambda = d_inv*d @@ -284,7 +284,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do i=1,4 w(i) = S_inv(i,l)*d_inv @@ -297,22 +297,22 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) S_inv(j,i) = S_inv(j,i)*lambda -z(i)*w(j) enddo enddo - + end BEGIN_TEMPLATE ! Version for mod(n,4) = 0 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -321,20 +321,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (LDS >= $n) integer :: i,j double precision :: zj, zj1, zj2, zj3 - + !DIR$ NOPREFETCH - !DIR$ SIMD NOVECREMAINDER + !OMP$ SIMD NOVECREMAINDER do i=1,$n u(i) = m(i) - S(i,l) enddo - + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-1,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo d_inv = 1.d0/d @@ -344,7 +344,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n,4 zj = 0.d0 @@ -353,26 +353,26 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) zj2 = zj2 + S_inv(i,j+2)*u(i) zj3 = zj3 + S_inv(i,j+3)*u(i) enddo - z(j ) = zj - z(j+1) = zj1 - z(j+2) = zj2 + z(j ) = zj + z(j+1) = zj1 + z(j+2) = zj2 z(j+3) = zj3 enddo !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n,4 zj = z(i ) zj1 = z(i+1) @@ -380,7 +380,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 @@ -435,15 +435,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 1 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -452,20 +452,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (LDS >= $n) integer :: i,j double precision :: zj, zj1, zj2, zj3 - + do i=1,$n u(i) = m(i) - S(i,l) enddo zj = 0.d0 !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) + !OMP$ SIMD REDUCTION(+:zj) do i=1,$n-1,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo zj = zj + S_inv($n,l)*u($n) - + d_inv = 1.d0/d d = d+zj lambda = d*d_inv @@ -473,7 +473,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-1,4 zj = 0.d0 @@ -482,7 +482,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n-1 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -498,19 +498,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-1 zj = zj + S_inv(i,$n)*u(i) enddo z($n) = zj + S_inv($n,$n)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-1,4 zj = z(i ) zj1 = z(i+1) @@ -518,14 +518,14 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n-1 S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 enddo - S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj + S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*zj1 S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*zj2 S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*zj3 @@ -534,7 +534,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = z($n) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER do i=1,$n S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*zj enddo @@ -588,15 +588,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 2 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -604,24 +604,24 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) !DIR$ ASSUME (LDS >= $n) integer :: i,j - + double precision :: zj, zj1, zj2, zj3 !DIR$ NOPREFETCH - !DIR$ SIMD NOVECREMAINDER + !OMP$ SIMD NOVECREMAINDER do i=1,$n u(i) = m(i) - S(i,l) enddo - + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-2,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo i=$n-1 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) d_inv = 1.d0/d d = d+zj @@ -630,7 +630,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-2,4 zj = 0.d0 @@ -638,7 +638,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 zj3 = 0.d0 !DIR$ VECTOR ALIGNED - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n-2 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -660,7 +660,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj1 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER do i=1,$n-2 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -671,19 +671,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-2,4 zj = z(i) zj1 = z(i+1) zj2 = z(i+2) zj3 = z(i+3) !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n-2 S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -704,7 +704,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = z(i) zj1= z(i+1) !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1) do j=1,$n-2 S_inv(j,i ) = S_inv(j,i )*lambda -zj*w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -761,15 +761,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 3 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -780,18 +780,18 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) double precision :: zj, zj1, zj2, zj3 - !DIR$ SIMD + !OMP$ SIMD do i=1,$n u(i) = m(i) - S(i,l) enddo - - zj = 0.d0 + + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-3,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo i=$n-2 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + S_inv(i+2,l)*u(i+2) @@ -804,7 +804,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-3,4 zj = 0.d0 @@ -812,7 +812,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 zj3 = 0.d0 !DIR$ VECTOR ALIGNED - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) do i=1,$n-3 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -839,7 +839,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2) do i=1,$n-3 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -856,12 +856,12 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-3,4 zj = z(i) zj1 = z(i+1) @@ -869,9 +869,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) do j=1,$n-3 - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 @@ -896,7 +896,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = z(i+2) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2) do j=1,$n S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 @@ -950,15 +950,15 @@ END_TEMPLATE subroutine det_update_general(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(LDS) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: lambda, d_inv double precision :: u(3840), z(3840), w(3840) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u @@ -976,19 +976,19 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) do i=1,n u(i) = m(i) - S(i,l) enddo - + zl = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zl) + !OMP$ SIMD REDUCTION(+:zl) do i=1,n zl = zl + S_inv(i,l)*u(i) enddo - d_inv = 1.d0/d + d_inv = 1.d0/d d = d+zl lambda = d*d_inv - + if ( dabs(lambda) < 1.d-3 ) then d = 0.d0 return @@ -1006,7 +1006,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) do i=1,n zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -1023,7 +1023,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) + !OMP$ SIMD REDUCTION(+:zj) do i=1,n zj = zj + S_inv(i,j)*u(i) enddo @@ -1031,14 +1031,14 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) enddo !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) @@ -1051,7 +1051,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) do j=1,n S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -1064,7 +1064,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj = z(i) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj) do j=1,n S_inv(j,i) = S_inv(j,i)*lambda -zj*w(j) enddo @@ -1146,9 +1146,9 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, det_alpha_value_curr ] &BEGIN_PROVIDER [ real, slater_matrix_alpha, (elec_alpha_num_8,elec_alpha_num) ] &BEGIN_PROVIDER [ double precision, slater_matrix_alpha_inv_det, (elec_alpha_num_8,elec_alpha_num) ] - + implicit none - + BEGIN_DOC ! det_alpha_value_curr : Value of the current alpha determinant ! @@ -1159,7 +1159,7 @@ END_PROVIDER ! ! slater_matrix_alpha_inv_det: Inverse of the alpha Slater matrix x determinant END_DOC - + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num), n_to_do_old, n_to_do @@ -1172,7 +1172,7 @@ END_PROVIDER !DIR$ VECTOR ALIGNED slater_matrix_alpha_inv_det = 0.d0 endif - + PROVIDE mo_value if (det_i /= det_alpha_order(1) ) then ! if (det_i == -1 ) then @@ -1187,7 +1187,7 @@ END_PROVIDER enddo ddet = 0.d0 - + ! n_to_do < 2*elec_alpha_num if (n_to_do < shiftl(elec_alpha_num,1)) then @@ -1219,13 +1219,13 @@ END_PROVIDER enddo endif - + else ddet = 0.d0 endif - + ! Avoid NaN if (ddet /= 0.d0) then continue @@ -1247,11 +1247,11 @@ END_PROVIDER enddo enddo call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet) - + endif ASSERT (ddet /= 0.d0) - - det_alpha_value_curr = ddet + + det_alpha_value_curr = ddet END_PROVIDER @@ -1268,24 +1268,24 @@ END_PROVIDER ! ! slater_matrix_beta_inv_det : Inverse of the beta Slater matrix x determinant END_DOC - + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num-mo_closed_num), n_to_do_old, n_to_do - + integer, save :: ifirst if (elec_beta_num == 0) then det_beta_value_curr = 0.d0 return endif - + if (ifirst == 0) then ifirst = 1 slater_matrix_beta = 0. slater_matrix_beta_inv_det = 0.d0 endif PROVIDE mo_value - + if (det_j /= det_beta_order(1)) then ! if (det_j == -1) then @@ -1297,7 +1297,7 @@ END_PROVIDER to_do(n_to_do) = k endif enddo - + ddet = 0.d0 if (n_to_do < shiftl(elec_beta_num,1)) then @@ -1329,13 +1329,13 @@ END_PROVIDER enddo endif - + else - + ddet = 0.d0 endif - + ! Avoid NaN if (ddet /= 0.d0) then continue @@ -1359,9 +1359,9 @@ END_PROVIDER call invert(slater_matrix_beta_inv_det,elec_beta_num_8,elec_beta_num,ddet) endif ASSERT (ddet /= 0.d0) - - det_beta_value_curr = ddet - + + det_beta_value_curr = ddet + END_PROVIDER BEGIN_PROVIDER [ integer, det_alpha_num_pseudo ] @@ -1383,15 +1383,15 @@ END_PROVIDER BEGIN_PROVIDER [ double precision , det_alpha_value, (det_alpha_num_8) ] &BEGIN_PROVIDER [ double precision , det_alpha_grad_lapl, (4,elec_alpha_num,det_alpha_num) ] &BEGIN_PROVIDER [ double precision , det_alpha_pseudo, (elec_alpha_num_8,det_alpha_num_pseudo) ] - + implicit none - + BEGIN_DOC ! Values of the alpha determinants ! Gradients of the alpha determinants ! Laplacians of the alpha determinants END_DOC - + integer :: j,i,k integer, save :: ifirst = 0 if (ifirst == 0) then @@ -1400,73 +1400,73 @@ END_PROVIDER det_alpha_grad_lapl = 0.d0 det_alpha_pseudo = 0.d0 endif - - + + do j=1,det_alpha_num - + det_i = det_alpha_order(j) if (j > 1) then TOUCH det_i endif - + det_alpha_value(det_i) = det_alpha_value_curr det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:) if (do_pseudo) then - det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) + det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) endif - + enddo - + det_i = det_alpha_order(1) SOFT_TOUCH det_i - + END_PROVIDER BEGIN_PROVIDER [ double precision, det_beta_value, (det_beta_num_8) ] &BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_alpha_num+1:elec_num,det_beta_num) ] &BEGIN_PROVIDER [ double precision, det_beta_pseudo, (elec_alpha_num+1:elec_num,det_beta_num_pseudo) ] - - + + implicit none - + BEGIN_DOC ! Values of the beta determinants ! Gradients of the beta determinants ! Laplacians of the beta determinants END_DOC - + integer :: j,i,k integer, save :: ifirst = 0 if (elec_beta_num == 0) then det_beta_value = 1.d0 return endif - + if (ifirst == 0) then ifirst = 1 det_beta_value = 0.d0 det_beta_grad_lapl = 0.d0 det_beta_pseudo = 0.d0 endif - + do j=1,det_beta_num - + det_j = det_beta_order(j) if (j > 1) then TOUCH det_j endif - + det_beta_value(det_j) = det_beta_value_curr det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:) if (do_pseudo) then - det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) + det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) endif enddo - + det_j = det_beta_order(1) SOFT_TOUCH det_j - + END_PROVIDER @@ -1474,34 +1474,34 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, psidet_inv ] &BEGIN_PROVIDER [ double precision, psidet_grad_lapl, (4,elec_num_8) ] &BEGIN_PROVIDER [ double precision, pseudo_non_local, (elec_num) ] - + implicit none BEGIN_DOC ! Value of the determinantal part of the wave function - + ! Gradient of the determinantal part of the wave function - + ! Laplacian of determinantal part of the wave function ! Non-local component of the pseudopotentials ! Regularized 1/psi = 1/(psi + 1/psi) END_DOC - + integer, save :: ifirst=0 if (ifirst == 0) then ifirst = 1 psidet_grad_lapl = 0.d0 endif - + double precision :: CDb(det_alpha_num_8) double precision :: CDb_i double precision :: DaC(det_beta_num_8) !DIR$ ATTRIBUTES ALIGN : 32 :: DaC,CDb ! C x D_beta - ! D_alpha^t x C - ! D_alpha^t x (C x D_beta) + ! D_alpha^t x C + ! D_alpha^t x (C x D_beta) integer :: i,j,k, l integer :: i1,i2,i3,i4,det_num4 @@ -1560,11 +1560,11 @@ END_PROVIDER enddo else - ! DaC(det_beta_num) = Trans( det_coef_matrix_dense(det_alpha_num,det_beta_num) ) + ! DaC(det_beta_num) = Trans( det_coef_matrix_dense(det_alpha_num,det_beta_num) ) ! @ det_alpha_value(det_alpha_num) call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1) - ! CDb(det_alpha_num) = det_coef_matrix_dense(det_alpha_num,det_beta_num) + ! CDb(det_alpha_num) = det_coef_matrix_dense(det_alpha_num,det_beta_num) ! @ det_beta_value(det_beta_num) call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1) @@ -1592,20 +1592,20 @@ END_PROVIDER ! Gradients ! --------- if(det_num .eq. 1) then - do i = 1, elec_alpha_num + do i = 1, elec_alpha_num psidet_grad_lapl(1,i) = det_alpha_grad_lapl_curr(1,i) * det_beta_value_curr psidet_grad_lapl(2,i) = det_alpha_grad_lapl_curr(2,i) * det_beta_value_curr psidet_grad_lapl(3,i) = det_alpha_grad_lapl_curr(3,i) * det_beta_value_curr psidet_grad_lapl(4,i) = det_alpha_grad_lapl_curr(4,i) * det_beta_value_curr enddo - do i = elec_alpha_num+1, elec_num + do i = elec_alpha_num+1, elec_num psidet_grad_lapl(1,i) = det_beta_grad_lapl_curr(1,i) * det_alpha_value_curr psidet_grad_lapl(2,i) = det_beta_grad_lapl_curr(2,i) * det_alpha_value_curr psidet_grad_lapl(3,i) = det_beta_grad_lapl_curr(3,i) * det_alpha_value_curr psidet_grad_lapl(4,i) = det_beta_grad_lapl_curr(4,i) * det_alpha_value_curr enddo else - ! psidet_grad_lapl(4,elec_alpha_num) = + ! psidet_grad_lapl(4,elec_alpha_num) = ! det_alpha_grad_lapl(4,elec_alpha_num,det_alpha_num) @ CDb(det_alpha_num) call dgemv('N',elec_alpha_num*4,det_alpha_num,1.d0, & det_alpha_grad_lapl, & @@ -1703,30 +1703,20 @@ END_PROVIDER , 1.d0, psi_svd_beta_unique(:,:,1), size(psi_svd_beta_unique,1), det_beta_value, 1 & , 0.d0, det_beta_value_SVD_unique, 1) - do mm = 1, 4 - call dgemm('N', 'N', elec_alpha_num, n_svd_coefs_unique, det_alpha_num, 1.d0 & - , det_alpha_grad_lapl(mm,:,:), size(det_alpha_grad_lapl,2) & - , psi_svd_alpha_unique(:,:,1), size(psi_svd_alpha_unique,1) & - , 0.d0, det_alpha_grad_lapl_SVD_unique(mm,:,:), size(det_alpha_grad_lapl_SVD_unique,2) ) - if (elec_beta_num /= 0) then - call dgemm('N', 'N', elec_beta_num, n_svd_coefs_unique, det_beta_num, 1.d0 & - , det_beta_grad_lapl(mm,:,:), size(det_beta_grad_lapl,2) & - , psi_svd_beta_unique(:,:,1), size(psi_svd_beta_unique,1) & - , 0.d0, det_beta_grad_lapl_SVD_unique(mm,:,:), size(det_beta_grad_lapl_SVD_unique,2) ) - endif - enddo + call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs_unique, det_alpha_num, 1.d0 & + , det_alpha_grad_lapl, 4*size(det_alpha_grad_lapl,2) & + , psi_svd_alpha_unique, size(psi_svd_alpha_unique,1) & + , 0.d0, det_alpha_grad_lapl_SVD_unique, 4*size(det_alpha_grad_lapl_SVD_unique,2) ) + if (elec_beta_num /= 0) then + call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs_unique, det_beta_num, 1.d0 & + , det_beta_grad_lapl, 4*size(det_beta_grad_lapl,2) & + , psi_svd_beta_unique, size(psi_svd_beta_unique,1) & + , 0.d0, det_beta_grad_lapl_SVD_unique, 4*size(det_beta_grad_lapl_SVD_unique,2) ) + endif END_PROVIDER - - - - - - - - BEGIN_PROVIDER [ logical, utilise_SVD ] &BEGIN_PROVIDER [ integer, n_svd_coefs ] implicit none @@ -1760,12 +1750,12 @@ END_PROVIDER ! !!! ! truncated SVD END_DOC - integer :: l + integer :: l do l = 1, n_svd_coefs psi_svd_coefs(l,1) = psi_svd_coefs_unique(l,1) psi_svd_alpha(:,l,1) = psi_svd_alpha_unique(:,l,1) psi_svd_beta (:,l,1) = psi_svd_beta_unique (:,l,1) - enddo + enddo END_PROVIDER @@ -1853,14 +1843,15 @@ END_PROVIDER ! !!! if (do_pseudo) then ! det_alpha_pseudo_SVD = det_alpha_pseudo @ psi_svd_alpha * psidet_inv_SVD - call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, psidet_inv_SVD & - , det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) & - , 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) ) + call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs , det_alpha_num, psidet_inv_SVD & + , det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(1,1,1) & + , size(psi_svd_alpha,1), 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) ) ! !!! if (elec_beta_num /= 0) then ! det_beta_pseudo_SVD = det_beta_pseudo @ psi_svd_beta * psidet_inv_SVD - call dgemm('N', 'N', elec_beta_num, n_svd_coefs , det_beta_num, psidet_inv_SVD & - , det_beta_pseudo(elec_alpha_num+1,1), size(det_beta_pseudo,1), psi_svd_beta(:,:,1), size(psi_svd_beta,1) & + call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs , det_beta_num, psidet_inv_SVD & + , det_beta_pseudo(elec_alpha_num+1,1), size(det_beta_pseudo,1) & + , psi_svd_beta(1,1,1), size(psi_svd_beta,1) & , 0.d0, det_beta_pseudo_SVD(elec_alpha_num+1,1), size(det_beta_pseudo_SVD,1) ) endif endif @@ -1890,15 +1881,15 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_SVD, (4, elec_alpha_num , n_svd_coefs) ] &BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_SVD , (4, elec_alpha_num+1:elec_num, n_svd_coefs) ] - - implicit none + + implicit none BEGIN_DOC ! !!! ! det_alpha_grad_lapl_SVD_{k} = sum_i U_{i,k} det_alpha_grad_lapl_{i} ! det_beta_grad_lapl_SVD_{k } = sum_j V_{j,k} det_beta_grad_lapl_{j } ! !!! END_DOC - + integer :: mm integer, save :: ifirst = 0 if (ifirst == 0) then @@ -1909,21 +1900,17 @@ END_PROVIDER ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ ! !!! - do mm = 1, 4 - ! !!! - call dgemm('N', 'N', elec_alpha_num, n_svd_coefs, det_alpha_num, 1.d0 & - , det_alpha_grad_lapl(mm,:,:), size(det_alpha_grad_lapl,2) & - , psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) & - , 0.d0, det_alpha_grad_lapl_SVD(mm,:,:), size(det_alpha_grad_lapl_SVD,2) ) - ! !!! - if (elec_beta_num /= 0) then - call dgemm('N', 'N', elec_beta_num, n_svd_coefs, det_beta_num, 1.d0 & - , det_beta_grad_lapl(mm,:,:), size(det_beta_grad_lapl,2) & - , psi_svd_beta(:,:,1), size(psi_svd_beta,1) & - , 0.d0, det_beta_grad_lapl_SVD(mm,:,:), size(det_beta_grad_lapl_SVD,2) ) - endif - ! !!! - enddo + call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs, det_alpha_num, 1.d0 & + , det_alpha_grad_lapl, 4*size(det_alpha_grad_lapl,2) & + , psi_svd_alpha, size(psi_svd_alpha,1) & + , 0.d0, det_alpha_grad_lapl_SVD, 4*size(det_alpha_grad_lapl_SVD,2) ) + ! !!! + if (elec_beta_num /= 0) then + call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs, det_beta_num, 1.d0 & + , det_beta_grad_lapl, 4*size(det_beta_grad_lapl,2) & + , psi_svd_beta, size(psi_svd_beta,1) & + , 0.d0, det_beta_grad_lapl_SVD, 4*size(det_beta_grad_lapl_SVD,2) ) + endif ! !!! ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ @@ -1951,7 +1938,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, psidet_grad_lapl_SVD, (4,elec_num) ] &BEGIN_PROVIDER [ double precision, pseudo_non_local_SVD, (elec_num) ] - + implicit none BEGIN_DOC ! !!! @@ -1962,8 +1949,8 @@ END_PROVIDER ! Non-local component of the pseudopotentials after SVD ! !!! ! for each electron: - ! for mm = 1, 2, 3 : grad_x Psi, grad_y Psi, grad_z Psi - ! for mm = 4: first term (only) of Laplacian Psi + ! for mm = 1, 2, 3 : grad_x Psi, grad_y Psi, grad_z Psi + ! for mm = 4: first term (only) of Laplacian Psi ! !!! END_DOC @@ -1985,7 +1972,7 @@ END_PROVIDER tmp = 0.d0 do l = 1, n_svd_coefs tmp = tmp + det_alpha_grad_lapl_SVD(mm,ee,l) * psi_svd_coefs(l,1) & - * det_beta_value_SVD(l) + * det_beta_value_SVD(l) enddo psidet_grad_lapl_SVD(mm,ee) = tmp enddo @@ -2012,7 +1999,7 @@ END_PROVIDER do ee = 1, elec_alpha_num tmp = 0.d0 do l = 1, n_svd_coefs - tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l) + tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l) enddo pseudo_non_local_SVD(ee) = tmp * psidet_inv_SVD enddo @@ -2086,13 +2073,13 @@ END_PROVIDER ! ! pair of | u_i v_j > to select from SVD ! ! !!! ! END_DOC -! integer :: l, i, j +! integer :: l, i, j ! do l = 1, n_svd_toselect ! i = psi_svd_alpha_numtoselect(l,1) ! j = psi_svd_beta_numtoselect (l,1) ! psi_svd_alpha_toselect(:,l,1) = psi_svd_alpha_unique(:,i,1) ! psi_svd_beta_toselect (:,l,1) = psi_svd_beta_unique (:,j,1) -! enddo +! enddo ! END_PROVIDER !____________________________________________________________________________________________________________________ @@ -2164,7 +2151,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) BEGIN_DOC ! Gradient of the current alpha determinant END_DOC - + integer :: i, j, k !DIR$ VECTOR ALIGNED do i=1,elec_alpha_num @@ -2173,7 +2160,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) det_alpha_grad_lapl_curr(3,i) = 0.d0 det_alpha_grad_lapl_curr(4,i) = 0.d0 enddo - + integer :: imo, imo2 ! ------- @@ -2184,7 +2171,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) ! do i=1,elec_alpha_num ! do k=1,4 ! det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + & -! mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j) +! mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j) ! enddo ! enddo ! enddo @@ -2234,13 +2221,13 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) do i=1,elec_alpha_num !DIR$ VECTOR ALIGNED do k=1,4 - det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) enddo enddo endif - - + + END_PROVIDER @@ -2249,7 +2236,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 BEGIN_DOC ! Gradient and Laplacian of the current beta determinant END_DOC - + integer :: i, j, k, l !DIR$ VECTOR ALIGNED @@ -2259,7 +2246,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 det_beta_grad_lapl_curr(3,i) = 0.d0 det_beta_grad_lapl_curr(4,i) = 0.d0 enddo - + integer :: imo, imo2 ! ------- @@ -2275,7 +2262,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 ! enddo ! enddo ! -! -------- +! -------- if (iand(elec_beta_num,1) == 0) then @@ -2289,10 +2276,10 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & - mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) enddo enddo enddo @@ -2309,10 +2296,10 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & - mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) enddo enddo i = elec_num @@ -2321,7 +2308,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) enddo enddo @@ -2336,7 +2323,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 enddo enddo - endif + endif END_PROVIDER From 248545e558effff89eb124b999baab562bb8db8d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 14:01:33 +0200 Subject: [PATCH 4/4] DIR -> OMP --- configure.sh | 4 +- src/det.irp.f | 340 +++++++++++++++++++++++++------------------------- src/mo.irp.f | 170 ++++++++++++------------- 3 files changed, 257 insertions(+), 257 deletions(-) diff --git a/configure.sh b/configure.sh index 535eea0..c48dfe7 100755 --- a/configure.sh +++ b/configure.sh @@ -39,11 +39,11 @@ export PATH="\${QMCCHEM_PATH}/bin:\${PATH}" export LD_LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LD_LIBRARY_PATH}" export LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LIBRARY_PATH}" export QMCCHEM_MPIRUN="mpirun" -export QMCCHEM_MPIRUN_FLAGS="--bind-to-core" +export QMCCHEM_MPIRUN_FLAGS="" export C_INCLUDE_PATH="\${QMCCHEM_PATH}/include:\${C_INCLUDE_PATH}" #export QMCCHEM_NIC=ib0 source \${QMCCHEM_PATH}/irpf90/bin/irpman #source \${QMCCHEM_PATH}/EZFIO/Bash/ezfio.sh -eval $(opam env) +eval \$(opam env) EOF diff --git a/src/det.irp.f b/src/det.irp.f index a7c0fca..91e2034 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -1,32 +1,32 @@ BEGIN_PROVIDER [ integer, det_i ] - + BEGIN_DOC ! Current running alpha determinant END_DOC det_i=det_alpha_order(1) - + END_PROVIDER BEGIN_PROVIDER [ integer, det_j ] - + BEGIN_DOC ! Current running beta determinant END_DOC det_j=det_beta_order(1) - + END_PROVIDER subroutine det_update(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(LDS) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + if (d == 0.d0) then return endif @@ -193,15 +193,15 @@ end subroutine det_update2(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(2) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,2) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,2) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + S(1,l) = m(1) S(2,l) = m(2) S_inv(1,1) = S(1,1) @@ -209,37 +209,37 @@ subroutine det_update2(n,LDS,m,l,S,S_inv,d) S_inv(2,1) = S(1,2) S_inv(2,2) = S(2,2) call invert2(S_inv,LDS,n,d) - + end subroutine det_update1(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(1) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,1) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,1) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + S(1,l) = m(1) S_inv(1,1) = S(1,1) call invert1(S_inv,LDS,n,d) - + end subroutine det_update3(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(3) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,3) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,3) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + integer :: i do i=1,3 S(i,l) = m(i) @@ -249,22 +249,22 @@ subroutine det_update3(n,LDS,m,l,S,S_inv,d) S_inv(2,i) = S(i,2) S_inv(3,i) = S(i,3) enddo - + call invert3(S_inv,LDS,n,d) - + end subroutine det_update4(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(4) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,4) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,4) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u(4), z(4), w(4), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u integer :: i,j @@ -276,7 +276,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) z(2) = S_inv(1,2)*u(1) + S_inv(2,2)*u(2) + S_inv(3,2)*u(3) + S_inv(4,2)*u(4) z(3) = S_inv(1,3)*u(1) + S_inv(2,3)*u(2) + S_inv(3,3)*u(3) + S_inv(4,3)*u(4) z(4) = S_inv(1,4)*u(1) + S_inv(2,4)*u(2) + S_inv(3,4)*u(3) + S_inv(4,4)*u(4) - + d_inv = 1.d0/d d = d+z(l) lambda = d_inv*d @@ -284,7 +284,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do i=1,4 w(i) = S_inv(i,l)*d_inv @@ -297,22 +297,22 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d) S_inv(j,i) = S_inv(j,i)*lambda -z(i)*w(j) enddo enddo - + end BEGIN_TEMPLATE ! Version for mod(n,4) = 0 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -321,20 +321,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (LDS >= $n) integer :: i,j double precision :: zj, zj1, zj2, zj3 - + !DIR$ NOPREFETCH - !DIR$ SIMD NOVECREMAINDER + !OMP$ SIMD NOVECREMAINDER do i=1,$n u(i) = m(i) - S(i,l) enddo - + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-1,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo d_inv = 1.d0/d @@ -344,7 +344,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n,4 zj = 0.d0 @@ -353,26 +353,26 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) zj2 = zj2 + S_inv(i,j+2)*u(i) zj3 = zj3 + S_inv(i,j+3)*u(i) enddo - z(j ) = zj - z(j+1) = zj1 - z(j+2) = zj2 + z(j ) = zj + z(j+1) = zj1 + z(j+2) = zj2 z(j+3) = zj3 enddo !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n,4 zj = z(i ) zj1 = z(i+1) @@ -380,7 +380,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 @@ -435,15 +435,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 1 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -452,20 +452,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (LDS >= $n) integer :: i,j double precision :: zj, zj1, zj2, zj3 - + do i=1,$n u(i) = m(i) - S(i,l) enddo zj = 0.d0 !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) + !OMP$ SIMD REDUCTION(+:zj) do i=1,$n-1,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo zj = zj + S_inv($n,l)*u($n) - + d_inv = 1.d0/d d = d+zj lambda = d*d_inv @@ -473,7 +473,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-1,4 zj = 0.d0 @@ -482,7 +482,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n-1 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -498,19 +498,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-1 zj = zj + S_inv(i,$n)*u(i) enddo z($n) = zj + S_inv($n,$n)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-1,4 zj = z(i ) zj1 = z(i+1) @@ -518,14 +518,14 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n-1 S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 enddo - S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj + S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*zj1 S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*zj2 S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*zj3 @@ -534,7 +534,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = z($n) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER do i=1,$n S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*zj enddo @@ -588,15 +588,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 2 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -604,24 +604,24 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) !DIR$ ASSUME (LDS >= $n) integer :: i,j - + double precision :: zj, zj1, zj2, zj3 !DIR$ NOPREFETCH - !DIR$ SIMD NOVECREMAINDER + !OMP$ SIMD NOVECREMAINDER do i=1,$n u(i) = m(i) - S(i,l) enddo - + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-2,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo i=$n-1 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) d_inv = 1.d0/d d = d+zj @@ -630,7 +630,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-2,4 zj = 0.d0 @@ -638,7 +638,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 zj3 = 0.d0 !DIR$ VECTOR ALIGNED - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER do i=1,$n-2 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -660,7 +660,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj1 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER do i=1,$n-2 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -671,19 +671,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-2,4 zj = z(i) zj1 = z(i+1) zj2 = z(i+2) zj3 = z(i+3) !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER do j=1,$n-2 S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -704,7 +704,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj = z(i) zj1= z(i+1) !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1) do j=1,$n-2 S_inv(j,i ) = S_inv(j,i )*lambda -zj*w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -761,15 +761,15 @@ BEGIN_TEMPLATE ! Version for mod(n,4) = 3 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m($n) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,$n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: u($n), z($n), w($n), lambda, d_inv !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN @@ -780,18 +780,18 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) double precision :: zj, zj1, zj2, zj3 - !DIR$ SIMD + !OMP$ SIMD do i=1,$n u(i) = m(i) - S(i,l) enddo - - zj = 0.d0 + + zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER do i=1,$n-3,4 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) enddo i=$n-2 zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + S_inv(i+2,l)*u(i+2) @@ -804,7 +804,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) d = 0.d0 return endif - + !DIR$ VECTOR ALIGNED do j=1,$n-3,4 zj = 0.d0 @@ -812,7 +812,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 zj3 = 0.d0 !DIR$ VECTOR ALIGNED - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) do i=1,$n-3 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -839,7 +839,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2) do i=1,$n-3 zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -856,12 +856,12 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,$n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + do i=1,$n-3,4 zj = z(i) zj1 = z(i+1) @@ -869,9 +869,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) do j=1,$n-3 - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 @@ -896,7 +896,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d) zj2 = z(i+2) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2) do j=1,$n S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 @@ -950,15 +950,15 @@ END_TEMPLATE subroutine det_update_general(n,LDS,m,l,S,S_inv,d) implicit none - + integer, intent(in) :: n,LDS ! Dimension of the vector real, intent(in) :: m(LDS) ! New vector integer, intent(in) :: l ! New position in S - + real,intent(inout) :: S(LDS,n) ! Slater matrix double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix double precision,intent(inout) :: d ! Det(S) - + double precision :: lambda, d_inv double precision :: u(3840), z(3840), w(3840) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u @@ -976,19 +976,19 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) do i=1,n u(i) = m(i) - S(i,l) enddo - + zl = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zl) + !OMP$ SIMD REDUCTION(+:zl) do i=1,n zl = zl + S_inv(i,l)*u(i) enddo - d_inv = 1.d0/d + d_inv = 1.d0/d d = d+zl lambda = d*d_inv - + if ( dabs(lambda) < 1.d-3 ) then d = 0.d0 return @@ -1006,7 +1006,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj3 = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) do i=1,n zj = zj + S_inv(i,j )*u(i) zj1 = zj1 + S_inv(i,j+1)*u(i) @@ -1023,7 +1023,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj = 0.d0 !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD REDUCTION(+:zj) + !OMP$ SIMD REDUCTION(+:zj) do i=1,n zj = zj + S_inv(i,j)*u(i) enddo @@ -1031,14 +1031,14 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) enddo !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) enddo - + !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(d_inv) + !OMP$ SIMD FIRSTPRIVATE(d_inv) do i=1,n w(i) = S_inv(i,l)*d_inv S(i,l) = m(i) @@ -1051,7 +1051,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj3 = z(i+3) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) do j=1,n S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) @@ -1064,7 +1064,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) zj = z(i) !DIR$ VECTOR ALIGNED !DIR$ NOPREFETCH - !DIR$ SIMD FIRSTPRIVATE(lambda,zj) + !OMP$ SIMD FIRSTPRIVATE(lambda,zj) do j=1,n S_inv(j,i) = S_inv(j,i)*lambda -zj*w(j) enddo @@ -1146,9 +1146,9 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, det_alpha_value_curr ] &BEGIN_PROVIDER [ real, slater_matrix_alpha, (elec_alpha_num_8,elec_alpha_num) ] &BEGIN_PROVIDER [ double precision, slater_matrix_alpha_inv_det, (elec_alpha_num_8,elec_alpha_num) ] - + implicit none - + BEGIN_DOC ! det_alpha_value_curr : Value of the current alpha determinant ! @@ -1159,7 +1159,7 @@ END_PROVIDER ! ! slater_matrix_alpha_inv_det: Inverse of the alpha Slater matrix x determinant END_DOC - + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num), n_to_do_old, n_to_do @@ -1172,7 +1172,7 @@ END_PROVIDER !DIR$ VECTOR ALIGNED slater_matrix_alpha_inv_det = 0.d0 endif - + PROVIDE mo_value if (det_i /= det_alpha_order(1) ) then ! if (det_i == -1 ) then @@ -1187,7 +1187,7 @@ END_PROVIDER enddo ddet = 0.d0 - + if (n_to_do < shiftl(elec_alpha_num,1)) then do while ( n_to_do > 0 ) @@ -1218,13 +1218,13 @@ END_PROVIDER enddo endif - + else ddet = 0.d0 endif - + ! Avoid NaN if (ddet /= 0.d0) then continue @@ -1246,11 +1246,11 @@ END_PROVIDER enddo enddo call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet) - + endif ASSERT (ddet /= 0.d0) - - det_alpha_value_curr = ddet + + det_alpha_value_curr = ddet END_PROVIDER BEGIN_PROVIDER [ double precision, det_beta_value_curr ] @@ -1266,24 +1266,24 @@ END_PROVIDER ! ! slater_matrix_beta_inv_det : Inverse of the beta Slater matrix x determinant END_DOC - + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num-mo_closed_num), n_to_do_old, n_to_do - + integer, save :: ifirst if (elec_beta_num == 0) then det_beta_value_curr = 0.d0 return endif - + if (ifirst == 0) then ifirst = 1 slater_matrix_beta = 0. slater_matrix_beta_inv_det = 0.d0 endif PROVIDE mo_value - + if (det_j /= det_beta_order(1)) then ! if (det_j == -1) then @@ -1295,7 +1295,7 @@ END_PROVIDER to_do(n_to_do) = k endif enddo - + ddet = 0.d0 if (n_to_do < shiftl(elec_beta_num,1)) then @@ -1327,13 +1327,13 @@ END_PROVIDER enddo endif - + else - + ddet = 0.d0 endif - + ! Avoid NaN if (ddet /= 0.d0) then continue @@ -1357,9 +1357,9 @@ END_PROVIDER call invert(slater_matrix_beta_inv_det,elec_beta_num_8,elec_beta_num,ddet) endif ASSERT (ddet /= 0.d0) - - det_beta_value_curr = ddet - + + det_beta_value_curr = ddet + END_PROVIDER BEGIN_PROVIDER [ integer, det_alpha_num_pseudo ] @@ -1381,15 +1381,15 @@ END_PROVIDER BEGIN_PROVIDER [ double precision , det_alpha_value, (det_alpha_num_8) ] &BEGIN_PROVIDER [ double precision , det_alpha_grad_lapl, (4,elec_alpha_num,det_alpha_num) ] &BEGIN_PROVIDER [ double precision , det_alpha_pseudo, (elec_alpha_num_8,det_alpha_num_pseudo) ] - + implicit none - + BEGIN_DOC ! Values of the alpha determinants ! Gradients of the alpha determinants ! Laplacians of the alpha determinants END_DOC - + integer :: j,i,k integer, save :: ifirst = 0 if (ifirst == 0) then @@ -1398,73 +1398,73 @@ END_PROVIDER det_alpha_grad_lapl = 0.d0 det_alpha_pseudo = 0.d0 endif - - + + do j=1,det_alpha_num - + det_i = det_alpha_order(j) if (j > 1) then TOUCH det_i endif - + det_alpha_value(det_i) = det_alpha_value_curr det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:) if (do_pseudo) then - det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) + det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) endif - + enddo - + det_i = det_alpha_order(1) SOFT_TOUCH det_i - + END_PROVIDER BEGIN_PROVIDER [ double precision, det_beta_value, (det_beta_num_8) ] &BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_alpha_num+1:elec_num,det_beta_num) ] &BEGIN_PROVIDER [ double precision, det_beta_pseudo, (elec_alpha_num+1:elec_num,det_beta_num_pseudo) ] - - + + implicit none - + BEGIN_DOC ! Values of the beta determinants ! Gradients of the beta determinants ! Laplacians of the beta determinants END_DOC - + integer :: j,i,k integer, save :: ifirst = 0 if (elec_beta_num == 0) then det_beta_value = 1.d0 return endif - + if (ifirst == 0) then ifirst = 1 det_beta_value = 0.d0 det_beta_grad_lapl = 0.d0 det_beta_pseudo = 0.d0 endif - + do j=1,det_beta_num - + det_j = det_beta_order(j) if (j > 1) then TOUCH det_j endif - + det_beta_value(det_j) = det_beta_value_curr det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:) if (do_pseudo) then - det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) + det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) endif enddo - + det_j = det_beta_order(1) SOFT_TOUCH det_j - + END_PROVIDER @@ -1472,34 +1472,34 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, psidet_inv ] &BEGIN_PROVIDER [ double precision, psidet_grad_lapl, (4,elec_num_8) ] &BEGIN_PROVIDER [ double precision, pseudo_non_local, (elec_num) ] - + implicit none BEGIN_DOC ! Value of the determinantal part of the wave function - + ! Gradient of the determinantal part of the wave function - + ! Laplacian of determinantal part of the wave function ! Non-local component of the pseudopotentials ! Regularized 1/psi = 1/(psi + 1/psi) END_DOC - + integer, save :: ifirst=0 if (ifirst == 0) then ifirst = 1 psidet_grad_lapl = 0.d0 endif - + double precision :: CDb(det_alpha_num_8) double precision :: CDb_i double precision :: DaC(det_beta_num_8) !DIR$ ATTRIBUTES ALIGN : 32 :: DaC,CDb ! C x D_beta - ! D_alpha^t x C - ! D_alpha^t x (C x D_beta) + ! D_alpha^t x C + ! D_alpha^t x (C x D_beta) integer :: i,j,k, l integer :: i1,i2,i3,i4,det_num4 @@ -1576,7 +1576,7 @@ END_PROVIDER psidet_value = psidet_value + det_beta_value(j) * DaC(j) enddo - + if (psidet_value == 0.d0) then call abrt(irp_here,'Determinantal component of the wave function is zero.') endif @@ -1595,7 +1595,7 @@ END_PROVIDER size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2), & DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1) endif - + if (do_pseudo) then call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv, & det_alpha_pseudo, size(det_alpha_pseudo,1), & @@ -1669,7 +1669,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) BEGIN_DOC ! Gradient of the current alpha determinant END_DOC - + integer :: i, j, k !DIR$ VECTOR ALIGNED do i=1,elec_alpha_num @@ -1678,7 +1678,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) det_alpha_grad_lapl_curr(3,i) = 0.d0 det_alpha_grad_lapl_curr(4,i) = 0.d0 enddo - + integer :: imo, imo2 ! ------- @@ -1688,7 +1688,7 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) ! imo = mo_list_alpha_curr(j) ! do i=1,elec_alpha_num ! do k=1,4 -! det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j) +! det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j) ! enddo ! enddo ! enddo @@ -1738,13 +1738,13 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) do i=1,elec_alpha_num !DIR$ VECTOR ALIGNED do k=1,4 - det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) enddo enddo endif - - + + END_PROVIDER @@ -1753,7 +1753,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 BEGIN_DOC ! Gradient and Laplacian of the current beta determinant END_DOC - + integer :: i, j, k, l !DIR$ VECTOR ALIGNED @@ -1763,7 +1763,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 det_beta_grad_lapl_curr(3,i) = 0.d0 det_beta_grad_lapl_curr(4,i) = 0.d0 enddo - + integer :: imo, imo2 ! ------- @@ -1779,7 +1779,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 ! enddo ! enddo ! -! -------- +! -------- if (iand(elec_beta_num,1) == 0) then @@ -1793,10 +1793,10 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & - mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) enddo enddo enddo @@ -1813,10 +1813,10 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & - mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) enddo enddo i = elec_num @@ -1825,7 +1825,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 do k=1,4 det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & - mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) enddo enddo @@ -1840,7 +1840,7 @@ BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1 enddo enddo - endif + endif END_PROVIDER diff --git a/src/mo.irp.f b/src/mo.irp.f index 38228b2..124890d 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -5,12 +5,12 @@ ! Number of Molecular orbitals END_DOC integer, external :: mod_align - + mo_num = maxval(present_mos) call iinfo(irp_here,'mo_num',mo_num) - + mo_num_8 = mod_align(mo_num) - + END_PROVIDER @@ -22,7 +22,7 @@ BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ] integer :: i, j real,allocatable :: buffer(:,:) allocate (buffer(ao_num,mo_tot_num)) - + buffer = 0. call get_mo_basis_mo_coef(buffer) do i=1,mo_tot_num @@ -35,7 +35,7 @@ BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ] call set_order(mo_coef_input(1,i),ao_nucl_sort_idx,ao_num) enddo deallocate(buffer) - + END_PROVIDER @@ -56,7 +56,7 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ] ! Molecular orbital coefficients END_DOC integer :: i, j - + do j=1,mo_num do i=1,ao_num_8 mo_coef(i,j) = mo_coef_input(i,j) @@ -68,10 +68,10 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ] mo_coef(i,j) = 0. enddo enddo - + ! Input MOs are not needed any more FREE mo_coef_input - + real :: f f = 1./mo_scale do j=1,mo_num @@ -80,7 +80,7 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ] mo_coef(i,j) *= f enddo enddo - + END_PROVIDER @@ -90,7 +90,7 @@ BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num_8,ao_num_8) ] ! Transpose of the Molecular orbital coefficients END_DOC call transpose(mo_coef,ao_num_8,mo_coef_transp,mo_num_8,ao_num_8,mo_num_8) - + END_PROVIDER @@ -102,7 +102,7 @@ END_PROVIDER ! orbital coefficients END_DOC integer :: i, j - + integer :: idx mo_coef_transp_sparsity = 0. do j=1,ao_num @@ -117,7 +117,7 @@ END_PROVIDER mo_coef_transp_sparsity += float(idx) enddo mo_coef_transp_sparsity *= 1./(mo_num*ao_num) - + END_PROVIDER @@ -133,7 +133,7 @@ BEGIN_PROVIDER [ real, mo_coef_transp_present, (num_present_mos_8,ao_num_8) ] mo_coef_transp_present(j,i) = mo_coef_transp(present_mos(j),i) enddo enddo - + END_PROVIDER @@ -143,24 +143,24 @@ END_PROVIDER &BEGIN_PROVIDER [ real, mo_grad_transp_z, (mo_num_8,elec_num) ] &BEGIN_PROVIDER [ real, mo_lapl_transp, (mo_num_8,elec_num) ] implicit none - + BEGIN_DOC ! Values, gradients, laplacians of the molecular orbitals ! ! Arrays are padded for efficiency END_DOC - + integer :: i, j, k, l, m - + PROVIDE primitives_reduced - + if (do_nucl_fitcusp) then PROVIDE nucl_fitcusp_param PROVIDE nucl_elec_dist_vec PROVIDE nucl_elec_dist_inv endif - - + + do i=1,elec_num if (i>1) then ao_elec = i @@ -180,7 +180,7 @@ END_PROVIDER mo_grad_transp_z(1,i), & mo_lapl_transp(1,i), & ao_num) - + else call sparse_full_mv(mo_coef_transp_present,num_present_mos_8, & ao_value_block(1),ao_num_8, & @@ -195,7 +195,7 @@ END_PROVIDER mo_grad_transp_z(1,i), & mo_lapl_transp(1,i), & ao_num) - + do j=num_present_mos,1,-1 mo_value_transp (present_mos(j),i) = mo_value_transp (j,i) mo_grad_transp_x(present_mos(j),i) = mo_grad_transp_x(j,i) @@ -207,10 +207,10 @@ END_PROVIDER endif enddo endif - + if (do_nucl_fitcusp) then real :: r, r2, r_inv, d, expzr, Z, Z2, a, b, c, phi, rx, ry, rz - + do k=1,nucl_num r = nucl_elec_dist(k,i) if (r > nucl_fitcusp_radius(k)) then @@ -235,14 +235,14 @@ END_PROVIDER enddo exit enddo ! k - + endif - + enddo ! i ao_elec = 1 SOFT_TOUCH ao_elec - - + + if (do_prepare) then real :: lambda, t ! Scale off-diagonal elements @@ -268,7 +268,7 @@ END_PROVIDER enddo enddo endif - + do i=1,mo_num do j=1,elec_num mo_value_transp(i,j) *= mo_cusp_rescale(i) @@ -285,10 +285,10 @@ BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ] BEGIN_DOC ! Values of the molecular orbitals END_DOC - + integer :: i,j integer, save :: ifirst = 0 - + if (ifirst == 0) then ifirst = 1 PROVIDE primitives_reduced @@ -296,7 +296,7 @@ BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ] mo_value = 0. endif call transpose(mo_value_transp(1,1),mo_num_8,mo_value,elec_num_8,mo_num,elec_num) - + END_PROVIDER @@ -304,14 +304,14 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, mo_grad_y, (elec_num_8,mo_num) ] &BEGIN_PROVIDER [ double precision, mo_grad_z, (elec_num_8,mo_num) ] implicit none - + BEGIN_DOC ! Gradients of the molecular orbitals END_DOC - + integer :: i,j integer, save :: ifirst = 0 - + if (ifirst == 0) then !DIR$ VECTOR ALIGNED mo_grad_x = 0.d0 @@ -327,7 +327,7 @@ END_PROVIDER call transpose_to_dp(mo_grad_transp_z(1,1),mo_num_8,mo_grad_z(1,1),elec_num_8,mo_num,elec_num) call transpose_to_dp(mo_grad_transp_x(1,1),mo_num_8,mo_grad_x(1,1),elec_num_8,mo_num,elec_num) - + END_PROVIDER BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ] @@ -335,7 +335,7 @@ BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ] BEGIN_DOC ! Laplacians of the molecular orbitals END_DOC - + integer :: i,j integer, save :: ifirst = 0 if (ifirst == 0) then @@ -345,7 +345,7 @@ BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ] mo_lapl = 0.d0 endif call transpose_to_dp(mo_lapl_transp(1,1),mo_num_8,mo_lapl,elec_num_8,mo_num,elec_num) - + END_PROVIDER @@ -358,21 +358,21 @@ END_PROVIDER integer :: i,j do j=1,mo_num do i=1,elec_alpha_num - mo_grad_lapl_alpha(1,i,j) = mo_grad_transp_x(j,i) - mo_grad_lapl_alpha(2,i,j) = mo_grad_transp_y(j,i) - mo_grad_lapl_alpha(3,i,j) = mo_grad_transp_z(j,i) - mo_grad_lapl_alpha(4,i,j) = mo_lapl_transp (j,i) + mo_grad_lapl_alpha(1,i,j) = mo_grad_transp_x(j,i) + mo_grad_lapl_alpha(2,i,j) = mo_grad_transp_y(j,i) + mo_grad_lapl_alpha(3,i,j) = mo_grad_transp_z(j,i) + mo_grad_lapl_alpha(4,i,j) = mo_lapl_transp (j,i) enddo enddo do j=1,mo_num do i=elec_alpha_num+1,elec_num - mo_grad_lapl_beta(1,i,j) = mo_grad_transp_x(j,i) - mo_grad_lapl_beta(2,i,j) = mo_grad_transp_y(j,i) - mo_grad_lapl_beta(3,i,j) = mo_grad_transp_z(j,i) - mo_grad_lapl_beta(4,i,j) = mo_lapl_transp (j,i) + mo_grad_lapl_beta(1,i,j) = mo_grad_transp_x(j,i) + mo_grad_lapl_beta(2,i,j) = mo_grad_transp_y(j,i) + mo_grad_lapl_beta(3,i,j) = mo_grad_transp_z(j,i) + mo_grad_lapl_beta(4,i,j) = mo_lapl_transp (j,i) enddo enddo - + END_PROVIDER BEGIN_PROVIDER [ double precision, mo_grad_lapl_transp, (4,mo_num,elec_num) ] @@ -389,7 +389,7 @@ BEGIN_PROVIDER [ double precision, mo_grad_lapl_transp, (4,mo_num,elec_num) ] mo_grad_lapl_transp(4,j,i) = mo_lapl_transp (j,i) enddo enddo - + END_PROVIDER @@ -404,18 +404,18 @@ END_PROVIDER BEGIN_PROVIDER [ integer, mo_tot_num ] - + BEGIN_DOC ! Total number of MOs in the EZFIO file END_DOC - + mo_tot_num = -1 call get_mo_basis_mo_num(mo_tot_num) if (mo_tot_num <= 0) then call abrt(irp_here,'Total number of MOs can''t be <0') endif call iinfo(irp_here,'mo_tot_num',mo_tot_num) - + END_PROVIDER @@ -432,16 +432,16 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ] END_DOC integer :: i, j, k, l real :: ao_value_at_nucl_no_S(ao_num) - + PROVIDE mo_fitcusp_normalization_before do k=1,nucl_num point(1) = nucl_coord(k,1) point(2) = nucl_coord(k,2) point(3) = nucl_coord(k,3) TOUCH point - + PROVIDE ao_value_p - + do i=1,ao_num if (ao_nucl(i) /= k) then ao_value_at_nucl_no_S(i) = ao_value_p(i) @@ -449,7 +449,7 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ] ao_value_at_nucl_no_S(i) = 0. endif enddo - + integer :: jj do jj=1,num_present_mos j = present_mos(jj) @@ -459,11 +459,11 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ] mo_value_at_nucl(j,k) = mo_value_at_nucl(j,k) + mo_coef(i,j)*ao_value_at_nucl_no_S(i) enddo enddo - + enddo FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p SOFT_TOUCH point - + END_PROVIDER @@ -474,15 +474,15 @@ END_PROVIDER BEGIN_DOC ! Values of the atomic orbitals without S components on atoms END_DOC - + integer :: i, j, k - + do k=1,nucl_num point(1) = nucl_coord(k,1) point(2) = nucl_coord(k,2) point(3) = nucl_coord(k,3)+ nucl_fitcusp_radius(k) TOUCH point - + do j=1,ao_num ao_value_at_fitcusp_radius(j,k) = ao_value_p(j) ao_grad_at_fitcusp_radius(j,k) = ao_grad_p(j,3) @@ -496,7 +496,7 @@ END_PROVIDER enddo FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p SOFT_TOUCH point - + END_PROVIDER BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_before, (mo_tot_num) ] @@ -560,7 +560,7 @@ BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) t = 0.d0 do j=1,ao_num if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then - t = t + mo_coef(j,i) * ao_value_p(j) + t = t + mo_coef(j,i) * ao_value_p(j) endif enddo t = t + nucl_fitcusp_param(1,i,k) + & @@ -601,7 +601,7 @@ END_PROVIDER ! Values of the molecular orbitals without S components on atoms END_DOC integer :: i, j, k, l - + do k=1,nucl_num do j=1,mo_num mo_value_at_fitcusp_radius(j,k) = 0.d0 @@ -615,7 +615,7 @@ END_PROVIDER enddo enddo enddo - + END_PROVIDER @@ -626,7 +626,7 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ] END_DOC integer :: i,k, niter character*(80) :: message - + nucl_fitcusp_param = 0.d0 do k=1,nucl_num double precision :: r, Z @@ -642,14 +642,14 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ] grad_phi = mo_grad_at_fitcusp_radius(i,k) phi = mo_value_at_fitcusp_radius(i,k) eta = mo_value_at_nucl(i,k) - + nucl_fitcusp_param(1,i,k) = -(R*(2.d0*eta*Z-6.d0*grad_phi)+lap_phi*R*R+6.d0*phi)/(2.d0*R*Z-6.d0) nucl_fitcusp_param(2,i,k) = (lap_phi*R*R*Z-6.d0*grad_phi*R*Z+6.d0*phi*Z+6.d0*eta*Z)/(2.d0*R*Z-6.d0) nucl_fitcusp_param(3,i,k) = -(R*(-5.d0*grad_phi*Z-1.5d0*lap_phi)+lap_phi*R*R*Z+3.d0*phi*Z+& 3.d0*eta*Z+6.d0*grad_phi)/(R*R*Z-3.d0*R) nucl_fitcusp_param(4,i,k) = (R*(-2.d0*grad_phi*Z-lap_phi)+0.5d0*lap_phi*R*R*Z+phi*Z+& eta*Z+3.d0*grad_phi)/(R*R*R*Z-3.d0*R*R) - + enddo enddo END_PROVIDER @@ -689,7 +689,7 @@ subroutine sparse_full_mv(A,LDA, & !DIR$ ASSUME_ALIGNED C3 : $IRP_ALIGN !DIR$ ASSUME_ALIGNED C4 : $IRP_ALIGN !DIR$ ASSUME_ALIGNED C5 : $IRP_ALIGN - + integer :: kao, kmax, kmax2, kmax3 integer :: i,j,k integer :: k_vec(8) @@ -698,9 +698,9 @@ subroutine sparse_full_mv(A,LDA, & real :: d21, d22, d23, d24, d25 real :: d31, d32, d33, d34, d35 real :: d41, d42, d43, d44, d45 - + ! LDC and LDA have to be factors of simd_sp - + ! IRP_IF NO_PREFETCH ! IRP_ELSE ! call MM_PREFETCH (A(1,indices(1)),3) @@ -709,7 +709,7 @@ subroutine sparse_full_mv(A,LDA, & ! call MM_PREFETCH (A(1,indices(4)),3) ! IRP_ENDIF - !DIR$ SIMD + !OMP$ SIMD do j=1,LDC C1(j) = 0. C2(j) = 0. @@ -727,35 +727,35 @@ subroutine sparse_full_mv(A,LDA, & k_vec(2) = indices(kao+1) k_vec(3) = indices(kao+2) k_vec(4) = indices(kao+3) - + d11 = B1(kao ) d21 = B1(kao+1) d31 = B1(kao+2) d41 = B1(kao+3) - + d12 = B2(kao ) d22 = B2(kao+1) d32 = B2(kao+2) d42 = B2(kao+3) - + d13 = B3(kao ) d23 = B3(kao+1) d33 = B3(kao+2) d43 = B3(kao+3) - + d14 = B4(kao ) d24 = B4(kao+1) d34 = B4(kao+2) d44 = B4(kao+3) - + d15 = B5(kao ) d25 = B5(kao+1) d35 = B5(kao+2) d45 = B5(kao+3) - + do k=0,LDA-1,$IRP_ALIGN/4 !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(d11,d21,d31,d41) + !OMP$ SIMD FIRSTPRIVATE(d11,d21,d31,d41) do j=1,$IRP_ALIGN/4 ! IRP_IF NO_PREFETCH ! IRP_ELSE @@ -767,18 +767,18 @@ subroutine sparse_full_mv(A,LDA, & C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 + A(j+k,k_vec(2))*d21& + A(j+k,k_vec(3))*d31 + A(j+k,k_vec(4))*d41 enddo - + !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(d12,d22,d32,d42,d13,d23,d33,d43) + !OMP$ SIMD FIRSTPRIVATE(d12,d22,d32,d42,d13,d23,d33,d43) do j=1,$IRP_ALIGN/4 C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12 + A(j+k,k_vec(2))*d22& + A(j+k,k_vec(3))*d32 + A(j+k,k_vec(4))*d42 C3(j+k) = C3(j+k) + A(j+k,k_vec(1))*d13 + A(j+k,k_vec(2))*d23& + A(j+k,k_vec(3))*d33 + A(j+k,k_vec(4))*d43 enddo - + !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(d14,d24,d34,d44,d15,d25,d35,d45) + !OMP$ SIMD FIRSTPRIVATE(d14,d24,d34,d44,d15,d25,d35,d45) do j=1,$IRP_ALIGN/4 C4(j+k) = C4(j+k) + A(j+k,k_vec(1))*d14 + A(j+k,k_vec(2))*d24& + A(j+k,k_vec(3))*d34 + A(j+k,k_vec(4))*d44 @@ -786,9 +786,9 @@ subroutine sparse_full_mv(A,LDA, & + A(j+k,k_vec(3))*d35 + A(j+k,k_vec(4))*d45 enddo enddo - + enddo - + do kao = kmax2+1, kmax3 k_vec(1) = indices(kao) d11 = B1(kao) @@ -799,7 +799,7 @@ subroutine sparse_full_mv(A,LDA, & !DIR$ VECTOR ALIGNED do k=0,LDA-1,$IRP_ALIGN/4 !DIR$ VECTOR ALIGNED - !DIR$ SIMD FIRSTPRIVATE(d11,d12,d13,d14,d15) + !OMP$ SIMD FIRSTPRIVATE(d11,d12,d13,d14,d15) do j=1,$IRP_ALIGN/4 C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12 @@ -809,7 +809,7 @@ subroutine sparse_full_mv(A,LDA, & enddo enddo enddo - + end