mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-10-06 16:25:58 +02:00
Compare commits
3 Commits
aee13ced42
...
217a828b15
Author | SHA1 | Date | |
---|---|---|---|
217a828b15 | |||
c4b7fda051 | |||
6778974ae5 |
@ -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
|
||||
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
|
||||
!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
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
pqsr = 0.d0
|
||||
if ( ( iand(p,1) == iand(r,1) ) .and. &
|
||||
( 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
|
||||
|
||||
pqrs = Kronecker_delta(mod(p,2),mod(r,2)) &
|
||||
* Kronecker_delta(mod(q,2),mod(s,2)) &
|
||||
* get_two_e_integral( &
|
||||
(p+1)/2, &
|
||||
(q+1)/2, &
|
||||
(r+1)/2, &
|
||||
(s+1)/2, &
|
||||
mo_two_e_integrals_in_map)
|
||||
if ( ( iand(p,1) == iand(s,1) ) .and. &
|
||||
( iand(q,1) == iand(r,1) ) ) then
|
||||
|
||||
pqsr = Kronecker_delta(mod(p,2),mod(s,2)) &
|
||||
* Kronecker_delta(mod(q,2),mod(r,2)) &
|
||||
* get_two_e_integral( &
|
||||
(p+1)/2, &
|
||||
(q+1)/2, &
|
||||
(s+1)/2, &
|
||||
(r+1)/2, &
|
||||
mo_two_e_integrals_in_map)
|
||||
! 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
|
||||
|
7
stable/amplitudes/EZFIO.cfg.save
Normal file
7
stable/amplitudes/EZFIO.cfg.save
Normal 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
1
stable/amplitudes/NEED
Normal file
@ -0,0 +1 @@
|
||||
determinants hartree_fock
|
6
stable/amplitudes/README.rst
Normal file
6
stable/amplitudes/README.rst
Normal file
@ -0,0 +1,6 @@
|
||||
==========
|
||||
amplitudes
|
||||
==========
|
||||
|
||||
Computes the amplitudes from a wave function.
|
||||
use the QP_STATE environment variable to select an excited state
|
@ -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
|
||||
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)
|
||||
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_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
|
||||
@ -182,21 +203,24 @@ subroutine run
|
||||
double precision :: e_cor
|
||||
double precision, external :: get_two_e_integral
|
||||
|
||||
e_cor = 0.d0
|
||||
e_cor = 0.d0
|
||||
do b=elec_alpha_num+1, mo_num
|
||||
do a=elec_alpha_num+1, mo_num
|
||||
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 + ( &
|
||||
t2_ab(i,j,a,b) + &
|
||||
t1_a(i,a) * t1_b(j,b) ) * &
|
||||
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)
|
||||
enddo
|
||||
enddo
|
Loading…
Reference in New Issue
Block a user