Merge pull request #294 from QuantumPackage/io_ccsd

Improved  I/O in CCSD

src/ao_two_e_ints/cholesky.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
  ! Number of Cholesky vectors in AO basis
  END_DOC
 
- cholesky_ao_num_guess = ao_num*ao_num / 2
+ cholesky_ao_num_guess = ao_num*ao_num 
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer, cholesky_ao_num ]
@@ -44,19 +44,12 @@ END_PROVIDER
    do m=0,9
      do l=1+m,ao_num,10
        !$OMP DO SCHEDULE(dynamic)
-       do j=1,l
+       do j=1,ao_num
          do k=1,ao_num
-           do i=1,min(k,j)
+           do i=1,ao_num
              if (ao_two_e_integral_zero(i,j,k,l)) cycle
              integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map)
              ao_integrals(i,k,j,l) = integral
-             ao_integrals(k,i,j,l) = integral
-             ao_integrals(i,k,l,j) = integral
-             ao_integrals(k,i,l,j) = integral
-             ao_integrals(j,l,i,k) = integral
-             ao_integrals(j,l,k,i) = integral
-             ao_integrals(l,j,i,k) = integral
-             ao_integrals(l,j,k,i) = integral
            enddo
          enddo
        enddo

src/ccsd/EZFIO.cfg

@@ -0,0 +1,11 @@
+[energy]
+type: double precision
+doc: CCSD energy
+interface: ezfio
+
+[energy_t]
+type: double precision
+doc: CCSD(T) energy
+interface: ezfio
+
+

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -135,8 +135,11 @@ subroutine run_ccsd_space_orb
   write(*,'(A15,1pE10.2,A3)')' Conv        = ', max_r
   print*,''
 
-  call write_t1(nO,nV,t1)
-  call write_t2(nO,nV,t2)
+  if (write_amplitudes) then
+    call write_t1(nO,nV,t1)
+    call write_t2(nO,nV,t2)
+    call ezfio_set_utils_cc_io_amplitudes('Read')
+  endif
 
   ! Deallocation
   if (cc_update_method == 'diis') then
@@ -147,6 +150,7 @@ subroutine run_ccsd_space_orb
 
   ! CCSD(T)
   double precision :: e_t
+  e_t = 0.d0
 
   if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then
 
@@ -182,8 +186,7 @@ subroutine run_ccsd_space_orb
     print*,''
   endif
 
-  print*,'Reference determinant:'
-  call print_det(det,N_int)
+  call save_energy(uncorr_energy + energy, e_t)
 
   deallocate(t1,t2)
 

src/ccsd/ccsd_spin_orb_sub.irp.f

@@ -269,8 +269,11 @@ subroutine run_ccsd_spin_orb
   write(*,'(A15,1pE10.2,A3)')' Conv        = ', max_r
   print*,''
 
-  call write_t1(nO,nV,t1)
-  call write_t2(nO,nV,t2)
+  if (write_amplitudes) then
+    call write_t1(nO,nV,t1)
+    call write_t2(nO,nV,t2)
+    call ezfio_set_utils_cc_io_amplitudes('Read')
+  endif
 
   ! Deallocate
   if (cc_update_method == 'diis') then
@@ -284,8 +287,9 @@ subroutine run_ccsd_spin_orb
   deallocate(v_ovoo,v_oovo)
   deallocate(v_ovvo,v_ovov,v_oovv)
   
+  double precision :: t_corr
+  t_corr = 0.d0
   if (cc_par_t .and. elec_alpha_num  +elec_beta_num > 2) then
-    double precision :: t_corr
     print*,'CCSD(T) calculation...'
     call wall_time(ta)
     !allocate(v_vvvo(nV,nV,nV,nO))
@@ -307,8 +311,8 @@ subroutine run_ccsd_spin_orb
     write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + t_corr, ' Ha'
     print*,''
   endif
-  print*,'Reference determinant:'
-  call print_det(det,N_int)
+
+  call save_energy(uncorr_energy + energy, t_corr)
   
   deallocate(f_oo,f_ov,f_vv,f_o,f_v)
   deallocate(v_ooov,v_vvoo,t1,t2)

src/ccsd/save_energy.irp.f

@@ -0,0 +1,13 @@
+subroutine save_energy(E,ET)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E, ET
+  call ezfio_set_ccsd_energy(E)
+  if (ET /= 0.d0) then
+    call ezfio_set_ccsd_energy_t(E+ET)
+  endif
+end
+
+

src/mo_two_e_ints/cholesky.irp.f

@@ -27,6 +27,8 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_ao_num, mo_num,
  double precision, allocatable :: buffer(:,:)
 
  print *, 'AO->MO Transformation of Cholesky vectors  .'
+
+ call set_multiple_levels_omp(.False.)
  !$OMP PARALLEL PRIVATE(i,j,k,buffer)
  allocate(buffer(mo_num,mo_num))
  !$OMP DO SCHEDULE(static)

src/utils/linear_algebra.irp.f

@@ -1831,7 +1831,7 @@ double precision, intent(in)     :: tol
 
 integer, dimension(:), allocatable          :: piv
 double precision, dimension(:), allocatable :: work
-character, parameter :: uplo = "U"
+character, parameter :: uplo = 'L'
 integer :: LDA
 integer :: info
 integer :: k, l, rank0
@@ -1848,14 +1848,14 @@ if (rank > rank0) then
 end if
 
 do k = 1, ndim
-  A(k+1:ndim, k) = 0.00D+0
+  A(k,k+1:ndim) = 0.00D+0
 end do
 ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
 ! to do the swapping in-place
 U(:,:) = 0.00D+0
 do k = 1, ndim
   l = piv(k)
-  U(l, 1:rank) = A(1:rank, k)
+  U(l, 1:rank) = A(k,1:rank)
 end do
 
 end subroutine pivoted_cholesky

src/utils_cc/EZFIO.cfg

@@ -46,17 +46,11 @@ doc: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation t
 interface: ezfio,ocaml,provider
 default: MP
 
-[cc_write_t1]
-type: logical
-doc: If true, it will write on disk the T1 amplitudes at the end of the calculation.
-interface: ezfio,ocaml,provider
-default: False
-
-[cc_write_t2]
-type: logical
-doc: If true, it will write on disk the T2 amplitudes at the end of the calculation.
-interface: ezfio,ocaml,provider
-default: False
+[io_amplitudes]
+type: Disk_access
+doc: Read/Write |CCSD| amplitudes from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
 
 [cc_par_t]
 type: logical

src/utils_cc/guess_t.irp.f

@@ -91,16 +91,17 @@ subroutine write_t1(nO,nV,t1)
   double precision, intent(in) :: t1(nO, nV)
 
   ! internal
-  integer                      :: i,a
+  integer                      :: i,a, iunit
+  integer, external :: getunitandopen
 
-  if (cc_write_t1) then
-    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+  if (write_amplitudes) then
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','w')
     do a = 1, nV
       do i = 1, nO
-         write(11,'(F20.12)') t1(i,a)
+         write(iunit,'(F20.12)') t1(i,a)
       enddo
     enddo
-    close(11)
+    close(iunit)
   endif
   
 end
@@ -120,20 +121,21 @@ subroutine write_t2(nO,nV,t2)
   double precision, intent(in) :: t2(nO, nO, nV, nV)
 
   ! internal
-  integer                      :: i,j,a,b
+  integer                      :: i,j,a,b, iunit
+  integer, external :: getunitandopen
 
-  if (cc_write_t2) then
-    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+  if (write_amplitudes) then
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','w')
     do b = 1, nV
       do a = 1, nV
         do j = 1, nO
           do i = 1, nO
-             write(11,'(F20.12)') t2(i,j,a,b)
+             write(iunit,'(F20.12)') t2(i,j,a,b)
           enddo
         enddo
       enddo
     enddo
-    close(11)
+    close(iunit)
   endif
   
 end
@@ -153,23 +155,19 @@ subroutine read_t1(nO,nV,t1)
   double precision, intent(out) :: t1(nO, nV)
 
   ! internal
-  integer                       :: i,a
+  integer                       :: i,a, iunit
   logical                       :: ok
+  integer, external :: getunitandopen
 
-  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
-  if (.not. ok) then
-     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
-     print*, 'Do a first calculation with cc_write_t1 = True'
-     print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
-     call abort
-  endif
-  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
-  do a = 1, nV
-    do i = 1, nO
-       read(11,'(F20.12)') t1(i,a)
+  if (read_amplitudes) then
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','r')
+    do a = 1, nV
+      do i = 1, nO
+         read(iunit,'(F20.12)') t1(i,a)
+      enddo
     enddo
-  enddo
-  close(11)
+    close(iunit)
+  endif
   
 end
 
@@ -188,26 +186,23 @@ subroutine read_t2(nO,nV,t2)
   double precision, intent(out) :: t2(nO, nO, nV, nV)
 
   ! internal
-  integer                       :: i,j,a,b
+  integer                       :: i,j,a,b, iunit
   logical                       :: ok
 
-  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
-  if (.not. ok) then
-     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
-     print*, 'Do a first calculation with cc_write_t2 = True'
-     print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
-     call abort
-  endif
-  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
-  do b = 1, nV
-    do a = 1, nV
-      do j = 1, nO
-        do i = 1, nO
-           read(11,'(F20.12)') t2(i,j,a,b)
+  integer, external :: getunitandopen
+
+  if (read_amplitudes) then
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','r')
+    do b = 1, nV
+      do a = 1, nV
+        do j = 1, nO
+          do i = 1, nO
+             read(iunit,'(F20.12)') t2(i,j,a,b)
+          enddo
         enddo
       enddo
     enddo
-  enddo
-  close(11)
+    close(iunit)
+  endif
   
 end