LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 19:13:29 +01:00
Code Releases Activity

Properly ordered natural MOs

Browse Source
This commit is contained in:
Anthony Scemama 2019-07-04 00:22:44 +02:00
parent 62ef1526a2
commit 932befb2bb
5 changed files with 41 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/casscf/bavard.irp.f
View File

@ -1,6 +1,6 @@
! -*- F90 -*-
BEGIN_PROVIDER [logical, bavard]
! bavard=.true.
bavard=.false.
bavard=.true.
! bavard=.false.
END_PROVIDER

2
src/casscf/bielec.irp.f
View File

@ -52,7 +52,7 @@ END_PROVIDER
BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_inact_orb+n_act_orb, mo_num)]
BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
BEGIN_DOC
! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
! indices are unshifted orbital numbers

48
src/casscf/bielec_natorb.irp.f
View File

@ -4,13 +4,13 @@
! indices are unshifted orbital numbers
END_DOC
implicit none
integer :: i,j,k,l,t,u,p,q,pp
integer :: i,j,k,l,t,u,p,q
double precision, allocatable :: f(:,:,:), d(:,:,:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,l,p,pp,d,f) &
!$OMP PRIVATE(j,k,l,p,d,f) &
!$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
!$OMP bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI)
@ -36,8 +36,7 @@
do k=1,n_core_inact_act_orb
do j=1,mo_num
do p=1,n_act_orb
pp=n_act_orb-p+1
bielec_PQxx_no(list_act(p),j,k,l)=d(pp,j,k)
bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k)
end do
end do
@ -54,9 +53,8 @@
d, n_act_orb)
do k=1,n_core_inact_act_orb
do p=1,n_act_orb
pp=n_act_orb-p+1
do j=1,mo_num
bielec_PQxx_no(j,list_act(p),k,l)=d(pp,j,k)
bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k)
end do
end do
end do
@ -83,10 +81,9 @@
0.d0, &
d, mo_num*mo_num)
do p=1,n_act_orb
pp=n_act_orb-p+1
do k=1,mo_num
do j=1,mo_num
bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,pp)
bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p)
end do
end do
end do
@ -110,10 +107,9 @@
0.d0, &
d, mo_num*mo_num)
do p=1,n_act_orb
pp=n_act_orb-p+1
do k=1,mo_num
do j=1,mo_num
bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp)
bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
end do
end do
end do
@ -133,11 +129,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
! indices are unshifted orbital numbers
END_DOC
implicit none
integer :: i,j,k,l,t,u,p,q,pp
integer :: i,j,k,l,t,u,p,q
double precision, allocatable :: f(:,:,:), d(:,:,:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,l,p,pp,d,f) &
!$OMP PRIVATE(j,k,l,p,d,f) &
!$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
!$OMP bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI)
@ -163,8 +159,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
do l=1,n_core_inact_act_orb
do k=1,n_core_inact_act_orb
do p=1,n_act_orb
pp=n_act_orb-p+1
bielec_PxxQ_no(list_act(p),k,l,j)=d(pp,k,l)
bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l)
end do
end do
end do
@ -193,8 +188,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
do l=1,n_core_inact_act_orb
do j=1,mo_num
do p=1,n_act_orb
pp=n_act_orb-p+1
bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(pp,j,l)
bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l)
end do
end do
end do
@ -221,10 +215,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
0.d0, &
d, mo_num*n_core_inact_act_orb)
do p=1,n_act_orb
pp=n_act_orb-p+1
do l=1,n_core_inact_act_orb
do j=1,mo_num
bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,pp)
bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,p)
end do
end do
end do
@ -248,10 +241,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
0.d0, &
d, mo_num*n_core_inact_act_orb)
do p=1,n_act_orb
pp=n_act_orb-p+1
do k=1,n_core_inact_act_orb
do j=1,mo_num
bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp)
bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
end do
end do
end do
@ -269,11 +261,11 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
! index p runs over the whole basis, t,u,v only over the active orbitals
END_DOC
implicit none
integer :: i,j,k,l,t,u,p,q,pp
integer :: i,j,k,l,t,u,p,q
double precision, allocatable :: f(:,:,:), d(:,:,:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,l,p,pp,d,f) &
!$OMP PRIVATE(j,k,l,p,d,f) &
!$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
!$OMP bielecCI_no,bielecCI,list_act,natorbsCI)
@ -298,8 +290,7 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
do k=1,n_act_orb
do j=1,n_act_orb
do p=1,n_act_orb
pp=n_act_orb-p+1
bielecCI_no(p,j,k,l)=d(pp,j,k)
bielecCI_no(p,j,k,l)=d(p,j,k)
end do
end do
@ -316,9 +307,8 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
d, n_act_orb)
do k=1,n_act_orb
do p=1,n_act_orb
pp=n_act_orb-p+1
do j=1,n_act_orb
bielecCI_no(j,p,k,l)=d(pp,j,k)
bielecCI_no(j,p,k,l)=d(p,j,k)
end do
end do
end do
@ -337,10 +327,9 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
d, n_act_orb*n_act_orb)
do p=1,n_act_orb
pp=n_act_orb-p+1
do k=1,n_act_orb
do j=1,n_act_orb
bielecCI_no(j,k,p,l)=d(j,k,pp)
bielecCI_no(j,k,p,l)=d(j,k,p)
end do
end do
end do
@ -363,10 +352,9 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
d, n_act_orb*n_act_orb)
do p=1,n_act_orb
pp=n_act_orb-p+1
do k=1,n_act_orb
do j=1,n_act_orb
bielecCI_no(j,k,l,list_act(p))=d(j,k,pp)
bielecCI_no(j,k,l,list_act(p))=d(j,k,p)
end do
end do
end do

1
src/casscf/casscf.irp.f
View File

@ -32,6 +32,7 @@ subroutine run
converged = dabs(energy_improvement) < thresh_scf
pt2_max = dabs(energy_improvement / pt2_relative_error)
call update_integrals
mo_coef = NewOrbs
call save_mos
call map_deinit(mo_integrals_map)

38
src/casscf/natorb.irp.f
View File

@ -11,7 +11,7 @@
end do
do i=1,n_act_orb
occnum(list_act(i))=occ_act(n_act_orb-i+1)
occnum(list_act(i))=occ_act(i)
end do
if (bavard) then
@ -31,8 +31,10 @@ END_PROVIDER
! Natural orbitals of CI
END_DOC
integer :: i, j
double precision :: Vt(n_act_orb,n_act_orb)
call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
! call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
call svd(D0tu, size(D0tu,1), natorbsCI,size(natorbsCI,1), occ_act, Vt, size(Vt,1),n_act_orb,n_act_orb)
if (bavard) then
write(6,*) ' found occupation numbers as '
@ -70,7 +72,7 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
BEGIN_DOC
! 4-index transformation of 2part matrices
END_DOC
integer :: i,j,k,l,p,q,pp
integer :: i,j,k,l,p,q
real*8 :: d(n_act_orb)
! index per index
@ -84,9 +86,8 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -103,9 +104,8 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -122,9 +122,8 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -141,9 +140,8 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -162,7 +160,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
BEGIN_DOC
! Transformed one-e integrals
END_DOC
integer :: i,j, p, pp, q
integer :: i,j, p, q
real*8 :: d(n_act_orb)
one_ints_no(:,:)=mo_one_e_integrals(:,:)
@ -172,9 +170,8 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -188,9 +185,8 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
@ -205,7 +201,7 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
BEGIN_DOC
! FCI natural orbitals
END_DOC
integer :: i,j, p, pp, q
integer :: i,j, p, q
real*8 :: d(n_act_orb)
NatOrbsFCI(:,:)=mo_coef(:,:)
@ -215,14 +211,18 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
d(p)=0.D0
end do
do p=1,n_act_orb
pp=n_act_orb-p+1
do q=1,n_act_orb
d(pp)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p)
d(p)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
NatOrbsFCI(j,list_act(p))=d(p)
end do
end do
! call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, &
! NatOrbsFCI, size(NatOrbsFCI,1), &
! Umat, size(Umat,1), 0.d0, &
! NewOrbs, size(NewOrbs,1))
END_PROVIDER