diff --git a/data/basis/cc-pcvdz b/data/basis/cc-pcvdz index d874fb06..76985d4a 100644 --- a/data/basis/cc-pcvdz +++ b/data/basis/cc-pcvdz @@ -991,4 +991,266 @@ D 1 1 1.3743000 1.0000000 D 1 1 0.0537000 1.00000000 -$END \ No newline at end of file + +COPPER +S 20 +1 5.430321E+06 7.801026E-06 +2 8.131665E+05 6.065666E-05 +3 1.850544E+05 3.188964E-04 +4 5.241466E+04 1.344687E-03 +5 1.709868E+04 4.869050E-03 +6 6.171994E+03 1.561013E-02 +7 2.406481E+03 4.452077E-02 +8 9.972584E+02 1.103111E-01 +9 4.339289E+02 2.220342E-01 +10 1.962869E+02 3.133739E-01 +11 9.104280E+01 2.315121E-01 +12 4.138425E+01 7.640920E-02 +13 1.993278E+01 1.103818E-01 +14 9.581891E+00 1.094372E-01 +15 4.234516E+00 1.836311E-02 +16 1.985814E+00 -6.043084E-04 +17 8.670830E-01 5.092245E-05 +18 1.813390E-01 -5.540730E-05 +19 8.365700E-02 3.969482E-05 +20 3.626700E-02 -1.269538E-05 +S 20 +1 5.430321E+06 -4.404706E-06 +2 8.131665E+05 -3.424801E-05 +3 1.850544E+05 -1.801238E-04 +4 5.241466E+04 -7.600455E-04 +5 1.709868E+04 -2.759348E-03 +6 6.171994E+03 -8.900970E-03 +7 2.406481E+03 -2.579378E-02 +8 9.972584E+02 -6.623861E-02 +9 4.339289E+02 -1.445927E-01 +10 1.962869E+02 -2.440110E-01 +11 9.104280E+01 -2.504837E-01 +12 4.138425E+01 2.852577E-02 +13 1.993278E+01 5.115874E-01 +14 9.581891E+00 4.928061E-01 +15 4.234516E+00 8.788437E-02 +16 1.985814E+00 -5.820281E-03 +17 8.670830E-01 2.013508E-04 +18 1.813390E-01 -5.182553E-04 +19 8.365700E-02 3.731503E-04 +20 3.626700E-02 -1.193171E-04 +S 20 +1 5.430321E+06 9.704682E-07 +2 8.131665E+05 7.549245E-06 +3 1.850544E+05 3.968892E-05 +4 5.241466E+04 1.677200E-04 +5 1.709868E+04 6.095101E-04 +6 6.171994E+03 1.978846E-03 +7 2.406481E+03 5.798049E-03 +8 9.972584E+02 1.534158E-02 +9 4.339289E+02 3.540484E-02 +10 1.962869E+02 6.702098E-02 +11 9.104280E+01 8.026945E-02 +12 4.138425E+01 -1.927231E-02 +13 1.993278E+01 -3.160129E-01 +14 9.581891E+00 -4.573162E-01 +15 4.234516E+00 1.550841E-01 +16 1.985814E+00 7.202872E-01 +17 8.670830E-01 3.885122E-01 +18 1.813390E-01 1.924326E-02 +19 8.365700E-02 -7.103807E-03 +20 3.626700E-02 3.272906E-03 +S 20 +1 5.430321E+06 -1.959354E-07 +2 8.131665E+05 -1.523472E-06 +3 1.850544E+05 -8.014808E-06 +4 5.241466E+04 -3.383992E-05 +5 1.709868E+04 -1.231191E-04 +6 6.171994E+03 -3.992085E-04 +7 2.406481E+03 -1.171900E-03 +8 9.972584E+02 -3.096141E-03 +9 4.339289E+02 -7.171993E-03 +10 1.962869E+02 -1.356621E-02 +11 9.104280E+01 -1.643989E-02 +12 4.138425E+01 4.107628E-03 +13 1.993278E+01 6.693964E-02 +14 9.581891E+00 1.028221E-01 +15 4.234516E+00 -4.422945E-02 +16 1.985814E+00 -2.031191E-01 +17 8.670830E-01 -2.230022E-01 +18 1.813390E-01 2.517975E-01 +19 8.365700E-02 5.650091E-01 +20 3.626700E-02 3.247243E-01 +S 20 +1 5.430321E+06 -7.508267E-07 +2 8.131665E+05 -5.972018E-06 +3 1.850544E+05 -3.039682E-05 +4 5.241466E+04 -1.340405E-04 +5 1.709868E+04 -4.615778E-04 +6 6.171994E+03 -1.601064E-03 +7 2.406481E+03 -4.330942E-03 +8 9.972584E+02 -1.265434E-02 +9 4.339289E+02 -2.586864E-02 +10 1.962869E+02 -5.835428E-02 +11 9.104280E+01 -5.132322E-02 +12 4.138425E+01 -1.908953E-02 +13 1.993278E+01 3.586116E-01 +14 9.581891E+00 3.885818E-01 +15 4.234516E+00 -3.057106E-01 +16 1.985814E+00 -2.069896E+00 +17 8.670830E-01 2.431774E+00 +18 1.813390E-01 -2.121974E-02 +19 8.365700E-02 -1.820251E+00 +20 3.626700E-02 1.434585E+00 +S 20 +1 5.430321E+06 -3.532229E-07 +2 8.131665E+05 -2.798812E-06 +3 1.850544E+05 -1.432517E-05 +4 5.241466E+04 -6.270946E-05 +5 1.709868E+04 -2.179490E-04 +6 6.171994E+03 -7.474316E-04 +7 2.406481E+03 -2.049271E-03 +8 9.972584E+02 -5.885203E-03 +9 4.339289E+02 -1.226885E-02 +10 1.962869E+02 -2.683147E-02 +11 9.104280E+01 -2.479261E-02 +12 4.138425E+01 -5.984746E-03 +13 1.993278E+01 1.557124E-01 +14 9.581891E+00 1.436683E-01 +15 4.234516E+00 8.374103E-03 +16 1.985814E+00 -7.460711E-01 +17 8.670830E-01 1.244367E-01 +18 1.813390E-01 1.510110E+00 +19 8.365700E-02 -3.477122E-01 +20 3.626700E-02 -9.774169E-01 +S 1 +1 3.626700E-02 1.000000E+00 +S 1 +1 0.0157200 1.0000000 +P 16 +1 2.276057E+04 4.000000E-05 +2 5.387679E+03 3.610000E-04 +3 1.749945E+03 2.083000E-03 +4 6.696653E+02 9.197000E-03 +5 2.841948E+02 3.266000E-02 +6 1.296077E+02 9.379500E-02 +7 6.225415E+01 2.082740E-01 +8 3.092964E+01 3.339930E-01 +9 1.575827E+01 3.324930E-01 +10 8.094211E+00 1.547280E-01 +11 4.046921E+00 2.127100E-02 +12 1.967869E+00 -1.690000E-03 +13 9.252950E-01 -1.516000E-03 +14 3.529920E-01 -2.420000E-04 +15 1.273070E-01 2.300000E-05 +16 4.435600E-02 -9.000000E-06 +P 16 +1 2.276057E+04 -1.500000E-05 +2 5.387679E+03 -1.310000E-04 +3 1.749945E+03 -7.550000E-04 +4 6.696653E+02 -3.359000E-03 +5 2.841948E+02 -1.208100E-02 +6 1.296077E+02 -3.570300E-02 +7 6.225415E+01 -8.250200E-02 +8 3.092964E+01 -1.398900E-01 +9 1.575827E+01 -1.407290E-01 +10 8.094211E+00 3.876600E-02 +11 4.046921E+00 3.426950E-01 +12 1.967869E+00 4.523100E-01 +13 9.252950E-01 2.770540E-01 +14 3.529920E-01 4.388500E-02 +15 1.273070E-01 -2.802000E-03 +16 4.435600E-02 1.152000E-03 +P 16 +1 2.276057E+04 5.000000E-06 +2 5.387679E+03 4.900000E-05 +3 1.749945E+03 2.780000E-04 +4 6.696653E+02 1.253000E-03 +5 2.841948E+02 4.447000E-03 +6 1.296077E+02 1.337000E-02 +7 6.225415E+01 3.046900E-02 +8 3.092964E+01 5.344700E-02 +9 1.575827E+01 5.263900E-02 +10 8.094211E+00 -1.688100E-02 +11 4.046921E+00 -1.794480E-01 +12 1.967869E+00 -2.095880E-01 +13 9.252950E-01 -3.963300E-02 +14 3.529920E-01 5.021300E-01 +15 1.273070E-01 5.811110E-01 +16 4.435600E-02 4.566600E-02 +P 16 +1 2.276057E+04 1.100000E-05 +2 5.387679E+03 9.600000E-05 +3 1.749945E+03 5.900000E-04 +4 6.696653E+02 2.484000E-03 +5 2.841948E+02 9.463000E-03 +6 1.296077E+02 2.645300E-02 +7 6.225415E+01 6.568900E-02 +8 3.092964E+01 1.027320E-01 +9 1.575827E+01 1.370410E-01 +10 8.094211E+00 -7.096100E-02 +11 4.046921E+00 -5.047080E-01 +12 1.967869E+00 -4.780560E-01 +13 9.252950E-01 9.428920E-01 +14 3.529920E-01 5.446990E-01 +15 1.273070E-01 -8.327660E-01 +16 4.435600E-02 -1.084160E-01 +P 16 +1 2.276057E+04 3.000000E-06 +2 5.387679E+03 2.500000E-05 +3 1.749945E+03 1.470000E-04 +4 6.696653E+02 6.560000E-04 +5 2.841948E+02 2.351000E-03 +6 1.296077E+02 7.004000E-03 +7 6.225415E+01 1.613100E-02 +8 3.092964E+01 2.777000E-02 +9 1.575827E+01 2.756700E-02 +10 8.094211E+00 -1.011500E-02 +11 4.046921E+00 -8.100900E-02 +12 1.967869E+00 -1.104090E-01 +13 9.252950E-01 -7.173200E-02 +14 3.529920E-01 1.879300E-01 +15 1.273070E-01 5.646290E-01 +16 4.435600E-02 4.070000E-01 +P 1 +1 4.435600E-02 1.000000E+00 +P 1 +1 0.0154500 1.0000000 +D 8 +1 1.738970E+02 2.700000E-03 +2 5.188690E+01 2.090900E-02 +3 1.934190E+01 8.440800E-02 +4 7.975720E+00 2.139990E-01 +5 3.398230E+00 3.359800E-01 +6 1.409320E+00 3.573010E-01 +7 5.488580E-01 2.645780E-01 +8 1.901990E-01 1.039720E-01 +D 8 +1 1.738970E+02 -3.363000E-03 +2 5.188690E+01 -2.607900E-02 +3 1.934190E+01 -1.082310E-01 +4 7.975720E+00 -2.822170E-01 +5 3.398230E+00 -3.471900E-01 +6 1.409320E+00 2.671100E-02 +7 5.488580E-01 4.920470E-01 +8 1.901990E-01 4.384220E-01 +D 8 +1 1.738970E+02 4.133000E-03 +2 5.188690E+01 3.308500E-02 +3 1.934190E+01 1.383360E-01 +4 7.975720E+00 3.901660E-01 +5 3.398230E+00 1.698420E-01 +6 1.409320E+00 -6.830180E-01 +7 5.488580E-01 -2.657970E-01 +8 1.901990E-01 8.380630E-01 +D 1 +1 1.901990E-01 1.000000E+00 +D 1 +1 0.0659100 1.0000000 +F 1 +1 5.082100E+00 1.000000E+00 +F 1 +1 1.279700E+00 1.000000E+00 +F 1 +1 0.4617200 1.0000000 +G 1 +1 3.483500E+00 1.0000000 +G 1 +1 1.4597900 1.0000000 +$END diff --git a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f index ad215b41..d2115d9e 100644 --- a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f @@ -224,7 +224,7 @@ subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_ double precision, allocatable :: analytical_j(:) resv(:) = 0.d0 - if(ao_overlap_abs(j,i).lt.1.d-12) then + if(ao_overlap_abs(j,i) .lt. 1.d-12) then return endif @@ -360,9 +360,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta, ASSERT(beta .gt. 0.d0) if(beta .lt. 1d-10) then - call overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points) - return endif @@ -379,19 +377,20 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta, A1_center(1:3) = nucl_coord(ao_nucl(i),1:3) A2_center(1:3) = nucl_coord(ao_nucl(j),1:3) - allocate (fact_g(n_points), G_center(n_points,3), analytical_j(n_points) ) + allocate(fact_g(n_points), G_center(n_points,3), analytical_j(n_points)) bg = beta * gama_inv dg = delta * gama_inv bdg = bg * delta - do ipoint=1,n_points + + do ipoint = 1, n_points + G_center(ipoint,1) = bg * B_center(1) + dg * D_center(ipoint,1) G_center(ipoint,2) = bg * B_center(2) + dg * D_center(ipoint,2) G_center(ipoint,3) = bg * B_center(3) + dg * D_center(ipoint,3) - fact_g(ipoint) = bdg * ( & - (B_center(1) - D_center(ipoint,1)) * (B_center(1) - D_center(ipoint,1)) & - + (B_center(2) - D_center(ipoint,2)) * (B_center(2) - D_center(ipoint,2)) & - + (B_center(3) - D_center(ipoint,3)) * (B_center(3) - D_center(ipoint,3)) ) + fact_g(ipoint) = bdg * ( (B_center(1) - D_center(ipoint,1)) * (B_center(1) - D_center(ipoint,1)) & + + (B_center(2) - D_center(ipoint,2)) * (B_center(2) - D_center(ipoint,2)) & + + (B_center(3) - D_center(ipoint,3)) * (B_center(3) - D_center(ipoint,3)) ) if(fact_g(ipoint) < 10d0) then fact_g(ipoint) = dexp(-fact_g(ipoint)) @@ -415,8 +414,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta, do ipoint = 1, n_points coef12f = coef12 * fact_g(ipoint) resv(ipoint) += coef12f * analytical_j(ipoint) - end do - + enddo enddo enddo diff --git a/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f b/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f index dcd1db66..54c2d95b 100644 --- a/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f +++ b/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f @@ -1,5 +1,9 @@ -double precision function overlap_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta) +! --- + +double precision function overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta) + BEGIN_DOC + ! ! Computes the following integral : ! ! .. math :: @@ -8,23 +12,25 @@ double precision function overlap_gauss_r12(D_center,delta,A_center,B_center,pow ! END_DOC - implicit none include 'constants.include.F' - double precision, intent(in) :: D_center(3), delta ! pure gaussian "D" - double precision, intent(in) :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B" - integer, intent(in) :: power_A(3),power_B(3) - double precision :: overlap_x,overlap_y,overlap_z,overlap + implicit none + double precision, intent(in) :: D_center(3), delta ! pure gaussian "D" + double precision, intent(in) :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B" + integer, intent(in) :: power_A(3),power_B(3) + + double precision :: overlap_x,overlap_y,overlap_z,overlap ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 ) - double precision :: A_new(0:max_dim,3)! new polynom - double precision :: A_center_new(3) ! new center - integer :: iorder_a_new(3) ! i_order(i) = order of the new polynom ==> should be equal to power_A - double precision :: alpha_new ! new exponent - double precision :: fact_a_new ! constant factor - double precision :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr - integer :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1 - dim1=100 - thr = 1.d-10 + double precision :: A_new(0:max_dim,3)! new polynom + double precision :: A_center_new(3) ! new center + integer :: iorder_a_new(3) ! i_order(i) = order of the new polynom ==> should be equal to power_A + double precision :: alpha_new ! new exponent + double precision :: fact_a_new ! constant factor + double precision :: accu, coefx, coefy, coefz, coefxy, coefxyz, thr + integer :: d(3), i, lx, ly, lz, iorder_tmp(3), dim1 + + dim1 = 100 + thr = 1.d-10 d(:) = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 ) == 0 overlap_gauss_r12 = 0.d0 @@ -38,17 +44,22 @@ double precision function overlap_gauss_r12(D_center,delta,A_center,B_center,pow coefx = A_new(lx,1)*fact_a_new if(dabs(coefx).lt.thr)cycle iorder_tmp(1) = lx + do ly = 0, iorder_a_new(2) - coefy = A_new(ly,2) + coefy = A_new(ly,2) coefxy = coefx * coefy - if(dabs(coefxy).lt.thr)cycle + if(dabs(coefxy) .lt. thr) cycle iorder_tmp(2) = ly + do lz = 0, iorder_a_new(3) - coefz = A_new(lz,3) + coefz = A_new(lz,3) coefxyz = coefxy * coefz - if(dabs(coefxyz).lt.thr)cycle + if(dabs(coefxyz) .lt. thr) cycle iorder_tmp(3) = lz - call overlap_gaussian_xyz(A_center_new,B_center,alpha_new,beta,iorder_tmp,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + call overlap_gaussian_xyz( A_center_new, B_center, alpha_new, beta, iorder_tmp, power_B & + , overlap_x, overlap_y, overlap_z, overlap, dim1) + accu += coefxyz * overlap enddo enddo @@ -159,11 +170,9 @@ subroutine overlap_gauss_r12_v(D_center, LD_D, delta, A_center, B_center, power_ maxab = maxval(power_A(1:3)) - allocate(A_new(n_points, 0:maxab, 3), A_center_new(n_points, 3), fact_a_new(n_points), iorder_a_new(3), overlap(n_points)) + allocate(A_new(n_points,0:maxab,3), A_center_new(n_points,3), fact_a_new(n_points), iorder_a_new(3), overlap(n_points)) - call give_explicit_poly_and_gaussian_v(A_new, maxab, A_center_new, & - alpha_new, fact_a_new, iorder_a_new, delta, alpha, d, power_A, & - D_center, LD_D, A_center, n_points) + call give_explicit_poly_and_gaussian_v(A_new, maxab, A_center_new, alpha_new, fact_a_new, iorder_a_new, delta, alpha, d, power_A, D_center, LD_D, A_center, n_points) rvec(:) = 0.d0 diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index 746593dc..c077dea1 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -182,6 +182,27 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3 enddo END_PROVIDER +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)] + + implicit none + integer :: i, j, ipoint + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + int2_grad1_u12_ao_t(ipoint,1,j,i) = int2_grad1_u12_ao(1,j,i,ipoint) + int2_grad1_u12_ao_t(ipoint,2,j,i) = int2_grad1_u12_ao(2,j,i,ipoint) + int2_grad1_u12_ao_t(ipoint,3,j,i) = int2_grad1_u12_ao(3,j,i,ipoint) + enddo + enddo + enddo + +END_PROVIDER + +! --- + BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, mo_num, mo_num, n_points_final_grid)] BEGIN_DOC diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f index c1c27f06..48fa84f7 100644 --- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -15,7 +15,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n character*(128) :: name_file three_body_ints_bi_ort = 0.d0 - print*,'Providing the three_body_ints_bi_ort ...' + print *, ' Providing the three_body_ints_bi_ort ...' call wall_time(wall0) name_file = 'six_index_tensor' @@ -71,7 +71,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -104,12 +104,11 @@ end subroutine give_integrals_3_body_bi_ort ! --- - subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -170,3 +169,39 @@ end subroutine give_integrals_3_body_bi_ort_old ! --- +subroutine give_integrals_3_body_bi_ort_ao(n, l, k, m, j, i, integral) + + BEGIN_DOC + ! + ! < n l k | -L | m j i > with a BI-ORTHONORMAL ATOMIC ORBITALS + ! + END_DOC + + implicit none + integer, intent(in) :: n, l, k, m, j, i + double precision, intent(out) :: integral + integer :: ipoint + double precision :: weight + + integral = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + + integral += weight * aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i) & + * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,l,j) & + + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,l,j) & + + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,l,j) ) + integral += weight * aos_in_r_array_transp(ipoint,l) * aos_in_r_array_transp(ipoint,j) & + * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,k,i) & + + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,k,i) & + + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,k,i) ) + integral += weight * aos_in_r_array_transp(ipoint,n) * aos_in_r_array_transp(ipoint,m) & + * ( int2_grad1_u12_ao_t(ipoint,1,l,j) * int2_grad1_u12_ao_t(ipoint,1,k,i) & + + int2_grad1_u12_ao_t(ipoint,2,l,j) * int2_grad1_u12_ao_t(ipoint,2,k,i) & + + int2_grad1_u12_ao_t(ipoint,3,l,j) * int2_grad1_u12_ao_t(ipoint,3,k,i) ) + + enddo + +end subroutine give_integrals_3_body_bi_ort_ao + +! --- diff --git a/src/bi_ortho_mos/bi_density.irp.f b/src/bi_ortho_mos/bi_density.irp.f index 0de8ce69..56f44da1 100644 --- a/src/bi_ortho_mos/bi_density.irp.f +++ b/src/bi_ortho_mos/bi_density.irp.f @@ -10,6 +10,7 @@ BEGIN_PROVIDER [double precision, TCSCF_bi_ort_dm_ao_alpha, (ao_num, ao_num) ] END_DOC call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0 & , mo_l_coef, size(mo_l_coef, 1), mo_r_coef, size(mo_r_coef, 1) & + !, mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & , 0.d0, TCSCF_bi_ort_dm_ao_alpha, size(TCSCF_bi_ort_dm_ao_alpha, 1) ) END_PROVIDER @@ -24,6 +25,7 @@ BEGIN_PROVIDER [ double precision, TCSCF_bi_ort_dm_ao_beta, (ao_num, ao_num) ] END_DOC call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & , mo_l_coef, size(mo_l_coef, 1), mo_r_coef, size(mo_r_coef, 1) & + !, mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & , 0.d0, TCSCF_bi_ort_dm_ao_beta, size(TCSCF_bi_ort_dm_ao_beta, 1) ) END_PROVIDER diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index d7d8fa7d..cb698fbb 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -1,12 +1,27 @@ +! --- BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta , (ao_num, ao_num) ] - use map_module - implicit none +&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta , (ao_num, ao_num) ] + BEGIN_DOC - ! Alpha and Beta Fock matrices in AO basis set + ! + ! 2-e part of alpha and beta Fock matrices (F^{a} & F^{b}) in AO basis set + ! + ! F^{a} = h + G^{a} + ! F^{b} = h + G^{b} + ! + ! where : + ! F^{a} = J^{a} + J^{b} - K^{a} ==> G_{ij}^{a} = \sum_{k,l} P_{kl} (kl|ij) - P_{kl}^{a} (ki|lj) + ! F^{b} = J^{a} + J^{b} - K^{b} ==> G_{ij}^{b} = \sum_{k,l} P_{kl} (kl|ij) - P_{kl}^{b} (ki|lj) + ! + ! and P_{kl} = P_{kl}^{a} + P_{kl}^{b} + ! END_DOC + use map_module + + implicit none + integer :: i,j,k,l,k1,r,s integer :: i0,j0,k0,l0 integer*8 :: p,q @@ -153,6 +168,8 @@ END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] &BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] implicit none diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 3226073d..f4123c85 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -68,20 +68,33 @@ subroutine create_guess endif end -subroutine run +! --- + +subroutine run() BEGIN_DOC -! Run SCF calculation + ! Run SCF calculation END_DOC use bitmasks implicit none - integer :: i_it, i, j, k - mo_label = 'Orthonormalized' - call Roothaan_Hall_SCF + PROVIDE scf_algorithm + + if(scf_algorithm .eq. "DIIS_MO") then + call Roothaan_Hall_SCF_MO() + elseif(scf_algorithm .eq. "DIIS_MODIF") then + call Roothaan_Hall_SCF_MODIF() + elseif(scf_algorithm .eq. "DIIS") then + call Roothaan_Hall_SCF() + elseif(scf_algorithm .eq. "Simple") then + call Roothaan_Hall_SCF_Simple() + else + print *, ' not implemented yet:', scf_algorithm + endif + call ezfio_set_hartree_fock_energy(SCF_energy) end diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index bb585f63..ca00b816 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -17,7 +17,7 @@ program debug_integ_jmu_modif PROVIDE mu_erf j1b_pen - call test_v_ij_u_cst_mu_j1b() +! call test_v_ij_u_cst_mu_j1b() ! call test_v_ij_erf_rk_cst_mu_j1b() ! call test_x_v_ij_erf_rk_cst_mu_j1b() ! call test_int2_u2_j1b2() @@ -31,6 +31,9 @@ program debug_integ_jmu_modif ! call test_u12_grad1_u12_j1b_grad1_j1b() ! !call test_gradu_squared_u_ij_mu() + !call test_vect_overlap_gauss_r12_ao() + call test_vect_overlap_gauss_r12_ao_with1s() + end ! --- @@ -595,7 +598,183 @@ subroutine test_u12_grad1_u12_j1b_grad1_j1b() print*, ' normalz = ', normalz return -end subroutine test_u12_grad1_u12_j1b_grad1_j1b, +end subroutine test_u12_grad1_u12_j1b_grad1_j1b ! --- +subroutine test_vect_overlap_gauss_r12_ao() + + implicit none + + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: expo_fit, r(3) + double precision, allocatable :: I_vec(:,:,:), I_ref(:,:,:), int_fit_v(:) + + double precision, external :: overlap_gauss_r12_ao + + print *, ' test_vect_overlap_gauss_r12_ao ...' + + provide mu_erf final_grid_points_transp j1b_pen + + expo_fit = expo_gauss_j_mu_x_2(1) + + ! --- + + allocate(int_fit_v(n_points_final_grid)) + allocate(I_vec(ao_num,ao_num,n_points_final_grid)) + + I_vec = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + + call overlap_gauss_r12_ao_v(final_grid_points_transp, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, n_points_final_grid) + + do ipoint = 1, n_points_final_grid + I_vec(j,i,ipoint) = int_fit_v(ipoint) + enddo + enddo + enddo + + ! --- + + allocate(I_ref(ao_num,ao_num,n_points_final_grid)) + + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = 1, ao_num + + I_ref(j,i,ipoint) = overlap_gauss_r12_ao(r, expo_fit, i, j) + enddo + enddo + enddo + + ! --- + + eps_ij = 1d-3 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + + i_exc = I_ref(i,j,ipoint) + i_num = I_vec(i,j,ipoint) + acc_ij = dabs(i_exc - i_num) + !acc_ij = dabs(i_exc - i_num) / dabs(i_exc) + if(acc_ij .gt. eps_ij) then + print *, ' problem in overlap_gauss_r12_ao_v on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + stop + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + enddo + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_vect_overlap_gauss_r12_ao + +! --- + +subroutine test_vect_overlap_gauss_r12_ao_with1s() + + implicit none + + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: expo_fit, r(3), beta, B_center(3) + double precision, allocatable :: I_vec(:,:,:), I_ref(:,:,:), int_fit_v(:) + + double precision, external :: overlap_gauss_r12_ao_with1s + + print *, ' test_vect_overlap_gauss_r12_ao_with1s ...' + + provide mu_erf final_grid_points_transp j1b_pen + + expo_fit = expo_gauss_j_mu_x_2(1) + beta = List_all_comb_b3_expo (2) + B_center(1) = List_all_comb_b3_cent(1,2) + B_center(2) = List_all_comb_b3_cent(2,2) + B_center(3) = List_all_comb_b3_cent(3,2) + + ! --- + + allocate(int_fit_v(n_points_final_grid)) + allocate(I_vec(ao_num,ao_num,n_points_final_grid)) + + I_vec = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + + call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, n_points_final_grid) + + do ipoint = 1, n_points_final_grid + I_vec(j,i,ipoint) = int_fit_v(ipoint) + enddo + enddo + enddo + + ! --- + + allocate(I_ref(ao_num,ao_num,n_points_final_grid)) + + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = 1, ao_num + + I_ref(j,i,ipoint) = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + enddo + enddo + enddo + + ! --- + + eps_ij = 1d-3 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + + i_exc = I_ref(i,j,ipoint) + i_num = I_vec(i,j,ipoint) + acc_ij = dabs(i_exc - i_num) + !acc_ij = dabs(i_exc - i_num) / dabs(i_exc) + if(acc_ij .gt. eps_ij) then + print *, ' problem in overlap_gauss_r12_ao_v on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + stop + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + enddo + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_vect_overlap_gauss_r12_ao + diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index e05492f0..a304324c 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -82,12 +82,77 @@ END_PROVIDER ! --- +BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) + ! + ! where r1 = r(ipoint) + ! + ! if J(r1,r2) = u12: + ! + ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1) + ! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] + ! = -int2_grad1_u12_ao(:,i,j,ipoint) + ! + ! if J(r1,r2) = u12 x v1 x v2 + ! + ! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ] + ! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ] + ! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) + ! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) + ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) + ! + ! + END_DOC + + implicit none + integer :: ipoint, i, j + double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 + + PROVIDE j1b_type + + if(j1b_type .eq. 3) then + + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + tmp0 = 0.5d0 * v_1b(ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) + + do j = 1, ao_num + do i = 1, ao_num + + tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) + + int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(1,i,j,ipoint) - tmp2 * tmp_x + int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint) - tmp2 * tmp_y + int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(3,i,j,ipoint) - tmp2 * tmp_z + enddo + enddo + enddo + + else + + int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao + + endif + +END_PROVIDER + +! --- BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)] BEGIN_DOC ! - ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij > + ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) . \grad_1 | ij > ! ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) ! @@ -99,11 +164,14 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, integer :: ipoint, i, j, k, l double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z double precision :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz + double precision :: ao_j_r, ao_l_r, ao_l_dx, ao_l_dy, ao_l_dz double precision, allocatable :: ac_mat(:,:,:,:) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) ac_mat = 0.d0 + ! --- + do ipoint = 1, n_points_final_grid weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) @@ -133,12 +201,47 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, enddo enddo enddo + + ! --- + + !do ipoint = 1, n_points_final_grid + ! weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) + + ! do l = 1, ao_num + ! ao_l_r = weight1 * aos_in_r_array_transp (ipoint,l) + ! ao_l_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,1) + ! ao_l_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,2) + ! ao_l_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,3) + + ! do j = 1, ao_num + ! ao_j_r = aos_in_r_array_transp(ipoint,j) + + ! tmp_x = ao_j_r * ao_l_dx - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,1) + ! tmp_y = ao_j_r * ao_l_dy - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,2) + ! tmp_z = ao_j_r * ao_l_dz - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,3) + + ! do i = 1, ao_num + ! do k = 1, ao_num + + ! contrib_x = int2_grad1_u12_ao(1,k,i,ipoint) * tmp_x + ! contrib_y = int2_grad1_u12_ao(2,k,i,ipoint) * tmp_y + ! contrib_z = int2_grad1_u12_ao(3,k,i,ipoint) * tmp_z + + ! ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z + ! enddo + ! enddo + ! enddo + ! enddo + !enddo + + ! --- do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + !tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) enddo enddo enddo diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f index 5188581a..a6f19c05 100644 --- a/src/scf_utils/diagonalize_fock.irp.f +++ b/src/scf_utils/diagonalize_fock.irp.f @@ -57,7 +57,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num) do i = elec_beta_num+1, elec_alpha_num F(i,i) += 0.5d0*level_shift enddo - do i = elec_alpha_num+1, mo_num F(i,i) += level_shift enddo diff --git a/src/scf_utils/diis.irp.f b/src/scf_utils/diis.irp.f index 713de1b3..00d4addb 100644 --- a/src/scf_utils/diis.irp.f +++ b/src/scf_utils/diis.irp.f @@ -1,3 +1,5 @@ +! --- + BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero ] implicit none BEGIN_DOC @@ -12,6 +14,8 @@ BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero ] END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_AO, (AO_num, AO_num)] implicit none BEGIN_DOC @@ -60,6 +64,8 @@ BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_AO, (AO_num, AO_num)] END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_MO, (mo_num, mo_num)] implicit none begin_doc @@ -69,6 +75,7 @@ BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_MO, (mo_num, mo_num)] FPS_SPF_Matrix_MO, size(FPS_SPF_Matrix_MO,1)) END_PROVIDER +! --- BEGIN_PROVIDER [ double precision, eigenvalues_Fock_matrix_AO, (AO_num) ] &BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_AO, (AO_num,AO_num) ] @@ -137,3 +144,107 @@ END_PROVIDER END_PROVIDER +! --- + +!BEGIN_PROVIDER [double precision, error_diis_Fmo, (ao_num, ao_num)] +! +! BEGIN_DOC +! ! +! ! error_diis_Fmo = (S x C) x [F_mo x \eta_occ - \eta_occ x F_mo] x (S x C).T +! ! +! ! \eta_occ is the matrix of occupation : \eta_occ = \eta_occ(alpha) + \eta_occ(beta) +! ! +! END_DOC +! +! implicit none +! integer :: i, j +! double precision, allocatable :: tmp(:,:) +! +! provide Fock_matrix_mo +! +! allocate(tmp(mo_num,mo_num)) +! tmp = 0.d0 +! +! ! F_mo x \eta_occ(alpha) - \eta_occ x F_mo(alpha) +! do j = 1, elec_alpha_num +! do i = elec_alpha_num + 1, mo_num +! tmp(i,j) = Fock_matrix_mo(i,j) +! enddo +! enddo +! do j = elec_alpha_num + 1, mo_num +! do i = 1, elec_alpha_num +! tmp(i,j) = -Fock_matrix_mo(i,j) +! enddo +! enddo +! +! ! F_mo x \eta_occ(beta) - \eta_occ x F_mo(beta) +! do j = 1, elec_beta_num +! do i = elec_beta_num + 1, mo_num +! tmp(i,j) += Fock_matrix_mo(i,j) +! enddo +! enddo +! do j = elec_beta_num + 1, mo_num +! do i = 1, elec_beta_num +! tmp(i,j) -= Fock_matrix_mo(i,j) +! enddo +! enddo +! +! call mo_to_ao(tmp, size(tmp, 1), error_diis_Fmo, size(error_diis_Fmo, 1)) +! +! deallocate(tmp) +! +!END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, error_diis_Fmo, (mo_num, mo_num)] + + BEGIN_DOC + ! + ! error_diis_Fmo = [F_mo x \eta_occ - \eta_occ x F_mo] + ! + ! \eta_occ is the matrix of occupation : \eta_occ = \eta_occ(alpha) + \eta_occ(beta) + ! + END_DOC + + implicit none + integer :: i, j + double precision, allocatable :: tmp(:,:) + + provide Fock_matrix_mo + + error_diis_Fmo = 0.d0 + + ! F_mo x \eta_occ(alpha) - \eta_occ x F_mo(alpha) + do j = 1, elec_alpha_num + do i = elec_alpha_num + 1, mo_num + error_diis_Fmo(i,j) += Fock_matrix_mo(i,j) + enddo + enddo + do j = elec_alpha_num + 1, mo_num + do i = 1, elec_alpha_num + error_diis_Fmo(i,j) -= Fock_matrix_mo(i,j) + enddo + enddo + + ! F_mo x \eta_occ(beta) - \eta_occ x F_mo(beta) + do j = 1, elec_beta_num + do i = elec_beta_num + 1, mo_num + error_diis_Fmo(i,j) += Fock_matrix_mo(i,j) + enddo + enddo + do j = elec_beta_num + 1, mo_num + do i = 1, elec_beta_num + error_diis_Fmo(i,j) -= Fock_matrix_mo(i,j) + enddo + enddo + + !allocate(tmp(ao_num,ao_num)) + !call mo_to_ao(error_diis_Fmo, size(error_diis_Fmo, 1), tmp, size(tmp, 1)) + !call ao_to_mo(tmp, size(tmp, 1), error_diis_Fmo, size(error_diis_Fmo, 1)) + !deallocate(tmp) + +END_PROVIDER + +! --- + diff --git a/src/scf_utils/rh_scf_mo.irp.f b/src/scf_utils/rh_scf_mo.irp.f new file mode 100644 index 00000000..5b70fb9c --- /dev/null +++ b/src/scf_utils/rh_scf_mo.irp.f @@ -0,0 +1,308 @@ +! --- + +subroutine Roothaan_Hall_SCF_MO() + + BEGIN_DOC + ! + ! Roothaan-Hall algorithm for SCF Hartree-Fock calculation + ! + END_DOC + + implicit none + + double precision :: energy_SCF, energy_SCF_previous, Delta_energy_SCF + double precision :: max_error_DIIS + double precision, allocatable :: Fock_matrix_DIIS(:,:,:), error_matrix_DIIS(:,:,:) + + integer :: iteration_SCF, dim_DIIS, index_dim_DIIS + + integer :: i, j + double precision :: level_shift_save + double precision, allocatable :: mo_coef_save(:,:) + + logical, external :: qp_stop + + PROVIDE ao_md5 mo_occ level_shift + + allocate( mo_coef_save(ao_num,mo_num) & + , Fock_matrix_DIIS (mo_num,mo_num,max_dim_DIIS) & + , error_matrix_DIIS(mo_num,mo_num,max_dim_DIIS) ) + + Fock_matrix_DIIS = 0.d0 + error_matrix_DIIS = 0.d0 + mo_coef_save = 0.d0 + + call write_time(6) + + print*,'energy of the guess = ',SCF_energy + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + ' N ', 'energy ', 'energy diff ', 'DIIS error ', 'Level shift ' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + +! Initialize energies and density matrices + energy_SCF_previous = SCF_energy + Delta_energy_SCF = 1.d0 + iteration_SCF = 0 + dim_DIIS = 0 + max_error_DIIS = 1.d0 + + +! +! Start of main SCF loop +! + PROVIDE Fock_matrix_mo error_diis_Fmo + + do while ( & + ( (max_error_DIIS > threshold_DIIS_nonzero) .or. & + (dabs(Delta_energy_SCF) > thresh_SCF) & + ) .and. (iteration_SCF < n_it_SCF_max) ) + + iteration_SCF += 1 + if(frozen_orb_scf) then + call initialize_mo_coef_begin_iteration + endif + + dim_DIIS = min(dim_DIIS+1, max_dim_DIIS) + + if( (scf_algorithm == 'DIIS_MO').and.(dabs(Delta_energy_SCF) > 1.d-6)) then + !if(scf_algorithm == 'DIIS_MO') then + + index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS) + 1 + do j = 1, mo_num + do i = 1, mo_num + Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_mo(i,j) + error_matrix_DIIS(i,j,index_dim_DIIS) = error_diis_Fmo(i,j) + enddo + enddo + + call extrapolate_Fock_matrix_mo(error_matrix_DIIS, Fock_matrix_DIIS, Fock_matrix_mo, size(Fock_matrix_mo, 1), iteration_SCF, dim_DIIS) + do i = 1, mo_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo + TOUCH Fock_matrix_mo fock_matrix_diag_mo + endif + + mo_coef = eigenvectors_Fock_matrix_mo + if(frozen_orb_scf) then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + + TOUCH mo_coef + + max_error_DIIS = maxval(Abs(error_diis_Fmo)) + + energy_SCF = SCF_energy + Delta_energy_SCF = energy_SCF - energy_SCF_previous + + if( (SCF_algorithm == 'DIIS_MO') .and. (Delta_energy_SCF > 0.d0) ) then + Fock_matrix_MO(1:mo_num,1:mo_num) = Fock_matrix_DIIS(1:mo_num,1:mo_num,index_dim_DIIS) + do i = 1, mo_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo + TOUCH Fock_matrix_mo fock_matrix_diag_mo + mo_coef = eigenvectors_Fock_matrix_mo + max_error_DIIS = maxval(Abs(error_diis_Fmo)) + energy_SCF = SCF_energy + Delta_energy_SCF = energy_SCF - energy_SCF_previous + endif + + level_shift_save = level_shift + mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num) + do while(Delta_energy_SCF > 0.d0) + mo_coef(1:ao_num,1:mo_num) = mo_coef_save(1:ao_num,1:mo_num) + if(level_shift <= .1d0) then + level_shift = 1.d0 + else + level_shift = level_shift * 3.0d0 + endif + TOUCH mo_coef level_shift + mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_mo(1:ao_num,1:mo_num) + if(frozen_orb_scf) then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + Delta_energy_SCF = SCF_energy - energy_SCF_previous + energy_SCF = SCF_energy + if(level_shift-level_shift_save > 40.d0) then + level_shift = level_shift_save * 4.d0 + SOFT_TOUCH level_shift + exit + endif + + dim_DIIS=0 + enddo + + level_shift = level_shift * 0.5d0 + SOFT_TOUCH level_shift + energy_SCF_previous = energy_SCF + +! Print results at the end of each iteration + + write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') & + iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS + + if(Delta_energy_SCF < 0.d0) then + call save_mos + endif + + if(qp_stop()) exit + enddo + +! +! End of Main SCF loop +! + + if(iteration_SCF < n_it_SCF_max) then + mo_label = 'Canonical' + endif + + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,*) + + if(.not.frozen_orb_scf)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo, size(Fock_matrix_mo, 1), size(Fock_matrix_mo, 2), mo_label, 1, .true.) + call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10) + call orthonormalize_mos + call save_mos + endif + + call write_double(6, energy_SCF, 'SCF energy') + + call write_time(6) + +end + +! --- + +subroutine extrapolate_Fock_matrix_mo(error_matrix_DIIS, Fock_matrix_DIIS, Fock_matrix_MO_, size_Fock_matrix_MO, iteration_SCF, dim_DIIS) + + BEGIN_DOC + ! Compute the extrapolated Fock matrix using the DIIS procedure + END_DOC + + implicit none + + integer,intent(inout) :: dim_DIIS + double precision,intent(in) :: Fock_matrix_DIIS(mo_num,mo_num,dim_DIIS), error_matrix_DIIS(mo_num,mo_num,dim_DIIS) + integer,intent(in) :: iteration_SCF, size_Fock_matrix_MO + double precision,intent(inout):: Fock_matrix_MO_(size_Fock_matrix_MO,mo_num) + + double precision,allocatable :: B_matrix_DIIS(:,:),X_vector_DIIS(:) + double precision,allocatable :: C_vector_DIIS(:) + + double precision,allocatable :: scratch(:,:) + integer :: i,j,k,l,i_DIIS,j_DIIS + double precision :: rcond, ferr, berr + integer, allocatable :: iwork(:) + integer :: lwork + + if(dim_DIIS < 1) then + return + endif + + allocate( & + B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), & + X_vector_DIIS(dim_DIIS+1), & + C_vector_DIIS(dim_DIIS+1), & + scratch(mo_num,mo_num) & + ) + + ! Compute the matrices B and X + B_matrix_DIIS(:,:) = 0.d0 + do j = 1, dim_DIIS + j_DIIS = min(dim_DIIS, mod(iteration_SCF-j, max_dim_DIIS) + 1) + + do i = 1, dim_DIIS + i_DIIS = min(dim_DIIS, mod(iteration_SCF-i, max_dim_DIIS) + 1) + + ! Compute product of two errors vectors + do l = 1, mo_num + do k = 1, mo_num + B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + error_matrix_DIIS(k,l,i_DIIS) * error_matrix_DIIS(k,l,j_DIIS) + enddo + enddo + + enddo + enddo + +! Pad B matrix and build the X matrix + + C_vector_DIIS(:) = 0.d0 + do i = 1, dim_DIIS + B_matrix_DIIS(i,dim_DIIS+1) = -1.d0 + B_matrix_DIIS(dim_DIIS+1,i) = -1.d0 + enddo + C_vector_DIIS(dim_DIIS+1) = -1.d0 + + deallocate(scratch) + +! Estimate condition number of B + double precision :: anorm + integer :: info + integer,allocatable :: ipiv(:) + double precision, allocatable :: AF(:,:) + double precision, external :: dlange + + lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5) + allocate(AF(dim_DIIS+1,dim_DIIS+1)) + allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) ) + allocate(scratch(lwork,1)) + scratch(:,1) = 0.d0 + + anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1)) + + AF(:,:) = B_matrix_DIIS(:,:) + call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + call dgecon( '1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + if(rcond < 1.d-14) then + dim_DIIS = 0 + return + endif + + ! solve the linear system C = B.X + + X_vector_DIIS = C_vector_DIIS + call dgesv(dim_DIIS+1 , 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv, X_vector_DIIS, size(X_vector_DIIS, 1), info) + + deallocate(scratch, AF, iwork) + + if(info < 0) then + stop 'bug in DIIS_MO' + endif + + ! Compute extrapolated Fock matrix + + + !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (mo_num > 200) + do j = 1, mo_num + do i = 1, mo_num + Fock_matrix_MO_(i,j) = 0.d0 + enddo + do k = 1, dim_DIIS + if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle + do i = 1, mo_num + ! FPE here + Fock_matrix_MO_(i,j) = Fock_matrix_MO_(i,j) + X_vector_DIIS(k) * Fock_matrix_DIIS(i,j,dim_DIIS-k+1) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + diff --git a/src/scf_utils/rh_scf_modif.irp.f b/src/scf_utils/rh_scf_modif.irp.f new file mode 100644 index 00000000..c63871f3 --- /dev/null +++ b/src/scf_utils/rh_scf_modif.irp.f @@ -0,0 +1,196 @@ +subroutine Roothaan_Hall_SCF_MODIF + +BEGIN_DOC +! Roothaan-Hall algorithm for SCF Hartree-Fock calculation +END_DOC + + implicit none + + double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF + double precision :: max_error_DIIS,max_error_DIIS_alpha,max_error_DIIS_beta + double precision, allocatable :: Fock_matrix_DIIS(:,:,:),error_matrix_DIIS(:,:,:) + + integer :: iteration_SCF,dim_DIIS,index_dim_DIIS + + integer :: i,j + logical, external :: qp_stop + double precision, allocatable :: mo_coef_save(:,:) + + PROVIDE ao_md5 mo_occ level_shift + + allocate(mo_coef_save(ao_num,mo_num), & + Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS), & + error_matrix_DIIS(ao_num,ao_num,max_dim_DIIS) & + ) + + Fock_matrix_DIIS = 0.d0 + error_matrix_DIIS = 0.d0 + mo_coef_save = 0.d0 + + call write_time(6) + + print*,'energy of the guess = ',SCF_energy + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + ' N ', 'energy ', 'energy diff ', 'DIIS error ', 'Level shift ' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + +! Initialize energies and density matrices + energy_SCF_previous = SCF_energy + Delta_energy_SCF = 1.d0 + iteration_SCF = 0 + dim_DIIS = 0 + max_error_DIIS = 1.d0 + + +! +! Start of main SCF loop +! + PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO + + do while ( & + ( (max_error_DIIS > threshold_DIIS_nonzero) .or. & + (dabs(Delta_energy_SCF) > thresh_SCF) & + ) .and. (iteration_SCF < n_it_SCF_max) ) + +! Increment cycle number + + iteration_SCF += 1 + if(frozen_orb_scf)then + call initialize_mo_coef_begin_iteration + endif + +! Current size of the DIIS space + + dim_DIIS = min(dim_DIIS+1,max_dim_DIIS) + + if( (scf_algorithm == 'DIIS_MODIF') .and. (dabs(Delta_energy_SCF) > 1.d-6) ) then + !if(scf_algorithm == 'DIIS_MODIF') then + + ! Store Fock and error matrices at each iteration + index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1 + do j=1,ao_num + do i=1,ao_num + Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO(i,j) + error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO(i,j) + enddo + enddo + + ! Compute the extrapolated Fock matrix + + call extrapolate_Fock_matrix( & + error_matrix_DIIS,Fock_matrix_DIIS, & + Fock_matrix_AO,size(Fock_matrix_AO,1), & + iteration_SCF,dim_DIIS & + ) + call ao_to_mo(Fock_matrix_AO, size(Fock_matrix_AO, 1), Fock_matrix_MO, size(Fock_matrix_MO, 1)) + do i = 1, mo_num + Fock_matrix_diag_MO(i) = Fock_matrix_MO(i,i) + enddo + TOUCH Fock_matrix_MO Fock_matrix_diag_MO + + !Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0 + !Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 + !TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta + endif + + MO_coef = eigenvectors_Fock_matrix_MO + if(frozen_orb_scf)then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + + TOUCH MO_coef + +! Calculate error vectors + + max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO)) + +! SCF energy + + energy_SCF = SCF_energy + Delta_energy_SCF = energy_SCF - energy_SCF_previous + if( (SCF_algorithm == 'DIIS_MODIF') .and. (Delta_energy_SCF > 0.d0) ) then + Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS(1:ao_num,1:ao_num,index_dim_DIIS) + call ao_to_mo(Fock_matrix_AO, size(Fock_matrix_AO, 1), Fock_matrix_MO, size(Fock_matrix_MO, 1)) + do i = 1, mo_num + Fock_matrix_diag_MO(i) = Fock_matrix_MO(i,i) + enddo + TOUCH Fock_matrix_MO Fock_matrix_diag_MO + + !Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0 + !Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 + !TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta + endif + + double precision :: level_shift_save + level_shift_save = level_shift + mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num) + do while (Delta_energy_SCF > 0.d0) + mo_coef(1:ao_num,1:mo_num) = mo_coef_save + if (level_shift <= .1d0) then + level_shift = 1.d0 + else + level_shift = level_shift * 3.0d0 + endif + TOUCH mo_coef level_shift + mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num) + if(frozen_orb_scf)then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + Delta_energy_SCF = SCF_energy - energy_SCF_previous + energy_SCF = SCF_energy + if (level_shift-level_shift_save > 40.d0) then + level_shift = level_shift_save * 4.d0 + SOFT_TOUCH level_shift + exit + endif + + dim_DIIS=0 + enddo + + level_shift = level_shift * 0.5d0 + SOFT_TOUCH level_shift + energy_SCF_previous = energy_SCF + +! Print results at the end of each iteration + + write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') & + iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS + + if (Delta_energy_SCF < 0.d0) then + call save_mos + endif + if (qp_stop()) exit + + enddo + + if (iteration_SCF < n_it_SCF_max) then + mo_label = 'Canonical' + endif +! +! End of Main SCF loop +! + + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,*) + + if(.not.frozen_orb_scf)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), & + size(Fock_matrix_mo,2),mo_label,1,.true.) + call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) + call orthonormalize_mos + call save_mos + endif + + call write_double(6, energy_SCF, 'SCF energy') + + call write_time(6) + +end + diff --git a/src/scf_utils/rh_scf_simple.irp.f b/src/scf_utils/rh_scf_simple.irp.f new file mode 100644 index 00000000..59b12749 --- /dev/null +++ b/src/scf_utils/rh_scf_simple.irp.f @@ -0,0 +1,130 @@ +subroutine Roothaan_Hall_SCF_Simple + +BEGIN_DOC +! Roothaan-Hall algorithm for SCF Hartree-Fock calculation +END_DOC + + implicit none + + integer :: iteration_SCF, dim_DIIS + double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF + double precision :: max_error_DIIS + + integer :: i,j + logical, external :: qp_stop + double precision, allocatable :: mo_coef_save(:,:) + + PROVIDE ao_md5 mo_occ level_shift + + allocate(mo_coef_save(ao_num,mo_num)) + + + dim_DIIS = 0 + mo_coef_save = 0.d0 + + call write_time(6) + + print*,'energy of the guess = ',SCF_energy + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + ' N ', 'energy ', 'energy diff ', 'DIIS error ', 'Level shift ' + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + +! Initialize energies and density matrices + energy_SCF_previous = SCF_energy + Delta_energy_SCF = 1.d0 + iteration_SCF = 0 + max_error_DIIS = 1.d0 + + do while ( & + ( (max_error_DIIS > threshold_DIIS_nonzero) .or. & + (dabs(Delta_energy_SCF) > thresh_SCF) & + ) .and. (iteration_SCF < n_it_SCF_max) ) + + iteration_SCF += 1 + if(frozen_orb_scf)then + call initialize_mo_coef_begin_iteration + endif + + MO_coef = eigenvectors_Fock_matrix_MO + if(frozen_orb_scf)then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH MO_coef + +! Calculate error vectors + max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO)) + +! SCF energy + + energy_SCF = SCF_energy + Delta_energy_SCF = energy_SCF - energy_SCF_previous + + double precision :: level_shift_save + level_shift_save = level_shift + mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num) + do while (Delta_energy_SCF > 0.d0) + mo_coef(1:ao_num,1:mo_num) = mo_coef_save + if (level_shift <= .1d0) then + level_shift = 1.d0 + else + level_shift = level_shift * 3.0d0 + endif + TOUCH mo_coef level_shift + mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num) + if(frozen_orb_scf)then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + Delta_energy_SCF = SCF_energy - energy_SCF_previous + energy_SCF = SCF_energy + if (level_shift-level_shift_save > 40.d0) then + level_shift = level_shift_save * 4.d0 + SOFT_TOUCH level_shift + exit + endif + + enddo + + level_shift = level_shift * 0.5d0 + SOFT_TOUCH level_shift + energy_SCF_previous = energy_SCF + +! Print results at the end of each iteration + + write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') & + iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS + + if(Delta_energy_SCF < 0.d0) then + call save_mos + endif + if(qp_stop()) exit + + enddo + + if (iteration_SCF < n_it_SCF_max) then + mo_label = 'Canonical' + endif + + write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & + '====','================','================','================','================' + write(6,*) + + if(.not.frozen_orb_scf)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), & + size(Fock_matrix_mo,2),mo_label,1,.true.) + call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) + call orthonormalize_mos + call save_mos + endif + + call write_double(6, energy_SCF, 'SCF energy') + + call write_time(6) + +end + diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 56a1ed8e..45522079 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -66,7 +66,8 @@ END_DOC dim_DIIS = min(dim_DIIS+1,max_dim_DIIS) - if ( (scf_algorithm == 'DIIS').and.(dabs(Delta_energy_SCF) > 1.d-6) ) then + if( (scf_algorithm == 'DIIS') .and. (dabs(Delta_energy_SCF) > 1.d-6)) then + !if(scf_algorithm == 'DIIS') then ! Store Fock and error matrices at each iteration index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1 diff --git a/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f index 60201f5f..eb812401 100644 --- a/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f +++ b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f @@ -14,21 +14,36 @@ program save_bitcpsileft_for_qmcchem e_ref = 0.d0 iunit = 13 - open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write') - call ezfio_has_fci_energy_pt2(exists) - - if(.not.exists) then - call ezfio_has_fci_energy(exists) + open(unit=iunit, file=trim(ezfio_filename)//'/simulation/e_ref', action='write') + call ezfio_has_fci_energy_pt2(exists) if(.not.exists) then - call ezfio_has_tc_scf_bitc_energy(exists) - if(exists) then - call ezfio_get_tc_scf_bitc_energy(e_ref) + + call ezfio_has_fci_energy(exists) + if(.not.exists) then + + call ezfio_has_cisd_energy(exists) + if(.not.exists) then + + call ezfio_has_tc_scf_bitc_energy(exists) + if(exists) then + call ezfio_get_tc_scf_bitc_energy(e_ref) + endif + + else + call ezfio_get_cisd_energy(e_ref) + endif + + else + call ezfio_get_fci_energy(e_ref) endif + + else + call ezfio_get_fci_energy_pt2(e_ref) endif - endif - write(iunit,*) e_ref + write(iunit,*) e_ref + close(iunit) end diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 0cbdb753..26d75ad4 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -158,7 +158,7 @@ default: 0. type: character*(32) doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] interface: ezfio,provider,ocaml -default: DIIS +default: Simple [im_thresh_tcscf] type: Threshold diff --git a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f new file mode 100644 index 00000000..048255f6 --- /dev/null +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -0,0 +1,483 @@ + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_cs ...' + call wall_time(ti) + + fock_3e_uhf_mo_cs = 0.d0 + + do a = 1, mo_num + do b = 1, mo_num + + do j = 1, elec_beta_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_cs(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - 2.d0 * I_bij_aji & + - 2.d0 * I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j, o + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_a ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - 2.d0 * I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + ! --- + + do j = o, elec_alpha_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + ! --- + + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j, o + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_b ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_iaj ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_jia ) + + enddo + enddo + + ! --- + + do j = o, elec_alpha_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( I_bij_aij & + - I_bij_aji ) + + enddo + enddo + + ! --- + + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + double precision, allocatable :: f_tmp(:,:) + + print *, ' PROVIDING fock_3e_uhf_ao_a ...' + call wall_time(ti) + + fock_3e_uhf_ao_a = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a) + + allocate(f_tmp(ao_num,ao_num)) + f_tmp = 0.d0 + + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + do mu = 1, ao_num + do nu = 1, ao_num + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_a * dm_dk_a * i_mugd_eknu & + + dm_ge_a * dm_dk_a * i_mugd_knue & + - dm_ge_a * dm_dk * i_mugd_enuk & + - dm_ge * dm_dk_a * i_mugd_kenu & + - dm_ge_a * dm_dk_a * i_mugd_nuke & + - dm_ge_b * dm_dk_b * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do mu = 1, ao_num + do nu = 1, ao_num + fock_3e_uhf_ao_a(mu,nu) += f_tmp(mu,nu) + enddo + enddo + !$OMP END CRITICAL + + deallocate(f_tmp) + !$OMP END PARALLEL + +! TODO +! !$OMP PARALLEL DEFAULT (NONE) & +! !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & +! !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & +! !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a) +! !$OMP DO +! do g = 1, ao_num +! do e = 1, ao_num +! dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) +! dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) +! dm_ge = dm_ge_a + dm_ge_b +! do d = 1, ao_num +! do k = 1, ao_num +! dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) +! dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) +! dm_dk = dm_dk_a + dm_dk_b +! do mu = 1, ao_num +! do nu = 1, ao_num +! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) +! fock_3e_uhf_ao_a(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & +! + dm_ge_a * dm_dk_a * i_mugd_eknu & +! + dm_ge_a * dm_dk_a * i_mugd_knue & +! - dm_ge_a * dm_dk * i_mugd_enuk & +! - dm_ge * dm_dk_a * i_mugd_kenu & +! - dm_ge_a * dm_dk_a * i_mugd_nuke & +! - dm_ge_b * dm_dk_b * i_mugd_nuke ) +! enddo +! enddo +! enddo +! enddo +! enddo +! enddo +! !$OMP END DO +! !$OMP END PARALLEL + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + double precision, allocatable :: f_tmp(:,:) + + print *, ' PROVIDING fock_3e_uhf_ao_b ...' + call wall_time(ti) + + fock_3e_uhf_ao_b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b) + + allocate(f_tmp(ao_num,ao_num)) + f_tmp = 0.d0 + + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + do mu = 1, ao_num + do nu = 1, ao_num + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_b * dm_dk_b * i_mugd_eknu & + + dm_ge_b * dm_dk_b * i_mugd_knue & + - dm_ge_b * dm_dk * i_mugd_enuk & + - dm_ge * dm_dk_b * i_mugd_kenu & + - dm_ge_b * dm_dk_b * i_mugd_nuke & + - dm_ge_a * dm_dk_a * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do mu = 1, ao_num + do nu = 1, ao_num + fock_3e_uhf_ao_b(mu,nu) += f_tmp(mu,nu) + enddo + enddo + !$OMP END CRITICAL + + deallocate(f_tmp) + !$OMP END PARALLEL + +! TODO +! !$OMP PARALLEL DO DEFAULT (NONE) & +! !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & +! !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & +! !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b) +! do g = 1, ao_num +! do e = 1, ao_num +! dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) +! dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) +! dm_ge = dm_ge_a + dm_ge_b +! do d = 1, ao_num +! do k = 1, ao_num +! dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) +! dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) +! dm_dk = dm_dk_a + dm_dk_b +! do mu = 1, ao_num +! do nu = 1, ao_num +! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) +! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) +! fock_3e_uhf_ao_b(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & +! + dm_ge_b * dm_dk_b * i_mugd_eknu & +! + dm_ge_b * dm_dk_b * i_mugd_knue & +! - dm_ge_b * dm_dk * i_mugd_enuk & +! - dm_ge * dm_dk_b * i_mugd_kenu & +! - dm_ge_b * dm_dk_b * i_mugd_nuke & +! - dm_ge_a * dm_dk_a * i_mugd_nuke ) +! enddo +! enddo +! enddo +! enddo +! enddo +! enddo +! !$OMP END PARALLEL DO + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti + +END_PROVIDER + +! --- + diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 2a08e469..5981791c 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -31,13 +31,22 @@ density_b = TCSCF_density_matrix_ao_beta (l,j) density = density_a + density_b + !! rho(l,j) * < k l| T | i j> + !two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i) + !! rho(l,j) * < k l| T | i j> + !two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i) + !! rho_a(l,j) * < l k| T | i j> + !two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i) + !! rho_b(l,j) * < l k| T | i j> + !two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) + ! rho(l,j) * < k l| T | i j> - two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i) + two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) ! rho(l,j) * < k l| T | i j> - two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i) - ! rho_a(l,j) * < l k| T | i j> + two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(k,i,l,j) + ! rho_a(l,j) * < k l| T | j i> two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i) - ! rho_b(l,j) * < l k| T | i j> + ! rho_b(l,j) * < k l| T | j i> two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) enddo @@ -84,13 +93,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] END_DOC implicit none + double precision, allocatable :: tmp(:,:) if(bi_ortho) then + !allocate(tmp(ao_num,ao_num)) + !tmp = Fock_matrix_tc_ao_alpha + !if(three_body_h_tc) then + ! tmp += fock_3e_uhf_ao_a + !endif + !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1)) + !deallocate(tmp) + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) - if(three_body_h_tc.and.elec_alpha_num == elec_beta_num) then - Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + if(three_body_h_tc) then + !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a endif else @@ -110,14 +129,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] END_DOC implicit none + double precision, allocatable :: tmp(:,:) if(bi_ortho) then - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & - , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) + !allocate(tmp(ao_num,ao_num)) + !tmp = Fock_matrix_tc_ao_beta + !if(three_body_h_tc) then + ! tmp += fock_3e_uhf_ao_b + !endif + !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1)) + !deallocate(tmp) - if(three_body_h_tc.and.elec_alpha_num == elec_beta_num) then - Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & + , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) + if(three_body_h_tc) then + !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b endif else diff --git a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f index b0345957..f73171a3 100644 --- a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f +++ b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f @@ -1,202 +1,286 @@ +! --- + BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo,contrib - fock_a_tot_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth(a,i) - fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i) - fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i) + + implicit none + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i) + enddo enddo - enddo + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo,contrib - fock_b_tot_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth(a,i) - fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i) - fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i) + + implicit none + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i) + enddo enddo - enddo + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo, contrib - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - double precision :: new - fock_cs_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = 1, elec_beta_num -! call contrib_3e_sss(a,i,j,k,contrib_sss) -! call contrib_3e_soo(a,i,j,k,contrib_soo) -! call contrib_3e_sos(a,i,j,k,contrib_sos) -! contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - ! negative terms :: exchange contrib - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) & - -1.5d0 * exch_13_int - exch_23_int - fock_cs_3e_bi_orth(a,i) += new + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_cs_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = 1, elec_beta_num + + !!call contrib_3e_sss(a,i,j,k,contrib_sss) + !!call contrib_3e_soo(a,i,j,k,contrib_soo) + !!call contrib_3e_sos(a,i,j,k,contrib_sos) + !!contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) -1.5d0 * exch_13_int - exch_23_int + + fock_cs_3e_bi_orth(a,i) += new + enddo + enddo enddo - enddo - enddo - enddo - fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth + + fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth END_PROVIDER +! --- BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo, contrib - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - double precision :: new - fock_a_tmp1_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = elec_beta_num + 1, elec_alpha_num - do k = 1, elec_beta_num - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) & - + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int) - enddo - enddo + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_beta_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int) + enddo + enddo + enddo enddo - enddo - fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + + fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss - fock_a_tmp2_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - do j = 1, elec_alpha_num - do k = elec_beta_num+1, elec_alpha_num - call contrib_3e_sss(a,i,j,k,contrib_sss) - fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_alpha_num + do k = elec_beta_num+1, elec_alpha_num + call contrib_3e_sss(a, i, j, k, contrib_sss) + + fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss + enddo + enddo enddo - enddo enddo - enddo + END_PROVIDER - - - +! --- BEGIN_PROVIDER [double precision, fock_b_tmp1_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int - double precision :: new - fock_b_tmp1_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = elec_beta_num+1, elec_alpha_num - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int - enddo - enddo + implicit none + integer :: i, a, j, k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_beta_num + do k = elec_beta_num+1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + + fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int + enddo + enddo + enddo enddo - enddo - fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho + + fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_soo - fock_b_tmp2_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - do j = elec_beta_num + 1, elec_alpha_num - do k = 1, elec_alpha_num - call contrib_3e_soo(a,i,j,k,contrib_soo) - fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo + + implicit none + integer :: i, a, j, k + double precision :: contrib_soo + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_soo(a, i, j, k, contrib_soo) + + fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo + enddo + enddo enddo - enddo enddo - enddo + END_PROVIDER -subroutine contrib_3e_sss(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the pure same spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - integral = direct_int + c_3_int + c_minus_3_int - ! negative terms :: exchange contrib - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - integral += - exch_13_int - exch_23_int - exch_12_int - integral = -integral +! --- + +subroutine contrib_3e_sss(a, i, j, k, integral) + + BEGIN_DOC + ! returns the pure same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + + PROVIDE mo_l_coef mo_r_coef + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + integral += - exch_13_int - exch_23_int - exch_12_int + + integral = -integral + end +! --- + subroutine contrib_3e_soo(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_23_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 - integral = direct_int - exch_23_int - integral = -integral + + BEGIN_DOC + ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_23_int + + PROVIDE mo_l_coef mo_r_coef + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + integral = direct_int - exch_23_int + + integral = -integral + end -subroutine contrib_3e_sos(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_13_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 - integral = direct_int - exch_13_int - integral = -integral +! --- + +subroutine contrib_3e_sos(a, i, j, k, integral) + + BEGIN_DOC + ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + PROVIDE mo_l_coef mo_r_coef + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + integral = direct_int - exch_13_int + + integral = -integral + end + +! --- + diff --git a/src/tc_scf/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f index 597c3e67..0312df5f 100644 --- a/src/tc_scf/rh_tcscf.irp.f +++ b/src/tc_scf/rh_tcscf.irp.f @@ -67,10 +67,9 @@ subroutine rh_tcscf() iteration_TCSCF += 1 if(iteration_TCSCF > n_it_TCSCF_max) then print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max - exit + stop endif - ! current size of the DIIS space dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) ! --- @@ -86,10 +85,7 @@ subroutine rh_tcscf() enddo enddo - ! Compute the extrapolated Fock matrix - call extrapolate_TC_Fock_matrix( e_DIIS, F_DIIS & - , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & - , iteration_TCSCF, dim_DIIS ) + call extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), iteration_TCSCF, dim_DIIS) Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot @@ -100,7 +96,6 @@ subroutine rh_tcscf() call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta , size(Fock_matrix_tc_ao_beta , 1) & , Fock_matrix_tc_mo_beta , size(Fock_matrix_tc_mo_beta , 1) ) TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta - endif ! --- @@ -121,9 +116,10 @@ subroutine rh_tcscf() ! --- - do while((dabs(delta_energy_tmp) > 0.1d0) .and. (iteration_TCSCF > 1)) -! print *, ' very big step : ', delta_energy_tmp -! print *, ' TC level shift = ', level_shift_TCSCF + do while((delta_gradie_tmp > 1.d-7) .and. (iteration_TCSCF > 1)) + !do while((dabs(delta_energy_tmp) > 0.5d0) .and. (iteration_TCSCF > 1)) + print *, ' very big or bad step : ', delta_energy_tmp, delta_gradie_tmp + print *, ' TC level shift = ', level_shift_TCSCF mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) @@ -139,7 +135,8 @@ subroutine rh_tcscf() mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) TOUCH mo_l_coef mo_r_coef - delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous if(level_shift_TCSCF - level_shift_save > 40.d0) then level_shift_TCSCF = level_shift_save * 4.d0 @@ -183,7 +180,7 @@ subroutine rh_tcscf() print *, ' 1-e TC energy = ', energy_TCSCF_1e print *, ' 2-e TC energy = ', energy_TCSCF_2e print *, ' 3-e TC energy = ', energy_TCSCF_3e - print *, ' |delta TC energy| = ', delta_energy_TCSCF + print *, ' |delta TC energy| = ', dabs(delta_energy_TCSCF) print *, ' TC gradient = ', gradie_TCSCF print *, ' delta TC gradient = ', delta_gradie_TCSCF print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF @@ -199,6 +196,9 @@ subroutine rh_tcscf() ! --- + print *, ' TCSCF DIIS converged !' + call print_energy_and_mos() + call write_time(6) deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS) diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 283ec2ae..fd11c48e 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -21,8 +21,11 @@ program tc_scf PROVIDE tcscf_algorithm if(tcscf_algorithm == 'DIIS') then call rh_tcscf() - else + elseif(tcscf_algorithm == 'Simple') then call simple_tcscf() + else + print *, ' not implemented yet', tcscf_algorithm + stop endif call minimize_tc_orb_angles() @@ -127,7 +130,7 @@ subroutine simple_tcscf() it += 1 if(it > n_it_tcscf_max) then print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max - exit + stop endif @@ -142,8 +145,10 @@ subroutine simple_tcscf() endif e_delta = dabs(TC_HF_energy - e_save) - print *, ' delta E = ', e_delta - print *, ' gradient = ', grad_non_hermit + print *, ' delta E = ', e_delta + print *, ' gradient = ', grad_non_hermit + print *, ' max TC DIIS error = ', maxval(abs(FQS_SQF_mo)) + !print *, ' gradient= ', grad_non_hermit_right !rho_new = TCSCF_bi_ort_dm_ao @@ -165,6 +170,8 @@ subroutine simple_tcscf() TOUCH mo_l_coef mo_r_coef call ezfio_set_tc_scf_bitc_energy(TC_HF_energy) + call test_fock_3e_uhf_mo() + print *, ' ***' print *, '' @@ -190,7 +197,7 @@ subroutine simple_tcscf() endif - print*,'Energy converged !' + print *, ' TCSCF Simple converged !' call print_energy_and_mos() deallocate(rho_old, rho_new) @@ -199,3 +206,64 @@ end subroutine simple_tcscf ! --- +subroutine test_fock_3e_uhf_mo() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + + thr_ih = 1d-12 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + !print *, ' difference on ', j, i + !print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + !print *, ' UHF : ', fock_3e_uhf_mo_a (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' norm_a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + !print *, ' difference on ', j, i + !print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + !print *, ' UHF : ', fock_3e_uhf_mo_b (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' norm_b = ', norm + print *, ' ' + + ! --- + +end subroutine test_fock_3e_uhf_mo() + diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index a81b09d5..ae88dac3 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -6,7 +6,7 @@ program test_ints implicit none - print *, 'starting ...' + print *, ' starting test_ints ...' my_grid_becke = .True. my_n_pt_r_grid = 30 @@ -14,6 +14,7 @@ program test_ints ! my_n_pt_r_grid = 15 ! small grid for quick debug ! my_n_pt_a_grid = 26 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + my_extra_grid_becke = .True. my_n_pt_r_extra_grid = 30 my_n_pt_a_extra_grid = 50 ! small extra_grid for quick debug @@ -44,8 +45,13 @@ program test_ints ! call test_tc_scf call test_int_gauss + !call test_fock_3e_uhf_ao() + call test_fock_3e_uhf_mo() + end +! --- + subroutine test_tc_scf implicit none integer :: i @@ -70,6 +76,8 @@ subroutine test_ao_tc_int_chemist ! provide tc_grad_and_lapl_ao_test end +! --- + subroutine routine_test_j1b implicit none integer :: i,icount,j @@ -332,13 +340,13 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2 double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: ints(:,:,:) allocate(ints(ao_num, ao_num, n_points_final_grid)) - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, ao_num - read(33,*)ints(j,i,ipoint) - enddo - enddo - enddo +! do ipoint = 1, n_points_final_grid +! do i = 1, ao_num +! do j = 1, ao_num +! read(33,*)ints(j,i,ipoint) +! enddo +! enddo +! enddo allocate(array(ao_num, ao_num, ao_num, ao_num)) array = 0.d0 @@ -563,10 +571,147 @@ subroutine routine_v_ij_u_cst_mu_j1b print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - end +! --- + +subroutine test_fock_3e_uhf_ao() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + double precision, allocatable :: fock_3e_uhf_ao_a_mo(:,:), fock_3e_uhf_ao_b_mo(:,:) + + thr_ih = 1d-7 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_ao_a fock_3e_uhf_ao_b + + ! --- + + allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a , size(fock_3e_uhf_ao_a , 1) & + , fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_a_mo(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_a_mo (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' ' + + deallocate(fock_3e_uhf_ao_a_mo) + + ! --- + + allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b , size(fock_3e_uhf_ao_b , 1) & + , fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_b_mo(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_b_mo (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' ' + + deallocate(fock_3e_uhf_ao_b_mo) + + ! --- + +end subroutine test_fock_3e_uhf_ao() + +! --- + +subroutine test_fock_3e_uhf_mo() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + + thr_ih = 1d-12 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_a (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' norm_a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_b (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' norm_b = ', norm + print *, ' ' + + ! --- + +end subroutine test_fock_3e_uhf_mo() + +! --- + subroutine test_total_grad_lapl implicit none integer :: i,j,ipoint,k,l @@ -702,3 +847,4 @@ subroutine test_int_gauss end + diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index f593cefb..5079daa7 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -48,7 +48,7 @@ end ! TODO remove dim -subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) +subroutine give_explicit_poly_and_gaussian(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim) BEGIN_DOC ! Transforms the product of @@ -65,19 +65,19 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, implicit none include 'constants.include.F' - integer, intent(in) :: dim - integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) - double precision, intent(in) :: alpha, beta ! exponents - double precision, intent(in) :: A_center(3) ! A center - double precision, intent(in) :: B_center (3) ! B center - double precision, intent(out) :: P_center(3) ! new center - double precision, intent(out) :: p ! new exponent - double precision, intent(out) :: fact_k ! constant factor - double precision, intent(out) :: P_new(0:max_dim,3)! polynomial - integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials + integer, intent(in) :: dim + integer, intent(in) :: a(3), b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) + double precision, intent(in) :: alpha, beta ! exponents + double precision, intent(in) :: A_center(3) ! A center + double precision, intent(in) :: B_center (3) ! B center + integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials + double precision, intent(out) :: P_center(3) ! new center + double precision, intent(out) :: p ! new exponent + double precision, intent(out) :: fact_k ! constant factor + double precision, intent(out) :: P_new(0:max_dim,3)! polynomial - double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) - integer :: n_new,i,j + integer :: n_new, i, j + double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b iorder(1) = 0 @@ -87,46 +87,46 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, P_new(0,2) = 0.d0 P_new(0,3) = 0.d0 !DIR$ FORCEINLINE - call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) - if (fact_k < thresh) then + call gaussian_product(alpha, A_center, beta, B_center, fact_k, p, P_center) + if(fact_k < thresh) then ! IF fact_k is too smal then: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef P_center = 0.d0 - p = 1.d+15 - fact_k = 0.d0 + p = 1.d+15 + fact_k = 0.d0 return endif !DIR$ FORCEINLINE - call recentered_poly2(P_a(0,1),A_center(1),P_center(1),a(1),P_b(0,1),B_center(1),P_center(1),b(1)) + call recentered_poly2(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1)) iorder(1) = a(1) + b(1) - do i=0,iorder(1) + do i = 0, iorder(1) P_new(i,1) = 0.d0 enddo - n_new=0 + n_new = 0 !DIR$ FORCEINLINE - call multiply_poly(P_a(0,1),a(1),P_b(0,1),b(1),P_new(0,1),n_new) + call multiply_poly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new) !DIR$ FORCEINLINE - call recentered_poly2(P_a(0,2),A_center(2),P_center(2),a(2),P_b(0,2),B_center(2),P_center(2),b(2)) + call recentered_poly2(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2)) iorder(2) = a(2) + b(2) - do i=0,iorder(2) + do i = 0, iorder(2) P_new(i,2) = 0.d0 enddo - n_new=0 + n_new = 0 !DIR$ FORCEINLINE - call multiply_poly(P_a(0,2),a(2),P_b(0,2),b(2),P_new(0,2),n_new) + call multiply_poly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new) !DIR$ FORCEINLINE - call recentered_poly2(P_a(0,3),A_center(3),P_center(3),a(3),P_b(0,3),B_center(3),P_center(3),b(3)) + call recentered_poly2(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3)) iorder(3) = a(3) + b(3) - do i=0,iorder(3) + do i = 0, iorder(3) P_new(i,3) = 0.d0 enddo - n_new=0 + n_new = 0 !DIR$ FORCEINLINE - call multiply_poly(P_a(0,3),a(3),P_b(0,3),b(3),P_new(0,3),n_new) + call multiply_poly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) end @@ -167,26 +167,33 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center, p, fact_k, io call gaussian_product_v(alpha, A_center, LD_A, beta, B_center, fact_k, p, P_center, n_points) - if ( ior(ior(b(1),b(2)),b(3)) == 0 ) then ! b == (0,0,0) - - lda = maxval(a) - ldb = 0 - allocate(P_a(n_points,0:lda,3), P_b(n_points,0:0,3)) - - call recentered_poly2_v0(P_a, lda, A_center, LD_A, P_center, a, P_b, B_center, P_center, n_points) + if(ior(ior(b(1), b(2)), b(3)) == 0) then ! b == (0,0,0) iorder(1:3) = a(1:3) + + lda = maxval(a) + allocate(P_a(n_points,0:lda,3)) + !ldb = 0 + !allocate(P_b(n_points,0:0,3)) + + !call recentered_poly2_v0(P_a, lda, A_center, LD_A, P_center, a, P_b, B_center, P_center, n_points) + call recentered_poly2_v0(P_a, lda, A_center, LD_A, P_center, a, n_points) + do ipoint = 1, n_points do xyz = 1, 3 - P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) + !P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) + P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) do i = 1, a(xyz) - P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) + !P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) + P_new(ipoint,i,xyz) = P_a(ipoint,i,xyz) enddo enddo enddo - return + deallocate(P_a) + !deallocate(P_b) + return endif lda = maxval(a) @@ -198,20 +205,27 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center, p, fact_k, io iorder(1:3) = a(1:3) + b(1:3) do xyz = 1, 3 - if (b(xyz) == 0) then + if(b(xyz) == 0) then + do ipoint = 1, n_points - P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) + !P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) + P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) do i = 1, a(xyz) - P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) + !P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) + P_new(ipoint,i,xyz) = P_a(ipoint,i,xyz) enddo enddo + else + do i = 0, iorder(xyz) do ipoint = 1, n_points P_new(ipoint,i,xyz) = 0.d0 enddo enddo + call multiply_poly_v(P_a(1,0,xyz), a(xyz), P_b(1,0,xyz), b(xyz), P_new(1,0,xyz), ldp, n_points) + endif enddo @@ -720,45 +734,57 @@ end subroutine recentered_poly2_v ! --- -subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, P_new2, x_B, x_Q, n_points) +!subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, P_new2, x_B, x_Q, n_points) +subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, n_points) BEGIN_DOC + ! ! Recenter two polynomials. Special case for b=(0,0,0) + ! + ! (x - A)^a (x - B)^0 = (x - P + P - A)^a (x - Q + Q - B)^0 + ! = (x - P + P - A)^a + ! END_DOC implicit none integer, intent(in) :: a(3), n_points, lda, LD_xA - double precision, intent(in) :: x_A(LD_xA,3) - double precision, intent(in) :: x_B(3) - double precision, intent(in) :: x_P(n_points,3), x_Q(n_points,3) - double precision, intent(out) :: P_new(n_points,0:lda,3), P_new2(n_points,3) + double precision, intent(in) :: x_A(LD_xA,3), x_P(n_points,3) + !double precision, intent(in) :: x_B(3), x_Q(n_points,3) + double precision, intent(out) :: P_new(n_points,0:lda,3) + !double precision, intent(out) :: P_new2(n_points,3) + integer :: i, j, k, l, xyz, ipoint, maxab(3) double precision :: fa - double precision, allocatable :: pows_a(:,:), pows_b(:,:) + double precision, allocatable :: pows_a(:,:) + !double precision, allocatable :: pows_b(:,:) double precision :: binom_func - maxab(1:3) = max(a(1:3),(/0,0,0/)) + maxab(1:3) = max(a(1:3), (/0,0,0/)) - allocate( pows_a(n_points,-2:maxval(maxab)+4), pows_b(n_points,-2:maxval(maxab)+4) ) + allocate(pows_a(n_points,-2:maxval(maxab)+4)) + !allocate(pows_b(n_points,-2:maxval(maxab)+4)) do xyz = 1, 3 - if (a(xyz)<0) cycle - do ipoint=1,n_points + if(a(xyz) < 0) cycle + + do ipoint = 1, n_points pows_a(ipoint,0) = 1.d0 pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz)) - pows_b(ipoint,0) = 1.d0 - pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) + !pows_b(ipoint,0) = 1.d0 + !pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) enddo - do i = 2,maxab(xyz) - do ipoint=1,n_points - pows_a(ipoint,i) = pows_a(ipoint,i-1)*pows_a(ipoint,1) - pows_b(ipoint,i) = pows_b(ipoint,i-1)*pows_b(ipoint,1) + + do i = 2, maxab(xyz) + do ipoint = 1, n_points + pows_a(ipoint,i) = pows_a(ipoint,i-1) * pows_a(ipoint,1) + !pows_b(ipoint,i) = pows_b(ipoint,i-1) * pows_b(ipoint,1) enddo enddo - do ipoint=1,n_points + + do ipoint = 1, n_points P_new (ipoint,0,xyz) = pows_a(ipoint,a(xyz)) - P_new2(ipoint,xyz) = pows_b(ipoint,0) + !P_new2(ipoint,xyz) = pows_b(ipoint,0) enddo do i = 1, min(a(xyz), 20) fa = binom_transp(a(xyz)-i, a(xyz)) @@ -775,11 +801,12 @@ subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, P_new2, x_B, x_Q, enddo !xyz - deallocate(pows_a, pows_b) + deallocate(pows_a) + !deallocate(pows_b) end subroutine recentered_poly2_v0 -!-- +! --- subroutine pol_modif_center(A_center, B_center, iorder, A_pol, B_pol) diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index c797c87e..cf417613 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -31,7 +31,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_ overlap_gaussian_x*= fact_p end +! --- +! TODO +! gaussian_product is called twice: in give_explicit_poly_and_gaussian and here subroutine overlap_gaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, overlap_y, overlap_z, overlap, dim) BEGIN_DOC @@ -45,51 +48,50 @@ subroutine overlap_gaussian_xyz(A_center, B_center, alpha, beta, power_A, power_ include 'constants.include.F' implicit none - integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials - double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions - double precision, intent(in) :: alpha,beta - integer,intent(in) :: power_A(3), power_B(3) ! power of the x1 functions - double precision, intent(out) :: overlap_x,overlap_y,overlap_z,overlap - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p - double precision :: F_integral_tab(0:max_dim) - integer :: iorder_p(3) - integer :: nmax - double precision :: F_integral + integer, intent(in) :: dim ! dimension maximum for the arrays representing the polynomials + integer, intent(in) :: power_A(3), power_B(3) ! power of the x1 functions + double precision, intent(in) :: A_center(3), B_center(3) ! center of the x1 functions + double precision, intent(in) :: alpha, beta + double precision, intent(out) :: overlap_x, overlap_y, overlap_z, overlap + integer :: i, nmax, iorder_p(3) + double precision :: P_new(0:max_dim,3), P_center(3), fact_p, p + double precision :: F_integral_tab(0:max_dim) + + double precision :: F_integral call give_explicit_poly_and_gaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim) - if(fact_p.lt.1d-20)then + if(fact_p .lt. 1d-20) then overlap_x = 1.d-10 overlap_y = 1.d-10 overlap_z = 1.d-10 - overlap = 1.d-10 + overlap = 1.d-10 return endif nmax = maxval(iorder_p) - do i = 0,nmax - F_integral_tab(i) = F_integral(i,p) + do i = 0, nmax + F_integral_tab(i) = F_integral(i, p) enddo overlap_x = P_new(0,1) * F_integral_tab(0) overlap_y = P_new(0,2) * F_integral_tab(0) overlap_z = P_new(0,3) * F_integral_tab(0) - integer :: i do i = 1,iorder_p(1) overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) enddo - call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1)) + call gaussian_product_x(alpha, A_center(1), beta, B_center(1), fact_p, p, P_center(1)) overlap_x *= fact_p - do i = 1,iorder_p(2) + do i = 1, iorder_p(2) overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) enddo - call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2)) + call gaussian_product_x(alpha, A_center(2), beta, B_center(2), fact_p, p, P_center(2)) overlap_y *= fact_p do i = 1,iorder_p(3) overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) enddo - call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3)) + call gaussian_product_x(alpha, A_center(3), beta, B_center(3), fact_p, p, P_center(3)) overlap_z *= fact_p overlap = overlap_x * overlap_y * overlap_z @@ -183,7 +185,7 @@ subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, powe double precision :: F_integral double precision, allocatable :: P_new(:,:,:), P_center(:,:), fact_p(:) - ldp = maxval( power_A(1:3) + power_B(1:3) ) + ldp = maxval(power_A(1:3) + power_B(1:3)) allocate(P_new(n_points,0:ldp,3), P_center(n_points,3), fact_p(n_points))