mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 06:05:44 -04:00
various little fixes. Basically the a3a was working
now the ij/bk symmetry is accounted for by the subroutine the abba=> abba - sign is fixed in code (needs consideration)
This commit is contained in:
Joop van Lenthe
parent
39f609f868
commit
b6a2c3cdd9
2 changed files with 95 additions and 69 deletions
|
|
@ -5,7 +5,8 @@ c ==========================================
|
|||
c CHECK
|
||||
c c/a symmetry
|
||||
c triangular symmetries (ij/ab in summation)
|
||||
c can be sorted by the factors in the occ's
|
||||
c are sorted in cct_uhf_t2 ......
|
||||
c can be (but are not)sorted by the factors in the occ's
|
||||
c ==========================================
|
||||
c
|
||||
c... SYMMETRY of integrals requires careful CONSIDERATION
|
||||
|
|
@ -24,7 +25,7 @@ c... (refers to matrix after possible transpose)
|
|||
c... ld : leading dimensions
|
||||
c
|
||||
*
|
||||
* $Id: cct_UHF.F,v 1.9 1999-07-28 00:42:13 d3e129 Exp $
|
||||
* $Id: cct_UHF.F,v 1.10 1999-09-06 13:36:37 jhvan Exp $
|
||||
*
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -141,18 +142,19 @@ c k m i n a b j c
|
|||
1 '<ij||an>(iajn)','<ij||an>(iajn)'/
|
||||
c======================================================================
|
||||
data a3a/'A3A'/
|
||||
c... this one already corrected for compensating ij symm (cases 1,2)
|
||||
c... cases contains the spin cases and the factor to multiply by
|
||||
c... occ contains the occupations (1=occ,2=virt) and the divide factor
|
||||
c... and (to be save) the inices to be weighted in the intermediate)
|
||||
c... and (to be save) the indices to be weighted in the intermediate)
|
||||
c... type contains types for intgrals/t's in first and second mult
|
||||
c 1 2 3 4 5 6 7 8
|
||||
c a e c f b k i j, fac
|
||||
data cases_a3a/ 1,1,1,1, 0,0, 1,1, 2,
|
||||
2 2,2,2,2, 0,0, 2,2, 2,
|
||||
data cases_a3a/ 1,1,1,1, 0,0, 1,1, 1,
|
||||
2 2,2,2,2, 0,0, 2,2, 1,
|
||||
3 1,2,1,2, 2,1, 1,2, 2,
|
||||
4 2,1,2,1, 1,2, 1,2, 2,
|
||||
5 1,1,2,2, 0,0, 1,2, 4/
|
||||
5 1,1,2,2, 0,0, 1,2, -4/
|
||||
c.. above since we seem to get abab while specifying abba
|
||||
c 5 1,1,2,2, 0,0, 1,2, 4/
|
||||
data occ_a3a / 2,2,2,2, 2,1, 1,1, 2, 1,3/
|
||||
data type_a3a /'<ba||ek>(bkea)','<ba||ek>(bkea)',
|
||||
1 'T(ij,ce)(ijec)','T(ij,ce)(ijec)'/
|
||||
|
|
@ -272,22 +274,34 @@ c r_a2b = cct_uhf_T2(a2b,occ_a2b,cases_a2b,ncase_a2b,type_a2b,
|
|||
c 1 r_a1b)
|
||||
c r_b2b = cct_uhf_T2(b2b,occ_b2b,cases_b2b,ncase_b2b,type_b2b,
|
||||
c 1 r_b1b)
|
||||
c print *,' a2a ',r_a2a
|
||||
c print *,' a2b ',r_a2b
|
||||
c print *,' b2b ',r_b2b
|
||||
c print *,' total term2 ',r_a2a+r_a2b+r_b2b
|
||||
c
|
||||
print *,' diagram ',A3A,ncase_a3a
|
||||
r_a3a = cct_uhf_T2(a3a,occ_a3a,cases_a3a,ncase_a3a,type_a3a,
|
||||
1 flop)
|
||||
print *,' r_a3a ', r_a3a
|
||||
c r_a3b = cct_uhf_T2(a3b,occ_a3b,cases_a3b,ncase_a3b,type_a3b,
|
||||
c 1 flop)
|
||||
c r_b3b = cct_uhf_T2(b3b,occ_b3b,cases_b3b,ncase_b3b,type_b3b,
|
||||
c 1 flop)
|
||||
print *,' a3a ',r_a3a
|
||||
c print *,' a3b ',r_a3b
|
||||
c print *,' b3b ',r_b3b
|
||||
c print *,' total term3 ',r_a3a+r_a3b+r_b3b
|
||||
stop
|
||||
c
|
||||
c r_a4a = cct_uhf_T2(a4a,occ_a4a,cases_a4a,ncase_a4a,type_a4a,
|
||||
c 1 flop)
|
||||
c r_a4b = cct_uhf_T2(a4b,occ_a4b,cases_a4b,ncase_a4b,type_a4b,
|
||||
c 1 flop)
|
||||
c r_b4b = cct_uhf_T2(b4b,occ_b4b,cases_b4b,ncase_b4b,type_b4b,
|
||||
c 1 flop)
|
||||
r_a4a = cct_uhf_T2(a4a,occ_a4a,cases_a4a,ncase_a4a,type_a4a,
|
||||
1 flop)
|
||||
r_a4b = cct_uhf_T2(a4b,occ_a4b,cases_a4b,ncase_a4b,type_a4b,
|
||||
1 flop)
|
||||
r_b4b = cct_uhf_T2(b4b,occ_b4b,cases_b4b,ncase_b4b,type_b4b,
|
||||
1 flop)
|
||||
print *,' a4a ',r_a4a
|
||||
print *,' a4b ',r_a4b
|
||||
print *,' b4b ',r_b4b
|
||||
print *,' total term4 ',r_a4a+r_a4b+r_b4b
|
||||
c
|
||||
CCT_UHF = .true.
|
||||
c
|
||||
|
|
@ -461,14 +475,15 @@ c character*20 type(2,2)
|
|||
c... cases contains the spin cases and the factor to multiply by
|
||||
c... occ contains the occupations (1=occ,2=virt) and the divide factor
|
||||
c... type contains types for intgrals/t's in first and second mult
|
||||
c (ij)t(ij,ce)*t(ij,af) * (bk) <ba||ek><bc||fk>
|
||||
c 1 2 3 4 5 6 7 8
|
||||
c a e c f b k i j fac
|
||||
c data cases/ 1,1,1,1, 0,0, 1,1, 1,
|
||||
c 1 2,2,2,2, 0,0, 2,2, 1,
|
||||
c 2 1,2,1,2, 2,1, 1,2, 2,
|
||||
c 3 2,1,2,1, 1,2, 1,2, 2,
|
||||
c 4 2,1,1,2, 0,0, 1,2, 4/
|
||||
c data occupation / 2,2,2,2, 1,2, 1,1, 2/
|
||||
c 4 1,1,2,2, 0,0, 1,2, 4/
|
||||
c data occupation / 2,2,2,2, 2,1, 1,1, 2/
|
||||
c data type /'<ba||ek>(bkea)','<ba||ek>(bkea)',
|
||||
c 1 'T(ij,ce)(ijec)','T(ij,ce)(ijec)'/
|
||||
c
|
||||
|
|
@ -503,6 +518,7 @@ c
|
|||
integer case,nbatch,nbcha,nbchc,size,size_1
|
||||
logical symmetry,a2a,a1a
|
||||
double precision factor,cct_sum,sum,cct_sum_a1a
|
||||
double precision deb1,deb2,deb3,deb4
|
||||
c
|
||||
sum = 0.0d0
|
||||
a2a = diagram.eq.'A2A'.or.diagram.eq.'A2B'.or.diagram.eq.'B2B'
|
||||
|
|
@ -543,6 +559,11 @@ c
|
|||
spinj = cases(8,case)
|
||||
c
|
||||
factor = cases(9,case)/(occupation(9)*1.0D0)
|
||||
c
|
||||
c... handle triangular factors
|
||||
c
|
||||
if (i.eq.j.and.spini.eq.spinj) factor = factor*2
|
||||
if (b.eq.k.and.spinb.eq.spink) factor = factor*2
|
||||
c
|
||||
eb = s_occ_virt(1,spine,e)
|
||||
fb = s_occ_virt(1,spinf,f)
|
||||
|
|
@ -644,7 +665,10 @@ c
|
|||
c
|
||||
call cct_scale_ac(g_TESULT,ab,ne,cb,nf)
|
||||
c
|
||||
sum = sum + cct_sum(g_RESULT,g_TESULT,na,nc,ne,nf) * factor
|
||||
c sum = sum + cct_sum(g_RESULT,g_TESULT,na,nc,ne,nf) * factor
|
||||
deb1 = cct_sum(g_RESULT,g_TESULT,na,nc,ne,nf) * factor
|
||||
print *,' case ',case,' intermediate ',deb1
|
||||
sum = sum + deb1
|
||||
c
|
||||
if (a1a.and.ab.eq.cb) then
|
||||
value_X1X = value_X1X +
|
||||
|
|
@ -749,8 +773,8 @@ c... figure occupations (1=occ,2-virt)
|
|||
occb = iocc(bb)
|
||||
occk = iocc(kb)
|
||||
c
|
||||
symmetric = (ab.eq.ae.and.cb.eq.ce.and.typeA.eq.typeB.and.
|
||||
1 spine.eq.spinf.and.spine.eq.spinf)
|
||||
symmetric = (ab.eq.cb.and.ae.eq.ce.and.typeA.eq.typeB.and.
|
||||
1 spine.eq.spinf)
|
||||
size = avail_ga
|
||||
if (.not.symmetric) size = size/2
|
||||
c
|
||||
|
|
@ -782,7 +806,9 @@ c
|
|||
1 call errquit('unable to crate global arrray matA',0)
|
||||
if (.not.symmetric) then
|
||||
if (.not.ga_create(MT_DBL,len,lcf,'matA(c,f)',0,0,g_MATB))
|
||||
1 call errquit('unable to crate global arrray matB',0)
|
||||
1 call errquit('unable to create global arrray matB',0)
|
||||
else
|
||||
g_MATB = g_MATA
|
||||
end if
|
||||
c
|
||||
bke = 0
|
||||
|
|
@ -798,16 +824,12 @@ c
|
|||
c
|
||||
c... note ;, the integrals are weighted for the summation indices
|
||||
c
|
||||
* print *,' mata '
|
||||
* call ga_print(g_MATA)
|
||||
* print *,' matb '
|
||||
* call ga_print(g_MATB)
|
||||
call ga_dgemm('T','N',lea,lcf,lbk,1.0d0,g_MATA,g_MATB,
|
||||
1 1.0d0,g_RESULT)
|
||||
c
|
||||
end do
|
||||
c
|
||||
if (symmetric) then
|
||||
if (.not.symmetric) then
|
||||
if (.not.ga_destroy(g_MATB)) call errquit('cant destroy MATB',0)
|
||||
end if
|
||||
if (.not.ga_destroy(g_MATA)) call errquit('cant destroy MATA',0)
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
logical function jvltest(rtdb)
|
||||
*
|
||||
* $Id: jvltest.F,v 1.8 1999-07-08 19:10:52 d3e129 Exp $
|
||||
* $Id: jvltest.F,v 1.9 1999-09-06 13:36:42 jhvan Exp $
|
||||
* calling routine for CCSD(T) trials
|
||||
*
|
||||
implicit none
|
||||
|
|
@ -76,34 +76,35 @@ c
|
|||
call get_T2(dbl_mb(k_t2),occ,virt)
|
||||
call get_EVALS(dbl_mb(k_orb),nbasis,occ,virt)
|
||||
c
|
||||
call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orb),nbasis,occ,virt)
|
||||
* call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orb),nbasis,occ,virt)
|
||||
c
|
||||
call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orb),nbasis,occ,virt)
|
||||
* call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orb),nbasis,occ,virt)
|
||||
c
|
||||
c... now make modified (homo-lumo) orbital energies and
|
||||
c... call all + cullen
|
||||
c
|
||||
if (.not. ma_push_get(mt_dbl, nbasis,'orbn',l_orbn, k_orbn))
|
||||
$ call errquit('orbn',0)
|
||||
call mod_EVALS(dbl_mb(k_orb),dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
* if (.not. ma_push_get(mt_dbl, nbasis,'orbn',l_orbn, k_orbn))
|
||||
* $ call errquit('orbn',0)
|
||||
* call mod_EVALS(dbl_mb(k_orb),dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
c
|
||||
print *,' ORBITAL ENERGIES TO HOMO-LUMO '
|
||||
delta = dbl_mb(k_orb+occ-1) - dbl_mb(k_orb+occ)
|
||||
print *,' delta e used ',delta
|
||||
delta = delta * 3.0d0
|
||||
* print *,' ORBITAL ENERGIES TO HOMO-LUMO '
|
||||
* delta = dbl_mb(k_orb+occ-1) - dbl_mb(k_orb+occ)
|
||||
* print *,' delta e used ',delta
|
||||
* delta = delta * 3.0d0
|
||||
c
|
||||
call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
* call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
c
|
||||
call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
* call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orbn),nbasis,occ,virt)
|
||||
c
|
||||
c... set orbital energies to 1 to get almlof to do cullen
|
||||
c... delta contains de delta e (*3)
|
||||
c
|
||||
call dfill(nbasis,1.0d0,dbl_mb(k_orbn),1)
|
||||
delta = 1.0d0
|
||||
call calc_almlof(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 delta,dbl_mb(k_orbn),
|
||||
2 nbasis,occ,virt,energy,.true.)
|
||||
|
|
@ -114,37 +115,37 @@ c call mod_EVALS(dbl_mb(k_orb),dbl_mb(k_orbn),nbasis,occ,virt)
|
|||
c call dcopy(nbasis,dbl_mb(k_orbn),1,dbl_mb(k_orb),1)
|
||||
c print *,' cw keep ORBITAL ENERGIES as HOMO-LUMO '
|
||||
c
|
||||
if (.not.rtdb_get_info(rtdb,'cct_almlof',i,n_almlof,date))
|
||||
1 then
|
||||
n_almlof = 1
|
||||
np_almlof(1) = 5
|
||||
else
|
||||
if (.not.rtdb_get(rtdb,'cct_almlof',MT_INT,nn_almlof,np_almlof))
|
||||
1 call errquit(' cct_almlof ',0)
|
||||
end if
|
||||
* if (.not.rtdb_get_info(rtdb,'cct_almlof',i,n_almlof,date))
|
||||
* 1 then
|
||||
* n_almlof = 1
|
||||
* np_almlof(1) = 5
|
||||
* else
|
||||
* if (.not.rtdb_get(rtdb,'cct_almlof',MT_INT,nn_almlof,np_almlof))
|
||||
* 1 call errquit(' cct_almlof ',0)
|
||||
* end if
|
||||
c
|
||||
do i=1,n_almlof
|
||||
call int_almlof(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orb),dbl_mb(k_orbn),
|
||||
2 nbasis,occ,virt,energy,np_almlof(i))
|
||||
print *,np_almlof(i),' point almlof integrated energy ',energy
|
||||
end do
|
||||
* do i=1,n_almlof
|
||||
* call int_almlof(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orb),dbl_mb(k_orbn),
|
||||
* 2 nbasis,occ,virt,energy,np_almlof(i))
|
||||
* print *,np_almlof(i),' point almlof integrated energy ',energy
|
||||
* end do
|
||||
c
|
||||
if (.not.rtdb_get_info(rtdb,'cct_laguer',i,n_laguer,date))
|
||||
1 then
|
||||
n_laguer = 1
|
||||
np_laguer(1) = 5
|
||||
else
|
||||
if (.not.rtdb_get(rtdb,'cct_laguer',MT_INT,nn_laguer,np_laguer))
|
||||
1 call errquit(' cct_laguer ',0)
|
||||
end if
|
||||
* if (.not.rtdb_get_info(rtdb,'cct_laguer',i,n_laguer,date))
|
||||
* 1 then
|
||||
* n_laguer = 1
|
||||
* np_laguer(1) = 5
|
||||
* else
|
||||
* if (.not.rtdb_get(rtdb,'cct_laguer',MT_INT,nn_laguer,np_laguer))
|
||||
* 1 call errquit(' cct_laguer ',0)
|
||||
* end if
|
||||
c
|
||||
do i=1,n_laguer
|
||||
call int_laguer(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
1 dbl_mb(k_orb),dbl_mb(k_orbn),
|
||||
2 nbasis,occ,virt,energy,np_laguer(i),rtdb)
|
||||
print *,np_laguer(i),' point laguer integrated energy ',energy
|
||||
end do
|
||||
* do i=1,n_laguer
|
||||
* call int_laguer(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2),
|
||||
* 1 dbl_mb(k_orb),dbl_mb(k_orbn),
|
||||
* 2 nbasis,occ,virt,energy,np_laguer(i),rtdb)
|
||||
* print *,np_laguer(i),' point laguer integrated energy ',energy
|
||||
* end do
|
||||
c
|
||||
c... free all core beyond rint
|
||||
c
|
||||
|
|
@ -784,6 +785,9 @@ c
|
|||
1 dbl_mb(k_w),dbl_mb(k_w+lenijab))
|
||||
if (.not.ma_pop_stack(l_w)) call errquit('pop b3b',0)
|
||||
c
|
||||
print *,' term_a3a *.5 ',term_a3a*0.5d0
|
||||
print *,' term_a3b *1, ',term_a3b
|
||||
print *,' term_b3b *.5 ',term_b3b*0.5d0
|
||||
term3 = term_a3a + 2.0d0*term_a3b + term_b3b
|
||||
c
|
||||
if (opr) then
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue