diff --git a/README.md b/README.md index b4b8912..04bec50 100644 --- a/README.md +++ b/README.md @@ -26,9 +26,13 @@ diffusion Monte Carlo algorithm. in Cloud environments (rance Grilles) coupled to supercomputers -*Warning*: QMC=Chem is under the GPLv2 license. Any modifications to or -software including (via compiler) GPL-licensed code must also be made available -under the GPL along with build & install instructions. +Warnings: +* QMC=Chem is under the GPLv2 license. Any modifications to or + software including (via compiler) GPL-licensed code must also be made available + under the GPL along with build & install instructions. +* Pseudopotentials are about to change in EZFIO database. Current calculations + will not be compatible with future versions + Requirements ------------ diff --git a/src/constants.F b/src/constants.F new file mode 100644 index 0000000..d4d35c8 --- /dev/null +++ b/src/constants.F @@ -0,0 +1,11 @@ + real, parameter :: pi=3.14159265359 + real, parameter :: two_pi=3.14159265359*2. + real, parameter :: sqpi=1.77245385091 + real, parameter :: sq3=1.7320508075688772 + real, parameter :: one_over_sqpi = 0.5641895835477563 + + double precision, parameter :: dpi=3.14159265359d0 + double precision, parameter :: dsqpi=1.77245385091d0 + double precision, parameter :: dsq3=1.7320508075688772d0 + double precision, parameter :: dtwo_pi=3.14159265359d0*2.d0 + double precision, parameter :: dfour_pi=3.14159265359d0*4.d0 diff --git a/src/det.irp.f b/src/det.irp.f new file mode 100644 index 0000000..4d5de55 --- /dev/null +++ b/src/det.irp.f @@ -0,0 +1,1551 @@ +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 + select case (n) + case default + call det_update_general(n,LDS,m,l,S,S_inv,d) + BEGIN_TEMPLATE + case ($n) + call det_update$n(n,LDS,m,l,S,S_inv,d) + SUBST [n] + 1;; + 2;; + 3;; + 4;; + 5;; + 6;; + 7;; + 8;; + 9;; + 10;; + 11;; + 12;; + 13;; + 14;; + 15;; + 16;; + 17;; + 18;; + 19;; + 20;; + 21;; + 22;; + 23;; + 24;; + 25;; + 26;; + 27;; + 28;; + 29;; + 30;; + 31;; + 32;; + 33;; + 34;; + 35;; + 36;; + 37;; + 38;; + 39;; + 40;; + 41;; + 42;; + 43;; + 44;; + 45;; + 46;; + 47;; + 48;; + 49;; + 50;; + END_TEMPLATE + end select +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) + S_inv(1,2) = S(2,1) + 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 + !DIR$ VECTOR ALIGNED + do i=1,3 + S(i,l) = m(i) + enddo + do i=1,3 + S_inv(1,i) = S(i,1) + 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 + u(1) = m(1) - S(1,l) + u(2) = m(2) - S(2,l) + u(3) = m(3) - S(3,l) + u(4) = m(4) - S(4,l) + z(1) = S_inv(1,1)*u(1) + S_inv(2,1)*u(2) + S_inv(3,1)*u(3) + S_inv(4,1)*u(4) + 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 + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + do i=1,4 + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,4 + !DIR$ VECTOR ALIGNED + do j=1,4 + 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 + !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN + !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) + !DIR$ ASSUME (LDS >= $n) + integer :: i,j + + !$OMP SIMD + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + !$OMP END SIMD + + z(l) = S_inv($n,l)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(l) = z(l) + S_inv(i,l)*u(i) + enddo + + d_inv = 1.d0/d + d = d+z(l) + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + do j=1,$n,4 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + z(j+2) = S_inv($n,j+2)*u($n) + z(j+3) = S_inv($n,j+3)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + z(j+3) = z(j+3) + S_inv(i,j+3)*u(i) + enddo + enddo + + !$OMP SIMD + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + !$OMP END SIMD + + do i=1,$n,4 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda -z(i )*w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -z(i+1)*w(j) + S_inv(j,i+2) = S_inv(j,i+2)*lambda -z(i+2)*w(j) + S_inv(j,i+3) = S_inv(j,i+3)*lambda -z(i+3)*w(j) + enddo + !$OMP END SIMD + enddo + +end + +SUBST [ n ] +8 ;; +12 ;; +16 ;; +20 ;; +24 ;; +28 ;; +32 ;; +36 ;; +40 ;; +44 ;; +48 ;; + +END_TEMPLATE + +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 + !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN + !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) + !DIR$ ASSUME (LDS >= $n) + integer :: i,j + + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + z(l) = S_inv($n,l)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(l) = z(l) + S_inv(i,l)*u(i) + enddo + + d_inv = 1.d0/d + d = d+z(l) + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + !DIR$ VECTOR ALIGNED + do j=1,$n-1,4 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + z(j+2) = S_inv($n,j+2)*u($n) + z(j+3) = S_inv($n,j+3)*u($n) + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + z(j+3) = z(j+3) + S_inv(i,j+3)*u(i) + enddo + enddo + + z($n) = S_inv($n,$n)*u($n) + do i=1,$n-1 + z($n) = z($n) + S_inv(i,$n)*u(i) + enddo + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n-1,4 + !$OMP SIMD + do j=1,$n-1 + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*z(i ) + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*z(i+1) + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*z(i+2) + S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*z(i+3) + enddo + !$OMP END SIMD + enddo + + do i=1,$n-1,4 + S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*z(i ) + S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*z(i+1) + S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*z(i+2) + S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*z(i+3) + enddo + + do i=1,$n + S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*z($n) + enddo + + +end + +SUBST [ n ] +5 ;; +9 ;; +13 ;; +17 ;; +21 ;; +25 ;; +29 ;; +33 ;; +37 ;; +41 ;; +45 ;; +49 ;; + +END_TEMPLATE + + +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 + !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN + !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) + !DIR$ ASSUME (LDS >= $n) + integer :: i,j + + !$OMP SIMD + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + !$OMP END SIMD + + z(l) = S_inv($n,l)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(l) = z(l) + S_inv(i,l)*u(i) + enddo + + d_inv = 1.d0/d + d = d+z(l) + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + do j=1,$n-2,4 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + z(j+2) = S_inv($n,j+2)*u($n) + z(j+3) = S_inv($n,j+3)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + z(j+3) = z(j+3) + S_inv(i,j+3)*u(i) + enddo + enddo + + j=$n-1 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + enddo + + !$OMP SIMD + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + !$OMP END SIMD + + do i=1,$n-2,4 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda -z(i )*w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -z(i+1)*w(j) + S_inv(j,i+2) = S_inv(j,i+2)*lambda -z(i+2)*w(j) + S_inv(j,i+3) = S_inv(j,i+3)*lambda -z(i+3)*w(j) + enddo + !$OMP END SIMD + enddo + + i=$n-1 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda -z(i )*w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -z(i+1)*w(j) + enddo + !$OMP END SIMD + +end + +SUBST [ n ] +6 ;; +10 ;; +14 ;; +18 ;; +22 ;; +26 ;; +30 ;; +34 ;; +38 ;; +42 ;; +46 ;; +50 ;; + +END_TEMPLATE + +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 + !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN + !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) + !DIR$ ASSUME (LDS >= $n) + integer :: i,j + + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + z(l) = S_inv($n,l)*u($n) + !DIR$ VECTOR ALIGNED + do i=1,$n-1 + z(l) = z(l) + S_inv(i,l)*u(i) + enddo + + + d_inv = 1.d0/d + d = d+z(l) + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + !DIR$ VECTOR ALIGNED + do j=1,$n-3,4 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + z(j+2) = S_inv($n,j+2)*u($n) + z(j+3) = S_inv($n,j+3)*u($n) + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + z(j+3) = z(j+3) + S_inv(i,j+3)*u(i) + enddo + enddo + + j=$n-2 + z(j ) = S_inv($n,j )*u($n) + z(j+1) = S_inv($n,j+1)*u($n) + z(j+2) = S_inv($n,j+2)*u($n) + do i=1,$n-1 + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + enddo + + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n-3,4 + !$OMP SIMD + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*z(i ) + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*z(i+1) + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*z(i+2) + S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*z(i+3) + enddo + !$OMP END SIMD + enddo + + i=$n-2 + !$OMP SIMD + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*z(i ) + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*z(i+1) + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*z(i+2) + enddo + !$OMP END SIMD + +end + +SUBST [ n ] +7 ;; +11 ;; +15 ;; +19 ;; +23 ;; +27 ;; +31 ;; +35 ;; +39 ;; +43 ;; +47 ;; + +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 + !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN + !DIR$ ASSUME (LDS >= n) + !DIR$ ASSUME (LDS <= 3840) + !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0) + integer :: i,j,n4 + + !$OMP SIMD + do i=1,n + u(i) = m(i) - S(i,l) + enddo + !$OMP END SIMD + + z(l) = 0.d0 + !DIR$ VECTOR ALIGNED + do i=1,n + z(l) = z(l) + S_inv(i,l)*u(i) + enddo + + d_inv = 1.d0/d + d = d+z(l) + lambda = d*d_inv + + if ( dabs(lambda) < 1.d-3 ) then + d = 0.d0 + return + endif + + n4 = iand(n,not(4)) + do j=1,n4,8 + z(j ) = 0.d0 + z(j+1) = 0.d0 + z(j+2) = 0.d0 + z(j+3) = 0.d0 + z(j+4) = 0.d0 + z(j+5) = 0.d0 + z(j+6) = 0.d0 + z(j+7) = 0.d0 + !DIR$ VECTOR ALIGNED + do i=1,n + z(j ) = z(j ) + S_inv(i,j )*u(i) + z(j+1) = z(j+1) + S_inv(i,j+1)*u(i) + z(j+2) = z(j+2) + S_inv(i,j+2)*u(i) + z(j+3) = z(j+3) + S_inv(i,j+3)*u(i) + z(j+4) = z(j+4) + S_inv(i,j+4)*u(i) + z(j+5) = z(j+5) + S_inv(i,j+5)*u(i) + z(j+6) = z(j+6) + S_inv(i,j+6)*u(i) + z(j+7) = z(j+7) + S_inv(i,j+7)*u(i) + enddo + enddo + + do j=n4+1,n + z(j) = 0.d0 + !DIR$ VECTOR ALIGNED + do i=1,n + z(j) = z(j) + S_inv(i,j)*u(i) + enddo + enddo + + !$OMP SIMD + do i=1,n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + !$OMP END SIMD + + do i=1,n + !DIR$ VECTOR ALIGNED + !$OMP SIMD aligned(S_inv,z) + do j=1,n + S_inv(j,i) = S_inv(j,i)*lambda -z(i)*w(j) + enddo + !$OMP END SIMD + enddo + +end + + + +subroutine bitstring_to_list( string, list, n_elements, Nint) + implicit none + BEGIN_DOC + ! Gives the inidices(+1) of the bits set to 1 in the bit string + END_DOC + integer, intent(in) :: Nint + integer*8, intent(in) :: string(Nint) + integer, intent(out) :: list(Nint*64) + integer, intent(out) :: n_elements + + integer :: i, ishift + integer*8 :: l + + n_elements = 0 + ishift = 2 + do i=1,Nint + l = string(i) + do while (l /= 0_8) + n_elements = n_elements+1 + list(n_elements) = ishift+popcnt(l-1_8) - popcnt(l) + l = iand(l,l-1_8) + enddo + ishift = ishift + 64 + enddo +end + + BEGIN_PROVIDER [ integer, mo_list_alpha_curr, (elec_alpha_num) ] +&BEGIN_PROVIDER [ integer, mo_list_alpha_prev, (elec_alpha_num) ] + implicit none + BEGIN_DOC + ! List of MOs in the current alpha determinant + END_DOC + integer :: l + if (det_i /= det_alpha_order(1)) then + mo_list_alpha_prev = mo_list_alpha_curr + else + mo_list_alpha_prev = 0 + endif + !DIR$ FORCEINLINE + call bitstring_to_list ( psi_det_alpha(1,det_i), mo_list_alpha_curr, l, N_int ) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, mo_list_beta_curr, (elec_beta_num) ] +&BEGIN_PROVIDER [ integer, mo_list_beta_prev, (elec_beta_num) ] + implicit none + BEGIN_DOC + ! List of MOs in the current beta determinant + END_DOC + integer :: l + if (elec_beta_num == 0) then + return + endif + if (det_j /= det_beta_order(1)) then + mo_list_beta_prev = mo_list_beta_curr + else + mo_list_beta_prev = 0 + endif + !DIR$ FORCEINLINE + call bitstring_to_list ( psi_det_beta(1,det_j), mo_list_beta_curr, l, N_int ) +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 + ! + ! det_alpha_value_curr : Slater matrix for the current alpha determinant. + ! 1st index runs over electrons and + ! 2nd index runs over MOs. + ! Built with 1st determinant + ! + ! 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 + integer, save :: ifirst + + if (ifirst == 0) then + ifirst = 1 + !DIR$ VECTOR ALIGNED + slater_matrix_alpha = 0. + !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 + + n_to_do = 0 + do k=1,elec_alpha_num + imo = mo_list_alpha_curr(k) + if ( imo /= mo_list_alpha_prev(k) ) then + n_to_do += 1 + to_do(n_to_do) = k + endif + enddo + + ddet = 0.d0 + + if (n_to_do < ishft(elec_alpha_num,1)) then + + do while ( n_to_do > 0 ) + ddet = det_alpha_value_curr + n_to_do_old = n_to_do + n_to_do = 0 + do l=1,n_to_do_old + k = to_do(l) + imo = mo_list_alpha_curr(k) + call det_update(elec_alpha_num, elec_alpha_num_8, & + mo_value(1,imo), & + k, & + slater_matrix_alpha, & + slater_matrix_alpha_inv_det, & + ddet) + if (ddet /= 0.d0) then + n_updates_det += 1.d0 + det_alpha_value_curr = ddet + else + n_to_do += 1 + to_do(n_to_do) = k + ddet = det_alpha_value_curr + endif + enddo + if (n_to_do == n_to_do_old) then + ddet = 0.d0 + exit + endif + enddo + + endif + + else + + ddet = 0.d0 + + endif + + ! Avoid NaN + if (ddet /= 0.d0) then + n_updated_det += 1.d0 + else + n_inverted_det += 1.d0 + do j=1,mo_closed_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(100) + do i=1,elec_alpha_num + slater_matrix_alpha(i,j) = mo_value(i,j) + slater_matrix_alpha_inv_det(j,i) = mo_value(i,j) + enddo + enddo + do k=mo_closed_num+1,elec_alpha_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(100) + do i=1,elec_alpha_num + slater_matrix_alpha(i,k) = mo_value(i,mo_list_alpha_curr(k)) + slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k)) + 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 +END_PROVIDER + + BEGIN_PROVIDER [ double precision, det_beta_value_curr ] +&BEGIN_PROVIDER [ real, slater_matrix_beta, (elec_beta_num_8,elec_beta_num) ] +&BEGIN_PROVIDER [ double precision, slater_matrix_beta_inv_det, (elec_beta_num_8,elec_beta_num) ] + BEGIN_DOC + ! det_beta_value_curr : Value of the current beta determinant + ! + ! slater_matrix_beta : Slater matrix for the current beta determinant. + ! 1st index runs over electrons and + ! 2nd index runs over MOs. + ! Built with 1st determinant + ! + ! 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 + + n_to_do = 0 + do k=mo_closed_num+1,elec_beta_num + imo = mo_list_beta_curr(k) + if ( imo /= mo_list_beta_prev(k) ) then + n_to_do += 1 + to_do(n_to_do) = k + endif + enddo + + ddet = 0.d0 + if (n_to_do < ishft(elec_beta_num,1)) then + + do while ( n_to_do > 0 ) + ddet = det_beta_value_curr + n_to_do_old = n_to_do + n_to_do = 0 + do l=1,n_to_do_old + k = to_do(l) + imo = mo_list_beta_curr(k) + call det_update(elec_beta_num, elec_beta_num_8, & + mo_value(elec_alpha_num+1,imo), & + k, & + slater_matrix_beta, & + slater_matrix_beta_inv_det, & + ddet) + if (ddet /= 0.d0) then + n_updates_det += 1.d0 + det_beta_value_curr = ddet + else + n_to_do += 1 + to_do(n_to_do) = k + ddet = det_beta_value_curr + endif + enddo + if (n_to_do == n_to_do_old) then + ddet = 0.d0 + exit + endif + enddo + + endif + + else + + ddet = 0.d0 + + endif + + ! Avoid NaN + if (ddet /= 0.d0) then + n_updated_det += 1.d0 + else + n_inverted_det += 1.d0 + do j=1,mo_closed_num + !DIR$ VECTOR UNALIGNED + !DIR$ LOOP COUNT (100) + do i=1,elec_beta_num + slater_matrix_beta(i,j) = mo_value(i+elec_alpha_num,j) + slater_matrix_beta_inv_det(j,i) = mo_value(i+elec_alpha_num,j) + enddo + enddo + do k=mo_closed_num+1,elec_beta_num + !DIR$ VECTOR UNALIGNED + !DIR$ LOOP COUNT (100) + do i=1,elec_beta_num + slater_matrix_beta(i,k) = mo_value(i+elec_alpha_num,mo_list_beta_curr(k)) + slater_matrix_beta_inv_det(k,i) = mo_value(i+elec_alpha_num,mo_list_beta_curr(k)) + enddo + enddo + 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 + +END_PROVIDER + + BEGIN_PROVIDER [ integer, det_alpha_num_pseudo ] +&BEGIN_PROVIDER [ integer, det_beta_num_pseudo ] + implicit none + BEGIN_DOC + ! Dimensioning of large arrays made smaller without pseudo + END_DOC + if (do_pseudo) then + det_alpha_num_pseudo = det_alpha_num + det_beta_num_pseudo = det_beta_num + else + det_alpha_num_pseudo = 1 + det_beta_num_pseudo = 1 + endif +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 + ifirst = 1 + det_alpha_value = 0.d0 + 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 + do i=1,elec_alpha_num + !$OMP SIMD + do k=1,4 + det_alpha_grad_lapl(k,i,det_i) = det_alpha_grad_lapl_curr(k,i) + enddo + !$OMP END SIMD + if (do_pseudo) then + det_alpha_pseudo(i,det_i) = det_alpha_pseudo_curr(i) + endif + enddo + + 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 + !DIR$ LOOP COUNT (200) + do i=elec_alpha_num+1,elec_num + !$OMP SIMD + do k=1,4 + det_beta_grad_lapl(k,i,det_j) = det_beta_grad_lapl_curr(k,i) + enddo + !$OMP END SIMD + if (do_pseudo) then + det_beta_pseudo(i,det_j) = det_beta_pseudo_curr(i) + endif + enddo + + enddo + + det_j = det_beta_order(1) + SOFT_TOUCH det_j + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, single_det_value ] +&BEGIN_PROVIDER [ double precision, single_det_grad, (elec_num_8,3) ] +&BEGIN_PROVIDER [ double precision, single_det_lapl, (elec_num) ] + BEGIN_DOC + ! Value of a single determinant wave function from the 1st determinant + END_DOC + det_i = 1 + det_j = 1 + integer :: i + single_det_value = det_alpha_value_curr * det_beta_value_curr + do i=1,elec_alpha_num + single_det_grad(i,1) = det_alpha_grad_lapl_curr(1,i) * det_beta_value_curr + single_det_grad(i,2) = det_alpha_grad_lapl_curr(2,i) * det_beta_value_curr + single_det_grad(i,3) = det_alpha_grad_lapl_curr(3,i) * det_beta_value_curr + single_det_lapl(i) = det_alpha_grad_lapl_curr(4,i) * det_beta_value_curr + enddo + do i=elec_alpha_num+1,elec_num + single_det_grad(i,1) = det_alpha_value_curr * det_beta_grad_lapl_curr(1,i) + single_det_grad(i,2) = det_alpha_value_curr * det_beta_grad_lapl_curr(2,i) + single_det_grad(i,3) = det_alpha_value_curr * det_beta_grad_lapl_curr(3,i) + single_det_lapl(i) = det_alpha_value_curr * det_beta_grad_lapl_curr(4,i) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, psidet_value ] +&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) + + integer :: i,j,k, l + integer :: i1,i2,i3,i4,det_num4 + integer :: j1,j2,j3,j4 + double precision :: f + + DaC = 0.d0 + CDb = 0.d0 + + det_num4 = iand(det_num,not(3)) + !DIR$ VECTOR ALIGNED + do k=1,det_num4,4 + i1 = det_coef_matrix_rows(k ) + i2 = det_coef_matrix_rows(k+1) + i3 = det_coef_matrix_rows(k+2) + i4 = det_coef_matrix_rows(k+3) + j1 = det_coef_matrix_columns(k ) + j2 = det_coef_matrix_columns(k+1) + j3 = det_coef_matrix_columns(k+2) + j4 = det_coef_matrix_columns(k+3) + if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then + f = det_beta_value (j1) + CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*f + CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f + CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f + CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f + + if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then + DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) & + + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) & + + det_coef_matrix_values(k+2)*det_alpha_value(i1+2) & + + det_coef_matrix_values(k+3)*det_alpha_value(i1+3) + else + DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) & + + det_coef_matrix_values(k+1)*det_alpha_value(i2) & + + det_coef_matrix_values(k+2)*det_alpha_value(i3) & + + det_coef_matrix_values(k+3)*det_alpha_value(i4) + endif + else + DaC(j1) = DaC(j1) + det_coef_matrix_values(k )*det_alpha_value(i1) + DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2) + DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3) + DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4) + CDb(i1) = CDb(i1) + det_coef_matrix_values(k )*det_beta_value (j1) + CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2) + CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3) + CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4) + endif + enddo + + do k=det_num4+1,det_num + i = det_coef_matrix_rows(k) + j = det_coef_matrix_columns(k) + DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i) + CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j) + enddo + + ! Value + ! ----- + + psidet_value = 0.d0 + !$OMP SIMD reduction(+:psidet_value) + do j=1,det_beta_num + psidet_value = psidet_value + det_beta_value(j) * DaC(j) + enddo + !$OMP END SIMD + + + if (psidet_value == 0.d0) then + call abrt(irp_here,'Determinantal component of the wave function is zero.') + endif + psidet_inv = 1.d0/psidet_value + + ! Gradients + ! --------- + + if (do_pseudo) then + do j=1,elec_num + psidet_grad_lapl(1:4,j) = 0.d0 + pseudo_non_local(j) = 0.d0 + enddo + + do i=1,det_alpha_num + do j=1,elec_alpha_num + !$OMP SIMD + do k=1,4 + psidet_grad_lapl(k,j) = psidet_grad_lapl(k,j) + det_alpha_grad_lapl(k,j,i)*CDb(i) + enddo + !$OMP END SIMD + pseudo_non_local(j) = pseudo_non_local(j) + det_alpha_pseudo(j,i)*CDb(i) + enddo + enddo + + do i=1,det_beta_num + do j=elec_alpha_num+1,elec_num + !$OMP SIMD + do k=1,4 + psidet_grad_lapl(k,j) = psidet_grad_lapl(k,j) + det_beta_grad_lapl(k,j,i)*DaC(i) + enddo + !$OMP END SIMD + pseudo_non_local(j) = pseudo_non_local(j) + det_beta_pseudo(j,i)*DaC(i) + enddo + enddo + + do j=1,elec_num + pseudo_non_local(j) = pseudo_non_local(j) * psidet_inv + enddo + + else + + do j=1,elec_num + psidet_grad_lapl(1:4,j) = 0.d0 + enddo + + do i=1,det_alpha_num + do j=1,elec_alpha_num + !$OMP SIMD + do k=1,4 + psidet_grad_lapl(k,j) = psidet_grad_lapl(k,j) + det_alpha_grad_lapl(k,j,i)*CDb(i) + enddo + !$OMP END SIMD + enddo + enddo + + do i=1,det_beta_num + do j=elec_alpha_num+1,elec_num + !$OMP SIMD + do k=1,4 + psidet_grad_lapl(k,j) = psidet_grad_lapl(k,j) + det_beta_grad_lapl(k,j,i)*DaC(i) + enddo + !$OMP END SIMD + enddo + enddo + + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, det_alpha_pseudo_curr, (elec_alpha_num) ] + implicit none + BEGIN_DOC +! Pseudopotential non-local contribution + END_DOC + integer :: i,j,l,m,k,n + integer :: imo,kk + double precision :: c + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + det_alpha_pseudo_curr = 0.d0 + endif + if (do_pseudo) then + do i=1,elec_alpha_num + det_alpha_pseudo_curr(i) = 0.d0 + do n=1,elec_alpha_num + imo = mo_list_alpha_curr(n) + c = slater_matrix_alpha_inv_det(i,n) + det_alpha_pseudo_curr(i) = & + det_alpha_pseudo_curr(i) + c*pseudo_mo_term(imo,i) + enddo + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, det_beta_pseudo_curr, (elec_alpha_num+1:elec_num) ] + implicit none + BEGIN_DOC +! Pseudopotential non-local contribution + END_DOC + integer :: i,j,l,m,k,n + integer :: imo,kk + double precision :: c + integer, save :: ifirst = 0 + if (elec_beta_num == 0) then + return + endif + if (ifirst == 0) then + ifirst = 1 + det_beta_pseudo_curr = 0.d0 + endif + if (do_pseudo) then + do i=elec_alpha_num+1,elec_num + det_beta_pseudo_curr(i) = 0.d0 + do n=1,elec_beta_num + imo = mo_list_beta_curr(n) + c = slater_matrix_beta_inv_det(i-elec_alpha_num,n) + det_beta_pseudo_curr(i) = & + det_beta_pseudo_curr(i) + c*pseudo_mo_term(imo,i) + enddo + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num) ] + implicit none + BEGIN_DOC + ! Gradient of the current alpha determinant + END_DOC + + integer :: i, j, k + !DIR$ VECTOR ALIGNED + do i=1,elec_alpha_num + det_alpha_grad_lapl_curr(1,i) = 0.d0 + det_alpha_grad_lapl_curr(2,i) = 0.d0 + det_alpha_grad_lapl_curr(3,i) = 0.d0 + det_alpha_grad_lapl_curr(4,i) = 0.d0 + enddo + + integer :: imo, imo2 + +! ------- +! The following code does the same as this: +! +! do j=1,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(k,i,imo)*slater_matrix_alpha_inv_det(i,j) +! enddo +! enddo +! enddo +! +! ------- + + if (iand(elec_alpha_num,1) == 0) then + + do j=1,elec_alpha_num,2 + imo = mo_list_alpha_curr(j ) + imo2 = mo_list_alpha_curr(j+1) + do i=1,elec_alpha_num,2 + !$OMP SIMD + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl(k,i ,imo )*slater_matrix_alpha_inv_det(i ,j ) & + + mo_grad_lapl(k,i ,imo2)*slater_matrix_alpha_inv_det(i ,j+1) + det_alpha_grad_lapl_curr(k,i+1) = det_alpha_grad_lapl_curr(k,i+1) + mo_grad_lapl(k,i+1,imo )*slater_matrix_alpha_inv_det(i+1,j ) & + + mo_grad_lapl(k,i+1,imo2)*slater_matrix_alpha_inv_det(i+1,j+1) + enddo + !$OMP END SIMD + enddo + enddo + + else + + do j=1,elec_alpha_num-1,2 + imo = mo_list_alpha_curr(j ) + imo2 = mo_list_alpha_curr(j+1) + do i=1,elec_alpha_num-1,2 + !$OMP SIMD + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl(k,i ,imo )*slater_matrix_alpha_inv_det(i ,j ) & + + mo_grad_lapl(k,i ,imo2)*slater_matrix_alpha_inv_det(i ,j+1) + det_alpha_grad_lapl_curr(k,i+1) = det_alpha_grad_lapl_curr(k,i+1) + mo_grad_lapl(k,i+1,imo )*slater_matrix_alpha_inv_det(i+1,j ) & + + mo_grad_lapl(k,i+1,imo2)*slater_matrix_alpha_inv_det(i+1,j+1) + enddo + !$OMP END SIMD + enddo + i=elec_alpha_num + !$OMP SIMD + do k=1,4 + det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl(k,i,imo )*slater_matrix_alpha_inv_det(i,j ) & + + mo_grad_lapl(k,i,imo2)*slater_matrix_alpha_inv_det(i,j+1) + enddo + !$OMP END SIMD + enddo + + j=elec_alpha_num + imo = mo_list_alpha_curr(j) + do i=1,elec_alpha_num + !$OMP SIMD + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) + enddo + !$OMP END SIMD + enddo + + endif + + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1:elec_num) ] + implicit none + BEGIN_DOC + ! Gradient and Laplacian of the current beta determinant + END_DOC + + integer :: i, j, k, l + + !DIR$ VECTOR ALIGNED + do i=elec_alpha_num+1,elec_num + det_beta_grad_lapl_curr(1,i) = 0.d0 + det_beta_grad_lapl_curr(2,i) = 0.d0 + det_beta_grad_lapl_curr(3,i) = 0.d0 + det_beta_grad_lapl_curr(4,i) = 0.d0 + enddo + + integer :: imo, imo2 + +! ------- +! The following code does the same as this: +! +! do j=1,elec_beta_num +! imo = mo_list_beta_curr(j) +! do i=elec_alpha_num+1,elec_num +! do k=1,4 +! det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& +! mo_grad_lapl(k,i,imo)*slater_matrix_beta_inv_det(i-elec_alpha_num,j) +! enddo +! enddo +! enddo +! +! -------- + + if (iand(elec_beta_num,1) == 0) then + + do j=1,elec_beta_num,2 + imo = mo_list_beta_curr(j ) + imo2 = mo_list_beta_curr(j+1) + !DIR$ LOOP COUNT (16) + do i=elec_alpha_num+1,elec_num,2 + l = i-elec_alpha_num + !$OMP SIMD + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl(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(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & + mo_grad_lapl(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + enddo + !$OMP END SIMD + enddo + enddo + + else + + do j=1,elec_beta_num-1,2 + imo = mo_list_beta_curr(j ) + imo2 = mo_list_beta_curr(j+1) + !DIR$ LOOP COUNT (16) + do i=elec_alpha_num+1,elec_num-1,2 + l = i-elec_alpha_num + !$OMP SIMD + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl(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(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & + mo_grad_lapl(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + enddo + !$OMP END SIMD + enddo + i = elec_num + l = elec_num-elec_alpha_num + !$OMP SIMD + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + enddo + !$OMP END SIMD + enddo + + j = elec_beta_num + imo = mo_list_beta_curr(j) + !DIR$ LOOP COUNT (16) + do i=elec_alpha_num+1,elec_num + l = i-elec_alpha_num + !$OMP SIMD + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl(k,i,imo)*slater_matrix_beta_inv_det(l,j) + enddo + !$OMP END SIMD + enddo + + endif + +END_PROVIDER + + diff --git a/src/electrons.irp.f b/src/electrons.irp.f new file mode 100644 index 0000000..f120b39 --- /dev/null +++ b/src/electrons.irp.f @@ -0,0 +1,333 @@ +BEGIN_PROVIDER [ double precision, xbrown, (elec_num_8,3) ] + BEGIN_DOC + ! Brownian step. Built in Brownian_step subroutine. + END_DOC + integer :: i,l + xbrown = 0.d0 +END_PROVIDER + + + BEGIN_PROVIDER [ integer, elec_alpha_num ] +&BEGIN_PROVIDER [ integer, elec_alpha_num_8 ] + implicit none + BEGIN_DOC + ! Number of alpha electrons + END_DOC + integer, external :: mod_align + elec_alpha_num = -1 + call get_electrons_elec_alpha_num(elec_alpha_num) + if (elec_alpha_num <= 0) then + call abrt(irp_here,'Number of alpha electrons should be > 0') + endif + elec_alpha_num_8 = mod_align(elec_alpha_num) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer, elec_beta_num ] +&BEGIN_PROVIDER [ integer, elec_beta_num_8 ] + implicit none + BEGIN_DOC + ! Number of beta electrons + END_DOC + integer, external :: mod_align + elec_beta_num = 0 + call get_electrons_elec_beta_num(elec_beta_num) + if (elec_beta_num < 0) then + call abrt(irp_here,'Number of beta electrons should be >= 0') + endif + elec_beta_num_8 = mod_align(elec_beta_num) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer, elec_num ] +&BEGIN_PROVIDER [ integer, elec_num_8 ] +&BEGIN_PROVIDER [ integer, elec_num_1_8 ] + implicit none + BEGIN_DOC + ! Number of electrons + END_DOC + integer, external :: mod_align + elec_num = elec_alpha_num + elec_beta_num + ASSERT ( elec_num > 0 ) + elec_num_8 = mod_align(elec_num) + elec_num_1_8 = mod_align(elec_num+1) + +END_PROVIDER + + + +BEGIN_PROVIDER [ real, elec_coord_full, (elec_num+1,3,walk_num) ] + implicit none + BEGIN_DOC + ! Electron coordinates of all walkers + ! Component (elec_num+1,1,walk_num) contains the length realized by the walker. + ! Initialized in init_walkers + END_DOC + integer :: i,k + real, allocatable :: buffer2(:,:,:) + if ( is_worker ) then + + call get_elec_coord_full(elec_coord_full,size(elec_coord_full,1)) + + else + + if (.not.do_prepare) then + allocate ( buffer2(elec_num+1,3,elec_coord_pool_size) ) + call get_electrons_elec_coord_pool(buffer2) + do k=1,walk_num + do i=1,elec_num+1 + elec_coord_full(i,1,k) = buffer2(i,1,k) + elec_coord_full(i,2,k) = buffer2(i,2,k) + elec_coord_full(i,3,k) = buffer2(i,3,k) + enddo + enddo + deallocate ( buffer2 ) + else + elec_coord_full = 0. + endif + + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, elec_coord_pool_size ] + implicit none + BEGIN_DOC + ! Size of the pool of electron coordinates + END_DOC + elec_coord_pool_size = walk_num_tot + call get_electrons_elec_coord_pool_size(elec_coord_pool_size) + call iinfo(irp_here,'elec_coord_pool_size',elec_coord_pool_size) + +END_PROVIDER + + + +BEGIN_PROVIDER [ real, elec_coord, (elec_num_1_8,3) ] + implicit none + BEGIN_DOC + ! Electron coordinates + END_DOC + integer :: i, k + elec_coord = 0. + do k=1,3 + do i=1,elec_num+1 + elec_coord(i,k) = elec_coord_full(i,k,walk_i) + enddo + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [ real, elec_coord_transp, (8,elec_num) + implicit none + BEGIN_DOC + ! Transposed array of elec_coord + END_DOC + integer :: i, k + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + elec_coord_transp = 0. + endif + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do i=1,elec_num + elec_coord_transp(1,i) = elec_coord(i,1) + elec_coord_transp(2,i) = elec_coord(i,2) + elec_coord_transp(3,i) = elec_coord(i,3) + enddo +END_PROVIDER + + + + BEGIN_PROVIDER [ real, elec_dist, (elec_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, elec_dist_vec_x, (elec_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, elec_dist_vec_y, ((-simd_sp+1):elec_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, elec_dist_vec_z, ((-2*simd_sp+1):elec_num_8,elec_num) ] + implicit none + BEGIN_DOC + ! Electron-electron distances + END_DOC + integer :: ie1, ie2, l + + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + !DIR$ VECTOR ALIGNED + elec_dist = 0. + !DIR$ VECTOR ALIGNED + elec_dist_vec_x = 0. + !DIR$ VECTOR ALIGNED + elec_dist_vec_y = 0. + !DIR$ VECTOR ALIGNED + elec_dist_vec_z = 0. + endif + + do ie2 = 1,elec_num + real :: x, y, z + x = elec_coord(ie2,1) + y = elec_coord(ie2,2) + z = elec_coord(ie2,3) + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do ie1 = 1,elec_num + elec_dist_vec_x(ie1,ie2) = elec_coord(ie1,1) - x + elec_dist_vec_y(ie1,ie2) = elec_coord(ie1,2) - y + elec_dist_vec_z(ie1,ie2) = elec_coord(ie1,3) - z + enddo + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do ie1 = 1,elec_num + elec_dist(ie1,ie2) = sqrt( & + elec_dist_vec_x(ie1,ie2)*elec_dist_vec_x(ie1,ie2) + & + elec_dist_vec_y(ie1,ie2)*elec_dist_vec_y(ie1,ie2) + & + elec_dist_vec_z(ie1,ie2)*elec_dist_vec_z(ie1,ie2) ) + enddo + enddo + +END_PROVIDER + + + + BEGIN_PROVIDER [ real, nucl_elec_dist, (nucl_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, nucl_elec_dist_vec, (3,nucl_num,elec_num) ] + implicit none + BEGIN_DOC + ! Electron-nucleus distances |r_elec - R_nucl| + END_DOC + integer :: i,j,l + integer, save :: ifirst = 0 + + if (ifirst == 0) then + ifirst = 1 + !DIR$ VECTOR ALIGNED + nucl_elec_dist = 0. + !DIR$ VECTOR ALIGNED + nucl_elec_dist_vec = 0. + endif + + do i = 1,elec_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j = 1,nucl_num + nucl_elec_dist_vec(1,j,i) = elec_coord_transp(1,i) - nucl_coord(j,1) + nucl_elec_dist_vec(2,j,i) = elec_coord_transp(2,i) - nucl_coord(j,2) + nucl_elec_dist_vec(3,j,i) = elec_coord_transp(3,i) - nucl_coord(j,3) + enddo + enddo + + do i = 1,elec_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j = 1,nucl_num + nucl_elec_dist(j,i) = (elec_coord(i,1) - nucl_coord(j,1)) & + * (elec_coord(i,1) - nucl_coord(j,1)) & + + (elec_coord(i,2) - nucl_coord(j,2)) & + * (elec_coord(i,2) - nucl_coord(j,2)) & + + (elec_coord(i,3) - nucl_coord(j,3)) & + * (elec_coord(i,3) - nucl_coord(j,3)) + nucl_elec_dist(j,i) = max(1.e-6,sqrt(nucl_elec_dist(j,i))) + enddo + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [ integer, elec_num_2, (2) ] + + BEGIN_DOC + ! Number of alpha and beta electrons in an array + END_DOC + + elec_num_2(1) = elec_alpha_num + elec_num_2(2) = elec_beta_num + +END_PROVIDER + + +BEGIN_PROVIDER [ integer, elec_spin, (elec_num) ] + implicit none + BEGIN_DOC + ! Electron spin. +1 for alpha and -1 for beta + END_DOC + integer :: i + do i=1,elec_alpha_num + elec_spin(i) = 1 + enddo + do i=elec_alpha_num+1,elec_num + elec_spin(i) = -1 + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [ real, elec_dist_inv, (elec_num_8,elec_num) ] + implicit none + BEGIN_DOC + ! 1/rij matrix + END_DOC + integer :: i,j + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + elec_dist_inv = 0. + endif + do i=1,elec_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do j=1,elec_num + elec_dist_inv(j,i) = 1./(elec_dist(j,i)+1.e-12) + enddo + elec_dist_inv(i,i) = 0. + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ real, nucl_elec_dist_inv, (nucl_num_8,elec_num) ] + implicit none + BEGIN_DOC + ! 1/rij matrix + END_DOC + integer :: i,j + do j=1,elec_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do i=1,nucl_num + nucl_elec_dist_inv(i,j) = 1./nucl_elec_dist(i,j) + enddo + enddo +END_PROVIDER + + +subroutine save_elec_coord_full + implicit none + BEGIN_DOC +! Save the electron coordinates to disk + END_DOC + integer :: i,k,l + real, allocatable :: buffer2(:,:,:) + + allocate ( buffer2(elec_num+1,3,elec_coord_pool_size) ) + k=0 + do l=1,elec_coord_pool_size + k = k+1 + if (k == walk_num+1) then + k=1 + endif + do i=1,elec_num+1 + buffer2(i,1,l) = elec_coord_full(i,1,k) + buffer2(i,2,l) = elec_coord_full(i,2,k) + buffer2(i,3,l) = elec_coord_full(i,3,k) + enddo + enddo + call ezfio_set_electrons_elec_coord_pool(buffer2) + deallocate ( buffer2 ) + +end + diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f new file mode 100644 index 0000000..c225411 --- /dev/null +++ b/src/ezfio_interface.irp.f @@ -0,0 +1,182 @@ +BEGIN_SHELL [ /usr/bin/python ] + +data = [ \ +("electrons_elec_coord_pool_size" , "integer" , "" ), +("electrons_elec_coord_pool" , "real" , "(elec_num+1,3,elec_coord_pool_size)" ), +("nuclei_nucl_num" , "integer" , "" ), +("nuclei_nucl_charge" , "real" , "(nucl_num)" ), +("nuclei_nucl_coord" , "real" , "(nucl_num,3)" ), +("nuclei_nucl_fitcusp_radius" , "real" , "(nucl_num)" ), +("mo_basis_mo_coef" , "real" , "(ao_num,mo_tot_num)" ), +("electrons_elec_fitcusp_radius" , "real" , "" ), +("electrons_elec_alpha_num" , "integer" , "" ), +("electrons_elec_beta_num" , "integer" , "" ), +("electrons_elec_walk_num" , "integer" , "" ), +("electrons_elec_walk_num_tot" , "integer" , "" ), +("ao_basis_ao_num" , "integer" , "" ), +("ao_basis_ao_prim_num" , "integer" , "(ao_num)" ), +("ao_basis_ao_nucl" , "integer" , "(ao_num)" ), +("ao_basis_ao_power" , "integer" , "(ao_num,3)" ), +("ao_basis_ao_expo" , "real" , "(ao_num,ao_prim_num_max)" ), +("ao_basis_ao_coef" , "real" , "(ao_num,ao_prim_num_max)" ), +("jastrow_jast_a_up_up" , "real" , "" ), +("jastrow_jast_a_up_dn" , "real" , "" ), +("jastrow_jast_b_up_up" , "real" , "" ), +("jastrow_jast_b_up_dn" , "real" , "" ), +("jastrow_jast_pen" , "real" , "(nucl_num)" ), +("jastrow_jast_eeN_e_a" , "real" , "" ), +("jastrow_jast_eeN_e_b" , "real" , "" ), +("jastrow_jast_eeN_N" , "real" , "(nucl_num)" ), +("jastrow_jast_core_a1" , "real" , "(nucl_num)" ), +("jastrow_jast_core_a2" , "real" , "(nucl_num)" ), +("jastrow_jast_core_b1" , "real" , "(nucl_num)" ), +("jastrow_jast_core_b2" , "real" , "(nucl_num)" ), +("jastrow_jast_type" , "character*(32)", "" ), +("simulation_stop_time" , "integer" , "" ), +("simulation_equilibration" , "logical" , "" ), +("simulation_block_time" , "integer" , "" ), +("simulation_time_step" , "real" , "" ), +("simulation_method" , "character*(32)", "" ), +("simulation_save_data" , "logical" , "" ), +("simulation_print_level" , "integer" , "" ), +("simulation_do_nucl_fitcusp" , "logical" , "" ), +("simulation_sampling" , "character*(32)", "" ), +("simulation_ci_threshold" , "double precision" , "" ), +("simulation_http_server" , "character*(128)", "" ), +("simulation_md5_key" , "character*(32)" , "" ), +("simulation_e_ref" , "double precision" , "" ), +("simulation_do_run" , "logical " , "" ), +("pseudo_do_pseudo" , "logical " , "" ), + +] + +data_no_set = [\ +("mo_basis_mo_tot_num" , "integer" , ""), +("mo_basis_mo_active_num" , "integer" , ""), +("mo_basis_mo_closed_num" , "integer" , ""), +("pseudo_ao_pseudo_grid" , "double precision" , "(ao_num,pseudo_lmax+pseudo_lmax+1,pseudo_lmax-0+1,nucl_num,pseudo_grid_size)"), +("pseudo_mo_pseudo_grid" , "double precision" , "(ao_num,pseudo_lmax+pseudo_lmax+1,pseudo_lmax-0+1,nucl_num,pseudo_grid_size)"), +("pseudo_pseudo_dz_k" , "double precision" , "(nucl_num,pseudo_klocmax)"), +("pseudo_pseudo_dz_kl" , "double precision" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"), +("pseudo_pseudo_grid_rmax" , "double precision" , ""), +("pseudo_pseudo_grid_size" , "integer" , ""), +("pseudo_pseudo_klocmax" , "integer" , ""), +("pseudo_pseudo_kmax" , "integer" , ""), +("pseudo_pseudo_lmax" , "integer" , ""), +("pseudo_pseudo_n_k" , "integer" , "(nucl_num,pseudo_klocmax)"), +("pseudo_pseudo_n_kl" , "integer" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"), +("pseudo_pseudo_v_k" , "double precision" , "(nucl_num,pseudo_klocmax)"), +("pseudo_pseudo_v_kl" , "double precision" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"), +("spindeterminants_n_det_alpha" , "integer" , ""), +("spindeterminants_n_det_beta" , "integer" , ""), +("spindeterminants_n_det" , "integer" , ""), +("spindeterminants_n_int" , "integer" , ""), +("spindeterminants_bit_kind" , "integer" , ""), +("spindeterminants_n_states" , "integer" , ""), +("spindeterminants_psi_det_alpha" , "integer*8" , "(N_int*bit_kind/8,det_alpha_num)"), +("spindeterminants_psi_det_beta" , "integer*8" , "(N_int*bit_kind/8,det_beta_num)"), +("spindeterminants_psi_coef_matrix_rows" , "integer" , "(det_num_input)"), +("spindeterminants_psi_coef_matrix_columns" , "integer" , "(det_num_input)"), +("spindeterminants_psi_coef_matrix_values" , "double precision" , "(det_num_input,N_states)"), + + +] + +def do_subst(t0,d): + t = t0 + t = t.replace("$X",d[0]) + t = t.replace("$T",d[1]) + t = t.replace("$D",d[2]) + if d[1].startswith("character"): + size = d[1].split("*")[1][1:-1] + u = "character" + elif d[1].startswith("double precision"): + u = d[1].replace(" ","_") + size = "1" + elif "*" in d[1]: + size = "1" + u = d[1].replace("*","") + else: + size = "1" + u = d[1] + t = t.replace("$U",u) + if d[2] == "": + t = t.replace("$S",size) + else: + if size == "1": + t = t.replace("$S","size(res)") + else: + t = t.replace("$S","%s*size(res)"%(size)) + provide = "" + tmp = d[2].replace('(','').replace(')','') + for i in "+-*/": + tmp = tmp.replace(i,',') + for i in tmp.split(','): + if ":" in i: + i = i.split(':')[1] + try: + eval(i+"+1") + except NameError: + provide += " PROVIDE "+i+"\n" + t = t.replace("$P",provide) + print t + +t0 = """ +subroutine get_$X(res) + implicit none + BEGIN_DOC +! Calls EZFIO subroutine to get $X + END_DOC + $T :: res$D + integer :: ierr, i + logical :: exists + PROVIDE ezfio_filename + $P + if (.not.is_worker) then + call ezfio_has_$X(exists) + if (exists) then + call ezfio_get_$X(res) + call ezfio_free_$X + else + call ezfio_set_$X(res) + endif + else + call zmq_ezfio_has('$X',exists) + if (exists) then + call zmq_ezfio_get_$U('$X',res,$S) + endif + endif + +end +""" + + +t1 = """ +subroutine get_$X(res) + implicit none + BEGIN_DOC +! Calls EZFIO subroutine to get $X + END_DOC + $T :: res$D + integer :: ierr + PROVIDE ezfio_filename + $P + if (.not.is_worker) then + call ezfio_get_$X(res) + call ezfio_free_$X + else + call zmq_ezfio_get_$U('$X',res,$S) + endif + +end +""" + +for i,d in enumerate(data): + do_subst(t0,d) + +for i,d in enumerate(data_no_set): + i += len(data) + do_subst(t1,d) + +END_SHELL + diff --git a/src/finish.irp.f b/src/finish.irp.f new file mode 100644 index 0000000..0d41eea --- /dev/null +++ b/src/finish.irp.f @@ -0,0 +1,22 @@ +subroutine abrt (here,message) + implicit none + character*(*) :: here + character*(*) :: message + write(0,*) '-------------------------' + write(0,*) 'Error in '//trim(here)//':' + write(0,*) '-------------------------' + write(0,*) trim(message)//'.' + write(0,*) '-------------------------' + if (is_worker) then + call worker_log(here,message) + call sleep(2) + endif + call finish() + stop +end + +subroutine finish() + implicit none + call ezfio_finish() +end + diff --git a/src/mo.irp.f b/src/mo.irp.f new file mode 100644 index 0000000..233137a --- /dev/null +++ b/src/mo.irp.f @@ -0,0 +1,693 @@ + BEGIN_PROVIDER [ integer, mo_num ] +&BEGIN_PROVIDER [ integer, mo_num_8 ] + implicit none + BEGIN_DOC +! 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 + + +BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ] + implicit none + BEGIN_DOC +! Molecular orbital coefficients read from the input file + END_DOC + 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 + do j=1,ao_num + mo_coef_input(j,i) = buffer(j,i) + enddo + do j=ao_num+1,ao_num_8 + mo_coef_input(j,i) = 0. + enddo + call set_order(mo_coef_input(1,i),ao_nucl_sort_idx,ao_num) + enddo + deallocate(buffer) + +END_PROVIDER + + + BEGIN_PROVIDER [ real, mo_scale ] +&BEGIN_PROVIDER [ real, mo_norm ] + implicit none + BEGIN_DOC +! Scaling factor for MOs to keep the determinant in a defined domain + END_DOC + mo_scale = 1.d0/(0.4d0*log(float(elec_num+1))) + mo_norm = mo_scale*mo_scale +END_PROVIDER + + +BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ] + implicit none + BEGIN_DOC +! 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) + enddo + enddo + do j =mo_num+1,mo_num_8 + !DIR$ VECTOR ALIGNED + do i=1,ao_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 + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (2000) + do i=1,ao_num_8 + mo_coef(i,j) *= f + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num_8,ao_num_8) ] + implicit none + BEGIN_DOC +! 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 + + + BEGIN_PROVIDER [ integer, mo_coef_transp_non_zero_idx, (0:mo_num,ao_num) ] +&BEGIN_PROVIDER [ real, mo_coef_transp_sparsity ] + implicit none + BEGIN_DOC +! Indices of the non-zero elements of the transpose of the Molecular +! orbital coefficients + END_DOC + integer :: i, j + + integer :: idx + mo_coef_transp_sparsity = 0. + do j=1,ao_num + idx = 0 + do i=1,mo_num + if (mo_coef_transp(i,j) /= 0.) then + idx += 1 + mo_coef_transp_non_zero_idx(idx,j) = i + endif + enddo + mo_coef_transp_non_zero_idx(0,j) = idx + mo_coef_transp_sparsity += float(idx) + enddo + mo_coef_transp_sparsity *= 1./(mo_num*ao_num) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, mo_coef_transp_present, (num_present_mos_8,ao_num_8) ] + implicit none + BEGIN_DOC +! mo_coef_transp without MOs absent in all determinants + END_DOC + integer :: i,j,n + mo_coef_transp_present = 0. + do i=1,ao_num + do j=1,num_present_mos + mo_coef_transp_present(j,i) = mo_coef_transp(present_mos(j),i) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ real, mo_value_transp, ((-simd_sp+1):mo_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, mo_grad_transp_x, ((-2*simd_sp+1):mo_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, mo_grad_transp_y, ((-3*simd_sp+1):mo_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, mo_grad_transp_z, ((-4*simd_sp+1):mo_num_8,elec_num) ] +&BEGIN_PROVIDER [ real, mo_lapl_transp, ((-5*simd_sp+1):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 + TOUCH ao_elec + endif + if (num_present_mos == mo_num) then + call sparse_full_mv(mo_coef_transp,mo_num_8, & + ao_value_block(1),ao_num_8, & + ao_grad_block_x(1), & + ao_grad_block_y(1), & + ao_grad_block_z(1), & + ao_lapl_block(1), & + ao_value_non_zero_idx(0), & + mo_value_transp(1,i),mo_num_8, & + mo_grad_transp_x(1,i), & + mo_grad_transp_y(1,i), & + 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, & + ao_grad_block_x(1), & + ao_grad_block_y(1), & + ao_grad_block_z(1), & + ao_lapl_block(1), & + ao_value_non_zero_idx(0), & + mo_value_transp(1,i),mo_num_8, & + mo_grad_transp_x(1,i), & + mo_grad_transp_y(1,i), & + 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) + mo_grad_transp_y(present_mos(j),i) = mo_grad_transp_y(j,i) + mo_grad_transp_z(present_mos(j),i) = mo_grad_transp_z(j,i) + mo_lapl_transp (present_mos(j),i) = mo_lapl_transp (j,i) + if (present_mos(j) == j) then + exit + 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 + cycle + endif + r_inv = nucl_elec_dist_inv(k,i) + !DIR$ LOOP COUNT (500) + do j=1,mo_num + mo_value_transp(j,i) = mo_value_transp(j,i) + nucl_fitcusp_param(1,j,k) +& + r * (nucl_fitcusp_param(2,j,k) + & + r * (nucl_fitcusp_param(3,j,k) + & + r * nucl_fitcusp_param(4,j,k) )) + mo_lapl_transp(j,i) = mo_lapl_transp(j,i) + & + nucl_fitcusp_param(2,j,k)*(r_inv+r_inv) + & + 6.*nucl_fitcusp_param(3,j,k) + & + r * 12.*nucl_fitcusp_param(4,j,k) + c = r_inv * (nucl_fitcusp_param(2,j,k) + & + r * (2.*nucl_fitcusp_param(3,j,k) + & + r * 3.*nucl_fitcusp_param(4,j,k) )) + mo_grad_transp_x(j,i) = mo_grad_transp_x(j,i) + nucl_elec_dist_vec(1,k,i)*c + mo_grad_transp_y(j,i) = mo_grad_transp_y(j,i) + nucl_elec_dist_vec(2,k,i)*c + mo_grad_transp_z(j,i) = mo_grad_transp_z(j,i) + nucl_elec_dist_vec(3,k,i)*c + 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 + t = prepare_walkers_t + do i=1,mo_num + !DIR$ LOOP COUNT (100) + do j=1,elec_alpha_num + if (i /= j) then + mo_value_transp(i,j) *= t + mo_grad_transp_x(i,j) *= t + mo_grad_transp_y(i,j) *= t + mo_grad_transp_z(i,j) *= t + mo_lapl_transp(i,j) *= t + endif + enddo + !DIR$ LOOP COUNT (100) + do j=1,elec_beta_num + if (i /= j) then + mo_value_transp(i,j+elec_alpha_num) *= t + mo_grad_transp_x(i,j+elec_alpha_num) *= t + mo_grad_transp_y(i,j+elec_alpha_num) *= t + mo_grad_transp_z(i,j+elec_alpha_num) *= t + mo_lapl_transp(i,j+elec_alpha_num) *= t + endif + enddo + enddo + endif + +END_PROVIDER + + +BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ] + implicit none + 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 + !DIR$ VECTOR ALIGNED + mo_value = 0. + endif + call transpose(mo_value_transp(1,1),mo_num_8+simd_sp,mo_value,elec_num_8,mo_num,elec_num) + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, mo_grad_x, (elec_num_8,mo_num) ] +&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 + !DIR$ VECTOR ALIGNED + mo_grad_y = 0.d0 + !DIR$ VECTOR ALIGNED + mo_grad_z = 0.d0 + ifirst = 1 + PROVIDE primitives_reduced + endif + ! Transpose x last for cache efficiency + call transpose_to_dp(mo_grad_transp_y(1,1),mo_num_8+3*simd_sp,mo_grad_y(1,1),elec_num_8,mo_num,elec_num) + call transpose_to_dp(mo_grad_transp_z(1,1),mo_num_8+4*simd_sp,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+2*simd_sp,mo_grad_x(1,1),elec_num_8,mo_num,elec_num) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_grad_lapl, (4,elec_num,mo_num) ] + implicit none + BEGIN_DOC +! Gradients and laplacian + END_DOC + integer :: i,j + do j=1,mo_num + do i=1,elec_num + mo_grad_lapl(1,i,j) = mo_grad_x(i,j) + mo_grad_lapl(2,i,j) = mo_grad_y(i,j) + mo_grad_lapl(3,i,j) = mo_grad_z(i,j) + mo_grad_lapl(4,i,j) = mo_lapl (i,j) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ] + implicit none + BEGIN_DOC +! Laplacians of the molecular orbitals + END_DOC + + integer :: i,j + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + PROVIDE primitives_reduced + !DIR$ VECTOR ALIGNED + mo_lapl = 0.d0 + endif + call transpose_to_dp(mo_lapl_transp(1,1),mo_num_8+5*simd_sp,mo_lapl,elec_num_8,mo_num,elec_num) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, prepare_walkers_t ] + implicit none + BEGIN_DOC +! prepare_walkers_t : scaling of the off-diagonal elements +! of the mo_value matrix + END_DOC + prepare_walkers_t = 1. +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_tot_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 + + +!----------------- +! Fit cusp +!----------------- + + +BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ] + implicit none + BEGIN_DOC +! Values of the molecular orbitals at the nucleus without the +! S components of the current nucleus + END_DOC + integer :: i, j, k, l + real :: ao_value_at_nucl_no_S(ao_num) + + 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 + + !DIR$ LOOP COUNT (2000) + do i=1,ao_num + if (ao_nucl(i) /= k) then + ao_value_at_nucl_no_S(i) = ao_value_p(i) + else + ao_value_at_nucl_no_S(i) = 0. + endif + enddo + + integer :: jj + do jj=1,num_present_mos + j = present_mos(jj) + mo_value_at_nucl(j,k) = 0. + !DIR$ VECTOR ALIGNED + do i=1,ao_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 + + + BEGIN_PROVIDER [ double precision, ao_value_at_fitcusp_radius, (ao_num_8,nucl_num) ] +&BEGIN_PROVIDER [ double precision, ao_grad_at_fitcusp_radius, (ao_num_8,nucl_num) ] +&BEGIN_PROVIDER [ double precision, ao_lapl_at_fitcusp_radius, (ao_num_8,nucl_num) ] + implicit none + 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 + + !DIR$ LOOP COUNT (2000) + 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) + ao_lapl_at_fitcusp_radius(j,k) = ao_lapl_p(j) + if ( (ao_nucl(j) /= k).or.(ao_power(j,4) >0) ) then + ao_value_at_fitcusp_radius(j,k) = 0. + ao_grad_at_fitcusp_radius(j,k) = 0. + ao_lapl_at_fitcusp_radius(j,k) = 0. + endif + 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 + + + BEGIN_PROVIDER [ double precision, mo_value_at_fitcusp_radius, (mo_num_8,nucl_num) ] +&BEGIN_PROVIDER [ double precision, mo_grad_at_fitcusp_radius, (mo_num_8,nucl_num) ] +&BEGIN_PROVIDER [ double precision, mo_lapl_at_fitcusp_radius, (mo_num_8,nucl_num) ] + implicit none + BEGIN_DOC +! 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 + mo_grad_at_fitcusp_radius(j,k) = 0.d0 + mo_lapl_at_fitcusp_radius(j,k) = 0.d0 + !DIR$ VECTOR ALIGNED + do i=1,ao_num + mo_value_at_fitcusp_radius(j,k) = mo_value_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_value_at_fitcusp_radius(i,k) + mo_grad_at_fitcusp_radius(j,k) = mo_grad_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_grad_at_fitcusp_radius(i,k) + mo_lapl_at_fitcusp_radius(j,k) = mo_lapl_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_lapl_at_fitcusp_radius(i,k) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ] + implicit none + BEGIN_DOC +! Parameters of the splines + END_DOC + integer :: i,k, niter + character*(80) :: message + + nucl_fitcusp_param = 0.d0 + do k=1,nucl_num + double precision :: r, Z + Z = nucl_charge(k) + if (Z < 1.d-2) then + ! Avoid dummy atoms + cycle + endif + R = nucl_fitcusp_radius(k) + !DIR$ LOOP COUNT (500) + do i=1,mo_num + double precision :: lap_phi, grad_phi, phi, eta + lap_phi = mo_lapl_at_fitcusp_radius(i,k) + 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 + + + +subroutine sparse_full_mv(A,LDA, & + B1,LDB, & + B2, B3, B4, B5, indices, & + C1,LDC,C2,C3,C4,C5,an) + implicit none + BEGIN_DOC +! Performs a vectorized product between a dense matrix (the MO coefficients +! matrix) and 5 sparse vectors (the value, gradients and laplacian of the AOs). + END_DOC + integer, intent(in) :: an,LDA,LDB,LDC + integer, intent(in) :: indices(0:LDB) + real, intent(in) :: A(LDA,an) + real, intent(in) :: B1(LDB) + real, intent(in) :: B2(LDB) + real, intent(in) :: B3(LDB) + real, intent(in) :: B4(LDB) + real, intent(in) :: B5(LDB) + real, intent(out) :: C1(LDC) + real, intent(out) :: C2(LDC) + real, intent(out) :: C3(LDC) + real, intent(out) :: C4(LDC) + real, intent(out) :: C5(LDC) + !DIR$ ASSUME_ALIGNED A : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED B1 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED B2 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED B3 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED B4 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED B5 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED C1 : $IRP_ALIGN + !DIR$ ASSUME_ALIGNED C2 : $IRP_ALIGN + !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) + !DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: k_vec + real :: d11, d12, d13, d14, d15 + 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 + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (256) + !$OMP SIMD + do j=1,LDC + C1(j) = 0. + C2(j) = 0. + C3(j) = 0. + C4(j) = 0. + C5(j) = 0. + enddo + !$OMP END SIMD + kmax2 = ishft(indices(0),-2) + kmax2 = ishft(kmax2,2) + kmax3 = indices(0) + !DIR$ LOOP COUNT (200) + do kao=1,kmax2,4 + k_vec(1) = indices(kao ) + 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) + + !DIR$ LOOP COUNT (256) + do k=0,LDA-1,$IRP_ALIGN/4 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$IRP_ALIGN/4 + 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 + 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 + enddo + !$OMP END SIMD + enddo + + 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) + + !DIR$ LOOP COUNT (256) + do k=0,LDA-1,$IRP_ALIGN/4 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$IRP_ALIGN/4 + 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 + 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 + enddo + !$OMP END SIMD + enddo + + d15 = B5(kao ) + d25 = B5(kao+1) + d35 = B5(kao+2) + d45 = B5(kao+3) + + !DIR$ LOOP COUNT (256) + do k=0,LDA-1,$IRP_ALIGN/4 + !DIR$ VECTOR ALIGNED + !$OMP SIMD + do j=1,$IRP_ALIGN/4 + C5(j+k) = C5(j+k) + A(j+k,k_vec(1))*d15 + A(j+k,k_vec(2))*d25& + + A(j+k,k_vec(3))*d35 + A(j+k,k_vec(4))*d45 + enddo + !$OMP END SIMD + enddo + + enddo + + !DIR$ LOOP COUNT (200) + do kao = kmax2+1, kmax3 + k_vec(1) = indices(kao) + d11 = B1(kao) + d12 = B2(kao) + d13 = B3(kao) + d14 = B4(kao) + d15 = B5(kao) + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (256) + do k=0,LDA-1,simd_sp + !DIR$ VECTOR ALIGNED + !$OMP SIMD + 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 + C3(j+k) = C3(j+k) + A(j+k,k_vec(1))*d13 + C4(j+k) = C4(j+k) + A(j+k,k_vec(1))*d14 + C5(j+k) = C5(j+k) + A(j+k,k_vec(1))*d15 + enddo + !$OMP END SIMD + enddo + enddo + +end + + + diff --git a/src/mo_point.irp.f b/src/mo_point.irp.f new file mode 100644 index 0000000..5a50a86 --- /dev/null +++ b/src/mo_point.irp.f @@ -0,0 +1,21 @@ +BEGIN_PROVIDER [ real, mo_value_p, (mo_tot_num) ] + implicit none + + BEGIN_DOC + ! Values of the molecular orbitals + END_DOC + + integer :: i, j, k + + do j=1,mo_num + mo_value_p(j) = 0. + enddo + do k=1,ao_num + do j=1,mo_num + mo_value_p(j) = mo_value_p(j)+mo_coef_transp(j,k)*ao_value_p(k) + enddo + enddo + +END_PROVIDER + + diff --git a/src/nuclei.irp.f b/src/nuclei.irp.f new file mode 100644 index 0000000..39d84a5 --- /dev/null +++ b/src/nuclei.irp.f @@ -0,0 +1,148 @@ + BEGIN_PROVIDER [ integer, nucl_num ] +&BEGIN_PROVIDER [ integer, nucl_num_8 ] + implicit none + BEGIN_DOC +! Number of nuclei + END_DOC + + nucl_num = -1 + call get_nuclei_nucl_num(nucl_num) + if (nucl_num <= 0) then + call abrt(irp_here,'Number of nuclei should be > 0') + endif + integer, external :: mod_align + nucl_num_8 = mod_align(nucl_num) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, nucl_charge, (nucl_num) ] + implicit none + BEGIN_DOC +! Nuclear charge + END_DOC + + nucl_charge = -1.d0 + call get_nuclei_nucl_charge(nucl_charge) + + integer :: i + do i=1,nucl_num + if (nucl_charge(i) < 0.) then + call abrt(irp_here,'Nuclear charges should be > 0') + endif + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ real, nucl_coord, (nucl_num_8,3) ] + implicit none + BEGIN_DOC +! Nuclear coordinates + END_DOC + + nucl_coord = 0. + real, allocatable :: buffer(:,:) + allocate (buffer(nucl_num,3)) + buffer = 0. + call get_nuclei_nucl_coord(buffer) + integer :: i,j + + do i=1,3 + do j=1,nucl_num + nucl_coord(j,i) = buffer(j,i) + enddo + enddo + deallocate(buffer) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, nucl_coord_transp, (8,nucl_num) + implicit none + BEGIN_DOC +! Transposed array of nucl_coord + END_DOC + integer :: i, k + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + nucl_coord_transp = 0. + endif + + !DIR$ VECTOR ALIGNED + do i=1,nucl_num + nucl_coord_transp(1,i) = nucl_coord(i,1) + nucl_coord_transp(2,i) = nucl_coord(i,2) + nucl_coord_transp(3,i) = nucl_coord(i,3) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ real, nucl_dist, (nucl_num_8,nucl_num) ] +&BEGIN_PROVIDER [ real, nucl_dist_vec_x, (nucl_num_8,nucl_num) ] +&BEGIN_PROVIDER [ real, nucl_dist_vec_y, (nucl_num_8,nucl_num) ] +&BEGIN_PROVIDER [ real, nucl_dist_vec_z, (nucl_num_8,nucl_num) ] + implicit none + BEGIN_DOC +! nucl_dist : Nucleus-nucleus distances : nucl_dist(i,j) = |R_i-R_j| + +! nucl_dist_vec : Nucleus-nucleus distances vectors + END_DOC + + integer :: ie1, ie2, l + integer,save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + nucl_dist = 0. + nucl_dist_vec_x = 0. + nucl_dist_vec_y = 0. + nucl_dist_vec_z = 0. + endif + + do ie2 = 1,nucl_num + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do ie1 = 1,nucl_num + nucl_dist_vec_x(ie1,ie2) = nucl_coord(ie1,1) - nucl_coord(ie2,1) + nucl_dist_vec_y(ie1,ie2) = nucl_coord(ie1,2) - nucl_coord(ie2,2) + nucl_dist_vec_z(ie1,ie2) = nucl_coord(ie1,3) - nucl_coord(ie2,3) + enddo + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do ie1 = 1,nucl_num + nucl_dist (ie1,ie2) = nucl_dist_vec_x(ie1,ie2)*nucl_dist_vec_x(ie1,ie2) +& + nucl_dist_vec_y(ie1,ie2)*nucl_dist_vec_y(ie1,ie2) + & + nucl_dist_vec_z(ie1,ie2)*nucl_dist_vec_z(ie1,ie2) + nucl_dist(ie1,ie2) = sqrt(nucl_dist (ie1,ie2)) + ASSERT (nucl_dist(ie1,ie2) > 0.) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, nucl_fitcusp_radius, (nucl_num) ] + implicit none + BEGIN_DOC +! Distance threshold for the fit + END_DOC + real :: def(nucl_num) + integer :: k + + if (.not.do_nucl_fitcusp) then + nucl_fitcusp_radius = 0.d0 + return + endif + do k=1,nucl_num + nucl_fitcusp_radius(k) = .5/nucl_charge(k) + enddo + call get_nuclei_nucl_fitcusp_radius(nucl_fitcusp_radius) + ! Avoid dummy atoms + do k=1,nucl_num + if (nucl_charge(k) < 5.d-1) then + nucl_fitcusp_radius(k) = 0. + endif + enddo + +END_PROVIDER + + diff --git a/src/point.irp.f b/src/point.irp.f new file mode 100644 index 0000000..5f46883 --- /dev/null +++ b/src/point.irp.f @@ -0,0 +1,39 @@ +BEGIN_PROVIDER [ real, point, (3) ] + implicit none + BEGIN_DOC + ! Coordinates of the current point + END_DOC + point(1) = 0. + point(2) = 0. + point(3) = 0. +END_PROVIDER + + +BEGIN_PROVIDER [ real, point_nucl_dist_vec, (nucl_num,3) ] + implicit none + BEGIN_DOC + ! Distance vector between the current point and the nuclei + END_DOC + integer :: k + do k=1,nucl_num + point_nucl_dist_vec(k,1) = point(1)-nucl_coord(k,1) + point_nucl_dist_vec(k,2) = point(2)-nucl_coord(k,2) + point_nucl_dist_vec(k,3) = point(3)-nucl_coord(k,3) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ real, point_nucl_dist, (nucl_num) ] + implicit none + BEGIN_DOC + ! Distance between the current point and the nuclei + END_DOC + integer :: k,l + do k=1,nucl_num + point_nucl_dist(k) = point_nucl_dist_vec(k,1)*point_nucl_dist_vec(k,1) +& + point_nucl_dist_vec(k,2)*point_nucl_dist_vec(k,2) + & + point_nucl_dist_vec(k,3)*point_nucl_dist_vec(k,3) + point_nucl_dist(k) = sqrt(point_nucl_dist(k)) + enddo +END_PROVIDER + diff --git a/src/prepare_walkers.irp.f b/src/prepare_walkers.irp.f new file mode 100644 index 0000000..d6f6622 --- /dev/null +++ b/src/prepare_walkers.irp.f @@ -0,0 +1,217 @@ +subroutine draw_init_points + implicit none + BEGIN_DOC +! Place randomly electron around nuclei + END_DOC + integer :: iwalk + logical, allocatable :: do_elec(:) + integer :: acc_num + + real, allocatable :: xmin(:,:) + + integer :: i, j, k, l, kk + + real :: norm + allocate (do_elec(elec_num), xmin(3,elec_num)) + xmin = -huge(1.) + norm = 0. + do i=1,elec_alpha_num + do j=1,ao_num + norm += mo_coef_transp(i,j)*mo_coef_transp(i,j) + enddo + enddo + norm = sqrt(norm/float(elec_alpha_num)) + call rinfo( irp_here, 'Norm : ', norm ) + call rinfo( irp_here, 'mo_scale: ' , mo_scale ) + mo_coef_transp = mo_coef_transp/norm + + double precision :: qmc_ranf + real :: mo_max + do i=1,elec_alpha_num + l=1 + xmin(1,i) = mo_coef_transp(i,1)*mo_coef_transp(i,1) - 0.001*qmc_ranf() + do j=2,ao_num + xmin(2,i) = mo_coef_transp(i,j)*mo_coef_transp(i,j) - 0.001*qmc_ranf() + if (xmin(2,i) > xmin(1,i) ) then + xmin(1,i) = xmin(2,i) + l = ao_nucl(j) + endif + enddo + xmin(1,i) = nucl_coord(l,1) + xmin(2,i) = nucl_coord(l,2) + xmin(3,i) = nucl_coord(l,3) + enddo + + call iinfo(irp_here, 'Det num = ', det_num ) + do k=1,elec_beta_num + i = k+elec_alpha_num + l=1 + xmin(1,i) = mo_coef_transp(k,1)*mo_coef_transp(k,1) - 0.001*qmc_ranf() + do j=2,ao_num + xmin(2,i) = mo_coef_transp(k,j)*mo_coef_transp(k,j) - 0.001*qmc_ranf() + if (xmin(2,i) > xmin(1,i) ) then + xmin(1,i) = xmin(2,i) + l = ao_nucl(j) + endif + enddo + xmin(1,i) = nucl_coord(l,1) + xmin(2,i) = nucl_coord(l,2) + xmin(3,i) = nucl_coord(l,3) + enddo + + call rinfo( irp_here, 'time step =', time_step ) + do iwalk=1,walk_num + call iinfo( irp_here, 'Generate initial positions for walker', iwalk ) + acc_num = 0 + do_elec = .True. + do while (acc_num < elec_num) + double precision :: gauss + real :: re_compute + re_compute = 0. + do while (re_compute < 1.e-6) + do i=1,elec_num + if (do_elec(i)) then + do l=1,3 + elec_coord(i,l) = xmin(l,i) + 2.*(0.5-qmc_ranf()) + enddo + endif + enddo + TOUCH elec_coord + re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num)) + enddo + + do i=1,elec_alpha_num + if (do_elec(i)) then + if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then + acc_num += 1 + do_elec(i) = .False. + endif + endif + enddo + + do i=1,elec_beta_num + if (do_elec(i+elec_alpha_num)) then + if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then + acc_num += 1 + do_elec(i+elec_alpha_num) = .False. + endif + endif + enddo + + + enddo + + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,iwalk) = elec_coord(i,l) + enddo + enddo + enddo + if (.not.is_worker) then + call ezfio_set_electrons_elec_coord_pool_size(walk_num) + call ezfio_set_electrons_elec_coord_pool(elec_coord_full) + endif + SOFT_TOUCH elec_coord elec_coord_full + deallocate (do_elec, xmin) + +end + + +subroutine run_prepare_walkers + implicit none + BEGIN_DOC +! Create starting points for walkers + END_DOC + include 'types.F' + integer :: istep, iwalk + integer :: i,j, l + + do iwalk=1,walk_num + do l=1,3 + do i=1,elec_num+1 + elec_coord(i,l) = elec_coord_full(i,l,iwalk) + enddo + enddo + TOUCH elec_coord + + double precision :: qmc_ranf, rcond, lambda + rcond = 100.d0 + lambda = 1.d0 + do while ( (rcond > 3.d0) .or. (rcond < -3.d0) ) + rcond = 0. + do i=1,elec_alpha_num + rcond += log(lambda*abs(mo_value_transp(i,i))) + enddo + do i=1,elec_beta_num + rcond += log(lambda*abs(mo_value_transp(i,elec_alpha_num+i))) + enddo + if (rcond > 2.d0) then + lambda = lambda/(1.d0+.1*qmc_ranf()) + endif + if (rcond< -2.d0) then + lambda = lambda*(1.d0+.1*qmc_ranf()) + endif + enddo + do i=1,ao_num + !DIR$ VECTOR ALIGNED + do j=1,mo_num_8 + mo_coef_transp(j,i) *= lambda + enddo + enddo + TOUCH mo_coef_transp + call iinfo (irp_here, 'Starting walker ', iwalk ) + + do istep=1,1000 + if (single_det_value == 0.d0) then + exit + endif + prepare_walkers_t = float(istep)/1000. + TOUCH prepare_walkers_t + rcond = log(abs(dble(single_det_value))) + real :: factor + rcond = log(abs(dble(single_det_value))) + integer :: icount + icount = 0 + do while ( (rcond > 10.d0) .or. (rcond < -10.d0) ) + icount += 1 + if (icount == 1000) then + exit + endif + if (rcond > 10.d0) then + factor = 1./(1.+.10) !*qmc_ranf()) + else if (rcond< -10.d0) then + factor = 1.+.10 !*qmc_ranf() + endif + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,mo_num_8 + mo_coef_transp(i,j) *= factor + enddo + enddo + TOUCH mo_coef_transp + + rcond = log(abs(dble(single_det_value))) + enddo + double precision :: p,q + logical :: accepted + real :: delta_x + accepted = .False. + do while (.not.accepted) + if (vmc_algo == t_Brownian) then + call brownian_step(p,q,accepted,delta_x) + else if (vmc_algo == t_Langevin) then + call langevin_step(p,q,accepted,delta_x) + endif + enddo + enddo + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,iwalk) = elec_coord(i,l) + enddo + enddo + call iinfo(irp_here, 'Walker done', iwalk) + TOUCH elec_coord_full + enddo + +end + diff --git a/src/properties.py b/src/properties.py new file mode 100755 index 0000000..ff180f2 --- /dev/null +++ b/src/properties.py @@ -0,0 +1,95 @@ +#!/usr/bin/python + +import string +import os + +properties = [] +dims = {} + +files = filter(lambda x: x.startswith("PROPERTIES") and \ + x.endswith("irp.f"), os.listdir(os.getcwd())) + +files = map(lambda x: 'PROPERTIES/'+x, filter(lambda x: x.endswith("irp.f"), os.listdir(os.getcwd()+'/PROPERTIES'))) + + +#files = filter(lambda x: x.endswith("irp.f"), os.listdir(os.getcwd())) + +for filename in files: + lines = [] + check_dims = False + file = open(filename,'r') + lines += file.readlines() + file.close() + for i,line in enumerate(lines): + if line.startswith("! PROPERTIES"): + lines = lines[i:] + break + + for line in map(lambda x: x.lower(),lines): + if line.lstrip().startswith('begin_provider'): + check_dims = False + buffer = line + buffer = buffer.split('[')[1] + buffer = buffer.split(']')[0] + buffer = buffer.split(',') + if (len(buffer) == 2): + buffer.append("") + else: + buffer = [ buffer[0], buffer[1], ','.join(buffer[2:]) ] + check_dims = True + buffer = map(string.strip,buffer) + properties.append(buffer) + current_prop = buffer[1] + elif check_dims: + if 'dimensions :' in line: + dims[current_prop] = line.split(':')[1].strip() + + + +def sq(item): + return [item[0], item[1]+"_2", item[2]] +properties_with_square = properties + map(sq,properties) + +for p in [ properties, properties_with_square ]: + def c(x,y): + if x[1] > y[1]: return 1 + if x[1] == y[1]: return 0 + if x[1] < y[1]: return -1 + p.sort(c) + +def namelist(): + buffer = "" + for p in properties: + buffer += "calc_"+p[1]+", &\n" + buffer = buffer[:-4] + result = " namelist /properties/"+buffer + return result + +def touch_all(): + print "TOUCH", + for p in properties: + print "calc_"+p[1], + print "" + +#file = open('../scripts/properties.py','w') +#print >>file,'properties = ',properties +#file.close() + +def make_dims(): + template = """ +BEGIN_PROVIDER [ integer, size_%(p)s ] + implicit none + BEGIN_DOC +! Size of %(p)s + END_DOC + if (calc_%(p)s) then + size_%(p)s = %(d)s + else + size_%(p)s = 1 + endif +END_PROVIDER +""" + for p in dims: + print template%{'p': p, 'd': dims[p]} + + diff --git a/src/pseudo.irp.f b/src/pseudo.irp.f new file mode 100644 index 0000000..b0999e1 --- /dev/null +++ b/src/pseudo.irp.f @@ -0,0 +1,383 @@ +BEGIN_PROVIDER [ double precision, pseudo_v_k , (nucl_num,pseudo_klocmax) ] + implicit none + BEGIN_DOC +! V_k + END_DOC + call get_pseudo_pseudo_v_k(pseudo_v_k) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pseudo_v_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ] + implicit none + BEGIN_DOC +! V_Kl + END_DOC + call get_pseudo_pseudo_v_kl(pseudo_v_kl) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_kmax ] + implicit none + BEGIN_DOC +! kmax + END_DOC + call get_pseudo_pseudo_kmax(pseudo_kmax) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pseudo_dz_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ] + implicit none + BEGIN_DOC +! Exponents in the non-local part of the pseudo-potential + END_DOC + call get_pseudo_pseudo_dz_kl(pseudo_dz_kl) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_klocmax ] + implicit none + BEGIN_DOC +! klocmax + END_DOC + call get_pseudo_pseudo_klocmax(pseudo_klocmax) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_lmax ] + implicit none + BEGIN_DOC +! Max value of l in the pseudo-potential + END_DOC + call get_pseudo_pseudo_lmax(pseudo_lmax) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_grid_size ] + implicit none + BEGIN_DOC +! Size of the QMC grid (number of points) + END_DOC + call get_pseudo_pseudo_grid_size(pseudo_grid_size) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pseudo_grid_rmax ] + implicit none + BEGIN_DOC +! Size of the QMC grid (max distance) + END_DOC + call get_pseudo_pseudo_grid_rmax(pseudo_grid_rmax) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, pseudo_non_loc_dim ] +&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_8 ] +&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_count, (nucl_num) ] + implicit none + BEGIN_DOC + ! Dimension of of pseudo_non_local arrays + END_DOC + pseudo_non_loc_dim = 0 + pseudo_non_loc_dim_count = 0 + integer :: k,l,m + do k=1,nucl_num + do l=0,pseudo_lmax + pseudo_non_loc_dim_count(k) += 2*l+1 + enddo + enddo + pseudo_non_loc_dim = sum(pseudo_non_loc_dim_count) + integer, external :: mod_align + pseudo_non_loc_dim_8 = mod_align(pseudo_non_loc_dim) +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_n_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ] + implicit none + BEGIN_DOC +! n_kl + END_DOC + call get_pseudo_pseudo_n_kl(pseudo_n_kl) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pseudo_dz_k , (nucl_num,pseudo_klocmax) ] + implicit none + BEGIN_DOC +! dz_k + END_DOC + call get_pseudo_pseudo_dz_k(pseudo_dz_k) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, pseudo_n_k , (nucl_num,pseudo_klocmax) ] + implicit none + BEGIN_DOC +! n_k + END_DOC + call get_pseudo_pseudo_n_k(pseudo_n_k) + +END_PROVIDER + +BEGIN_PROVIDER [ logical, do_pseudo ] + implicit none + BEGIN_DOC +! Using pseudo potential integral of not + END_DOC + call get_pseudo_do_pseudo(do_pseudo) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ] + implicit none + BEGIN_DOC + ! Pseudopotential grid points + END_DOC + call get_pseudo_ao_pseudo_grid(ao_pseudo_grid) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ] + implicit none + BEGIN_DOC + ! Pseudopotential grid points + END_DOC + call get_pseudo_mo_pseudo_grid(mo_pseudo_grid) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, mo_pseudo_grid_scaled, (pseudo_non_loc_dim_8,ao_num,pseudo_grid_size) ] + implicit none + BEGIN_DOC + ! Pseudopotential grid points + END_DOC + integer :: i,k,l,m,kk,n + double precision :: c + c = 1.d0/mo_scale + do n=1,pseudo_grid_size + do i=1,ao_num + kk = 0 + do k=1,nucl_num + do l=0,pseudo_lmax + do m=-l,l + kk = kk+1 + mo_pseudo_grid_scaled(kk,i,n) = c * mo_pseudo_grid(i,m,l,k,n) + enddo + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, v_pseudo_local, (elec_num) ] + implicit none + BEGIN_DOC + ! Local component of the pseudo-potential + END_DOC + integer :: i,k,l + double precision :: r, rn, alpha + do l=1,elec_num + v_pseudo_local(l) = 0.d0 + do k=1,pseudo_klocmax + do i=1,nucl_num + r = nucl_elec_dist(i,l) + alpha = pseudo_dz_k(i,k)*r*r + if (alpha > 20.d0) then + cycle + endif + select case (pseudo_n_k(i,k)) + case (-2) + rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l) + case (-1) + rn = nucl_elec_dist_inv(i,l) + case (0) + rn = 1.d0 + case (1) + rn = r + case (2) + rn = r*r + case default + rn = r**pseudo_n_k(i,k) + end select + v_pseudo_local(l) += pseudo_v_k(i,k) * exp(-alpha) * rn + enddo + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, v_pseudo_non_local, (pseudo_non_loc_dim_8,elec_num) ] + implicit none + BEGIN_DOC + ! Non-Local component of the pseudo-potential + END_DOC + integer :: i,k,j,l + integer :: kk,m + double precision :: r, rn, r2, alpha + double precision :: tmp(0:pseudo_lmax) + v_pseudo_non_local = 0.d0 + do j=1,elec_num + + kk = 0 + do i=1,nucl_num + r = nucl_elec_dist(i,j) + r2 = r*r + tmp(0:pseudo_lmax) = 0.d0 + do k=1,pseudo_kmax + do l=0,pseudo_lmax + alpha = pseudo_dz_kl(i,k,l)*r2 + if (alpha > 20.d0) then + cycle + endif + select case (pseudo_n_kl(i,k,l)) + case (0) + rn = 1.d0 + case (-1) + rn = nucl_elec_dist_inv(i,j) + case (1) + rn = r + case (2) + rn = r*r + case (-2) + rn = nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j) + case default + rn = r**pseudo_n_kl(i,k,l) + end select + tmp(l) += pseudo_v_kl(i,k,l) * exp(-alpha) * rn + enddo + enddo + do l=0,pseudo_lmax + do m=-l,l + kk += 1 + v_pseudo_non_local(kk,j) = tmp(l) + enddo + enddo + enddo + + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ] + implicit none + BEGIN_DOC +! \sum < Ylm | MO > x Ylm(r) x V_nl(r) + END_DOC + integer :: ii,i,j,l,m,k,n,kk + double precision :: r, dr_inv + double precision :: tmp(pseudo_non_loc_dim_8,mo_num) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: tmp + integer, save :: ifirst = 0 + + if (ifirst == 0) then + pseudo_mo_term = 0.d0 + ifirst = 1 + endif + + PROVIDE pseudo_ylm present_mos mo_pseudo_grid_scaled pseudo_non_loc_dim_count + dr_inv = dble(pseudo_grid_size)/pseudo_grid_rmax + do j=1,elec_num + tmp = 0.d0 + kk=0 + do k=1,nucl_num + r = nucl_elec_dist(k,j) + n = 1 + int(0.5+r*dr_inv) + if (n<=pseudo_grid_size) then + do ii=1,num_present_mos + i = present_mos(ii) + !DIR$ LOOP COUNT(4) + do l=kk+1,kk+pseudo_non_loc_dim_count(k) + tmp(l,i) = mo_pseudo_grid_scaled (l,i,n) * pseudo_ylm(l,j) + enddo + enddo + endif + kk = kk+pseudo_non_loc_dim_count(k) + enddo + do ii=1,num_present_mos + i = present_mos(ii) + pseudo_mo_term(i,j) = 0.d0 + do kk=1,pseudo_non_loc_dim + pseudo_mo_term(i,j) = pseudo_mo_term(i,j) + tmp(kk,i) * v_pseudo_non_local(kk,j) + enddo + enddo + enddo +END_PROVIDER + + +real function ylm(l,m,x,y,z,r_inv) + implicit none + include 'constants.F' + BEGIN_DOC +! Y(l,m) evaluated at x,y,z + END_DOC + integer, intent(in) :: l,m + real, intent(in) :: x,y,z,r_inv + + ylm = 0.5 * one_over_sqpi + select case(l) + + case(0) + continue + + case(1) + select case (m) + case (-1) + ylm = ylm * sq3 * y * r_inv + case (0) + ylm = ylm * sq3 * z * r_inv + case (1) + ylm = ylm * sq3 * x * r_inv + end select + +! case(2) +! select case (m) +! case(-2) +! ylm = ylm * sqrt(15.) * x * y * r_inv * r_inv +! case(-1) +! ylm = ylm * sqrt(15.) * y * z * r_inv * r_inv +! case(0) +! ylm = 0.5 * ylm * sqrt(15.) * (2.*z*z - x*x - y*y) * r_inv * r_inv +! case(1) +! ylm = ylm * sqrt(15.) * z * x * r_inv * r_inv +! case(2) +! ylm = 0.5 * ylm * sqrt(15.) * (x*x - y*y) * r_inv * r_inv +! end select + + case default + stop 'problem in Ylm of pseudo' + + end select + +end + + + +BEGIN_PROVIDER [ double precision, pseudo_ylm, (pseudo_non_loc_dim_8,elec_num) ] + implicit none + BEGIN_DOC + ! Y(l,m) evaluated for every electron position centered on every nuclei + END_DOC + + integer :: i,j,l,m,kk + real, external :: ylm + integer, save :: ifirst = 0 + + if (ifirst == 0) then + pseudo_ylm = 0.d0 + ifirst = 1 + endif + + do j=1,elec_num + kk = 0 + do i=1,nucl_num + do l=0,pseudo_lmax + do m=-l,l + kk = kk+1 + pseudo_ylm(kk,j) = ylm(l,m, & + nucl_elec_dist_vec(1,i,j), & + nucl_elec_dist_vec(2,i,j), & + nucl_elec_dist_vec(3,i,j), & + nucl_elec_dist_inv(i,j)) + enddo + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/psi.irp.f b/src/psi.irp.f new file mode 100644 index 0000000..791f89d --- /dev/null +++ b/src/psi.irp.f @@ -0,0 +1,102 @@ + BEGIN_PROVIDER [ double precision, psi_value ] + implicit none + BEGIN_DOC + ! Value of the wave function + END_DOC + + psi_value = psidet_value*jast_value + if (psi_value == 0.d0) then + call abrt(irp_here,"Value of the wave function is 0.") + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_value_inv ] + implicit none + BEGIN_DOC +! 1./psi_value + END_DOC + psi_value_inv = 1.d0/psi_value +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_value_inv2 ] + implicit none + BEGIN_DOC + ! 1./(psi_value)**2 + END_DOC + psi_value_inv2 = psi_value_inv*psi_value_inv +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, psi_grad_x, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_y, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_z, (elec_num_8) ] + implicit none + BEGIN_DOC +! Gradients of the wave function + END_DOC + + integer :: j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_grad_x(j) = psi_grad_psi_inv_x(j)*psi_value + psi_grad_y(j) = psi_grad_psi_inv_y(j)*psi_value + psi_grad_z(j) = psi_grad_psi_inv_z(j)*psi_value + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, psi_lapl, (elec_num_8) ] + implicit none + BEGIN_DOC + ! Laplacian of the wave function + END_DOC + + integer :: i, j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_lapl(j) = jast_value*(psidet_grad_lapl(4,j) + psidet_value*jast_lapl_jast_inv(j) + 2.d0*(& + psidet_grad_lapl(1,j)*jast_grad_jast_inv_x(j) + & + psidet_grad_lapl(2,j)*jast_grad_jast_inv_y(j) + & + psidet_grad_lapl(3,j)*jast_grad_jast_inv_z(j) )) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_y, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_z, (elec_num_8) ] + implicit none + BEGIN_DOC +! grad(psi)/psi + END_DOC + + integer :: j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_grad_psi_inv_x(j) = psidet_grad_lapl(1,j)*psidet_inv + jast_grad_jast_inv_x(j) + psi_grad_psi_inv_y(j) = psidet_grad_lapl(2,j)*psidet_inv + jast_grad_jast_inv_y(j) + psi_grad_psi_inv_z(j) = psidet_grad_lapl(3,j)*psidet_inv + jast_grad_jast_inv_z(j) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, psi_lapl_psi_inv, (elec_num_8) ] + implicit none + BEGIN_DOC + ! (Laplacian psi) / psi + END_DOC + + integer :: i, j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_lapl_psi_inv(j) = psi_lapl(j)*psi_value_inv + enddo + +END_PROVIDER + diff --git a/src/simulation.irp.f b/src/simulation.irp.f new file mode 100644 index 0000000..3283dda --- /dev/null +++ b/src/simulation.irp.f @@ -0,0 +1,350 @@ +BEGIN_PROVIDER [ logical, is_worker ] + implicit none + BEGIN_DOC + ! True if the process is a worker process + END_DOC + is_worker = .False. +END_PROVIDER + + +BEGIN_PROVIDER [ integer, walk_num_tot ] + implicit none + BEGIN_DOC + ! Total number of walkers + END_DOC + + walk_num_tot = 1000 + call get_electrons_elec_walk_num_tot(walk_num_tot) + walk_num_tot = max(walk_num,walk_num_tot) + call iinfo(irp_here,'walk_num', walk_num_tot) + + if (walk_num_tot <= 0) then + call abrt(irp_here,'Total number of walkers should be > 0') + endif +END_PROVIDER + + + BEGIN_PROVIDER [ integer, walk_num ] +&BEGIN_PROVIDER [ integer, walk_num_8 ] + implicit none + BEGIN_DOC + ! Number of walkers + END_DOC + + walk_num = 100 + call get_electrons_elec_walk_num(walk_num) + call iinfo(irp_here,'walk_num', walk_num) + + if (walk_num <= 0) then + call abrt(irp_here,'Number of walkers should be > 0') + endif + integer :: mod_align + walk_num_8 = mod_align(walk_num) +END_PROVIDER + + +BEGIN_PROVIDER [ logical, do_equilibration ] + implicit none + BEGIN_DOC + ! Equilibrate walkers + END_DOC + + do_equilibration = .True. + if (.not.do_prepare) then + call get_simulation_equilibration(do_equilibration) + endif + call iinfo(irp_here,'equilibration', do_equilibration) + +END_PROVIDER + + +BEGIN_PROVIDER [ logical, do_prepare ] + implicit none + BEGIN_DOC + ! If true, prepare new walkers + END_DOC + do_prepare = .False. + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, block_time ] + implicit none + BEGIN_DOC + ! Wall time requested to realize one block + END_DOC + block_time = 30.d0 + integer :: block_time_int + call get_simulation_block_time(block_time_int) + + if (block_time<= 1) then + call abrt(irp_here,'Block time should be > 1s') + endif + double precision, external :: qmc_ranf + block_time = dble(block_time_int) + qmc_ranf() + call dinfo(irp_here,'block_time',block_time) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, stop_time ] + implicit none + BEGIN_DOC + ! Termination condition of the run + END_DOC + stop_time = 3600*24 + call get_simulation_stop_time(stop_time) + call iinfo(irp_here,'stop_time',stop_time) + + if (stop_time<= 1) then + call abrt(irp_here,'Stop time should be > 1s') + endif +END_PROVIDER + + + BEGIN_PROVIDER [ real, time_step ] +&BEGIN_PROVIDER [ real, time_step_inv ] +&BEGIN_PROVIDER [ double precision, dtime_step ] + implicit none + BEGIN_DOC + ! time_step : The time step of the random walk + END_DOC + + time_step = 0.0 + call get_simulation_time_step(time_step) + call rinfo(irp_here,'time_step',time_step) + + if (time_step <= 0.) then + call abrt(irp_here,'Time step should be > 0') + endif + dtime_step = dble(time_step) + time_step_inv = 1./time_step +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, time_step_sq ] +&BEGIN_PROVIDER [ double precision, time_step_exp ] +&BEGIN_PROVIDER [ double precision, time_step_exp_sq ] +&BEGIN_PROVIDER [ double precision, time_step_exp_sq_sq ] + implicit none + BEGIN_DOC + ! + ! time_step_sq : sqrt(time_step) + ! + ! time_step_exp = exp(-time_step) + ! + ! time_step_exp_sq = sqrt(time_step_exp) + ! + ! time_step_exp_sq_sq = sqrt(time_step_exp_sq) + END_DOC + + time_step_sq = sqrt(dble(time_step)) + time_step_exp = exp(-dble(time_step)) + time_step_exp_sq = sqrt(time_step_exp) + time_step_exp_sq_sq = sqrt(time_step_exp_sq) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, qmc_method ] + implicit none + include 'types.F' + BEGIN_DOC + ! qmc_method : Calculation method. Can be t_VMC, t_DMC + END_DOC + character*(32) :: method + method = types(t_VMC) + call get_simulation_method(method) + + if (method == types(t_VMC)) then + qmc_method = t_VMC + else if (method == types(t_DMC)) then + qmc_method = t_DMC + else + call abrt(irp_here, 'Method should be ( VMC | DMC )') + endif + + call cinfo(irp_here,'qmc_method',trim(method)) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, events_num ] + BEGIN_DOC + ! Number of Monte Carlo events to average + END_DOC + events_num = real(walk_num)*real(step_num) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, walk_i ] + BEGIN_DOC + ! Current walker + END_DOC + walk_i = 1 +END_PROVIDER + + + BEGIN_PROVIDER [ real, accepted_num ] +&BEGIN_PROVIDER [ real, rejected_num ] + BEGIN_DOC + ! Number of accepted steps + ! Number of rejected steps + END_DOC + accepted_num= 0. + rejected_num= 0. +END_PROVIDER + + +BEGIN_PROVIDER [ logical, save_data ] + implicit none + BEGIN_DOC + ! If true, the updated simulation data is saved for restart. + END_DOC + save_data = .True. + call get_simulation_save_data(save_data) + call linfo(irp_here,'save_data',save_data) + +END_PROVIDER + + +real function accep_rate() + if ((accepted_num+rejected_num) > 0.) then + accep_rate = accepted_num/(accepted_num+rejected_num) + else + accep_rate = 0. + endif +end + + +logical function first_step() + first_step = (accepted_num+rejected_num < 1.) +end + + +subroutine accep_reset + FREE accepted_num + FREE rejected_num +end + + + + +BEGIN_PROVIDER [ integer, print_level ] + + BEGIN_DOC + ! Level of verbosity for standard output printing + END_DOC + print_level = 1 + call get_simulation_print_level(print_level) + +END_PROVIDER + + +BEGIN_PROVIDER [ character*(64), hostname] + implicit none + BEGIN_DOC + ! Name of the current host + END_DOC + call HOSTNM(hostname) +END_PROVIDER + + +BEGIN_PROVIDER [ logical, do_nucl_fitcusp ] + implicit none + BEGIN_DOC + ! If true, do the fit of the electron-nucleus cusp + END_DOC + do_nucl_fitcusp = .True. + call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp) + call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, vmc_algo ] + implicit none + include 'types.F' + BEGIN_DOC + ! Type of VMC algorithm: Brownian, MTM or Langevin + END_DOC + character*(32) :: Sampling + + Sampling = types(t_Langevin) + call get_simulation_sampling(Sampling) + + vmc_algo = 0 + if (Sampling == types(t_Brownian)) then + vmc_algo = t_Brownian + else if (Sampling == types(t_Langevin)) then + vmc_algo = t_Langevin + if (qmc_method == t_DMC) then + vmc_algo = t_Brownian + endif + else if (Sampling == types(t_MTM)) then + vmc_algo = t_MTM + else + call abrt(irp_here,'Sampling should be (Brownian|Langevin|MTM|Read)') + endif + call cinfo(irp_here,'Sampling',Sampling) + + ASSERT (vmc_algo > 0) +END_PROVIDER + + +BEGIN_PROVIDER [ character*(512), ezfio_filename ] + implicit none + BEGIN_DOC + ! Name of the ezfio file. + ! Defined in init_ezfio_filename + END_DOC + integer :: command_argument_count + if (command_argument_count() == 0) then + ezfio_filename = 'NOT_SET' + call ezfio_set_file(ezfio_filename) + else + call get_command_argument(1,ezfio_filename) + if (.not.is_worker) then + call ezfio_set_file(ezfio_filename) + endif + endif +END_PROVIDER + +BEGIN_PROVIDER [ character*(128), http_server ] + implicit none + BEGIN_DOC + ! Address of the data server + END_DOC + integer :: command_argument_count + if (command_argument_count() > 1) then + call get_command_argument(2,http_server) + else + call get_simulation_http_server(http_server) + endif +END_PROVIDER + + +subroutine read_do_run(do_run) + implicit none + integer, intent(out) :: do_run + BEGIN_DOC + ! Read the do_run variable from the ezfio directory + END_DOC + include 'types.F' + do_run = t_Stopping + call get_simulation_do_run(do_run) +end + + + +BEGIN_PROVIDER [ character*(32), md5_key ] + implicit none + BEGIN_DOC + ! Digest of the input + END_DOC + + md5_key = '' + call get_simulation_md5_key(md5_key) + + if (md5_key == '') then + call abrt(irp_here,'MD5 key of input is absent') + endif +END_PROVIDER + diff --git a/src/types.F b/src/types.F new file mode 100644 index 0000000..cc4ba67 --- /dev/null +++ b/src/types.F @@ -0,0 +1,34 @@ + integer, parameter :: t_Brownian = 3 + integer, parameter :: t_Langevin = 4 + integer, parameter :: t_MTM = 5 + integer, parameter :: t_Read = 6 + + integer, parameter :: t_VMC = 7 + integer, parameter :: t_DMC = 8 + + integer, parameter :: t_Simple = 11 + integer, parameter :: t_None = 12 + integer, parameter :: t_Core = 14 + + integer, parameter :: t_Stopped = 0 + integer, parameter :: t_Queued = 1 + integer, parameter :: t_Running = 2 + integer, parameter :: t_Stopping = 3 + + character*(32) :: types(15) = & + (/ '', & + '', & + 'Brownian', & + 'Langevin', & + '', & + '', & + 'VMC ', & + 'DMC ', & + ' ', & + '', & + 'Simple ', & + 'None ', & + ' ', & + 'Core ', & + ' '/) + diff --git a/src/wf.irp.f b/src/wf.irp.f new file mode 100644 index 0000000..b973e9a --- /dev/null +++ b/src/wf.irp.f @@ -0,0 +1,370 @@ +BEGIN_PROVIDER [ integer, i_state ] + implicit none + BEGIN_DOC + ! Current state + END_DOC + i_state = 1 +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_int ] + implicit none + BEGIN_DOC + ! Number of 64-bit integers needed to represent determinants as binary strings + END_DOC + call get_spindeterminants_n_int(N_int) +END_PROVIDER + +BEGIN_PROVIDER [ integer, bit_kind ] + implicit none + BEGIN_DOC + ! Number of octets per integer storing determinants + END_DOC + call get_spindeterminants_bit_kind(bit_kind) + ASSERT (bit_kind == 8) +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_states ] + implicit none + BEGIN_DOC + ! Number of states in EZFIO file + END_DOC + call get_spindeterminants_n_states(N_states) +END_PROVIDER + +BEGIN_PROVIDER [ integer, det_num_input ] + implicit none + BEGIN_DOC + ! Number of Det_a x Det_b products in input file + END_DOC + call get_spindeterminants_n_det(det_num_input) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, det_alpha_norm, (det_alpha_num) ] +&BEGIN_PROVIDER [ double precision, det_beta_norm, (det_beta_num) ] + implicit none + BEGIN_DOC + ! Norm of the alpha and beta spin determinants in the wave function: + ! + ! ||Da||_i \sum_j C_{ij}**2 + END_DOC + + integer :: i,j,k + double precision :: f + + det_alpha_norm = 0.d0 + det_beta_norm = 0.d0 + do k=1,det_num + i = det_coef_matrix_rows(k) + j = det_coef_matrix_columns(k) + f = det_coef_matrix_values(k)*det_coef_matrix_values(k) + det_alpha_norm(i) += f + det_beta_norm(j) += f + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, det_coef_matrix_values, (det_num_input) ] +&BEGIN_PROVIDER [ integer, det_coef_matrix_rows, (det_num_input) ] +&BEGIN_PROVIDER [ integer, det_coef_matrix_columns, (det_num_input) ] + implicit none + BEGIN_DOC + ! det_coef_matrix in sparse storage (Coordinate format for sparse BLAS) + END_DOC + double precision, allocatable :: buffer(:,:) + allocate (buffer(det_num_input,N_states)) + call get_spindeterminants_psi_coef_matrix_rows(det_coef_matrix_rows) + call get_spindeterminants_psi_coef_matrix_columns(det_coef_matrix_columns) + call get_spindeterminants_psi_coef_matrix_values(buffer) + det_coef_matrix_values(:) = buffer(:,i_state) + deallocate(buffer) +END_PROVIDER + +BEGIN_PROVIDER [ integer, det_num ] + implicit none + BEGIN_DOC + ! Number of Det_a x Det_b products. The determinant basis set is reduced with + ! the CI threshold + END_DOC + integer :: i,j,k,l + double precision :: f + + double precision :: d_alpha(det_alpha_num), d_beta (det_beta_num) + integer :: i_alpha(det_alpha_num), i_beta(det_beta_num) + integer :: iorder(max(det_alpha_num,det_beta_num)) + double precision :: t, norm + + t = ci_threshold + + ! Compute the norm of the alpha and beta determinants + d_alpha = 0.d0 + d_beta = 0.d0 + do k=1,det_num_input + i = det_coef_matrix_rows(k) + j = det_coef_matrix_columns(k) + f = det_coef_matrix_values(k)*det_coef_matrix_values(k) + d_alpha(i) += f + d_beta (j) += f + enddo + t = min(t, maxval(d_alpha)) + t = min(t, maxval(d_beta)) + + ! Reorder alpha determinants + do i=1,det_alpha_num + iorder(i) = i + if (d_alpha(i) < t) then + i_alpha(i) = det_alpha_num+i + else + i_alpha(i) = i + endif + enddo + call isort(i_alpha,iorder,det_alpha_num) + + i=det_alpha_num + do while (i > 0) + if (i_alpha(i) <= det_alpha_num) then + det_alpha_num = i + exit + else + i = i-1 + endif + enddo + + do i=1,det_alpha_num + psi_det_alpha(:,i) = psi_det_alpha(:,iorder(i)) + i_alpha(iorder(i)) = i + enddo + + ! Reorder beta determinants + do i=1,det_beta_num + iorder(i) = i + if (d_beta(i) < t) then + i_beta(i) = det_beta_num+i + else + i_beta(i) = i + endif + enddo + call isort(i_beta,iorder,det_beta_num) + + + i=det_beta_num + do while (i > 0) + if (i_beta(i) <= det_beta_num) then + det_beta_num = i + exit + else + i = i-1 + endif + enddo + + do i=1,det_beta_num + psi_det_beta(:,i) = psi_det_beta(:,iorder(i)) + i_beta(iorder(i)) = i + enddo + + + ! Apply the threshold to the wave function + l = 1 + norm = 0.d0 + do k=1,det_num_input + i = det_coef_matrix_rows(k) + j = det_coef_matrix_columns(k) + det_coef_matrix_rows(l) = i_alpha(i) + det_coef_matrix_columns(l) = i_beta(j) + det_coef_matrix_values(l) = det_coef_matrix_values(k) + if ( (d_alpha(i) >= t).and.(d_beta(j) >= t) ) then + l = l+1 + norm += det_coef_matrix_values(k)*det_coef_matrix_values(k) + endif + enddo + det_num = l-1 + norm = 1.d0/dsqrt(norm) + do k=1,det_num + det_coef_matrix_values(k) *= norm + enddo + + SOFT_TOUCH det_alpha_num det_beta_num det_coef_matrix_values det_coef_matrix_rows det_coef_matrix_columns psi_det_beta psi_det_alpha + +END_PROVIDER + + BEGIN_PROVIDER [ integer, det_alpha_num ] +&BEGIN_PROVIDER [ integer, det_beta_num ] + implicit none + BEGIN_DOC + ! Number of alpha and beta determinants + END_DOC + call get_spindeterminants_n_det_alpha(det_alpha_num) + call get_spindeterminants_n_det_beta(det_beta_num) +END_PROVIDER + + BEGIN_PROVIDER [ integer, det_alpha_num_8 ] +&BEGIN_PROVIDER [ integer, det_beta_num_8 ] + implicit none + BEGIN_DOC + ! Number of alpha and beta determinants + END_DOC + integer :: mod_align + det_alpha_num_8 = max(4,mod_align(det_alpha_num)) ! + det_beta_num_8 = max(4,mod_align(det_beta_num)) ! Used in 4x unrolling +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ci_threshold ] + + implicit none + BEGIN_DOC + ! Threshold on absolute value of the CI coefficients of the wave functioE + END_DOC + ci_threshold = 0.d0 + call get_simulation_ci_threshold(ci_threshold) + call rinfo(irp_here,'ci_threshold',ci_threshold) + +END_PROVIDER + +BEGIN_PROVIDER [ integer*8, psi_det_alpha, (N_int,det_alpha_num) ] + implicit none + BEGIN_DOC + ! Alpha determinants + END_DOC + call get_spindeterminants_psi_det_alpha(psi_det_alpha) +END_PROVIDER + +BEGIN_PROVIDER [ integer*8, psi_det_beta, (N_int,det_beta_num) ] + implicit none + BEGIN_DOC + ! Beta determinants + END_DOC + call get_spindeterminants_psi_det_beta(psi_det_beta) +END_PROVIDER + +BEGIN_PROVIDER [ integer, present_mos, (mo_tot_num) ] +&BEGIN_PROVIDER [ integer, num_present_mos ] +&BEGIN_PROVIDER [ integer, num_present_mos_8 ] +&BEGIN_PROVIDER [ integer, mo_closed_num ] + implicit none + BEGIN_DOC + ! List of used MOs to build the wf in the CI expansion + END_DOC + integer*8 :: tmp_det(N_int) + integer :: i,k + integer, external :: mod_align + PROVIDE det_num + + num_present_mos = mo_tot_num + do i=1,mo_tot_num + present_mos(i) = i + enddo + +!--- + present_mos = 0 + tmp_det = 0_8 + do i=1,det_alpha_num + do k=1,N_int + tmp_det(k) = ior(tmp_det(k),psi_det_alpha(k,i)) + enddo + enddo + do i=1,det_beta_num + do k=1,N_int + tmp_det(k) = ior(tmp_det(k),psi_det_beta(k,i)) + enddo + enddo + call bitstring_to_list(tmp_det,present_mos,num_present_mos,N_int) +!--- + + num_present_mos_8 = mod_align(num_present_mos) + + integer :: list(mo_tot_num), n + logical :: good + + list = present_mos + mo_closed_num = elec_beta_num + do n=1,elec_beta_num + call list_to_bitstring(tmp_det,present_mos,n,N_int) + do k=1,N_int + if (tmp_det(k) == 0_8) then + exit + endif + good = .True. + do i=1,det_alpha_num + if (iand(tmp_det(k),psi_det_alpha(k,i)) /= tmp_det(k)) then + good = .False. + exit + endif + enddo + if (good) then + do i=1,det_beta_num + if (iand(tmp_det(k),psi_det_beta(k,i)) /= tmp_det(k)) then + good = .False. + exit + endif + enddo + endif + if (.not.good) then + exit + endif + enddo + if (.not.good) then + mo_closed_num = n-1 + exit + endif + enddo +END_PROVIDER + +subroutine list_to_bitstring( string, list, n_elements, Nint) + implicit none + BEGIN_DOC + ! Returns the physical string "string(N_int,2)" from the array of + ! occupations "list(N_int*64,2) + END_DOC + integer, intent(in) :: Nint + integer*8, intent(out) :: string(Nint) + integer, intent(in) :: list(Nint*64) + integer, intent(in) :: n_elements + + + integer :: i, j + integer :: ipos, iint + + ! + ! <== ipos ==> + ! | + ! v + !string :|------------------------|-------------------------|------------------------| + ! <==== 64 ====> <==== 64 ====> <==== 64 ====> + ! { iint } { iint } { iint } + ! + + string = 0_8 + + do i=1,n_elements + iint = ishft(list(i)-1,-6) + 1 + ipos = list(i)-ishft((iint-1),6)-1 + string(iint) = ibset( string(iint), ipos ) + enddo + +end + + +BEGIN_PROVIDER [ integer, det_alpha_order, (det_alpha_num) ] + implicit none + BEGIN_DOC + ! Order in which to compute the alhpa determinants + END_DOC + integer :: i +! double precision :: tmp(det_alpha_num) + do i=1,det_alpha_num + det_alpha_order(i) = i + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, det_beta_order, (det_beta_num) ] + implicit none + BEGIN_DOC + ! Order in which to compute the beta determinants + END_DOC + integer :: i + do i=1,det_beta_num + det_beta_order(i) = i + enddo +END_PROVIDER +