From c4b7fda0514904a81c9e2a17f7159e858864e0cc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 10 Sep 2019 16:35:14 +0200 Subject: [PATCH] Created amplitudes module --- devel/cc/amplitude_guess.irp.f | 32 ++++++------------- devel/cc/form_r2.irp.f | 4 +++ stable/amplitudes/EZFIO.cfg.save | 7 ++++ stable/amplitudes/NEED | 1 + stable/amplitudes/README.rst | 5 +++ .../extract_amplitudes.irp.f | 27 +++++++++------- 6 files changed, 42 insertions(+), 34 deletions(-) create mode 100644 stable/amplitudes/EZFIO.cfg.save create mode 100644 stable/amplitudes/NEED create mode 100644 stable/amplitudes/README.rst rename stable/{utilities => amplitudes}/extract_amplitudes.irp.f (92%) diff --git a/devel/cc/amplitude_guess.irp.f b/devel/cc/amplitude_guess.irp.f index 354efd5..d284134 100644 --- a/devel/cc/amplitude_guess.irp.f +++ b/devel/cc/amplitude_guess.irp.f @@ -19,7 +19,6 @@ BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ] read(iunit,*,err=10) i, a, amplitude i = 2*i-1 a = 2*a-1 - spin_occ_num -print '(I4,I4,F20.10,X,F20.10)', i,a,amplitude,t1_guess(i,a) t1_guess(i,a) = amplitude enddo 10 continue @@ -27,7 +26,6 @@ print '(I4,I4,F20.10,X,F20.10)', i,a,amplitude,t1_guess(i,a) read(iunit,*,end=20) i, a, amplitude i = 2*i a = 2*a - spin_occ_num -print '(I4,I4,F20.10,X,F20.10)', i,a,amplitude,t1_guess(i,a) t1_guess(i,a) = amplitude enddo 20 continue @@ -54,7 +52,6 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir integer :: i, j, a, b double precision :: amplitude - t2_guess(:,:,:,:) = 0.d0 iunit = getunitandopen('t2','r') read(iunit,*) do @@ -63,54 +60,43 @@ 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 -print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) - t2_guess(i,j,a,b) = amplitude +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 -print *, '' do read(iunit,*,err=20) i, j, a, b, amplitude i = 2*i j = 2*j a = 2*a - spin_occ_num b = 2*b - spin_occ_num -print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) - t2_guess(i,j,a,b) = amplitude +print 100, i,j,a,b,t2_guess(i,j,a,b) , amplitude + t2_guess(i,j,a,b) = amplitude enddo 20 continue -print *, '' do read(iunit,*,end=30) i, j, a, b, amplitude i = 2*i-1 j = 2*j a = 2*a-1 - spin_occ_num b = 2*b - spin_occ_num -print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) +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 '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) +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) - print *, 'Non-zero amplitudes:' - do i=1,spin_occ_num - do j=1,spin_occ_num - do a=1,spin_vir_num - do b=1,spin_vir_num - if (dabs(t2_guess(i,j,a,b)) > 1.d-16) then - print *, i,j,a,b,t2_guess(i,j,a,b) - endif - enddo - enddo - enddo - enddo else if (cc_guess == 2) then call random_number(t2_guess) t2_guess *= 1.d-3 diff --git a/devel/cc/form_r2.irp.f b/devel/cc/form_r2.irp.f index f869300..d974132 100644 --- a/devel/cc/form_r2.irp.f +++ b/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 diff --git a/stable/amplitudes/EZFIO.cfg.save b/stable/amplitudes/EZFIO.cfg.save new file mode 100644 index 0000000..66ef730 --- /dev/null +++ b/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) + + diff --git a/stable/amplitudes/NEED b/stable/amplitudes/NEED new file mode 100644 index 0000000..d3d4d2c --- /dev/null +++ b/stable/amplitudes/NEED @@ -0,0 +1 @@ +determinants diff --git a/stable/amplitudes/README.rst b/stable/amplitudes/README.rst new file mode 100644 index 0000000..06d99fd --- /dev/null +++ b/stable/amplitudes/README.rst @@ -0,0 +1,5 @@ +========== +amplitudes +========== + +Computes the amplitudes from a wave function. diff --git a/stable/utilities/extract_amplitudes.irp.f b/stable/amplitudes/extract_amplitudes.irp.f similarity index 92% rename from stable/utilities/extract_amplitudes.irp.f rename to stable/amplitudes/extract_amplitudes.irp.f index e025bb8..2a201ea 100644 --- a/stable/utilities/extract_amplitudes.irp.f +++ b/stable/amplitudes/extract_amplitudes.irp.f @@ -84,10 +84,10 @@ subroutine run endif if ( (s1 == 1).and.(s2 == 1) ) then - t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0 - t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 - t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 - t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 + 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 * 0.5d0 t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 @@ -95,10 +95,10 @@ subroutine run 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 * 0.5d0 - t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 - t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 - t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 + 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 @@ -118,6 +118,8 @@ subroutine run 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_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 @@ -202,13 +204,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 + 0.5d0*( & + 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 + 1.0d0 * ( & + 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)