diff --git a/make.config.gfortran b/make.config.gfortran index ac6352b..6e5b0e6 100644 --- a/make.config.gfortran +++ b/make.config.gfortran @@ -7,7 +7,7 @@ CPU_TYPE="-mavx" ## FORTRAN compiler FC="gfortran -ffree-line-length-none" FCFLAGS="-O2 -g ${CPU_TYPE}" -LIB="-lblas -llapack -lpthread" +LIB="-lblas -llapack -lpthread -lqmckl" ## IRPF90 IRPF90="${QMCCHEM_PATH}/bin/irpf90" diff --git a/src/det.irp.f b/src/det.irp.f index 71488b2..6218438 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -20,7 +20,8 @@ END_PROVIDER END_PROVIDER -subroutine det_update(n,LDS,m,l,S,S_inv,d) + subroutine det_update(n,LDS,m,l,S,S_inv,d) + use qmckl implicit none integer, intent(in) :: n,LDS ! Dimension of the vector @@ -30,169 +31,245 @@ subroutine det_update(n,LDS,m,l,S,S_inv,d) 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 :: zl, lambda, d_inv 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 + +!! -- START QMCKL/SMW ADDITIONS -- !! +!! !! + integer :: i, j + integer (qmckl_exit_code) :: rc + integer (qmckl_context) :: context + integer(kind=8) :: ddim, nupdates, updates_index(1) + real(c_double) :: updates(n,1), breakdown + double precision :: S_inv_sq(n,n) + + S_inv = S_inv / d + + ! open(unit = 2000, file = "Slater_old.dat") + ! open(unit = 3000, file = "Slater_old_inv.dat") + ! do i=1,n + ! do j=1,n + ! write(2000,"(E23.15, 1X)", advance="no") S(j,i) ! write transpose for Octave + ! write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) + ! end do + ! write(2000,*) + ! write(3000,*) + ! end do + ! flush(2000) + ! flush(3000) + ! close(2000) + ! close(3000) + + do i=1,n + updates(i,1) = m(i) - S(i,l) ! transform repl. upds. into additive upds. + S(i,l) = m(i) ! update S with repl. upds + end do + + zl = 0 + do i=1,n + zl = zl + S_inv(i,l) * updates(i,1) + end do + + d_inv = 1.d0/d + d = d + zl + lambda = d * d_inv + if ( dabs(lambda) < 1.d-3 ) then + d = 0.d0 + return + endif + + context = qmckl_context_create() + ddim = n + nupdates = 1 + updates_index(1) = l + breakdown = 1e-3 + S_inv_sq = S_inv(1:n,1:n) + rc = qmckl_sherman_morrison_splitting(context, ddim, nupdates, updates, updates_index, breakdown, S_inv_sq) + S_inv = S_inv_sq(1:n,1:n) + rc = qmckl_context_destroy(context) + + ! open(unit = 4000, file = "Slater.dat") + ! open(unit = 5000, file = "Slater_inv.dat") + ! do i=1,n + ! do j=1,n + ! write(4000,"(E23.15, 1X)", advance="no") S(j,i) ! write transpose for Octave + ! write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j) + ! end do + ! write(4000,*) + ! write(5000,*) + ! end do + ! flush(4000) + ! flush(5000) + ! close(4000) + ! close(5000) + + S_inv = S_inv * d + + ! 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 QMCKL/SMW ADDITIONS -- !! end subroutine det_update2(n,LDS,m,l,S,S_inv,d) @@ -430,7 +507,7 @@ SUBST [ n ] END_TEMPLATE -BEGIN_TEMPLATE +BEGIN_TEMPLATE ! Version for mod(n,4) = 1 subroutine det_update$n(n,LDS,m,l,S,S_inv,d) implicit none @@ -1186,7 +1263,7 @@ END_PROVIDER BEGIN_DOC ! det_alpha_value_curr : Value of the current alpha determinant ! - ! det_alpha_value_curr : Slater matrix for 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 @@ -1200,6 +1277,7 @@ END_PROVIDER double precision :: tmp_inv(elec_alpha_num_8) real :: tmp_det(elec_alpha_num_8) integer, save :: ifirst + logical :: file_exists !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: tmp_inv, tmp_det if (ifirst == 0) then @@ -1211,12 +1289,12 @@ END_PROVIDER endif PROVIDE mo_value - if (det_i /= det_alpha_order(1) ) then + if (det_i /= det_alpha_order(1) ) then ! alpha determinant order changes n_to_do = 0 do k=1,elec_alpha_num imo = mo_list_alpha_curr(k) - if ( imo /= mo_list_alpha_prev(k) ) then + if ( imo /= mo_list_alpha_prev(k) ) then ! mo for electron k has changed n_to_do += 1 to_do(n_to_do) = k endif @@ -1288,14 +1366,14 @@ END_PROVIDER ddet = 0.d0 if (n_to_do < shiftl(elec_alpha_num,1)) then - do while ( n_to_do > 0 ) - ddet = det_alpha_value_curr - n_to_do_old = n_to_do + ddet = det_alpha_value_curr ! remember value of det_alpha_value_curr + n_to_do_old = n_to_do ! remember n_to_do value n_to_do = 0 do l=1,n_to_do_old - k = to_do(l) - imo = mo_list_alpha_curr(k) + k = to_do(l) ! select electron to change + imo = mo_list_alpha_curr(k) ! select mo to change + ! write(*,*) "k, imo, mo_value(1,imo) = ", k, imo, mo_value(1,imo) call det_update(elec_alpha_num, elec_alpha_num_8, & mo_value(1,imo), & k, & @@ -1344,12 +1422,13 @@ END_PROVIDER slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k)) enddo enddo + ! write(*,*) "FIRST TIME OR ALL FAILED; DO LAPACK" 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 ] @@ -1467,6 +1546,7 @@ END_PROVIDER ddet = det_beta_value_curr n_to_do_old = n_to_do n_to_do = 0 + loopcount = 0 do l=1,n_to_do_old k = to_do(l) imo = mo_list_beta_curr(k) diff --git a/src/qmckl b/src/qmckl index de03986..8535e3c 160000 --- a/src/qmckl +++ b/src/qmckl @@ -1 +1 @@ -Subproject commit de03986bda2be207377875ed5a0852cb721b86b9 +Subproject commit 8535e3c0b5a4f4194ca0400b348774c3cc7c74d0