diff --git a/.gitignore b/.gitignore index 65d0183..d36ee77 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ *.o *.mod -cppSherman-Morrison -fSherman-Morrison -.vscode \ No newline at end of file +.vscode +cMaponiA3_test +fMaponiA3_test +QMCChem_dataset_test + diff --git a/Helpers.hpp b/Helpers.hpp index 176e97c..d0e6ac5 100644 --- a/Helpers.hpp +++ b/Helpers.hpp @@ -119,3 +119,15 @@ T matDet(T **A, unsigned int M) { } delete [] temp; } + + +template +bool is_identity(T *A, unsigned int M, double tolerance) { + for (unsigned int i = 0; i < M; i++) { + for (unsigned int j = 0; j < M; j++) { + if (i==j && fabs(A[i*M+j]-1) > tolerance) return false; + if (i!=j && fabs(A[i*M+j]) > tolerance) return false; + } + } + return true; +} diff --git a/Makefile b/Makefile index e913743..4d2890d 100644 --- a/Makefile +++ b/Makefile @@ -1,44 +1,56 @@ ## Used compilers -CXX = icpc -FC = ifort +H5CXX = h5c++ +CXX = clang++ +FC = flang ## Compiler flags -CXXFLAGS = -O0 #-debug full -traceback -FFLAGS = -O0 #-debug full -traceback -# ARCH = -xCORE-AVX2 +CXXFLAGS = -O0 +FFLAGS = -O0 -## Deps & objs for the C++ stuff -cppDEPS = cppmain.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp -cppOBJ = cppmain.o SM_MaponiA3.o +## Deps & objs for C++ cMaponiA3_test +cMaponiA3_testDEP = cMaponiA3_test.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp +cMaponiA3_testOBJ = cMaponiA3_test.o SM_MaponiA3.o -## Deps & objs for the Fortran stuff -fDEPS = fmain.f90 SM_MaponiA3_mod.f90 -fOBJ = SM_MaponiA3.o SM_MaponiA3_mod.o fmain.o -fLIBS = -lstdc++ +## Deps & objs for Fortran fMaponiA3_test +fMaponiA3_testDEP = fMaponiA3_test.f90 SM_MaponiA3_mod.f90 +fMaponiA3_testOBJ = SM_MaponiA3.o SM_MaponiA3_mod.o fMaponiA3_test.o +fMaponiA3_testLIB = -lstdc++ -## Compile recipes for C++ stuff -%.o: %.cpp $(cppDEPS) +## Deps & objs for Fortran QMCChem_dataset_test +QMCChem_dataset_testDEP = QMCChem_dataset_test.f90 SM_MaponiA3_mod.f90 Utils_mod.f90 +QMCChem_dataset_testOBJ = SM_MaponiA3.o Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o +QMCChem_dataset_testLIB = -lstdc++ + +## Compile recipes for C++ cMaponiA3_test +%.o: %.cpp $(cMaponiA3_testDEP) $(CXX) $(ARCH) $(CXXFLAGS) -c -o $@ $< -## Compile recepies for Fortran stuff -%.o: %.f90 $(fDEPS) +## Compile recepies for Fortran fMaponiA3_test +%.o: %.f90 $(fMaponiA3_testDEP) $(FC) $(ARCH) $(FFLAGS) -c -o $@ $< ## Build tagets .PHONY: all clean distclean -all: cppSherman-Morrison fSherman-Morrison +all: cMaponiA3_test fMaponiA3_test QMCChem_dataset_test tests/test clean: @rm -vf *.o *.mod distclean: clean - @rm -vf cppSherman-Morrison fSherman-Morrison + @rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test ## Linking the C++ example program -cppSherman-Morrison: $(cppOBJ) +cMaponiA3_test: $(cMaponiA3_testOBJ) $(CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ ## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' -fSherman-Morrison: $(fOBJ) - $(FC) $(ARCH) $(FFLAGS) $(fLIBS) -o $@ $^ +fMaponiA3_test: $(fMaponiA3_testOBJ) + $(FC) $(ARCH) $(FFLAGS) $(fMaponiA3_testLIB) -o $@ $^ + +## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' +QMCChem_dataset_test: $(QMCChem_dataset_testOBJ) + $(FC) $(ARCH) $(FFLAGS) $(QMCChem_dataset_testLIB) -o $@ $^ + +tests/test: tests/test.cpp SM_MaponiA3.o + $(H5CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ diff --git a/QMCChem_dataset_test.f90 b/QMCChem_dataset_test.f90 new file mode 100644 index 0000000..d53bd01 --- /dev/null +++ b/QMCChem_dataset_test.f90 @@ -0,0 +1,66 @@ +program QMCChem_dataset_test + use Sherman_Morrison, only : MaponiA3 + use Utils, only : Read_dataset + use, intrinsic :: iso_c_binding, only : c_int, c_double + implicit none + + integer :: i, j, col !! Iterators + integer :: cycle_id, dim, n_updates + integer(c_int), dimension(:), allocatable :: Updates_index + real(c_double), dimension(:,:), allocatable :: S, S_inv, Updates + + call Read_dataset("test.dataset.dat", & + cycle_id, & + dim, & + n_updates, & + S, & + S_inv, & + Updates_index, & + Updates) + + !! Write current S and S_inv to file for check in Octave + open(unit = 2000, file = "Slater_old.dat") + open(unit = 3000, file = "Slater_inv_old.dat") + do i=1,dim + do j=1,dim + write(2000,"(E23.15, 1X)", advance="no") S(i,j) + write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) + end do + write(2000,*) + write(3000,*) + end do + call flush(2000) + call flush(3000) + close(2000) + close(3000) + + !! Send S, S_inv and Updates to MaponiA3 algo + call MaponiA3(S, S_inv, dim, n_updates, Updates, Updates_index) + + !! Update S itself + do j=1,n_updates + do i=1,dim + col = Updates_index(j) + S(i,col) = S(i,col) + Updates(i,col) + end do + end do + + !! Write new S and S_inv to file for check in Octave + open(unit = 2000, file = "Slater_new.dat") + open(unit = 3000, file = "Slater_inv_new.dat") + do i=1,dim + do j=1,dim + write(2000,"(E23.15, 1X)", advance="no") S(i,j) + write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) + end do + write(2000,*) + write(3000,*) + end do + call flush(2000) + call flush(3000) + + close(2000) + close(3000) + + deallocate(S, S_inv, Updates, Updates_index) +end program diff --git a/SM_MaponiA3.cpp b/SM_MaponiA3.cpp index 3a54730..f24a831 100644 --- a/SM_MaponiA3.cpp +++ b/SM_MaponiA3.cpp @@ -28,7 +28,7 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, } // Calculate all the y0k in M^2 multiplications instead of M^3 - for (k = 1; k < M + 1; k++) { + for (k = 1; k < N_updates + 1; k++) { for (i = 1; i < M + 1; i++) { ylk[0][k][i] = Slater_inv[(i - 1) * M + (i - 1)] * Updates[(i - 1) * M + (k - 1)]; @@ -92,10 +92,11 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, } extern "C" { -void MaponiA3_f(double **linSlater0, double **linSlater_inv, unsigned int *Dim, - unsigned int *N_updates, double **linUpdates, - unsigned int **Updates_index) { - MaponiA3(*linSlater0, *linSlater_inv, *Dim, *N_updates, *linUpdates, - *Updates_index); -} + void MaponiA3_f(double **linSlater0, double **linSlater_inv, unsigned int *Dim, + unsigned int *N_updates, double **linUpdates, + unsigned int **Updates_index) { + MaponiA3(*linSlater0, *linSlater_inv, *Dim, + *N_updates, *linUpdates, + *Updates_index); + } } diff --git a/Utils_mod.f90 b/Utils_mod.f90 new file mode 100644 index 0000000..d933dd8 --- /dev/null +++ b/Utils_mod.f90 @@ -0,0 +1,47 @@ +module Utils + implicit none + contains + subroutine Read_dataset(filename, cycle_id, dim, n_updates, & + S, S_inv, Updates_index, Updates) + use, intrinsic :: iso_c_binding, only : c_int, c_double + implicit none + + character (len = *), intent(in) :: filename + integer, intent(inout) :: cycle_id, dim, n_updates + integer(c_int), allocatable, intent(inout) :: Updates_index(:) + real(c_double), allocatable, intent(inout) :: S(:,:), S_inv(:,:) + real(c_double), allocatable, intent(inout) :: Updates(:,:) + integer :: i, j + character (len = 32) :: ignore + + !! Start of reading the dataset from file + open(unit = 1000, file = filename) + read(1000,*) + read(1000,*) ignore, cycle_id + read(1000,*) ignore, dim + read(1000,*) ignore,n_updates + + allocate(Updates_index(n_updates), S(dim,dim), & + S_inv(dim,dim), Updates(dim,n_updates)) + + !! Read S and S_inv + read(1000,*) + do i=1,dim + do j=1,dim + read(1000,*) ignore, ignore, S(i,j), S_inv(i,j) + end do + end do + + !! Read the updates Updates and Updates_index + do j=1,n_updates + read(1000,*) ignore, Updates_index(j) + do i=1,dim + read(1000,*) ignore, Updates(i,j) + end do + end do + + read(1000,*) ignore + close(1000) + !! End of reading the dataset from file + end subroutine Read_dataset +end module Utils \ No newline at end of file diff --git a/cppmain.cpp b/cMaponiA3_test.cpp similarity index 100% rename from cppmain.cpp rename to cMaponiA3_test.cpp diff --git a/det.irp.f b/det.irp.f new file mode 100644 index 0000000..10a0968 --- /dev/null +++ b/det.irp.f @@ -0,0 +1,1966 @@ +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;; + 51;; + 52;; + 53;; + 54;; + 55;; + 56;; + 57;; + 58;; + 59;; + 60;; + 61;; + 62;; + 63;; + 64;; + 65;; + 66;; + 67;; + 68;; + 69;; + 70;; + 71;; + 72;; + 73;; + 74;; + 75;; + 76;; + 77;; + 78;; + 79;; + 80;; + 81;; + 82;; + 83;; + 84;; + 85;; + 86;; + 87;; + 88;; + 89;; + 90;; + 91;; + 92;; + 93;; + 94;; + 95;; + 96;; + 97;; + 98;; + 99;; + 100;; + 101;; + 102;; + 103;; + 104;; + 105;; + 106;; + 107;; + 108;; + 109;; + 110;; + 111;; + 112;; + 113;; + 114;; + 115;; + 116;; + 117;; + 118;; + 119;; + 120;; + 121;; + 122;; + 123;; + 124;; + 125;; + 126;; + 127;; + 128;; + 129;; + 130;; + 131;; + 132;; + 133;; + 134;; + 135;; + 136;; + 137;; + 138;; + 139;; + 140;; + 141;; + 142;; + 143;; + 144;; + 145;; + 146;; + 147;; + 148;; + 149;; + 150;; + 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 + 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 + + !DIR$ VECTOR ALIGNED + 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 + double precision :: zj, zj1, zj2, zj3 + + !DIR$ NOPREFETCH + !DIR$ SIMD NOVECREMAINDER + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + zj = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + do i=1,$n-1,4 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + enddo + + d_inv = 1.d0/d + d = d+zj + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + !DIR$ VECTOR ALIGNED + do j=1,$n,4 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + zj3 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + do i=1,$n + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + zj3 = zj3 + S_inv(i,j+3)*u(i) + enddo + z(j ) = zj + z(j+1) = zj1 + z(j+2) = zj2 + z(j+3) = zj3 + enddo + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n,4 + zj = z(i ) + zj1 = z(i+1) + zj2 = z(i+2) + zj3 = z(i+3) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 + S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 + enddo + enddo + +end + +SUBST [ n ] +8 ;; +12 ;; +16 ;; +20 ;; +24 ;; +28 ;; +32 ;; +36 ;; +40 ;; +44 ;; +48 ;; +52 ;; +56 ;; +60 ;; +64 ;; +68 ;; +72 ;; +76 ;; +80 ;; +84 ;; +88 ;; +92 ;; +96 ;; +100 ;; +104 ;; +108 ;; +112 ;; +116 ;; +120 ;; +124 ;; +128 ;; +132 ;; +136 ;; +140 ;; +144 ;; +148 ;; + +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 + double precision :: zj, zj1, zj2, zj3 + + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + zj = 0.d0 + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) + do i=1,$n-1,4 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + enddo + zj = zj + S_inv($n,l)*u($n) + + d_inv = 1.d0/d + d = d+zj + lambda = d*d_inv + if (dabs(lambda) < 1.d-6) then + d = 0.d0 + write(502,"('#BREAKDOWN_OCCURED')") + return + end if + + !DIR$ VECTOR ALIGNED + do j=1,$n-1,4 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + zj3 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + do i=1,$n-1 + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + zj3 = zj3 + S_inv(i,j+3)*u(i) + enddo + z(j ) = zj + S_inv($n,j )*u($n) + z(j+1) = zj1 + S_inv($n,j+1)*u($n) + z(j+2) = zj2 + S_inv($n,j+2)*u($n) + z(j+3) = zj3 + S_inv($n,j+3)*u($n) + enddo + + zj = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + do i=1,$n-1 + zj = zj + S_inv(i,$n)*u(i) + enddo + z($n) = zj + S_inv($n,$n)*u($n) + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n-1,4 + zj = z(i ) + zj1 = z(i+1) + zj2 = z(i+2) + zj3 = z(i+3) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + do j=1,$n-1 + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 + S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 + enddo + S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj + S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*zj1 + S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*zj2 + S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*zj3 + enddo + + zj = z($n) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER + do i=1,$n + S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*zj + enddo + + +end + +SUBST [ n ] +5 ;; +9 ;; +13 ;; +17 ;; +21 ;; +25 ;; +29 ;; +33 ;; +37 ;; +41 ;; +45 ;; +49 ;; +53 ;; +57 ;; +61 ;; +65 ;; +69 ;; +73 ;; +77 ;; +81 ;; +85 ;; +89 ;; +93 ;; +97 ;; +101 ;; +105 ;; +109 ;; +113 ;; +117 ;; +121 ;; +125 ;; +129 ;; +133 ;; +137 ;; +141 ;; +145 ;; +149 ;; + +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 + + double precision :: zj, zj1, zj2, zj3 + !DIR$ NOPREFETCH + !DIR$ SIMD NOVECREMAINDER + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + zj = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + do i=1,$n-2,4 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + enddo + i=$n-1 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + + d_inv = 1.d0/d + d = d+zj + lambda = d*d_inv + if (dabs(lambda) < 1.d-3) then + d = 0.d0 + return + endif + + !DIR$ VECTOR ALIGNED + do j=1,$n-2,4 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + zj3 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER + do i=1,$n-2 + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + zj3 = zj3 + S_inv(i,j+3)*u(i) + enddo + z(j ) = zj + S_inv($n-1,j )*u($n-1) + z(j ) = z(j ) + S_inv($n,j )*u($n) + z(j+1) = zj1 + S_inv($n-1,j+1)*u($n-1) + z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) + z(j+2) = zj2 + S_inv($n-1,j+2)*u($n-1) + z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) + z(j+3) = zj3 + S_inv($n-1,j+3)*u($n-1) + z(j+3) = z(j+3) + S_inv($n,j+3)*u($n) + enddo + + j=$n-1 + zj = 0.d0 + zj1 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER + do i=1,$n-2 + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + enddo + z(j ) = zj + S_inv($n-1,j )*u($n-1) + z(j ) = z(j ) + S_inv($n,j )*u($n) + z(j+1) = zj1 + S_inv($n-1,j+1)*u($n-1) + z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n-2,4 + zj = z(i) + zj1 = z(i+1) + zj2 = z(i+2) + zj3 = z(i+3) + !DIR$ VECTOR ALIGNED + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER + do j=1,$n-2 + S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) + S_inv(j,i+2) = S_inv(j,i+2)*lambda -zj2*w(j) + S_inv(j,i+3) = S_inv(j,i+3)*lambda -zj3*w(j) + enddo + S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj *w($n-1) + S_inv($n ,i ) = S_inv($n ,i )*lambda -zj *w($n ) + S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) + S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) + S_inv($n-1,i+2) = S_inv($n-1,i+2)*lambda -zj2*w($n-1) + S_inv($n ,i+2) = S_inv($n ,i+2)*lambda -zj2*w($n ) + S_inv($n-1,i+3) = S_inv($n-1,i+3)*lambda -zj3*w($n-1) + S_inv($n ,i+3) = S_inv($n ,i+3)*lambda -zj3*w($n ) + enddo + + i=$n-1 + zj = z(i) + zj1= z(i+1) + !DIR$ VECTOR ALIGNED + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1) + do j=1,$n-2 + S_inv(j,i ) = S_inv(j,i )*lambda -zj*w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) + enddo + S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj*w($n-1) + S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) + S_inv($n ,i ) = S_inv($n ,i )*lambda -zj*w($n ) + S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) + +end + +SUBST [ n ] +6 ;; +10 ;; +14 ;; +18 ;; +22 ;; +26 ;; +30 ;; +34 ;; +38 ;; +42 ;; +46 ;; +50 ;; +54 ;; +58 ;; +62 ;; +66 ;; +70 ;; +74 ;; +78 ;; +82 ;; +86 ;; +90 ;; +94 ;; +98 ;; +102 ;; +106 ;; +110 ;; +114 ;; +118 ;; +122 ;; +126 ;; +130 ;; +134 ;; +138 ;; +142 ;; +146 ;; +150 ;; + +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 + + double precision :: zj, zj1, zj2, zj3 + + !DIR$ SIMD + do i=1,$n + u(i) = m(i) - S(i,l) + enddo + + zj = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER + do i=1,$n-3,4 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & + + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) + enddo + i=$n-2 + zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + S_inv(i+2,l)*u(i+2) + + + d_inv = 1.d0/d + d = d+zj + 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 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + zj3 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + do i=1,$n-3 + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + zj3 = zj3 + S_inv(i,j+3)*u(i) + enddo + z(j ) = zj + S_inv($n-2,j )*u($n-2) + z(j ) = z(j ) + S_inv($n-1,j )*u($n-1) + z(j ) = z(j ) + S_inv($n,j )*u($n) + z(j+1) = zj1 + S_inv($n-2,j+1)*u($n-2) + z(j+1) = z(j+1) + S_inv($n-1,j+1)*u($n-1) + z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) + z(j+2) = zj2 + S_inv($n-2,j+2)*u($n-2) + z(j+2) = z(j+2) + S_inv($n-1,j+2)*u($n-1) + z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) + z(j+3) = zj3 + S_inv($n-2,j+3)*u($n-2) + z(j+3) = z(j+3) + S_inv($n-1,j+3)*u($n-1) + z(j+3) = z(j+3) + S_inv($n,j+3)*u($n) + enddo + + j=$n-2 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2) + do i=1,$n-3 + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + enddo + z(j ) = zj + S_inv($n-2,j )*u($n-2) + z(j ) = z(j ) + S_inv($n-1,j )*u($n-1) + z(j ) = z(j ) + S_inv($n,j )*u($n) + z(j+1) = zj1 + S_inv($n-2,j+1)*u($n-2) + z(j+1) = z(j+1) + S_inv($n-1,j+1)*u($n-1) + z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) + z(j+2) = zj2 + S_inv($n-2,j+2)*u($n-2) + z(j+2) = z(j+2) + S_inv($n-1,j+2)*u($n-1) + z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) + do i=1,$n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,$n-3,4 + zj = z(i) + zj1 = z(i+1) + zj2 = z(i+2) + zj3 = z(i+3) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + do j=1,$n-3 + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 + S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 + enddo + S_inv($n-2,i ) = S_inv($n-2,i )*lambda -zj *w($n-2) + S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj *w($n-1) + S_inv($n ,i ) = S_inv($n ,i )*lambda -zj *w($n ) + S_inv($n-2,i+1) = S_inv($n-2,i+1)*lambda -zj1*w($n-2) + S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) + S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) + S_inv($n-2,i+2) = S_inv($n-2,i+2)*lambda -zj2*w($n-2) + S_inv($n-1,i+2) = S_inv($n-1,i+2)*lambda -zj2*w($n-1) + S_inv($n ,i+2) = S_inv($n ,i+2)*lambda -zj2*w($n ) + S_inv($n-2,i+3) = S_inv($n-2,i+3)*lambda -zj3*w($n-2) + S_inv($n-1,i+3) = S_inv($n-1,i+3)*lambda -zj3*w($n-1) + S_inv($n ,i+3) = S_inv($n ,i+3)*lambda -zj3*w($n ) + enddo + + i=$n-2 + zj = z(i) + zj1 = z(i+1) + zj2 = z(i+2) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2) + do j=1,$n + S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj + S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 + S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 + enddo + + +end + +SUBST [ n ] +7 ;; +11 ;; +15 ;; +19 ;; +23 ;; +27 ;; +31 ;; +35 ;; +39 ;; +43 ;; +47 ;; +51 ;; +55 ;; +59 ;; +63 ;; +67 ;; +71 ;; +75 ;; +79 ;; +83 ;; +87 ;; +91 ;; +95 ;; +99 ;; +103 ;; +107 ;; +111 ;; +115 ;; +119 ;; +123 ;; +127 ;; +131 ;; +135 ;; +139 ;; +143 ;; +147 ;; + +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) + !DIR$ ASSUME (n>150) + + integer :: i,j,n4 + double precision :: zl + + !DIR$ NOPREFETCH + do i=1,n + u(i) = m(i) - S(i,l) + enddo + + zl = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zl) + do i=1,n + zl = zl + S_inv(i,l)*u(i) + enddo + + d_inv = 1.d0/d + d = d+zl + lambda = d*d_inv + + if ( dabs(lambda) < 1.d-3 ) then + d = 0.d0 + endif + + double precision :: zj, zj1, zj2, zj3 + + n4 = iand(n,not(3)) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + do j=1,n4,4 + zj = 0.d0 + zj1 = 0.d0 + zj2 = 0.d0 + zj3 = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) + do i=1,n + zj = zj + S_inv(i,j )*u(i) + zj1 = zj1 + S_inv(i,j+1)*u(i) + zj2 = zj2 + S_inv(i,j+2)*u(i) + zj3 = zj3 + S_inv(i,j+3)*u(i) + enddo + z(j ) = zj + z(j+1) = zj1 + z(j+2) = zj2 + z(j+3) = zj3 + enddo + + do j=n4+1,n + zj = 0.d0 + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD REDUCTION(+:zj) + do i=1,n + zj = zj + S_inv(i,j)*u(i) + enddo + z(j ) = zj + enddo + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) + do i=1,n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(d_inv) + do i=1,n + w(i) = S_inv(i,l)*d_inv + S(i,l) = m(i) + enddo + + do i=1,n4,4 + zj = z(i) + zj1 = z(i+1) + zj2 = z(i+2) + zj3 = z(i+3) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) + do j=1,n + S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) + S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) + S_inv(j,i+2) = S_inv(j,i+2)*lambda -zj2*w(j) + S_inv(j,i+3) = S_inv(j,i+3)*lambda -zj3*w(j) + enddo + enddo + + do i=n4+1,n + zj = z(i) + !DIR$ VECTOR ALIGNED + !DIR$ NOPREFETCH + !DIR$ SIMD FIRSTPRIVATE(lambda,zj) + do j=1,n + S_inv(j,i) = S_inv(j,i)*lambda -zj*w(j) + enddo + 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 ) + if (l /= elec_alpha_num) then + stop 'error in number of alpha electrons' + endif + +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 ) + if (l /= elec_beta_num) then + stop 'error in number of beta electrons' + endif +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 + ! + ! slater_matrix_alpha : 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 * 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 + integer, save :: cycle_id=0 + + !! Some usefull formats for output + 10001 format ('#START_PACKET') + 10008 format ('#CYCLE_ID: ', I4) + 10000 format ('#SLATER_MATRIX_DIM: ', I3) + 10002 format ('#NUPDATES: ', I2) + 10003 format ('#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet') + 10004 format ('(',I0.2,',',I0.2,')',2(2X,E23.15)) + 10005 format ('#COL_UPDATE_INDEX: ', I2) + 10006 format ('#COL_UPDATE_COMP_(',I0.2,'): ', E23.15) + 10007 format ('#END_PACKET',/) + + open (unit = 501, file = "dataset.dat") !! slightly cleaner output + open (unit = 502, file = "dataset.fulltrace.dat") + + if (ifirst == 0) then !! If this is the first time we enter this subroutine + 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 + ! write(*,*) "det_i: ", det_i +! if (det_i == -1 ) then + + !! determin number of updates, the updates and + 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 + + !! Write number of alpha electrons, number of updates to do, slater and slater_inv to file unit 10000 + write(501,10001) + write(501,10008) cycle_id + write(501,10000) elec_alpha_num + write(501,10002) n_to_do + write(501,10003) + do i=1,elec_alpha_num + do j=1,elec_alpha_num + write(501,10004) i, j, slater_matrix_alpha(i,j), slater_matrix_alpha_inv_det(i,j) / det_alpha_value_curr + end do + end do + + !! write all the updates to file unit 10000 + do l=1,n_to_do + k = to_do(l) + imo = mo_list_alpha_curr(k) + write(501,10005) k + do i=1,elec_alpha_num + write(501,10006) i, mo_value(i, imo) + end do + end do + +! print n_to_do (number of updates to do) +! to_do (array of the columns that need to be swapped) +! mo_list_alpha_curr (list of orbitals to build the current determinant) +! mo_list_alpha_prev (list of the previous determinant) +! +! slater_matrix_alpha (slater matrix that needs to be inverted. This is the initial matrix) +! slater_matrix_alpha_inv_det = inverse of the slater matrix divided by the determinant (ddet: l. 1216) +! +! +! +! 1 2 3 4 +! -------- +! 2 4 6 8 +! 2 4 10 8 +! +! mo_list(:) (2,4,10,8) +! mo_list_prev(:) (2,4,6,8) +! n_todo 1 +! to_do (3) +! +! +! 1 2 3 4 +! -------- +! 2 4 6 8 +! 2 5 10 8 +! +! mo_list(:) (2,5,10,8) +! mo_list_prev(:) (2,4,6,8) +! n_todo 2 +! to_do (2,3) + + + ddet = 0.d0 + + if (n_to_do < shiftl(elec_alpha_num,1)) then + + write(502,10001) + write(502,10008) cycle_id + write(502,10000) elec_alpha_num + write(502,10002) n_to_do + + 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) + + write(502,10003) + do i=1,elec_alpha_num + do j=1,elec_alpha_num + write(502,10004) i, j, slater_matrix_alpha(i,j), slater_matrix_alpha_inv_det(i,j) / ddet + end do + end do + write(502,10005) k + do i=1,elec_alpha_num + write(502,10006) i, mo_value(i,imo) + end do + + 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 + 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 + + write(501,10007) + write(502,10007) + cycle_id = cycle_id + 1 + !close(501) !! Close file + !close(502) !! Close file + !stop + + else + + ddet = 0.d0 + + endif + + + ! Avoid NaN + if (ddet /= 0.d0) then + continue + else + 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 < shiftl(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 + 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 + continue + else + 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 + det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:) + if (do_pseudo) then + det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) + endif + + enddo + + det_i = det_alpha_order(1) + SOFT_TOUCH det_i + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, det_beta_value, (det_beta_num_8) ] +&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_alpha_num+1:elec_num,det_beta_num) ] +&BEGIN_PROVIDER [ double precision, det_beta_pseudo, (elec_alpha_num+1:elec_num,det_beta_num_pseudo) ] + + + implicit none + + BEGIN_DOC + ! Values of the beta determinants + ! Gradients of the beta determinants + ! Laplacians of the beta determinants + END_DOC + + integer :: j,i,k + integer, save :: ifirst = 0 + if (elec_beta_num == 0) then + det_beta_value = 1.d0 + return + endif + + if (ifirst == 0) then + ifirst = 1 + det_beta_value = 0.d0 + det_beta_grad_lapl = 0.d0 + det_beta_pseudo = 0.d0 + endif + + do j=1,det_beta_num + + det_j = det_beta_order(j) + if (j > 1) then + TOUCH det_j + endif + + det_beta_value(det_j) = det_beta_value_curr + det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:) + if (do_pseudo) then + det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) + endif + + enddo + + det_j = det_beta_order(1) + SOFT_TOUCH det_j + +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 + + if (det_num < shiftr(det_alpha_num*det_beta_num,2)) then + + 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 + + else + + call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & + size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1) + call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & + size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1) + + endif + + ! Value + ! ----- + + psidet_value = 0.d0 + do j=1,det_beta_num + psidet_value = psidet_value + det_beta_value(j) * DaC(j) + enddo + + + if (psidet_value == 0.d0) then + call abrt(irp_here,'Determinantal component of the wave function is zero.') + endif + psidet_inv = 1.d0/psidet_value + + ! Gradients + ! --------- + + call dgemv('N',elec_alpha_num*4,det_alpha_num,1.d0, & + det_alpha_grad_lapl, & + size(det_alpha_grad_lapl,1)*size(det_alpha_grad_lapl,2), & + CDb, 1, 0.d0, psidet_grad_lapl, 1) + if (elec_beta_num /= 0) then + call dgemv('N',elec_beta_num*4,det_beta_num,1.d0, & + det_beta_grad_lapl(1,elec_alpha_num+1,1), & + size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2), & + DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1) + endif + + if (do_pseudo) then + call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv, & + det_alpha_pseudo, size(det_alpha_pseudo,1), & + CDb, 1, 0.d0, pseudo_non_local, 1) + if (elec_beta_num /= 0) then + call dgemv('N',elec_beta_num,det_beta_num,psidet_inv, & + det_beta_pseudo, size(det_beta_pseudo,1), & + DaC, 1, 0.d0, pseudo_non_local(elec_alpha_num+1), 1) + endif + 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_alpha(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 + !DIR$ VECTOR ALIGNED + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo )*slater_matrix_alpha_inv_det(i ,j ) & + + mo_grad_lapl_alpha(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_alpha(k,i+1,imo )*slater_matrix_alpha_inv_det(i+1,j ) & + + mo_grad_lapl_alpha(k,i+1,imo2)*slater_matrix_alpha_inv_det(i+1,j+1) + enddo + 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 + !DIR$ VECTOR ALIGNED + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo )*slater_matrix_alpha_inv_det(i ,j ) & + + mo_grad_lapl_alpha(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_alpha(k,i+1,imo )*slater_matrix_alpha_inv_det(i+1,j ) & + + mo_grad_lapl_alpha(k,i+1,imo2)*slater_matrix_alpha_inv_det(i+1,j+1) + enddo + enddo + i=elec_alpha_num + !DIR$ VECTOR ALIGNED + do k=1,4 + det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo )*slater_matrix_alpha_inv_det(i,j ) & + + mo_grad_lapl_alpha(k,i,imo2)*slater_matrix_alpha_inv_det(i,j+1) + enddo + enddo + + j=elec_alpha_num + imo = mo_list_alpha_curr(j) + do i=1,elec_alpha_num + !DIR$ VECTOR ALIGNED + do k=1,4 + det_alpha_grad_lapl_curr(k,i ) = det_alpha_grad_lapl_curr(k,i ) + mo_grad_lapl_alpha(k,i ,imo)*slater_matrix_alpha_inv_det(i ,j) + enddo + 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_alpha(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 + !DIR$ VECTOR ALIGNED + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& + mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + enddo + 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 + !DIR$ VECTOR ALIGNED + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +& + mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j ) + & + mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) + enddo + enddo + i = elec_num + l = elec_num-elec_alpha_num + !DIR$ VECTOR ALIGNED + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j ) + & + mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) + enddo + enddo + + j = elec_beta_num + imo = mo_list_beta_curr(j) + do i=elec_alpha_num+1,elec_num + l = i-elec_alpha_num + !DIR$ VECTOR ALIGNED + do k=1,4 + det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +& + mo_grad_lapl_beta(k,i,imo)*slater_matrix_beta_inv_det(l,j) + enddo + enddo + + endif + +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 + diff --git a/fmain.f90 b/fMaponiA3_test.f90 similarity index 98% rename from fmain.f90 rename to fMaponiA3_test.f90 index 5fdc1d0..b8c519c 100644 --- a/fmain.f90 +++ b/fMaponiA3_test.f90 @@ -10,7 +10,7 @@ program Interface_test real(c_double), dimension(:,:), allocatable :: A0_inv dim = 3 - N_updates = dim + N_updates = 3 allocate(Ar_index(dim), A(dim,dim), A0(dim,dim), Ar(dim,dim), A0_inv(dim,dim)) !! Initialize A with M=3 and fill acc. to Eq. (17) from paper diff --git a/test.dataset.dat b/test.dataset.dat new file mode 100644 index 0000000..fce4db0 --- /dev/null +++ b/test.dataset.dat @@ -0,0 +1,513 @@ +#START_PACKET +#CYCLE_ID: 13 +#SLATER_MATRIX_DIM: 21 +#NUPDATES: 3 +#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet +(01,01) -0.123240642249584E+00 0.311919238733319E+02 +(01,02) 0.150952935218811E+00 0.391211616261786E+02 +(01,03) 0.871514901518822E-01 0.230468185418855E+02 +(01,04) -0.150895461440086E+00 0.385544395586345E+02 +(01,05) -0.871202647686005E-01 0.233591113858616E+02 +(01,06) -0.123183816671371E+00 -0.312005871898339E+02 +(01,07) 0.100359566509724E+00 -0.751758654676229E+00 +(01,08) -0.135056734085083E+00 -0.440356424219148E+01 +(01,09) -0.744080543518066E-01 -0.157135293238612E+01 +(01,10) -0.120210982859135E+00 0.311608268850865E+01 +(01,11) -0.757721215486526E-01 0.269147143852639E+01 +(01,12) 0.525853671133518E-01 -0.295675853033276E+01 +(01,13) 0.102125488221645E+00 0.695722472661512E+01 +(01,14) -0.307190138846636E-01 -0.656548036745546E+01 +(01,15) 0.202103536576033E-01 0.357345834272637E+01 +(01,16) -0.211704343557358E+00 -0.127226560851769E+00 +(01,17) -0.152789223939180E-01 0.211510373691347E+01 +(01,18) -0.396802611649036E-01 0.225402683839516E+01 +(01,19) 0.156803593039513E+00 0.138095041843308E+00 +(01,20) 0.273103892803192E+00 0.144755926374198E+00 +(01,21) -0.666790679097176E-01 -0.129650521353452E+01 +(02,01) -0.383377403020859E+00 0.748881862107261E+02 +(02,02) 0.469519078731537E+00 0.913799202443389E+02 +(02,03) -0.271040737628937E+00 0.511079244308798E+02 +(02,04) 0.469401627779007E+00 0.908925159786654E+02 +(02,05) -0.271042704582214E+00 0.524296869618370E+02 +(02,06) 0.383267551660538E+00 -0.743781724291279E+02 +(02,07) 0.116775326430798E+00 -0.308858192791382E+01 +(02,08) -0.936452969908714E-01 -0.883914805747968E+01 +(02,09) 0.369134694337845E-01 -0.251520306199258E+01 +(02,10) -0.204805508255959E-01 0.120501905090525E+02 +(02,11) -0.173394829034805E-01 0.968749784123027E+01 +(02,12) -0.197960838675499E+00 -0.128360153025071E+02 +(02,13) 0.962643101811409E-01 0.124602848636645E+02 +(02,14) 0.274471193552017E+00 -0.108315097729706E+02 +(02,15) 0.145023554563522E+00 0.538637730088244E+01 +(02,16) -0.179225616157055E-01 -0.155954216668878E+01 +(02,17) 0.136946782469749E+00 0.733071759611877E+01 +(02,18) -0.174913182854652E+00 0.106785970043980E+02 +(02,19) -0.128193218261003E-01 -0.526404570077434E+00 +(02,20) 0.225567314773798E-01 0.126274276019119E+01 +(02,21) -0.171542761381716E-02 -0.176385334707397E+01 +(03,01) -0.109934650361538E+00 0.981074335673308E+01 +(03,02) 0.393143796827644E-03 0.142283742204400E+02 +(03,03) 0.155300989747047E+00 0.106061203769129E+02 +(03,04) -0.378105294657871E-03 0.142205894966935E+02 +(03,05) 0.154906913638115E+00 0.106658385090178E+02 +(03,06) 0.109414055943489E+00 -0.999705856360834E+01 +(03,07) 0.162110984325409E+00 -0.235530921594611E+01 +(03,08) -0.104241691529751E+00 -0.397822961500586E+01 +(03,09) -0.199010908603668E+00 -0.615755307417915E+00 +(03,10) -0.178463548421860E+00 0.302954990045915E+01 +(03,11) 0.115249522030354E+00 0.260774312202084E+01 +(03,12) -0.751317143440247E-01 -0.324474908269315E+01 +(03,13) -0.180085301399231E-01 0.506349143477028E+01 +(03,14) 0.776397660374641E-01 -0.403331084954255E+01 +(03,15) -0.162087708711624E+00 0.216159755322281E+01 +(03,16) 0.115154005587101E+00 -0.356236831434911E+00 +(03,17) -0.111836232244968E+00 0.229287231445069E+01 +(03,18) 0.215133726596832E+00 0.345423911992016E+01 +(03,19) -0.165389522910118E+00 -0.522840315767084E-01 +(03,20) -0.160914268344641E-01 -0.798660202084652E+00 +(03,21) -0.184740740805864E-01 -0.669484792656382E-01 +(04,01) -0.505007863044739E+00 -0.509168796576331E+00 +(04,02) -0.618503987789154E+00 -0.626061168716399E+00 +(04,03) 0.357068866491318E+00 0.106282405724139E+00 +(04,04) 0.618444144725800E+00 0.187208027678628E+00 +(04,05) -0.357080936431885E+00 -0.361241794759462E+00 +(04,06) -0.504969716072083E+00 -0.152171018401898E+00 +(04,07) 0.778413712978363E-01 0.786005474778943E-02 +(04,08) 0.487910546362400E-01 0.335536466772946E-01 +(04,09) -0.115782171487808E-01 0.204150929795949E-01 +(04,10) -0.530001707375050E-01 -0.431873457736125E-01 +(04,11) 0.214109313674271E-02 -0.266642681570363E-01 +(04,12) -0.194653034210205E+00 0.354230568719380E-01 +(04,13) -0.116052091121674E+00 -0.662672792520103E-01 +(04,14) -0.260605067014694E+00 0.541727654659900E-01 +(04,15) -0.137633457779884E+00 -0.215053141567270E-01 +(04,16) 0.422291792929173E-01 -0.313441408257335E-02 +(04,17) 0.133351176977158E+00 -0.360061392608368E-01 +(04,18) -0.170262515544891E+00 -0.364432600501005E-01 +(04,19) -0.304439514875412E-01 -0.281329679904227E-02 +(04,20) 0.533868409693241E-01 0.147042058872862E-01 +(04,21) 0.596583448350430E-02 0.983501095240267E-02 +(05,01) -0.947399914264679E+00 -0.515300150381945E+00 +(05,02) 0.440122239524499E-04 -0.414148002449120E+00 +(05,03) -0.133992493152618E+01 -0.483927352922599E+00 +(05,04) 0.392828042095061E-04 -0.414236968594294E+00 +(05,05) 0.134024262428284E+01 0.892894843202583E-02 +(05,06) -0.947782099246979E+00 0.163447351211741E+00 +(05,07) -0.155572161078453E+00 -0.263796606276840E-01 +(05,08) -0.790303274989128E-01 -0.422167712445766E-01 +(05,09) -0.178044617176056E+00 -0.244456092073667E-01 +(05,10) 0.141150355339050E+00 0.343720228179847E-01 +(05,11) -0.857104435563087E-01 0.188162844929908E-01 +(05,12) 0.108973577618599E+00 -0.227730288498809E-01 +(05,13) 0.100100971758366E-01 0.389866683213947E-01 +(05,14) 0.567631311714649E-01 -0.313257678215315E-01 +(05,15) -0.203381881117821E+00 0.166784616908835E-01 +(05,16) 0.391849316656590E-01 -0.748798889931418E-02 +(05,17) 0.140034869313240E+00 0.141492354958890E-01 +(05,18) -0.172185719013214E+00 0.166134715418269E-01 +(05,19) 0.583750754594803E-01 -0.569910669664292E-02 +(05,20) -0.916111283004284E-03 -0.139449114927599E-01 +(05,21) -0.292119309306145E-01 -0.277054057678700E-02 +(06,01) -0.351941259577870E-02 0.320422300739978E+02 +(06,02) -0.429157400503755E-02 0.396796238151428E+02 +(06,03) 0.200041988864541E-02 0.225744499599300E+02 +(06,04) 0.343890837393701E-02 0.385319291283625E+02 +(06,05) -0.242624944075942E-02 0.228620731935291E+02 +(06,06) -0.279381335712969E-02 -0.315787444812651E+02 +(06,07) 0.133628234267235E+00 -0.201688119637744E+01 +(06,08) 0.164535552263260E+00 -0.691868674131645E+01 +(06,09) -0.303669814020395E-01 -0.337527478391271E+01 +(06,10) 0.335082374513149E-01 0.737932786002596E+01 +(06,11) -0.115301951766014E+00 0.532046905755991E+01 +(06,12) -0.643811970949173E-01 -0.704646020919806E+01 +(06,13) 0.168905649334192E-01 0.116798006516197E+02 +(06,14) -0.131173312664032E+00 -0.103900446585121E+02 +(06,15) -0.228150542825460E-01 0.367492551215626E+01 +(06,16) -0.140278100967407E+00 0.541292498306359E+00 +(06,17) 0.121666260063648E+00 0.688933313780047E+01 +(06,18) -0.169033091515303E-01 0.643920432296654E+01 +(06,19) 0.515771955251694E-01 0.640069090766210E+00 +(06,20) -0.180635631084442E+00 -0.275193635264416E+01 +(06,21) 0.195857465267181E+00 -0.188561818645725E+01 +(07,01) -0.520476046949625E-02 0.236539170987535E+02 +(07,02) 0.580822164192796E-02 0.295003969224203E+02 +(07,03) -0.399843230843544E-02 0.174888381981832E+02 +(07,04) 0.576473958790302E-02 0.294531380129501E+02 +(07,05) -0.271833199076355E-02 0.175044078873038E+02 +(07,06) 0.426918268203735E-02 -0.241830822724836E+02 +(07,07) 0.147665038704872E+00 0.314282931518621E+01 +(07,08) -0.117447517812252E+00 0.403784215759159E+01 +(07,09) 0.147420719265938E+00 0.169155742879292E+01 +(07,10) 0.135893806815147E+00 -0.554066860884136E+01 +(07,11) 0.422651544213295E-01 -0.525713150347206E+01 +(07,12) -0.695310905575752E-01 0.394727380633284E+01 +(07,13) -0.168554261326790E-01 -0.641137465954821E+01 +(07,14) 0.921243280172348E-01 0.423388578775427E+01 +(07,15) 0.117380328476429E+00 0.306527904979264E-01 +(07,16) 0.138013571500778E+00 0.202300380704869E+01 +(07,17) -0.623580329120159E-01 -0.459569791294739E+01 +(07,18) -0.136343047022820E+00 -0.468425326686841E+01 +(07,19) 0.131682857871056E+00 0.121902516832242E+01 +(07,20) -0.132132634520531E+00 -0.707896318054816E+00 +(07,21) 0.119739077985287E+00 0.441694177169657E+00 +(08,01) -0.241280291229486E-01 -0.964925469634583E+02 +(08,02) 0.295925457030535E-01 -0.118149667679186E+03 +(08,03) 0.161907169967890E-01 -0.672201839913842E+02 +(08,04) -0.280807800590992E-01 -0.117411734713197E+03 +(08,05) -0.170978084206581E-01 -0.686625184338945E+02 +(08,06) -0.229465700685978E-01 0.961995663618476E+02 +(08,07) 0.170763313770294E+00 0.451976509920424E+01 +(08,08) -0.246933296322823E+00 0.100355852033043E+02 +(08,09) -0.118787530809641E-01 0.362241589185164E+01 +(08,10) -0.233876928687096E-01 -0.130069484167442E+02 +(08,11) -0.235260799527168E+00 -0.110760346601338E+02 +(08,12) -0.413223467767239E-01 0.139329830045770E+02 +(08,13) 0.334977582097054E-01 -0.153873537547839E+02 +(08,14) 0.146476894617081E+00 0.146103537661863E+02 +(08,15) 0.599658722057939E-02 -0.765612032914714E+01 +(08,16) -0.962516665458679E-01 0.886407137879243E+00 +(08,17) 0.239439934492111E+00 -0.705337777438639E+01 +(08,18) -0.289920177310705E-01 -0.116989846863909E+02 +(08,19) 0.406377874314785E-01 0.592924574459365E-01 +(08,20) 0.131046742200851E+00 -0.301336277255331E+00 +(08,21) 0.136250451207161E+00 0.221249720574218E+01 +(09,01) -0.887189060449600E-02 0.408696880275784E+02 +(09,02) -0.451377010904253E-03 0.498549177991145E+02 +(09,03) -0.123341614380479E-01 0.283099915956156E+02 +(09,04) -0.489434518385679E-03 0.498508177346691E+02 +(09,05) 0.119013069197536E-01 0.288351506046625E+02 +(09,06) -0.825028307735920E-02 -0.408181390328127E+02 +(09,07) 0.120959408581257E+00 0.274581588439313E+01 +(09,08) 0.899301841855049E-01 0.442598366847707E+01 +(09,09) 0.183146566152573E+00 0.258341971052896E+01 +(09,10) -0.167249992489815E+00 -0.343796285403762E+01 +(09,11) 0.150660663843155E+00 -0.176879400649229E+01 +(09,12) 0.621877647936344E-01 0.220820305885640E+01 +(09,13) 0.775983557105064E-01 -0.388855121748206E+01 +(09,14) -0.421430803835392E-01 0.313381413329812E+01 +(09,15) -0.270811710506678E-01 -0.180792375354553E+01 +(09,16) 0.361108258366585E-01 0.773423930388479E+00 +(09,17) 0.712988025043160E-03 -0.124926313426792E+01 +(09,18) 0.177535593509674E+00 -0.149560597825620E+01 +(09,19) 0.527953356504440E-01 0.612080363385255E+00 +(09,20) 0.141406338661909E-01 0.160255194437205E+01 +(09,21) -0.784571282565594E-02 0.260363897304434E+00 +(10,01) 0.852987250254955E-05 0.301411990266088E+03 +(10,02) 0.136617054522503E-04 0.368201194080180E+03 +(10,03) 0.202008941414533E-04 0.207873209289172E+03 +(10,04) 0.554962934984360E-05 0.367592383804457E+03 +(10,05) -0.211104361369507E-04 0.212546863836424E+03 +(10,06) 0.204913503694115E-04 -0.300694117440840E+03 +(10,07) 0.147012909874320E-01 0.516692929993857E+01 +(10,08) -0.299783907830715E-01 0.478244480753043E+01 +(10,09) 0.106510026380420E-01 0.190200983103669E+01 +(10,10) 0.321145392954350E-01 0.479668913780619E+01 +(10,11) -0.327010415494442E-01 -0.406712268774822E+01 +(10,12) 0.469501204788685E-01 -0.129331403067363E+02 +(10,13) -0.509594380855560E-01 -0.101318538315566E+02 +(10,14) -0.599967837333679E-01 -0.435269015938419E+01 +(10,15) -0.323531590402126E-01 0.642934255142023E+01 +(10,16) -0.460793916136026E-02 0.455735933529251E+00 +(10,17) -0.224983040243387E-01 -0.346483670488126E+01 +(10,18) 0.588687174022198E-01 0.108965933547002E+02 +(10,19) -0.269168405793607E-02 0.146020551701604E+01 +(10,20) 0.760632846504450E-02 0.935891778893215E+01 +(10,21) -0.926463399082422E-02 -0.109012901159488E+01 +(11,01) -0.162329792510718E-02 -0.439722468868773E+02 +(11,02) -0.155836262274534E-02 -0.536247043860828E+02 +(11,03) -0.135748169850558E-02 -0.303963712026272E+02 +(11,04) -0.150117883458734E-02 -0.528578708468503E+02 +(11,05) -0.450650841230527E-03 -0.306542964087078E+02 +(11,06) 0.917027471587062E-03 0.432088828086063E+02 +(11,07) 0.112307444214821E+00 0.284895856343348E+01 +(11,08) 0.859148427844048E-01 0.612570733710553E+01 +(11,09) 0.108720205724239E+00 0.319395569891390E+01 +(11,10) -0.932282283902168E-01 -0.756710064416565E+01 +(11,11) 0.284656975418329E-01 -0.555260155688967E+01 +(11,12) -0.560129657387733E-01 0.729138222067640E+01 +(11,13) -0.110948169603944E-01 -0.982617231657246E+01 +(11,14) -0.682298317551613E-01 0.839812251519917E+01 +(11,15) 0.883226171135902E-01 -0.269194027969023E+01 +(11,16) -0.131774798035622E+00 -0.226683480127203E+01 +(11,17) -0.403642915189266E-01 -0.641562127865956E+01 +(11,18) 0.920902267098427E-01 -0.680007273581139E+01 +(11,19) -0.128290548920631E+00 -0.197494901501065E+01 +(11,20) -0.114749923348427E+00 -0.173653424919703E+00 +(11,21) -0.100844800472260E+00 0.163116606552273E+01 +(12,01) -0.122732308227569E-03 -0.152518665905767E+03 +(12,02) 0.161253090482205E-03 -0.186570947710702E+03 +(12,03) -0.120117227197625E-03 -0.106601315763174E+03 +(12,04) 0.204894982744008E-03 -0.187054534721548E+03 +(12,05) -0.675647606840357E-04 -0.108517073248987E+03 +(12,06) 0.143623954500072E-03 0.153127156022849E+03 +(12,07) 0.312066711485386E-01 -0.123830418453637E+02 +(12,08) -0.353118106722832E-01 -0.304682312306106E+02 +(12,09) 0.467785857617855E-01 -0.686119827740334E+01 +(12,10) 0.669700726866722E-01 0.285727839324843E+02 +(12,11) 0.159731749445200E-01 0.272695375357948E+02 +(12,12) 0.522553734481335E-01 -0.145329020460164E+02 +(12,13) -0.301160980015993E-01 0.394519469582405E+02 +(12,14) -0.422964245080948E-01 -0.262032227418349E+02 +(12,15) -0.366332642734051E-01 0.440123540181984E+01 +(12,16) 0.137751828879118E-01 -0.480152374806063E+01 +(12,17) -0.348058715462685E-01 0.206844374870974E+02 +(12,18) 0.687717227265239E-02 0.176681540870106E+02 +(12,19) 0.186093822121620E-01 -0.322685179905062E+01 +(12,20) -0.138963768258691E-01 -0.693072145710399E+01 +(12,21) 0.271552260965109E-01 -0.299148354256703E+01 +(13,01) -0.116050858050585E-01 -0.974433821352575E+02 +(13,02) 0.380535630029044E-05 -0.118894186687256E+03 +(13,03) 0.164038110524416E-01 -0.678145724931819E+02 +(13,04) -0.270099262706935E-04 -0.120113847111271E+03 +(13,05) 0.164625588804483E-01 -0.697800218501777E+02 +(13,06) 0.116816693916917E-01 0.982709046226804E+02 +(13,07) 0.839531794190407E-01 0.183794147361240E+01 +(13,08) -0.290049593895674E-01 -0.211184933123850E+01 +(13,09) -0.165680781006813E+00 -0.341701670805002E+01 +(13,10) -0.581433475017548E-01 -0.617525955449609E+00 +(13,11) 0.208601281046867E+00 -0.107290136247085E+01 +(13,12) 0.179246842861176E+00 0.576311903496641E+01 +(13,13) -0.210931494832039E+00 0.348525089218890E+01 +(13,14) 0.144276302307844E-01 -0.182340900435565E+01 +(13,15) 0.231533348560333E+00 -0.196481152900671E+01 +(13,16) -0.211453773081303E-01 0.515398790870046E+01 +(13,17) 0.206030920147896E+00 0.612830413789551E+01 +(13,18) 0.661853179335594E-01 -0.330426561798086E+01 +(13,19) 0.334390364587307E-01 0.266628566506736E+01 +(13,20) 0.252219894900918E-02 -0.631641611962441E+01 +(13,21) 0.526576638221741E-01 -0.833732776936585E+01 +(14,01) -0.437273271381855E-02 -0.359695109817143E+02 +(14,02) -0.462692696601152E-02 -0.441431137577791E+02 +(14,03) -0.348892970941961E-02 -0.256068810893284E+02 +(14,04) -0.462776422500610E-02 -0.442923373719554E+02 +(14,05) -0.185345404315740E-02 -0.255716486578541E+02 +(14,06) 0.320841209031641E-02 0.361682912515179E+02 +(14,07) 0.167358413338661E+00 -0.106016872374744E+01 +(14,08) 0.120894066989422E+00 -0.188238211287414E+01 +(14,09) 0.146612733602524E+00 -0.994102032235108E+00 +(14,10) -0.114104241132736E+00 0.223591971782469E+01 +(14,11) 0.352907255291939E-01 0.159655429090737E+01 +(14,12) -0.118334278464317E+00 -0.205361173619823E+01 +(14,13) -0.409047538414598E-02 0.264356759617923E+01 +(14,14) -0.129458680748940E+00 -0.219743789511244E+01 +(14,15) 0.150047942996025E+00 0.750051299535086E+00 +(14,16) 0.104290992021561E+00 0.595419527031237E+00 +(14,17) -0.437837503850460E-01 0.164954543747536E+01 +(14,18) 0.142482638359070E+00 0.187932831454260E+01 +(14,19) 0.945448875427246E-01 0.525623994704737E+00 +(14,20) 0.911889076232910E-01 0.178475542825685E+00 +(14,21) 0.619743093848228E-01 -0.364538154174523E+00 +(15,01) -0.218217610381544E-02 -0.310917736032053E+01 +(15,02) -0.264662108384073E-02 -0.377923704249770E+01 +(15,03) 0.152970384806395E-02 -0.184422147463024E+01 +(15,04) 0.263672322034836E-02 -0.324061827036610E+01 +(15,05) -0.152716226875782E-02 -0.209529573424705E+01 +(15,06) -0.216715293936431E-02 0.269519062926751E+01 +(15,07) 0.769277736544609E-01 0.508196426964923E+00 +(15,08) 0.140253871679306E+00 0.117105367866777E+01 +(15,09) -0.757992193102837E-01 -0.908058618311773E-01 +(15,10) 0.182921692728996E+00 0.164352591385766E+00 +(15,11) -0.116774775087833E+00 -0.788766372658534E+00 +(15,12) 0.197751015424728E+00 0.116885552810645E+01 +(15,13) 0.239971160888672E+00 -0.215046535598994E+00 +(15,14) 0.236179351806641E+00 0.144881662476397E+01 +(15,15) 0.138329446315765E+00 0.178339724976908E+00 +(15,16) -0.386442616581917E-02 -0.116071749607031E+00 +(15,17) -0.119596652686596E+00 -0.909498295398762E+00 +(15,18) 0.228947326540947E+00 0.610646926717143E-01 +(15,19) 0.306998658925295E-02 -0.153530058087171E+00 +(15,20) -0.562079949304461E-02 0.286291057924642E+00 +(15,21) -0.900598522275686E-02 0.297948743912008E+00 +(16,01) -0.428344011306763E-01 -0.149205064262646E+02 +(16,02) -0.119414471555501E-03 -0.182142441709329E+02 +(16,03) 0.605285130441189E-01 -0.974811729836851E+01 +(16,04) 0.105399143649265E-03 -0.174196586571617E+02 +(16,05) 0.603122040629387E-01 -0.983403637118013E+01 +(16,06) 0.425829663872719E-01 0.148648771225544E+02 +(16,07) 0.147546008229256E+00 -0.335646648007986E+00 +(16,08) 0.401769354939461E-01 0.201470033855699E+01 +(16,09) -0.189309567213058E+00 0.104369831287472E+01 +(16,10) 0.683742910623550E-01 -0.121294114124208E+01 +(16,11) 0.125439286231995E+00 -0.829328960956144E+00 +(16,12) -0.422803349792957E-01 0.689052934581414E+00 +(16,13) -0.458180941641331E-01 -0.383226269716630E+01 +(16,14) -0.296942126005888E-01 0.373098537300082E+01 +(16,15) -0.109751708805561E+00 -0.170586149813142E+01 +(16,16) -0.193516850471497E+00 -0.193030177533155E+01 +(16,17) -0.665694326162338E-01 -0.168714597279047E+01 +(16,18) -0.817999318242073E-01 -0.547045141065854E+00 +(16,19) 0.279022544622421E+00 0.178195066309647E+01 +(16,20) -0.160081461071968E-01 0.313524660082812E+00 +(16,21) 0.335453636944294E-01 0.541301151089991E+00 +(17,01) -0.921052098274231E+00 -0.313781087163451E+02 +(17,02) 0.112806463241577E+01 -0.378487347233482E+02 +(17,03) -0.651386320590973E+00 -0.214243947749121E+02 +(17,04) 0.112832391262054E+01 -0.376461131371173E+02 +(17,05) -0.651349246501923E+00 -0.219732301247688E+02 +(17,06) 0.921293020248413E+00 0.311686483844811E+02 +(17,07) -0.129972562193871E+00 0.126800002733523E+01 +(17,08) 0.182197511196136E+00 0.365548308626798E+01 +(17,09) -0.253000296652317E-02 0.103660756397868E+01 +(17,10) -0.357133313082159E-02 -0.498166767375460E+01 +(17,11) 0.184675022959709E+00 -0.400092670381805E+01 +(17,12) 0.551109127700329E-01 0.531406946549384E+01 +(17,13) 0.231322161853313E-01 -0.514788968725777E+01 +(17,14) -0.134687021374702E+00 0.447991231100838E+01 +(17,15) -0.352683709934354E-02 -0.223730164941427E+01 +(17,16) 0.955944135785103E-01 0.638760991508856E+00 +(17,17) -0.204981788992882E+00 -0.302527057763186E+01 +(17,18) -0.310501866042614E-01 -0.441444840698545E+01 +(17,19) 0.732915475964546E-01 0.213306065587418E+00 +(17,20) -0.121698610484600E+00 -0.520607531935427E+00 +(17,21) 0.487866103649139E-01 0.730688376616387E+00 +(18,01) 0.845294289320009E-05 0.369665247269374E+03 +(18,02) 0.200470458366908E-04 0.451735937919445E+03 +(18,03) -0.500704288697307E-06 0.258293555787007E+03 +(18,04) -0.117598438009736E-04 0.455925831388285E+03 +(18,05) 0.122016963359783E-04 0.265784280467111E+03 +(18,06) 0.177297315531177E-04 -0.372418765344440E+03 +(18,07) 0.166253838688135E-01 -0.225545830162084E+01 +(18,08) -0.616515334695578E-02 0.124672280645443E+02 +(18,09) -0.348300337791443E-01 0.112465777171278E+02 +(18,10) -0.131978942081332E-01 -0.945248700625405E+00 +(18,11) 0.483418852090836E-01 0.332314552507453E+01 +(18,12) 0.501377508044243E-01 -0.139535378216261E+02 +(18,13) -0.548814050853252E-01 -0.195743707104671E+02 +(18,14) -0.183727266266942E-02 0.870092400349449E+01 +(18,15) 0.718188360333443E-01 0.100152523582761E+02 +(18,16) -0.103118065744638E-01 -0.171444008663923E+02 +(18,17) 0.654948502779007E-01 -0.218953533683405E+02 +(18,18) 0.109703140333295E-01 0.713189479375495E+01 +(18,19) 0.177329909056425E-01 -0.119032365190316E+02 +(18,20) 0.323858810588717E-02 0.251995942921836E+02 +(18,21) 0.499825514853001E-01 0.312135303090508E+02 +(19,01) -0.682211592793465E-01 0.765967706327448E+01 +(19,02) -0.430850574048236E-03 0.877594356269220E+01 +(19,03) 0.962435081601143E-01 0.415620530045576E+01 +(19,04) 0.422304205130786E-03 0.850982221207004E+01 +(19,05) 0.957677885890007E-01 0.445399360802617E+01 +(19,06) 0.675818175077438E-01 -0.745326318445295E+01 +(19,07) 0.202333241701126E+00 0.370654267814242E+01 +(19,08) 0.746065378189087E-01 0.561632690325829E+01 +(19,09) -0.223254755139351E+00 0.100203495234038E+01 +(19,10) 0.119545228779316E+00 -0.408091666676030E+01 +(19,11) 0.991457328200340E-01 -0.360166368714016E+01 +(19,12) -0.162345185875893E+00 0.392121124117465E+01 +(19,13) 0.185034181922674E-01 -0.644936945503097E+01 +(19,14) -0.627368614077568E-01 0.453191069282601E+01 +(19,15) -0.280645161867142E+00 -0.207986659450949E+01 +(19,16) 0.505855195224285E-01 0.914629091127449E+00 +(19,17) -0.182546600699425E+00 -0.370670813302335E+01 +(19,18) -0.147651746869087E+00 -0.476882671010381E+01 +(19,19) -0.701333060860634E-01 -0.148427101363476E+01 +(19,20) 0.668072421103716E-02 0.220237913196195E+01 +(19,21) 0.869733467698097E-02 0.117938679293896E+01 +(20,01) -0.449361046776175E-02 0.902376000084666E+00 +(20,02) 0.220308941788971E-02 -0.801193767960690E+00 +(20,03) 0.513966614380479E-02 -0.307272834790146E+01 +(20,04) -0.225867680273950E-02 -0.110857606565975E+01 +(20,05) 0.263794837519526E-02 -0.320007272165156E+01 +(20,06) 0.947628752328455E-03 -0.121510907580541E+01 +(20,07) 0.152274832129478E+00 -0.703746631901387E-01 +(20,08) -0.117567442357540E+00 -0.221727238903005E+01 +(20,09) -0.206277355551720E+00 -0.190035138613005E+01 +(20,10) -0.207551613450050E+00 0.118705706035607E+01 +(20,11) 0.124761372804642E+00 0.218395870949059E+01 +(20,12) 0.863072555512190E-02 -0.225848673886085E+01 +(20,13) -0.714888703078032E-02 0.183569641575782E+01 +(20,14) 0.456247404217720E-01 -0.103488813996630E+01 +(20,15) -0.708496421575546E-01 -0.190399084333949E+00 +(20,16) 0.107937427237630E-01 -0.149848769988273E+00 +(20,17) -0.103561900556087E+00 0.263149954374569E+00 +(20,18) 0.195673510432243E+00 0.226312056314707E+01 +(20,19) -0.138991530984640E-01 -0.215379472713573E+00 +(20,20) -0.738288741558790E-02 -0.717177684354636E+00 +(20,21) 0.119809703901410E-01 0.395015106415659E+00 +(21,01) -0.111714076995850E+01 -0.157006272972466E+01 +(21,02) 0.136832118034363E+01 -0.165695410559176E+01 +(21,03) 0.789960503578186E+00 -0.102464959222303E+01 +(21,04) -0.136850726604462E+01 -0.197490477262258E+01 +(21,05) -0.790138483047485E+00 -0.123963009795365E+01 +(21,06) -0.111746346950531E+01 0.127940329351346E+01 +(21,07) -0.158849164843559E+00 -0.115726516828563E-01 +(21,08) 0.157852247357368E+00 0.278837130191972E+00 +(21,09) 0.170709997415543E+00 0.100161817596371E+00 +(21,10) 0.214689865708351E+00 -0.726523579353106E-01 +(21,11) -0.167143829166889E-01 -0.678547976954667E-01 +(21,12) -0.241510886698961E-01 0.377502144275716E-01 +(21,13) -0.105082295835018E+00 -0.447407742013355E+00 +(21,14) 0.101131219416857E-01 0.419175912148910E+00 +(21,15) 0.525299832224846E-01 -0.233384037060628E+00 +(21,16) -0.143686249852180E+00 -0.418753822037995E-02 +(21,17) 0.132168933749199E+00 -0.856999606537821E-01 +(21,18) -0.693865939974785E-01 -0.714580334063162E-02 +(21,19) 0.104032047092915E+00 -0.161645225161624E-01 +(21,20) 0.185707733035088E+00 -0.856160510487734E-02 +(21,21) -0.296475943177938E-01 0.957054267749156E-01 +#COL_UPDATE_INDEX: 12 +#COL_UPDATE_COMP_(01): 0.102125488221645E+00 +#COL_UPDATE_COMP_(02): 0.962643101811409E-01 +#COL_UPDATE_COMP_(03): -0.180085301399231E-01 +#COL_UPDATE_COMP_(04): -0.116052091121674E+00 +#COL_UPDATE_COMP_(05): 0.100100971758366E-01 +#COL_UPDATE_COMP_(06): 0.168905649334192E-01 +#COL_UPDATE_COMP_(07): -0.168554261326790E-01 +#COL_UPDATE_COMP_(08): 0.334977582097054E-01 +#COL_UPDATE_COMP_(09): 0.775983557105064E-01 +#COL_UPDATE_COMP_(10): -0.509594380855560E-01 +#COL_UPDATE_COMP_(11): -0.110948169603944E-01 +#COL_UPDATE_COMP_(12): -0.301160980015993E-01 +#COL_UPDATE_COMP_(13): -0.210931494832039E+00 +#COL_UPDATE_COMP_(14): -0.409047538414598E-02 +#COL_UPDATE_COMP_(15): 0.239971160888672E+00 +#COL_UPDATE_COMP_(16): -0.458180941641331E-01 +#COL_UPDATE_COMP_(17): 0.231322161853313E-01 +#COL_UPDATE_COMP_(18): -0.548814050853252E-01 +#COL_UPDATE_COMP_(19): 0.185034181922674E-01 +#COL_UPDATE_COMP_(20): -0.714888703078032E-02 +#COL_UPDATE_COMP_(21): -0.105082295835018E+00 +#COL_UPDATE_INDEX: 13 +#COL_UPDATE_COMP_(01): -0.727039668709040E-02 +#COL_UPDATE_COMP_(02): -0.329451821744442E-01 +#COL_UPDATE_COMP_(03): 0.230633586645126E+00 +#COL_UPDATE_COMP_(04): 0.322856456041336E-01 +#COL_UPDATE_COMP_(05): 0.188164815306664E+00 +#COL_UPDATE_COMP_(06): 0.976034477353096E-01 +#COL_UPDATE_COMP_(07): 0.124376058578491E+00 +#COL_UPDATE_COMP_(08): -0.241902604699135E+00 +#COL_UPDATE_COMP_(09): -0.214029267430305E+00 +#COL_UPDATE_COMP_(10): -0.164923388510942E-01 +#COL_UPDATE_COMP_(11): -0.796175077557564E-01 +#COL_UPDATE_COMP_(12): 0.552074126899242E-01 +#COL_UPDATE_COMP_(13): 0.794003605842590E-01 +#COL_UPDATE_COMP_(14): -0.999749153852463E-01 +#COL_UPDATE_COMP_(15): 0.135955372825265E-01 +#COL_UPDATE_COMP_(16): -0.876737236976624E-01 +#COL_UPDATE_COMP_(17): 0.210409909486771E+00 +#COL_UPDATE_COMP_(18): 0.177621953189373E-01 +#COL_UPDATE_COMP_(19): -0.151124328374863E+00 +#COL_UPDATE_COMP_(20): 0.247122272849083E+00 +#COL_UPDATE_COMP_(21): -0.161797225475311E+00 +#COL_UPDATE_INDEX: 21 +#COL_UPDATE_COMP_(01): -0.817385390400887E-01 +#COL_UPDATE_COMP_(02): 0.108162150718272E-02 +#COL_UPDATE_COMP_(03): 0.693925749510527E-02 +#COL_UPDATE_COMP_(04): 0.493063451722264E-02 +#COL_UPDATE_COMP_(05): -0.334013439714909E-01 +#COL_UPDATE_COMP_(06): 0.323923602700233E-01 +#COL_UPDATE_COMP_(07): 0.120032556355000E+00 +#COL_UPDATE_COMP_(08): -0.119067849591374E-01 +#COL_UPDATE_COMP_(09): -0.334841385483742E-01 +#COL_UPDATE_COMP_(10): 0.125436363741755E-01 +#COL_UPDATE_COMP_(11): -0.145120620727539E+00 +#COL_UPDATE_COMP_(12): -0.101416520774364E-01 +#COL_UPDATE_COMP_(13): -0.723919421434402E-01 +#COL_UPDATE_COMP_(14): 0.134314149618149E+00 +#COL_UPDATE_COMP_(15): -0.125767393037677E-01 +#COL_UPDATE_COMP_(16): 0.361425074515864E-03 +#COL_UPDATE_COMP_(17): -0.351854860782623E-01 +#COL_UPDATE_COMP_(18): -0.506495498120785E-01 +#COL_UPDATE_COMP_(19): -0.291761625558138E-01 +#COL_UPDATE_COMP_(20): -0.932136957999319E-03 +#COL_UPDATE_COMP_(21): -0.500183776021004E-01 +#END_PACKET \ No newline at end of file diff --git a/tests/convert-to-h5.py b/tests/convert-to-h5.py new file mode 100644 index 0000000..32068a4 --- /dev/null +++ b/tests/convert-to-h5.py @@ -0,0 +1,46 @@ +import h5py +import numpy as np +from parse import parse + +def rl(rf): + return " ".join(rf.readline().split()) + + +with h5py.File('datasets.hdf5', 'w') as f: + with open('dataset.dat', 'r') as rf: + while(1): + line = rl(rf) + if not line or not line.startswith('#START_PACKET'): + break + cycle_id = parse('#CYCLE_ID: {:d}', rl(rf))[0] + slater_matrix_dim = parse('#SLATER_MATRIX_DIM: {:d}', rl(rf))[0] + nupdates = parse('#NUPDATES: {:d}', rl(rf))[0] + assert(rf.readline().startswith('#SLATER_MATRIX')) + + # Read matrices + slater_matrix = np.zeros((slater_matrix_dim,slater_matrix_dim)) + slater_inverse = np.zeros((slater_matrix_dim,slater_matrix_dim)) + for i in range(slater_matrix_dim*slater_matrix_dim): + res = parse('({i:d},{j:d}) {sla:e} {inv:e}', rl(rf)) + slater_matrix[res['i']-1, res['j']-1] = res['sla'] + slater_inverse[res['i']-1, res['j']-1] = res['inv'] + + # Read updates + col_update_index = np.zeros(nupdates, dtype='i') + updates = np.zeros((nupdates, slater_matrix_dim)) + for n in range(nupdates): + col_update_index[n] = parse('#COL_UPDATE_INDEX: {:d}', rl(rf))[0] + for i in range(slater_matrix_dim): + res = parse('#COL_UPDATE_COMP_({i:d}): {x:e}', rl(rf)) + updates[n][res['i']-1] = res['x'] + + assert(rf.readline().startswith('#END_PACKET')) + rf.readline() + + cycle = f.create_group('cycle_{}'.format(cycle_id)) + cycle.create_dataset("slater_matrix_dim", data=slater_matrix_dim) + cycle.create_dataset("nupdates", data=nupdates) + cycle.create_dataset("slater_matrix", data=slater_matrix, compression='gzip') + cycle.create_dataset("slater_inverse", data=slater_inverse, compression='gzip') + cycle.create_dataset("col_update_index", data=col_update_index) + cycle.create_dataset("updates", data=updates, compression='gzip') diff --git a/tests/test.cpp b/tests/test.cpp new file mode 100644 index 0000000..0e4ded0 --- /dev/null +++ b/tests/test.cpp @@ -0,0 +1,89 @@ +#include +#include +#include "hdf5/serial/hdf5.h" +#include "H5Cpp.h" + +#include "../SM_MaponiA3.hpp" +#include "../Helpers.hpp" + +using namespace H5; +#define DEBUG 1 + +const H5std_string FILE_NAME( "datasets.hdf5" ); + +void read_int(H5File file, std::string key, unsigned int * data) { + DataSet ds = file.openDataSet(key); + ds.read(data, PredType::STD_U32LE); + ds.close(); +} + +void read_double(H5File file, std::string key, double * data) { + DataSet ds = file.openDataSet(key); + ds.read(data, PredType::IEEE_F64LE); + ds.close(); +} + + +int test_cycle(H5File file, int cycle) { + + /* Read the data */ + + std::string group = "cycle_" + std::to_string(cycle); + + unsigned int dim, nupdates; + read_int(file, group + "/slater_matrix_dim", &dim); + read_int(file, group + "/nupdates", &nupdates); + + double * slater_matrix = new double[dim*dim]; + read_double(file, group + "/slater_matrix", slater_matrix); + + double * slater_inverse = new double[dim*dim]; + read_double(file, group + "/slater_inverse", slater_inverse); + + unsigned int * col_update_index = new unsigned int[nupdates]; + read_int(file, group + "/col_update_index", col_update_index); + + double * updates = new double[nupdates*dim]; + read_double(file, group + "/updates", updates); + + /* Test */ +#ifdef DEBUG + showMatrix(slater_matrix, dim, "Slater"); +#endif + + MaponiA3(slater_matrix, slater_inverse, dim, nupdates, updates, col_update_index); + +#ifdef DEBUG + showMatrix(slater_inverse, dim, "Inverse"); +#endif + + double * res = matMul(slater_matrix, slater_inverse, dim); + bool ok = is_identity(res, dim, 1.0e-8); + +#ifdef DEBUG + showMatrix(res, dim, "Result"); +#endif + + delete [] res, updates, col_update_index, slater_matrix, slater_inverse; + + return ok; +} + +int main(int argc, char **argv) { + if (argc != 2) { + std::cerr << "usage: ./test_dataset " << std::endl; + return 1; + } + int cycle = std::stoi(argv[1]); + H5File file(FILE_NAME, H5F_ACC_RDONLY); + + bool ok = test_cycle(file, 21); + + if (ok) { + std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl; + } else { + std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl; + } + return ok; +} +