C C $Id: diat.F,v 1.3 1998/07/16 16:39:36 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE 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, . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC,resinter) c . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC) c resinter = 1 if interresidue, 0 if intraresidue interaction between I and J C C ILEVEL CODES: 0 ==> DO ENNIJ, HIJ, REPIJ, EENUC C 1 ==> DO ENNIJ, REPIJ, EENUC C 2 ==> DO ENNIJ, HIJ, REPIJ C 3 ==> DO ENNIJ, REPIJ C 4 ==> DO REPIJ C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" DIMENSION AGI(4),BGI(4),CGI(4),AGJ(4),BGJ(4),CGJ(4), . HIJ(9,9),REPIJ(10,10),EENUC(2,10) integer resinter logical lok C C LOCAL ARRAYS: C LOGICAL ORBSI,ORBSJ DIMENSION SIJ(9,9),IREPUL(10) DATA IREPUL /1,2,5,3,6,8,4,7,9,10/ SAVE IREPUL C C do i=1,10 do j=1,10 REPIJ(j,i) = 0.0d0 enddo enddo do i=1,9 do j=1,9 HIJ(j,i) = 0.0d0 SIJ(j,i) = 0.0d0 enddo enddo ENNIJ = 0.0d0 ORBSI = NORBSI.GT.0 ORBSJ = NORBSJ.GT.0 C C GET DIATOMIC OVERLAPS AND ASSIGN OFF-DIAGONAL BLOCK OF C 1-ELECTRON MATRIX (NON-SPARKLE ATOMS ONLY). C IF(ORBSI.AND.ORBSJ.AND.(ILEVEL.EQ.0.OR.ILEVEL.EQ.2))THEN CALL OVERLP(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,XI,YI,ZI, . IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,XJ,YJ,ZJ, . RIJ,SIJ) BETAI = BSI DO 20 I=1,NORBSI IF(I.EQ.2)THEN BETAI = BPI ELSEIF(I.EQ.5)THEN BETAI = BDI ENDIF BETAJ = BSJ DO 10 J=1,NORBSJ IF(J.EQ.2)THEN BETAJ = BPJ ELSEIF(J.EQ.5)THEN BETAJ = BDJ ENDIF HIJ(I,J) = 0.5D0*(BETAI+BETAJ)*SIJ(I,J) 10 CONTINUE 20 CONTINUE ENDIF C C GET ELECTRON-ELECTRON REPULSIONS. TREAT SPARKLES AS MONOPOLES. C NIMAX = MAX(NORBSI,1) NJMAX = MAX(NORBSJ,1) CALL REPUL(IAI,NIMAX,AI0,AI1,AI2,DI1,DI2,XI,YI,ZI, . IAJ,NJMAX,AJ0,AJ1,AJ2,DJ1,DJ2,XJ,YJ,ZJ, . RIJ,REPIJ) C C WE ARE DONE IF ONLY REPIJ HAS BEEN REQUESTED. C IF(ILEVEL.EQ.4) RETURN C C ASSEMBLE ELECTRON-NUCLEAR ATTRACTION INTEGRALS. C C ELECTRONS ON I ATTRACTED TO CORE OF J: C IF(ILEVEL.EQ.0.OR.ILEVEL.EQ.1)THEN NI = (NORBSI*(NORBSI+1))/2 IF(NI.GT.0)THEN DO 130 K=1,NI EENUC(1,K) = -ZNUCJ*REPIJ(IREPUL(K),1) 130 CONTINUE ENDIF C NJ = (NORBSJ*(NORBSJ+1))/2 C C ELECTRONS ON J ATTRACTED TO CORE OF I: C IF(NJ.GT.0)THEN DO 140 K=1,NJ EENUC(2,K) = -ZNUCI*REPIJ(1,IREPUL(K)) 140 CONTINUE ENDIF ENDIF C C CORE-CORE REPULSION: C IAMAX = MAX(IAI,IAJ) IAMIN = MIN(IAI,IAJ) ZIJ = ZNUCI*ZNUCJ REPSS = REPIJ(1,1) EIJ = ZIJ*REPSS ARGI = ACI*RIJ IF(ARGI.GT.25.0D0)THEN EXPI = 0.0D0 ELSE EXPI = EXP(-ARGI) ENDIF ARGJ = ACJ*RIJ IF(ARGJ.GT.25.0D0)THEN EXPJ = 0.0D0 ELSE EXPJ = EXP(-ARGJ) ENDIF IF(IAMIN.EQ.1.AND.(IAMAX.EQ.7.OR.IAMAX.EQ.8))THEN C C N-H OR O-H INTERACTIONS: C IF(IAI.EQ.1)THEN HTERM = EXPI XTERM = EXPJ ELSE HTERM = EXPJ XTERM = EXPI ENDIF ENNIJ = EIJ*(1.0D0 + RIJ*XTERM + HTERM) ELSE C C REGULAR INTERACTION: C ENNIJ = EIJ*(1.0D0 + EXPI + EXPJ) ENDIF C C ADD GAUSSIANS FOR AM1 OR PM3 (NON-SPARKLE ATOMS ONLY). C IF(AM1)THEN GTERM = 0.0D0 IF(ORBSI)THEN DO 150 IGAUS=1,4 ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2 IF(ARGI.LT.25.0D0)THEN GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI) ENDIF 150 CONTINUE ENDIF IF(ORBSJ)THEN DO 160 JGAUS=1,4 ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2 IF(ARGJ.LT.25.0D0)THEN GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ) ENDIF 160 CONTINUE ENDIF ENNIJ = ENNIJ + GTERM*ZIJ/RIJ ELSEIF(RM1)THEN GTERM = 0.0D0 IF(ORBSI)THEN DO 151 IGAUS=1,4 ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2 IF(ARGI.LT.25.0D0)THEN GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI) ENDIF 151 CONTINUE ENDIF IF(ORBSJ)THEN DO 161 JGAUS=1,4 ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2 IF(ARGJ.LT.25.0D0)THEN GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ) ENDIF 161 CONTINUE ENDIF ENNIJ = ENNIJ + GTERM*ZIJ/RIJ ELSEIF(AM1D)THEN GTERM = 0.0D0 IF(ORBSI)THEN DO 152 IGAUS=1,4 ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2 IF(ARGI.LT.25.0D0)THEN GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI) ENDIF 152 CONTINUE ENDIF IF(ORBSJ)THEN DO 162 JGAUS=1,4 ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2 IF(ARGJ.LT.25.0D0)THEN GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ) ENDIF 162 CONTINUE ENDIF ENNIJ = ENNIJ + GTERM*ZIJ/RIJ ELSEIF(PM3PDDG)THEN GTERM = 0.0D0 IF(ORBSI)THEN DO IGAUS=1,4 ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2 IF(ARGI.LT.25.0D0)THEN GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI) ENDIF END DO ENDIF IF(ORBSJ)THEN DO JGAUS=1,4 ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2 IF(ARGJ.LT.25.0D0)THEN GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ) ENDIF END DO ENDIF ENNIJ = ENNIJ + GTERM*ZIJ/RIJ DNAPNB = 1.0d0/(ZNUCI+ZNUCJ) DO IGAUS = 1,2 DO JGAUS = 1,2 TERM1 = ZNUCI*PA7(IGAUS,IAI) + ZNUCJ*PA7(JGAUS,IAJ) TERM2 = RIJ - DA7(IGAUS,IAI) - DA7(JGAUS,IAJ) TERM3 = DNAPNB*TERM1*EXP(-10.0d0*TERM2*TERM2) ENNIJ = ENNIJ + TERM3 END DO END DO ELSEIF(PM3)THEN GTERM = 0.0D0 IF(ORBSI)THEN DO 180 IGAUS=1,2 ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2 IF(ARGI.LT.25.0D0)THEN GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI) ENDIF 180 CONTINUE ENDIF IF(ORBSJ)THEN DO 200 JGAUS=1,2 ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2 IF(ARGJ.LT.25.0D0)THEN GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ) ENDIF 200 CONTINUE ENDIF c c *** Modification : GM - 01 fev 2000 c *** Change : divcon.1.0.C011 c *** Object : introduce PIF c (PM3 Parametrized Intermolecular Function for water) c IF (PIF) then ! intermolecular bond ! use PIF term ! three possibilities: ! 1) H2O-H2O -> gpif (interres = 1) ! 2) solute-H2O -> gwalid (interres = 2) ! 3) solute-solute' -> gwalid (except no C-C parameters) ! -> not implemented if (resinter.eq.0) then ! intramolecular bond ! use standard PM3 term ! write(*,'(a,2i5)') "intramolecular-pm3", interres, resinter ENNIJ = ENNIJ + GTERM*ZIJ/RIJ elseif (resinter.eq.2) then ! intramolecular bond ! solute-water ! write(*,'(a,2i5)') "solute-water-gwalid", interres, resinter call gwalid(iai,iaj,rij,ewalid,dewalid) ENNIJ = ENNIJ + ewalid elseif (resinter.eq.1) then ! intramolecular bond ! water-water ! write(*,'(a,2i5)') "water-water-gpif", interres, resinter call gpif(iai,iaj,rij,epif,depif,lok) if (lok) then ENNIJ = ENNIJ + epif else write(*,'("WARNING: Use standard interaction between",i5, & " and",i5)') iai, iaj ENNIJ = ENNIJ + GTERM*ZIJ/RIJ endif else continue endif ELSEIF (MAIS) THEN call gmais(iai,iaj,rij,emais,demais) ENNIJ = ENNIJ + emais ELSE ENNIJ = ENNIJ + GTERM*ZIJ/RIJ ENDIF ENDIF RETURN END