Single-ref amplitudes OK

devel/cc/amplitude_guess.irp.f

@@ -30,6 +30,9 @@ BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ]
     enddo
     20 continue
     close(iunit)
+ else if (cc_guess == 2) then
+   call random_number(t1_guess)
+   t1_guess *= 1.d-3
  endif
 END_PROVIDER
 
@@ -57,6 +60,8 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir
       j = 2*j-1
       a = 2*a-1 - spin_occ_num
       b = 2*b-1 - spin_occ_num
+!100 format (4(I3,X), 2(F20.10,X))
+!print 100,  i,j,a,b,t2_guess(i,j,a,b) , amplitude 
       t2_guess(i,j,a,b) = amplitude
     enddo
     10 continue
@@ -66,6 +71,7 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir
       j = 2*j  
       a = 2*a   - spin_occ_num
       b = 2*b   - spin_occ_num
+!print 100,  i,j,a,b,t2_guess(i,j,a,b) , amplitude 
       t2_guess(i,j,a,b) = amplitude
     enddo
     20 continue
@@ -75,10 +81,25 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir
       j = 2*j
       a = 2*a-1 - spin_occ_num
       b = 2*b   - spin_occ_num
+!print 100,  i,j,a,b,t2_guess(i,j,a,b) , amplitude 
       t2_guess(i,j,a,b) = amplitude 
+!print 100,  i,j,a,b,t2_guess(i,j,b,a) , -amplitude 
+      t2_guess(i,j,b,a) = -amplitude 
+
+      i = i+1
+      j = j-1
+      a = a+1
+      b = b-1
+!print 100,  i,j,a,b,t2_guess(i,j,a,b) , amplitude 
+      t2_guess(i,j,a,b) = amplitude 
+!print 100,  i,j,a,b,t2_guess(i,j,b,a) , -amplitude 
+      t2_guess(i,j,b,a) = -amplitude 
     enddo
     30 continue
     close(iunit)
+ else if (cc_guess == 2) then
+   call random_number(t2_guess)
+   t2_guess *= 1.d-3
  endif
 
 END_PROVIDER

devel/cc/form_r2.irp.f

@@ -29,6 +29,9 @@ subroutine form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
 
   r2(:,:,:,:) = OOVV(:,:,:,:)
 
+  !$OMP PARALLEL DO DEFAULT(NONE) &
+  !$OMP SHARED(nO,nV,r2,aoooo,bvvvv,gvv,goo,tau,OVOO,OVVV,t1,t2,hovvo,OVVO) &
+  !$OMP PRIVATE(i,j,a,b,k,l,c,d)
   do b=1,nV
     do a=1,nV
       do j=1,nO
@@ -130,5 +133,6 @@ subroutine form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
 
 end subroutine form_r2

devel/cc/spatial_to_spin_ERI.irp.f

@@ -21,24 +21,34 @@ BEGIN_PROVIDER [ double precision, dbERI, (spin_mo_num,spin_mo_num,spin_mo_num,s
       do r=1,spin_mo_num
         do q=1,spin_mo_num
           do p=1,spin_mo_num
-
+            pqrs = 0.d0
-            pqrs = Kronecker_delta(mod(p,2),mod(r,2)) &
+            pqsr = 0.d0
-                 * Kronecker_delta(mod(q,2),mod(s,2)) &
+            if ( ( iand(p,1) == iand(r,1) )  .and. &
-                 * get_two_e_integral(                &
+                 ( iand(q,1) == iand(s,1) ) ) then
+!            pqrs = Kronecker_delta(mod(p,2),mod(r,2)) &
+!                 * Kronecker_delta(mod(q,2),mod(s,2)) &
+!                 * get_two_e_integral(                &
+              pqrs = get_two_e_integral(                &
                       (p+1)/2,                         &
                       (q+1)/2,                         &
                       (r+1)/2,                         &
                       (s+1)/2,                         &
                       mo_two_e_integrals_in_map)
+            endif
 
-            pqsr = Kronecker_delta(mod(p,2),mod(s,2)) &
+            if ( ( iand(p,1) == iand(s,1) )  .and. &
-                 * Kronecker_delta(mod(q,2),mod(r,2)) &
+                 ( iand(q,1) == iand(r,1) ) ) then
-                 * get_two_e_integral(                &
+
+!              pqsr = Kronecker_delta(mod(p,2),mod(s,2)) &
+!                   * Kronecker_delta(mod(q,2),mod(r,2)) &
+!                   * get_two_e_integral(                &
+                pqsr = get_two_e_integral(                &
                       (p+1)/2,                         &
                       (q+1)/2,                         &
                       (s+1)/2,                         &
                       (r+1)/2,                         &
                       mo_two_e_integrals_in_map)
+            endif
 
             dbERI(p,q,r,s) = pqrs - pqsr
           enddo

stable/amplitudes/EZFIO.cfg.save

@@ -0,0 +1,7 @@
+[t1_amplitudes]
+type: double precision
+doc: Amplitudes for the single-excitation operator
+interface: ezfio,provider
+size: (mo_basis.mo_num,mo_basis.mo_num)
+
+

stable/amplitudes/NEED

@@ -0,0 +1 @@
+determinants hartree_fock

stable/amplitudes/README.rst

@@ -0,0 +1,6 @@
+==========
+amplitudes
+==========
+
+Computes the amplitudes from a wave function.
+use the QP_STATE environment variable to select an excited state

stable/amplitudes/extract_amplitudes.irp.f

@@ -22,8 +22,13 @@ subroutine run
   integer                        :: exc(0:2,2,2), degree
   integer                        :: h1,h2,p1,p2,s1,s2, istate
   integer                        :: i,j,k,l,a,b,c,d
+  character*(32)                 :: buffer
 
+  call getenv('QP_STATE',buffer)
   istate = 1
+  read(buffer,*,err=5,end=5) istate
+5 continue
+  call write_int(6,istate,'State for amplitudes')
 
   allocate (                                                         &
       t1_a(elec_alpha_num,mo_num),                                   &
@@ -85,14 +90,23 @@ subroutine run
 
         if ( (s1 == 1).and.(s2 == 1) ) then
           t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
+          t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm
+          t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
+          t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm
         else if ( (s1 == 1).and.(s2 == 2) ) then
-          t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
+          t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0
+          t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
         else if ( (s1 == 2).and.(s2 == 1) ) then
-          t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
+          print *,  irp_here, ': Bug!'
+          stop -1
         else if ( (s1 == 2).and.(s2 == 2) ) then
           t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
+          t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm
+          t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
+          t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm
         else
            print *,  irp_here, ': Bug!'
+           stop -2
         endif
         
     end select
@@ -103,9 +117,16 @@ subroutine run
     do a=elec_alpha_num+1, mo_num
       do j=1,elec_alpha_num
         do i=1,elec_alpha_num
+          if ( (i == j).or.(a == b)) then
+            t2_aa(i,j,a,b) = 0.d0
+            t2_bb(i,j,a,b) = 0.d0
+          else
             t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b)
             t2_bb(i,j,a,b) -= t1_b(i,a)*t1_b(j,b)
-          t2_ab(i,j,a,b) -= t1_a(i,a)*t1_b(j,b)
+            t2_aa(i,j,a,b) += t1_a(i,b)*t1_a(j,a)
+            t2_bb(i,j,a,b) += t1_b(i,b)*t1_b(j,a)
+          endif
+          t2_ab(i,j,a,b) -= 0.5d0* (t1_a(i,a)*t1_b(j,b) + t1_b(i,a)*t1_a(j,b))
         enddo
       enddo
     enddo
@@ -188,13 +209,16 @@ subroutine run
       do j=1,elec_alpha_num
         do i=1,elec_alpha_num
           ! The t1 terms are zero because HF
-          e_cor = e_cor + (                                                 &
+          e_cor = e_cor + 0.25d0*(                                   &
               t2_aa(i,j,a,b) + t2_bb(i,j,a,b) +                      &
               t1_a(i,a) * t1_a(j,b)           +                      &
-              t1_b(i,a) * t1_b(j,b) ) * (                            &
+              t1_b(i,a) * t1_b(j,b)           -                      &
+              t1_a(i,b) * t1_a(j,a)           -                      &
+              t1_b(i,b) * t1_b(j,a)                                  &
+              ) * (                                                  &
               get_two_e_integral(i,j,a,b,mo_integrals_map) -         &
               get_two_e_integral(i,j,b,a,mo_integrals_map) )
-          e_cor = e_cor + (                                                 &
+          e_cor = e_cor + 1.0d0  * (                                 &
               t2_ab(i,j,a,b) +                                       &
               t1_a(i,a) * t1_b(j,b) ) *                              &
               get_two_e_integral(i,j,a,b,mo_integrals_map)