1      subroutine tce_mrcc_create_cas1(rtdb)
2        implicit none
3#include "global.fh"
4#include "rtdb.fh"
5#include "mafdecls.fh"
6#include "sym.fh"
7#include "util.fh"
8#include "stdio.fh"
9#include "errquit.fh"
10#include "tce.fh"
11#include "tce_mrcc.fh"
12#include "tce_main.fh"
13#include "geom.fh"
14
15        integer iactel,iactorb
16        integer rtdb
17ckbn gf -2      integer nodezero,nprocs
18      logical nodezero
19        character*4 irrepname
20        integer i,j,irrep
21c        integer nref,iref
22        integer iref
23        integer k,n
24        integer id,ia,ib
25        integer mrcc_factorial
26        external mrcc_factorial
27        integer itotal
28        integer iacounts(maxorb)
29        integer ibcounts(maxorb)
30        integer iapos(maxorb+1)
31        integer ibpos(maxorb+1)
32c        integer itotreal,idbocc(10240,maxorb)
33        integer itotreal
34        character*10240 s
35c        character*maxref s
36        integer iactshift
37        character*8 t1
38        character*3 t2
39        integer ic
40        integer l_idbocc,k_idbocc
41c        integer l_sidbocc,k_sidbocc
42
43        nodezero = (ga_nodeid().eq.0)
44
45        do k=1,10240
46          s(k:k)=''
47        enddo
48
49        t1 = 'bwcc:ref'
50
51      if (.not.rtdb_get(rtdb,'mrcc:iactel',mt_int,1,iactel))
52     2   call  errquit('tce_mrcc_cas',2,UNKNOWN_ERR)
53      if (.not.rtdb_get(rtdb,'mrcc:iactorb',mt_int,1,iactorb))
54     2    call  errquit('tce_mrcc_cas',3,UNKNOWN_ERR)
55
56
57      if (.not.ma_push_get(mt_int,maxorb*maxref,"idbocc",
58     +  l_idbocc,k_idbocc))
59     +  call errquit("tce_mrcc_energy: MA problem",0,MA_ERR)
60c       do i=0,maxorb*10240
61c        int_mb(k_idbocc + i-1 ) = 0
62c       enddo
63c      call ifill(maxorb*maxref,0,int_mb(k_idbocc),1)
64
65c      if (.not.ma_push_get(mt_int,maxorb*maxref,"sidbocc",
66c      if (.not.ma_push_get(mt_byte,maxorb*maxref,"sidbocc",
67c     +  l_sidbocc,k_sidbocc))
68c     +  call errquit("tce_mrcc_energy: MA problem",0,MA_ERR)
69c      call ifill(maxorb*maxref,0,int_mb(k_sidbocc),1)
70
71        if(nodezero) then
72
73          write(LuOut,"(/,'Active elec :',I4)")iactel
74          write(LuOut,"('Active orbs :',I4,/)")iactorb
75
76          if(mod(iactel,2).eq.1)
77     1 write(LuOut,"('Odd number of elec.',I4,/)")iactel/2
78
79        endif
80
81        itotal = 0
82
83        do i=1,min(iactel/2,iactorb)
84          id=iactel/2-i+1
85          ia=i-1+mod(iactel,2)
86          ib=i-1
87          itotal = itotal +
88     1 mrcc_factorial(iactorb)/(mrcc_factorial(id)*
89     4 mrcc_factorial(ia)*mrcc_factorial(ib)*
90     5 mrcc_factorial(iactorb-id-ia-ib))
91
92c        if(nodezero)
93c     1  write(LuOut,"('Category ',I5,I5,I5,' has ',I8,
94c     2 ' possibilities')")id,ia,ib,
95c     3 mrfactorial(iactorb)/(mrfactorial(id)*
96c     4 mrfactorial(ia)*mrfactorial(ib)*
97c     5 mrfactorial(iactorb-id-ia-ib))
98
99        enddo
100
101
102        do k=1,iactorb
103            iacounts(k) = 0
104            ibcounts(k) = 0
105          if(k.le.id) then
106            iacounts(k) = 1
107            ibcounts(k) = 1
108          else
109           if(k.le.(id+ia)) iacounts(k) = 1
110           if(k.le.(id+ib)) ibcounts(k) = 1
111          endif
112        enddo
113
114        if(nodezero) then
115          write(LuOut,"(/,'Alpha ',100I1)")(iacounts(k),k=1,iactorb)
116          write(LuOut,"('Beta  ',100I1,/)")(ibcounts(k),k=1,iactorb)
117        endif
118
119        do k=1,id+ib
120          ibpos(k) = k
121        enddo
122
123          ibpos(id+ib+1) = iactorb+1
124
125        itotreal = 0
126 8746   continue
127
128        do i=1,id+ib
129             if(i.gt.1) then
130               if(ibpos(id+ib-i+1).lt.(iactorb-i+1)) then
131             ibpos(id+ib-i+1)=ibpos(id+ib-i+1)+1
132                 do n=1,i-1
133                   ibpos(id+ib-n+1)=ibpos(id+ib-i+1)+i-n
134                 enddo
135               goto 8746
136               else
137               goto 8745
138               endif
139             endif
140         do k=ibpos(id+ib-i+1),ibpos(id+ib-i+2)
141           if(k.eq.(ibpos(id+ib-i+2))) then
142            if((ibpos(id+ib-i+1)+1).lt.(ibpos(id+ib-i+2))) then
143             ibpos(id+ib-i+1)=ibpos(id+ib-i+1)+1
144             if(i.gt.1) goto 8746
145            endif
146            goto 8745
147           endif
148           do n=1,iactorb
149            ibcounts(n) = 0
150           enddo
151           do n=1,id+ib
152            if(n.eq.(id+ib+1-i)) then
153            ibcounts(k) = 1
154            else
155            ibcounts(ibpos(n)) = 1
156            endif
157           enddo
158c           if(nodezero)write(LuOut,"(100I1)")(ibcounts(n),n=1,iactorb)
159           itotreal = itotreal + 1
160           do n=1,iactorb
161c             idbocc(itotreal,n)=ibcounts(n)
162            int_mb(k_idbocc+(n-1)*maxorb+itotreal-1)=ibcounts(n)
163           enddo
164         enddo
165 8745    continue
166        enddo
167
168corg        if(ibpos(1).lt.(iactorb-id-ib)) goto 8746
169         do i=1,(nocc(1)+nocc(2)-iactel)/2
170           s(i:i)='2'
171c           byte_mb(k_sidbocc+(i-1)*maxorb+i-1)='2'
172c       write(*,'(A,A,A)') "s",s(i:i),byte_mb(k_sidbocc+(i-1)*maxorb+i-1)
173c       write(*,'(A,A)') "s",byte_mb(k_sidbocc+(i-1)*maxorb+i-1)
174c           write(*,'(A,A,A)') "s",s
175         enddo
176
177         iactshift = (nocc(1)+nocc(2)-iactel)/2
178
179c         write(6,*)iactshift,nocc(1),nocc(2)
180
181         ic = 0
182
183         do i=1,itotreal
184          do k=1,itotreal
185          irrep = 0
186           do n=1,iactorb
187c          if((idbocc(i,n)+idbocc(k,n)).eq.2) then
188           if(((int_mb(k_idbocc+(n-1)*maxorb+i-1))+
189     +            int_mb(k_idbocc+(n-1)*maxorb+k-1)).eq.2) then
190           s(n+iactshift:n+iactshift)='2'
191c         byte_mb(k_sidbocc+((n+iactshift)-1)*maxorb+(n+iactshift)-1)='2'
192c             elseif(idbocc(i,n).eq.1) then
193             elseif(int_mb(k_idbocc+(n-1)*maxorb+i-1).eq.1) then
194               s(n+iactshift:n+iactshift)='a'
195c         byte_mb(k_sidbocc+((n+iactshift)-1)*maxorb+(n+iactshift)-1)='a'
196               irrep = ieor(irrep,int_mb(k_irs(1)+n+iactshift-1))
197c             elseif(idbocc(k,n).eq.1) then
198             elseif((int_mb(k_idbocc+(n-1)*maxorb+k-1)).eq.1) then
199               s(n+iactshift:n+iactshift)='b'
200c         byte_mb(k_sidbocc+((n+iactshift)-1)*maxorb+(n+iactshift)-1)='b'
201               irrep = ieor(irrep,int_mb(k_irs(2)+n+iactshift-1))
202             else
203               s(n+iactshift:n+iactshift)='0'
204c         byte_mb(k_sidbocc+((n+iactshift)-1)*maxorb+(n+iactshift)-1)='0'
205             endif
206           enddo
207           call sym_irrepname(geom,irrep+1,irrepname)
208        if(targetsym.eq.irrepname) then
209
210          ic = ic + 1
211
212             if(nodezero)write(LuOut,*)s(1:iactorb+iactshift),' ',
213     1 ic
214c      if(nodezero) write(LuOut,'(A,A,I5)')
215c     +byte_mb(k_sidbocc+((iactorb+iactshift)-1)*maxorb+(1)-1),'ii ',ic
216
217      write(t2,"(I3.3)")ic
218      if (.not.rtdb_cput(rtdb,t1//t2,1,s(1:iactorb+iactshift)))
219c      if (.not.rtdb_cput(rtdb,t1//t2,1,byte_mb(k_sidbocc)))
220     1   call errquit('tce_mrcc_cas: failed writing to rtdb',0,
221     2   RTDB_ERR)
222
223           endif
224
225           enddo
226         enddo
227
228        if(nodezero)
229     1  write(LuOut,"(/,'Total no. of references:',I5,/)")ic
230
231      if (.not.rtdb_put(rtdb,'bwcc:nref',mt_int,1,ic))
232     1   call errquit('tce_mrcc_input: failed writing to rtdb',0,
233     2   RTDB_ERR)
234
235c       this routine has to be rewritten at some point
236
237ckbn    check whether number of references go beyond maxref
238        if(ic .gt. maxref)  then
239        if(nodezero)
240     +   write(LuOut,'(A,I5,A,I5)') "Number of references: ",ic,
241     +   " greater than maximum number of references: " , maxref
242        if(nodezero) call util_flush(LuOut)
243        call errquit('tce_mrcc_create_cas: nref > maxref',3,RTDB_ERR)
244        endif
245
246c      if (.not.ma_pop_stack(l_sidbocc))
247c     +    call errquit("tce_mrcc_energy: MA problem",615,MA_ERR)
248
249      if (.not.ma_pop_stack(l_idbocc))
250     +    call errquit("tce_mrcc_energy: MA problem",615,MA_ERR)
251
252        return
253        end
254
255        integer function mrcc_factorial(input)
256         implicit none
257#include "errquit.fh"
258         integer input,res
259         integer i
260
261c         if(input .lt.  0)
262c     + call errquit("mrcc_factorial: cas input",30,INPUT_ERR)
263c
264c         if(input .gt. 20)
265c     + call errquit("mrcc_factorial: cas input",30,INPUT_ERR)
266
267
268         res = 1
269         do i=1,input
270           res = res*i
271         enddo
272
273         if(res .le. 0)
274     + call errquit("mrcc_factorial: cas input",30,INPUT_ERR)
275         mrcc_factorial = res
276
277        end
278
279
280c $Id$
281