scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-10-05 07:45:59 +02:00
Code Issues Releases Wiki Activity

Created amplitudes module

Browse Source
This commit is contained in:
Anthony Scemama 2019-09-10 16:35:14 +02:00
parent 6778974ae5
commit c4b7fda051
6 changed files with 42 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
devel/cc/amplitude_guess.irp.f
View File

@ -19,7 +19,6 @@ BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ]
read(iunit,*,err=10) i, a, amplitude read(iunit,*,err=10) i, a, amplitude
i = 2*i-1 i = 2*i-1
a = 2*a-1 - spin_occ_num 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 t1_guess(i,a) = amplitude
enddo enddo
10 continue 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 read(iunit,*,end=20) i, a, amplitude
i = 2*i i = 2*i
a = 2*a - spin_occ_num 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 t1_guess(i,a) = amplitude
enddo enddo
20 continue 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 integer :: i, j, a, b
double precision :: amplitude double precision :: amplitude
t2_guess(:,:,:,:) = 0.d0
iunit = getunitandopen('t2','r') iunit = getunitandopen('t2','r')
read(iunit,*) read(iunit,*)
do do
@ -63,54 +60,43 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir
j = 2*j-1 j = 2*j-1
a = 2*a-1 - spin_occ_num a = 2*a-1 - spin_occ_num
b = 2*b-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) 100 format (4(I3,X), 2(F20.10,X))
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 enddo
10 continue 10 continue
print *, ''
do do
read(iunit,*,err=20) i, j, a, b, amplitude read(iunit,*,err=20) i, j, a, b, amplitude
i = 2*i i = 2*i
j = 2*j j = 2*j
a = 2*a - spin_occ_num a = 2*a - spin_occ_num
b = 2*b - 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 t2_guess(i,j,a,b) = amplitude
enddo enddo
20 continue 20 continue
print *, ''
do do
read(iunit,*,end=30) i, j, a, b, amplitude read(iunit,*,end=30) i, j, a, b, amplitude
i = 2*i-1 i = 2*i-1
j = 2*j j = 2*j
a = 2*a-1 - spin_occ_num a = 2*a-1 - spin_occ_num
b = 2*b - 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 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 t2_guess(i,j,b,a) = -amplitude
i = i+1 i = i+1
j = j-1 j = j-1
a = a+1 a = a+1
b = b-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 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 t2_guess(i,j,b,a) = -amplitude
enddo enddo
30 continue 30 continue
close(iunit) 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 else if (cc_guess == 2) then
call random_number(t2_guess) call random_number(t2_guess)
t2_guess *= 1.d-3 t2_guess *= 1.d-3

4
devel/cc/form_r2.irp.f
View File

@ -29,6 +29,9 @@ subroutine form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
r2(:,:,:,:) = OOVV(:,:,:,:) 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 b=1,nV
do a=1,nV do a=1,nV
do j=1,nO 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 end do
end do end do
!$OMP END PARALLEL DO
end subroutine form_r2 end subroutine form_r2

7
stable/amplitudes/EZFIO.cfg.save Normal file
View File

@ -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)

1
stable/amplitudes/NEED Normal file
View File

@ -0,0 +1 @@
determinants

5
stable/amplitudes/README.rst Normal file
View File

@ -0,0 +1,5 @@
==========
amplitudes
==========
Computes the amplitudes from a wave function.

27
stable/utilities/extract_amplitudes.irp.f → stable/amplitudes/extract_amplitudes.irp.f
View File

@ -84,10 +84,10 @@ subroutine run
endif endif
if ( (s1 == 1).and.(s2 == 1) ) then 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,p1,p2) += phase * psi_coef(k,istate) * norm
t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm
t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm
else if ( (s1 == 1).and.(s2 == 2) ) then else if ( (s1 == 1).and.(s2 == 2) ) then
t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0 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 t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
@ -95,10 +95,10 @@ subroutine run
print *, irp_here, ': Bug!' print *, irp_here, ': Bug!'
stop -1 stop -1
else if ( (s1 == 2).and.(s2 == 2) ) then 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,p1,p2) += phase * psi_coef(k,istate) * norm
t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm
t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm
else else
print *, irp_here, ': Bug!' print *, irp_here, ': Bug!'
stop -2 stop -2
@ -118,6 +118,8 @@ subroutine run
else else
t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b) 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_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 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)) 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
@ -202,13 +204,16 @@ subroutine run
do j=1,elec_alpha_num do j=1,elec_alpha_num
do i=1,elec_alpha_num do i=1,elec_alpha_num
! The t1 terms are zero because HF ! 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) + & t2_aa(i,j,a,b) + t2_bb(i,j,a,b) + &
t1_a(i,a) * t1_a(j,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,a,b,mo_integrals_map) - &
get_two_e_integral(i,j,b,a,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) + & t2_ab(i,j,a,b) + &
t1_a(i,a) * t1_b(j,b) ) * & t1_a(i,a) * t1_b(j,b) ) * &
get_two_e_integral(i,j,a,b,mo_integrals_map) get_two_e_integral(i,j,a,b,mo_integrals_map)