diff --git a/src/casscf/bavard.irp.f b/src/casscf/bavard.irp.f index 402e67ec..a9797712 100644 --- a/src/casscf/bavard.irp.f +++ b/src/casscf/bavard.irp.f @@ -1,6 +1,6 @@ ! -*- F90 -*- BEGIN_PROVIDER [logical, bavard] -! bavard=.true. - bavard=.false. + bavard=.true. +! bavard=.false. END_PROVIDER diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f index daf3f68b..1c6d9e6b 100644 --- a/src/casscf/bielec.irp.f +++ b/src/casscf/bielec.irp.f @@ -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 diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f index cb09be3e..9968530c 100644 --- a/src/casscf/bielec_natorb.irp.f +++ b/src/casscf/bielec_natorb.irp.f @@ -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 diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 4270a9fb..1b77cf43 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -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) diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f index dcbcd413..c84b4862 100644 --- a/src/casscf/natorb.irp.f +++ b/src/casscf/natorb.irp.f @@ -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