C C $Id: fock.F,v 1.7 1998/07/16 16:39:47 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE FOCK(ITER,ECORE1) C C BUG FIX: 7/26/99 -- S. DIXON. WHEN SPARKLES APPEAR AFTER THE FIRST C ATOM, THEN THE IJBOND POINTER CAN MISIDENTIFY ENTRIES IN THE BONDED C PAIRLIST. THIS WAS FIXED BY MAKING SURE THAT ATOMS I AND J ARE BOTH C VALID CANDIDATES FOR THE PAIRLIST BEFORE USING THE IJBOND POINTER. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" #ifdef MPI #include "mpif.h" integer ijk #endif C C LOCAL: C LOGICAL ORBSI,ORBSJ,ORBSIJ,ITER1,SHARE1,SHARE2 DIMENSION AGI(4),BGI(4),CGI(4),AGJ(4),BGJ(4),CGJ(4),HIJ(9,9), . REPIJ(10,10),EENUC(2,10),PII(9,9),PJJ(9,9),PIJ(9,9) real*8 DXREPIJ(10,10),DYREPIJ(10,10),DZREPIJ(10,10) real*8 DXHIJ(9,9),DYHIJ(9,9),DZHIJ(9,9) integer resinter DIMENSION IREPUL(4,4),LQUANT(9) DATA IREPUL /1,2,3,4, . 2,5,6,7, . 3,6,8,9, . 4,7,9,10/ DATA LQUANT /0,1,1,1,2,2,2,2,2/ SAVE IREPUL,LQUANT C C C INITIALIZE FOCK MATRIX. C IIMAX = IIMAT(NATOMS+1)-1 IPMAX = IP1(NATOMS+1) IJMAX = IJMAT(IPMAX)-1 DO 10 II=1,IIMAX FDIAG(II) = 0.0D0 10 CONTINUE DO 20 IJ=1,IJMAX FDIAT(IJ) = 0.0D0 20 CONTINUE C C INITIALIZE 1-ELECTRON MATRIX IF THIS IS THE FIRST SCF ITERATION. C ITER1 = ITER.EQ.1 IF(ITER1)THEN ecore1 = 0.0 C JV First init everything to zero, then go back and do diagonal of hdiag C Can leave this in serial version, but done this way for parallel version do II=1,IIMAX HDIAG(II) = 0.0D0 enddo do IJ=1,IJMAX HDIAT(IJ) = 0.0D0 enddo C JV Must set whole arrays to zero before doing any of this in parallel C because after first cycle arrays will have garabage all over the place C Want everything zero, but diagonal done only by pieces C Can't have all PE's sum filled in diagonal again #ifdef MPI DO 50 I=myid+1,NATOMS,nproc #else DO 50 I=1,NATOMS #endif IAI = IATNUM(I) NORBSI = NATORB(IAI) C C SKIP IF IT'S NOT A REAL ATOM. C IF(NORBSI.EQ.0) GO TO 50 IJ = IIMAT(I) DO 40 II=1,NORBSI LI = LQUANT(II) DO 30 JJ=1,II IF(JJ.EQ.II) HDIAG(IJ) = UCORE(LI,IAI) IJ = IJ + 1 30 CONTINUE 40 CONTINUE 50 CONTINUE ENDIF C================ SOLVATION CALCULATION =================== C V.G. updated 10/12/98 C CALCULATE 1 ELECTRON INTEGRAL MATRIX ELEMENTS DUE TO THE C INTERACTION BETWEEN THE SURFACE CHARGES, INDUCED BY C SOLVENT REACTION FIELD, AND SOLUTE CORES and ELECTRONS C #ifdef SCRF_IS_ON IF(SOLVCALC) THEN call pbfock(ierror,ecore1) if(ierror.gt.0) then c write(iout,'(" ERROR IN SUBROUTINE ")') RETURN endif ENDIF #endif C ============= END of SOLVATION CALCULATION =============== C C FOCK MATRIX CONTRIBUTIONS FROM 1-CENTER, 2-ELECTRON TERMS. C C jv Would it be possible to split this up by atoms C jv then just do a global sum on fdiag, like DENSUB ? #ifdef MPI DO 80 I=myid+1,NATOMS,nproc #else DO 80 I=1,NATOMS #endif IAI = IATNUM(I) NORBSI = NATORB(IAI) IF(NORBSI.EQ.0) GO TO 80 GSSI = GSS(IAI) ISS = IIMAT(I) PSS = PDIAG(ISS) FDIAG(ISS) = FDIAG(ISS) + 0.5D0*PSS*GSSI IF(NORBSI.EQ.1) GO TO 80 GPPI = GPP(IAI) GSPI = GSP(IAI) GP2I = GP2(IAI) HSPI = HSP(IAI) HGPLUS = HSPI + GSPI GGPLUS = GPPI + GP2I GGMINU = GPPI - GP2I C C INDEXING SCHEME FOR LOWER TRIANGLE OF DIAGONAL BLOCK: C C IIMAT(I) ISS C IIMAT(I)+1 IIMAT(I)+2 --\ IPXS IPXPX C IIMAT(I)+3 IIMAT(I)+4 IIMAT(I)+5 --/ IPYS IPYPX IPYPY C IIMAT(I)+6 IIMAT(I)+7 IIMAT(I)+8 IIMAT(I)+9 IPZS IPZPX IPZPY IPZPZ C IPXS = ISS + 1 PPXS = PDIAG(IPXS) IPXPX = IPXS + 1 PPXPX = PDIAG(IPXPX) IPYS = IPXPX + 1 PPYS = PDIAG(IPYS) IPYPX = IPYS + 1 PPYPX = PDIAG(IPYPX) IPYPY = IPYPX + 1 PPYPY = PDIAG(IPYPY) IPZS = IPYPY + 1 PPZS = PDIAG(IPZS) IPZPX = IPZS + 1 PPZPX = PDIAG(IPZPX) IPZPY = IPZPX + 1 PPZPY = PDIAG(IPZPY) IPZPZ = IPZPY + 1 PPZPZ = PDIAG(IPZPZ) C PXY = PDIAG(IPXPX) + PDIAG(IPYPY) PXZ = PDIAG(IPXPX) + PDIAG(IPZPZ) PYZ = PDIAG(IPYPY) + PDIAG(IPZPZ) PXYZ = PXY + PDIAG(IPZPZ) FDIAG(ISS) = FDIAG(ISS) + PXYZ*(GSPI - 0.5D0*HSPI) C HTERM = 2.0D0*HSPI - 0.5D0*HGPLUS FDIAG(IPXS) = FDIAG(IPXS) + PPXS*HTERM FDIAG(IPYS) = FDIAG(IPYS) + PPYS*HTERM FDIAG(IPZS) = FDIAG(IPZS) + PPZS*HTERM C PSSGH = PSS*(GSPI - 0.5D0*HSPI) GHALF = GPPI*0.5D0 GTERM = GP2I - 0.25D0*GGMINU FDIAG(IPXPX) = FDIAG(IPXPX) + PSSGH + PPXPX*GHALF + PYZ*GTERM FDIAG(IPYPY) = FDIAG(IPYPY) + PSSGH + PPYPY*GHALF + PXZ*GTERM FDIAG(IPZPZ) = FDIAG(IPZPZ) + PSSGH + PPZPZ*GHALF + PXY*GTERM C GTERM = GGMINU - 0.25D0*GGPLUS FDIAG(IPYPX) = FDIAG(IPYPX) + PPYPX*GTERM FDIAG(IPZPX) = FDIAG(IPZPX) + PPZPX*GTERM FDIAG(IPZPY) = FDIAG(IPZPY) + PPZPY*GTERM 80 CONTINUE C C OFF-DIAGONAL BLOCKS AND DIATOMIC CONTRIBUTIONS TO DIAGONAL BLOCKS: C IF(NATOMS.EQ.1) RETURN C IJFULL = 1 IJBOND = 1 C C IJFULL IS A POINTER TO POSITION IN LIST OF ALL POSSIBLE ATOM PAIRS. C IJBOND IS A POINTER TO POSTIION IN BONDED ATOM PAIRLIST IPAIR. C c jv Think about doing just atoms of the subsystem here - but maybe can't c if they require interaction with all other atoms #ifdef MPI DO 1000 I=myid+2,NATOMS,nproc C JV Make IJBOND start at the beginning of the pairlist for this atom IJBOND = IP1(I) #else DO 1000 I=2,NATOMS #endif II0 = IIMAT(I) IAI = IATNUM(I) IF(IAI.EQ.0) GO TO 1000 NORBSI = NATORB(IAI) ORBSI = NORBSI.GT.0 NREPI = MAX((NORBSI*(NORBSI+1))/2,1) AI0 = AL(0,IAI) AI1 = AL(1,IAI) AI2 = AL(2,IAI) DI1 = DL(1,IAI) DI2 = DL(2,IAI) BSI = BETA(0,IAI) BPI = BETA(1,IAI) BDI = BETA(2,IAI) IF(ITER1)THEN NPQI = NQUANT(IAI) XISI = EXPNT(0,IAI) XIPI = EXPNT(1,IAI) XIDI = EXPNT(2,IAI) ACI = ACORE(IAI) DO 90 IGAUS=1,4 AGI(IGAUS) = AGAUS(IGAUS,IAI) BGI(IGAUS) = BGAUS(IGAUS,IAI) CGI(IGAUS) = CGAUS(IGAUS,IAI) 90 CONTINUE ENDIF ZNUCI = ZCHG(IAI) C C TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM I. C IF(ORBSI)THEN IJ = IIMAT(I) DO 110 IP=1,NORBSI DO 100 JP=1,IP PII(IP,JP) = PDIAG(IJ) PII(JP,IP) = PDIAG(IJ) IJ = IJ + 1 100 CONTINUE 110 CONTINUE ENDIF XI = XYZ(1,I) YI = XYZ(2,I) ZI = XYZ(3,I) C C GET I-J INTERSECTION IN FOCK MATRIX. THE DIATOMIC BLOCK FDIAT C WILL ONLY BE FILLED IF I AND J ARE A "BONDED" PAIR STORED C IN IPAIR. HOWEVER I AND J ALWAYS AFFECT THE DIAGONAL BLOCK C FDIAG OF THE OTHER ATOM. C c#ifdef MPI_NOMC c DO 500 J=myid+1,I-1,nproc cC JV Must set IJBOND to start of pairtlist for atom I in outer loop cC before using with increment in 500 loop c#else DO 500 J=1,I-1 c#endif IAJ = IATNUM(J) IF(IAJ.EQ.0) GO TO 500 resinter = interres(i,j) NORBSJ = NATORB(IAJ) ORBSJ = NORBSJ.GT.0 ORBSIJ = ORBSI.AND.ORBSJ C C SEE IF (I,J) IS THE NEXT PAIR IN THE BONDED PAIRLIST AND C ASSIGN "SHARE" STATUS. C C ********************************************************************** C BUG FIX - AN ADDITIONAL CONDITION WAS ADDED: C C (PREVIOUS) (FIXED) C IF(IPAIR(IJBOND).EQ.J) --> IF(ORBSIJ.AND.IPAIR(IJBOND).EQ.J) C C WHEN ATOM I IS A SPARKLE, THEN (I,J) IS NOT A VALID CANDIDATE C FOR THE BONDED PAIRLIST, SO THE IPAIR(IJBOND).EQ.J TEST IS C NOT VALID. THE ORBSIJ CONDITION WAS ADDED TO MAKE SURE A C COINCIDENTAL MATCH DOESN'T THROW OFF THE SHARE1/SHARE2 STATUS C AND THE FOCK INDEXING SCHEME. C ********************************************************************** C IF(ORBSIJ.AND.IPAIR(IJBOND).EQ.J)THEN SHARE1 = NSHARE(1,IJBOND).NE.0 SHARE2 = .TRUE. ELSE SHARE1 = .FALSE. SHARE2 = .FALSE. ENDIF C C SHARE1 INDICATES WHETHER OR NOT ATOMS I AND J EVER APPEAR IN THE C SAME SUBSYSTEM AS NON-BUFFER ATOMS. IF NOT, THEN THE I-J BLOCK C OF PDIAT WILL BE ZERO. THIS INFORMATION MAY BE USED TO SKIP OVER C THE 2-ELECTRON PORTION OF FDIAT. C C SHARE2 INDICATES WHETHER THE ATOMS EVER APPEAR WITH ANY STATUS C IN THE SAME SUBSYSTEMS. IF NOT, THEN THE I-J BLOCK OF HDIAT C WILL BE IGNORED AND THERE IS NO NEED TO COMPUTE OVERLAP INTEGRALS. C C NOTE THAT BOTH SHARE1 AND SHARE2 ARE AFFECTED IF THE CUTBOND C OPTION IS USED. IN THAT CASE, ATOMS I AND J MUST BE IN THE C SAME SUBSYSTEM AND SEPARATED BY NO MORE THAN THE CUTOFF TO C BE CONSIDERED. C NREPJ = MAX((NORBSJ*(NORBSJ+1))/2,1) AJ0 = AL(0,IAJ) AJ1 = AL(1,IAJ) AJ2 = AL(2,IAJ) DJ1 = DL(1,IAJ) DJ2 = DL(2,IAJ) BSJ = BETA(0,IAJ) BPJ = BETA(1,IAJ) BDJ = BETA(2,IAJ) IF(ITER1)THEN NPQJ = NQUANT(IAJ) XISJ = EXPNT(0,IAJ) XIPJ = EXPNT(1,IAJ) XIDJ = EXPNT(2,IAJ) ACJ = ACORE(IAJ) DO 120 IGAUS=1,4 AGJ(IGAUS) = AGAUS(IGAUS,IAJ) BGJ(IGAUS) = BGAUS(IGAUS,IAJ) CGJ(IGAUS) = CGAUS(IGAUS,IAJ) 120 CONTINUE ENDIF IF(PBC)THEN CALL PBCXYZ(I,J,XJ,YJ,ZJ) ELSE XJ = XYZ(1,J) YJ = XYZ(2,J) ZJ = XYZ(3,J) ENDIF ZNUCJ = ZCHG(IAJ) C C TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM J. C IF(ORBSJ)THEN IJ = IIMAT(J) DO 140 IP=1,NORBSJ DO 130 JP=1,IP PJJ(IP,JP) = PDIAG(IJ) PJJ(JP,IP) = PDIAG(IJ) IJ = IJ + 1 130 CONTINUE 140 CONTINUE ENDIF C C TEMPORARILY STORE OFF-DIAGONAL BLOCK OF DENSITY MATRIX FOR I-J C INTERSECTION. C IF(ORBSIJ.AND.SHARE1)THEN IJ = IJMAT(IJBOND) DO 170 IP=1,NORBSI DO 160 JP=1,NORBSJ PIJ(IP,JP) = PDIAT(IJ) IJ = IJ + 1 160 CONTINUE 170 CONTINUE ENDIF C C INTEGRAL COMPUTATION C -------------------- C C FIRST SCF ITERATION: COMPUTE 2-ELECTRON INTEGRALS REPIJ, C 1-ELECTRON INTEGRALS EENUC, AND C CORE-CORE REPULSIONS ENNIJ. C COMPUTE DIATOMIC 1-ELECTRON INTEGRALS C HIJ ONLY IF ATOMS I AND J SHARE A C SUBSYSTEM (ANY STATUS). C C C SUBSEQUENT ITERATIONS: IF DIRECT CALCULATION, COMPUTE ONLY REPIJ. C IF NOT DIRECT, THEN COMPUTE NO INTEGRALS. C IF(ITER1.OR.DIRECT)THEN IF(.NOT.DIRECT)THEN RIJ = RPAIR(IJFULL) ELSE RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2) ENDIF ILEVEL = 0 IF(.NOT.ITER1)THEN C C DIRECT CALCULATION, AND NOT FIRST ITERATION. GET ONLY C TWO-ELECTRON INTEGRALS. C ILEVEL = 4 ELSEIF(.NOT.SHARE2)THEN C C FIRST ITERATION, BUT NOT BONDED. NEED ONLY CORE-ELECTRON C INTEGRALS EENUC AND CORE-CORE REPULSION ENNIJ. C ILEVEL = 1 ENDIF CALL DIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1, . DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI, . IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1, . DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ, c . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,DXREPIJ,DYREPIJ,DZREPIJ, c . DXHIJ,DYHIJ,DZHIJ, c . DXENNIJ,DYENNIJ,DZENNIJ, c . EENUC,resinter) . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC,resinter) c . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC) C C STORE REPULSION INTEGRALS IN GLOBAL ARRAY IF IT'S NOT A C DIRECT CALCULATION. C IF(.NOT.DIRECT)THEN IJR = IJREP(IJFULL) DO 220 IR=1,NREPI DO 210 JR=1,NREPJ EEREP(IJR) = REPIJ(IR,JR) IJR = IJR + 1 210 CONTINUE 220 CONTINUE ENDIF C C ON FIRST SCF ITERATION, STORE 1-ELECTRON MATRIX AND C CORE-CORE REPULSIONS. C IF(ITER1)THEN C C 1-ELECTRON (DIAGONAL, I): C IF(ORBSI)THEN IJ = IIMAT(I) DO 230 II=1,NREPI HDIAG(IJ) = HDIAG(IJ) + EENUC(1,II) IJ = IJ + 1 230 CONTINUE ENDIF C C 1-ELECTRON (DIAGONAL, J): C IF(ORBSJ)THEN IJ = IIMAT(J) DO 240 JJ=1,NREPJ HDIAG(IJ) = HDIAG(IJ) + EENUC(2,JJ) IJ = IJ + 1 240 CONTINUE ENDIF C C 1-ELECTRON (OFF-DIAGONAL, IJ): C IF(ORBSIJ.AND.SHARE2)THEN IJ = IJMAT(IJBOND) DO 260 II=1,NORBSI DO 250 JJ=1,NORBSJ HDIAT(IJ) = HIJ(II,JJ) IJ = IJ + 1 250 CONTINUE 260 CONTINUE ENDIF C C CORE-CORE REPULSIONS: C IF(.NOT.DIRECT) ENUCLR(IJFULL) = ENNIJ ECORE1 = ECORE1 + ENNIJ ENDIF ELSE C C IT'S NOT THE FIRST SCF ITERATION AND IT'S NOT A DIRECT C CALCULATION. EXTRACT REPULSION INTEGRALS FROM GLOBAL C ARRAY. C IJR = IJREP(IJFULL) DO 310 IR=1,NREPI DO 300 JR=1,NREPJ REPIJ(IR,JR) = EEREP(IJR) IJR = IJR + 1 300 CONTINUE 310 CONTINUE ENDIF C C DETERMINE OFF-DIAGONAL BLOCK OF FOCK MATRIX FOR I-J INTERSECTION. C BLOCK WILL BE ZERO, HOWEVER, IF ATOMS I AND J DO NOT SHARE A C SUBSYSTEM AS NON-BUFFER ATOMS. C IF(ORBSIJ.AND.SHARE1)THEN IJ = IJMAT(IJBOND) DO 350 IG=1,NORBSI DO 340 JG=1,NORBSJ GIJ = 0.0D0 DO 330 IP=1,NORBSI IR = IREPUL(IG,IP) DO 320 JP=1,NORBSJ JR = IREPUL(JG,JP) GIJ = GIJ + PIJ(IP,JP)*REPIJ(IR,JR) 320 CONTINUE 330 CONTINUE FDIAT(IJ) = FDIAT(IJ) - 0.5D0*GIJ IJ = IJ + 1 340 CONTINUE 350 CONTINUE ENDIF C C ADD 2-CENTER 2-ELECTRON TERMS ARISING FROM ATOM J TO DIAGONAL C BLOCK OF ATOM I. C IJ = II0 DO 400 IG=1,NORBSI DO 390 JG=1,IG JJG = IREPUL(JG,JG) IJG = IREPUL(IG,JG) GIJ = 0.0D0 DO 380 IP=1,NORBSJ DO 370 JP=1,NORBSJ IJP = IREPUL(IP,JP) GIJ = GIJ + REPIJ(IJG,IJP)*PJJ(IP,JP) 370 CONTINUE 380 CONTINUE FDIAG(IJ) = FDIAG(IJ) + GIJ IJ = IJ + 1 390 CONTINUE 400 CONTINUE C C ADD 2-CENTER 2-ELECTRON TERMS ARISING FROM ATOM I TO DIAGONAL C OF ATOM J. C IJ = IIMAT(J) DO 480 IG=1,NORBSJ DO 470 JG=1,IG JJG = IREPUL(JG,JG) IJG = IREPUL(IG,JG) GIJ = 0.0D0 DO 460 IP=1,NORBSI DO 450 JP=1,NORBSI IJP = IREPUL(IP,JP) GIJ = GIJ + REPIJ(IJP,IJG)*PII(IP,JP) 450 CONTINUE 460 CONTINUE FDIAG(IJ) = FDIAG(IJ) + GIJ IJ = IJ + 1 470 CONTINUE 480 CONTINUE C C UPDATE POINTERS FOR FULL AND BOND-BASED PAIRLIST. C IJFULL = IJFULL + 1 IF(SHARE2) IJBOND = IJBOND + 1 500 CONTINUE 1000 CONTINUE #ifdef MPI C JV Sum FDIAG and FDIAT globally call mpi_allreduce(FDIAG,TMP_MPI,IIMAX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IIMAX fdiag(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(FDIAT,TMP_MPI,IJMAX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IJMAX fdiat(ijk) = TMP_MPI(ijk) enddo C JV Do this only if ITER1 also if(ITER1) then call mpi_allreduce(HDIAG,TMP_MPI,IIMAX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IIMAX hdiag(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(HDIAT,TMP_MPI,IJMAX, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,IJMAX hdiat(ijk) = TMP_MPI(ijk) enddo C JV Sum ECORE from all PE's - Only on first ITER (first subsys) call mpi_allreduce(ECORE1,TMP_MPI(1),1, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) ECORE1 = TMP_MPI(1) endif #endif C C C ADD IN 1-ELECTRON MATRIX. C DO 1010 II=1,IIMAX FDIAG(II) = FDIAG(II) + HDIAG(II) 1010 CONTINUE DO 1020 IJ=1,IJMAX FDIAT(IJ) = FDIAT(IJ) + HDIAT(IJ) 1020 CONTINUE RETURN END