lukasplevac/easyconfigs-it4i
1
0
Fork 0
mirror of https://code.it4i.cz/sccs/easyconfigs-it4i.git synced 2025-04-16 11:48:06 +01:00
Code Issues Packages Projects Releases Wiki Activity
easyconfigs-it4i/m/MRCC/mrcc_files/xmrcc.f
Lukas Krupcik eb1ea94130 merge master
2024-07-25 10:27:17 +02:00

5742 lines
213 KiB
Fortran
Executable File
Raw Blame History

************************************************************************
program xmrcc
************************************************************************
* Date last modified: 01/04/2005 *
* Author: Mihaly Kallay *
* Version: 459 *
* Description: This program solves the (MR)CC/CI/(left and right *
* hand) EOMCC equations for any excitation level. *
* Routines called: getmemr8,getmem,gconfig,intinp,initialt,mpdenom,*
* ccit,gcoupl *
* Compile: g77 -O4 -o mrcc mrcc.f lambda.f mem.f xalloc.c combin.f *
* -lblaspii -llapack (-lessl on AIX; -ldxml on DEC) *
* References: J. Chem. Phys. 115, 2945 (2001); 117, 980 (2002); *
* 119, 2991 (2003). *
************************************************************************
#include "MRCCCOMMON"
integer icyc,nactoo,nactvo,convero,l2mapmem,itcase2,ntcase,nn,kk
integer i,j,k,l,m,n,nnb,ii,jj,itol,maxbuf,spsht,intadd,dbladd
integer rmem,j1mem,j2mem,j3mem,j4mem,j5mem,j6mem,ioffs
real*8 mimem,mamem,inmem,remem,almem,almemo,tomem,tomemo,coef
character*4 localcc,dfbasis_cor,ovirt
character*5 scftype
character*8 iface,hamilton
character*13 rohftype
character*16 c16
logical lopt,llow,lll,log1,lfock,lcvs,lf12,dof12
C
call mrccini
icyc=0
#if defined (OMP)
C$OMP PARALLEL
C$OMP& DEFAULT(SHARED)
xyzsiz1=OMP_GET_NUM_THREADS()
xyzsize=xyzsiz1-1
ompmem=0.d0
C$OMP END PARALLEL
#endif
tomemo=0.d0
inquire(file='fort.24',exist=log1)
if(.not.log1) goto 1005
open(scrfile6,status='old',file='fort.24',form='unformatted')
rewind(scrfile6)
read(scrfile6,end=1005) icyc,nactoo,nactvo,convero,almem,almemo,
$remem,lopt,llow,nacto,nactv,nocc,nvirt,conver,pert,tomem,tomemo
close(scrfile6,status='delete')
1005 continue
C
llow=.false.
open(intfile ,status='unknown',access='direct',recl=irecln)
open(mpfile ,status='unknown',access='direct',recl=irecln)
open(diisfile,status='unknown',access='direct',recl=irecln)
open(errfile ,status='unknown',access='direct',recl=irecln)
open(scrfile2,status='unknown',access='direct',recl=irecln)
open(scrfile3,status='unknown',access='direct',recl=irecln)
open(scrfile4,status='unknown',access='direct',recl=irecln)
open(strfile ,status='unknown',access='direct',recl=irecln)
open(inp ,status='unknown',file='fort.55')
open(inpfile ,status='unknown',file='fort.56')
recmax=idnint(2.d0*dble(2**30)-10.d0) !No 2Gbyte limit
c recmax=idnint((2.d0*dble(2**30)-1.d0)/dble(irecln))!2Gb limit
c recmax=idnint(dble(2**30)/dble(irecln))
C Default values of options
mrop=2
nsing=1
ntrip=0
l=0
calc=1
dens=0
conver=0
isym=1
isy1=1
isy2=1
diag=0
iclsh=1
k=1
ihf=1
ndoub=0
nacto=0
nactv=0
nroot=1
mr=.false.
eomcc=.false.
eomroute=.false.
left=.false.
hess=.false.
leftroute=.false.
densroute=.false.
ptroute=.false.
ccroute=.false.
t1route=.false.
t2route=.false.
l1route=.false.
l2route=.false.
l3route=.false.
d1route=.false.
d2route=.false.
dbroute=.false.
ss1route=.false.
ss2route=.false.
ss3route=.false.
pertroute=0
pert=0
nmr=.false.
eomgrad=.false.
zroute=.false.
osc=.false.
dipm=.false.
ptsym(1)=1
sacc=0
col=.false.
itol=8
maxex=0
ptfreq=0.d0
fdpci=.false.
dboc=0
lcc=.false.
ssmrcc=0
lfvo=.false.
qcc=0
rel=0
rest=0
lcvs=.false.
lf12=.false.
C Try to read options
rewind(inpfile)
read(inpfile,*,end=1001) mrop,nsing,ntrip,rest,calc,dens,conver,
$isym,diag,iclsh,k,ihf,ndoub,nacto,nactv,itol,maxex,sacc,ptfreq,
$dboc,memory
1001 continue
inquire(file='MINP',exist=lll)
if(lll) then
call getkey('localcc',7,localcc,4)
call getkey('dfbasis_cor',11,dfbasis_cor,4)
call getkey('iface',5,iface,8)
call getkey('hamilton',8,hamilton,8)
call getkey('ovirt',5,ovirt,4)
call getkey('rohftype',8,rohftype,13)
call getkey('scftype',7,scftype,5)
call getkey('cvs',3,c16,16)
lcvs=trim(c16).ne.'off'
call getkey('calc',4,c16,16)
lf12=dof12(c16)
c lfock=localcc.ne.'off '.or.dfbasis_cor.ne.'none'.or.
c $(iface.eq.'dirac '.and.hamilton.eq.'x2cmmf ')
lfock=ovirt.eq.'off '.and.dens.eq.0.and.
$ .not.(scftype.eq.'rohf '.and.rohftype.eq.'standard ')
else
lfock=.false.
endif
icyc=icyc+1
tol=10.d0**dble(-itol)
if(calc.ge.2.and.calc.le.8) then
if(mrop.le.2) calc=1
endif
C MR options
mrop2=max(mrop,2)
if(mod(nactv,2).eq.0) then
nactva=nactv/2
nactvb=nactv/2
else
nactva=(nactv-1)/2
nactvb=(nactv-1)/2+1
endif
if(mod(nacto,2).eq.0) then
nactoa=nacto/2
nactob=nacto/2
else
nactoa=(nacto-1)/2+1
nactob=(nacto-1)/2
endif
C Multiplicities
nroot=nsing+ntrip+ndoub
call getmem(nroot+2,icore(1),multmem)
do m=1,nsing
icore(multmem+m-1)=1
enddo
do m=nsing+1,nsing+ntrip
icore(multmem+m-1)=3
enddo
do m=nsing+ntrip+1,nroot
icore(multmem+m-1)=2
enddo
multip=icore(multmem)
if(iclsh.eq.1) then
cs=.true.
else
cs=.false.
endif
if(k.eq.1) then
spatial=.true.
else
spatial=.false.
endif
calc=iabs(calc)
if(calc.eq.20) then
calc=1
qcc=1
eomcc=.true.
left=.true.
endif
if(calc.eq.9) then
calc=3
lfvo=.true.
endif
if(dboc.ne.0) dens=4
if(calc.gt.0.and.nroot.gt.1) then
eomcc=.true.
if(dens.ne.0) eomgrad=.true.
endif
if(calc.ge.2.and.calc.le.8) then
pert=calc-1
calc=1
endif
if(calc.eq.10) then
calc=0
lcc=.true.
endif
if(calc.gt.10) then
ssmrcc=calc-10
calc=1
open(ssfile,status='old',file='fort.31',form='formatted')
rewind(ssfile)
read(ssfile,*) ss1route,ss2route,ss3route,i,coef
if(ss1route) calc=0
endif
if(dens.lt.0) col=.true.
dens=iabs(dens)
if((dens.gt.0.and.calc.eq.1).or.pert.eq.3) then
eomcc=.true.
left=.true.
endif
if(dens.gt.2.and.dens.le.4) then
hess=.true.
dens=dens-2
endif
if(dens.eq.5.or.dens.eq.6) then
osc=.true.
c dens=1
dens=dens-4
endif
if(dens.eq.7.or.dens.eq.8) then
hess=.true.
third=.true.
dens=dens-6
endif
lc3=left.or.dens.gt.0.or.pert.ne.0
rewind(inp)
read(inp,*,end=1002) m
if(m.le.0) then
if(.not.hess) then
write(iout,*) 'Corrupted input file!'
call mrccend(1)
endif
ptroute=.true.
else
rewind(inp)
read(inp,*,end=1002) nbasis,nocc
read(inp,*,end=1002) (m,i=1,nbasis)
read(inp,*,end=1002) m
rewind(inp)
if(m.eq.-2.or.m.eq.-3) rel=iabs(m)-1
endif
if(rel.gt.0) then
nactva=nactv
nactvb=0
nactoa=nacto
nactob=0
endif
fdpci=calc.eq.0.and.ptroute.and.ptfreq.ne.0.d0
c if(dboc.ne.0.and.calc.eq.0.and.ptroute) then
c dboc=1
c else
c dboc=0
c endif
read(inp,*,end=1002) nbasis,nocc
nbasis=2*nbasis
nvirt=nbasis-nocc
if(mod(nocc,2).ne.0) then
nal=(nocc+1)/2
else
nal=nocc/2
endif
nbe=nocc-nal
C Occupancies and active orbitals
call getmem(nbasis+2,icore(1),noccmem)
call getmem(nbasis+2,icore(1),nactmem)
call getmem(nbasis+2,icore(1),iordmem)
call ifillzero(icore(noccmem),nbasis+1)
call ifillzero(icore(nactmem),nbasis+1)
do m=1,min(nal,nbe)
icore(noccmem+m)=2
enddo
do m=nbe+1,nal
icore(noccmem+m)=1
enddo
do m=nal-max(nactoa,nactob)+1,nal+max(nactva,nactvb)
icore(nactmem+m)=1
enddo
read(inpfile,*,end=1013)
read(inpfile,*,end=1013) (icore(noccmem+m),m=1,nbasis/2)
if(rel.gt.0) then
do m=1,nbasis/2
if(icore(noccmem+m).ne.1.and.icore(noccmem+m).ne.0) then
write(iout,*) 'Invalid occupation numbers! '
call mrccend(1)
endif
enddo
endif
nal=0
nbe=0
do m=1,nbasis/2
if(icore(noccmem+m).eq.2) then
nal=nal+1
nbe=nbe+1
endif
if(iabs(icore(noccmem+m)).eq.1) then
if(icore(noccmem+m).eq.1) nal=nal+1
if(icore(noccmem+m).eq.-1) nbe=nbe+1
if(cs) then
write(iout,*) 'Invalid option! Open shell system!'
call mrccend(1)
endif
endif
enddo
nocc=nal+nbe
read(inpfile,*,end=1013) (icore(nactmem+m),m=1,nbasis/2)
nactva=0
nactvb=0
nactoa=0
nactob=0
do m=1,nbasis/2
if(icore(nactmem+m).eq.1) then
if(icore(noccmem+m).eq.2) then
nactoa=nactoa+1
nactob=nactob+1
endif
if(icore(noccmem+m).eq.1) then
nactoa=nactoa+1
nactvb=nactvb+1
endif
if(icore(noccmem+m).eq.-1) then
nactob=nactob+1
nactva=nactva+1
endif
if(icore(noccmem+m).eq.0) then
nactva=nactva+1
nactvb=nactvb+1
endif
endif
enddo
if(rel.ne.0) then
nactvb=0
nactob=0
endif
nactv=nactva+nactvb
nacto=nactoa+nactob
1013 continue
ms2=nal-nbe
C
if(nactv+nacto.gt.0) mr=.true.
c if(nactv*nacto.gt.0) mr=.true.
if((ndoub.gt.0.and.cs).or.(cs.and.ms2.ne.0).or.
$(cs.and.nsing.eq.0.and.calc.eq.1)) then
write(iout,*) 'Invalid multiplicity!'
call mrccend(1)
endif
nvirtal=nbasis/2-nal
nvirtbe=nbasis/2-nbe
if(rel.gt.0) then
nvirtbe=0
nactvb=0
nactob=0
endif
nvirt=nvirtal+nvirtbe
nelec=nal+nbe
op=min(mrop+min(nactv,nacto),nocc,nvirt)
if(maxex.gt.0) op=min(op,maxex)
if(qcc.eq.1) then
op=min(2*op,nocc,nvirt)
mrop=op
mrop2=max(mrop,2)
endif
nactmax=min(max(nactva,nactvb,nactoa,nactob),op)
op2=op+2
if(op.gt.max(nocc,nvirt)) then
write(iout,*) 'Invalid excitation level!'
call mrccend(1)
endif
rhf=cs.and.spatial
minvirtal=min(op,nvirtal,nal)
minoccal=min(op,nvirtal,nal)
minvirtbe=min(op,nvirtbe,nbe)
minoccbe=min(op,nvirtbe,nbe)
minvirtal2=max(minvirtal,2)
minvirtbe2=max(minvirtbe,2)
minoccal2=max(minoccal,2)
minoccbe2=max(minoccbe,2)
if(rel.gt.0) then
minvirtbe=0
minoccbe=0
endif
oo=max(op,minvirtal2,minvirtbe2,minoccal2,minoccbe2)
lfvo=lfvo.or.(pert.eq.3.and.op.ge.4).or.(ihf.eq.0.and..not.
$(mod(op,2).eq.0.and.(pert.eq.1.or.pert.eq.4)))
if(icyc.eq.1) then
lopt=nacto.eq.0.and.nactv.eq.0
nactoo=nacto
nactvo=nactv
convero=conver
endif
write(iout,1003)
1003 format(1x,70('*'))
if(calc.eq.0) then
if(mr) then
write(iout,*) 'MRCI(',mrop,') calculation '
else
write(iout,*) 'CI(',op,') calculation '
endif
else
if(lcc) then
conver=0
if(mr) then
write(iout,*) 'MRLCC(',mrop,') calculation '
else
write(iout,*) 'LCC(',op,') calculation '
endif
else
if(mr) then
write(iout,*) 'MRCC(',mrop,') calculation '
else
if(pert.eq.1) then
write(iout,"(' CC(',i1,')[',i1,'] calculation ')") op-1,op
else if(pert.eq.2) then
write(iout,"(' CC(',i1,')(',i1,') calculation ')") op-1,op
else if(pert.eq.3) then
write(iout,"(' CC(',i1,')(',i1,')_L calculation ')")
$op-1,op
else if(pert.eq.4) then
write(iout,"(' CC(',i1,')-1a calculation ')") op
else if(pert.eq.5) then
write(iout,"(' CC(',i1,')-1b calculation ')") op
else if(pert.eq.6) then
write(iout,"(' CC',i1,' calculation ')") op
else if(pert.eq.7) then
write(iout,"(' CC(',i1,')-3 calculation ')") op
else
write(iout,*) 'CC(',op,') calculation '
endif
endif
endif
endif
write(iout,*)
call getmem(op+1,icore(1),ntcase)
rewind(gfile)
read(gfile,*,end=1000) i,j,k,maxmem,l
read(gfile,*,end=1000) nimed,wspcmax,iimed,(icore(ntcase+m),
$m=0,op),ccimed,eomimed,leftimed,densimed,l3imed,d2imed,nimed2,
$ds2imed,dt2imed
if(ptroute.and.calc.eq.0) nimed2=d2imed+ccimed
nimed3=nimed2+3*iimed
if(i.ne.op.or.j.ne.nvirt.or.k.ne.nocc.or.l.ne.ihf) goto 1008
maxmem=min(maxmem,dble(memory))
call getmem((oo+1)*(2*nactmax+1)**2,icore(1),itcase)
call getmem(2*(oo+1)*(2*nactmax+1)**2,icore(1),itcase2)
call ifillzero(icore(itcase),(oo+1)*(2*nactmax+1)**2)
call ifillzero(icore(itcase2),2*(oo+1)*(2*nactmax+1)**2)
do m=0,op
do n=1,icore(ntcase+m)
read(gfile,*) i1,i2,ii
icore(itcase+m*(2*nactmax+1)**2+i2*(2*nactmax+1)+i1)=ii
icore(itcase2+m*2*(2*nactmax+1)**2+(ii-1)*2 )=i1
icore(itcase2+m*2*(2*nactmax+1)**2+(ii-1)*2+1)=i2
enddo
enddo
write(iout,*)
write(iout,132) maxmem
132 format(' Allocation of',f6.1, ' Mbytes of memory...')
maxmem=maxmem*dble(twoto20)/dble(ifltln)
maxcor=idnint(maxmem)
maxcor=maxcor+ibufln
call getmemr8(maxcor,icore(1),imem)
c write(6,*) 'getmem imem',imem
eex=0.d0
maxmem=0.d0
nopmax=0.d0
write(iout,*) 'Number of spinorbitals:',nbasis
write(iout,*) 'Number of alpha electrons: ',nal
write(iout,*) 'Number of beta electrons: ',nbe
if(mr) then
write(iout,*) 'Number of active particles:',nactv
write(iout,*) 'Number of active holes: ',nacto
endif
if(nroot.eq.1) write(iout,*) 'Spin multiplicity:',multip
write(iout,1014) ms2*0.5d0
1014 format(' z-component of spin: ',f4.1)
if(isym.eq.0) then
write(iout,*) 'Spatial symmetry is not used.'
else
if(isym.gt.10) then
isy1=(isym-mod(isym,10))/10
isy2=isym-isy1*10
write(iout,*) 'Spatial symmetry: ',isy1,isy2
isym=isy2
else
isy1=isym
isy2=isym
if(eomcc) isy1=1
write(iout,*) 'Spatial symmetry: ',isym
endif
endif
write(iout,2015) tol
2015 format(' Convergence criterion: ',1pe8.1)
C Initialize symmetry variables
call getmem(5*nbasis+1,icore(1),mosymmem)
call symini(icore(mosymmem),icore(mosymmem+nbasis),nbasis,
$icore(noccmem+1),icore(nactmem+1))
if(isym.gt.nir.or.(calc.eq.1.and.isym.ne.1.and.(.not.eomcc)).or.
$isy1.gt.nir.or.isy2.gt.nir) then
write(iout,*) 'Invalid spatial symmetry!'
call mrccend(1)
endif
read(inp,*)
C Number of configurations and occupation graphs
call getmem(4*nir*(oo+1)*(nactmax+1)+(oo+1)**3*(nactmax+1)**4,
$icore(1),iconfmem)
call getmem((oo+1)**2*(nactmax+1)**4,icore(1),itrecmem)
call getmem((oo+1)**2*(nactmax+1)**4,icore(1),itarecmem)
eex=eex-3.d0*dble((oo+1)**2*(nactmax+1)**4-(oo+1)**2*(nactmax+1)
$**4)/dble(iintfp)
call getmem((oo+1)*(2*nactmax+1)**2,icore(1),itrec2mem)
if((eomcc.and.isym.ne.1).or.(ptroute.and.ptsym(1).ne.isym).or.
$(calc.eq.0.and.isy1.ne.isy2)) then
call getmem((oo+1)**2*(nactmax+1)**4,icore(1),econfmem)
call getmem((oo+1)**2*(nactmax+1)**4,icore(1),iearecmem)
call getmem((oo+1)**2*(nactmax+1)**4,icore(1),ierecmem)
eex=eex-3.d0*dble((oo+1)**2*(nactmax+1)**4-(oo+1)**2*(nactmax+1)
$**4)/dble(iintfp)
call getmem((oo+1)*(2*nactmax+1)**2,icore(1),ierec2mem)
else
econfmem=iconfmem+4*nir*(oo+1)*(nactmax+1)
iearecmem=itarecmem
ierecmem=itrecmem
ierec2mem=itrec2mem
endif
call getmem(4*nir*(oo+1)*(nactmax+1),icore(1),istrmem)
call getmem(4*nir*(oo+1)*(nactmax+1),icore(1),istrecmem)
call getmem(2*4*(oo+1)*(nactmax+1),icore(1),itransmem)
call getmem(4*(oo+1)*(nactmax+1),icore(1),iarcmem)
call getmem(2*(oo+2)**3+10,icore(1),itspmem)
call getmem(oo+2,icore(1),iucvmem)
call getmem(oo+5,icore(1),iindmem)
call getmem(2*(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2,icore(1),
$isympmem)
call getmem(nir**4,icore(1),itamem)
call ifillzero(icore(itamem),nir**4)
open(tafile,status='unknown',access='direct',recl=nir**4*iintln)
iimem=intadd(imem)
call gconfig(icore(iconfmem),
$icore(iconfmem+nir*(oo+1)*(nactmax+1)),
$icore(iconfmem+2*nir*(oo+1)*(nactmax+1)),
$icore(iconfmem+3*nir*(oo+1)*(nactmax+1)),
$icore(iconfmem+4*nir*(oo+1)*(nactmax+1)),
$icore(iarcmem),icore(itspmem),oo,nir,
$icore(istrmem),icore(istrmem+nir*(oo+1)*(nactmax+1)),
$icore(istrmem+2*nir*(oo+1)*(nactmax+1)),
$icore(istrmem+3*nir*(oo+1)*(nactmax+1)),
$icore(mosymmem+nbasis),nbasis,icore(isympmem),
$icore(isympmem+(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2),
$icore(itamem),icore(itarecmem),nactmax,icore(econfmem),
$icore(iearecmem),icore(itrecmem),icore(ierecmem),icore(itrec2mem),
$icore(ierec2mem),icore(itcase),icore(ntcase),
$icore(istrecmem),icore(istrecmem+nir*(oo+1)*(nactmax+1)),
$icore(istrecmem+2*nir*(oo+1)*(nactmax+1)),
$icore(istrecmem+3*nir*(oo+1)*(nactmax+1)),
$icore(itransmem),icore(itransmem+2*(oo+1)*(nactmax+1)),
$icore(itransmem+4*(oo+1)*(nactmax+1)),
$icore(itransmem+6*(oo+1)*(nactmax+1)))
C Coupling coefficients
call getmem(4*(oo+1)**2*nir**2*(nactmax+1)**2,icore(1),icoupmem)
call gcoupl(icore(icoupmem),icore(iconfmem),
$icore(iconfmem+nir*(oo+1)*(nactmax+1)),
$icore(iconfmem+2*nir*(oo+1)*(nactmax+1)),
$icore(iconfmem+3*nir*(oo+1)*(nactmax+1)),icore(iarcmem),oo,nir,
$icore(istrmem),icore(istrmem+nir*(oo+1)*(nactmax+1)),
$icore(istrmem+2*nir*(oo+1)*(nactmax+1)),
$icore(istrmem+3*nir*(oo+1)*(nactmax+1)),nactmax,icore(iimem),
$icore(istrecmem),icore(istrecmem+nir*(oo+1)*(nactmax+1)),
$icore(istrecmem+2*nir*(oo+1)*(nactmax+1)),
$icore(istrecmem+3*nir*(oo+1)*(nactmax+1)),
$icore(itransmem),icore(itransmem+2*(oo+1)*(nactmax+1)),
$icore(itransmem+4*(oo+1)*(nactmax+1)),
$icore(itransmem+6*(oo+1)*(nactmax+1)))
C Initialize cluster amplitudes
if(ss3route) then
iimem=intadd(imem)
rmem=iimem
iimem=iimem+nbasis/2
read(ssfile,*) (icore(rmem-1+i),i=1,nbasis/2)
j1mem=iimem
iimem=iimem+(nactmax+1)**4*(oo+1)**2
j2mem=iimem
j3mem=iimem
iimem=iimem+2*4*(oo+1)*(nactmax+1)
j4mem=iimem
j5mem=iimem
j6mem=iimem
c
ioffs=dbladd(iimem)
call cconv(icore(econfmem),icore(ierecmem),icore(iconfmem),oo,
$dcore(ioffs),nir,icore(isympmem),
$icore(isympmem+(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2),nnb,
$nactmax,icore(noccmem+1),icore(iabsmem),icore(istrecmem),ioffs,
$icore(rmem),icore(nactmem+1),icore(j1mem),icore(j2mem),
$icore(mosymmem),icore(iarcmem+4*(oo+1)*(nactmax+1)),icore(j3mem),
$icore(itarecmem),icore(itamem),icore(j4mem),icore(j5mem),
$icore(j6mem))
endif
C Integral lists and intermediates
nnb=nbasis/2
wspcmax=wspcmax*12
call getmem(nimed3+op+wspcmax+19,icore(1),intrecmem)
call getmem(iimed+1,icore(1),denrecmem)
call ifillzero(icore(intrecmem),op2)
intrecmem=intrecmem+op+1
call getmem(leftimed-eomimed+2,icore(1),l2mapmem)
call getmem(nimed3+op+18,icore(1),itypamem)
call getmem(nimed3+op+18,icore(1),itypbmem)
call getmem(nimed3+op+wspcmax+18,icore(1),isymamem)
call getmem(nimed3+op+wspcmax+18,icore(1),isymbmem)
itypamem=itypamem+op+1
itypbmem=itypbmem+op+1
isymamem=isymamem+op+1
isymbmem=isymbmem+op+1
do i=-op,-1
icore(itypamem+i-1)=2
icore(itypbmem+i-1)=2
enddo
call getmem(16*(nimed3+14)+10,icore(1),imedmem)
call getmem((wspcmax+1)*nimed3,icore(1),wspcbmem)
call getmem((wspcmax+1)*(nimed3+op+1),icore(1),wspcamem)
wspcamem=wspcamem+(op+1)*(wspcmax+1)
i=4
if(calc.ne.0.and.(left.or.dens.ne.0)) i=5
call getmem(max(op+2,wspcmax+1,12*nir**i),icore(1),isamem) !Not lower than 6*nir**4 (intinp)
call getmem(14*nir**4,icore(1),iswmem)
call getmem(iimed+1,icore(1),int2mem)
call getmem(8*nnb+1,icore(1),iabsmem)
call readsp(icore(wspcbmem),icore(wspcamem),oo,wspcmax,
$icore(int2mem),iimed,icore(imedmem),icore(intrecmem),
$icore(itypamem),icore(itypbmem),icore(l2mapmem))
jj=1
if(.not.spatial.or.rel.eq.2) jj=2
ii=maxcor-jj*nnb**2
if(lfock) ii=maxcor-2*jj*nnb**2
maxbuf=min(nnb,(ii-mod(ii,nnb**3))/nnb**3)
if(maxbuf.eq.0) then
write(iout,*) 'Insufficient memory!'
call mrccend(1)
endif
nopmax=max(nopmax,dble(nnb**4+jj*nnb**2))
maxmem=max(maxmem,dble(nnb**3+jj*nnb**2))
C Fock-matrix or one-electron integrals?
if(lfock) then
nopmax=max(nopmax,dble(nnb**4+2*jj*nnb**2))
maxmem=max(maxmem,dble(nnb**3+2*jj*nnb**2))
endif
if(rel.eq.2) then
nopmax=max(nopmax,dble(nnb**4+jj*nnb**2)+dble(nnb+1+2*nnb**2)/
$iintfp)
maxmem=max(maxmem,dble(nnb**3+jj*nnb**2)+dble(nnb+1+2*nnb**2)/
$iintfp)
endif
call intinp(dcore(imem),dcore(imem),
$dcore(imem),dcore(imem),maxbuf,nnb,oo,
$icore(iconfmem),icore(intrecmem),icore(imedmem),icore(mosymmem),
$nir,icore(isympmem),
$icore(isympmem+(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2),
$icore(istrmem),icore(isamem),nactmax,icore(wspcbmem),
$icore(wspcamem),wspcmax,icore(int2mem),icore(noccmem+1),
$icore(denrecmem),icore(iabsmem),icore(itypamem),icore(itypbmem),
$icore(isymamem),icore(isymbmem),icore(nactmem+1),icore(iordmem),
$icore(iordmem+nnb))
C Read CVS orbitals
if(lcvs) maxmem=max(maxmem,dble(nocc))
C Calculate many-body denominators
spsht=0
if(.not.spatial) spsht=nnb+2*nnb**2
if(eomcc.and.(.not.ptroute)) calc=0 !szemet, diag H elements
ii=econfmem
jj=ierecmem
kk=iearecmem
347 continue
call mpdenom(icore(ii),icore(jj),icore(iconfmem),icore(intrecmem),
$oo,dcore(imem),dcore(imem),dcore(imem),
$dcore(imem),dcore(imem),nir,
$icore(isympmem),
$icore(isympmem+(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2),
$icore(istrmem),dcore(imem),dcore(imem),
$dcore(imem),dcore(imem),dcore(imem),
$nnb,icore(wspcbmem),icore(wspcamem),wspcmax,nactmax,
$icore(imedmem),icore(int2mem),icore(noccmem+1),icore(iabsmem),
$icore(istrecmem),imem+spsht+nnb+2*nnb**2)
C Initalize CI coefficients
if(eomcc.and.calc.eq.0) then
ii=iconfmem+4*nir*(oo+1)*(nactmax+1)
jj=itrecmem
calc=1
goto 347
endif
if(calc.eq.0.and.isy1.ne.isy2.and.isym.eq.isy2) then
ii=iconfmem+4*nir*(oo+1)*(nactmax+1)
jj=itrecmem
kk=itarecmem
isym=isy1
goto 347
endif
if(sacc.eq.0.and.(pert.eq.0)) close(strfile,status='delete')
close(scrfile1,status='delete')
open(scrfile1,status='unknown',access='direct',recl=irecln)
C CC/CI iteration
c write(6,*) 'ccit imem',imem
call getmem(max(2*nbasis,nir**6+nir**5),icore(1),iwamem)
call getmemr8(3*ndiis**2+ndiis+5,icore(1),idiismem)
call ccit(icore(iconfmem),icore(iconfmem+4*nir*(oo+1)*(nactmax+1))
$,icore(itrecmem),icore(intrecmem),icore(imedmem),dcore(imem),oo,
$icore(istrmem),icore(isympmem),
$icore(isympmem+(oo+1)**2*nir*(nir+1)*2*(nactmax+1)**2),
$icore(itamem),icore(itarecmem),icore(iwamem),icore(isamem),
$icore(iswmem),icore(wspcbmem),icore(wspcamem),nactmax,
$icore(int2mem),icore(ntcase),icore(multmem),icore(iearecmem),
$icore(econfmem),icore(ierecmem),wspcmax,icore(denrecmem),
$icore(itrec2mem),icore(ierec2mem),icore(itcase2),icore(itypamem),
$icore(itypbmem),icore(isymamem),icore(isymbmem),icore(l2mapmem),
$imem)
C F12 szemet
c if(lf12) then
c maxmem=maxmem+(nvirtal*nvirtal*nvirtal*nal*nal*nal)
c $ +(nvirtbe*nvirtbe*nvirtbe*nbe*nbe*nbe)
c $ +(nvirtal*nvirtal*nvirtbe*nal*nal*nbe)
c $ +(nvirtal*nvirtbe*nvirtbe*nal*nbe*nbe)
c endif
C
maxmem=max(maxmem,dble(ndiis+3*ndiis**2))
nopmax=nopmax+dble(ibufln)
nopmax=max(nopmax,maxmem)
mimem=8.d0*maxmem/dble(twoto20)
mamem=8.d0*nopmax/dble(twoto20)
inmem=8.d0*eex/dble(twoto20)
#if defined (OMP)
c inmem=inmem+max(0.d0,8.d0*(ompmem-maxmem)/dble(twoto20))
mimem=mimem+max(0.d0,8.d0*(ompmem-maxmem)/dble(twoto20))
mamem=mamem+max(0.d0,8.d0*(ompmem-maxmem)/dble(twoto20))
#endif
remem=dble(memory)-inmem
almem=min(max(mimem,remem),mamem)
tomem=mimem+inmem
if(.not.log1) then
almemo=almem
tomemo=tomem
nactoo=nacto
nactvo=nactv
endif
if(tomem.le.tomemo) then
c if(almem.le.almemo) then
llow=.true.
almemo=almem
tomemo=tomem
nactoo=nacto
nactvo=nactv
convero=conver
endif
write(6,*)
write(6,*) 'Memory requirements /Mbyte/: '
write(6,*) ' Minimal Optimal'
write(6,6473) 'Real*8: ',mimem,mamem
write(6,6473) 'Integer:',inmem
write(6,6473) 'Total: ',mimem+inmem,mamem+inmem
c write(6,"(f20.12)") mimem+inmem,mamem+inmem
write(6,*)
6473 format(1x,a8,1x,2(f14.4,1x))
C
call updateg(icore(wspcbmem),icore(wspcamem),oo,wspcmax,
$icore(int2mem),iimed,icore(imedmem),icore(intrecmem),
$icore(itypamem),icore(itypbmem),icore(l2mapmem),icore(itspmem),
$itcase2,ntcase,almem)
C
close(mpfile,status='delete')
close(diisfile,status='delete')
close(errfile,status='delete')
close(scrfile1,status='delete')
close(scrfile2,status='delete')
close(scrfile3,status='delete')
close(scrfile4,status='delete')
close(tafile,status='delete')
if(.not.hess) then
do i=0,nfile
close(intfile+i,status='delete')
enddo
endif
if(left) close(scrfile6,status='delete')
if(ptroute.or.eomgrad) then
close(t1file,status='delete')
close(l2file,status='delete')
endif
if(dipm) close(dipfile,status='delete')
C
open(scrfile6,status='unknown',file='fort.24',form='unformatted')
rewind(scrfile6)
write(scrfile6) icyc,nactoo,nactvo,convero,almem,almemo,remem,
$lopt,llow,nacto,nactv,nocc,nvirt,conver,pert,tomem,tomemo
close(scrfile6)
C
call mrccend(0)
1000 continue
write(iout,*) 'Execute goldstone before this program!'
1002 continue
write(iout,*) 'Input file does not exist!'
call mrccend(1)
1008 continue
write(iout,*) 'Formula tape is too old!'
write(iout,*) 'Execute goldstone before this program!'
call mrccend(1)
end
C
************************************************************************
subroutine intinp(vint,fa,fb,v,maxinc,nbmax,nmax,nstr,intrec,imed,
$mosym,nnir,isympv,isympo,istr,iwa,nactm,wspcb,wspca,wsmax,intn,
$noccup,denrec,absind,itypa,itypb,imedsyma,imedsymb,nactd,iorda,
$iordb)
************************************************************************
* This subroutine reads and sorts integrals and constructs correlation *
* energy operator (temporary) *
************************************************************************
#include "MRCCCOMMON"
integer nnir,nbmax,nactm,i,j,k,l,ii,absind(2*nbmax,-3:0),ire,nspc
integer nmax,irec,mosym(*),jj,kk,nstr,wsmax,iw,intn(*),intrec(*)
integer isc,imed(16,1),intl,kkk,noccup(*),istr,ispc,irecold,maxinc
integer isympv,isympo,iwa,wslen,ssym,ndist,idist,il1,il11,il2,ix
integer wspcb(0:wsmax,13),wspca(0:wsmax,13),denrec(13),ifile
real*8 vint(nbmax,nbmax,nbmax,maxinc),v(*),itol,sum
real*8 fa(nbmax,nbmax),fb(nbmax,nbmax)
integer ncn,ncnt,imedsyma(*),itypa(*),itypb(*),imedsymb(*),diprec
integer nactd(*),iorda(*),iordb(*)
parameter(itol=10e-20)
C Numbering of orbitals
C Virtual alpha
ii=0
do i=1,nbasis/2
if(noccup(i).le.0.and.nactd(i).eq.1) then
ii=ii+1
absind(ii,-3)=i
endif
enddo
do i=1,nbasis/2
if(noccup(i).le.0.and.nactd(i).eq.0) then
ii=ii+1
absind(ii,-3)=i
endif
enddo
C Virtual beta
ii=0
do i=1,nbasis/2
if((noccup(i).eq.1.or.noccup(i).eq.0).and.nactd(i).eq.1) then
ii=ii+1
absind(ii,-2)=i
endif
enddo
do i=1,nbasis/2
if((noccup(i).eq.1.or.noccup(i).eq.0).and.nactd(i).eq.0) then
ii=ii+1
absind(ii,-2)=i
endif
enddo
C Occupied alpha
ii=0
do i=1,nbasis/2
if(noccup(i).ge.1.and.nactd(i).eq.0) then
ii=ii+1
absind(ii,-1)=i
endif
enddo
do i=1,nbasis/2
if(noccup(i).ge.1.and.nactd(i).eq.1) then
ii=ii+1
absind(ii,-1)=i
endif
enddo
C Occupied beta
ii=0
do i=1,nbasis/2
if((noccup(i).eq.-1.or.noccup(i).eq.2).and.nactd(i).eq.0) then
ii=ii+1
absind(ii,0)=i
endif
enddo
do i=1,nbasis/2
if((noccup(i).eq.-1.or.noccup(i).eq.2).and.nactd(i).eq.1) then
ii=ii+1
absind(ii,0)=i
endif
enddo
C Initialize wsfile
if(nroot.gt.1.and.((calc.eq.0.and.isy1.ne.isy2).or.(left.and.
$(isym.ne.1.or.eomgrad)))) then
wslen=iintln*(2*(wsmax+1)+3)
open(wsfile,status='unknown',access='direct',recl=wslen)
if(calc.eq.0) then
il1=iimed
il2=ccimed
else
il1=eomimed
il2=nimed
endif
imedsyma(il1+1:il2)=0
imedsymb(il1+1:il2)=0
do i=il1+1,il2
call wswrite(intrec(i),imedsyma(i),imedsymb(i),wspca(0,i),
$wspcb(0,i),wsmax+1,i-il1)
enddo
endif
C Length of intermediates and structure of intermediate file
ptfst=0
irecold=0
irec=0
C Integral list
c hsym,rsym,lsym,csym
call imedlen(1,isym,isym,ptsym(1),0,iimed,0,itypa,itypb,wspca,
$wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
ccfst=irec
C CC intermediates
call imedlen(1,isym,isym,ptsym(1),iimed,ccimed,0,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
leftfst=irec
C Right-hand EOMCC intermediates
if(eomcc) then
call imedlen(1,isym,isym,isym,ccimed,eomimed,0,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
C Left-hand EOMCC intermediates
if(left) then
irecold=max(irec,irecold)
ssym=isym
if(.not.eomgrad) irec=leftfst
if(eomgrad.and.osc) ssym=1
call imedlen(1,isym,ssym,isym,eomimed,leftimed,0,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
if(eomgrad.and.osc.and.nroot.gt.2) leftfst=irec
C Density-matrix intermediates
if(dens.gt.0.and.(.not.eomgrad)) then
ssym=isym
kk=isym
if(calc.eq.0) then
ssym=ptsym(1)
if(isy1.ne.isy2) then
kk=isy1
ssym=isy2
endif
endif
call imedlen(1,kk,ssym,kk,leftimed,densimed,0,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
ptfst=irec
if(ptroute) then
C T1 intermediates, <K|exp(-T) dH/dx exp(T)|0>
call imedlen(mult(isym,ptsym(1)),ptsym,isym,isym,d2imed,t1imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
C L2 intermediates, <0|(1+L) exp(-T) dH/dx exp(T)|K>
call imedlen(mult(isym,ptsym(1)),ptsym,isym,isym,d1imed,l2imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
C T2 intermediates, <K|[exp(-T) H exp(T), dT/dx]|0>
irecold=max(irec,irecold)
irec=ptfst
call imedlen(1,ptsym,isym,isym,t1imed,t2imed,d2imed,itypa,itypb
$,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
if(ptroute.or.eomgrad) then
C L3 intermediates, <0|(1+L)[exp(-T) H exp(T), dT/dx]|K>
ssym=isym
kk=isym
if(ptroute) ssym=ptsym(1)
if(eomgrad.and.osc) kk=1
call imedlen(1,ssym,kk,isym,densimed,l3imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
C D2 intermediates, <0|(1+L)[exp(-T) {p+q-} exp(T), dT/dx]|0>
call imedlen(1,ssym,kk,isym,l3imed,d2imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
if(ptroute) then
C L1 intermediates, <0|dL/dx exp(-T) H exp(T)|K>
irecold=max(irec,irecold)
irec=ptfst
call imedlen(1,ptsym,ptsym,isym,t2imed,l1imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
C D1 intermediates, <0|dL/dx exp(-T) {p+q-} exp(T)|0>
call imedlen(1,ptsym,ptsym,isym,l1imed,d1imed,
$d2imed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,
$nactm,irec,intrec,nstr,denrec,ptfst,intn,imedsyma,imedsymb)
endif
C Integral lists from diagrams 12 and 13
if(ptroute) then
ptsh=d2imed
call diag12(0,i,12, 3,intn,intrec,wspca,denrec,wsmax,imed,
$.true.)
call diag12(0,i,13,11,intn,intrec,wspca,denrec,wsmax,imed,
$.true.)
call diag12(d2imed,drecmax,12, 3,intn,intrec,wspca,denrec,wsmax,
$imed,.true.)
call diag12(d2imed,drecmax,13,11,intn,intrec,wspca,denrec,wsmax,
$imed,.true.)
else
ptsh=0
call diag12(0,drecmax,12, 3,intn,intrec,wspca,denrec,wsmax,imed,
$.true.)
call diag12(0,drecmax,13,11,intn,intrec,wspca,denrec,wsmax,imed,
$.true.)
endif
C
call wszero(wsmax,wspca,0,nimed2)
call wszero(wsmax,wspcb,0,nimed2)
C Ground-state (tot. symmetric) lh EOMCC intermediates
if(nroot.gt.1.and.left.and.(isym.ne.1.or.eomgrad)) then
irecold=max(irec,irecold)
irec=leftfst
call wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
ssym=1
if(eomgrad.and.osc) ssym=isym
C Left-hand EOMCC intermediates for tot.sym. L amplitudes
call imedlen(1,1,ssym,isym,eomimed,leftimed,0,itypa,itypb,wspca,
$wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,nstr,
$denrec,ptfst,intn,imedsyma,imedsymb)
C EOMCC density-matrix intermediates
if(eomgrad)
$call imedlen(1,1,ssym,isym,leftimed,densimed,0,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
call wszero(wsmax,wspca,eomimed,nimed2)
call wszero(wsmax,wspcb,eomimed,nimed2)
if(eomgrad.and.osc) call wschange(wspca,wspcb,wsmax,intrec,
$imedsyma,imedsymb)
endif
C Ground-state CI "intermediates"
if(calc.eq.0.and.isy1.ne.isy2) then
irecold=max(irec,irecold)
irec=ccfst
call wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
call imedlen(1,isy1,isy1,isy1,iimed,ccimed,0,itypa,itypb,wspca,
$wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,ptfst,intn,imedsyma,imedsymb)
call wszero(wsmax,wspca,iimed,ccimed)
call wszero(wsmax,wspcb,iimed,ccimed)
endif
C Initialize integral file
irec=max(irec,irecold)
maxrec=irec
sum=dble(irec)*dble(ifltln)*dble(ibufln)/dble(twoto20)
write(iout,113) sum
113 format(' Length of intermediate file (Mbytes): ',f9.1)
j=mod(irec,recmax)
nfile=(irec-j)/recmax
jj=intfile
do k=1,nfile
open(intfile+k,status='unknown',access='direct',recl=irecln)
enddo
if(.not.ptroute) ptfst=0
l=mod(ptfst,recmax)
ifile=(ptfst-l)/recmax
C Reading and sorting integrals
insym=1
if(ptroute) insym=mult(isym,ptsym(1))
C Try to read dipole integrals
if(osc) then
dipm=.true.
if(calc.ne.0) ptsh=nimed2
open(dipfile,status='unknown',access='direct',recl=irecln)
C Length of dipole moment integrals
diprec=0
irecold=0
do i=1,3
call imedlen(denssym,1,1,1,nimed+(i-1)*iimed,nimed+i*iimed,
$nimed+(i-1)*iimed,itypa,itypb,wspca,wspcb,wsmax,imed,nmax,nnir,
$isympv,isympo,nactm,diprec,intrec,nstr,denrec,irecold,intn,
$imedsyma,imedsymb)
call diag12(nimed+(i-1)*iimed,drecmax,12, 3,intn,intrec,wspca,
$denrec,wsmax,imed,.true.)
call diag12(nimed+(i-1)*iimed,drecmax,13,11,intn,intrec,wspca,
$denrec,wsmax,imed,.true.)
irecold=diprec
enddo
call wszero(wsmax,wspca,nimed2,nimed3)
call wszero(wsmax,wspcb,nimed2,nimed3)
C Initialize dipole integral file
write(iout,"(' Length of dipole integral file (Mbytes): ',f9.1)")
$dble(diprec)*dble(ifltln)*dble(ibufln)/dble(twoto20)
endif
C
return
end
C
************************************************************************
subroutine gconfig(nvstral,nvstrbe,nostral,nostrbe,nconf,iarc,tsp,
$nmax,nnir,ivstral,ivstrbe,iostral,iostrbe,mosym,nb,isympv,isympo,
$ita,tarec,nactm,econf,earec,trec,erec,trec2,erec2,tcase,ntcase,
$ivrecal,ivrecbe,iorecal,iorecbe,itrecva,itrecvb,itrecoa,itrecob)
************************************************************************
* This subroutine generates virtual and occupied excitation operator *
* strings and their lexical indices. *
* For references see e.g. J. Phys. A 18, 3283 (1985); Adv. Quant. Chem.*
* 34, 143 (1999). *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nb,ii,mosym(nb,4),nact,nactm,nacta,nactb,ntcase
integer nvstral(nnir,0:nactm,0:nmax),nvstrbe(nnir,0:nactm,0:nmax)
integer nostral(nnir,0:nactm,0:nmax),nostrbe(nnir,0:nactm,0:nmax)
integer ivstral(nnir,0:nactm,0:nmax),ivstrbe(nnir,0:nactm,0:nmax)
integer iostral(nnir,0:nactm,0:nmax),iostrbe(nnir,0:nactm,0:nmax)
integer iactv,iacto,iactva,iactvb,minva,minvb
integer nconf,iarc(0:nactm,0:nmax,4),iactoa,econf,tarec,earec,ssym
integer tsp(0:2*(nmax+1)**2,0:nmax),nex,i,j,k,m,isum,istr2mem,ir1
integer isymv,isymo,isymva,isymvb,isymoa,isymob,irv,iro,iactob
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 ita(nnir,nnir,nnir,nnir),trec,erec,trec2,erec2,tcase
integer ivrecal(nnir,0:nactm,0:nmax),ivrecbe(nnir,0:nactm,0:nmax)
integer iorecal(nnir,0:nactm,0:nmax),iorecbe(nnir,0:nactm,0:nmax)
integer itrecva(2,0:nactm,0:nmax),itrecvb(2,0:nactm,0:nmax)
integer itrecoa(2,0:nactm,0:nmax),itrecob(2,0:nactm,0:nmax)
C
write(iout,*) 'Construction of occupation graphs...'
do i=0,op
read(gfile,*) tsp(0,i),(tsp(j,i),j=1,2*tsp(0,i))
enddo
i=(oo+1)*nir*(nactm+1)
call ifillzero(nvstral,i)
call ifillzero(nvstrbe,i)
call ifillzero(nostral,i)
call ifillzero(nostrbe,i)
i=(oo+1)**2*nir*(nir+1)*2*(nactm+1)**2
call ifillzero(isympv,i)
call ifillzero(isympo,i)
i=2*(nmax+1)*(nactm+1)
call ifillzero(itrecva,i)
call ifillzero(itrecvb,i)
call ifillzero(itrecoa,i)
call ifillzero(itrecob,i)
istrc=1
minva=max(min(op,nvirtal,nal),2)
minvb=max(min(op,nvirtbe,nbe),2)
C Virtual alpha strings
do nex=0,minva
do nact=max(0,nex-mrop2,nex-nvirtal+nactva),min(nactva,nex)
call getmem((nvirtal+1)*(nex+1),icore(1),iarc(nact,nex,1))
call occgraph(nvirtal,nex,icore(iarc(nact,nex,1)),
$nvstral(1,nact,nex),ivstral(1,nact,nex),
$mosym(1,1),nactva,nact,.true.,ivrecal(1,nact,nex),
$itrecva(1,nact,nex))
if(nex.le.2) then
do ir=1,nir
call getmem(nvstral(ir,nact,nex)*nex,icore(1),istr2mem)
do i=0,nvstral(ir,nact,nex)*nex-1
icore(istr2mem+i)=icore(ivstral(ir,nact,nex)+i)
enddo
ivstral(ir,nact,nex)=istr2mem
enddo
endif
enddo
enddo
C Virtual beta strings
if(cs.and.nvirtal.eq.nvirtbe.and.minvirtal2.eq.minvirtbe2.and.
$nactva.eq.nactvb) then
do nex=0,minvb
do nact=max(0,nex-mrop2),min(nactvb,nex)
iarc(nact,nex,2)=iarc(nact,nex,1)
itrecvb(1,nact,nex)=itrecva(1,nact,nex)
itrecvb(2,nact,nex)=itrecva(2,nact,nex)
do ir=1,nir
nvstrbe(ir,nact,nex)=nvstral(ir,nact,nex)
ivstrbe(ir,nact,nex)=ivstral(ir,nact,nex)
ivrecbe(ir,nact,nex)=ivrecal(ir,nact,nex)
enddo
enddo
enddo
else
do nex=0,minvb
do nact=max(0,nex-mrop2,nex-nvirtbe+nactvb),min(nactvb,nex)
call getmem((nvirtbe+1)*(nex+1),icore(1),iarc(nact,nex,2))
call occgraph(nvirtbe,nex,icore(iarc(nact,nex,2)),
$nvstrbe(1,nact,nex),ivstrbe(1,nact,nex),
$mosym(1,2),nactvb,nact,.true.,ivrecbe(1,nact,nex),
$itrecvb(1,nact,nex))
if(nex.le.2) then
do ir=1,nir
call getmem(nvstrbe(ir,nact,nex)*nex,icore(1),istr2mem)
do i=0,nvstrbe(ir,nact,nex)*nex-1
icore(istr2mem+i)=icore(ivstrbe(ir,nact,nex)+i)
enddo
ivstrbe(ir,nact,nex)=istr2mem
enddo
endif
enddo
enddo
endif
C Symmetry pairs for virtual strings
do i=0,minva
do nacta=max(0,i-mrop2),min(nactva,i)
do j=0,minvb
do nactb=max(0,j-mrop2),min(nactvb,j)
if(i+j-mrop2.le.nacta+nactb) then
do ir=1,nir
isum=0
do irv=1,nir
if(nvstral(isympair(ir,irv,1),nacta,i).gt.0.and.
$ nvstrbe(isympair(ir,irv,2),nactb,j).gt.0) then
isum=isum+1
isympv(isum,ir,nacta,nactb,i,j,1)=isympair(ir,irv,1)
isympv(isum,ir,nacta,nactb,i,j,2)=isympair(ir,irv,2)
endif
enddo
isympv(0,ir,nacta,nactb,i,j,1)=isum
enddo
endif
enddo
enddo
enddo
enddo
C Occupied alpha strings
do nex=0,minoccal2
do nact=max(0,nex-mrop2,nex-nal+nactoa),min(nactoa,nex)
call getmem((nal+1)*(nex+1),icore(1),iarc(nact,nex,3))
call occgraph(nal,nex,icore(iarc(nact,nex,3)),
$nostral(1,nact,nex),iostral(1,nact,nex),
$mosym(1,3),nactoa,nact,.false.,iorecal(1,nact,nex),
$itrecoa(1,nact,nex))
if(nex.le.2) then
do ir=1,nir
call getmem(nostral(ir,nact,nex)*nex,icore(1),istr2mem)
do i=0,nostral(ir,nact,nex)*nex-1
icore(istr2mem+i)=icore(iostral(ir,nact,nex)+i)
enddo
iostral(ir,nact,nex)=istr2mem
enddo
endif
enddo
enddo
C Occupied beta strings
if(cs.and.nal.eq.nbe.and.minoccal2.eq.minoccbe2.and.
$nactoa.eq.nactob) then
do nex=0,minoccbe2
do nact=max(0,nex-mrop2),min(nactob,nex)
iarc(nact,nex,4)=iarc(nact,nex,3)
itrecob(1,nact,nex)=itrecoa(1,nact,nex)
itrecob(2,nact,nex)=itrecoa(2,nact,nex)
do ir=1,nir
nostrbe(ir,nact,nex)=nostral(ir,nact,nex)
iostrbe(ir,nact,nex)=iostral(ir,nact,nex)
iorecbe(ir,nact,nex)=iorecal(ir,nact,nex)
enddo
enddo
enddo
else
do nex=0,minoccbe2
do nact=max(0,nex-mrop2,nex-nbe+nactob),min(nactob,nex)
call getmem((nbe+1)*(nex+1),icore(1),iarc(nact,nex,4))
call occgraph(nbe,nex,icore(iarc(nact,nex,4)),
$nostrbe(1,nact,nex),iostrbe(1,nact,nex),
$mosym(1,4),nactob,nact,.false.,iorecbe(1,nact,nex),
$itrecob(1,nact,nex))
if(nex.le.2) then
do ir=1,nir
call getmem(nostrbe(ir,nact,nex)*nex,icore(1),istr2mem)
do i=0,nostrbe(ir,nact,nex)*nex-1
icore(istr2mem+i)=icore(iostrbe(ir,nact,nex)+i)
enddo
iostrbe(ir,nact,nex)=istr2mem
enddo
endif
enddo
enddo
endif
C Symmetry pairs for occupied strings
do i=0,minoccal2
do nacta=max(0,i-mrop2),min(nactoa,i)
do j=0,minoccbe2
do nactb=max(0,j-mrop2),min(nactob,j)
if(i+j-mrop2.le.nacta+nactb) then
do ir=1,nir
isum=0
do irv=1,nir
if(nostral(isympair(ir,irv,1),nacta,i).gt.0.and.
$ nostrbe(isympair(ir,irv,2),nactb,j).gt.0) then
isum=isum+1
isympo(isum,ir,nacta,nactb,i,j,1)=isympair(ir,irv,1)
isympo(isum,ir,nacta,nactb,i,j,2)=isympair(ir,irv,2)
endif
enddo
isympo(0,ir,nacta,nactb,i,j,1)=isum
enddo
endif
enddo
enddo
enddo
enddo
C Number of configurations
ii=0
trecmax=0
ssym=isym
if(hess) then
if(ptroute) then
ssym=ptsym(1)
write(iout,*) 'Number of perturbed amplitudes: '
endif
do ir1=1,nir
if(ir1.ne.ssym)
$call nconfig(nvstral,nvstrbe,nostral,nostrbe,nconf,tsp,nmax,
$nnir,isympv,isympo,ita,tarec,nactm,i,ir1,trec,trec2,tcase,ntcase,
$.false.)
enddo
endif
call nconfig(nvstral,nvstrbe,nostral,nostrbe,econf,tsp,nmax,nnir,
$isympv,isympo,ita,earec,nactm,ii,ssym,erec,erec2,tcase,ntcase,
$.true.)
if(eomcc.and.isym.ne.1.or.(ptroute.and.ptsym(1).ne.isym.or.
$(calc.eq.0.and.isy1.ne.isy2))) then
write(iout,*) 'Number of ground state excitations: '
ssym=1
if(calc.eq.0) then
ssym=isym
if(isy1.ne.isy2) ssym=isy1
endif
call nconfig(nvstral,nvstrbe,nostral,nostrbe,nconf,tsp,nmax,
$nnir,isympv,isympo,ita,tarec,nactm,ii,ssym,trec,trec2,tcase,ntcase
$,.true.)
endif
C
return
end
C
************************************************************************
subroutine occgraph(norb,nex,idrt,nstr,istr,mosym,nactvo,
$nact,lvirt,irec,itrecv)
************************************************************************
* This subroutine generates occupation graph of excitation operators *
* and generates states in symmetry order *
************************************************************************
#include "MRCCCOMMON"
integer norb,nex,idrt(0:norb,0:nex),nstr(*),istr(*)
integer mosym(*),i,j,ii,m,n,nactvo,nact,irec(*),itrecv(*)
logical lvirt
C Vertex weights
do i=0,norb
do j=0,nex
idrt(i,j)=0
enddo
enddo
if(lvirt) then
m=nactvo-nact
do i=0,nactvo
ii=min(i,nact)
do j=max(0,i-m,ii-nact),ii
if(j.eq.0) then
idrt(i,j)=1
else
idrt(i,j)=idrt(i-1,j-1)+idrt(i-1,j)
endif
enddo
enddo
m=(norb-nactvo)-(nex-nact)
do i=0,(norb-nactvo)
ii=min(i,(nex-nact))
do j=max(0,i-m,ii-(nex-nact)),ii
if(j+nact.eq.0) then
idrt(i+nactvo,j+nact)=1
else
idrt(i+nactvo,j+nact)=idrt(i+nactvo-1,j+nact-1)+
$ idrt(i+nactvo-1,j+nact)
endif
enddo
enddo
else
m=(norb-nactvo)-(nex-nact)
do i=0,norb-nactvo
ii=min(i,nex-nact)
do j=max(0,i-m,ii-(nex-nact)),ii
if(j.eq.0) then
idrt(i,j)=1
else
idrt(i,j)=idrt(i-1,j-1)+idrt(i-1,j)
endif
enddo
enddo
m=nactvo-nact
do i=0,nactvo
ii=min(i,nact)
do j=max(0,i-m,ii-nact),ii
if(j+(nex-nact).eq.0) then
idrt(i+(norb-nactvo),j+(nex-nact))=1
else
idrt(i+(norb-nactvo),j+(nex-nact))=
$idrt(i+(norb-nactvo)-1,j+(nex-nact)-1)+
$idrt(i+(norb-nactvo)-1,j+(nex-nact))
endif
enddo
enddo
endif
n=idrt(norb,nex)
C Generate strings
istr(1)=iimem+n
iivirt=0
iiocc=0
istrv=1
call opstr1(norb,nex,icore(iucvmem),idrt,icore(istr(1)))
call ifillzero(nstr,nir)
call symanal(icore(iimem),icore(istr(1)),nex,n,mosym,nstr,
$icore(iimem+n*(nex+1)))
itrecv(1)=istrc
itrecv(2)=n
call putint(strfile,itrecv(1),icore(iimem),n)
call irecupd(istrc,n)
do i=1,nir
if(i.gt.1) istr(i)=istr(i-1)+nex*nstr(i-1)
m=nex*nstr(i)
irec(i)=istrc
call putint(strfile,irec(i),icore(istr(i)),m)
call irecupd(istrc,m)
enddo
C
return
end
C
************************************************************************
subroutine irecupd(ircm,n)
************************************************************************
* Calculate records for an integer direct access file *
************************************************************************
#include "MRCCCOMMON"
integer irc,ircm,n,i,k,m
C
i=mod(n*iintln,ifltln)
if(i.ne.0) i=ifltln-i
m=(n*iintln+i)/ifltln
k=mod(m,ibufln)
if(k.ne.0) k=ibufln-k
ircm=ircm+(m+k)/ibufln
C
return
end
C
************************************************************************
subroutine symanal(istrsym,istr,nex,n,mosym,nstr,istrsym2)
************************************************************************
* This subroutine analyses the symmetry of operator strings *
************************************************************************
#include "MRCCCOMMON"
integer nex,n,istrsym(*),istr(nex,2*n),mosym(*),nstr(*)
integer iisym,i,j,isum,istrsym2(*)
maxmem=max(maxmem,2.d0*dble(n+nex*n))
C Symmetry of strings
do i=1,n
iisym=1
do j=1,nex
iisym=mult(iisym,mosym(istr(j,i)))
enddo
istrsym(i)=iisym
nstr(iisym)=nstr(iisym)+1
enddo
C Symmetry order
do i=1,n
do j=1,nex
istr(j,n+i)=istr(j,i)
enddo
enddo
isum=0
do ir=1,nir
do i=1,n
if(istrsym(i).eq.ir) then
isum=isum+1
do j=1,nex
istr(j,isum)=istr(j,n+i)
enddo
endif
enddo
enddo
do i=1,n
istrsym2(i)=istrsym(i)
enddo
do ir=1,nir
isum=0
do i=1,n
if(istrsym2(i).eq.ir) then
isum=isum+1
istrsym(i)=isum
endif
enddo
enddo
C
nr5=nir**5
nr4=nir**4
nr3=nir**3
nr2=nir**2
nrr=-nr3-nr2-nir
C
return
end
C
************************************************************************
subroutine opstr1(norb,nex,iucv,idrt,istr)
************************************************************************
* This subroutine generates operator strings *
************************************************************************
#include "MRCCCOMMON"
integer norb,nex,iucv(*),i,idrt(0:norb,0:nex),istr(nex,*)
C
if(iivirt.eq.nex) then
do i=1,nex
istr(i,istrv)=iucv(i)
enddo
else
iiocc=iiocc+1
if(idrt(iiocc,iivirt).gt.0) then
call opstr2(norb,nex,iucv,idrt,istr)
endif
if(idrt(iiocc,iivirt+1).gt.0) then
iivirt=iivirt+1
iucv(iivirt)=iiocc
istrv=istrv+idrt(iiocc-1,iivirt)
call opstr2(norb,nex,iucv,idrt,istr)
istrv=istrv-idrt(iiocc-1,iivirt)
iivirt=iivirt-1
endif
iiocc=iiocc-1
endif
C
return
end
C
************************************************************************
subroutine opstr2(norb,nex,iucv,idrt,istr)
************************************************************************
* This subroutine generates operator strings *
************************************************************************
#include "MRCCCOMMON"
integer norb,nex,iucv,idrt,istr
C
call opstr1(norb,nex,iucv,idrt,istr)
C
return
end
C
************************************************************************
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)
************************************************************************
* 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 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 iactv,iacto,iactva,iactoa,iactvb,iactob,noccup(*)
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
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,dbladd
real*8 f(*),faaal(*),faabe(*),fiial(*),fiibe(*)
real*8 hval(nbmax),kmal(nbmax,nbmax),jmat(nbmax,nbmax)
real*8 hvbe(nbmax),kmbe(nbmax,nbmax)
C
rewind(scrfile1)
C Loop over excitations
ssym=isym
if(calc.eq.1) ssym=1
if(ptroute) ssym=ptsym(1)
do nex=1,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
nn=nconf(iactva,iactvb,iactoa,iactob,i1,nex)
if(nn.gt.0) then
nampvirtal=i1
nampvirtbe=nex-nampvirtal
nampoccal=i1
nampoccbe=nex-nampoccal
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)
C
if(calc.eq.0) maxmem=max(maxmem,dble(ioffs-imem)+dble(nn)+
$dble(dbladd(n1*nampvirtal+n2*nampvirtbe+n3*nampoccal+n4*nampoccbe)
$)+dble(nbmax*n2*n4))
C
ii=ii+n3*n4*k
enddo
enddo
enddo
endif
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine ccit(nstr,nconf,trec,intrec,imed,v,nmax,istr,isympv,
$isympo,ita,tarec,iwa,isa,isw,wspcb,wspca,nactm,intn,ntcase,vmult,
$earec,econf,erec,wsmax,denrec,trec2,erec2,tcase,itypa,itypb,
$imedsyma,imedsymb,l2map,ioffs)
************************************************************************
* This subroutine controls CC iteration *
************************************************************************
#include "MRCCCOMMON"
integer wsmax,nmax,nstr,trec,wspcb(0:wsmax,1),wspca(0:wsmax,1),jj
integer nactm,imed,ii,k,istr,isympv,isympo,ita,tarec,iwa,j,tcase
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ioffs
integer intn(*),ntcase,vmult(*),earec,econf,erec,isw,i,intrec(*)
integer denrec,trec2,erec2,itypa,itypb,imedsyma(*),imedsymb(*)
integer l2map,isa(*)
real*8 v(*),norm,test,sum,eene
logical lll,lscr,log1,log2
C
if(ptroute.or.eomgrad) then
open(t1file,status='unknown',access='direct',recl=irecln)
open(l2file,status='unknown',access='direct',recl=irecln)
endif
if(hess) then
open(sdfile ,status='unknown',access='direct',recl=irecln)
endif
if(left.and.nroot.gt.1.and.isym.ne.1) then
open(strfile,status='unknown',access='direct',recl=irecln)
endif
irecmax=0
do iroot=1,nroot
if(eomgrad.and.osc.and.iroot.eq.3) call mrccend(1)
leftroute=.false.
multip=vmult(iroot)
test=1.d0
if(ptroute) then
C <K|exp(-T) dH/dx exp(T)|0> term
write(iout,*)
if(calc.eq.1) then
write(iout,*)
$'Calculation of the <K|exp(-T) dH/dx exp(T)|0> term...'
else
write(iout,*)
$'Calculation of the <K|dH/dx C|0> - dE/Dx <K|C|0> term...'
endif
call flush(iout)
t1route=.true.
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax,
$norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
t1route=.false.
if(calc.eq.1) then
C <0|(1+L) exp(-T) dH/dx exp(T)|K> term
write(iout,*)
write(iout,*)
$'Calculation of the <0|(1+L) exp(-T) dH/dx exp(T)|K> term...'
call flush(iout)
leftroute=.true.
l2route=.true.
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax
$,norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
leftroute=.false.
l2route=.false.
C
eomroute=.true.
t2route=.true.
endif
endif
lscr=.true.
2018 continue
nit=1
write(iout,*)
write(iout,122)
call flush(iout)
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax,
$norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
norm=dsqrt(norm)
call flush(iout)
122 format(1x,70('='))
C Construct intermediates for perturbation corrections
if(pert.gt.0.and.pert.le.3.and..not.leftroute) then
pertroute=1
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax,
$norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
pertroute=0
endif
C Perturbation corrections
if((pert.gt.0.and.pert.le.2).or.(pert.eq.3.and.leftroute)) then
write(iout,*)
if(pert.eq.3) then
leftroute=.false.
eomroute=.false.
endif
call pertcorr(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax,
$norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
endif
C Calculate constant terms
if((ptroute.and.lscr).or.(eomgrad.and.leftroute.and.iroot.gt.1
$.and.(.not.zroute).and.(.not.osc)).or.(eomgrad.and.osc.and.
$iroot.gt.1.and.eomroute)) then
if(eomgrad.and.iroot.gt.2)
$call wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
eomroute=.false.
t2route=.false.
if(calc.eq.1) then
C <0|(1+L)[exp(-T) H exp(T), dT/dx]|K> term
write(iout,*)
write(iout,*)
$'Calculation of the <0|(1+L)[exp(-T) H exp(T), dT/dx]|K> term...'
call flush(iout)
leftroute=.true.
l3route=.true.
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax
$,norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
l3route=.false.
C <0|(1+L)[exp(-T) {p+q-} exp(T), dT/dx]|0> term
write(iout,*)
write(iout,*)
$'Calculation of the <0|(1+L)[exp(-T){p+q-}exp(T), dT/dx]|0> term..
$.'
call flush(iout)
densroute=.true.
d2route=.true.
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax
$,norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
d2route=.false.
densroute=.false.
if(ptroute) then
l1route=.true.
endif
if(eomgrad) then
if(isym.ne.1.and.(.not.osc)) then
close(strfile,status='delete')
endif
call wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
if(iroot.eq.nroot) close(wsfile,status='delete')
if(osc) then
leftroute=.true.
else
zroute=.true.
endif
endif
lscr=.false.
goto 2018
endif
endif
l1route=.false.
C Calculate density-matrix
if(dens.gt.0.and.(((.not.lscr).and.(iroot.gt.1.or.
$(.not.eomgrad))).or.calc.eq.0).and.(iroot.gt.1.or.(.not.osc)))
$then
write(iout,*)
call flush(iout)
densroute=.true.
log1=eomroute
log2=leftroute
eomroute=.false.
leftroute=.true.
if(eomgrad.and.((.not.osc).or.zroute)) d1route=.true.
zroute=.false.
if(ptroute.and.calc.ne.0) then
d1route=.true.
ii=2*(naopt-1)*trecmax+trecmax+1
endif
call step(nstr,nconf,trec,intrec,imed,icore(icoupmem),v,nmax,
$norm,isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,
$wsmax,intn,nactmax,ntcase,earec,econf,erec,denrec,trec2,erec2
$,tcase,itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
densroute=.false.
eomroute=log1
leftroute=log2
write(iout,*)
if(eomgrad.and.osc.and.(.not.d1route)) then
zroute=.true.
goto 2018
endif
d1route=.false.
if(fdpci) then
fdpci=.false.
ptfreq=-ptfreq
goto 2018
endif
endif
if(ptroute.and.calc.eq.1) then
return
endif
if(eomcc.and.(lscr.or.(.not.hess))) then
eomroute=.true.
endif
if(left.and.lscr.and.(iroot.gt.1.or.(.not.eomgrad).or.osc))then
C Left-hand EOMCC (or Lambda) calculations
leftroute=.true.
eomroute=.false.
lscr=.false.
goto 2018
endif
irecmax=irecmax+trecmax
C Retrieve non-tot.symmetric spin-cases for higher lh EOMCC roots
if(nroot.gt.1.and.iroot.eq.1.and.((calc.eq.0.and.isy1.ne.isy2)
$.or.(left.and.(isym.ne.1.or.eomgrad).and.(.not.osc)))) then
call wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
if(calc.eq.0) then
isym=isy2
ptsh=nimed2
endif
endif
enddo
130 format(' Total CI energy /au/: ',25f18.12)
121 format(' Total CC energy /au/: ',25f18.12)
133 format(' Total EOMCC energy /au/: ',25f18.12)
132 format(' Excitation energy /eV/: ',25f9.3)
127 format(' Total LCC energy /au/: ',25f14.8)
126 format(' Total MP(2) energy /au/: ',25f14.8)
128 format(' Correlation energy /au/: ',25f14.8)
134 format(' Derivative of correlation energy /au/:',25f14.8)
C
return
end
C
************************************************************************
subroutine step(nstr,nconf,trec,intrec,imed,icmem,v,nmax,norm,
$isympv,isympo,ita,tarec,iwa,isa,isw,istr,test,wspcb,wspca,wsmax,
$intn,nactm,ntcase,earec,econf,erec,denrec,trec2,erec2,tcase,
$itypa,itypb,imedsyma,imedsymb,l2map,ioffs)
************************************************************************
* This subroutine performs one CC iteration step *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,ita(*),istr,nstr,i,k,l,m,n,icmem,ii,jj,l2map
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ioffs
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ic1,j
integer econf(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),file1
integer intrec(*),isympv,isympo,isw,iwa(*),file,wsmax,intn(*),isa
integer iadd0(0:wsmax),iadd1(0:wsmax),nn,iactvb,iactob,ntcase(0:1)
integer iadd2(0:wsmax),iadd3(0:wsmax),iactv,iacto,iactva,iactoa
integer n1,n2,n3,n4,nn1,nn2,nn3,nn4,imed,ifi,earec,wspc3(0:wsmax)
integer nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1,nointold1
integer nvirtold1,noccold1,namp1,ical1,isig11,isave1,nptype
integer nvintnew2,nointnew2,nvirtnew2,noccnew2,nvintold2,nointold2
integer nvirtold2,noccold2,namp2,ical2,isig12,isave2,nlist
integer nvintnew3,nointnew3,nvirtnew3,noccnew3,nvintold3,nointold3
integer nvirtold3,noccold3,namp3,ical3,isig13,isave3,ssym
integer nvintnew4,nointnew4,nvirtnew4,noccnew4,nvintold4,nointold4
integer nvirtold4,noccold4,namp4,ical4,isig14,isave4,imedsyma(*)
integer nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,iadds(0:wsmax)
integer nvintold1a,nointold1a,nvirtold1a,noccold1a,itypb(*),icl1
integer nvintnew2a,nointnew2a,nvirtnew2a,noccnew2a,file2
integer nvintold2a,nointold2a,nvirtold2a,noccold2a,itypa(*)
integer nvintnew3a,nointnew3a,nvirtnew3a,noccnew3a,imedsymb(*)
integer nvintold3a,nointold3a,nvirtold3a,noccold3a
integer nvintnew4a,nointnew4a,nvirtnew4a,noccnew4a
integer nvintold4a,nointold4a,nvirtold4a,noccold4a
integer wspcb(0:wsmax,13),wspca(0:wsmax,13),tarec,denrec(*),drec
integer trec2((2*nactm+1)**2,0:nmax),erec2((2*nactm+1)**2,0:nmax)
integer tcase(2,(2*nactm+1)**2,0:nmax)
real*8 v(*),norm,sum,c0,c1,denom,e1,chc(2,2),a1,test,sum1
character*5 eomcheck
logical lll,llog,log1,log2,log3,log4,ldir
C File structure
lamfile=cfile
ampfile=tfile
if(calc.eq.0.and.ptroute.and.(.not.t1route).and.(.not.densroute))
$ampfile=ptfile
if(leftroute.or.densroute) then
if(ptroute) then
if(l1route.or.d1route) lamfile=plfile
if(calc.eq.0) lamfile=ptfile
else
lamfile=tfile
ampfile=cfile
if(osc.and.(l3route.or.d2route)) lamfile=lfile
endif
endif
c0=1.d0
C Symmetries of lambda amplitudes
j=isym
if(l1route.or.d1route.or.(ptroute.and.calc.eq.0)) j=ptsym(1)
if(calc.eq.1.and.(iroot.eq.1.or.(zroute.and.(.not.osc)).or.
$(d1route.and.(.not.osc)).or.((l3route.or.d2route).and.osc)).and.
$eomgrad) j=1
if(calc.eq.0.and.osc) j=isy2
do i=-op,0
imedsyma(i)=j
imedsymb(i)=j
enddo
C Initialize scratch file
sum=0.d0
ic1=calc
if((calc.eq.0.or.eomroute).and.(.not.densroute).and.
$((.not.ptroute).or.calc.ne.1)) then
if(nit.eq.1.or.diag.eq.0.or.diag.eq.2) then
file=tfile
else
file=diisfile
endif
if(ptroute) file=ampfile
if((econf(0,0,0,0,0,0).eq.1.and.(calc.ne.0.or.isy1.eq.isy2.or.
$isym.eq.isy2)).or.(nconf(0,0,0,0,0,0).eq.1.and.calc.eq.0.and.
$isy1.ne.isy2.and.isym.eq.isy1).or.(t1route.and.ptsym(1).eq.1))
$then
if(t1route.and.nconf(0,0,0,0,0,0).eq.0) c0=0.d0
else
c0=0.d0
endif
else
file=tfile
if((ptroute.and.leftroute).or.(osc.and.l3route)) file=lamfile
ibuf(1)=1.d0
c0=1.d0
if(densroute.and.calc.eq.0) then
if(ptroute.and.nconf(0,0,0,0,0,0).eq.0) c0=0.d0
endif
if(leftroute) then
ibuf(1)=0.d0
if((iroot.gt.1.and.((.not.zroute).or.osc).and.(.not.d1route)
$.and.(.not.densroute).and.(.not.l3route)).or.l1route) then
c0=0.d0
ic1=0
endif
if(d1route.and.((.not.eomgrad).or.osc)) c0=0.d0
if(eomgrad.and.(l3route.or.d2route).and.(.not.osc)) c0=0.d0
endif
endif
c1=c0
ibuf(1)=0.d0
k=trecmax
nn=op
if(pert.ne.0) then
if(pertroute.eq.0.and.pert.le.3) nn=op-1
if(pertroute.gt.2.and.pertroute.le.op+2) nn=pertroute-4
if(pertroute.gt.op+2) nn=pertroute-4-op
endif
if(sacc.lt.0) nn=op !PT
if(sacc.gt.0) nn=2 !Miriam
if(densroute) then
k=drecmax
nn=iimed
endif
C Reading formula tape
llog=((.not.leftroute).or.(iroot.eq.1.and.nit.eq.1).or.l2route.or.
$l3route.or.(eomgrad.and.iroot.ne.1.and.nit.eq.1.and.(.not.zroute)
$.and.(.not.osc))).and.(.not.l1route)
rewind(ftfile)
if(eomroute.and.calc.ne.0 ) call readlab('EOMCC')
if(leftroute.and.calc.ne.0) call readlab('LEFT ')
if(densroute) call readlab('DENS ')
if(l3route) call readlab('DLEFT')
if(d2route) call readlab('DDENS')
do i=1,nn
if(pert.ge.4.and.i.eq.op) pertroute=1
if(ptroute.or.eomroute.or.(leftroute.and.iroot.gt.1.and.
$(((.not.zroute).and.(.not.d1route).and.(.not.l3route)).or.osc)))
$then
ifi=ierecmem+i*((nactm+1)**4)*(nmax+1)
else
ifi=itrecmem+i*((nactm+1)**4)*(nmax+1)
endif
if(densroute) then
ii=1
drec=denrec(i)
else
ii=ntcase(i)
endif
do nptype=1,ii
if((pertroute.eq.1.and.i.ne.op).or.
$ (pertroute.gt.2.and.i.lt.nn)) then
nn1=0
call skip1(nn1,ftfile)
else
if(.not.densroute) then
if(ptroute.or.eomroute.or.(leftroute.and.iroot.gt.1.and.
$(((.not.zroute).and.(.not.l3route)).or.osc)).or.(calc.eq.0.and.
$isy1.ne.isy2.and.isym.eq.isy2)) then
drec=erec2(nptype,i)
else
drec=trec2(nptype,i)
endif
endif
read(ftfile,*) n1
call ifillzero(iadd0,wsmax+1)
if((ptroute.and.(eomroute.or.l1route.or.d1route.or.(calc.eq.0
$.and.(.not.t1route).and.(.not.densroute)))).or.(osc.and.(zroute
$.or.d1route))) call iaddini(econf,nmax,nactm,i,iadd0,tcase,nptype,
$wspca,wsmax)
if((.not.osc).and.(zroute.or.(eomgrad.and.d1route)))
$call iaddini(nconf,nmax,nactm,i,iadd0,tcase,nptype,wspca,wsmax)
do nn1=1,n1
call readft(nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1,
$nointold1,nvirtold1,noccold1,namp1,ical1,isig11,nlist,isave1,
$nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,nvintold1a,nointold1a,
$nvirtold1a,noccold1a,itypa,1,l2map,log1,intn)
ldir=.true.
if(densroute) then
icl1=i+ptsh
else
icl1=ical1
endif
if(isave1.gt.0.and.(llog.or.itypa(ical1).ne.0).and.log1) then
read(ftfile,*) n2
call ifillzero(iadd1,wsmax+1)
do nn2=1,n2
call readft(nvintnew2,nointnew2,nvirtnew2,noccnew2,
$nvintold2,nointold2,nvirtold2,noccold2,namp2,ical2,isig12,nlist,
$isave2,nvintnew2a,nointnew2a,nvirtnew2a,noccnew2a,nvintold2a,
$nointold2a,nvirtold2a,noccold2a,itypa,2,l2map,log2,intn)
if(isave2.eq.0.and.nn2.eq.1) ldir=.false.
if(isave2.gt.0.and.(llog.or.itypa(ical2).ne.0).and.log2)
$then
read(ftfile,*) n3
call ifillzero(iadd2,wsmax+1)
do nn3=1,n3
call readft(nvintnew3,nointnew3,nvirtnew3,noccnew3,
$nvintold3,nointold3,nvirtold3,noccold3,namp3,ical3,isig13,nlist,
$isave3,nvintnew3a,nointnew3a,nvirtnew3a,noccnew3a,nvintold3a,
$nointold3a,nvirtold3a,noccold3a,itypa,3,l2map,log3,intn)
if(isave3.gt.0.and.(llog.or.itypa(ical3).ne.0).and.
$log3) then
read(ftfile,*) n4
call ifillzero(iadd3,wsmax+1)
do nn4=1,n4
call readft(nvintnew4,nointnew4,nvirtnew4,noccnew4
$,nvintold4,nointold4,nvirtold4,noccold4,namp4,ical4,isig14,nlist,
$isave4,nvintnew4a,nointnew4a,nvirtnew4a,noccnew4a,nvintold4a,
$nointold4a,nvirtold4a,noccold4a,itypa,4,l2map,log4,intn)
if(isave4.gt.0.and.(llog.or.itypa(ical4).ne.0)
$.and.log4) then
call conin1(trec,nstr,nmax,icmem,nvintnew4,nointnew4,nvirtnew4,
$noccnew4,nvintold4,nointold4,nvirtold4,noccold4,namp4,isig14,
$scrfile4,0,intfile,intrec(nlist),iadd3,iadd3,.true.,v,isave4,
$intrec(ical4),isympv,isympo,nir,nconf,ita,tarec,iwa,isa,isw,
$iwa(nir**6+1),file,wspcb(0,ical3),wspca(0,nlist),wspca(0,ical4),
$wsmax,nvintnew4a,nointnew4a,nvirtnew4a,noccnew4a,nvintold4a,
$nointold4a,nvirtold4a,noccold4a,nactmax,earec,erec,econf,4,imed,
$wspca,intrec,itypa(ical4),itypb(ical4),imedsyma(ical4),
$imedsymb(ical4),ioffs)
else
call imedret(scrfile4,0,wspcb(0,ical3),intfile,intrec(ical4),
$isig14,iadd3,v,nmax,2,1,c0,wspcb(0,ical3),wspca(0,ical4),nactmax,
$0,0,ical4,nconf,trec,econf,erec,nvintnew4a+nvirtnew4a,
$nointnew4a+noccnew4a,wsmax,4,isave4,itypa(ical4),imedsyma(ical4))
endif
enddo
call conin1(trec,nstr,nmax,icmem,nvintnew3,nointnew3,nvirtnew3,
$noccnew3,nvintold3,nointold3,nvirtold3,noccold3,namp3,isig13,
$scrfile3,0,scrfile4,0,iadd2,iadd3,.false.,v,isave3,intrec(ical3),
$isympv,isympo,nir,nconf,ita,tarec,iwa,isa,isw,iwa(nir**6+1),file,
$wspcb(0,ical2),wspcb(0,ical3),wspca(0,ical3),wsmax,nvintnew3a,
$nointnew3a,nvirtnew3a,noccnew3a,nvintold3a,nointold3a,nvirtold3a,
$noccold3a,nactmax,earec,erec,econf,3,imed,wspca,intrec,
$itypa(ical3),itypb(ical3),imedsyma(ical3),imedsymb(ical3),ioffs)
else
call imedret(scrfile3,0,wspcb(0,ical2),intfile,intrec(ical3)
$,isig13,iadd2,v,nmax,2,1,c0,wspcb(0,ical2),wspca(0,ical3),nactmax,
$0,0,ical3,nconf,trec,econf,erec,nvintnew3a+nvirtnew3a,
$nointnew3a+noccnew3a,wsmax,3,isave3,itypa(ical3),imedsyma(ical3))
endif
enddo
if(nvintnew2+nvirtnew2.ne.op.or..not.ldir.or.pertroute.ne.1)
$call conin1(trec,nstr,nmax,icmem,nvintnew2,nointnew2,nvirtnew2,
$noccnew2,nvintold2,nointold2,nvirtold2,noccold2,namp2,isig12,
$scrfile2,0,scrfile3,0,iadd1,iadd2,.false.,v,isave2,intrec(ical2),
$isympv,isympo,nir,nconf,ita,tarec,iwa,isa,isw,iwa(nir**6+1),file,
$wspcb(0,ical1),wspcb(0,ical2),wspca(0,ical2),wsmax,nvintnew2a,
$nointnew2a,nvirtnew2a,noccnew2a,nvintold2a,nointold2a,nvirtold2a,
$noccold2a,nactmax,earec,erec,econf,2,imed,wspca,intrec,
$itypa(ical2),itypb(ical2),imedsyma(ical2),imedsymb(ical2),ioffs)
else
call imedret(scrfile2,0,wspcb(0,ical1),intfile,intrec(ical2)
$,isig12,iadd1,v,nmax,2,1,c0,wspcb(0,ical1),wspca(0,ical2),nactmax,
$0,0,ical2,nconf,trec,econf,erec,nvintnew2a+nvirtnew2a,
$nointnew2a+noccnew2a,wsmax,2,isave2,itypa(ical2),imedsyma(ical2))
endif
enddo
if(pertroute.ne.1.and.(pert.lt.4.or.namp1.ne.op))
$call conin1(trec,nstr,nmax,icmem,nvintnew1,nointnew1,nvirtnew1,
$noccnew1,nvintold1,nointold1,nvirtold1,noccold1,namp1,isig11,
$scrfile1,drec,scrfile2,0,iadd0,iadd1,.false.,v,isave1,
$intrec(ical1),isympv,isympo,nir,nconf,ita,tarec,iwa,isa,isw,
$iwa(nir**6+1),file,wspca(0,icl1),wspcb(0,ical1),wspca(0,ical1),
$wsmax,nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,nvintold1a,
$nointold1a,nvirtold1a,noccold1a,nactmax,earec,erec,econf,1,imed,
$wspca,intrec,itypa(ical1),itypb(ical1),imedsyma(ical1),
$imedsymb(ical1),ioffs)
else
if(leftroute) then
if(densroute) then
call imedret(scrfile1,drec,icore(ifi),intfile,intrec(ical1),isig11
$,iadd0,v,nmax,2,ic1,c1,wspca(0,icl1),wspca(0,ical1),nactmax,
c $,iadd0,v,nmax,2,ic1,c1,wspca(0,i+ptsh),wspca(0,ical1),nactmax,
$nvirtnew1a,noccnew1a,ical1,nconf,trec,econf,erec,nvintnew1a+
$nvirtnew1a,nointnew1a+noccnew1a,wsmax,1,isave1,itypa(ical1)
$,imedsyma(ical1))
else
if((iroot.eq.1.or.(zroute.and.(.not.osc)).or.(l3route
$.and.osc)).and.(itypa(ical1).eq.0.or.itypa(ical1).eq.3)) then
call rospc(wspca(0,ical1),wspc3,nvirtold1a,noccold1a,wsmax)
call imedret(scrfile1,0,icore(ifi),intfile,intrec(ical1),isig11,
$iadd0,v,nmax,0,ic1,c0,trec,wspc3,nactmax,nvintnew1a,
$nointnew1a,ical1,nconf,trec,econf,erec,nvintnew1a+nvirtnew1a,
$nointnew1a+noccnew1a,wsmax,1,isave1,itypa(ical1),imedsyma(ical1))
endif
if(itypa(ical1).ne.0.and.itypa(ical1).ne.3) then
call ifillzero(iadds,wsmax+1)
call imedret(scrfile4,0,icore(ifi),intfile,intrec(ical1),isig11
$,iadds,v,nmax,2,1,c0,wspca(0,ical1),wspca(0,ical1),nactmax,
$nvirtnew1a,noccnew1a,ical1,nconf,trec,econf,erec,nvintnew1a+
$nvirtnew1a,nointnew1a+noccnew1a,wsmax,1,isave1,itypa(ical1)
$,imedsyma(ical1))
call anti1(nconf,econf,trec,erec,nstr,nmax,v,nir,isympv,isympo,
$ita,tarec,earec,nactm,icmem,wspca(0,ical1),nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvirtnew1a,noccnew1a,iadd0,iadds)
endif
endif
else
if((pert.eq.0).or.log1) then
call imedret(scrfile1,0,icore(ifi),intfile,intrec(ical1),isig11,
$iadd0,v,nmax,0,ic1,c0,trec,wspca(0,ical1),nactmax,nvirtnew1a,
$noccnew1a,ical1,nconf,trec,econf,erec,nvintnew1a+nvirtnew1a,
$nointnew1a+noccnew1a,wsmax,1,isave1,itypa(ical1),imedsyma(ical1))
endif
endif
endif
enddo
endif
read(ftfile,*)
enddo
enddo
#if defined (MPI)
call psend(scrfile1,nconf,nmax,nactm,trec,v,0,wspca,wsmax,denrec,
$econf,erec)
if(d1route.and.dboc.eq.1.and.calc.eq.1.and..not.dbroute)
$call psend(t1file,nconf,nmax,nactm,trec,v,0,wspca,wsmax,denrec,
$econf,erec)
#endif
C Perturbative methods
if(pert.ge.4) then
pertroute=2
call coninp(trec,nstr,nmax,icore(icoupmem),v(nbasis+1),isympv,
$isympo,nir,nconf,ita,tarec,iwa,isa,iwa(nir**6+1),tfile,wspca,
$wspcb,nactmax,imedsyma,imedsymb,icore(istrecmem),icore(itransmem),
$icore(iarcmem),wsmax,intrec,ntcase,tcase,v,v(nvirtal+1),
$v(nvirt+1),v(nvirt+nal+1),sum,sum,nvirtal,nvirtbe,ioffs+nbasis+1)
pertroute=0
endif
C
if(pertroute.ne.0.or.t1route.or.l2route.or.l3route) return
if(densroute) return
C
norm=0.d0
if((calc.eq.0.or.eomroute).and.diag.eq.1.and.(.not.ptroute)) then
else
if((calc.ne.0.and.(.not.eomroute)).or.ptroute) then
C DIIS extrapolation
ndit=nit-ifirstdiis+1
c ndit=-1
if(ndit.gt.0.and.dsqrt(norm).gt.1.d0*tol) then
if(ptroute.or.(leftroute.and.iroot.gt.1.and.((.not.zroute).or.osc
$))) then
call diis(dcore(idiismem),dcore(idiismem+ndiis**2+1),
$dcore(idiismem+2*ndiis**2+1),dcore(idiismem+3*ndiis**2+1),ndiis,
$v,econf,erec,nmax,nactmax,test,file2)
else
call diis(dcore(idiismem),dcore(idiismem+ndiis**2+1),
$dcore(idiismem+2*ndiis**2+1),dcore(idiismem+3*ndiis**2+1),ndiis,
$v,nconf,trec,nmax,nactmax,test,file2)
endif
endif
else
C Davidson diagonalization
if(calc.eq.0.and.isy1.ne.isy2.and.isym.eq.isy1) then
call david(dcore(idiismem),dcore(idiismem+ndiis**2+1),
$dcore(idiismem+2*ndiis**2+1),dcore(idiismem+3*ndiis**2+1),ndiis,
$v,nconf,trec,nmax,nactmax,test,norm,file2)
else
call david(dcore(idiismem),dcore(idiismem+ndiis**2+1),
$dcore(idiismem+2*ndiis**2+1),dcore(idiismem+3*ndiis**2+1),ndiis,
$v,econf,erec,nmax,nactmax,test,norm,file2)
endif
endif
endif
if(calc.eq.0.or.eomroute.or.(leftroute.and.iroot.gt.1.and.
$((.not.zroute).or.osc)).or.ptroute) then
ssym=isym
file1=tfile
if(ptroute) then
ssym=ptsym(1)
if(t2route.or.calc.eq.0) file1=ptfile
if(l1route) file1=plfile
endif
if(calc.eq.0.and.isy1.ne.isy2.and.isym.eq.isy1) then
call symmetrize(nconf,trec,nstr,nmax,nir,isympv,isympo,ita,
$tarec,v,file1,nactmax,ssym) !szemet
else
call symmetrize(econf,erec,nstr,nmax,nir,isympv,isympo,ita,
$earec,v,file1,nactmax,ssym) !szemet
endif
else
call symmetrize(nconf,trec,nstr,nmax,nir,isympv,isympo,ita,tarec
$,v,tfile,nactmax,1)
endif
C
return
end
C
************************************************************************
subroutine gcoupl(icmem,nvstral,nvstrbe,nostral,nostrbe,iarc,nmax,
$nnir,ivstral,ivstrbe,iostral,iostrbe,nactm,iscr,ivrecal,ivrecbe,
$iorecal,iorecbe,itrva,itrvb,itroa,itrob)
************************************************************************
* This subroutine generates coupling coefficient for excitation *
* operator strings *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,iarc(0:nactm,0:nmax,4),ii,i,j,k,minvo
integer nvstral(nnir,0:nactm,0:nmax),nvstrbe(nnir,0:nactm,0:nmax)
integer nostral(nnir,0:nactm,0:nmax),nostrbe(nnir,0:nactm,0:nmax)
integer icmem(nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,4),ir1,ir2
integer ivstral(nnir,0:nactm,0:nmax),ivstrbe(nnir,0:nactm,0:nmax)
integer iostral(nnir,0:nactm,0:nmax),iostrbe(nnir,0:nactm,0:nmax)
integer iact,iact1,iact2,iscr(*)
integer ivrecal(nnir,0:nactm,0:nmax),ivrecbe(nnir,0:nactm,0:nmax)
integer iorecal(nnir,0:nactm,0:nmax),iorecbe(nnir,0:nactm,0:nmax)
integer itrva(2,0:nactm,0:nmax),itrvb(2,0:nactm,0:nmax)
integer itroa(2,0:nactm,0:nmax),itrob(2,0:nactm,0:nmax)
integer n1,n2,n3,nn1,nn2
C
write(iout,*) 'Calculation of coupling coefficients...'
minvo=2
C Virtual alpha strings
if(calc.eq.1.or.dens.gt.0) minvo=minvirtal
do i=0,minvo
do j=0,minvirtal
if(i+j.le.minvirtal2) then
i1=i+j
i2=i1-1
do iact1=max(0,i-mrop2),min(nactva,i)
do iact2=max(0,j-mrop2),min(nactva,j)
iact=iact1+iact2
if(iact.le.nactva.and.i1-iact.le.mrop2.and.
$(pert.eq.0.or.i1.lt.op.or.iact.ge.1)) then
do ir1=1,nir
do ir2=1,nir
n1=nvstral(ir1,iact1,i)
n2=nvstral(ir2,iact2,j)
n3=itrva(2,iact,i1)
nn1=n1*i
nn2=n2*j
call getint(strfile,ivrecal(ir1,iact1,i),iscr(1),nn1)
call getint(strfile,ivrecal(ir2,iact2,j),iscr(nn1+1),nn2)
call getint(strfile,itrva(1,iact,i1),iscr(nn1+nn2+1),n3)
call gcoupl1(icmem(ir1,ir2,iact1,iact2,i,j,1),n1,n2,nvirtal,
$icore(iarc(iact,i1,1)),i,j,iscr(1),iscr(nn1+1),iscr(nn1+nn2+1),
$icore(iindmem),icore(iucvmem),iscr(nn1+nn2+n3+1),ii)
maxmem=max(maxmem,dble(nn1+nn2+n3+ii))
enddo
enddo
endif
enddo
enddo
endif
enddo
enddo
C Virtual beta strings
if(calc.eq.1.or.dens.gt.0) minvo=minvirtbe
if(cs.and.nvirtal.eq.nvirtbe.and.minvirtal.eq.minvirtbe.and.
$nactva.eq.nactvb) then
do i=0,minvo
do j=0,minvirtbe
if(i+j.le.minvirtbe2) then
do iact1=max(0,i-mrop2),min(nactvb,i)
do iact2=max(0,j-mrop2),min(nactvb,j)
do ir1=1,nir
do ir2=1,nir
icmem(ir1,ir2,iact1,iact2,i,j,2)=
$ icmem(ir1,ir2,iact1,iact2,i,j,1)
enddo
enddo
enddo
enddo
endif
enddo
enddo
else
do i=0,minvo
do j=0,minvirtbe
if(i+j.le.minvirtbe2) then
i1=i+j
i2=i1-1
do iact1=max(0,i-mrop2),min(nactvb,i)
do iact2=max(0,j-mrop2),min(nactvb,j)
iact=iact1+iact2
if(iact.le.nactvb.and.i1-iact.le.mrop2.and.
$(pert.eq.0.or.i1.lt.op.or.iact.ge.1)) then
do ir1=1,nir
do ir2=1,nir
n1=nvstrbe(ir1,iact1,i)
n2=nvstrbe(ir2,iact2,j)
n3=itrvb(2,iact,i1)
nn1=n1*i
nn2=n2*j
call getint(strfile,ivrecbe(ir1,iact1,i),iscr(1),nn1)
call getint(strfile,ivrecbe(ir2,iact2,j),iscr(nn1+1),nn2)
call getint(strfile,itrvb(1,iact,i1),iscr(nn1+nn2+1),n3)
call gcoupl1(icmem(ir1,ir2,iact1,iact2,i,j,2),n1,n2,nvirtbe,
$icore(iarc(iact,i1,2)),i,j,iscr(1),iscr(nn1+1),iscr(nn1+nn2+1),
$icore(iindmem),icore(iucvmem),iscr(nn1+nn2+n3+1),ii)
maxmem=max(maxmem,dble(nn1+nn2+n3+ii))
enddo
enddo
endif
enddo
enddo
endif
enddo
enddo
endif
C Occupied alpha strings
if(calc.eq.1.or.dens.gt.0) minvo=minoccal
do i=0,minvo
do j=0,minoccal
if(i+j.le.minoccal2) then
i1=i+j
i2=i1-1
do iact1=max(0,i-mrop2),min(nactoa,i)
do iact2=max(0,j-mrop2),min(nactoa,j)
iact=iact1+iact2
if(iact.le.nactoa.and.i1-iact.le.mrop2) then
do ir1=1,nir
do ir2=1,nir
n1=nostral(ir1,iact1,i)
n2=nostral(ir2,iact2,j)
n3=itroa(2,iact,i1)
nn1=n1*i
nn2=n2*j
call getint(strfile,iorecal(ir1,iact1,i),iscr(1),nn1)
call getint(strfile,iorecal(ir2,iact2,j),iscr(nn1+1),nn2)
call getint(strfile,itroa(1,iact,i1),iscr(nn1+nn2+1),n3)
call gcoupl1(icmem(ir1,ir2,iact1,iact2,i,j,3),n1,n2,nal,
$icore(iarc(iact,i1,3)),i,j,iscr(1),iscr(nn1+1),iscr(nn1+nn2+1),
$icore(iindmem),icore(iucvmem),iscr(nn1+nn2+n3+1),ii)
maxmem=max(maxmem,dble(nn1+nn2+n3+ii))
enddo
enddo
endif
enddo
enddo
endif
enddo
enddo
C Occupied beta strings
if(calc.eq.1.or.dens.gt.0) minvo=minoccbe
if(cs.and.nal.eq.nbe.and.minoccal.eq.minoccbe.and.
$nactoa.eq.nactob) then
do i=0,minvo
do j=0,minoccbe
if(i+j.le.minoccbe2) then
do iact1=max(0,i-mrop2),min(nactob,i)
do iact2=max(0,j-mrop2),min(nactob,j)
do ir1=1,nir
do ir2=1,nir
icmem(ir1,ir2,iact1,iact2,i,j,4)=
$ icmem(ir1,ir2,iact1,iact2,i,j,3)
enddo
enddo
enddo
enddo
endif
enddo
enddo
else
do i=0,minvo
do j=0,minoccbe
if(i+j.le.minoccbe2) then
i1=i+j
i2=i1-1
do iact1=max(0,i-mrop2),min(nactob,i)
do iact2=max(0,j-mrop2),min(nactob,j)
iact=iact1+iact2
if(iact.le.nactob.and.i1-iact.le.mrop2) then
do ir1=1,nir
do ir2=1,nir
n1=nostrbe(ir1,iact1,i)
n2=nostrbe(ir2,iact2,j)
n3=itrob(2,iact,i1)
nn1=n1*i
nn2=n2*j
call getint(strfile,iorecbe(ir1,iact1,i),iscr(1),nn1)
call getint(strfile,iorecbe(ir2,iact2,j),iscr(nn1+1),nn2)
call getint(strfile,itrob(1,iact,i1),iscr(nn1+nn2+1),n3)
call gcoupl1(icmem(ir1,ir2,iact1,iact2,i,j,4),n1,n2,nbe,
$icore(iarc(iact,i1,4)),i,j,iscr(1),iscr(nn1+1),iscr(nn1+nn2+1),
$icore(iindmem),icore(iucvmem),iscr(nn1+nn2+n3+1),ii)
maxmem=max(maxmem,dble(nn1+nn2+n3+ii))
enddo
enddo
endif
enddo
enddo
endif
enddo
enddo
endif
C
return
end
C
************************************************************************
subroutine gcoupl1(icmem,nstr1,nstr2,nv,iarc,nex1,nex2,iistr1,
$iistr2,itr,iind,iucv,coup,ii)
************************************************************************
* This subroutine generates coupling coefficient for excitation *
* operator strings *
************************************************************************
#include "MRCCCOMMON"
integer nstr1,nstr2,coup(*),nex1,nex2,nv,a,ii,isum,icmem
integer iarc(0:nv,0:nex1+nex2),iistr1(nex1,nstr1)
integer itr(*),iind(*),iucv(*),i,j,k,iistr2(nex2,nstr2)
logical lll
C
ii=2*nstr2
do istr2=1,nstr2
isum=0
coup(istr2)=ii
do i=1,nex2
iucv(i)=iistr2(i,istr2)
enddo
do istr1=1,nstr1
do i=1,nex2
iind(nex1+i)=iucv(i)
enddo
do i=1,nex1
iind(i)=iistr1(i,istr1)
enddo
isig1=1
1011 continue
lll=.false.
do i=1,i2
j=i+1
if(iind(i).gt.iind(j)) then
lll=.true.
isig1=-isig1
k=iind(j)
iind(j)=iind(i)
iind(i)=k
endif
enddo
if(lll) goto 1011
do i=1,i2
if(iind(i).eq.iind(i+1)) then
goto 1012
endif
enddo
istrv=1
do i=1,i1
istrv=istrv+iarc(iind(i)-1,i)
enddo
isum=isum+1
ii=ii+1
coup(ii)=istr1
ii=ii+1
coup(ii)=isig1*itr(istrv)
1012 continue
enddo
coup(nstr2+istr2)=isum
enddo
call getmem(ii,icore(1),icmem)
call icp(coup,icore(icmem),ii)
C
return
end
C
************************************************************************
subroutine diis(bmat,bvec,invbmat,cvec,ndmax,v,nconf,trec
$,nmax,nactm,test,file)
************************************************************************
* DIIS extrapolation *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,nact,nactpl,iact,ii,i,j,ntrecact,file
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer ifi,ndmax,iactv,iacto,iactva,iactoa,iactvb,iactob
integer ipiv(ndmax)
real*8 bmat(ndmax,ndmax),bvec(ndmax**2),invbmat(ndmax,ndmax)
real*8 cvec(ndmax),sum,sum1,v(*),norm,test
C
if(calc.eq.0.or.(leftroute.and.iroot.gt.1.and.isym.ne.1.and.
$((.not.zroute).or.osc)).or.(ptroute.and.ptsym(1).ne.1)) then
ifi=0
else
ifi=1
endif
nact=min(ndiis,ndit)
nactpl=mod(ndit,ndiis)
if(nactpl.eq.0) nactpl=ndiis
ntrecact=(nactpl-1)*trecmax
C Update the scalar product matrix (B matrix)
do i=ifi,op
do iactv=max(0,i-mrop),min(nactv,i)
do iacto=max(0,i-mrop),min(nacto,i)
do i1=0,i
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
ii=nconf(iactva,iactvb,iactoa,iactob,i1,i)
if(ii.gt.0) then
maxmem=max(maxmem,2.d0*dble(ii))
c if(2*ii.gt.maxcor) write(6,*) 'M2 ',8.d0*dble(2*ii)/dble(twoto20)
endif
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine conini(trec,nstr,nmax,icmem,nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvintold1,nointold1,nvirtold1,noccold1,namp1,
$isig11,file1,irec1,file2,irec2,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,nconf,ita,tarec,iwa,isa,isw,iwan,file3,wspc1,
$wspc2,wspca,wsmax,nvintnewact,nointnewact,nvirtnewact,noccnewact,
$nvintoldact,nointoldact,nvirtoldact,noccoldact,nactm,earec,erec,
$econf,itypa,itypb,nnewsym,noldsym,ilev,ioffs)
************************************************************************
* This subroutine initializes variables for contractions and performs *
* transpositions if it is necessary *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,i,j,k,m,iamprec,ii,irec3,ilev,ioffs
integer nstr(nnir,0:nactm,0:nmax,4),isa(*),n1,n2,n3,n4,nmem,wsmax
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),file3
integer icmem(nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,4),inewadd
integer nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1,nointold1
integer nvirtold1,noccold1,namp1,isig11,isig12,iadda(0:wsmax)
integer file1,file2,irec1,irec2,ircold,ircnew,isave,intrec,nsyma
integer iadd(0:1),ntoldlen,namplen,ia1,ia2,ntoldleno,file4
integer nvintoldal,nvintoldbe,nsymw,nointoldal,nointoldbe,isw(*)
integer nvirtoldbelen,noccoldbelen,incsum(nnir),incold(nnir)
integer nsumvirtallen,nsumvirtbelen,nsumoccallen,nsumoccbelen,jj
integer nvirtnewbelen,noccnewbelen,nampvirtlen,kk,ir1,isymi,isyma
integer nampocclen,ntampvirtal,ntampvirtbe,ntampoccal,ntampoccbe
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 ira,nvirtnewsym,noccnewsym,irav,nvirtnewsyma,nvirtnewsymb
integer irao,noccnewsyma,noccnewsymb,nsumsym,ntampsym,noldsym
integer ntampsymv,ntampsymo,irtv,ntampsymva,ntampsymvb,noccoldlen
integer nvirtoldsyma,nvirtoldsymb,nvirtoldsym,irto,itypa,itypb
integer ntampsymoa,ntampsymob,noccoldsyma,irt,nvirtoldlen
integer tarec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer ita(nnir,nnir,nnir,nnir),noccoldsymb,noccoldsym,ntnewlen
integer iwa(nnir,nnir,nnir,nnir,nnir,nnir),incnew(nnir),ioldadd
integer iwan(nnir,nnir,nnir,nnir,nnir),nampsym,j1,j2,j3,j4,n
integer wspc1(0:1),wspc2(0:1),wspca(0:1),iw,iaddo(0:1),nnewsym
integer nvintnewact,nointnewact,nvirtnewact,noccnewact
integer nvintoldact,nointoldact,nvirtoldact,noccoldact
integer nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact
integer nvirtoldalact,noccoldalact,nvirtoldbeact,noccoldbeact
integer nampvirtactmax,nampoccactmax,nampvirtactmin,nampoccactmin
integer ntampvirtact,ntampoccact,nsumvirtact,nsumoccact
integer nampvirtact,nampoccact,nampvirtalact,nampvirtbeact
integer nampoccalact,nampoccbeact,nsumvirtalact,nsumvirtbeact
integer ntampvirtalact,ntampvirtbeact,nvintnewalact,nvintnewbeact
integer nvirtnewalact,nvirtnewbeact,noccnewalact,noccnewbeact
integer nsumoccalact,ntampoccalact,nsumoccbeact,ntampoccbeact
integer earec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer nointnewalact,nointnewbeact,isgnlt,nvirtnewallen
integer econf(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 nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,xyzomp
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe
integer nvintnewal,nvintnewbe,nointnewal,nointnewbe
integer nvirtnewal,nvirtnewbe,noccnewal,noccnewbe
integer nsumvirt,nsumocc,nvirtoldal,noccoldal,nvirtoldbe,noccoldbe
c
integer nvintnew,nointnew,nintnew,nvirtnew,noccnew,nvintold
integer nointold,nintold,nvirtold,noccold
c
real*8 v(*)
logical lll,lcalc,lcalcn,lcalcs,lcalcsn,ladd,let
C Initialize variables
#if defined (OMP)
xyzomp=xyzsiz1
#else
xyzomp=1
#endif
if(isave.eq.2) call ifillzero(iadda,wsmax+1)
nvintnew=nvintnew1
nointnew=nointnew1
nintnew=nvintnew+nointnew
nvirtnew=nvirtnew1
noccnew=noccnew1
C
nvintold=nvintold1
nointold=nointold1
nintold=nvintold+nointold
nvirtold=nvirtold1
noccold=noccold1
C
let=leftroute.and.(namp1.gt.0.or.namp1.lt.-iimed-op).and.
$(itypb.eq.1.or.itypb.eq.4).and.(itypa.eq.1.or.itypa.eq.4)
lcalc=((calc.ne.1.or.(calc.eq.1.and.nintnew.eq.0))
$.and.(.not.let) !(let.and.nvirtnew+noccnew.eq.0))
$.and.multip.ne.2.and.nal.eq.nbe.and.nactva.eq.nactvb.and.
$nactoa.eq.nactob.and.rhf)
c szemet eleje
c lcalc=.false.
c szemet vege
lcalcn=(.not.lcalc)
C
isgnlt=1
if(let) isgnlt=-1
namp=iabs(namp1)
file4=file3
irec3=0
if(ptroute) then
if(eomroute.and.namp1.lt.0) file4=ptfile
if(leftroute) then
if(namp1.gt.0) then
file4=ampfile
else if(namp1.ge.-op) then
file4=lamfile
else
file4=ptfile
namp=namp-iimed-op
endif
endif
else
if((eomroute.or.leftroute).and.namp1.gt.0) file4=cfile
if((eomgrad.and.l3route).and.namp1.lt.-iimed-op) then
namp=iabs(namp1)-iimed-op
irec3=(iroot-1)*trecmax
file4=cfile
endif
endif
nampsym=mult(nnewsym,noldsym)
C
nsumvirt=nvintold-nvintnew
nsumocc=nointold-nointnew
nsum=nsumvirt+nsumocc
ntampvirtact=iabs(nvirtnewact-nvirtoldact)
ntampoccact=iabs(noccnewact-noccoldact)
nampvirtactmax=min(nactv,nsumvirt+ntampvirtact)
nampoccactmax=min(nacto,nsumocc+ntampoccact)
nampvirtactmin=max(namp-mrop,ntampvirtact+nvintoldact-nvintnewact)
nampoccactmin=max(namp-mrop,ntampoccact+nointoldact-nointnewact)
if(lc3) then !c3
nampvirtactmax=nampvirtactmin
nampoccactmax=nampoccactmin
endif
C Loop over spin cases of intermediates
do iw=1,wspc2(0)
nvintoldal=wspc2((iw-1)*12+1)
nointoldal=wspc2((iw-1)*12+2)
nvirtoldal=wspc2((iw-1)*12+3)
noccoldal=wspc2((iw-1)*12+4)
nvintoldalact=wspc2((iw-1)*12+5)
nointoldalact=wspc2((iw-1)*12+6)
nvintoldbeact=wspc2((iw-1)*12+7)
nointoldbeact=wspc2((iw-1)*12+8)
nvirtoldalact=wspc2((iw-1)*12+9)
noccoldalact=wspc2((iw-1)*12+10)
nvirtoldbeact=nvirtoldact-nvirtoldalact
noccoldbeact=noccoldact-noccoldalact
ircold=wspc2((iw-1)*12+11)
ntoldlen=wspc2((iw-1)*12+12)
ntoldleno=ntoldlen
nvintoldbe=nvintold-nvintoldal
nointoldbe=nointold-nointoldal
nvirtoldbe=nvirtold-nvirtoldal
noccoldbe=noccold-noccoldal
if(ladd) then
ia1=1
else
call fwspc(iaddo,nvintoldal,nointoldal,nvirtoldal,noccoldal,
$nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact,
$nvirtoldalact,noccoldalact,ia1,ia2)
endif
if(ia1.gt.0) then
do nampvirtal=0,namp
nampvirtbe=namp-nampvirtal
nampoccal=nampvirtal
nampoccbe=namp-nampoccal
do nampvirtact=nampvirtactmin,nampvirtactmax
nsumvirtact=nampvirtact-ntampvirtact
do nampoccact=nampoccactmin,nampoccactmax
nsumoccact=nampoccact-ntampoccact
do nampvirtalact=max(0,nampvirtact-nactvb),
$min(nactva,nampvirtact,nampvirtal)
nampvirtbeact=nampvirtact-nampvirtalact
do nampoccalact=max(0,nampoccact-nactob),
$min(nactoa,nampoccact,nampoccal)
nampoccbeact=nampoccact-nampoccalact
if((namp1.gt.0.or.(leftroute.and.(iroot.eq.1.or.(zroute.and.
$(.not.osc)).or.(l3route.and.osc.and.namp1.lt.0.and.namp1.ge.-op))
$.and.(.not.l1route).and.namp1.ge.-op)).and.(calc.ne.0
$.or.t1route.or.(calc.eq.0.and.isy1.ne.isy2.and.isym.eq.isy1)))then
iamprec=trec(nampvirtalact,nampvirtbeact,
$nampoccalact,nampoccbeact,nampvirtal,namp)
namplen=nconf(nampvirtalact,nampvirtbeact,nampoccalact,
$nampoccbeact,nampvirtal,namp)
else
iamprec=erec(nampvirtalact,nampvirtbeact,
$nampoccalact,nampoccbeact,nampvirtal,namp)
namplen=econf(nampvirtalact,nampvirtbeact,nampoccalact,
$nampoccbeact,nampvirtal,namp)
endif
if(iamprec.gt.0.and.
$(lcalcn.or.nampvirtal.lt.nampvirtbe.or.(nampvirtal.eq.nampvirtbe
$.and.nampvirtalact.lt.nampvirtbeact).or.(nampvirtal.eq.nampvirtbe
$.and.nampvirtalact.eq.nampvirtbeact.and.nampoccalact.le.
$nampoccbeact))) then
do nsumvirtal=max(0,isign(nvirtoldal+nampvirtal,isgnlt)-minvirtal
$,nsumvirt-min(nvintoldbe,nampvirtbe)),
$min(nvintoldal,nsumvirt,nampvirtal,
$minvirtbe-isign(nvirtoldbe+nampvirtbe,isgnlt)+nsumvirt)
ntampvirtal=nampvirtal-nsumvirtal
nsumvirtbe=nsumvirt-nsumvirtal
ntampvirtbe=nampvirtbe-nsumvirtbe
nvintnewal=nvintoldal-nsumvirtal
nvintnewbe=nvintoldbe-nsumvirtbe
nvirtnewal=nvirtoldal+isign(ntampvirtal,isgnlt)
nvirtnewbe=nvirtoldbe+isign(ntampvirtbe,isgnlt)
do nsumvirtalact=max(0,nsumvirtact-min(nsumvirtbe,nvintoldbeact
$,nampvirtbeact),nvintoldalact-nvintnewal),min(nsumvirtal,
$nvintoldalact,nsumvirtact,nampvirtalact,nactva,
$nvintnewbe-nvintoldbeact+nsumvirtact)
nsumvirtbeact=nsumvirtact-nsumvirtalact
ntampvirtalact=nampvirtalact-nsumvirtalact
ntampvirtbeact=nampvirtbeact-nsumvirtbeact
nvintnewalact=nvintoldalact-nsumvirtalact
nvintnewbeact=nvintoldbeact-nsumvirtbeact
nvirtnewalact=nvirtoldalact+isign(ntampvirtalact,isgnlt)
nvirtnewbeact=nvirtoldbeact+isign(ntampvirtbeact,isgnlt)
if(nvirtnewbe.ge.0.and.nvirtnewbeact.ge.0.and.
$nvintnewalact.le.nvintnewal.and.nvintnewbeact.le.nvintnewbe.and.
$nvirtnewalact.le.nvirtnewal.and.nvirtnewbeact.le.nvirtnewbe.and.
$ntampvirtalact.le.ntampvirtal.and.ntampvirtbeact.le.ntampvirtbe)
$then
do nsumoccal=max(0,isign(noccoldal+nampoccal,isgnlt)-minoccal
$,nsumocc-min(nointoldbe,nampoccbe)),
$min(nointoldal,nsumocc,nampoccal,
$minoccbe-isign(noccoldbe+nampoccbe,isgnlt)+nsumocc)
ntampoccal=nampoccal-nsumoccal
nsumoccbe=nsumocc-nsumoccal
ntampoccbe=nampoccbe-nsumoccbe
nointnewal=nointoldal-nsumoccal
nointnewbe=nointoldbe-nsumoccbe
noccnewal=noccoldal+isign(ntampoccal,isgnlt)
noccnewbe=noccoldbe+isign(ntampoccbe,isgnlt)
do nsumoccalact=max(0,nampoccalact-ntampoccal,
$nsumoccact-min(nsumoccbe,nointoldbeact,nampoccbeact),
$nointoldalact-nointnewal),min(nsumoccal,nointoldalact,nsumoccact,
$nampoccalact,nactoa,ntampoccbe-nampoccbeact+nsumoccact,
$nointnewbe-nointoldbeact+nsumoccact)
nsumoccbeact=nsumoccact-nsumoccalact
ntampoccalact=nampoccalact-nsumoccalact
ntampoccbeact=nampoccbeact-nsumoccbeact
nointnewalact=nointoldalact-nsumoccalact
nointnewbeact=nointoldbeact-nsumoccbeact
noccnewalact=noccoldalact+isign(ntampoccalact,isgnlt)
noccnewbeact=noccoldbeact+isign(ntampoccbeact,isgnlt)
call fwspc(wspc1,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,ircnew,ntnewlen)
C
if(ircnew.gt.0.and.noccnewbe.ge.0.and.noccnewbeact.ge.0
$.and.(lcalcn.or.nampvirtal.ne.nampvirtbe.or.nampvirtalact.ne.
$nampvirtbeact.or.nampoccalact.ne.nampoccbeact.or.
$nvirtnewal.lt.nvirtnewbe.or.(nvirtnewal.eq.nvirtnewbe.and.
$nvirtnewalact.lt.nvirtnewbeact).or.(nvirtnewal.eq.nvirtnewbe.and.
$nvirtnewalact.eq.nvirtnewbeact.and.noccnewalact.le.noccnewbeact))
$) then
lcalcs=lcalc.and.nampvirtal.eq.nampvirtbe.and.nampvirtalact.eq.
$nampvirtbeact.and.nampoccalact.eq.nampoccbeact.and.
$nvirtnewal.eq.nvirtnewbe.and.nvirtnewalact.eq.nvirtnewbeact.and.
$noccnewalact.eq.noccnewbeact
c lcalcs=.false. !szemet
lcalcsn=.not.lcalcs
C Transposition of old intermediate
call transposition(nmax,nnir,nactm,nstr,icmem,isympv,isympo,iwa,
$iwan,v,isig12,ntoldlen,ntoldleno,namp1,file2,irec2,ircold,ita,
$ntnewlen,incnew,noldsym,nnewsym,nvintoldal,nvintoldbe,nointoldal,
$nointoldbe,nvintoldalact,nvintoldbeact,nointoldalact,nointoldbeact
$,nvintnewalact,nvintnewbeact,nointnewalact,nointnewbeact,
$nsumvirtalact,nsumvirtbeact,nsumoccalact,nsumoccbeact,
$nvirtoldalact,nvirtoldbeact,noccoldalact,noccoldbeact,
$nvirtnewalact,nvirtnewbeact,noccnewalact,noccnewbeact,incsum,
$.true.,iwa,incold,nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,
$nvintnewal,nvintnewbe,nointnewal,nointnewbe,nvirtnewal,nvirtnewbe,
$noccnewal,noccnewbe,nintnew,nintold,nvintold,nointold,nvirtoldal,
$nvirtoldbe,noccoldal,noccoldbe,.true.,.not.let,0,i,i,0,ioffs,
$.false.,xyzomp)
C Addresses of summation indices
do ir=1,nir
call tlength(nmax,nstr,nnir,isympv,isympo,nactm,jj,ir,
$nsumvirtalact,nsumvirtbeact,nsumvirtal,nsumvirtbe,
$nsumoccalact,nsumoccbeact,nsumoccal,nsumoccbe,isa,1,0)
enddo
C Read amplitudes
if((namp1.gt.0.or.(leftroute.and.(iroot.eq.1.or.(zroute.and.
$(.not.osc)).or.(l3route.and.osc.and.namp1.lt.0.and.namp1.ge.-op))
$.and.(.not.l1route).and.namp1.ge.-op)).and.(calc.ne.0
$.or.t1route.or.(calc.eq.0.and.isy1.ne.isy2.and.isym.eq.isy1)))then
read(tafile,rec=tarec(nampvirtalact,nampvirtbeact,nampoccalact,
$nampoccbeact,nampvirtal,namp)) ita
else
read(tafile,rec=earec(nampvirtalact,nampvirtbeact,nampoccalact,
$nampoccbeact,nampvirtal,namp)) ita
endif
nmem=ntnewlen+ntoldlen+namplen+1
C Loop over symmetry cases of cluster amplitudes
do ir=1,nir
nsumsym=isympair(nampsym,ir,1)
ntampsym=isympair(nampsym,ir,2)
do irt=1,nir
ntampsymv=isympair(ntampsym,irt,1)
ntampsymo=isympair(ntampsym,irt,2)
do irtv=1,isympv(0,ntampsymv,ntampvirtalact,ntampvirtbeact,
$ntampvirtal,ntampvirtbe,1)
ntampsymva=isympv(irtv,ntampsymv,ntampvirtalact,
$ntampvirtbeact,ntampvirtal,ntampvirtbe,1)
ntampsymvb=isympv(irtv,ntampsymv,ntampvirtalact,
$ntampvirtbeact,ntampvirtal,ntampvirtbe,2)
do irto=1,isympo(0,ntampsymo,ntampoccalact,ntampoccbeact,
$ntampoccal,ntampoccbe,1)
ntampsymoa=isympo(irto,ntampsymo,ntampoccalact,
$ntampoccbeact,ntampoccal,ntampoccbe,1)
ntampsymob=isympo(irto,ntampsymo,ntampoccalact,
$ntampoccbeact,ntampoccal,ntampoccbe,2)
C Loop over symmetry cases of summation indices
call sumsym(nstr,nmax,icmem,isympv,isympo,nnir,ita,nactm,isa,
$nsumsym,nsyma,nsumlen,nsumvirtalact,nsumvirtbeact,nsumvirtallen,
$nsumvirtbelen,ntampsymva,ntampsymvb,nampvirtalact,nampvirtbeact,
$nampvirtlen,nsumoccalact,nsumoccbeact,nsumoccallen,nsumoccbelen,
$ntampsymoa,ntampsymob,nampoccalact,nampoccbeact,nampocclen,
$ntnewlen,ntoldlen,ntampvirtal,ntampvirtalact,ntampvirtbe,
$ntampvirtbeact,ntampoccal,ntampoccalact,ntampoccbe,ntampoccbeact,1
$,nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,
$nampvirtal,nampvirtbe,nampoccal,nampoccbe)
nsumlen=incsum(nsumsym)
if(nsumlen.gt.0) then
C Loop over symmetry cases of new intermediate
nsymw=0
ii=0
do ir1=1,nir
isymi=isympair(nnewsym,ir1,1)
isyma=isympair(nnewsym,ir1,2)
do ira=1,nir
nvirtnewsym=isympair(isyma,ira,1)
noccnewsym=isympair(isyma,ira,2)
do irav=1,isympv(0,nvirtnewsym,nvirtnewalact,nvirtnewbeact,
$nvirtnewal,nvirtnewbe,1)
nvirtnewsyma=isympv(irav,nvirtnewsym,nvirtnewalact,nvirtnewbeact,
$nvirtnewal,nvirtnewbe,1)
nvirtnewallen=nstr(nvirtnewsyma,nvirtnewalact,nvirtnewal,1)
nvirtnewsymb=isympv(irav,nvirtnewsym,nvirtnewalact,nvirtnewbeact,
$nvirtnewal,nvirtnewbe,2)
nvirtnewbelen=nstr(nvirtnewsymb,nvirtnewbeact,nvirtnewbe,2)
nvirtnewlen=nvirtnewallen*nvirtnewbelen
nvirtoldsyma=mult(nvirtnewsyma,ntampsymva)
nvirtoldallen=nstr(nvirtoldsyma,nvirtoldalact,nvirtoldal,1)
nvirtoldsymb=mult(nvirtnewsymb,ntampsymvb)
nvirtoldbelen=nstr(nvirtoldsymb,nvirtoldbeact,nvirtoldbe,2)
nvirtoldsym=mult(nvirtoldsyma,nvirtoldsymb)
nvirtoldlen=nvirtoldallen*nvirtoldbelen
do irao=1,isympo(0,noccnewsym,noccnewalact,noccnewbeact,
$noccnewal,noccnewbe,1)
noccnewsyma=isympo(irao,noccnewsym,noccnewalact,noccnewbeact,
$noccnewal,noccnewbe,1)
noccnewallen=nstr(noccnewsyma,noccnewalact,noccnewal,3)
noccnewsymb=isympo(irao,noccnewsym,noccnewalact,noccnewbeact,
$noccnewal,noccnewbe,2)
noccnewbelen=nstr(noccnewsymb,noccnewbeact,noccnewbe,4)
noccnewlen=noccnewallen*noccnewbelen
noccoldsyma=mult(noccnewsyma,ntampsymoa)
noccoldallen=nstr(noccoldsyma,noccoldalact,noccoldal,3)
noccoldsymb=mult(noccnewsymb,ntampsymob)
noccoldbelen=nstr(noccoldsymb,noccoldbeact,noccoldbe,4)
noccoldsym=mult(noccoldsyma,noccoldsymb)
noccoldlen=noccoldallen*noccoldbelen
C
nintnewlen=incnew(isymi)
nintoldlen=nintnewlen*nsumlen
if(nintnewlen.gt.0.and.nintoldlen.gt.0.and.(lcalcsn.or.
$nvirtnewsyma.lt.nvirtnewsymb.or.(nvirtnewsyma.eq.nvirtnewsymb.and.
$noccnewsyma.le.noccnewsymb))) then
nsymw=nsymw+1
ioldadd=ntnewlen+iwa(isymi,nsumsym,nvirtoldsyma,nvirtoldsymb,
$noccoldsyma,noccoldsymb)
inewadd=iwan(isymi,nvirtnewsyma,nvirtnewsymb,noccnewsyma,
$noccnewsymb)
ii=ii+1
isw(ii)=inewadd
ii=ii+1
isw(ii)=noccnewallen*nvirtnewbelen*nvirtnewallen*nintnewlen
ii=ii+1
isw(ii)=nvirtnewbelen*nvirtnewallen*nintnewlen
ii=ii+1
isw(ii)=nvirtnewallen*nintnewlen
ii=ii+1
isw(ii)=nintnewlen
ii=ii+1
isw(ii)=ioldadd
ii=ii+1
isw(ii)=noccoldallen*nvirtoldbelen*nvirtoldallen*nintoldlen
ii=ii+1
isw(ii)=nvirtoldbelen*nvirtoldallen*nintoldlen
ii=ii+1
isw(ii)=nvirtoldallen*nintoldlen
ii=ii+1
isw(ii)=nintoldlen
ii=ii+1
if(let) then
isw(ii)=icmem(nvirtnewsyma,ntampsymva,nvirtnewalact,
$ntampvirtalact,nvirtnewal,ntampvirtal,1)-1
ii=ii+1
isw(ii)=icmem(nvirtnewsymb,ntampsymvb,nvirtnewbeact,
$ntampvirtbeact,nvirtnewbe,ntampvirtbe,2)-1
ii=ii+1
isw(ii)=icmem(noccnewsyma,ntampsymoa,noccnewalact,ntampoccalact,
$noccnewal,ntampoccal,3)-1
ii=ii+1
isw(ii)=icmem(noccnewsymb,ntampsymob,noccnewbeact,ntampoccbeact,
$noccnewbe,ntampoccbe,4)-1
else
isw(ii)=icmem(nvirtoldsyma,ntampsymva,nvirtoldalact,
$ntampvirtalact,nvirtoldal,ntampvirtal,1)-1
ii=ii+1
isw(ii)=icmem(nvirtoldsymb,ntampsymvb,nvirtoldbeact,
$ntampvirtbeact,nvirtoldbe,ntampvirtbe,2)-1
ii=ii+1
isw(ii)=icmem(noccoldsyma,ntampsymoa,noccoldalact,ntampoccalact,
$noccoldal,ntampoccal,3)-1
ii=ii+1
isw(ii)=icmem(noccoldsymb,ntampsymob,noccoldbeact,ntampoccbeact,
$noccoldbe,ntampoccbe,4)-1
endif
endif
enddo
enddo
enddo
enddo
C Contraction
maxmem=max(maxmem,dble(nmem)+dble(nsumlen))
#if defined (OMP)
ompmem=max(ompmem,dble(nmem)+dble(nsumlen)+
$ dble(xyzsize)*dble(ntnewlen+nsumlen))
#endif
c if(nmem+nsumlen.gt.maxcor) then
c write(6,"(100i3)") nvintnew,nointnew,nvirtnew,noccnew,
c $nvintnewact,nointnewact,nvirtnewact,noccnewact
c write(6,"(10f10.2)") 8.d0*dble(nmem)/dble(twoto20),
c $8.d0*dble(ntnewlen)/dble(twoto20),
c $8.d0*dble(ntoldlen)/dble(twoto20),
c $8.d0*dble(namplen)/dble(twoto20)
c write(6,*) 'M3 ',8.d0*dble(nmem+nsumlen)/dble(twoto20)
c endif
endif
enddo
enddo
enddo
enddo
C Sign of intermediate
if(let) then
if(mod(nsumvirtbe*ntampvirtal+nsumoccbe*ntampoccal+
$nvirtnewal*ntampvirtbe+noccnewal*ntampoccbe+
$nvirtnewal*ntampvirtal+noccnewal*ntampoccal+
$nvirtnewbe*ntampvirtbe+noccnewbe*ntampoccbe,2).ne.0)
$isig12=-isig12
else
if(mod(nsumvirtbe*ntampvirtal+nsumoccbe*ntampoccal+
$nvirtoldbe*ntampvirtal+noccoldbe*ntampoccal,2).ne.0)
$isig12=-isig12
endif
C Save intermediate
call savenew(nmem,ntnewlen,wspca,isave,iadd,nvintnewalact,
$nointnewalact,nvintnewbeact,nointnewbeact,nvirtnewalact,
$noccnewalact,intrec,v,iadda,isig11,isig12,irec1,ircnew,file1,
$lcalc,nampvirtalact,nampvirtbeact,nampoccalact,nampoccbeact,
$nvirtnewbeact,noccnewbeact,wspc1,nstr,nmax,isympv,isympo,ita,
$earec,tarec,nampvirtal,nampvirtbe,nampoccal,nampoccbe,nvintnewal,
$nvintnewbe,nointnewal,nointnewbe,nvirtnewal,nvirtnewbe,noccnewal,
$noccnewbe)
endif
enddo
enddo
endif
enddo
enddo
endif
enddo
enddo
enddo
enddo
enddo
endif !ia1.gt.0
enddo !iw
if(isave.eq.2)
$ call imedret(file1,irec1,wspc1,intfile,intrec,isig11,iadd,v,
$nmax,2,1,1.d0,wspc1,wspca,nactmax,0,0,1,nconf,trec,econf,erec,i,i,
$wsmax,0,0,1,nnewsym)
#if defined (MPI)
if(ilev.eq.1) call imedsync(intfile,intrec,wspca,v)
#endif
C
return
end
C
************************************************************************
subroutine readsp(wspcb,wspca,nmax,wsmax,intn,nim,imed,intrec,
$itypa,itypb,l2map)
************************************************************************
* This subroutine reads spin cases for intermediates and cluster *
* amplitudes *
************************************************************************
#include "MRCCCOMMON"
integer nmax,wmax,wsmax,i,j,k,i5,i6,i7,i8,i9,i10,nim,itr,ii
integer wspcb(0:wsmax,1),wspca(0:wsmax,1),intn(nim),imed(16,1)
integer intrec(*),itypa(*),itypb(*),l2map(*)
C
read(gfile,*) (intn(i),i=1,nim)
i1med=0
do i=1,iimed
if(intn(i).le.3) i1med=max(i1med,i)
enddo
call ifillzero(wspcb,(wspcmax+1)*nimed3)
call ifillzero(wspca,(wspcmax+1)*nimed3)
do i=1,nimed
read(gfile,*) wmax,itr
wspcb(0,itr)=wmax
do j=1,wmax
read(gfile,*) i1,i2,i3,i4,i5,i6,i7,i8,i9,i10
wspcb((j-1)*12+1,itr)=i1
wspcb((j-1)*12+2,itr)=i2
wspcb((j-1)*12+3,itr)=i3
wspcb((j-1)*12+4,itr)=i4
wspcb((j-1)*12+5,itr)=i5
wspcb((j-1)*12+6,itr)=i6
wspcb((j-1)*12+7,itr)=i7
wspcb((j-1)*12+8,itr)=i8
wspcb((j-1)*12+9,itr)=i9
wspcb((j-1)*12+10,itr)=i10
enddo
read(gfile,*) wmax,itr
wspca(0,itr)=wmax
do j=1,wmax
read(gfile,*) i1,i2,i3,i4,i5,i6,i7,i8,i9,i10
wspca((j-1)*12+1,itr)=i1
wspca((j-1)*12+2,itr)=i2
wspca((j-1)*12+3,itr)=i3
wspca((j-1)*12+4,itr)=i4
wspca((j-1)*12+5,itr)=i5
wspca((j-1)*12+6,itr)=i6
wspca((j-1)*12+7,itr)=i7
wspca((j-1)*12+8,itr)=i8
wspca((j-1)*12+9,itr)=i9
wspca((j-1)*12+10,itr)=i10
enddo
enddo
C Definition of itypa/itypb
C 0: IT
C 1: LIT
C 2: LT
C 3: ITR
C 4: LITR
C 5: LTR
C 6: IC
C 7: LC
C Read intermediates
do i=1,nimed
read(gfile,*) (imed(j,i),j=1,16),intrec(i),itypa(i),itypb(i)
enddo
if(calc.eq.0) then
C CI case
C ccimed=eomimed=leftimed,
C densimed=l3imed=d2imed=nimed=nimed2
C
C Numbering of intermediates
C 1,iimed: integral lists
C iimed+1,ccimed: <K|H C|0> or <K|H dC/dx|0> mx elements
C ccimed+1,densimed: (derivative) density matrix elements
C densimed+1,t1imed: T1, <K|dH/dx C|0>
C nimed2+1,nimed3: dipole moment integrals
t1imed=nimed
t2imed=nimed
l1imed=nimed
l2imed=nimed
d1imed=nimed
if(ptroute) t1imed=d2imed+ccimed
do i=iimed+1,nimed
if(itypa(i).eq.0) itypa(i)=6
if(itypa(i).eq.2) itypa(i)=7
enddo
C T1 intermediates
call derimed(wspcb,wspca,wsmax,imed,intrec,itypa,itypb,d2imed,
$t1imed)
else
C CC case
C Numbering of intermediates
C 1,iimed: integral lists
C iimed+1,ccimed: CC intemediates
C ccimed+1,eomimed: right-hand EOMCC intemediates
C eomimed+1,leftimed: left-hand EOMCC intemediates
C leftimed+1,densimed: density matrix intemediates
C densimed+1,l3imed: L3, <0|(1+L)[exp(-T) H exp(T), dT/dx]|K>
C l3imed+1,d2imed: D2, <0|(1+L)[exp(-T) {p+q-} exp(T), dT/dx]|0>
C d2imed+1,t1imed: T1, <K|exp(-T) dH/dx exp(T)|0>
C t1imed+1,t2imed: T2, <K|[exp(-T) H exp(T), dT/dx]|0>
C t2imed+1,l1imed: L1, <0|dL/dx exp(-T) H exp(T)|K>
C l1imed+1,d1imed: D1, <0|dL/dx exp(-T) {p+q-} exp(T)|K>
C d1imed+1,l2imed: L2, <0|(1+L) exp(-T) dH/dx exp(T)|K>
C nimed2+1,nimed3: dipole moment integrals
if(ptroute) then
t1imed=d2imed+ccimed
t2imed=d2imed+eomimed
l1imed=d2imed+leftimed
d1imed=d2imed+densimed
l2imed=nimed2
else
t1imed=nimed
t2imed=nimed
l1imed=nimed
l2imed=nimed
d1imed=nimed
endif
C T1, T2, L1, D1 intermediates
call derimed(wspcb,wspca,wsmax,imed,intrec,itypa,itypb,d2imed,
$d1imed)
C L2 intermediates
if(ptroute) then
ii=d1imed
do i=eomimed+1,leftimed
if(itypa(i).eq.1) then
ii=ii+1
l2map(i-eomimed)=ii
do j=1,16
imed(j,ii)=imed(j,i)
enddo
intrec(ii)=intrec(i)
itypa(ii)=itypa(i)
itypb(ii)=itypb(i)
do j=0,wsmax
wspcb(j,ii)=wspcb(j,i)
wspca(j,ii)=wspca(j,i)
enddo
endif
enddo
endif
endif
C Dipole moment integrals
if(.not.hess) then
do i=1,3
call derimed(wspcb,wspca,wsmax,imed,intrec,itypa,itypb,
$nimed+(i-1)*iimed,nimed+i*iimed)
enddo
endif
C
return
end
C
************************************************************************
subroutine derimed(wspcb,wspca,wsmax,imed,intrec,itypa,itypb,ifst,
$iend)
************************************************************************
* Intermediates in derivative calculation *
************************************************************************
#include "MRCCCOMMON"
integer wmax,wsmax,wspcb(0:wsmax,1),wspca(0:wsmax,1),imed(16,1)
integer intrec(*),itypa(*),itypb(*),i,j,ifst,iend
C
do i=ifst+1,iend
do j=1,16
imed(j,i)=imed(j,i-ifst)
enddo
intrec(i)=intrec(i-ifst)
itypa(i)=itypa(i-ifst)
itypb(i)=itypb(i-ifst)
do j=0,wsmax
wspcb(j,i)=wspcb(j,i-ifst)
wspca(j,i)=wspca(j,i-ifst)
enddo
enddo
C
return
end
C
************************************************************************
subroutine imedret(file1,irec1,tre1,file2,irec2,isg,iadd,v,nmax,n,
$ic,c0,wspc1,wspc2,nactm,nact1,nact2,ical,nconf,trec,econf,erec,
$ivact,ioact,wsmax,ilev1,isave,ityp,isyma)
************************************************************************
* This subroutine retrieves intermediates *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,n,file1,file2,irec1,irec2,isg,i,ic,file3,ical,j
integer tre1(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax),wsmax
integer wspc1(0:*),wspc2(0:*),ire1,ire2,ile1,iw,iadd(0:*),ia1,ia2
integer nact1,nact2,i5,i6,i7,i8,i9,i10,ivact,ioact,ilev1,ilev
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),isave
integer econf(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),ityp
integer isyma
real*8 v(10),c0
C
if(densroute.and.ilev1.eq.1.and.isyma.ne.denssym) return
if((leftroute.and.(iroot.gt.1.or.nit.gt.1.or.l1route).and.
$isave.gt.0.and.ityp.eq.0.and.ilev1.lt.4.and.(.not.l2route).and.
$(.not.l3route)).or.(l2route.and.ityp.eq.2.and.isave.gt.0)) then
ilev=ilev1
call skip1(ilev,ftfile)
endif
if(ilev1.eq.1.and.ityp.eq.0.and.l1route) return
file3=file2
if(ical.lt.0) then
file3=lamfile
i=iabs(ical)
if(iroot.gt.1.or.l1route.or.d1route) then
call tspcase(0,0,i,i,0,0,ivact,ioact,wspc2,wsmax,
$econf(0,0,0,0,0,i),erec(0,0,0,0,0,i),nmax,nactm)
else
call tspcase(0,0,i,i,0,0,ivact,ioact,wspc2,wsmax,
$nconf(0,0,0,0,0,i),trec(0,0,0,0,0,i),nmax,nactm)
endif
endif
if(ic.eq.0.and.c0.eq.0.d0) return
if(n.eq.0) then
do iw=1,wspc2(0)
i3=wspc2((iw-1)*12+3)
i4=wspc2((iw-1)*12+4)
i9=wspc2((iw-1)*12+9)
i10=wspc2((iw-1)*12+10)
ire2=wspc2((iw-1)*12+11)
ile1=wspc2((iw-1)*12+12)
ire1=tre1(i9,nact1-i9,i10,nact2-i10,i3)
if(ile1.gt.0.and.ire1.gt.0) then
maxmem=max(maxmem,dble(ile1))
c if(ile1.gt.maxcor) write(6,*) 'M4 ',8.d0*dble(ile1)/
c $dble(twoto20)
call fwspc(iadd,0,0,i3,i4,0,0,0,0,i9,i10,ia1,ia2)
if(ia1.gt.0) then
else
call wiadd(iadd,0,0,i3,i4,0,0,0,0,i9,i10,ile1)
endif
endif
enddo
else
do iw=1,wspc1(0)
i1=wspc1((iw-1)*12+1)
i2=wspc1((iw-1)*12+2)
i3=wspc1((iw-1)*12+3)
i4=wspc1((iw-1)*12+4)
i5=wspc1((iw-1)*12+5)
i6=wspc1((iw-1)*12+6)
i7=wspc1((iw-1)*12+7)
i8=wspc1((iw-1)*12+8)
i9=wspc1((iw-1)*12+9)
i10=wspc1((iw-1)*12+10)
ire1=wspc1((iw-1)*12+11)
call fwspc(wspc2,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,ire2,ile1)
if(ile1.gt.0) then
maxmem=max(maxmem,dble(ile1))
c if(ile1.gt.maxcor) write(6,*) 'M4 ',8.d0*dble(ile1)/dble(twoto20)
call fwspc(iadd,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,ia1,ia2)
if(ia1.gt.0) then
else
call wiadd(iadd,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,ile1)
endif
endif
enddo
endif
C
return
end
C
************************************************************************
subroutine symini(mosym,mosymvo,nb,noccup,nactd)
************************************************************************
* This subroutine initializes symmetry variables *
************************************************************************
#include "MRCCCOMMON"
integer i,j,nb,mosym(*),mosymvo(nb,4),noccup(*),ii
integer nactd(*),jj,kk
C Initialize multiplication table
do i=1,8
do j=1,8
mult(i,j)=multpg(i,j)
enddo
enddo
C Read symmetry labels
read(inp,*) (mosym(i),i=1,nbasis/2)
if(ptroute) then
read(inp,*) ptsym(1),nopt(1) !szemet
read(inp,*) (i,j=1,nir) !szemet
endif
C Read multiplication table for double groups
if(rel.eq.2) then
read(inp,*) i
read(inp,*) nir
do i=1,nir
read(inp,*) (mult(i,j),j=1,nir)
enddo
endif
denssym=1
nir=1
if(isym.ne.0) then
do i=1,nbasis/2
nir=max(nir,mosym(i))
enddo
if(nir.gt.2) then
if(nir.le.4) then
nir=4
else
nir=8
endif
endif
do i=1,nbasis/2
mosym(nbasis/2+i)=mosym(i)
enddo
refsym=1
do i=1,nbasis/2
if(iabs(noccup(i)).eq.1) refsym=mult(refsym,mosym(i))
enddo
if(refsym.ne.isym.and.eomcc)
$write(iout,*) 'Spatial symmetry of ground state: ',refsym
C Complex conjugate of irreps
do ir=1,nir
do i=1,ir
c write(6,*) ir,i,mult(ir,i)
if(mult(ir,i).eq.1) then
iconj(ir)=i
iconj(i)=ir
endif
enddo
enddo
if(eomcc) then
isy1=1
else
isy1=mult(iconj(refsym),isy1)
endif
isy2=mult(iconj(refsym),isy2)
isym=mult(iconj(refsym),isym)
if(ptroute) then
denssym=ptsym(1)
ptsym(1)=mult(isym,ptsym(1))
write(iout,*) 'Symmetry of perturbation: ',ptsym(1)
else
ptsym(1)=isym
if(osc) denssym=mult(iconj(isy1),isy2)
endif
else
iconj(1)=1
isym=1
isy1=1
isy2=1
ptsym(1)=1
refsym=1
do i=1,nbasis
mosym(i)=1
enddo
endif
C Orbital symmetries
C Virtual alpha
ii=0
do i=1,nbasis/2
if(noccup(i).le.0.and.nactd(i).eq.1) then
ii=ii+1
mosymvo(ii,1)=mosym(i)
endif
enddo
do i=1,nbasis/2
if(noccup(i).le.0.and.nactd(i).eq.0) then
ii=ii+1
mosymvo(ii,1)=mosym(i)
endif
enddo
C Virtual beta
ii=0
do i=1,nbasis/2
if((noccup(i).eq.1.or.noccup(i).eq.0).and.nactd(i).eq.1) then
ii=ii+1
mosymvo(ii,2)=mosym(i)
endif
enddo
do i=1,nbasis/2
if((noccup(i).eq.1.or.noccup(i).eq.0).and.nactd(i).eq.0) then
ii=ii+1
mosymvo(ii,2)=mosym(i)
endif
enddo
C Occupied alpha
ii=0
do i=1,nbasis/2
if(noccup(i).ge.1.and.nactd(i).eq.0) then
ii=ii+1
mosymvo(ii,3)=mosym(i)
endif
enddo
do i=1,nbasis/2
if(noccup(i).ge.1.and.nactd(i).eq.1) then
ii=ii+1
mosymvo(ii,3)=mosym(i)
endif
enddo
C Occupied beta
ii=0
do i=1,nbasis/2
if((noccup(i).eq.-1.or.noccup(i).eq.2).and.nactd(i).eq.0) then
ii=ii+1
mosymvo(ii,4)=mosym(i)
endif
enddo
do i=1,nbasis/2
if((noccup(i).eq.-1.or.noccup(i).eq.2).and.nactd(i).eq.1) then
ii=ii+1
mosymvo(ii,4)=mosym(i)
endif
enddo
C Generate symmetry pairs
do ir=1,nir
ii=0
jj=0
kk=0
do i=1,nir
do j=1,nir
if(mult(i,j).eq.ir) then
ii=ii+1
isympair(ir,ii,1)=i
isympair(ir,ii,2)=j
c write(6,"(a4,10i2)") 'isym',ir,i,j
endif
if(mult(iconj(i),j).eq.ir) then
jj=jj+1
csympair(ir,jj,1)=i
csympair(ir,jj,2)=j
c write(6,"(a4,10i2)") 'csym',ir,i,j
endif
if(mult(i,iconj(j)).eq.ir) then
kk=kk+1
dsympair(ir,kk,1)=i
dsympair(ir,kk,2)=j
endif
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine symmetrize(nconf,trec,nstr,nmax,nnir,isympv,
$isympo,ita,tarec,v,ifile,nactm,ssym)
************************************************************************
* This subroutine transposes CI coefficients *
************************************************************************
#include "MRCCCOMMON"
integer nnir,isymv,isymo,isymva,isymvb,isymoa,isymob,irv,iro,ii,n
integer nmax,nactm,i,j,k,nex,l,iactob,ssym,nvirtnewallen
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer nstr(nnir,0:nactm,0:nmax,4),ifile,n2,kk,ll,irav,irao
integer isympv(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2),jj,n1
integer isympo(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2),m,ira
integer nvstral,nvstrbe,nostral,nostrbe,navirtal,navirtbe
integer ivstral,ivstrbe,iostral,iostrbe,naoccal,naoccbe
integer nvirtnewsym,noccnewsym,nvirtnewsyma,nvirtnewsymb
integer nvirtnewbelen,noccnewsyma,noccnewsymb,noccnewbelen
integer ita(nnir,nnir,nnir,nnir),iactv,iacto,iactva,iactoa,iactvb
integer tarec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe
real*8 v(*),sum
C
if(calc.ne.0) return !szemet
if(multip.eq.2.or.nal.ne.nbe.or.nactva.ne.nactvb.or.
$nactoa.ne.nactob.or.(.not.rhf)) return
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
nampvirtal=i1
nampvirtbe=nex-nampvirtal
nampoccal=i1
nampoccbe=nex-nampoccal
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
ii=nconf(iactva,iactvb,iactoa,iactob,i1,nex)
if(ii.gt.0) then
maxmem=max(maxmem,2.d0*dble(ii))
c if(2*ii.gt.maxcor) write(6,*) 'M5 ',8.d0*dble(2*ii)/dble(twoto20)
endif
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine fwspc(wspc,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,ire,ile)
************************************************************************
* Give back the record number and the length of an intermediate *
************************************************************************
#include "MRCCCOMMON"
integer wspc(0:1),n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,ire,ile,i
C
i=0
9879 continue
if(wspc(i+11).eq.0) then
ire=0
ile=0
return
endif
i=i+1
if(wspc(i).ne.n1) then
i=i+11
goto 9879
endif
i=i+1
if(wspc(i).ne.n2) then
i=i+10
goto 9879
endif
i=i+1
if(wspc(i).ne.n3) then
i=i+9
goto 9879
endif
i=i+1
if(wspc(i).ne.n4) then
i=i+8
goto 9879
endif
i=i+1
if(wspc(i).ne.n5) then
i=i+7
goto 9879
endif
i=i+1
if(wspc(i).ne.n6) then
i=i+6
goto 9879
endif
i=i+1
if(wspc(i).ne.n7) then
i=i+5
goto 9879
endif
i=i+1
if(wspc(i).ne.n8) then
i=i+4
goto 9879
endif
i=i+1
if(wspc(i).ne.n9) then
i=i+3
goto 9879
endif
i=i+1
if(wspc(i).ne.n10) then
i=i+2
goto 9879
endif
ire=wspc(i+1)
ile=wspc(i+2)
C
return
end
C
************************************************************************
subroutine wiadd(iadd,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,ilen)
************************************************************************
* Write in iadd array *
************************************************************************
implicit none
integer iadd(0:12),n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,i,ilen
C
if(ilen.eq.0) return
i=iadd(0)
iadd(0)=i+1
i=i*12
iadd(i+1)=n1
iadd(i+2)=n2
iadd(i+3)=n3
iadd(i+4)=n4
iadd(i+5)=n5
iadd(i+6)=n6
iadd(i+7)=n7
iadd(i+8)=n8
iadd(i+9)=n9
iadd(i+10)=n10
iadd(i+11)=ilen
C
return
end
C
************************************************************************
subroutine david(bmat,bvec,invbmat,cvec,ndmax,v,nconf,trec,nmax,
$nactm,test,norm,file)
************************************************************************
* Davidson extrapolation *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,nact,nactpl,iact,ii,i,j,ntrecact,file
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer ndmax,iactv,iacto,iactva,iactoa,iactvb,iactob,n,nn,m,k
real*8 bmat(ndmax,ndmax),bvec(ndmax,ndmax),invbmat(ndmax,ndmax)
real*8 cvec(ndmax),sum,sum1,v(*),norm,test,norm2,denom
logical lll
C
nact=min(ndiis,nit)
nactpl=mod(nit,ndiis)
if(nactpl.eq.0) nactpl=ndiis
ntrecact=(nactpl-1)*trecmax
C Update the scalar product matrix (H matrix)
do i=0,op
do iactv=max(0,i-mrop),min(nactv,i)
do iacto=max(0,i-mrop),min(nacto,i)
do i1=0,i
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
ii=nconf(iactva,iactvb,iactoa,iactob,i1,i)
if(ii.gt.0) then
maxmem=max(maxmem,2.d0*dble(ii))
c if(2*ii.gt.maxcor) write(6,*) 'M6 ',8.d0*dble(2*ii)/dble(twoto20)
endif
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine nconfig(nvstral,nvstrbe,nostral,nostrbe,nconf,tsp,nmax,
$nnir,isympv,isympo,ita,tarec,nactm,ii,iisym,trec,trec2,tcase,
$ntcase,lll)
************************************************************************
* This subroutine calculates the number of excitations *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,ii,nact,nactm,iisym,nma
integer nvstral(nnir,0:nactm,0:nmax),nvstrbe(nnir,0:nactm,0:nmax)
integer nostral(nnir,0:nactm,0:nmax),nostrbe(nnir,0:nactm,0:nmax)
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe
integer tsp(0:2*(nmax+1)**2,0:nmax),nex,i,j,k,m,isum,iactoa
integer isymv,isymo,isymva,isymvb,isymoa,isymob,irv,iro,iactob
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 ita(nnir,nnir,nnir,nnir),iactv,iacto,iactva,iactvb
integer tarec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer tcase(0:2*nactm,0:2*nactm,0:nmax),ntcase(0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec2((2*nactm+1)**2,0:nmax)
logical lll
C
nma=op
if(pert.ne.0) nma=op-1
ntot=0
call ifillzero(nconf,(nactm+1)**4*(oo+1)**2)
do i=0,nma
isum=0
do iactv=max(0,i-mrop),min(nactv,i)
do iacto=max(0,i-mrop),min(nacto,i)
if(tcase(iactv,iacto,i).gt.0) then
do iactva=max(0,iactv-nactvb),min(nactva,iactv)
iactvb=iactv-iactva
do iactoa=max(0,iacto-nactob),min(nactoa,iacto)
iactob=iacto-iactoa
do j=1,tsp(0,i)
nampvirtal=tsp(2*j-1,i)
nampvirtbe=i-nampvirtal
nampoccal=tsp(2*j,i)
nampoccbe=i-nampoccal
if(lll) then
ii=ii+1
tarec(iactva,iactvb,iactoa,iactob,nampvirtal,i)=ii
endif
m=0
do ir=1,nir
isymv=csympair(iisym,ir,1)
isymo=csympair(iisym,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)
k=nvstral(isymva,iactva,nampvirtal)*
$ nvstrbe(isymvb,iactvb,nampvirtbe)
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)
if(lll) ita(isymva,isymvb,isymoa,isymob)=m+1
m=m+k*nostral(isymoa,iactoa,nampoccal)*
$ nostrbe(isymob,iactob,nampoccbe)
enddo
enddo
enddo
c call slength(nmax,nvstral,nnir,isympv,isympo,nactm,m,iisym,
c $iactva,iactvb,nampvirtal,nampvirtbe,iactoa,iactob,nampoccal,
c $nampoccbe)
if(lll) write(tafile,rec=ii) ita
nconf(iactva,iactvb,iactoa,iactob,nampvirtal,i)=m
isum=isum+m
enddo
enddo
enddo
endif
enddo
enddo
if(lll) write(iout,*) 'Number of',i,'-fold excitations: ',isum
ntot=ntot+isum
enddo
if(lll) write(iout,*) 'Total number of configurations: ',ntot
C Records in amplitude file
call ifillzero(trec,(oo+1)**2*(nactm+1)**4)
j=1
do nex=0,nma
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
isum=nconf(iactva,iactvb,iactoa,iactob,i1,nex)
if(isum.gt.0) then
trec(iactva,iactvb,iactoa,iactob,i1,nex)=j
k=mod(isum,ibufln)
if(k.ne.0) k=ibufln-k
j=j+(isum+k)/ibufln
endif
enddo
enddo
enddo
enddo
enddo
enddo
trecmax=max(trecmax,j-1) !???
do i=0,nma
do k=1,ntcase(i)
trec2(k,i)=trecmax
enddo
do iactv=max(0,i-mrop),min(nactv,i)
do iacto=max(0,i-mrop),min(nacto,i)
do i1=0,i
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
j=trec(iactva,iactvb,iactoa,iactob,i1,i)
isum=nconf(iactva,iactvb,iactoa,iactob,i1,i)
if(isum.gt.0) trec2(tcase(iactv,iacto,i),i)=
$min(trec2(tcase(iactv,iacto,i),i),j-1)
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine ilength(nmax,nstr,imed,nnir,isympv,isympo,nactm,wspc,m,
$ssym)
************************************************************************
* This subroutine calculates the length of intermediates *
************************************************************************
#include "MRCCCOMMON"
integer nnir,nactm,m,nmax,n1,n2,wspc(*),imed(*),isymi,isyma
integer nstr(nnir,0:nactm,0:nmax,4),ssym,nactoccal,nactvirtbe
integer nvintoldal,nvintoldbe,nointoldal,nointoldbe,nactoccbe
integer nactvintal,nactointal,nactvintbe,nactointbe,nactvirtal
integer nvirtoldal,noccoldal,nvirtoldbe,noccoldbe,isympv,isympo
C
nvintoldal=wspc(1)
nointoldal=wspc(2)
nvirtoldal=wspc(3)
noccoldal=wspc(4)
nvintoldbe=imed(1)-nvintoldal
nointoldbe=imed(2)-nointoldal
nvirtoldbe=imed(3)-nvirtoldal
noccoldbe=imed(4)-noccoldal
nactvintal=wspc(5)
nactointal=wspc(6)
nactvintbe=wspc(7)
nactointbe=wspc(8)
nactvirtal=wspc(9)
nactoccal=wspc(10)
nactvirtbe=imed(11)-nactvirtal
nactoccbe=imed(12)-nactoccal
m=0
do ir=1,nir
isymi=isympair(ssym,ir,1)
isyma=isympair(ssym,ir,2)
call slength(nmax,nstr,nnir,isympv,isympo,nactm,n1,isymi,
$nactvintal,nactvintbe,nvintoldal,nvintoldbe,nactointal,nactointbe,
$nointoldal,nointoldbe,dsympair)
call slength(nmax,nstr,nnir,isympv,isympo,nactm,n2,isyma,
$nactvirtal,nactvirtbe,nvirtoldal,nvirtoldbe,nactoccal,nactoccbe,
$noccoldal,noccoldbe,csympair)
m=m+n1*n2
enddo
C
return
end
C
************************************************************************
subroutine slength(nmax,nstr,nnir,isympv,isympo,nactm,isum,ssym,
$nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8,isp)
************************************************************************
* This subroutine calculates the length of indices running over strings*
************************************************************************
#include "MRCCCOMMON"
integer nnir,nactm,isum,nmax,nn,nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8
integer nstr(nnir,0:nactm,0:nmax,4),irao,isymaoa,isymaob,ssym
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 ira,isymav,isymao,irav,isymava,isymavb
integer isp(nirmax,nirmax,2)
C
isum=0
do ira=1,nir
isymav=isp(ssym,ira,1)
isymao=isp(ssym,ira,2)
do irav=1,isympv(0,isymav,nn1,nn2,nn3,nn4,1)
isymava=isympv(irav,isymav,nn1,nn2,nn3,nn4,1)
isymavb=isympv(irav,isymav,nn1,nn2,nn3,nn4,2)
nn=nstr(isymava,nn1,nn3,1)*nstr(isymavb,nn2,nn4,2)
do irao=1,isympo(0,isymao,nn5,nn6,nn7,nn8,1)
isymaoa=isympo(irao,isymao,nn5,nn6,nn7,nn8,1)
isymaob=isympo(irao,isymao,nn5,nn6,nn7,nn8,2)
isum=isum+nstr(isymaoa,nn5,nn7,3)*nstr(isymaob,nn6,nn8,4)*nn
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine conin1(trec,nstr,nmax,icmem,nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvintold1,nointold1,nvirtold1,noccold1,namp1,
$isig11,file1,irec1,file2,irec2,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,nconf,ita,tarec,iwa,isa,isw,iwan,file3,wspc1,
$wspc2,wspca,wsmax,nvintnewact,nointnewact,nvirtnewact,noccnewact,
$nvintoldact,nointoldact,nvirtoldact,noccoldact,nactm,earec,
$erec,econf,ilev,imed,wspci,intrec1,itypa,itypb,inewsym,ioldsym,
$ioffs)
************************************************************************
* This subroutine initializes variables for contractions *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,file3,wsmax,iadd(0:*),isa(*),isw(*),i
integer nstr(nnir,0:nactm,0:nmax,4),wspc3(0:wsmax),ilev,ii,j,file4
integer wspc1(0:12),wspc2(0:12),wspca(0:12),iaddo(0:*)
integer iaddn(0:wsmax)
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ioffs
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),file5
integer nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1,nointold1
integer nvirtold1,noccold1,namp1,isig11,isig12,iadds(0:wsmax)
integer file1,file2,irec1,irec2,ircold,ircnew,isave,intrec,icmem
integer isympv,isympo,itypa,itypb,inewsym,ioldsym,nointnewbeact
integer tarec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer ita(nnir,nnir,nnir,nnir),iwan(nnir,nnir,nnir,nnir,nnir)
integer earec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer nointnewalact,iwa(nnir,nnir,nnir,nnir,nnir,nnir)
integer econf(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 nvintnewact,nointnewact,nvirtnewact,noccnewact,imed(16,1)
integer nvintoldact,nointoldact,nvirtoldact,noccoldact,intrec1(*)
integer nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact
integer nvirtoldalact,noccoldalact,nvirtoldbeact,noccoldbeact
integer wspci(0:wsmax,13),ire11,irec3,file6
real*8 v(*)
logical ladd,lll
C
ire11=irec1
if(leftroute.and.iroot.gt.1.and.itypa.eq.0.and.namp1.gt.0.and.
$ilev.eq.1.and.(((.not.zroute).and.(.not.d1route).and.(.not.osc))
$.or.(osc.and.(.not.l3route).and.(.not.d2route)))) return
C
if(leftroute.and.ilev.eq.1.and.nvintnew1+nointnew1.gt.0.and.
$(.not.densroute)) then
lll=.true.
file5=scrfile4
do i=0,wsmax
iadds(i)=iadd(i)
enddo
call ifillzero(iadd,wsmax+1)
ire11=0
else
lll=.false.
file5=file1
endif
if(leftroute.and.namp1.lt.0.and.namp1.ge.-iimed-op) then
if(namp1.ge.-op) then
C Contract an L vertex with conventional intermediates at the end
if(ilev.eq.1) then
call rospc(iaddo,iaddn,nvirtoldact,noccoldact,wsmax)
call rospc(wspc2,wspc3,nvirtoldact,noccoldact,wsmax)
call conini(trec,nstr,nmax,icmem,nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvirtold1,noccold1,nvintold1,nointold1,namp1,
$isig11,file5,ire11,file2,irec2,iadd,iaddn,ladd,v,isave,intrec,
$isympv,isympo,nnir,nconf,ita,tarec,iwa,isa,isw,iwan,file3,wspc1,
$wspc3,wspca,wsmax,nvintnewact,nointnewact,nvirtnewact,noccnewact,
$nvirtoldact,noccoldact,nvintoldact,nointoldact,nactm,earec,
$erec,econf,itypa,itypb,inewsym,ioldsym,ilev,ioffs)
else
C Contract an L vertex with conventional intermediates not at the end
call conlam(trec,nstr,nmax,icmem,nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvintold1,nointold1,nvirtold1,noccold1,namp1,
$isig11,file5,ire11,file2,irec2,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,nconf,ita,tarec,iwa,isa,isw,iwan,file3,wspc1,
$wspc2,wspca,wsmax,nvintnewact,nointnewact,nvirtnewact,noccnewact,
$nvintoldact,nointoldact,nvirtoldact,noccoldact,nactm,earec,
$erec,econf,inewsym,ioldsym)
endif
else
C Contract an integral list with a L-T intermediate
i=iabs(namp1+op)
if(l2route) i=i+ptsh
call conint(nstr,nmax,icmem,nvintnew1,nointnew1,nvirtnew1,
$noccnew1,nvintold1,nointold1,nvirtold1,noccold1,isig11,
$file5,ire11,file2,intrec1(i),iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,ita,iwa,isa,isw,iwan,wspc1,wspc2,wspca,wsmax,
$nvintnewact,nointnewact,nvirtnewact,noccnewact,nvintoldact,
$nointoldact,nvirtoldact,noccoldact,nactm,imed(1,i),imed(2,i),
$imed(3,i),imed(4,i),imed(9,i),imed(10,i),imed(11,i),imed(12,i),
$wspci(0,i),.false.,nconf,trec,tarec,iwan,lll.or.ilev.ne.1,
$inewsym,ioldsym,ampfile,0,ioffs)
endif
else
if(leftroute.and.namp1.gt.0.and.ilev.eq.1.and.itypa.eq.0)
$ire11=1
if(nvintnew1.gt.nvintold1.or.nointnew1.gt.nointold1.or.
$densroute) then
C Contract a T vertex with a L-T intermediate
file4=file2
if(ilev.eq.4) then
file4=lamfile
if((iroot.eq.1.or.(eomgrad.and.(d1route.or.zroute).and.
$(.not.osc)).or.(osc.and.(l3route.or.d2route))).and.(.not.l1route)
$.and.((.not.d1route).or.eomgrad)) then
call tspcase(0,0,nvirtold1,noccold1,0,0,nvirtoldact,noccoldact,
$wspc2,wsmax,nconf(0,0,0,0,0,nvirtold1),
$trec(0,0,0,0,0,nvirtold1),nmax,nactm)
else
call tspcase(0,0,nvirtold1,noccold1,0,0,nvirtoldact,noccoldact,
$wspc2,wsmax,econf(0,0,0,0,0,nvirtold1),
$erec(0,0,0,0,0,nvirtold1),nmax,nactm)
endif
endif
if((l3route.or.d2route).and.namp1.lt.-iimed-op) then
namp=iabs(namp1)-iimed-op
if(eomgrad) then
irec3=(iroot-1)*trecmax
file6=cfile
else
irec3=0
file6=ptfile
endif
call conint(nstr,nmax,icmem,nvintnew1,nointnew1,nvirtnew1,
$noccnew1,nvintold1,nointold1,nvirtold1,noccold1,isig11,
$file5,ire11,file4,intrec1,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,ita,iwa,isa,isw,iwan,wspc1,wspc2,wspca,wsmax,
$nvintnewact,nointnewact,nvirtnewact,noccnewact,nvintoldact,
$nointoldact,nvirtoldact,noccoldact,nactm,nvintold1,nointold1,
$nvirtold1,noccold1,nvintoldact,nointoldact,nvirtoldact,noccoldact,
$wspc3,.true.,econf(0,0,0,0,0,namp),erec(0,0,0,0,0,namp),
$earec(0,0,0,0,0,namp),iwan,.true.,inewsym,ioldsym,file6,irec3,
$ioffs)
else
call conint(nstr,nmax,icmem,nvintnew1,nointnew1,nvirtnew1,
$noccnew1,nvintold1,nointold1,nvirtold1,noccold1,isig11,
$file5,ire11,file4,intrec1,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,ita,iwa,isa,isw,iwan,wspc1,wspc2,wspca,wsmax,
$nvintnewact,nointnewact,nvirtnewact,noccnewact,nvintoldact,
$nointoldact,nvirtoldact,noccoldact,nactm,nvintold1,nointold1,
$nvirtold1,noccold1,nvintoldact,nointoldact,nvirtoldact,noccoldact,
$wspc3,.true.,nconf(0,0,0,0,0,namp1),trec(0,0,0,0,0,namp1),
$tarec(0,0,0,0,0,namp1),iwan,.true.,inewsym,ioldsym,ampfile,0,
$ioffs)
endif
else
C Contract a T vertex with a conventional or a L-T-I intermediate
call conini(trec,nstr,nmax,icmem,nvintnew1,nointnew1,
$nvirtnew1,noccnew1,nvintold1,nointold1,nvirtold1,noccold1,namp1,
$isig11,file5,ire11,file2,irec2,iadd,iaddo,ladd,v,isave,intrec,
$isympv,isympo,nnir,nconf,ita,tarec,iwa,isa,isw,iwan,file3,wspc1,
$wspc2,wspca,wsmax,nvintnewact,nointnewact,nvirtnewact,noccnewact,
$nvintoldact,nointoldact,nvirtoldact,noccoldact,nactm,earec,
$erec,econf,itypa,itypb,inewsym,ioldsym,ilev,ioffs)
endif
endif
if(lll) then
call anti1(nconf,econf,trec,erec,nstr,nmax,
$v,nnir,isympv,isympo,ita,tarec,earec,nactm,icmem,wspc1,
$nvintnew1,nointnew1,nvirtnew1,noccnew1,nvirtnewact,noccnewact,
$iadds,iadd)
do i=0,wsmax
iadd(i)=iadds(i)
enddo
endif
C
return
end
C
************************************************************************
subroutine rospc(iadd,iaddn,nactv,nacto,wsmax)
************************************************************************
* Reorder iadd arrays *
************************************************************************
implicit none
integer wsmax,iadd(0:12),iaddn(0:12),nactv,nacto,i,j
C
call ifillzero(iaddn,wsmax+1)
iaddn(0)=iadd(0)
do j=0,iadd(0)-1
i=12*j
iaddn(i+1)=iadd(i+3)
iaddn(i+2)=iadd(i+4)
iaddn(i+3)=iadd(i+1)
iaddn(i+4)=iadd(i+2)
iaddn(i+5)=iadd(i+9)
iaddn(i+6)=iadd(i+10)
iaddn(i+7)=nactv-iadd(i+9)
iaddn(i+8)=nacto-iadd(i+10)
iaddn(i+9)=iadd(i+5)
iaddn(i+10)=iadd(i+6)
iaddn(i+11)=iadd(i+11)
iaddn(i+12)=iadd(i+12)
enddo
C
return
end
C
************************************************************************
subroutine wszero(wsmax,wspca,ifi,ila)
************************************************************************
* Zeros redundant spin cases *
************************************************************************
#include "MRCCCOMMON"
integer n,i,j,k,ii,wsmax,wspca(0:wsmax,13),iw,ifi,ila
C
do i=ifi+1,ila
ii=0
do iw=1,wspca(0,i)
if(wspca((iw-ii-1)*12+11,i)*wspca((iw-ii-1)*12+12,i).eq.0)
$then
do k=iw-ii,wspca(0,i)-ii
do j=1,12
wspca((k-1)*12+j,i)=wspca(k*12+j,i)
enddo
enddo
ii=ii+1
endif
enddo
wspca(0,i)=wspca(0,i)-ii
call ifillzero(wspca(wspca(0,i)*12+1,i),ii*12)
enddo
C
return
end
C
************************************************************************
subroutine wswrite(is1,is2,is3,wsa,wsb,wsm,ire)
************************************************************************
* Writes in wsfile *
************************************************************************
#include "MRCCCOMMON"
integer is1,is2,is3,wsm,wsa(wsm),wsb(wsm),ire,i,iscr(2*wsm+3)
C
iscr(1)=is1
iscr(2)=is2
iscr(3)=is3
do i=1,wsm
iscr(3+i)=wsa(i)
enddo
do i=1,wsm
iscr(3+wsm+i)=wsb(i)
enddo
write(wsfile,rec=ire) iscr
C
return
end
C
************************************************************************
subroutine wsch(is1,is2,is3,wsa,wsb,wsm,ire)
************************************************************************
* Cahnges the content of wsfile *
************************************************************************
#include "MRCCCOMMON"
integer is1,is2,is3,wsm,wsa(wsm),wsb(wsm),ire,i,j,iscr(2*wsm+3)
C
read(wsfile,rec=ire) iscr
j=is1
is1=iscr(1)
iscr(1)=j
j=is2
is2=iscr(2)
iscr(2)=j
j=is3
is3=iscr(3)
iscr(3)=j
do i=1,wsm
j=wsa(i)
wsa(i)=iscr(3+i)
iscr(3+i)=j
enddo
do i=1,wsm
j=wsb(i)
wsb(i)=iscr(3+wsm+i)
iscr(3+wsm+i)=j
enddo
write(wsfile,rec=ire) iscr
C
return
end
C
************************************************************************
subroutine icp(a,b,n)
************************************************************************
* Copies a to b *
************************************************************************
implicit none
integer n,i
integer a(*),b(*)
C
do i=1,n
b(i)=a(i)
enddo
C
return
end
C
************************************************************************
subroutine nlength(nmax,nstr,nnir,isympv,isympo,nactm,nn1,nn2,nn3,
$nn4,nn5,nn6,nn7,nn8,incnew,inewadd,nnewsym,iwan)
************************************************************************
* This subroutine calculates the length of indices running over strings*
************************************************************************
#include "MRCCCOMMON"
integer nnir,nactm,isum,nmax,nn,nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8
integer nstr(nnir,0:nactm,0:nmax,4),irao,isymaoa,isymaob,j1,j2
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 ira,isymav,isymao,irav,isymava,isymavb,incnew(8),inewadd
integer iwan(nnir,nnir,nnir,nnir,nnir),nnewsym,isymi,isyma
C
inewadd=1
do ir=1,nir
isymi=isympair(nnewsym,ir,1)
isyma=isympair(nnewsym,ir,2)
j1=incnew(isymi)
do ira=1,nir
isymav=csympair(isyma,ira,1)
isymao=csympair(isyma,ira,2)
do irav=1,isympv(0,isymav,nn1,nn2,nn3,nn4,1)
isymava=isympv(irav,isymav,nn1,nn2,nn3,nn4,1)
isymavb=isympv(irav,isymav,nn1,nn2,nn3,nn4,2)
nn=nstr(isymava,nn1,nn3,1)*nstr(isymavb,nn2,nn4,2)
do irao=1,isympo(0,isymao,nn5,nn6,nn7,nn8,1)
isymaoa=isympo(irao,isymao,nn5,nn6,nn7,nn8,1)
isymaob=isympo(irao,isymao,nn5,nn6,nn7,nn8,2)
j2=nstr(isymaoa,nn5,nn7,3)*nstr(isymaob,nn6,nn8,4)*nn
iwan(isymi,isymava,isymavb,isymaoa,isymaob)=inewadd
inewadd=inewadd+j1*j2
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine transposition(nmax,nnir,nactm,nstr,icmem,isympv,isympo,
$iwa,iwan,v,isig12,ntoldlen,ntoldleno,namp1,file2,irec2,ircold,ita,
$ntnewlen,incnew,noldsym,nnewsym,nvintoldal,nvintoldbe,nointoldal,
$nointoldbe,nvintoldalact,nvintoldbeact,nointoldalact,nointoldbeact
$,nvintnewalact,nvintnewbeact,nointnewalact,nointnewbeact,
$nsumvirtalact,nsumvirtbeact,nsumoccalact,nsumoccbeact,
$nvirtoldalact,nvirtoldbeact,noccoldalact,noccoldbeact,
$nvirtnewalact,nvirtnewbeact,noccnewalact,noccnewbeact,incsum,lll,
$iwa2,incold,nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,nvintnewal,
$nvintnewbe,nointnewal,nointnewbe,nvirtnewal,nvirtnewbe,noccnewal,
$noccnewbe,nintnew,nintold,nvintold,nointold,nvirtoldal,nvirtoldbe,
$noccoldal,noccoldbe,lsm,liw,ipt,imm,iadd,lcalc,ioffs,ll4,xyz)
************************************************************************
* This subroutine performs transpositions *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,isum,i,j,k,m,itin(8),itrp,j1,j2,j3,j4,iii
integer nvintoldal,nvintoldbe,nointoldal,nointoldbe,iscrn(nmax+2)
integer nstr(nnir,0:nactm,0:nmax,4),albeo(0:4),alben(0:4),nsumnew
integer icmem(nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,4),in(4),kk
integer isympv(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2),namp1
integer isympo(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2),file2
integer ita(nnir,nnir,nnir,nnir),iscro(nmax+2),is1,is2,is3,is4,jj
integer iwa(nnir,nnir,nnir,nnir,nnir,nnir),ntoldlen,ntoldleno,incn
integer iwan(nnir,nnir,nnir,nnir,nnir),irec2,ircold,ntnewlen,nmem
integer isymi,incnew(nnir),incsum(nnir),incold(nnir),ioldadd,ifn
integer inewadd,ir1,noldsym,isyma,inco,iri,iriv,irio,ira,irav,irao
integer nintoldsym,nintoldsymv,nintoldsymo,nintoldsymva,isumsym(4)
integer nintoldsymvb,nintoldsymoa,nintoldsymob,nvirtnewsym,isig12
integer noccnewsym,nvirtnewsyma,nvirtnewsymb,nvirtnewbelen,itinlen
integer noccnewsyma,noccnewsymb,noccnewbelen,nsumsym,irs,irsv,irso
integer nvintoldallen,nvintoldbelen,nointoldallen,nointoldbelen
integer nsumsymv,nsumsymva,nsumsymvb,nsumvirtallen,nsumvirtbelen
integer nsumsymo,nsumsymoa,nsumsymob,nsumoccallen,nsumoccbelen,ifo
integer nintnewsymv,nintnewsymo,nintnewsymva,nintnewsymvb,nnewsym
integer nintnewsymoa,nintnewsymob,nvintnewallen,nvintnewbelen,ir2
integer nointnewallen,nointnewbelen,isumact(4),ism1,ism2,ism3,ism4
integer nvintoldalact,nvintoldbeact,nointoldalact,nointoldbeact
integer nvintnewalact,nvintnewbeact,nointnewalact,nointnewbeact
integer nsumvirtalact,nsumvirtbeact,nsumoccalact,nsumoccbeact,ipt
integer nvirtoldalact,nvirtoldbeact,noccoldalact,noccoldbeact,incs
integer nvirtnewalact,nvirtnewbeact,noccnewalact,noccnewbeact,iadd
integer n1,n2,n3,n4,isa1,isa2,isa3,isa4,iwa2(nnir,nnir,nnir),nnnn
integer nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,nvirtoldlen,xyz
integer noccoldlen,isumind(8),noccoldal,noccoldbe,incint(nnir)
integer nvintnewal,nvintnewbe,nointnewal,nointnewbe,nvirtnewallen
integer nvirtnewal,nvirtnewbe,noccnewal,noccnewbe,in2(4),imm(0:1)
integer nintnew,nintold,nvintold,nointold,nvirtoldal,nvirtoldbe
integer itrpadd(nnir,nnir),icdc,lcalc,ioffs
real*8 v(*)
logical lll,log1(4),lsm,liw,incore,ll4
nsum=nsumvirtal+nsumvirtbe+nsumoccal+nsumoccbe
C Shall we do transposition?
call ifillzero(in,4)
call ifillzero(albeo,5)
call ifillzero(alben,5)
if(nintnew.gt.0) then
isum=0
do j=1,nointoldbe
isum=isum+1
itin(isum)=0
enddo
if(nointoldbe.gt.1) itin(isum)=isum-1
do j=1,nointoldal
isum=isum+1
itin(isum)=-1
enddo
if(nointoldal.gt.1) itin(isum)=isum-1
do j=1,nvintoldbe
isum=isum+1
itin(isum)=-2
enddo
if(nvintoldbe.gt.1) itin(isum)=isum-1
do j=1,nvintoldal
isum=isum+1
itin(isum)=-3
enddo
if(nvintoldal.gt.1) itin(isum)=isum-1
C
isum=0
do j=1,nsumoccbe
isum=isum+1
isumind(isum)=0
enddo
if(nsumoccbe.gt.1) isumind(isum)=isum-1
do j=1,nsumoccal
isum=isum+1
isumind(isum)=-1
enddo
if(nsumoccal.gt.1) isumind(isum)=isum-1
do j=1,nsumvirtbe
isum=isum+1
isumind(isum)=-2
enddo
if(nsumvirtbe.gt.1) isumind(isum)=isum-1
do j=1,nsumvirtal
isum=isum+1
isumind(isum)=-3
enddo
if(nsumvirtal.gt.1) isumind(isum)=isum-1
C
do j=1,nointnewbe
isum=isum+1
isumind(isum)=0
enddo
if(nointnewbe.gt.1) isumind(isum)=isum-1
do j=1,nointnewal
isum=isum+1
isumind(isum)=-1
enddo
if(nointnewal.gt.1) isumind(isum)=isum-1
do j=1,nvintnewbe
isum=isum+1
isumind(isum)=-2
enddo
if(nvintnewbe.gt.1) isumind(isum)=isum-1
do j=1,nvintnewal
isum=isum+1
isumind(isum)=-3
enddo
if(nvintnewal.gt.1) isumind(isum)=isum-1
itrp=0
do i=1,nintold
if(itin(i).ne.isumind(i)) itrp=nintold
enddo
else
itrp=0
endif
if(nintold.gt.1.and.nintnew.gt.0.and.nir.gt.1) itrp=nintold
C Sign of the spin case (not each sp. case of the integrals are stored)
isig12=1
if(nintnew.gt.0.and.(nointold.eq.2.or.nvintold.eq.2)) then
do i=1,nintold-nsum
alben(i)=isumind(nsum+i)
enddo
do i=1,nsum
alben(nintold-nsum+i)=isumind(i)
enddo
do i=1,nintold
albeo(i)=itin(i)
if(albeo(i).gt.0) albeo(i)=albeo(i-1)
if(alben(i).gt.0) alben(i)=alben(i-1)
enddo
if(nointold.eq.2) then
isum=0
do i=1,nintold
if(albeo(i).ge.-1) then
isum=isum+1
iscro(isum)=mod(albeo(i),2)
endif
enddo
isum=0
do i=1,nintold
if(alben(i).ge.-1) then
isum=isum+1
iscrn(isum)=mod(alben(i),2)
endif
enddo
if(iscro(1).ne.iscrn(1)) isig12=-isig12
endif
if(nvintold.eq.2) then
isum=0
do i=1,nintold
if(albeo(i).lt.-1) then
isum=isum+1
iscro(isum)=mod(albeo(i),2)
endif
enddo
isum=0
do i=1,nintold
if(alben(i).lt.-1) then
isum=isum+1
iscrn(isum)=mod(alben(i),2)
endif
enddo
if(iscro(1).ne.iscrn(1)) isig12=-isig12
endif
endif
C Indices to be transposed
if(itrp.ne.0) then
nsumnew=nsum
do i=2,nsum
if(isumind(i).gt.0) nsumnew=nsumnew-1
enddo
do i=2,nintold
j=isumind(i)
if(j.gt.0) then
m=isumind(i-1)
isumind(i-1)=isumind(i-1)-4
itrp=itrp-1
do k=1,nintold
if(itin(k).eq.m) itin(k)=m-4
enddo
endif
enddo
m=0
do i=2,nintold
if(itin(i-1).lt.-3.and.itin(i).gt.0) m=m+1
itin(i)=itin(i+m)
enddo
m=0
do i=2,nintold
if(isumind(i-1).lt.-3.and.isumind(i).gt.0) m=m+1
isumind(i)=isumind(i+m)
enddo
do i=2,itrp
if(isumind(i).gt.0) isumind(i)=i-1
if(itin(i).gt.0) itin(i)=i-1
enddo
C
do i=1,itrp
albeo(i)=itin(i)
alben(i)=isumind(i)
if(albeo(i).gt.0) albeo(i)=albeo(i-1)
if(alben(i).gt.0) alben(i)=alben(i-1)
log1(i)=.true.
enddo
do i=nsumnew,1,-1
j=itrp+1
2015 continue
j=j-1
if(log1(j).and.alben(i).eq.albeo(j)) then
log1(j)=.false.
in(j)=i
else
goto 2015
endif
enddo
do i=itrp,nsumnew+1,-1
j=itrp+1
2016 continue
j=j-1
if(log1(j).and.alben(i).eq.albeo(j)) then
log1(j)=.false.
in(j)=i
else
goto 2016
endif
enddo
endif
if(itrp.gt.1) then
is1=isumind(in(1))
is2=isumind(in(2))
endif
if(itrp.gt.2) is3=isumind(in(3))
if(itrp.gt.3) is4=isumind(in(4))
C
C Length of new free indices
do isymi=1,nir
call slength(nmax,nstr,nnir,isympv,isympo,nactm,incs,isymi,
$nvintnewalact,nvintnewbeact,nvintnewal,nvintnewbe,nointnewalact,
$nointnewbeact,nointnewal,nointnewbe,dsympair)
incnew(isymi)=incs
enddo
C Memory addresses of new intermediates
if(ipt.eq.0)
$call nlength(nmax,nstr,nnir,isympv,isympo,nactm,nvirtnewalact,
$nvirtnewbeact,nvirtnewal,nvirtnewbe,noccnewalact,noccnewbeact,
$noccnewal,noccnewbe,incnew,inewadd,nnewsym,iwan)
if(lsm) ntnewlen=inewadd-1
C Length of summation indices
do isymi=1,nir
call slength(nmax,nstr,nnir,isympv,isympo,nactm,incs,isymi,
$nsumvirtalact,nsumvirtbeact,nsumvirtal,nsumvirtbe,nsumoccalact,
$nsumoccbeact,nsumoccal,nsumoccbe,dsympair)
incsum(isymi)=incs
enddo
C Length of old fixed indices
do isymi=1,nir
call slength(nmax,nstr,nnir,isympv,isympo,nactm,incs,isymi,
$nvirtoldalact,nvirtoldbeact,nvirtoldal,nvirtoldbe,noccoldalact,
$noccoldbeact,noccoldal,noccoldbe,csympair)
incold(isymi)=incs
enddo
ntoldlen=ntoldleno !If the intermediate was skipped
C
if(ipt.gt.0) then
call fwspd(imm,nvintoldal,nointoldal,nvirtoldal,noccoldal,
$nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact,
$nvirtoldalact,noccoldalact,in,ipt,iii)
if(iii.ne.0) then
iadd=-iii
lcalc=0
return
endif
if(lcalc.eq.1) then
if(nvintoldal+nointoldal+nvirtoldal+noccoldal.eq.0.or.
$nvintoldbe+nointoldbe+nvirtoldbe+noccoldbe.eq.0) then
call fwspd(imm,nvintoldbe,nointoldbe,nvirtoldbe,noccoldbe,
$nvintoldbeact,nointoldbeact,nvintoldalact,nointoldalact,
$nvirtoldbeact,noccoldbeact,in,ipt,jj)
if(jj.ne.0) then
iadd=-jj
return
endif
else if(nintold.eq.2) then
in2(1)=in(2)
in2(2)=in(1)
in2(3)=0
in2(4)=0
call fwspd(imm,nvintoldal,nointoldal,nvirtoldal,noccoldal,
$nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact,
$nvirtoldalact,noccoldalact,in2,ipt,jj)
if(jj.ne.0) then
iadd=-jj
lcalc=2
return
endif
endif
lcalc=0
endif
call diadd(imm,nvintoldal,nointoldal,nvirtoldal,noccoldal,
$nvintoldalact,nointoldalact,nvintoldbeact,nointoldbeact,
$nvirtoldalact,noccoldalact,in,ipt,iadd)
endif
C Length of old integral indices
do isymi=1,nir
call slength(nmax,nstr,nnir,isympv,isympo,nactm,inco,isymi,
$nvintoldalact,nvintoldbeact,nvintoldal,nvintoldbe,nointoldalact,
$nointoldbeact,nointoldal,nointoldbe,dsympair)
incint(isymi)=inco
enddo
C Address of old intermediate in case of lambda contracted at the end
if(leftroute.and.namp1.lt.0.and.namp1.ge.-op) then
iii=0
do ir1=1,nir
isymi=isympair(noldsym,ir1,1)
isyma=isympair(noldsym,ir1,2)
itrpadd(isymi,isyma)=iii
iii=iii+incold(isymi)*incint(isyma)
enddo
endif
C Incore or out-of-core algorithm
if(itrp.gt.1) then
inewadd=0
do ir1=1,nir
nintoldsym=isympair(noldsym,ir1,1)
isyma= isympair(noldsym,ir1,2)
isum=0
do ir2=1,nir
isymi= isympair(nintoldsym,ir2,1)
nsumsym=isympair(nintoldsym,ir2,2)
isum=isum+incnew(isymi)*incsum(nsumsym)
enddo
inewadd=inewadd+isum*incold(isyma)
enddo
if(ntoldlen+inewadd.le.maxcor) then
incore=.true.
else
incore=.false.
endif
else
incore=.true.
endif
if(itrp.gt.1) then
nopmax=max(nopmax,dble(ntoldlen))
else
nopmax=max(nopmax,dble(ntnewlen)+dble(ntoldlen))
endif
C Transposition
iii=1
inewadd=1
nnnn=0
C Loop over symmetry cases of old intermediates
do ir1=1,nir
nintoldsym=isympair(noldsym,ir1,1)
isyma= isympair(noldsym,ir1,2)
C Set up increments for old integral indices
call tlength(nmax,nstr,nnir,isympv,isympo,nactm,inco,nintoldsym,
$nvintoldalact,nvintoldbeact,nvintoldal,nvintoldbe,nointoldalact,
$nointoldbeact,nointoldal,nointoldbe,ita,1,0)
nnnn=max(nnnn,inco)
C Transpose free and fixed labels in case of lambda vertex
if(leftroute.and.namp1.lt.0.and.namp1.ge.-op)
$maxmem=max(maxmem,dble(ntoldlen))
C Loop over target indices
do ir2=1,nir
isymi= isympair(nintoldsym,ir2,1)
nsumsym=isympair(nintoldsym,ir2,2)
c do nsumsym=1,nir
c isymi=mult(nsumsym,nintoldsym)
if(.not.lll) iwa2(isymi,nsumsym,isyma)=inewadd
ioldadd=iii
do ira=1,nir
nvirtnewsym=csympair(isyma,ira,1)
noccnewsym=csympair(isyma,ira,2)
do irav=1,isympv(0,nvirtnewsym,nvirtoldalact,nvirtoldbeact,
$nvirtoldal,nvirtoldbe,1)
nvirtnewsyma=isympv(irav,nvirtnewsym,nvirtoldalact,nvirtoldbeact,
$nvirtoldal,nvirtoldbe,1)
nvirtnewallen=nstr(nvirtnewsyma,nvirtoldalact,nvirtoldal,1)
nvirtnewsymb=isympv(irav,nvirtnewsym,nvirtoldalact,nvirtoldbeact,
$nvirtoldal,nvirtoldbe,2)
nvirtnewbelen=nstr(nvirtnewsymb,nvirtoldbeact,nvirtoldbe,2)
nvirtoldlen=nvirtnewallen*nvirtnewbelen
do irao=1,isympo(0,noccnewsym,noccoldalact,noccoldbeact,
$noccoldal,noccoldbe,1)
noccnewsyma=isympo(irao,noccnewsym,noccoldalact,noccoldbeact,
$noccoldal,noccoldbe,1)
noccnewallen=nstr(noccnewsyma,noccoldalact,noccoldal,3)
noccnewsymb=isympo(irao,noccnewsym,noccoldalact,noccoldbeact,
$noccoldal,noccoldbe,2)
noccnewbelen=nstr(noccnewsymb,noccoldbeact,noccoldbe,4)
noccoldlen=noccnewallen*noccnewbelen
if(lll) then
if(liw) then
iwa(isymi,nvirtnewsyma,nvirtnewsymb,noccnewsyma,
$noccnewsymb,nsumsym)=inewadd
else
iwa(isymi,nsumsym,nvirtnewsyma,nvirtnewsymb,noccnewsyma,
$noccnewsymb)=inewadd
endif
endif
C Increments for new indices
incs=incnew(isymi)
incn=incs*incsum(nsumsym)
j4=incn*nvirtoldlen*noccoldlen
jj=inewadd+j4
nmem=ntoldlen+jj
C Loop over summation indices
kk=0
do irs=1,nir
nsumsymv=dsympair(nsumsym,irs,1)
nsumsymo=dsympair(nsumsym,irs,2)
do irsv=1,isympv(0,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,1)
nsumsymva=isympv(irsv,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,1)
nsumvirtallen=nstr(nsumsymva,nsumvirtalact,nsumvirtal,1)
nsumsymvb=isympv(irsv,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,2)
nsumvirtbelen=nstr(nsumsymvb,nsumvirtbeact,nsumvirtbe,2)
j1=nsumvirtallen*nsumvirtbelen
do irso=1,isympo(0,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,1)
nsumsymoa=isympo(irso,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,1)
nsumoccallen=nstr(nsumsymoa,nsumoccalact,nsumoccal,3)
nsumsymob=isympo(irso,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,2)
nsumoccbelen=nstr(nsumsymob,nsumoccbeact,nsumoccbe,4)
j2=j1*nsumoccallen*nsumoccbelen
C Loop over new integral indices
ifn=kk
do iri=1,nir
nintnewsymv=dsympair(isymi,iri,1)
nintnewsymo=dsympair(isymi,iri,2)
do iriv=1,isympv(0,nintnewsymv,nvintnewalact,
$nvintnewbeact,nvintnewal,nvintnewbe,1)
nintnewsymva=isympv(iriv,nintnewsymv,nvintnewalact,nvintnewbeact,
$nvintnewal,nvintnewbe,1)
nvintnewallen=nstr(nintnewsymva,nvintnewalact,nvintnewal,1)
nintnewsymvb=isympv(iriv,nintnewsymv,nvintnewalact,nvintnewbeact,
$nvintnewal,nvintnewbe,2)
nvintnewbelen=nstr(nintnewsymvb,nvintnewbeact,nvintnewbe,2)
nintoldsymva=mult(nintnewsymva,nsumsymva)
nvintoldallen=nstr(nintoldsymva,nvintoldalact,nvintoldal,1)
nintoldsymvb=mult(nintnewsymvb,nsumsymvb)
nvintoldbelen=nstr(nintoldsymvb,nvintoldbeact,nvintoldbe,2)
j3=nvintnewallen*nvintnewbelen
do irio=1,isympo(0,nintnewsymo,nointnewalact,
$nointnewbeact,nointnewal,nointnewbe,1)
nintnewsymoa=isympo(irio,nintnewsymo,nointnewalact,nointnewbeact,
$nointnewal,nointnewbe,1)
nointnewallen=nstr(nintnewsymoa,nointnewalact,nointnewal,3)
nintnewsymob=isympo(irio,nintnewsymo,nointnewalact,nointnewbeact,
$nointnewal,nointnewbe,2)
nointnewbelen=nstr(nintnewsymob,nointnewbeact,nointnewbe,4)
nintoldsymoa=mult(nintnewsymoa,nsumsymoa)
nointoldallen=nstr(nintoldsymoa,nointoldalact,nointoldal,3)
nintoldsymob=mult(nintnewsymob,nsumsymob)
nointoldbelen=nstr(nintoldsymob,nointoldbeact,nointoldbe,4)
C
if(itrp.gt.1.and.incore) then
ifo=ita(nintoldsymva,nintoldsymvb,nintoldsymoa,nintoldsymob)
nintoldlen=nsumvirtallen*nsumvirtbelen*nsumoccallen*nsumoccbelen
$*nvintnewallen*nvintnewbelen*nointnewallen*nointnewbelen
itinlen=nvintoldallen*nvintoldbelen*nointoldallen*nointoldbelen
isum=0
if(nsumoccbe.gt.0) then
isum=isum+1
isumsym(isum)=nsumsymob
isumact(isum)=nsumoccbeact
endif
if(nsumoccal.gt.0) then
isum=isum+1
isumsym(isum)=nsumsymoa
isumact(isum)=nsumoccalact
endif
if(nsumvirtbe.gt.0) then
isum=isum+1
isumsym(isum)=nsumsymvb
isumact(isum)=nsumvirtbeact
endif
if(nsumvirtal.gt.0) then
isum=isum+1
isumsym(isum)=nsumsymva
isumact(isum)=nsumvirtalact
endif
C
if(nointnewbe.gt.0) then
isum=isum+1
isumsym(isum)=nintnewsymob
isumact(isum)=nointnewbeact
endif
if(nointnewal.gt.0) then
isum=isum+1
isumsym(isum)=nintnewsymoa
isumact(isum)=nointnewalact
endif
if(nvintnewbe.gt.0) then
isum=isum+1
isumsym(isum)=nintnewsymvb
isumact(isum)=nvintnewbeact
endif
if(nvintnewal.gt.0) then
isum=isum+1
isumsym(isum)=nintnewsymva
isumact(isum)=nvintnewalact
endif
C
ism1=isumsym(in(1))
ism2=isumsym(in(2))
isa1=isumact(in(1))
isa2=isumact(in(2))
n1=nstr(ism1,isa1,idint(dble(-is1)/4.d0)+1,mod(is1,4)+4)
n2=nstr(ism2,isa2,idint(dble(-is2)/4.d0)+1,mod(is2,4)+4)
if(itrp.eq.3) then
ism3=isumsym(in(3))
isa3=isumact(in(3))
n3=nstr(ism3,isa3,idint(dble(-is3)/4.d0)+1,mod(is3,4)+4)
endif
if(itrp.eq.4) then
ism3=isumsym(in(3))
ism4=isumsym(in(4))
isa3=isumact(in(3))
isa4=isumact(in(4))
n3=nstr(ism3,isa3,idint(dble(-is3)/4.d0)+1,mod(is3,4)+4)
n4=nstr(ism4,isa4,idint(dble(-is4)/4.d0)+1,mod(is4,4)+4)
endif
ifn=ifn+j3*nointnewallen*nointnewbelen
endif
enddo
enddo
enddo
kk=kk+j2*incs
enddo
enddo
enddo
inewadd=jj
ioldadd=ioldadd+inco*nvirtoldlen*noccoldlen
enddo
enddo
enddo
enddo
iii=iii+incold(isyma)*inco
enddo
C Skip old intermediate
if(itrp.gt.1.or.(leftroute.and.namp1.lt.0.and.namp1.ge.-op)) then
if(itrp.gt.1) then
nopmax=max(nopmax,dble(inewadd)+dble(ntoldlen))
maxmem=max(maxmem,dble(inewadd)+dble(nnnn+ibufln))
c if(inewadd+nnnn+ibufln.gt.maxcor)
c $write(6,*) 'M9 ',8.d0*dble(inewadd+nnnn+ibufln)/dble(twoto20)
endif
inewadd=inewadd-1
ntoldlen=inewadd
endif
C
return
end
C
************************************************************************
subroutine savenew(nmem,ntnewlen,wspca,isave,iadd,nvintnewalact,
$nointnewalact,nvintnewbeact,nointnewbeact,nvirtnewalact,
$noccnewalact,intrec,v,iadda,isig11,isig12,irec1,ircnew,file1,
$lcalc,nampvirtalact,nampvirtbeact,nampoccalact,nampoccbeact,
$nvirtnewbeact,noccnewbeact,wspc1,nstr,nmax,isympv,isympo,ita,
$earec,tarec,nampvirtal,nampvirtbe,nampoccal,nampoccbe,nvintnewal,
$nvintnewbe,nointnewal,nointnewbe,nvirtnewal,nvirtnewbe,noccnewal,
$noccnewbe)
************************************************************************
* This subroutine save the new intermediates after the contraction *
************************************************************************
#include "MRCCCOMMON"
integer nmem,ntnewlen,wspca,isave,i,j,k,ia1,ia2,iadd,ssym
integer nvintnewalact,nointnewalact,nvirtnewbeact,noccnewbeact
integer nvintnewbeact,nointnewbeact,nvirtnewalact,noccnewalact
integer nampvirtalact,nampvirtbeact,nampoccalact,nampoccbeact
integer intrec,iadda,isig11,isig12,irec1,ircnew,file1,wspc1
integer nstr,nmax,isympv,isympo,ita,earec,tarec
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe
integer nvintnewal,nvintnewbe,nointnewal,nointnewbe
integer nvirtnewal,nvirtnewbe,noccnewal,noccnewbe
real*8 v(*)
logical lcalc
C
nmem=ntnewlen+1
if(isave.eq.2) then
call fwspc(wspca,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,j,k)
j=intrec+j
call fwspc(iadda,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,ia1,ia2)
if(ia1.gt.0) then
else
call wiadd(iadda,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,ntnewlen)
endif
else
isig12=isig11*isig12
call fwspc(iadd,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,ia1,ia2)
if(ia1.gt.0) then
else
call wiadd(iadd,nvintnewal,nointnewal,nvirtnewal,noccnewal,
$nvintnewalact,nointnewalact,nvintnewbeact,nointnewbeact,
$nvirtnewalact,noccnewalact,ntnewlen)
endif
if(lcalc.and.(nampvirtal.ne.nampvirtbe.or.nampvirtalact.ne.
$nampvirtbeact.or.nampoccalact.ne.nampoccbeact.or.
$nvirtnewal.ne.nvirtnewbe.or.nvirtnewalact.ne.nvirtnewbeact.or.
$noccnewalact.ne.noccnewbeact)) then
maxmem=max(maxmem,2.d0*dble(ntnewlen))
c if(2*ntnewlen.gt.maxcor) write(6,*) 'M10 ',8.d0*dble(2*ntnewlen)/
c $dble(twoto20)
call fwspc(wspc1,nvintnewbe,nointnewbe,nvirtnewbe,noccnewbe,
$nvintnewbeact,nointnewbeact,nvintnewalact,nointnewalact,
$nvirtnewbeact,noccnewbeact,ircnew,k)
call fwspc(iadd,nvintnewbe,nointnewbe,nvirtnewbe,noccnewbe,
$nvintnewbeact,nointnewbeact,nvintnewalact,nointnewalact,
$nvirtnewbeact,noccnewbeact,ia1,ia2)
if(ia1.gt.0) then
else
call wiadd(iadd,nvintnewbe,nointnewbe,nvirtnewbe,noccnewbe,
$nvintnewbeact,nointnewbeact,nvintnewalact,nointnewalact,
$nvirtnewbeact,noccnewbeact,ntnewlen)
endif
endif
endif
C
return
end
C
************************************************************************
subroutine sumsym(nstr,nmax,icmem,isympv,isympo,nnir,ita,nactm,isa
$,nsumsym,nsyma,nnsumlen,nsumvirtalact,nsumvirtbeact,nsumvirtallen,
$nsumvirtbelen,ntampsymva,ntampsymvb,nampvirtalact,nampvirtbeact,
$nampvirtlen,nsumoccalact,nsumoccbeact,nsumoccallen,nsumoccbelen,
$ntampsymoa,ntampsymob,nampoccalact,nampoccbeact,nampocclen,
$ntnewlen,ntoldlen,ntampvirtal,ntampvirtalact,ntampvirtbe,
$ntampvirtbeact,ntampoccal,ntampoccalact,ntampoccbe,ntampoccbeact,
$nslen,nsumvirtal,nsumvirtbe,nsumoccal,nsumoccbe,nampvirtal,
$nampvirtbe,nampoccal,nampoccbe)
************************************************************************
* Symmetry cases of summation indices *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,isa(*),n4,nsumsym,nsyma,ii,nnsumlen,nslen
integer nstr(nnir,0:nactm,0:nmax,4),ita(nnir,nnir,nnir,nnir)
integer icmem(nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,4),ntoldlen
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 irs,irsv,irso,nsumsymv,nsumsymo,nsumsymva,nsumsymvb
integer nsumvirtalact,nsumvirtbeact,nsumvirtallen,nsumvirtbelen
integer nampsymva,nampsymvb,ntampsymva,ntampsymvb,nampvirtalact
integer nampvirtbelen,nampvirtbeact,nampvirtlen,nsumoccalact
integer nsumoccbeact,nsumsymoa,nsumsymob,nsumoccallen,nsumoccbelen
integer nampsymoa,nampsymob,ntampsymoa,ntampsymob,nampoccalact
integer nampoccbelen,nampoccbeact,nampocclen,itadd,ntnewlen
integer ntampvirtal,ntampvirtalact,ntampvirtbe,ntampvirtbeact
integer ntampoccal,ntampoccalact,ntampoccbe,ntampoccbeact
integer nampvirtallen,nampoccallen,nsumvirtal,nsumvirtbe,nsumoccal
integer nsumoccbe,nampvirtal,nampvirtbe,nampoccal,nampoccbe
C Loop over symmetry cases of summation indices
nnsumlen=0
ii=0
nsyma=0
do irs=1,nir
nsumsymv=isympair(nsumsym,irs,1)
nsumsymo=isympair(nsumsym,irs,2)
do irsv=1,isympv(0,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,1)
nsumsymva=isympv(irsv,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,1)
nsumvirtallen=nstr(nsumsymva,nsumvirtalact,nsumvirtal,1)
nsumsymvb=isympv(irsv,nsumsymv,nsumvirtalact,nsumvirtbeact,
$nsumvirtal,nsumvirtbe,2)
nsumvirtbelen=nstr(nsumsymvb,nsumvirtbeact,nsumvirtbe,2)
nampsymva=mult(nsumsymva,ntampsymva)
nampvirtallen=nstr(nampsymva,nampvirtalact,nampvirtal,1)
nampsymvb=mult(nsumsymvb,ntampsymvb)
nampvirtbelen=nstr(nampsymvb,nampvirtbeact,nampvirtbe,2)
nampvirtlen=nampvirtallen*nampvirtbelen
do irso=1,isympo(0,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,1)
nsumsymoa=isympo(irso,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,1)
nsumoccallen=nstr(nsumsymoa,nsumoccalact,nsumoccal,3)
nsumsymob=isympo(irso,nsumsymo,nsumoccalact,nsumoccbeact,
$nsumoccal,nsumoccbe,2)
nsumoccbelen=nstr(nsumsymob,nsumoccbeact,nsumoccbe,4)
nampsymoa=mult(nsumsymoa,ntampsymoa)
nampoccallen=nstr(nampsymoa,nampoccalact,nampoccal,3)
nampsymob=mult(nsumsymob,ntampsymob)
nampoccbelen=nstr(nampsymob,nampoccbeact,nampoccbe,4)
nampocclen=nampoccallen*nampoccbelen
n4=nsumvirtallen*nsumvirtbelen*nsumoccallen*nsumoccbelen
nnsumlen=nnsumlen+n4
itadd=ntnewlen+ntoldlen+
$ita(nampsymva,nampsymvb,nampsymoa,nampsymob)-1
nsyma=nsyma+1
ii=ii+1
isa(ii)=icmem(nsumsymva,ntampsymva,nsumvirtalact,ntampvirtalact,
$nsumvirtal,ntampvirtal,1)-1
ii=ii+1
isa(ii)=icmem(nsumsymvb,ntampsymvb,nsumvirtbeact,ntampvirtbeact,
$nsumvirtbe,ntampvirtbe,2)-1
ii=ii+1
isa(ii)=icmem(nsumsymoa,ntampsymoa,nsumoccalact,ntampoccalact,
$nsumoccal,ntampoccal,3)-1
ii=ii+1
isa(ii)=icmem(nsumsymob,ntampsymob,nsumoccbeact,ntampoccbeact,
$nsumoccbe,ntampoccbe,4)-1
ii=ii+1
isa(ii)=nampvirtallen*nampvirtbelen*nampoccallen*nslen
ii=ii+1
isa(ii)=nampvirtallen*nampvirtbelen*nslen
ii=ii+1
isa(ii)=nampvirtallen*nslen
ii=ii+1
isa(ii)=itadd
ii=ii+1
isa(ii)=nsumvirtallen*nsumvirtbelen*nsumoccallen*nslen
ii=ii+1
isa(ii)=nsumvirtallen*nsumvirtbelen*nslen
ii=ii+1
isa(ii)=nsumvirtallen*nslen
ii=ii+1
isa(ii)=n4*nslen
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine tlength(nmax,nstr,nnir,isympv,isympo,nactm,isum,ssym,
$nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8,ita,inc,isft)
************************************************************************
* This subroutine calculates the length of indices running over strings*
************************************************************************
#include "MRCCCOMMON"
integer nnir,nactm,isum,nmax,nn,nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8
integer nstr(nnir,0:nactm,0:nmax,4),irao,isymaoa,isymaob,ssym,inc
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 ira,isymav,isymao,irav,isymava,isymavb,isft
integer ita(nnir,nnir,nnir,nnir)
C
isum=isft
do ira=1,nir
isymav=dsympair(ssym,ira,1)
isymao=dsympair(ssym,ira,2)
do irav=1,isympv(0,isymav,nn1,nn2,nn3,nn4,1)
isymava=isympv(irav,isymav,nn1,nn2,nn3,nn4,1)
isymavb=isympv(irav,isymav,nn1,nn2,nn3,nn4,2)
nn=nstr(isymava,nn1,nn3,1)*nstr(isymavb,nn2,nn4,2)
do irao=1,isympo(0,isymao,nn5,nn6,nn7,nn8,1)
isymaoa=isympo(irao,isymao,nn5,nn6,nn7,nn8,1)
isymaob=isympo(irao,isymao,nn5,nn6,nn7,nn8,2)
ita(isymava,isymavb,isymaoa,isymaob)=isum*inc
isum=isum+nstr(isymaoa,nn5,nn7,3)*nstr(isymaob,nn6,nn8,4)*nn
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine add2(v1,v2,n,cc)
************************************************************************
* v1 = v1 + cc * v2 *
************************************************************************
#include "MRCCCOMMON"
integer n,i
real*8 v1(n),v2(n),cc
C
if(cc.eq.0.d0) return
maxmem=max(maxmem,2.d0*dble(n))
c if(2*n.gt.maxcor) write(6,*) 'M12 ',8.d0*dble(2*n)/dble(twoto20)
C
return
end
C
************************************************************************
subroutine intsym(nmax,nnir,nactm,nstr,istr,isympv,isympo,iwa,
$nsym,isymi,nactvintal,nactvintbe,nactointal,nactointbe,nvintnewal,
$nvintnewbe,nointnewal,nointnewbe)
************************************************************************
* Symmetry cases for integral lists *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nnir,nactm,isum,ii,iri,iriv,irio,isymiva,isymivb
integer nstr(nnir,0:nactm,0:nmax,-3:0),isymioa,isymiob,nsym,isymi
integer istr(nnir,0:nactm,0:nmax,-3:0),nactvintal,nactvintbe
integer nactointal,nactointbe,isymiv,isymio
integer isympv(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2),iwa(*)
integer isympo(0:nnir,nnir,0:nactm,0:nactm,0:nmax,0:nmax,2)
integer nvintnewal,nvintnewbe,nointnewal,nointnewbe
C
ii=0
nsym=0
do iri=1,nir
isymiv=dsympair(isymi,iri,1)
isymio=dsympair(isymi,iri,2)
do iriv=1,isympv(0,isymiv,nactvintal,nactvintbe,nvintnewal,
$nvintnewbe,1)
isymiva=isympv(iriv,isymiv,nactvintal,nactvintbe,nvintnewal,
$nvintnewbe,1)
isymivb=isympv(iriv,isymiv,nactvintal,nactvintbe,nvintnewal,
$nvintnewbe,2)
do irio=1,isympo(0,isymio,nactointal,nactointbe,nointnewal,
$nointnewbe,1)
isymioa=isympo(irio,isymio,nactointal,nactointbe,nointnewal,
$nointnewbe,1)
isymiob=isympo(irio,isymio,nactointal,nactointbe,nointnewal,
$nointnewbe,2)
nsym=nsym+1
if(nointnewbe.gt.0) then
ii=ii+1
iwa(ii)=nstr(isymiob,nactointbe,nointnewbe,0)
ii=ii+1
iwa(ii)=istr(isymiob,nactointbe,nointnewbe,0)-1
endif
if(nointnewal.gt.0) then
ii=ii+1
iwa(ii)=nstr(isymioa,nactointal,nointnewal,-1)
ii=ii+1
iwa(ii)=istr(isymioa,nactointal,nointnewal,-1)-1
endif
if(nvintnewbe.gt.0) then
ii=ii+1
iwa(ii)=nstr(isymivb,nactvintbe,nvintnewbe,-2)
ii=ii+1
iwa(ii)=istr(isymivb,nactvintbe,nvintnewbe,-2)-1
endif
if(nvintnewal.gt.0) then
ii=ii+1
iwa(ii)=nstr(isymiva,nactvintal,nvintnewal,-3)
ii=ii+1
iwa(ii)=istr(isymiva,nactvintal,nvintnewal,-3)-1
endif
enddo
enddo
enddo
C
return
end
C
************************************************************************
subroutine readft(nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1
$,nointold1,nvirtold1,noccold1,namp1,ical1,isig11,nlist,isave1,
$nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,nvintold1a,nointold1a,
$nvirtold1a,noccold1a,ityp,ii,l2map,log1,intn)
************************************************************************
* Read formula tape *
************************************************************************
#include "MRCCCOMMON"
integer nvintnew1,nointnew1,nvirtnew1,noccnew1,nvintold1,nointold1
integer nvirtold1,noccold1,namp1,ical1,isig11,isave1,nvintnew1a
integer nointnew1a,nvirtnew1a,noccnew1a,nvintold1a,nointold1a,ilst
integer nvirtold1a,noccold1a,ityp(*),nlist,ii,l2map(*),intn(*),ilv
logical log1
C
log1=.true.
if(ii.lt.4) then
read(ftfile,*) nvintnew1,nointnew1,nvirtnew1,noccnew1,
$nvintold1,nointold1,nvirtold1,noccold1,namp1,ical1,isig11,isave1,
$nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,
$nvintold1a,nointold1a,nvirtold1a,noccold1a
else
read(ftfile,*) nvintnew1,nointnew1,nvirtnew1,noccnew1,
$nvintold1,nointold1,nvirtold1,noccold1,namp1,ical1,isig11,nlist,
$isave1,nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,
$nvintold1a,nointold1a,nvirtold1a,noccold1a
endif
c write(6,"(a1,12i3,a1,8i3)")
c $'#',nvintnew1,nointnew1,nvirtnew1,noccnew1,
c $nvintold1,nointold1,nvirtold1,noccold1,namp1,ical1,isig11,isave1,
c $' ',nvintnew1a,nointnew1a,nvirtnew1a,noccnew1a,
c $nvintold1a,nointold1a,nvirtold1a,noccold1a
C
if(((eomgrad.and.(.not.osc)).and.l3route.and.ii.eq.1.and.
$(ityp(ical1).eq.0.or.ityp(ical1).eq.3)).or.
$((pert.ge.1.and.pert.le.3).and.pertroute.eq.0.and.namp1.eq.op).or.
$(pertroute.gt.2.and.pertroute.le.op+2.and.namp1.ne.pertroute-2)
$.or.(pertroute.gt.op+2.and.(namp1.ne.pertroute-3-op.or.
$nvintold1+nointold1+nvirtold1+noccold1.lt.3))) then
log1=.false.
if(isave1.ne.0) then
ilv=ii
call skip1(ilv,ftfile)
endif
endif
C
if(ptroute) then
if(t1route.or.(eomroute.and.ityp(ical1).ne.0))
$ical1=ical1+d2imed
if(t1route.or.l2route) nlist=nlist+d2imed
if(leftroute) then
if(l1route) then
if(ii.eq.1.and.ityp(ical1).eq.0) log1=.false.
if(ityp(ical1).ne.0.and.(ityp(ical1).ne.2.or.ical1.gt.0))
$ical1=ical1+d2imed
endif
if(l2route) then
if(ityp(ical1).eq.0) ical1=ical1+d2imed
if(ityp(ical1).eq.1) ical1=l2map(ical1-eomimed)
if(ityp(ical1).eq.2) log1=.false.
endif
if(l3route) then
if(ityp(ical1).eq.3.and.ical1.le.densimed)
$ical1=ical1+d2imed
endif
endif
if(d1route.and.(ityp(ical1).ne.2.or.ical1.gt.0))
$ical1=ical1+d2imed
if(nmr) then
ilst=ical1
if(ii.eq.4) ilst=nlist
if(ilst.gt.d2imed.and.ilst.le.d2imed+iimed) then
if(intn(ilst-d2imed).ge.12) isig11=-isig11
endif
endif
endif
C
C
return
end
C
************************************************************************
subroutine diag12(itsh,drecm,ii1,ii2,intn,intrec,wspca,denrec,
$wsmax,imed,ld)
************************************************************************
* Integral lists from diagrams 12 and 13 *
************************************************************************
#include "MRCCCOMMON"
integer i,j,ii,jj,isum,itsh,drecm,ii1,ii2,intn(*),intrec(*),wsmax
integer wspca(0:wsmax,13),denrec(*),imed(16,1),iw,nactvirtal
integer nactoccal,nactvirtbe,nactoccbe,ire,ile,k
logical ld
C
do i=itsh+1,itsh+iimed
if(intn(i-itsh).eq.ii1) then
do j=itsh+1,itsh+iimed
jj=j-itsh
if(intn(jj).eq.ii2.and.imed(11,i).eq.imed(9,j).and.
$imed(12,i).eq.imed(10,j)) then
intrec(i)=intrec(j)
isum=0
do iw=1,wspca(0,i)
nactvirtal=wspca((iw-1)*12+9,i)
nactoccal=wspca((iw-1)*12+10,i)
nactvirtbe=imed(11,i)-nactvirtal
nactoccbe=imed(12,i)-nactoccal
call fwspc(wspca(0,j),wspca((iw-1)*12+3,i),wspca((iw-1)*12+4,i),0,
$0,nactvirtal,nactoccal,nactvirtbe,nactoccbe,0,0,ire,ile)
if(ire*ile.gt.0) then
wspca((iw-1)*12+11,i)=ire
wspca((iw-1)*12+12,i)=ile
k=mod(ile,ibufln)
if(k.ne.0) k=ibufln-k
isum=isum+(ile+k)/ibufln
endif
enddo
if(ld.and.(dens.eq.2.or.(dens.eq.1.and.ii1.eq.12))) then
denrec(i-itsh)=drecm
if(ii2.eq.3.or.intn(jj).eq.intn(jj+1)) then
drecm=drecm+intrec(j+1)-intrec(j)
else
drecm=drecm+isum
endif
endif
endif
enddo
endif
enddo
C
return
end
C
************************************************************************
subroutine iaddini(nconf,nmax,nactm,i,iadd,tcase,nptype,wspc,
$wsmax)
************************************************************************
* Initialize iadd vector *
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,nact,nactpl,ii,i,iadd,wsmax,iw,nptype
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer iactv,iacto,iactva,iactoa,iactvb,iactob,wspc(0:wsmax,13)
integer tcase(2,(2*nactm+1)**2,0:nmax)
C
if(densroute) then
ii=i+ptsh
do iw=1,wspc(0,ii)
call wiadd(iadd,wspc((iw-1)*12+1,ii),wspc((iw-1)*12+2,ii),
$wspc((iw-1)*12+3,ii),wspc((iw-1)*12+4,ii),wspc((iw-1)*12+5,ii),
$wspc((iw-1)*12+6,ii),wspc((iw-1)*12+7,ii),wspc((iw-1)*12+8,ii),
$wspc((iw-1)*12+9,ii),wspc((iw-1)*12+10,ii),wspc((iw-1)*12+12,ii))
enddo
else
iactv=tcase(1,nptype,i)
iacto=tcase(2,nptype,i)
do i1=0,i
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
ii=nconf(iactva,iactvb,iactoa,iactob,i1,i)
if(ii.gt.0) then
call wiadd(iadd,0,0,i1,i1,0,0,0,0,iactva,iactoa,ii)
endif
enddo
enddo
enddo
endif
C
return
end
C
************************************************************************
subroutine imedlen(hsym,rsym,lsym,csym,ifi,ila,pts,itypa,itypb,
$wspca,wspcb,wsmax,imed,nmax,nnir,isympv,isympo,nactm,irec,intrec,
$nstr,denrec,irecpt,intn,imedsyma,imedsymb)
************************************************************************
* Lengrth of intermediates after contraction *
************************************************************************
#include "MRCCCOMMON"
integer hsym,rsym,lsym,csym,imsym(0:9),itypa(*),itypb(*),nmax,pts
integer nnir,iw,nstr,imed(16,1),irecpt,ifi,ila,imedsym
integer isympv,isympo,nactm,i,j,k,l,m,ii,irec,intrec(*),wsmax
integer wspca(0:wsmax,13),wspcb(0:wsmax,13),denrec(*),tsym
integer intn(*),imedsyma(*),imedsymb(*)
C
tsym=1
imsym(0)=mult(hsym,tsym)
imsym(1)=mult(lsym,mult(hsym,tsym))
imsym(2)=mult(lsym,tsym)
imsym(3)=mult(hsym,mult(tsym,rsym))
imsym(4)=mult(mult(lsym,hsym),mult(tsym,rsym))
imsym(5)=mult(lsym,mult(tsym,rsym))
imsym(6)=mult(hsym,csym)
imsym(7)=mult(lsym,csym)
C
do i=ifi+1,ila
imedsym=imsym(itypa(i))
imedsyma(i)=imedsym
ii=0
if(i-pts.gt.0.and.i-pts.le.iimed) ii=intn(i-pts)
if(ii.ne.12.and.ii.ne.13) then
j=1
do iw=1,wspca(0,i)
call ilength(nmax,nstr,imed(1,i),nnir,isympv,isympo,nactm,
$wspca((iw-1)*12+1,i),m,imedsym)
wspca((iw-1)*12+11,i)=j
wspca((iw-1)*12+12,i)=m
k=mod(m,ibufln)
if(k.ne.0) k=ibufln-k
j=j+(m+k)/ibufln
enddo
if(intrec(i).ne.0) then
intrec(i)=irec
irec=irec+j-1
if(ii.ne.0.and.(dens.eq.2.or.(dens.eq.1.and.ii.le.3))) then
denrec(i-pts)=intrec(i)-irecpt
drecmax=irec-irecpt
endif
endif
endif
C
imedsym=imsym(itypb(i))
imedsymb(i)=imedsym
j=1
do iw=1,wspcb(0,i)
call ilength(nmax,nstr,imed(5,i),nnir,isympv,isympo,nactm,
$wspcb((iw-1)*12+1,i),m,imedsym)
wspcb((iw-1)*12+11,i)=j
wspcb((iw-1)*12+12,i)=m
k=mod(m,ibufln)
if(k.ne.0) k=ibufln-k
j=j+(m+k)/ibufln
enddo
enddo
C
return
end
C
************************************************************************
subroutine readlab(lab)
************************************************************************
* Find labels in formula file *
************************************************************************
#include "MRCCCOMMON"
character*5 lab,eomcheck
C
eomcheck=' '
do while(eomcheck.ne.lab)
read(ftfile,'(a5)') eomcheck
enddo
C
return
end
C
************************************************************************
subroutine wschange(wspca,wspcb,wsmax,intrec,imedsyma,imedsymb)
************************************************************************
* Cahnges the content of wsfile *
************************************************************************
#include "MRCCCOMMON"
integer i,ii,jj,wsmax,wspca(0:wsmax,13),wspcb(0:wsmax,13)
integer intrec(*),imedsyma(*),imedsymb(*),il1,il2
C
if(calc.eq.0) then
il1=iimed
il2=ccimed
else
il1=eomimed
il2=nimed
endif
do i=il1+1,il2
call wsch(intrec(i),imedsyma(i),imedsymb(i),wspca(0,i),
$wspcb(0,i),wsmax+1,i-il1)
enddo
C
return
end
C
************************************************************************
subroutine updateg(wspcb,wspca,nmax,wsmax,intn,nim,imed,intrec,
$itypa,itypb,l2map,tsp,itcase2,ntcase,almem)
************************************************************************
* This subroutine reads spin cases for intermediates and cluster *
* amplitudes *
************************************************************************
#include "MRCCCOMMON"
integer nmax,wmax,wsmax,i,j,k,i5,i6,i7,i8,i9,i10,nim,itr,ii
integer wspcb(0:wsmax,1),wspca(0:wsmax,1),intn(nim),imed(16,1)
integer intrec(*),itypa(*),itypb(*),l2map(*),iii,jjj,kkk,lll
integer tsp(0:2*(nmax+1)**2,0:nmax),n,m,itcase2,ntcase
real*8 sum,almem
C
rewind(gfile)
read(gfile,*) iii,jjj,kkk,sum,lll
read(gfile,*) nimed,wspcmax,iimed,(icore(ntcase+m),
$m=0,op),ccimed,eomimed,leftimed,densimed,l3imed,d2imed,nimed2,
$ds2imed,dt2imed
i=0
j=0
do m=0,op
do n=1,icore(ntcase+m)
read(gfile,*) i1,i2,ii
icore(itcase+i)=ii
i=i+1
icore(itcase2+j)=i1
j=j+1
icore(itcase2+j)=i2
j=j+1
enddo
enddo
do i=0,op
read(gfile,*) tsp(0,i),(tsp(j,i),j=1,2*tsp(0,i))
enddo
read(gfile,*) (intn(i),i=1,nim)
call ifillzero(wspcb,(wspcmax+1)*nimed2)
call ifillzero(wspca,(wspcmax+1)*nimed2)
do i=1,nimed
read(gfile,*) wmax,itr
wspcb(0,itr)=wmax
do j=1,wmax
read(gfile,*) i1,i2,i3,i4,i5,i6,i7,i8,i9,i10
wspcb((j-1)*12+1,itr)=i1
wspcb((j-1)*12+2,itr)=i2
wspcb((j-1)*12+3,itr)=i3
wspcb((j-1)*12+4,itr)=i4
wspcb((j-1)*12+5,itr)=i5
wspcb((j-1)*12+6,itr)=i6
wspcb((j-1)*12+7,itr)=i7
wspcb((j-1)*12+8,itr)=i8
wspcb((j-1)*12+9,itr)=i9
wspcb((j-1)*12+10,itr)=i10
enddo
read(gfile,*) wmax,itr
wspca(0,itr)=wmax
do j=1,wmax
read(gfile,*) i1,i2,i3,i4,i5,i6,i7,i8,i9,i10
wspca((j-1)*12+1,itr)=i1
wspca((j-1)*12+2,itr)=i2
wspca((j-1)*12+3,itr)=i3
wspca((j-1)*12+4,itr)=i4
wspca((j-1)*12+5,itr)=i5
wspca((j-1)*12+6,itr)=i6
wspca((j-1)*12+7,itr)=i7
wspca((j-1)*12+8,itr)=i8
wspca((j-1)*12+9,itr)=i9
wspca((j-1)*12+10,itr)=i10
enddo
enddo
do i=1,nimed
read(gfile,*) (imed(j,i),j=1,16),intrec(i),itypa(i),itypb(i)
enddo
C Update gfile
close(gfile)
rewind(gfile)
write(gfile,*) iii,jjj,kkk,almem,lll
write(gfile,*) nimed,wspcmax,iimed,(icore(ntcase+m),
$m=0,op),ccimed,eomimed,leftimed,densimed,l3imed,d2imed,nimed2,
$ds2imed,dt2imed
i=0
j=0
do m=0,op
do n=1,icore(ntcase+m)
ii=icore(itcase+i)
i=i+1
i1=icore(itcase2+j)
j=j+1
i2=icore(itcase2+j)
j=j+1
write(gfile,*) i1,i2,ii
enddo
enddo
do i=0,op
write(gfile,*) tsp(0,i),(tsp(j,i),j=1,2*tsp(0,i))
enddo
write(gfile,*) (intn(i),i=1,nim)
do i=1,nimed
wmax=wspcb(0,i)
itr=i
write(gfile,*) wmax,itr
do j=1,wmax
write(gfile,*) (wspcb((j-1)*12+k,itr),k=1,10)
enddo
wmax=wspca(0,i)
itr=i
write(gfile,*) wmax,itr
do j=1,wmax
write(gfile,*) (wspca((j-1)*12+k,itr),k=1,10)
enddo
enddo
do i=1,nimed
write(gfile,*) (imed(j,i),j=1,8),' ',(imed(j,i),j=9,16),' ',
$intrec(i),' ',itypa(i),itypb(i)
enddo
write(gfile,*)
C
return
end
C
#if defined (MPI)
************************************************************************
subroutine psend(file,nconf,nmax,nactm,trec,v,irec,wspc,wsmax,
$denrec,econf,erec)
************************************************************************
* Synchronizing a file an all nodes
************************************************************************
#include "MRCCCOMMON"
integer nmax,nactm,i,k,l,m,n,nn,iw,wsmax,wspc(0:wsmax,13),ii,mm
integer denrec(*)
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer econf(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 iactvb,iactob,iactv,iacto,iactva,iactoa,file,irec
real*8 v(*)
C
if(densroute) then
do i=ptsh+1,ptsh+iimed
ii=i-ptsh
do iw=1,wspc(0,i)
nn=wspc((iw-1)*12+12,i)
mm=denrec(ii)+wspc((iw-1)*12+11,i)
maxmem=max(maxmem,dble(nn)+2.d0*dble(mpibfl))
enddo
enddo
else
do i=0,op
do iactv=max(0,i-mrop),min(nactv,i)
do iacto=max(0,i-mrop),min(nacto,i)
do i1=0,i
nn=0
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(isym.eq.isy1.or.(eomcc.and.iroot.eq.1).or.
$zroute.or.(eomgrad.and.l3route)) then
nn=nconf(iactva,iactvb,iactoa,iactob,i1,i)
mm= trec(iactva,iactvb,iactoa,iactob,i1,i)
else
nn=econf(iactva,iactvb,iactoa,iactob,i1,i)
mm= erec(iactva,iactvb,iactoa,iactob,i1,i)
endif
if(nn.gt.0) then
maxmem=max(maxmem,dble(nn)+2.d0*dble(mpibfl))
endif
enddo
enddo
enddo
enddo
enddo
enddo
endif
C
return
end
C
************************************************************************
subroutine imedsync(file,irec,wspc,v)
************************************************************************
* Synchronize intermediates in files *
************************************************************************
#include "MRCCCOMMON"
integer file,irec,wspc(0:1),iw,ire,ile
real*8 v(*)
C
do iw=1,wspc(0)
ile=wspc((iw-1)*12+12)
maxmem=max(maxmem,dble(ile)+2.d0*dble(mpibfl))
enddo
C
return
end
C
#endif
************************************************************************
subroutine cconv(nconf,trec,nstr,nmax,f,nnir,isympv,isympo,nbmax,
$nactm,noccup,absind,irec,ioffs,roccup,nactd,jconf,jrec,mosym,
$iarcrec,itr,tarec,ita,jstr,absinv,mosymvo)
************************************************************************
* Converts CI coefficients corresponding to different Fermi-vacua *
************************************************************************
#include "MRCCCOMMON"
integer*1 isign,isign1
integer nnir,nmax,nactm,i,j,k,nex,ii,jj,nn,ssym,absind,roccup(*)
integer nconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),ioffs
integer jconf(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax),jn
integer n1,n2,n3,n4,ire,ile,nbmax,l,kk,nactd(*),jampvirtal,jex
integer trec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer jrec(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),nva,nvb,noa,nob,jampvirtbe
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,iob(nactm)
integer nn1,nn2,nn3,nn4,iva(nactm),ivb(nactm),ioa(nactm),mosym(*)
integer nactvintal,nactointal,nactvintbe,nactointbe,nactvirtal
integer nactoccal,nactvirtbe,nactoccbe,irec(nnir,0:nactm,0:nmax,4)
integer nampvirtal,nampvirtbe,nampoccal,nampoccbe,intadd,jmem
integer dbladd,jmva,jmvb,jmoa,jmob,jtva,jtvb,jtoa,jtob
integer jactva,jactvb,jactoa,jactob,iarcrec(0:nactm,0:nmax,4)
integer tarec(0:nactm,0:nactm,0:nactm,0:nactm,0:nmax,0:nmax)
integer itr(2,0:nactm,0:nmax,4),ita(nnir,nnir,nnir,nnir)
integer jstr(nnir,0:nactm,0:nmax,4),absinv(2*nbmax,-3:0)
integer mosymvo(2*nbmax,4)
real*8 f(*)
C Operator string generating the ref det from the absolute Fermi-vacuum
nva=0
nvb=0
noa=0
nob=0
do i=1,nbmax
if(nactd(i).eq.1.and.noccup(i).ne.roccup(i)) then
if((noccup(i).eq.1.or.noccup(i).eq.2).and.
$ (roccup(i).eq.-1.or.roccup(i).eq.0)) then
nva=nva+1
iva(nva)=i
endif
if((noccup(i).eq.-1.or.noccup(i).eq.2).and.
$ (roccup(i).eq.0.or.roccup(i).eq.1)) then
nvb=nvb+1
ivb(nvb)=i
endif
if((noccup(i).eq.-1.or.noccup(i).eq.0).and.
$ (roccup(i).eq.1.or.roccup(i).eq.2)) then
noa=noa+1
ioa(noa)=i
endif
if((noccup(i).eq.0.or.noccup(i).eq.1).and.
$ (roccup(i).eq.-1.or.roccup(i).eq.2)) then
nob=nob+1
iob(nob)=i
endif
endif
enddo
C Read important arrays of the ref det
open(stfile,status='unknown',access='direct',recl=irecln)
rewind(snfile)
read(snfile) jconf
read(snfile)
read(snfile)
read(snfile) itr
read(snfile)
read(snfile)
read(snfile)
read(snfile)
C Loop over excitations
ssym=isym
if(calc.eq.1) ssym=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
nn=nconf(iactva,iactvb,iactoa,iactob,i1,nex)
if(nn.gt.0) then
nampvirtal=i1
nampvirtbe=nex-nampvirtal
nampoccal=i1
nampoccbe=nex-nampoccal
C
do jactva=max(0,iactva-noa),min(nactva,iactva+nva)
do jactvb=max(0,iactvb-nob),min(nactvb,iactvb+nvb)
jampvirtal=nampvirtal-iactva+jactva
jampvirtbe=nampvirtbe-iactvb+jactvb
jactoa=jampvirtal-nampoccal+iactoa
jactob=jampvirtbe-nampoccbe+iactob
jex=jampvirtal+jampvirtbe
if(jex.le.op.and.jampvirtal.ge.0.and.jampvirtbe.ge.0.and.
$jampvirtal.le.min(nvirtal,nal).and.jampvirtbe.le.min(nvirtbe,nbe)
$.and.jactoa.ge.max(0,iactoa-nva).and.
$jactoa.le.min(nactoa,iactoa+noa).and.jactob.ge.max(0,iactob-nvb)
$.and.jactob.le.min(nactob,iactob+nob)) then
jn=jconf(jactva,jactvb,jactoa,jactob,jampvirtal,jex)
if(jn.gt.0) then
C
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)
maxmem=max(maxmem,dble(nn)+dble(jn)+
$dble((nvirtal+1)*(jampvirtal+1)+(nvirtbe+1)*(jampvirtbe+1)+
$(nal+1)*(jampvirtal+1)+(nbe+1)*(jampvirtbe+1)+
$itr(2,jactva,jampvirtal,1)+itr(2,jactvb,jampvirtbe,2)+
$itr(2,jactoa,jampvirtal,3)+itr(2,jactob,jampvirtbe,4)+
$n1*nampvirtal+n2*nampvirtbe+n3*nampoccal+n4*nampoccbe+
$nbasis/2+(nactmax+1)**4*(oo+1)**2+
$(nactmax+1)**4*(oo+1)**2+2*4*(oo+1)*(nactmax+1)+
$4*nir*(oo+1)*(nactmax+1)+8*nbasis)/dble(iintfp))
ii=ii+n3*n4*k
enddo
enddo
enddo
C
endif
endif
enddo
enddo
C
endif
enddo
enddo
enddo
enddo
enddo
enddo
C
return
end
C
Reference in New Issue View Git Blame Copy Permalink