1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18 subroutine dalcip(nela,nelb,ndet,c,ne,trou,part,nd,nm, 19 * inact,iact,nact,lfermi,nisht,nsym,LUN98,mult) 20#include "implicit.h" 21#include "priunit.h" 22 23! ** Before 2011, the nevpt2 code was not compiled unless PGF90, GFORTRAN or IFORT compiler. 24! ** Now code is compiled for all compilers, if your compiler cannot handle 25! ** the code then do NOT #define GO_ON_COMP in dalcip.F (below) and koopro4.F /Jan. 2011 26 27! #if !(defined(VAR_IFORT)||defined(VAR_PGF90)||defined(VAR_GFORTRAN)) 28! return 29! end 30! #else 31#define GO_ON_COMP 32! #endif 33#ifdef GO_ON_COMP 34 dimension c(*),ne(*),trou(*),part(*),nd(*) 35c dimension c(*),ast(nela,*),bst(nelb,*),ne(*),trou(*),part(*),nd(*) 36c integer*1 ast,bst 37 integer*2 ne,trou,part,ispin(600),iorb(600) 38 integer astr(1600),bstr(1600),iecci(0:100) 39 character*1 zwspin(600),zplus,zmoins 40 dimension inact(8),iact(8) 41c integer tri 42c tri(i,j)=max(i,j)*(max(i,j)-1)/2+min(i,j) 43 zplus='+' 44 zmoins='-' 45 ngelo=0 46 ina=0 47 nocc=lfermi+nisht 48 do i=1,nsym 49 ina=ina+inact(i) 50 enddo 51 norb=nm 52 nd(1)=0 53 do i=1,norb 54 ispin(i)=0 55 zwspin(i)=zmoins 56 ispin(i+norb)=1 57 zwspin(i+norb)=zplus 58 iorb(i)=i 59 iorb(i+norb)=i 60 enddo 61c--renzo 62 if(mod((nela-nelb),2).ne.0)then !per un doppietto, quadrupletto,... 63 if(nela.ne.nelb)then 64 ispin(norb+norb+1)=1 65 zwspin(norb+norb+1)=zplus 66 iorb(norb+norb+1)=norb+1 67 endif 68 endif 69 rewind LUN98 70 do idet=1,ndet 71 read(LUN98)c(idet),(astr(i),i=1,nela),(bstr(i),i=1,nelb) 72c remco 73c if (dabs(c(idet)).gt.0.10) then 74c write(2,'(F10.6,5x,14I3)')c(idet),(astr(jj),jj=1,nela) 75c write(2,'(15x,14I3)')(bstr(jj),jj=1,nelb) 76c endif 77c remco 78c if(nela.eq.nelb)then 79 if(mult.eq.1)then 80 call bupa(astr,bstr,norb,nela,nelb,ne,trou,part,nd,ina, 81 $ idet,isign) 82c elseif (mod((nela-nelb),2).eq.1)then !doublets, quartets,... 83 elseif (mod(mult,2).eq.0)then !doublets, quartets,... 84 call bupa2(astr,bstr,norb,nela,nelb,ne,trou,part,nd,ina, 85 $ idet,isign,lfermi) 86 else !triplets, quintets,... 87 call bupa3(astr,bstr,norb,nela,nelb,ne,trou,part,nd,ina, 88 $ idet,isign) 89 endif 90 if(isign.eq.-1)c(idet)=-c(idet) 91 nec=ne(idet) 92 if(idet.lt.ndet)nd(idet+1)=nd(idet)+nec 93 enddo 94 call ascend(nd,ne,trou,part,c,ndet,nm,iecci,nexmax,nocc,norb 95 $ ,mult,iorb,ispin) 96c if(nela.eq.nelb)then 97c if(mult.eq.1)then 98c call orderfe3(nexmax,iecci,ndet,c,nd,ne,trou,part,nm,ina,ast, 99c $ bst,nela) 100c endif 101 return 102 end 103c--------------------------------------------------------------- 104 subroutine bupa2(astr,bstr,norb,nela,nelb,ne,trou,part,nd, 105 $ ina,idet,isign,nacto) 106 dimension astr(*),bstr(*),trou(*),part(*),ne(*),nd(*) 107 dimension c(1600) 108 integer*2 trou,part,ne,istri(20),jstri(20) 109 integer astr,bstr,c 110 logical zhole 111 ip=0 112 irest=0 113 ni=nd(idet) 114 do i=1,20 115 istri(i)=0 116 jstri(i)=0 117 enddo 118 do i=1,nelb 119 if(bstr(i).gt.nacto)then 120 ip=ip+1 121 part(ni+ip)=bstr(i)+ina 122 else 123 irest=irest+1 124 c(irest)=bstr(i) 125 endif 126 enddo 127 it=0 128 do i=1,nacto 129 zhole=.true. 130 do j=1,irest 131 if(c(j).eq.i)then 132 zhole=.false. 133 goto 1 134 endif 135 enddo 136 1 continue 137 if(zhole)then 138 it=it+1 139 trou(ni+it)=i+ina 140 endif 141 enddo 142 irest=0 143 do i=1,nela 144 if(astr(i).gt.nacto)then 145 ip=ip+1 146 part(ni+ip)=astr(i)+ina+norb 147 else 148 irest=irest+1 149 c(irest)=astr(i) 150 endif 151 enddo 152 if(nelb.lt.nela)then 153 ip=ip+1 154 part(ni+ip)=norb+norb+1 !practically infinity 155 endif 156 do i=1,nacto 157 zhole=.true. 158 do j=1,irest 159 if(c(j).eq.i)then 160 zhole=.false. 161 goto 2 162 endif 163 enddo 164 2 continue 165 if(zhole)then 166 it=it+1 167 trou(ni+it)=i+ina+norb 168 endif 169 enddo 170 nec=it 171 ne(idet)=nec 172 do i=1,nacto 173 istri(i)=i !vacuum state definition 174 jstri(i)=i !vacuum state definition 175 enddo 176 imore=0 177 do i=1,nec 178 it=trou(ni+i) 179 ip=part(ni+i) 180 if(it.gt.norb)then 181 if(ip.gt.norb)then 182 istri(it-norb-ina)=ip-norb-ina 183 else 184 istri(it-norb-ina)=0 185 endif 186 else 187 if(ip.gt.norb)then 188 jstri(it-ina)=0 189 imore=imore+1 190 istri(nacto+imore)=ip-norb-ina !for triplets 191 else 192 jstri(it-ina)=ip-ina 193 endif 194 endif 195 enddo 196 isign=1 197 do i=1,nacto+imore !this is for triplets (and for doublets a la cipsi) 198 do j=i+1,nacto+imore !same remark a s above 199 if(istri(i).gt.istri(j))isign=-isign 200 enddo 201 enddo 202 do i=1,nacto 203 do j=i+1,nacto 204 if(jstri(i).gt.jstri(j))isign=-isign 205 enddo 206 enddo 207 return 208 end 209c********************************************** 210 subroutine bupa3(astr,bstr,norb,nela,nelb,ne,trou,part,nd, 211 $ ina,idet,isign) 212c--routine for triplets (high spin) 213 dimension astr(*),bstr(*),trou(*),part(*),ne(*),nd(*) 214 dimension c(1600) 215 integer*2 trou,part,ne,istri(20),jstri(20) 216 integer astr,bstr,c 217 logical zhole 218 ip=0 219 irest=0 220 ni=nd(idet) 221 nel=(nela+nelb)/2 222 do i=1,nelb 223 if(bstr(i).gt.nel)then 224 ip=ip+1 225 part(ni+ip)=bstr(i)+ina 226 else 227 irest=irest+1 228 c(irest)=bstr(i) 229 endif 230 enddo 231 it=0 232 do i=1,nel 233 zhole=.true. 234 do j=1,irest 235 if(c(j).eq.i)then 236 zhole=.false. 237 goto 1 238 endif 239 enddo 240 1 continue 241 if(zhole)then 242 it=it+1 243 trou(ni+it)=i+ina 244 endif 245 enddo 246 irest=0 247 do i=1,nela 248 if(astr(i).gt.nel)then 249 ip=ip+1 250 part(ni+ip)=astr(i)+ina+norb 251 else 252 irest=irest+1 253 c(irest)=astr(i) 254 endif 255 enddo 256 do i=1,nel 257 zhole=.true. 258 do j=1,irest 259 if(c(j).eq.i)then 260 zhole=.false. 261 goto 2 262 endif 263 enddo 264 2 continue 265 if(zhole)then 266 it=it+1 267 trou(ni+it)=i+ina+norb 268 endif 269 enddo 270 nec=it 271 ne(idet)=nec 272 do i=1,nel 273 istri(i)=i !vacuum state definition 274 jstri(i)=i !vacuum state definition 275 enddo 276 imore=0 277 do i=1,nec 278 it=trou(ni+i) 279 ip=part(ni+i) 280 if(it.gt.norb)then 281 if(ip.gt.norb)then 282 istri(it-norb-ina)=ip-norb-ina 283 else 284 istri(it-norb-ina)=0 285 endif 286 else 287 if(ip.gt.norb)then 288 jstri(it-ina)=0 289 imore=imore+1 290 istri(nel+imore)=ip-norb-ina !for triplets 291 else 292 jstri(it-ina)=ip-ina 293 endif 294 endif 295 enddo 296 isign=1 297 do i=1,nel+imore !this is for triplets (and for doublets a la cipsi) 298 do j=i+1,nel+imore !same remark a s above 299 if(istri(i).gt.istri(j))isign=-isign 300 enddo 301 enddo 302 do i=1,nel 303 do j=i+1,nel 304 if(jstri(i).gt.jstri(j))isign=-isign 305 enddo 306 enddo 307 return 308 end 309c************************************************************** 310 subroutine bupa(astr,bstr,norb,nela,nelb,ne,trou,part,nd, 311 $ ina,idet,isign) 312 dimension astr(*),bstr(*),trou(*),part(*),ne(*),nd(*) 313 dimension c(1600) 314 integer*2 trou,part,ne,istri(20),jstri(20) 315 integer astr,bstr,c 316 logical zhole 317 ip=0 318 irest=0 319 ni=nd(idet) 320 do i=1,nelb 321 if(bstr(i).gt.nela)then 322 ip=ip+1 323 part(ni+ip)=bstr(i)+ina 324 else 325 irest=irest+1 326 c(irest)=bstr(i) 327 endif 328 enddo 329 it=0 330 do i=1,nela 331 zhole=.true. 332 do j=1,irest 333 if(c(j).eq.i)then 334 zhole=.false. 335 goto 1 336 endif 337 enddo 338 1 continue 339 if(zhole)then 340 it=it+1 341 trou(ni+it)=i+ina 342 endif 343 enddo 344 irest=0 345 do i=1,nela 346 if(astr(i).gt.nela)then 347 ip=ip+1 348 part(ni+ip)=astr(i)+ina+norb 349 else 350 irest=irest+1 351 c(irest)=astr(i) 352 endif 353 enddo 354 do i=1,nela 355 zhole=.true. 356 do j=1,irest 357 if(c(j).eq.i)then 358 zhole=.false. 359 goto 2 360 endif 361 enddo 362 2 continue 363 if(zhole)then 364 it=it+1 365 trou(ni+it)=i+ina+norb 366 endif 367 enddo 368 nec=it 369 ne(idet)=nec 370 do i=1,nela 371 istri(i)=i !vacuum state definition 372 jstri(i)=i !vacuum state definition 373 enddo 374 do i=1,nec 375 it=trou(ni+i) 376 ip=part(ni+i) 377 if(it.gt.norb)then 378 istri(it-norb-ina)=ip-norb-ina 379 else 380 jstri(it-ina)=ip-ina 381 endif 382 enddo 383 isign=1 384 do i=1,nela 385 do j=i+1,nela 386 if(istri(i).gt.istri(j))isign=-isign 387 enddo 388 enddo 389 do i=1,nela 390 do j=i+1,nela 391 if(jstri(i).gt.jstri(j))isign=-isign 392 enddo 393 enddo 394 return 395 end 396c**************************************************************** 397 subroutine orderfe(ifirst,ilast,c,ne,nd,trou,part,npos) 398#include "implicit.h" 399 DIMENSION c(*) 400 integer*2 ne(*),trou(*),part(*) 401 dimension nd(*),npos(*) 402 do idet=ifirst,ilast 403 nec=ne(idet) 404 inu=idet-ifirst+1 405 if(npos(inu).eq.inu) goto 1 406 do jdet=idet+1,ilast 407 jnu=jdet-ifirst+1 408 if(npos(jnu).eq.inu)then 409 call swapbp(nec,idet,jdet,nd,trou,part) 410 cdum=c(idet) 411 c(idet)=c(jdet) 412 c(jdet)=cdum 413 ndum=npos(inu) 414 npos(inu)=npos(jnu) 415 npos(jnu)=ndum 416 endif 417 enddo 418 1 continue 419 enddo 420 return 421 end 422c---------------------------------------------------------- 423 subroutine giveoccij(iocc,m,nd,ne,trou,part,nocc,norb) 424 integer*2 ne(*),trou(*),part(*) 425 integer*1 iocc(*) 426 integer nd(*) 427 do i=1,nocc 428 iocc(i)=2 429 enddo 430 do i=nocc+1,norb 431 iocc(i)=0 432 enddo 433 nstart=nd(m) 434 do i=1,ne(m) 435 ib=trou(nstart+i) 436 ip=part(nstart+i) 437 if(ib.gt.norb)ib=ib-norb 438 if(ip.gt.norb)ip=ip-norb 439 iocc(ib)=iocc(ib)-1 440 if (ip.le.norb) then 441 iocc(ip)=iocc(ip)+1 442 endif 443 enddo 444 return 445 end 446C------------------------------------------------------------------ 447 integer function ndiffij(iocc,jocc,norb) 448 integer*1 iocc(*),jocc(*) 449 ndiffij=0 450 do i=1,norb 451 ia=iocc(i)-jocc(i) 452 if(ia.ne.0)then 453 ndiffij=1 454 return 455 endif 456 enddo 457 return 458 end 459c************************************************************* 460 subroutine swapbp(nec,idet,jdet,nd,trou,part) 461 dimension nd(*),trou(*),part(*) 462 dimension troud(100),partd(100) 463 integer*2 trou,part,troud,partd 464 ndi=nd(idet) 465 do i=1,nec 466 troud(i)=trou(ndi+i) 467 partd(i)=part(ndi+i) 468 enddo 469 ndj=nd(jdet) 470 do i=1,nec 471 trou(ndi+i)=trou(ndj+i) 472 part(ndi+i)=part(ndj+i) 473 trou(ndj+i)=troud(i) 474 part(ndj+i)=partd(i) 475 enddo 476 return 477 end 478c****************************************************************** 479 subroutine givepos(iel,n,k,npos) 480c 481c n=numero di orbitali singoli occupati 482c k=numero di elettroni spaiati alpha 483c iel(i),i=1,k = posizione dell'i-esimo elettrone 484c nelle n caselle 485c 486 integer iel 487 dimension iel(*) 488 logical*1 zsz0 489 490 zsz0=n.eq.k*2 491 492 if (zsz0) then 493 ntot=1 494 do i=k+1,n 495 ntot=ntot*i 496 enddo 497 do i=1,k 498 ntot=ntot/i 499 enddo 500 endif 501 502 if (k.eq.1) then 503 npos=iel(1) 504 goto 110 505 endif 506 507 npos=1 508 kpos=0 509 kposold=0 510 do ijk=1,k-1 511 kpos=iel(ijk)-ijk-kposold 512 if (kpos.eq.0) goto 111 513 nogg=k-ijk 514 do kl=ijk+kposold,iel(ijk)-1 515 ncasel=n-kl 516 ndistr=1 517 do l=nogg+1,ncasel 518 ndistr=ndistr*l 519 enddo 520 do l=1,ncasel-nogg 521 ndistr=ndistr/l 522 enddo 523 npos=npos+ndistr 524 enddo 525 kposold=kposold+kpos 526 111 continue 527 enddo 528 npos=npos+iel(k)-iel(k-1)-1 529 530 531 110 if (zsz0) then 532 if (npos*2.le.ntot) then 533 npos=npos*2-1 534 else 535 npos=(ntot-npos)*2+2 536 endif 537 endif 538 539 return 540 end 541c--------------------------------------------------------------- 542 subroutine ascend(nd,ne,trou,part,c,ndet,nm,iecci,nexmax,nocc,norb 543 $ ,mult,iorb,ispin) 544#include "implicit.h" 545 DIMENSION c(*) 546 dimension nd(*),ne(*),trou(*),part(*),iorb(*),ispin(*) 547 integer*2 ne,trou,part,ibu(100),ipa(100),iorb,ispin 548 integer iecci(0:100) 549 integer*1 o 550 integer*1 occu 551 allocatable occu(:,:),iel(:),npos(:) 552 nexmax=0 553 do i=0,100 554 iecci(i)=0 555 enddo 556 do i=1,ndet 557 if(ne(i).gt.nexmax) nexmax=ne(i) 558 enddo 559 LUN15=-13300 560 CALL GPOPEN(LUN15,'scr15','UNKNOWN',' ','UNFORMATTED',IDUMMY, 561 * .FALSE.) 562 do iec=0,nexmax 563 do idet=1,ndet 564 nec=ne(idet) 565 if(nec.eq.iec)then 566 iecci(iec)=iecci(iec)+1 567 ndi=nd(idet) 568 write(lun15)c(idet),nec,(trou(i),part(i),i=ndi+1,ndi+nec) 569 endif 570 enddo 571 enddo 572 rewind lun15 573 nd(1)=0 574 idet1=1 575 allocate(occu(norb,ndet)) 576 allocate(iel(norb)) 577c allocate(npos(norb)) !wrong should be maximum of a deter space 578 allocate(npos(2000)) !covers more than 12-open shell 579 do iec=0,nexmax 580 do idet=idet1,idet1+iecci(iec)-1 581 read(lun15)c(idet),nec,(ibu(i),ipa(i),i=1,nec) 582 ne(idet)=nec 583 if(idet.lt.ndet)nd(idet+1)=nd(idet)+nec 584 ndi=nd(idet) 585 do i=1,nec 586 trou(ndi+i)=ibu(i) 587 part(ndi+i)=ipa(i) 588 enddo 589 call giveoccij(occu(1,idet),idet,nd,ne,trou,part,nocc,norb) 590 enddo 591 ifirst=idet1 592 do idet=idet1,idet1+iecci(iec)-1 593 nec=ne(idet) 594 inext=idet+1 595 if(idet.eq.ifirst)then 596 nsing=0 597 do i=1,norb 598 if(occu(i,idet).eq.1)nsing=nsing+1 599 enddo 600 isz2=mult-1 601 nha=(nsing-isz2)/2 602 nde=noverk(nsing,nha) 603 ilast=ifirst+nde-1 604 ifound=1 605 endif 606 if(idet.gt.ifirst.and.idet.le.ilast)goto 999 607 do jdet=idet+1,idet1+iecci(iec)-1 608 if(ndiff(occu(1,idet),occu(1,jdet),norb).eq.0)then 609 ifound=ifound+1 610 if(jdet.ne.inext)then 611 call swapbp(nec,inext,jdet,nd,trou 612 $ ,part) 613 cdum=c(inext) 614 c(inext)=c(jdet) 615 c(jdet)=cdum 616 do io=1,norb 617 o=occu(io,inext) 618 occu(io,inext)=occu(io,jdet) 619 occu(io,jdet)=o 620 enddo 621 endif 622 inext=inext+1 623 endif 624 if(ifound.eq.nde)goto 999 625 enddo 626 999 if(idet.eq.ilast)then 627 do kdet=ifirst,ilast 628 isingle=0 629 indalf=0 630 do k=1,norb 631 if(occu(k,kdet).eq.1)then 632 isingle=isingle+1 633 ndk=nd(kdet) 634 do i=1,ne(kdet) 635 ispb=ispin(trou(i+ndk)) 636 iorbb=iorb(trou(i+ndk)) 637 ispp=ispin(part(i+ndk)) 638 iorbp=iorb(part(i+ndk)) 639 if((iorbp.eq.k.and.ispp.eq.1).or. 640 $ (iorbb.eq.k.and.ispb.eq.0))then 641 indalf=indalf+1 642 iel(indalf)=isingle 643 if(indalf.gt.norb)then 644 print*,'ERROR: indalf .gt. norb',indalf,norb 645 CALL QUIT('givepos: indalf .gt. norb') 646 endif 647 endif 648 enddo 649 endif 650 enddo 651 if((kdet-ifirst+1).gt.2000)then 652 print*,'givepos: arg of npos too large' 653 CALL QUIT('givepos: arg of npos too large') 654 endif 655 if (indalf.gt.0) then 656 call givepos(iel,isingle,indalf,npos(kdet-ifirst+1)) 657 else 658 npos(kdet-ifirst+1) = -999 999 999 659 end if 660 enddo 661 call orderfe(ifirst,ilast,c,ne,nd,trou,part,npos) 662 ifirst=ilast+1 663 endif 664 enddo 665 idet1=idet1+iecci(iec) 666 enddo 667 CALL GPCLOSE(LUN15,'DELETE') 668c close(15,status='DELETE') 669 return 670 end 671c**************************************************************** 672 integer function noverk(n,k) 673 if(k.eq.0)then 674 noverk=1 675 return 676 elseif(k.eq.1)then 677 noverk=n 678 return 679 else 680 num=n 681 iden=1 682 do i=2,k 683 num=num*(n-i+1) 684 iden=iden*i 685 enddo 686 noverk=num/iden 687 return 688 endif 689 end 690#endif 691