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/external/qp2-dependencies b/external/qp2-dependencies index 242151e0..90ee61f5 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 242151e03d1d6bf042387226431d82d35845686a +Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0 diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f index 4dd87a60..f01ed5ba 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -515,4 +515,3 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p END_PROVIDER ! --- - diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 8d4fff57..262301e0 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -106,3 +106,26 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_ao +[n_pts_charge] +type: integer +doc: Number of point charges to be added to the potential +interface: ezfio +default: 0 + +[pts_charge_z] +type: double precision +doc: Charge associated to each point charge +interface: ezfio +size: (ao_one_e_ints.n_pts_charge) + +[pts_charge_coord] +type: double precision +doc: Coordinate of each point charge. +interface: ezfio +size: (ao_one_e_ints.n_pts_charge,3) + +[point_charges] +type: logical +doc: If |true|, point charges (see ao_one_e_ints/write_pt_charges.py) are added to the one-electron potential +interface: ezfio,provider,ocaml +default: False diff --git a/src/ao_one_e_ints/point_charges.irp.f b/src/ao_one_e_ints/point_charges.irp.f new file mode 100644 index 00000000..c038458d --- /dev/null +++ b/src/ao_one_e_ints/point_charges.irp.f @@ -0,0 +1,272 @@ + +! --- + + +BEGIN_PROVIDER [ integer, n_pts_charge ] + implicit none + BEGIN_DOC +! Number of point charges to be added to the potential + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + + call ezfio_has_ao_one_e_ints_n_pts_charge(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..' + call ezfio_get_ao_one_e_ints_n_pts_charge(n_pts_charge) + else + print *, 'ao_one_e_ints/n_pts_charge not found in EZFIO file' + stop 1 + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( n_pts_charge, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read n_pts_charge with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pts_charge_z, (n_pts_charge) ] + + BEGIN_DOC + ! Charge associated to each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_ao_one_e_ints_pts_charge_z(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_z ] <<<<< ..' + call ezfio_get_ao_one_e_ints_pts_charge_z(pts_charge_z) + IRP_IF MPI + call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + endif + + else + + integer :: i + do i = 1, n_pts_charge + pts_charge_z(i) = 0.d0 + enddo + + endif + print*,'Point charges ' + do i = 1, n_pts_charge + print*,'i,pts_charge_z(i)',i,pts_charge_z(i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, pts_charge_coord, (n_pts_charge,3) ] + + BEGIN_DOC + ! Coordinates of each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_ao_one_e_ints_pts_charge_coord(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + double precision, allocatable :: buffer(:,:) + allocate (buffer(n_pts_charge,3)) + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_coord ] <<<<< ..' + call ezfio_get_ao_one_e_ints_pts_charge_coord(buffer) + integer :: i,j + do i=1,3 + do j=1,n_pts_charge + pts_charge_coord(j,i) = buffer(j,i) + enddo + enddo + deallocate(buffer) + IRP_IF MPI + call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + endif + + else + + do i = 1, n_pts_charge + pts_charge_coord(i,:) = 0.d0 + enddo + + endif + print*,'Coordinates for the point charges ' + do i = 1, n_pts_charge + write(*,'(I3,X,3(F16.8,X))'),i,pts_charge_coord(i,1:3) + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)] + + BEGIN_DOC + ! Point charge-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + ! + ! These integrals also contain the pseudopotential integrals. + END_DOC + + implicit none + integer :: num_A, num_B, power_A(3), power_B(3) + integer :: i, j, k, l, n_pt_in, m + double precision :: alpha, beta + double precision :: A_center(3),B_center(3),C_center(3) + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_integrals_pt_chrg = 0.d0 + +! if (read_ao_integrals_pt_chrg) then +! +! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals read from disk' +! +! else + +! if(use_cosgtos) then +! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos +! +! do j = 1, ao_num +! do i = 1, ao_num +! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j) +! enddo +! enddo +! +! else + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,c1,n_pt_in) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,& + !$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_num + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + + double precision :: c, c1 + c = 0.d0 + + do k = 1, n_pts_charge + double precision :: Z + Z = pts_charge_z(k) + + C_center(1:3) = pts_charge_coord(k,1:3) + + c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & + , alpha, beta, C_center, n_pt_in ) + + c = c + Z * c1 + + enddo + ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) & + + ao_coef_normalized_ordered_transp(l,j) & + * ao_coef_normalized_ordered_transp(m,i) * c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + +! endif + + +! IF(do_pseudo) THEN +! ao_integrals_pt_chrg += ao_pseudo_integrals +! ENDIF + +! endif + + +! if (write_ao_integrals_pt_chrg) then +! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals written to disk' +! endif + +END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 928053ad..20e299af 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -104,6 +104,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] IF(do_pseudo) THEN ao_integrals_n_e += ao_pseudo_integrals ENDIF + IF(point_charges) THEN + ao_integrals_n_e += ao_integrals_pt_chrg + ENDIF endif diff --git a/src/ao_one_e_ints/write_pt_charges.py b/src/ao_one_e_ints/write_pt_charges.py new file mode 100755 index 00000000..d4b6d251 --- /dev/null +++ b/src/ao_one_e_ints/write_pt_charges.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python +import os +import sys + +# First argument is the EZFIO file +# It reads a file EZFIO_point_charges.xyz written in this way: +# charge x y z (Angstrom) +# for all charges + + +def zip_in_ezfio(ezfio,tmp): + tmpzip=tmp+".gz" + cmdzip="gzip -c "+tmp+" > "+tmpzip + os.system(cmdzip) + os.system("rm "+tmp) + cmdmv="mv "+tmpzip+" "+EZFIO+"/ao_one_e_ints/"+tmpzip + os.system(cmdmv) + +def mv_in_ezfio(ezfio,tmp): + cmdmv="mv "+tmp+" "+EZFIO+"/ao_one_e_ints/"+tmp + os.system(cmdmv) + + +# Getting the EZFIO +EZFIO=sys.argv[1] +EZFIO=EZFIO.replace("/", "") +print(EZFIO) + +# Reading the point charges and convert the Angstrom geometry in Bohr for QP +f = open(EZFIO+'_point_charges.xyz','r') +lines = f.readlines() +convert_angs_to_bohr=1.8897259885789233 + +n_charges=0 +coord_x=[] +coord_y=[] +coord_z=[] +charges=[] +for line in lines: + data = line.split() + if(len(data)>0): + n_charges += 1 + charges.append(str(data[0])) + coord_x.append(str(convert_angs_to_bohr*float(data[1]))) + coord_y.append(str(convert_angs_to_bohr*float(data[2]))) + coord_z.append(str(convert_angs_to_bohr*float(data[3]))) + +# Write the file containing the number of charges and set in EZFIO folder +tmp="n_pts_charge" +fncharges = open(tmp,'w') +fncharges.write(" "+str(n_charges)+'\n') +fncharges.close() +mv_in_ezfio(EZFIO,tmp) + +# Write the file containing the charges and set in EZFIO folder +tmp="pts_charge_z" +fcharges = open(tmp,'w') +fcharges.write(" 1\n") +fcharges.write(" "+str(n_charges)+'\n') +for i in range(n_charges): + fcharges.write(charges[i]+'\n') +fcharges.close() +zip_in_ezfio(EZFIO,tmp) + +# Write the file containing the charge coordinates and set in EZFIO folder +tmp="pts_charge_coord" +fcoord = open(tmp,'w') +fcoord.write(" 2\n") +fcoord.write(" "+str(n_charges)+' 3\n') +#fcoord.write(" "+' 3 '+str(n_charges)+' \n') +for i in range(n_charges): + fcoord.write(' '+coord_x[i]+'\n') +for i in range(n_charges): + fcoord.write(' '+coord_y[i]+'\n') +for i in range(n_charges): + fcoord.write(' '+coord_z[i]+'\n') +fcoord.close() +zip_in_ezfio(EZFIO,tmp) + diff --git a/src/bi_ortho_mos/bi_density.irp.f b/src/bi_ortho_mos/bi_density.irp.f index 90fe9634..2dad9485 100644 --- a/src/bi_ortho_mos/bi_density.irp.f +++ b/src/bi_ortho_mos/bi_density.irp.f @@ -15,6 +15,7 @@ BEGIN_PROVIDER [double precision, TCSCF_bi_ort_dm_ao_alpha, (ao_num, ao_num) ] 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 @@ -35,6 +36,7 @@ BEGIN_PROVIDER [ double precision, TCSCF_bi_ort_dm_ao_beta, (ao_num, ao_num) ] 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/determinants/filter_connected.irp.f b/src/determinants/filter_connected.irp.f index 6110eb89..2c9d7a49 100644 --- a/src/determinants/filter_connected.irp.f +++ b/src/determinants/filter_connected.irp.f @@ -96,7 +96,6 @@ subroutine filter_not_connected(key1,key2,Nint,sze,idx) idx(0) = l-1 end - subroutine filter_connected(key1,key2,Nint,sze,idx) use bitmasks implicit none diff --git a/src/determinants/sparse_mat.irp.f b/src/determinants/sparse_mat.irp.f new file mode 100644 index 00000000..889bbeba --- /dev/null +++ b/src/determinants/sparse_mat.irp.f @@ -0,0 +1,164 @@ + use bitmasks + +subroutine filter_connected_array(key1,key2,ld,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC + integer, intent(in) :: Nint, ld,sze + integer(bit_kind), intent(in) :: key1(Nint,2,ld) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) +! print*,degree_x2 + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DIR$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +! print*,'idx(0) = ',idx(0) +end + + + BEGIN_PROVIDER [ integer, n_sparse_mat] +&BEGIN_PROVIDER [ integer, n_connected_per_det, (N_det)] +&BEGIN_PROVIDER [ integer, n_max_connected_per_det] + implicit none + BEGIN_DOC +! n_sparse_mat = total number of connections in the CI matrix +! +! n_connected_per_det(i) = number of connected determinants to the determinant psi_det(1,1,i) +! +! n_max_connected_per_det = maximum number of connected determinants + END_DOC + integer, allocatable :: idx(:) + allocate(idx(0:N_det)) + integer :: i + n_sparse_mat = 0 + do i = 1, N_det + call filter_connected_array(psi_det_sorted,psi_det_sorted(1,1,i),psi_det_size,N_int,N_det,idx) + n_connected_per_det(i) = idx(0) + n_sparse_mat += idx(0) + enddo + n_max_connected_per_det = maxval(n_connected_per_det) +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), connected_det_per_det, (N_int,2,n_max_connected_per_det,N_det)] +&BEGIN_PROVIDER [ integer(bit_kind), list_connected_det_per_det, (n_max_connected_per_det,N_det)] + implicit none + BEGIN_DOC +! connected_det_per_det(:,:,j,i) = jth connected determinant to the determinant psi_det(:,:,i) +! +! list_connected_det_per_det(j,i) = index of jth determinant in psi_det which is connected to psi_det(:,:,i) + END_DOC + integer, allocatable :: idx(:) + allocate(idx(0:N_det)) + integer :: i,j + do i = 1, N_det + call filter_connected_array(psi_det_sorted,psi_det_sorted(1,1,i),psi_det_size,N_int,N_det,idx) + do j = 1, idx(0) + connected_det_per_det(1:N_int,1:2,j,i) = psi_det_sorted(1:N_int,1:2,idx(j)) + list_connected_det_per_det(j,i) = idx(j) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, sparse_h_mat, (n_max_connected_per_det, N_det)] + implicit none + BEGIN_DOC +! sparse matrix format +! +! sparse_h_mat(j,i) = matrix element between the jth connected determinant and psi_det(:,:,i) + END_DOC + integer :: i,j + double precision :: hij + do i = 1, N_det + do j = 1, n_connected_per_det(i) + call i_H_j(psi_det(1,1,i),connected_det_per_det(1,1,j,i),N_int,hij) + sparse_h_mat(j,i) = hij + enddo + enddo + +END_PROVIDER + diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index 81a8fe71..ff3d11f3 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -299,6 +299,7 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_ END_PROVIDER +! --- ! --- BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ] 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 9aef436f..854789bd 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -340,6 +340,7 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, 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/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/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index eb8fa8be..fabc3d14 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -136,36 +136,12 @@ doc: nb of Gaussians used to fit Jastrow fcts interface: ezfio,provider,ocaml default: 20 -[max_dim_diis_tcscf] -type: integer -doc: Maximum size of the DIIS extrapolation procedure -interface: ezfio,provider,ocaml -default: 15 - -[threshold_diis_tcscf] -type: Threshold -doc: Threshold on the convergence of the DIIS error vector during a TCSCF calculation. If 0. is chosen, the square root of thresh_tcscf will be used. -interface: ezfio,provider,ocaml -default: 0. - -[level_shift_tcscf] -type: Positive_float -doc: Energy shift on the virtual MOs to improve TCSCF convergence -interface: ezfio,provider,ocaml -default: 0. - [tcscf_algorithm] type: character*(32) doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] interface: ezfio,provider,ocaml default: Simple -[im_thresh_tcscf] -type: Threshold -doc: Thresholds on the Imag part of energy -interface: ezfio,provider,ocaml -default: 1.e-7 - [test_cycle_tc] type: logical doc: If |true|, the integrals of the three-body jastrow are computed with cycles @@ -184,3 +160,27 @@ doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph interface: ezfio,provider,ocaml default: 1.e-6 +[max_dim_diis_tcscf] +type: integer +doc: Maximum size of the DIIS extrapolation procedure +interface: ezfio,provider,ocaml +default: 15 + +[threshold_diis_tcscf] +type: Threshold +doc: Threshold on the convergence of the DIIS error vector during a TCSCF calculation. If 0. is chosen, the square root of thresh_tcscf will be used. +interface: ezfio,provider,ocaml +default: 0. + +[level_shift_tcscf] +type: Positive_float +doc: Energy shift on the virtual MOs to improve TCSCF convergence +interface: ezfio,provider,ocaml +default: 0. + +[im_thresh_tcscf] +type: Threshold +doc: Thresholds on the Imag part of energy +interface: ezfio,provider,ocaml +default: 1.e-7 + diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 7403049c..6796666d 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -41,6 +41,15 @@ !! rho_b(l,j) * < l k| T | i j> !two_e_tc_non_hermit_integral_seq_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) + !! 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_seq_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) ! rho(l,j) * < k l| T | i j> diff --git a/src/tc_scf/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f new file mode 100644 index 00000000..0312df5f --- /dev/null +++ b/src/tc_scf/rh_tcscf.irp.f @@ -0,0 +1,336 @@ +! --- + +subroutine rh_tcscf() + + BEGIN_DOC + ! + ! Roothaan-Hall algorithm for TC-SCF calculation + ! + END_DOC + + implicit none + + integer :: i, j + integer :: iteration_TCSCF, dim_DIIS, index_dim_DIIS + double precision :: energy_TCSCF, energy_TCSCF_1e, energy_TCSCF_2e, energy_TCSCF_3e, gradie_TCSCF + double precision :: energy_TCSCF_previous, delta_energy_TCSCF + double precision :: gradie_TCSCF_previous, delta_gradie_TCSCF + double precision :: max_error_DIIS_TCSCF + double precision :: level_shift_save + double precision :: delta_energy_tmp, delta_gradie_tmp + double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:) + double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:) + + logical, external :: qp_stop + + + !PROVIDE ao_md5 mo_occ + PROVIDE level_shift_TCSCF + + allocate( mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num) & + , F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), e_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF) ) + + F_DIIS = 0.d0 + e_DIIS = 0.d0 + mo_l_coef_save = 0.d0 + mo_r_coef_save = 0.d0 + + call write_time(6) + + ! --- + ! Initialize energies and density matrices + + energy_TCSCF_previous = TC_HF_energy + energy_TCSCF_1e = TC_HF_one_e_energy + energy_TCSCF_2e = TC_HF_two_e_energy + energy_TCSCF_3e = 0.d0 + if(three_body_h_tc) then + energy_TCSCF_3e = diag_three_elem_hf + endif + gradie_TCSCF_previous = grad_non_hermit + delta_energy_TCSCF = 1.d0 + delta_gradie_TCSCF = 1.d0 + iteration_TCSCF = 0 + dim_DIIS = 0 + max_error_DIIS_TCSCF = 1.d0 + + ! --- + + ! Start of main SCF loop + + PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot + + do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. & + !(dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. & + (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF)) ) + + iteration_TCSCF += 1 + if(iteration_TCSCF > n_it_TCSCF_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + stop + endif + + dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) + + ! --- + + if((tcscf_algorithm == 'DIIS') .and. (dabs(delta_energy_TCSCF) > 1.d-6)) then + + ! store Fock and error matrices at each iteration + index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1 + do j = 1, ao_num + do i = 1, ao_num + F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j) + e_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao(i,j) + enddo + enddo + + 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 + !TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta + + 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) ) + 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 + + ! --- + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + 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 + + ! --- + + ! calculate error vectors + max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo)) + + ! --- + + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous + + ! --- + + 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) + + if(level_shift_TCSCF <= .1d0) then + level_shift_TCSCF = 1.d0 + else + level_shift_TCSCF = level_shift_TCSCF * 3.0d0 + endif + TOUCH mo_l_coef mo_r_coef level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + 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_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 + SOFT_TOUCH level_shift_TCSCF + exit + endif + + dim_DIIS = 0 + enddo +! print *, ' very big step : ', delta_energy_tmp +! print *, ' TC level shift = ', level_shift_TCSCF + + ! --- + + level_shift_TCSCF = 0.d0 + !level_shift_TCSCF = level_shift_TCSCF * 0.5d0 + SOFT_TOUCH level_shift_TCSCF + + gradie_TCSCF = grad_non_hermit + energy_TCSCF = TC_HF_energy + energy_TCSCF_1e = TC_HF_one_e_energy + energy_TCSCF_2e = TC_HF_two_e_energy + energy_TCSCF_3e = 0.d0 + if(three_body_h_tc) then + energy_TCSCF_3e = diag_three_elem_hf + endif + delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous + delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous + + energy_TCSCF_previous = energy_TCSCF + gradie_TCSCF_previous = gradie_TCSCF + + + level_shift_save = level_shift_TCSCF + mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) + mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) + + + print *, ' iteration = ', iteration_TCSCF + print *, ' total TC energy = ', energy_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| = ', 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 + print *, ' TC DIIS dim = ', dim_DIIS + print *, ' TC level shift = ', level_shift_TCSCF + print *, ' ' + + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + + if(qp_stop()) exit + enddo + + ! --- + + 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) + +end + +! --- + +subroutine extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, F_ao, size_F_ao, iteration_TCSCF, dim_DIIS) + + BEGIN_DOC + ! + ! Compute the extrapolated Fock matrix using the DIIS procedure + ! + ! e = \sum_i c_i e_i and \sum_i c_i = 1 + ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1) + ! + END_DOC + + implicit none + + integer, intent(in) :: iteration_TCSCF, size_F_ao + integer, intent(inout) :: dim_DIIS + double precision, intent(in) :: F_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(in) :: e_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(inout) :: F_ao(size_F_ao,ao_num) + + double precision, allocatable :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:) + + integer :: i, j, k, l, i_DIIS, j_DIIS + integer :: lwork + double precision :: rcond, ferr, berr + integer, allocatable :: iwork(:) + double precision, allocatable :: scratch(:,:) + + 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(ao_num,ao_num) ) + + ! Compute the matrices B and X + B_matrix_DIIS(:,:) = 0.d0 + do j = 1, dim_DIIS + j_DIIS = min(dim_DIIS, mod(iteration_TCSCF-j, max_dim_DIIS_TCSCF)+1) + + do i = 1, dim_DIIS + i_DIIS = min(dim_DIIS, mod(iteration_TCSCF-i, max_dim_DIIS_TCSCF)+1) + + ! Compute product of two errors vectors + do l = 1, ao_num + do k = 1, ao_num + B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + e_DIIS(k,l,i_DIIS) * e_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 + integer :: info + double precision :: anorm + 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 + + 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 TC-DIIS' + endif + + ! Compute extrapolated Fock matrix + + !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) + do j = 1, ao_num + do i = 1, ao_num + F_ao(i,j) = 0.d0 + enddo + do k = 1, dim_DIIS + if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle + do i = 1,ao_num + ! FPE here + F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + +! --- + diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 187750ff..deaf8d82 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -73,4 +73,3 @@ subroutine create_guess() end subroutine create_guess ! --- - diff --git a/src/tc_scf/tc_scf_dm.irp.f b/src/tc_scf/tc_scf_dm.irp.f index 4750199c..90719f47 100644 --- a/src/tc_scf/tc_scf_dm.irp.f +++ b/src/tc_scf/tc_scf_dm.irp.f @@ -10,7 +10,6 @@ BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num else TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta endif - END_PROVIDER ! --- @@ -25,7 +24,6 @@ BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_nu else TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha endif - END_PROVIDER diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index 6abeddf1..a14c4126 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -62,12 +62,6 @@ subroutine test_tc_scf integer :: i ! provide int2_u_grad1u_x_j1b2_test provide x_v_ij_erf_rk_cst_mu_j1b_test -! do i = 1, ng_fit_jast -! print*,expo_gauss_1_erf_x_2(i),coef_gauss_1_erf_x_2(i) -! enddo -! provide tc_grad_square_ao_test -! provide tc_grad_and_lapl_ao_test -! provide int2_u_grad1u_x_j1b2_test ! provide x_v_ij_erf_rk_cst_mu_j1b_test ! print*,'TC_HF_energy = ',TC_HF_energy ! print*,'grad_non_hermit = ',grad_non_hermit @@ -1006,3 +1000,4 @@ end ! --- +>>>>>>> 92a4e33f8a21717cab0c0e4f8412ed6903afb04a diff --git a/src/tools/print_hmat.irp.f b/src/tools/print_hmat.irp.f new file mode 100644 index 00000000..48001e44 --- /dev/null +++ b/src/tools/print_hmat.irp.f @@ -0,0 +1,90 @@ +program print_h_mat + implicit none + BEGIN_DOC + ! program that prints out the CI matrix in sparse form + END_DOC + read_wf = .True. + touch read_wf + call print_wf_dets + call print_wf_coef + call sparse_mat + call full_mat + call test_sparse_mat +end + +subroutine print_wf_dets + implicit none + integer :: i,j + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.wf_det' + i_unit_output = getUnitAndOpen(output,'w') + write(i_unit_output,*)N_det,N_int + do i = 1, N_det + write(i_unit_output,*)psi_det_sorted(1:N_int,1,i) + write(i_unit_output,*)psi_det_sorted(1:N_int,2,i) + enddo +end + +subroutine print_wf_coef + implicit none + integer :: i,j + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.wf_coef' + i_unit_output = getUnitAndOpen(output,'w') + write(i_unit_output,*)N_det,N_states + do i = 1, N_det + write(i_unit_output,*)psi_coef_sorted(i,1:N_states) + enddo +end + +subroutine sparse_mat + implicit none + integer :: i,j + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.hmat_sparse' + i_unit_output = getUnitAndOpen(output,'w') + do i = 1, N_det + write(i_unit_output,*)i,n_connected_per_det(i) + do j =1, n_connected_per_det(i) + write(i_unit_output,*)list_connected_det_per_det(j,i),sparse_h_mat(j,i) + enddo + enddo +end + + +subroutine full_mat + implicit none + integer :: i,j + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.hmat_full' + i_unit_output = getUnitAndOpen(output,'w') + do i = 1, N_det + do j = i, N_det + write(i_unit_output,*)i,j,H_matrix_all_dets(j,i) + enddo + enddo +end + + +subroutine test_sparse_mat + implicit none + integer :: i,j + double precision, allocatable :: eigvec(:,:), eigval(:), hmat(:,:) + allocate(eigval(N_det), eigvec(N_det,N_det),hmat(N_det,N_det)) + hmat = 0.d0 + do i = 1, N_det + do j =1, n_connected_per_det(i) + hmat(list_connected_det_per_det(j,i),i) = sparse_h_mat(j,i) + enddo + enddo + call lapack_diag(eigval,eigvec,hmat,N_det,N_det) + print*,'The two energies should be the same ' + print*,'eigval(1) = ',eigval(1) + print*,'psi_energy= ',CI_electronic_energy(1) + + +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))