Print QCint in binary files

src/modules/quick_calculated_module.f90

@@ -167,6 +167,11 @@ module quick_calculated_module
       ! Mulliken charge and Lowdin charge
       double precision,dimension(:),allocatable :: Mulliken,Lowdin
 
+      ! ERI index for 8 fold symmetry
+      integer, dimension(:,:,:,:), allocatable :: iarray2
+
+      ! ERI in 8 fold symmetry
+      double precision,dimension(:),allocatable :: aoint2e 
 
    end type quick_qm_struct_type
 
@@ -222,6 +227,10 @@ subroutine allocate_quick_qm_struct(self)
       integer nelec
       integer idimA
       integer nelecb
+      !For QC output
+      integer i,j,k,l
+      integer ij,kl,ijkl
+      integer npair,ns8
 
       type (quick_qm_struct_type) self
 
@@ -247,6 +256,39 @@ subroutine allocate_quick_qm_struct(self)
       if(.not. allocated(self%Mulliken)) allocate(self%Mulliken(natom))
       if(.not. allocated(self%Lowdin)) allocate(self%Lowdin(natom))
 
+      if (quick_method%QCint) then
+         if(.not. allocated(self%iarray2)) allocate(self%iarray2(nbasis,nbasis,nbasis,nbasis))
+         self%iarray2=0.d0
+         ijkl=0
+         ij=0
+         do i=1,nbasis
+            do j=1,i
+               ij=ij+1
+               kl=0
+               do k=1,i
+                  do l=1,k
+                     kl=kl+1
+                     if (ij .ge. kl) then
+                        ijkl=ijkl+1
+                        self%iarray2(i,j,k,l)=ijkl
+                        self%iarray2(i,j,l,k)=ijkl
+                        self%iarray2(j,i,l,k)=ijkl
+                        self%iarray2(j,i,k,l)=ijkl
+                        self%iarray2(k,l,i,j)=ijkl
+                        self%iarray2(k,l,j,i)=ijkl
+                        self%iarray2(l,k,j,i)=ijkl
+                        self%iarray2(l,k,i,j)=ijkl
+                     endif
+                  enddo
+               enddo
+            enddo
+         enddo
+
+         npair=(nbasis*(nbasis+1))/2
+         ns8=(ns8*(ns8+1))/2
+         if(.not. allocated(self%aoint2e)) allocate(self%aoint2e(ns8))
+      endif
+
       ! if 1st order derivation, which is gradient calculation is requested
       if (quick_method%grad) then
          if(.not. allocated(self%gradient)) allocate(self%gradient(3*natom))
@@ -410,6 +452,11 @@ subroutine deallocate_quick_qm_struct(self)
       if (allocated(self%Mulliken)) deallocate(self%Mulliken)
       if (allocated(self%Lowdin)) deallocate(self%Lowdin)
 
+      if (quick_method%QCint) then
+         if (allocated(self%iarray2)) deallocate(self%iarray2)
+         if (allocated(self%aoint2e)) deallocate(self%aoint2e)
+      endif
+
       ! if 1st order derivation, which is gradient calculation is requested
       if (quick_method%grad) then
          if (allocated(self%gradient)) deallocate(self%gradient)
@@ -568,6 +615,11 @@ subroutine init_quick_qm_struct(self)
 
       call zeroMatrix(self%denseOld,nbasis)
 
+      ! For QC output
+      if (quick_method%QCint) then
+         self%aoint2e=0.0d0
+      endif
+
       ! if 1st order derivation, which is gradient calculation is requested
       if (quick_method%grad) then
          call zeroVec(self%gradient,3*natom)

src/modules/quick_method_module.f90

@@ -79,7 +79,7 @@ module quick_method_module
 
         ! this part is Div&Con options
         logical :: BEoff = .false.
-        logical :: Qint = .false.
+        logical :: QCint = .false.
         logical :: OWNfrag = .false.
 
         ! this is DFT grid
@@ -253,7 +253,7 @@ subroutine broadcast_quick_method(self, ierr)
             call MPI_BCAST(self%DIVCON,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%DCMP2only,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%BEoff,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
-            call MPI_BCAST(self%Qint,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%QCint,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%OWNfrag,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%DNCRB,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%DNCRB2,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
@@ -503,7 +503,7 @@ subroutine print_quick_method(self,io,ierr)
 
             ! additional DIVCON options
             if (self%BEoff)  write(io,'(" NO ELIMINATION OF EMBEDDED SUBSYSTEM IN DIVCON ")')
-            if (self%Qint)   write(io,'(" PRINT DIVCON FRAGMENT INTEGRALS IN QISKIT FORMAT ")') 
+            if (self%QCint)   write(io,'(" PRINT DIVCON FRAGMENT INTEGRALS IN QISKIT FORMAT ")') 
 
             ! computing cycles
             write(io,'(" MAX SCF CYCLES = ",i6)') self%iscf
@@ -672,7 +672,7 @@ subroutine read_quick_method(self,keywd,ierr)
             if (index(keyWD,'EXTCHARGES').ne.0) self%EXTCHARGES=.true.
             if (index(keyWD,'FORCE').ne.0)      self%grad=.true.
             if (index(keyWD,'BEOFF').ne.0)      self%BEoff=.true.
-            if (index(keyWD,'QINT').ne.0)       self%Qint=.true.
+            if (index(keyWD,'QCINT').ne.0)       self%QCint=.true.
 
             ! The BEoff keyword is redundant when OWNFRAG is used.
             ! OWNFRAG sets radii of both buffer1 and buffer2 to zero,
@@ -928,7 +928,7 @@ subroutine init_quick_method(self,ierr)
             self%DIVCON = .false.      ! Div&Con
             self%DCMP2only = .false.   ! Canonical HF and Div&Con for MP2 only
             self%BEoff = .false.       ! Turns off bEliminate in Div&Con
-            self%Qint = .false.        ! Prints 1e and 2e integrals in Qiskit format
+            self%QCint = .false.        ! Prints 1e and 2e integrals in Qiskit format
             self%OWNfrag = .false.     ! User defined fragmentation based on CORE.DNC and BUFFER.DNC files
 
             self%ifragbasis = 1        ! =2.residue basis,=1.atom basis(DEFUALT),=3 non-h atom basis

src/modules/quick_scf_module.f90

@@ -113,8 +113,12 @@ subroutine scf(ierr)
 
      logical :: done
      integer, intent(inout) :: ierr
-     integer :: jscf
+     integer :: jscf,ierr1
      done=.false.
+
+     if (quick_method%QCint) then
+        call quick_open(iIntFile,intFileName,'U','U','A',.false.,ierr1)
+     end if
 
      !-----------------------------------------------------------------
      ! The purpose of this subroutine is to perform scf cycles.  At this
@@ -206,6 +210,9 @@ subroutine electdiis(jscf,ierr)
 
      double precision :: c_coords(3),c_zeta,c_chg
 
+     ! For QC output
+     integer :: npair,ns8,ijkl, ijkl2
+
      !---------------------------------------------------------------------------
      ! The purpose of this subroutine is to utilize Pulay's accelerated
      ! scf convergence as detailed in J. Comp. Chem, Vol 3, #4, pg 566-60, 1982.
@@ -356,6 +363,21 @@ subroutine electdiis(jscf,ierr)
 
         call scf_operator(deltaO)
 
+        if (quick_method%QCint) then
+        if (jscf .eq. 2) then
+           write(iintfile)'1E INTEGRALS'
+           do ijkl=1, nbasis
+              write(iintfile)(quick_qm_struct%oneElecO(ijkl,ijkl2),ijkl2=1,nbasis)
+           enddo
+           write(iintfile)'2E INTEGRALS'
+           npair = (nbasis*(nbasis+1))/2
+           ns8 = (npair*(npair+1))/2
+           do ijkl=1,ns8
+              write(iintfile)quick_qm_struct%aoint2e(ijkl)
+           enddo
+        endif
+        endif
+
         quick_qm_struct%denseOld(:,:) = quick_qm_struct%dense(:,:)
 
         if (quick_method%debug)  call debug_SCF(jscf)

src/subs/hrr.f90

@@ -22,6 +22,9 @@ subroutine hrrwhole
    common /xiaostore/store
    common /hrrstore/II,JJ,KK,LL,NBI1,NBI2,NBJ1,NBJ2,NBK1,NBK2,NBL1,NBL2
 
+   !For QC output
+   integer ijkl
+
    select case (IJKLtype)
 
    case (0,10,1000,1010)
@@ -282,7 +285,13 @@ subroutine hrrwhole
 
    end select
    111 continue
-!  write(*,*) IJKLtype,mpirank, iii,jjj,kkk,lll,Y
+
+   if (quick_method%QCint) then
+      ijkl=quick_qm_struct%iarray2(iii,jjj,kkk,lll)
+      quick_qm_struct%aoint2e(ijkl)=Y
+print*,IJKL,quick_qm_struct%aoint2e(ijkl)
+   endif
+
 end subroutine hrrwhole