subroutine setprtvec implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension itmpc(maxsub), itmpb1(maxsub), itmpb2(maxsub) logical notb1, notb2 nprsub = 0 do i=1,maxsub i2 = 2*i i1 = i2-1 eprsub(i1) = 9.99e19 eprsub(i2) = -9.99e19 enddo C-RDC write(iout,'(//" FINAL EIGENVECTORS:", C-RDC $ /" ------------------")') do id=1,nprvec id2=2*id id1=id2-1 ia1 = iprvec(id1) ia2 = iprvec(id2) nc = 0 nb1 = 0 nb2 = 0 do 100 is=1,nsub notb1 = .true. notb2 = .true. j1 = iatom1(is) j2 = iatom1(is+1)-1 do j=j1,j2 if (iatoms(j).gt.ia2) goto 100 C . Since core atoms have the highest priority, C . we keep on looping over all atoms until a core C . atom is found. We don't have to check for buffer2 C . atoms if a buffer1 atom was already found, since C . buffer1 atoms have a higher priority than buffer2 C . atoms. if ((iatoms(j).ge.ia1).and.(iatoms(j).le.ia2)) then if (iabuff(j).eq.0) then nc = nc + 1 itmpc(nc) = is goto 100 elseif (notb1.and.(iabuff(j).eq.1)) then nb1 = nb1 + 1 itmpb1(nb1) = is notb1 = .false. notb2 = .false. elseif (notb2.and.(iabuff(j).eq.2)) then nb2 = nb2 + 1 itmpb2(nb2) = is notb2 = .false. endif endif enddo 100 continue C . some subsystems might now be listed twice (or more), C . but for the PRTVEC we only want to access each subsystem C . once (at most). C . so do some sorting now, to get a unique list: C . core has the highest priority, then inner buffer, then C . outer buffer C . if the same subsystem occurs in two or three groups, only C . the entry in the group with highest priority will survive, C . e.g. subsystem 6 is both in core (itmpc) and outer buffer C . (itmpb2), then the entry in the outer buffer will be deleted C . A. Compare outer and inner buffer kb1 = 1 kb2 = 1 110 if (kb2.le.nb2) then 120 if (itmpb1(kb1).lt.itmpb2(kb2)) then kb1 = kb1 + 1 if (kb1.gt.nb1) goto 150 goto 120 elseif (itmpb1(kb1).gt.itmpb2(kb2)) then kb2 = kb2 +1 goto 110 else do i=kb2,nb2-1 itmpb2(i) = itmpb2(i+1) enddo nb2 = nb2 - 1 goto 110 endif endif C . B. Compare outer buffer and core 150 kc = 1 kb2 = 1 210 if (kb2.le.nb2) then 220 if (itmpc(kc).lt.itmpb2(kb2)) then kc = kc + 1 if (kc.gt.nc) goto 250 goto 220 elseif (itmpc(kc).gt.itmpb2(kb2)) then kb2 = kb2 +1 goto 210 else do i=kb2,nb2-1 itmpb2(i) = itmpb2(i+1) enddo nb2 = nb2 - 1 goto 210 endif endif C . C. Compare inner buffer and core 250 kc = 1 kb1 = 1 310 if (kb1.le.nb1) then 320 if (itmpc(kc).lt.itmpb1(kb1)) then kc = kc + 1 if (kc.gt.nc) goto 350 goto 320 elseif (itmpc(kc).gt.itmpb1(kb1)) then kb1 = kb1 +1 goto 310 else do i=kb1,nb1-1 itmpb1(i) = itmpb1(i+1) enddo nb1 = nb1 - 1 goto 310 endif endif 350 continue C-RDC 350 write(iout,'(//" (",i3," ) EIGENVECTORS FOR ATOMS ",I6," - ", C-RDC $ I6," AS CORE ATOMS ", C-RDC $ /" ARE OBTAINED FROM THE EIGENVECTORS OF SUBSYSTEMS")') C-RDC $ id, iprvec(id1),iprvec(id2) C-RDC write(iout,'(8(1x,i6))') (itmpc(i),i=1,nc) C-RDC write(iout,'(" EIGENVECTORS FOR ATOMS ",I6," - ", C-RDC $ I6," AS CORE OR INNER BUFFER ATOMS ", C-RDC $ /" ARE OBTAINED FROM THE EIGENVECTORS OF THE ", C-RDC $ "SUBSYSTEMS ABOVE AND SUBSYSTEMS")') C-RDC $ iprvec(id1),iprvec(id2) if (nb1.eq.0) then C-RDC write(iout,'(" ---")') else C-RDC write(iout,'(8(1x,i6))') (itmpb1(i),i=1,nb1) endif C-RDC write(iout,'(" EIGENVECTORS FOR ATOMS ",I6," - ", C-RDC $ I6," AS CORE, INNER OR OUTER BUFFER ATOMS ", C-RDC $ /" ARE OBTAINED FROM THE EIGENVECTORS OF ALL ", C-RDC $ "THESE SUBSYSTEM, PLUS",/" SUBSYSTEMS")') C-RDC $ iprvec(id1),iprvec(id2) if (nb2.eq.0) then C-RDC write(iout,'(" ---")') else C-RDC write(iout,'(8(1x,i6))') (itmpb2(i),i=1,nb2) endif C . now copy these to the list of subsystems that need to be printed. C . make sure the list is unique and find the right energy interval. do i=1,nb1 nc = nc + 1 itmpc(nc) = itmpb1(i) enddo do i=1,nb2 nc = nc + 1 itmpc(nc) = itmpb2(i) enddo call bsort1(nc,itmpc) if (nprsub.eq.0) goto 450 kpr = 1 kc = 1 410 if (kc.le.nc) then 420 if (iprsub(kpr).lt.itmpc(kc)) then kpr = kpr + 1 if (kpr.gt.nprsub) goto 450 goto 420 elseif (iprsub(kpr).gt.itmpc(kc)) then kc = kc +1 goto 410 else C . find the greater energy interval kpr2 = 2*kpr kpr1 = kpr2-1 if (efprvec(id)) then C . find the Fermi energy for these subsystem do icr=1,ncores i1 = isubend(icr-1)+1 i2 = isubend(icr) if ((iprsub(kpr).ge.i1) $ .and.(iprsub(kpr).le.i2)) then eprsub(kpr1) = min(eprsub(kpr1), $ efermi(icr)-eprvec(id1)) eprsub(kpr2) = max(eprsub(kpr2), $ efermi(icr)+eprvec(id1)) goto 430 endif enddo else eprsub(kpr1) = min(eprsub(kpr1),eprvec(id1)) eprsub(kpr2) = max(eprsub(kpr2),eprvec(id2)) endif C . make the subsystem list unique 430 do i=kc,nc-1 itmpc(i) = itmpc(i+1) enddo nc = nc - 1 goto 410 endif endif 450 do 460 i=1,nc nprsub = nprsub + 1 nprsub2 = 2*nprsub nprsub1 = nprsub2-1 iprsub(nprsub) = itmpc(i) if (efprvec(id)) then C . find the Fermi energy for these subsystem do icr=1,ncores i1 = isubend(icr-1)+1 i2 = isubend(icr) if ((iprsub(nprsub).ge.i1) $ .and.(iprsub(nprsub).le.i2)) then eprsub(nprsub1) = min(eprsub(nprsub1), $ efermi(icr)-eprvec(id1)) eprsub(nprsub2) = max(eprsub(nprsub2), $ efermi(icr)+eprvec(id1)) goto 460 endif enddo else eprsub(nprsub1) = eprvec(id1) eprsub(nprsub2) = eprvec(id2) endif 460 continue call bsort2(nprsub,iprsub,eprsub) enddo if (nprvec.gt.1) then C-RDC write(iout,'(/" CONTRACTED PRTVEC PARAMETERS TO SAVE ON ", C-RDC $ "SCF ITERATIONS.")') endif end