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