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