1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
|
subroutine mpdenom(nconf,trec,nstr,intrec,nmax,faaal,faabe,fiial,
$fiibe,f,nnir,isympv,isympo,istr,hval,kmal,jmat,hvbe,kmbe,nbmax,
$wspcb,wspca,wsmax,nactm,imed,intn,noccup,absind,irec,ioffs,mpd,
$ltpr,tprtol,erec,econf)
************************************************************************
* This subroutine calculates MP denominators *
************************************************************************
#include "MRCCCOMMON"
integer nnir,nmax,nactm,i,j,k,intrec(*),nex,ii,jj,nn,ssym,absind
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ioffs
integer econf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer n1,n2,n3,n4,ire,ile,nbmax,wsmax,intn(*),imed(16,1),l,kk
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer erec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer iactv,iacto,iactva,iactoa,iactvb,iactob,noccup(nbmax)
integer nstr(nnir,0:nactm,0:nmax,4),istr(nnir,0:nactm,0:nmax,4)
integer isympv(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2)
integer isympo(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2)
integer isymv,isymo,isymva,isymvb,isymoa,isymob,irv,iro,dbladd
integer wspcb(0:wsmax,13),wspca(0:wsmax,13),nn1,nn2,nn3,nn4
integer nactvintal,nactointal,nactvintbe,nactointbe,nactvirtal
integer nactoccal,nactvirtbe,nactoccbe,irec(nnir,0:nactm,0:nmax,4)
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe,nref,intadd
real*8 f(*),faaal(*),faabe(*),fiial(*),fiibe(*),tprtol
real*8 hval(nbmax),kmal(nbmax,nbmax),jmat(nbmax,nbmax)
real*8 hvbe(nbmax),kmbe(nbmax,nbmax)
logical mpd,ltpr
C
if(mpd) then
rewind(scrfile1)
call dfillzero(ibuf,ibufln)
write(mpfile,rec=trecmax+1) ibuf
if((calc.eq.0.or.conver.eq.1).and.(.not.ptroute).and.
$rest.ge.0) then
write(iout,*)
$'Calculation of diagonal elements of Hamiltonian...'
do i=1,nbmax
read(scrfile1) (kmal(i,j),j=1,nbmax)
enddo
if(.not.spatial) then
do i=1,nbmax
read(scrfile1) (kmbe(i,j),j=1,nbmax)
enddo
endif
read(scrfile1) jmat
read(scrfile1) hval
if(.not.spatial) read(scrfile1) hvbe
call dfillzero(ibuf,ibufln)
write(mpfile,rec=1) ibuf
else
write(iout,*) 'Calculation of MP denominators...'
C Read diagonal Fock-matrix elements
call rfock(faaal,faabe,fiial,fiibe,f,intn,intrec,imed,wspca,
$wsmax,nstr,istr,nnir,nmax,nactm)
endif
else
if(ltpr) then
write(iout,*)
if(calc.eq.1.and.iroot.eq.1) then
write(iout,*) 'Dominant cluster amplitudes'
else
write(iout,*) 'Dominant CI coefficients'
endif
write(iout,"(' Printing threshold: ',1pe9.2)") tprtol
write(iout,*)
endif
rewind(scfile)
nref=0
endif
C Loop over excitations
ssym=isym
if(calc.eq.1.and.(.not.ltpr.or.iroot.eq.1)) ssym=1
if(ptroute) ssym=ptsym(1)
do nex=0,op
do iactv=max(0,nex-mrop),min(nactv,nex)
do iacto=max(0,nex-mrop),min(nacto,nex)
do i1=0,nex
do iactva=max(0,iactv-nactvb),min(nactva,i1,iactv)
iactvb=iactv-iactva
do iactoa=max(0,iacto-nactob),min(nactoa,i1,iacto)
iactob=iacto-iactoa
if(ltpr.and.iroot.gt.1) then
nn=econf(iactva,iactvb,iactoa,iactob,i1,nex)
else
nn=nconf(iactva,iactvb,iactoa,iactob,i1,nex)
endif
if(nn.gt.0) then
nampvirtal=i1
nampvirtbe=nex-nampvirtal
nampoccal=i1
nampoccbe=nex-nampoccal
if(.not.mpd) then
if(ltpr.and.iroot.gt.1.or.(calc.eq.0.and.nroot.gt.1))
$then
c call getlst(tfile,erec(iactva,iactvb,iactoa,iactob,i1,
c $nex),f(1),nn)
call getlst(cfile,(iroot-1)*trecmax+
$erec(iactva,iactvb,iactoa,iactob,i1,nex),f(1),nn)
else
call getlst(tfile,trec(iactva,iactvb,iactoa,iactob,i1,
$nex),f(1),nn)
endif
endif
ii=1
do ir=1,nir
isymv=csympair(ssym,ir,1)
isymo=csympair(ssym,ir,2)
do irv=1,isympv(0,isymv,iactva,iactvb,nampvirtal,
$nampvirtbe,1)
isymva=isympv(irv,isymv,iactva,iactvb,nampvirtal,nampvirtbe,1)
isymvb=isympv(irv,isymv,iactva,iactvb,nampvirtal,nampvirtbe,2)
n1=nstr(isymva,iactva,nampvirtal,1)
n2=nstr(isymvb,iactvb,nampvirtbe,2)
k=n1*n2
do iro=1,isympo(0,isymo,iactoa,iactob,nampoccal,
$nampoccbe,1)
isymoa=isympo(iro,isymo,iactoa,iactob,nampoccal,nampoccbe,1)
isymob=isympo(iro,isymo,iactoa,iactob,nampoccal,nampoccbe,2)
n3=nstr(isymoa,iactoa,nampoccal,3)
n4=nstr(isymob,iactob,nampoccbe,4)
call strread(nn1,nn2,nn3,nn4,n1,n2,n3,n4,nampvirtal,nampvirtbe,
$nampoccal,nampoccbe,isymva,isymvb,isymoa,isymob,iactva,iactvb,
$iactoa,iactob,nnir,nactm,nmax,irec,jj,ioffs+nn)
if(mpd.or.ltpr) then
kk=intadd(jj)
jj=dbladd(kk+nbmax)
call mpdenom1(nampvirtal,nampvirtbe,nampoccal,nampoccbe,n1,n2,n3,
$n4,faaal,faabe,fiial,fiibe,icore(nn1),icore(nn2),icore(nn3),
$icore(nn4),f(ii),hval,hvbe,kmal,kmbe,jmat,nbmax,f(ii),noccup,
$dcore(jj),absind,icore(kk),ltpr,tprtol)
else
if(nampvirtal.eq.iactva.and.nampvirtbe.eq.iactvb.and.
$ nampoccal.eq.iactoa.and.nampoccbe.eq.iactob) !mrszemet
$call mrrefdet(nampvirtal,nampvirtbe,nampoccal,nampoccbe,n1,n2,n3,
$n4,icore(nn1),icore(nn2),icore(nn3),icore(nn4),f(ii),nbmax,
$absind,noccup,nref)
endif
ii=ii+n3*n4*k
enddo
enddo
enddo
if(mpd) call putlst(mpfile,
$trec(iactva,iactvb,iactoa,iactob,i1,nex),f(1),nn)
endif
enddo
enddo
enddo
enddo
enddo
enddo
if(ss1route.and..not.mpd.and..not.ltpr) then
rewind(ssfile)
write(ssfile,*) nref
endif
call szemet(mpfile,nconf,nmax,nactm,trec,f,0,wspca,wsmax,i,
$nconf,trec)
C
return
end
C
| |