C C $Id: mosub.F,v 1.11 1998/07/16 16:40:06 jjv5 Exp $ C C------------------------------------------------------------------------ SUBROUTINE MOSUBfdmx(GOODF,ITER,PSEUDO,pvecnow,IERROR) C C DETERMINES THE MOLECULAR ORBITALS FOR ALL SUBSYSTEMS. C C GOODF = FLAG FOR A GOOD ESTIMATE OF THE FERMI ENERGY. WHEN C THE DENSITY MATRIX HAS JUST BEEN INITIALIZED TO C DIAGONAL FORM, THEN WE DON'T HAVE A GOOD ESTIMATE C OF EFERMI. IN THIS CASE, TWO ROUNDS OF DIAGONALIZATIONS C WILL BE NECESSARY TO DETERMINE THE GLOBAL DENSITY C MATRIX. C C ITER = SCF ITERATION NUMBER: 1,2,... C C PSEUDO = LOGICAL FLAG DO CARRY OUT PSEUDO-DIAGONALIZATION C C RETURNED ERROR CODES: C C IERROR = 0 --> OKAY C IERROR = 1 --> (1) QR ITERATION IN DIAG DIDN'T CONVERGE, OR C (2) FERMI ENERGY IN DOFERM DIDN'T CONVERGE C C NOTE: THIS ROUTINE USES THE /WORK/ COMMON BLOCK FOR TEMPORARY C STORAGE. ANYTHING IN FF, EVEC, OR WW BELOW WILL BE DESTROYED. C C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" LOGICAL GOODF,PSEUDO, firstfdmx, pvecnow #ifdef MPI_IS_ON #include "mpif.h" integer ijk #endif #ifdef CRAY C JV Arrays for SSYEVX LAPACK call and timers integer IWORK(MSVAL),IFAIL(MSVAL),M double precision v1,v2,v3,v4,v5,v6 #endif C C LOGICAL FIRST DATA FIRST /.TRUE./ SAVE FIRST C C BOLTZMANN CONSTANT IN eV/Kelvin: C DATA BOLTZ /8.617335408D-5/ SAVE BOLTZ C BETAH = 1.0D0/(BOLTZ*TEMPK) ifdmx = iunit(28) C C #ifdef MPI_NOMC C JV Init eval in case we changed subsetting and they C JV have garbage from other subsystems - OK for serial version IORBMX = IORBPT(NSUB+1)-1 DO I=1,IORBMX EVAL(I) = 0.0 ENDDO CC JV Balance by number of orbitals per subsytem as initial guess call balance_bynorbs CCJV Init timers for next load balancing (by completion time) do i = 1,nsub time_subs(i) = 0.0 enddo #endif C WHEN NO RELIABLE ESTIMATE OF EFERMI EXISTS, THEN TWO DIAGONALIZATION C PASSES ARE NECESSARY. ON THE FIRST PASS, THE EIGENVECTORS ARE C USED TO GET THE QUANTITIES EVECSQ, WHICH ARE IN TURN USED TO C DETERMINE THE FERMI ENERGY. THEN THE EIGENVECTORS ARE DETERMINED C AGAIN, AND THE GLOBAL DENSITY MATRIX IS CONSTRUCTED. THIS APPROACH C AVOIDS THE NECESSITY OF STORING ALL THE SUBSYSTEM EIGENVECTORS C SIMULTANEOUSLY. UNDER NORMAL CIRCUMSTANCES, MULTIPLE PASSES ARE C ONLY NECESSARY ON THE FIRST CALL AT PROGRAM STARTUP. ON LATER C SCF ITERATIONS AND FOR LATER GEOMETRIES, THE FERMI ENERGY FROM THE C PREVIOUS CALL MAY BE USED. C IF(GOODF.OR.STAND.OR.NSUB.EQ.1)THEN NPASS = 1 ELSE NPASS = 2 ENDIF firstfdmx = .true. DO 1000 IPASS=1,NPASS nfull = 0 nfroz = 0 C C HAVE DIAG COMPUTE EIGENVALUES ONLY ON THE FIRST PASS. C DOEVAL = IPASS.EQ.1 C C ASSEMBLE FOCK MATRIX FOR EACH SUBSYSTEM BY ACCESSING THE CORRECT C BLOCKS OF THE GLOBAL FOCK MATRIX. C TOLERA = 1.0D-16 #ifdef MPI_NOMC DO 500 ijk= 1,my_numsubs if (iter .eq. 1) call etimer(watch(15)) k = my_subs(ijk) #else DO 500 K=1,NSUB #endif C C DIAGONAL BLOCKS FIRST. C if (first.or.notfroz(k)) then NORBSK = IORBPT(K+1)-IORBPT(K) IIORB = 0 IJLOC0 = 0 DO 100 I=IATOM1(K),IATOM1(K+1)-1 IATM = IATOMS(I) NORBSI = NATORB(IATNUM(IATM)) IJGLB = IIMAT(IATM) DO 20 IF=1,NORBSI IJLOC = IJLOC0 + IF + IIORB DO 10 JF=1,IF FF(IJLOC) = FDIAG(IJGLB) IJGLB = IJGLB + 1 IJLOC = IJLOC + NORBSK 10 CONTINUE 20 CONTINUE IIORB = IIORB + NORBSI IJLOC0 = IIORB*NORBSK 100 CONTINUE C C NOW OFF-DIAGONAL BLOCKS. C NPAIRS = IP1(NATOMS+1)-1 NATMS = IATOM1(K+1) - IATOM1(K) IF(NATMS.GT.1)THEN IATM0 = IATOM1(K) - 1 IATM1 = IATOMS(IATOM1(K)) II = NATORB(IATNUM(IATM1)) DO 200 I=2,NATMS IATM = IATOMS(IATM0+I) NORBSI = NATORB(IATNUM(IATM)) JJ = 0 IJLOC0 = 0 DO 180 J=1,I-1 JATM = IATOMS(IATM0+J) C C GET POSITION IN PAIRLIST OF (IATM,JATM) PAIR. C CALL IJFIND(NPAIRS,IATM,JATM,IJADDR) NORBSJ = NATORB(IATNUM(JATM)) C C COPY GLOBAL DIATOMIC BLOCK FDIAT TO FF, BUT ONLY IF IATM C AND JATM ARE CLOSE ENOUGH TO BOND. OTHERWISE, SET THE C BLOCK OF FF TO ZERO. C IF(IJADDR.NE.0)THEN IJGLB = IJMAT(IJADDR) DO 130 IF=1,NORBSI IJLOC = IJLOC0 + II + IF DO 120 JF=1,NORBSJ FF(IJLOC) = FDIAT(IJGLB) IJLOC = IJLOC + NORBSK IJGLB = IJGLB + 1 120 CONTINUE 130 CONTINUE ELSE DO 150 IF=1,NORBSI IJLOC = IJLOC0 + II + IF DO 140 JF=1,NORBSJ FF(IJLOC) = 0.0D0 IJLOC = IJLOC + NORBSK 140 CONTINUE 150 CONTINUE ENDIF JJ = JJ + NORBSJ IJLOC0 = JJ*NORBSK 180 CONTINUE II = II + NORBSI 200 CONTINUE ENDIF NK = NORBSK nk2 = nk*nk NOCC = NELEC/2 K1 = IORBPT(K) if (wrtscr) then C-RDC write(iscr,'(" full ",i5,$)') k C-RDC if(mod(k,5).eq.0.or.k.eq.nsub) write(iscr,'(6X)') endif nfull = nfull + 1 C . either a pseudo diagonalization or a full diagonalization IF(PSEUDO)THEN CALL DIAGP(NK,NOCC,FF,EVAL(K1),WW,VV,EVEC) ELSE #ifdef CRAY C JV Use DSYEV for normal machines, SSYEV for T3E call SSYEVX('V','A','L', NK,FF,NK,0,0,0,0,1.0D-7, + M,EVAL(K1),EVEC,NK,WW,MSORB2,IWORK,IFAIL,IERROR ) #else CALL DIAG(NK,FF,NK,TOLERA,WW,EVAL(K1), & IDEGEN,EVEC,IERROR) #endif IF(IERROR.NE.0)THEN C-RDC WRITE(IOUT,'(/" ERROR IN DIAG -- NO CONVERGENCE IN", C-RDC . " QR ITERATION")') RETURN ENDIF ENDIF else nfroz = nfroz + 1 endif C . see if we have a frozen subsystem here if ((ipass.eq.npass).and.(.not.notfroz(k))) then if (first) then C . save eigenvectors to file if (firstfdmx) then call opnfil(28,ierror) if (ierror.ne.0) return firstfdmx = .false. endif C-RDC write(ifdmx,'(3i8)') k, nk, nk2 C-RDC write(ifdmx,'(4E32.25)') (evec(iii),iii=1,nk2) else C . read eigenvectors from file if (firstfdmx) then call opnfil(28,ierror) if (ierror.ne.0) return firstfdmx = .false. endif read(ifdmx,'(3i8)') kk, nk, nk2 read(ifdmx,'(4E32.25)') (evec(iii),iii=1,nk2) norbsk = nk if (wrtscr) then C-RDC write(iscr,'(" froz ",i5,$)') k C-RDC if(mod(k,5).eq.0.or.k.eq.nsub) write(iscr,'(6X)') endif endif endif C C COMPUTE GLOBAL DENSITY MATRIX FOR STANDARD CALCULATION, C OR, IF MULTIPLE SUBSYSTEMS AND FINAL PASS, ADD SUBSYSTEM C CONTRIBUTION TO GLOBAL DENSITY MATRIX. IF, FOR SOME C STRANGE REASON, THERE IS ONLY ONE SUBSYSTEM, AND THE USER C HAS NOT SPECIFIED 'STANDARD', THEN THE DENSITY MATRIX C CANNOT BE DETERMINED UNTIL AFTER THE CALL TO DOFERM. C IF(STAND)THEN CALL DENFUL(NORBSK,EVEC) ELSEIF(IPASS.EQ.NPASS.AND.NSUB.GT.1)THEN C C IF AUTOMATIC SUBSETTING IS BEING USED, THEN THE POINTERS C IN THE ARRAY FERMI CAN CHANGE FROM ONE GEOMETRY TO THE C NEXT. TO BE SAFE, RECOMPUTE THE FERMI OCCUPATION NUMBERS C ON THE FIRST SCF ITERATION OF ALL GEOMETRIES AFTER THE C INITIAL ONE. C if (setch) then do i=1,ncores if (k.le.isubend(i)) then efermix = efermi(i) goto 201 endif enddo else efermix = efermi(1) endif 201 IF(AUTOSUB.AND..NOT.FIRST.AND.ITER.EQ.1)THEN DO 300 I=IORBPT(K),IORBPT(K+1)-1 DELTAE = EVAL(I) - EFERMIx ARG = BETAH*DELTAE IF(ARG.GT.50.0D0)THEN FERMI(I) = 0.0D0 ELSEIF(ARG.LT.-50.0D0)THEN FERMI(I) = 1.0D0 ELSE FERMI(I) = 1.0D0/(1.0D0 + EXP(ARG)) ENDIF 300 CONTINUE ENDIF CALL DENSUB(K,NORBSK,EVEC) ENDIF C C FOR THE FIRST PASS OF A D&C CALCULATION, FILL IN THIS C SUBSYSTEM'S CONTRIBUTION TO EVECSQ. WHEN THE WHOLE C EVECSQ ARRAY HAS BEEN ASSIGNED, THE FERMI ENERGY MAY BE C DETERMINED. C IF(IPASS.EQ.1.AND..NOT.STAND)THEN CALL ESQR(K,NORBSK,EVEC) C JV Moved call to DOFERM outside 500 loop C C IF THIS IS THE LAST SUBSYSTEM, THEN WE CAN DETERMINE C THE FERMI ENERGY. C ENDIF C C ADD TO GLOBAL DENSITY MATRIX FOR THE STRANGE SITUATION WHERE C 'STANDARD' HAS NOT BEEN DEFINED, BUT THERE IS ONLY ONE C SUBSYSTEM. C IF(.NOT.STAND.AND.NSUB.EQ.1.AND.IPASS.EQ.NPASS)THEN CALL DENSUB(K,NORBSK,EVEC) ENDIF if ((ipass.eq.npass).and.pvecnow) then call wrtvecsub(k) endif #ifdef MPI_NOMC C JV Time susbystems for balancing only on first iteration if (iter .eq. 1) then call etimer(watch(16)) time_subs(k) = watch(16)-watch(15) endif #endif 500 CONTINUE IF(IPASS.EQ.1.AND..NOT.STAND)THEN C JV Init all fermi values to zero in case we rearranged subsystems C JV OK to leave this for serial version IORBMX = IORBPT(NSUB+1)-1 do ijk = 1,IORBMX fermi(ijk) = 0.0 enddo CALL DOFERM(IERROR) IF(IERROR.NE.0)THEN C-RDC WRITE(IOUT,'(/" ERROR IN DOFERM -- NO CONVERGENCE", C-RDC . " IN FERMI ENERGY DETERMINATION")') ENDIF #ifdef MPI_NOMC C-RDC C JV Gather fermi and eval so we can rebalance and still keep all info if ( iter .eq. 1 ) then call mpi_allreduce(fermi,TMP_MPI,IORBMX, + MPI_REAL8,MPI_SUM,commsebomd,ier) do ijk = 1,IORBMX fermi(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(eval,TMP_MPI,IORBMX, + MPI_REAL8,MPI_SUM,commsebomd,ier) do ijk = 1,IORBMX eval(ijk) = TMP_MPI(ijk) enddo C C JV Load balance based on time once we've done a single pass call mpi_allreduce(time_subs,TMP_MPI,NSUB, + MPI_REAL8,MPI_SUM,commsebomd,ier) do i = 1,NSUB time_subs(i) = TMP_MPI(i) enddo C call balance_bytime endif #endif ENDIF 1000 CONTINUE #ifdef MPI_NOMC c JV Sum PDIAG and PDIAT from all PE's C NUM_PDIAG = IIMAT(NATOMS+1)-1 IPMAX = IP1(NATOMS+1) NUM_PDIAT = IJMAT(IPMAX)-1 CC call mpi_allreduce(PDIAG,TMP_MPI,NUM_PDIAG, + MPI_REAL8,MPI_SUM,commsebomd,ier) do ijk = 1,NUM_PDIAG pdiag(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(PDIAT,TMP_MPI,NUM_PDIAT, + MPI_REAL8,MPI_SUM,commsebomd,ier) do ijk = 1,NUM_PDIAT pdiat(ijk) = TMP_MPI(ijk) enddo C #endif C C CC C THE FIRST CALL TO MOSUB IS OVER: C FIRST = .FALSE. close(ifdmx) RETURN END