C C $Id: mosub.F,v 1.11 1998/07/16 16:40:06 jjv5 Exp $ C C------------------------------------------------------------------------ SUBROUTINE MOSUB(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,pvecnow #ifdef MPI #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 BOLZMANN CONSTANT IN eV/Kelvin: C DATA BOLTZ /8.617335408D-5/ SAVE BOLTZ C BETAH = 1.0D0/(BOLTZ*TEMPK) nfull = 0 nfroz = 0 iflag = 0 ! TODO: check and remove need for iflag in this subroutine C IIMAX = IIMAT(NATOMS+1)-1 C #ifdef MPI 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.0d0 EVECSQ(I) = 0.0d0 ENDDO if ((iter.eq.1.) .and. (iflag.eq.0)) then DO I=1,IORBMX FERMI(I) = 0.0d0 ENDDO endif 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 DO 1000 IPASS=1,NPASS 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 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 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)+FPMEDIAG(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)+FPMEDIAT(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 NOCC = NELEC/2 K1 = IORBPT(K) ! if (wrtscr) then ! write(iscr,'(" full ",i5,$)') k ! 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 WRITE(iout,'(/" ERROR IN DIAG -- NO CONVERGENCE IN", . " QR ITERATION")') RETURN 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 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 #ifdef MPI 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_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IORBMX fermi(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(eval,TMP_MPI,IORBMX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IORBMX eval(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(evecsq,TMP_MPI,IORBMX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IORBMX evecsq(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_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) ! do i = 1,NSUB ! time_subs(i) = TMP_MPI(i) ! enddo ! call balance_bytime ! endif #endif IF(IPASS.EQ.1.AND..NOT.STAND .and. (iflag .eq. 0))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 WRITE(0,'(/" ERROR IN DOFERM -- NO CONVERGENCE", . " IN FERMI ENERGY DETERMINATION")') ENDIF ENDIF 1000 CONTINUE #ifdef MPI 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_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,NUM_PDIAG pdiag(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(PDIAT,TMP_MPI,NUM_PDIAT, + MPI_DOUBLE_PRECISION,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. RETURN END