C C $Id: gcart.F,v 1.8 1998/07/16 16:39:50 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE GCART(GXYZ) C C SUBROUTINE TO COMPUTE CARTESIAN ENERGY GRADIENT USING VARIATIONAL C FINITE-DIFFERENCE DERIVATIVES. RETURNS GRADIENT IN KCAL/ANGSTROM. C SEE SUBROUTINE FOCK FOR MORE EXTENSIVE COMMENTS ON ENERGY TERMS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" #include "constants.h" DIMENSION GXYZ(3,MAXATM) C C LOCAL VARIABLES: C LOGICAL ORBSI,ORBSJ,ORBSIJ,SHARE1,SHARE2 integer resinter DIMENSION AGI(4),BGI(4),CGI(4),AGJ(4),BGJ(4),CGJ(4),HIJ(9,9), . REPIJ(10,10),EENUC(2,10),DHII(9,9),DHJJ(9,9),DHIJ(9,9), . DFII(9,9),DFJJ(9,9),DFIJ(9,9),DREPIJ(10,10),PII(9,9), . PJJ(9,9),PIJ(9,9),REPIJ0(10,10),DCART(2),SCALE(2) DIMENSION DXREPIJ(10,10),DYREPIJ(10,10),DZREPIJ(10,10) DIMENSION DXHIJ(9,9),DYHIJ(9,9),DZHIJ(9,9) DIMENSION DXREPIJdum(10,10),DYREPIJdum(10,10),DZREPIJdum(10,10) DIMENSION DXHIJdum(9,9),DYHIJdum(9,9),DZHIJdum(9,9) c DIMENSION DXEENUC(2,10),DYEENUC(2,10),DZEENUC(2,10) DIMENSION IREPUL(4,4) DATA IREPUL /1,2,3,4, . 2,5,6,7, . 3,6,8,9, . 4,7,9,10/ SAVE IREPUL C call etimer(t1) DDXYZ = 5.0D-6 DO 20 I=1,NATOMS DO 10 J=1,3 GXYZ(J,I) = 0.0D0 10 CONTINUE 20 CONTINUE IF(NATOMS.EQ.1) RETURN C C LOOP OVER ALL POSSIBLE PAIRS OF ATOMS. 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 DO 1000 I=2,NATOMS IAI = IATNUM(I) IF(IAI.EQ.0) GO TO 1000 C C UNPERTURBED ATOM I COORDINATES: C XI = XYZ(1,I) YI = XYZ(2,I) ZI = XYZ(3,I) C NORBSI = NATORB(IAI) ORBSI = NORBSI.GT.0 NREPI = MAX((NORBSI*(NORBSI+1))/2,1) NPQI = NQUANT(IAI) XISI = EXPNT(0,IAI) XIPI = EXPNT(1,IAI) XIDI = EXPNT(2,IAI) 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) ACI = ACORE(IAI) DO 40 IGAUS=1,4 AGI(IGAUS) = AGAUS(IGAUS,IAI) BGI(IGAUS) = BGAUS(IGAUS,IAI) CGI(IGAUS) = CGAUS(IGAUS,IAI) 40 CONTINUE ZNUCI = ZCHG(IAI) C C TEMPORARILY STORE DIAGONAL BLOCK OF DENSITY MATRIX FOR ATOM I. C IF(ORBSI)THEN IJ = IIMAT(I) DO 70 II=1,NORBSI DO 60 JJ=1,II PII(II,JJ) = PDIAG(IJ) PII(JJ,II) = PDIAG(IJ) IJ = IJ + 1 60 CONTINUE 70 CONTINUE ENDIF C C GET INTERSECTION WITH ALL OTHER ATOMS J < I. C DO 900 J=1,I-1 IAJ = IATNUM(J) IF(IAJ.EQ.0) GO TO 900 resinter = interres(i,j) NORBSJ = NATORB(IAJ) ORBSJ = NORBSJ.GT.0 ORBSIJ = ORBSI.AND.ORBSJ IF(IPAIR(IJBOND).EQ.J)THEN SHARE1 = NSHARE(1,IJBOND).NE.0 SHARE2 = .TRUE. ELSE SHARE1 = .FALSE. SHARE2 = .FALSE. ENDIF NREPJ = MAX((NORBSJ*(NORBSJ+1))/2,1) NPQJ = NQUANT(IAJ) XISJ = EXPNT(0,IAJ) XIPJ = EXPNT(1,IAJ) XIDJ = EXPNT(2,IAJ) 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) ACJ = ACORE(IAJ) DO 90 IGAUS=1,4 AGJ(IGAUS) = AGAUS(IGAUS,IAJ) BGJ(IGAUS) = BGAUS(IGAUS,IAJ) CGJ(IGAUS) = CGAUS(IGAUS,IAJ) 90 CONTINUE 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 110 II=1,NORBSJ DO 100 JJ=1,II PJJ(II,JJ) = PDIAG(IJ) PJJ(JJ,II) = PDIAG(IJ) IJ = IJ + 1 100 CONTINUE 110 CONTINUE ENDIF C C TEMPORARILY STORE I-J OFF-DIAGONAL BLOCK OF DENSITY MATRIX. C IF(ORBSIJ.AND.SHARE1)THEN IJ = IJMAT(IJBOND) DO 130 II=1,NORBSI DO 120 JJ=1,NORBSJ PIJ(II,JJ) = PDIAT(IJ) IJ = IJ + 1 120 CONTINUE 130 CONTINUE ENDIF C C IF THIS IS A DIRECT CALCULATION, THEN WE NEED TO COMPUTE C CORE-CORE AND TWO-ELECTRON INTEGRALS CORRESPONDING TO C THE UNPERTURBED ATOM I COORDINATES X0,Y0,Z0. OTHERWISE C THEY CAN BE RETRIEVED FROM THE APPROPRIATE GLOBAL ARRAYS. C c IF(DIRECT.AND..NOT.CENTRL)THEN c ILEVEL = 3 ilevel=2 RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2) 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, . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,DXREPIJ,DYREPIJ,DZREPIJ, . DXHIJ,DYHIJ,DZHIJ, . DXENNIJ,DYENNIJ,DZENNIJ, . EENUC,resinter) c . RIJ0,ILEVEL,ENNIJ0,HIJ,REPIJ0,EENUC) c ELSEIF(.NOT.DIRECT)THEN c IJ = IJREP(IJFULL) c DO 220 II=1,NREPI c DO 210 JJ=1,NREPJ c REPIJ0(II,JJ) = EEREP(IJ) c IJ = IJ + 1 c 210 CONTINUE c 220 CONTINUE c ENNIJ0 = ENUCLR(IJFULL) c ENDIF DO 800 ICART=1,3 C C SHIFT EACH CARTESIAN COORDINATE OF ATOM I BY A SMALL AMOUNT C AND COMPUTE THE CORRESPONDING CHANGES IN DIATOMIC INTERACTIONS C BETWEEN ATOMS I AND J. C XYZ0 = XYZ(ICART,I) XYZ(ICART,I) = XYZ0 + DDXYZ XI = XYZ(1,I) YI = XYZ(2,I) ZI = XYZ(3,I) C C COMPUTE CORE-CORE, TWO-ELECTRON, AND POSSIBLY DIATOMIC C ONE-ELECTRON INTEGRALS FOR I-J INTERSECTION WITH PERTURBED C COORDINATE FOR ATOM I. CORE-ELECTRON MAY BE SKIPPED C BECAUSE THEY ARE REFORMED FROM TWO-ELECTRON INTEGRALS. C RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2) ILEVEL = 2 C C SKIP HIJ IF PIJ IS ZERO. C IF(.NOT.SHARE1) ILEVEL = 3 C 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, . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,DXREPIJdum,DYREPIJdum,DZREPIJdum, . DXHIJdum,DYHIJdum,DZHIJdum, . DXENNIJdum,DYENNIJdum,DZENNIJdum, c . DXEENUC,DYEENUC,DZEENUC, . EENUC,resinter) c . RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC) c DO 800 ICART=1,3 c write(*,*) c do ii=1,10 c do jj=1,10 c write(*,99)'D',icart,'REPIJ(',jj,',',ii,') =', c & drepij(jj,ii) c enddo c enddo c 99 format(A1,I1,A6,I2,A1,I2,A3,2x,f15.8) c c c write(*,*) c do ii=1,9 c do jj=1,9 c write(*,98)'D',icart,'HIJ(',jj,',',ii,') =', c & dhij(jj,ii) c enddo c enddo c 98 format(A1,I1,A4,I2,A1,I2,A3,2x,f15.8) c write(*,*) c write(*,97)'D',icart,'ENNIJ =',DENNIJ c 97 format(A1,I1,A7,2x,f15.8) C C DO REVERSE STEP IF THE USER WANTS A CENTRAL DIFFERENCE. C c IF(CENTRL)THEN XYZ(ICART,I) = XYZ0 - DDXYZ XI = XYZ(1,I) YI = XYZ(2,I) ZI = XYZ(3,I) RIJOLD = RIJ RIJ = DSQRT((XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2) 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, . RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,DXREPIJdum,DYREPIJdum,DZREPIJdum, . DXHIJdum,DYHIJdum,DZHIJdum, . DXENNIJdum,DYENNIJdum,DZENNIJdum, . EENUC,resinter) c . RIJ,ILEVEL,DENNIJ,DHIJ,DREPIJ,EENUC) if (icart.eq.1) then do iii=1,10 do jjj=1,10 DREPIJ(jjj,iii)= DXREPIJ(jjj,iii)/Bohr2Ang enddo enddo c do iii=1,9 c do jjj=1,9 c DHIJ(jjj,iii)= DXHIJ(jjj,iii)/Bohr2Ang c enddo c enddo DENNIJ= DXENNIJ/Bohr2Ang elseif (icart.eq.2) then do iii=1,10 do jjj=1,10 DREPIJ(jjj,iii)= DYREPIJ(jjj,iii)/Bohr2Ang enddo enddo c do iii=1,9 c do jjj=1,9 c DHIJ(jjj,iii)= DYHIJ(jjj,iii)/Bohr2Ang c enddo c enddo DENNIJ=DYENNIJ/Bohr2Ang elseif (icart.eq.3) then do iii=1,10 do jjj=1,10 DREPIJ(jjj,iii)= DZREPIJ(jjj,iii)/Bohr2Ang enddo enddo c do iii=1,9 c do jjj=1,9 c DHIJ(jjj,iii)= DZHIJ(jjj,iii)/Bohr2Ang c enddo c enddo DENNIJ= DZENNIJ/Bohr2Ang endif C C COMBINE DIFFERENCE QUANTITIES FROM FORWARD AND REVERSE C STEPS TO GET CHANGES IN OFF-DIAGONAL BLOCK OF 1-ELECTRON C MATRIX AND ELECTRON REPULSIONS. C c DENNIJ = (ENNIJ - DENNIJ)*0.5D0/ddxyz IF(ORBSIJ.AND.SHARE1)THEN DO 250 II=1,NORBSI DO 240 JJ=1,NORBSJ DHIJ(II,JJ) = (HIJ(II,JJ) - DHIJ(II,JJ))*0.5D0/ddxyz DFIJ(II,JJ) = DHIJ(II,JJ) 240 CONTINUE 250 CONTINUE ENDIF C C NEED DREPIJ EVEN IF SPARKLES TO CALCULATE CHANGES C IN ELECTRON-CORE ATTRACTIONS. C c DO 270 II=1,NREPI c DO 260 JJ=1,NREPJ c DREPIJ(II,JJ) = (REPIJ(II,JJ) - DREPIJ(II,JJ))*0.5D0/ddxyz c 260 CONTINUE c 270 CONTINUE c XYZ(ICART,I) = XYZ0 + DDXYZ c XI = XYZ(1,I) c YI = XYZ(2,I) c ZI = XYZ(3,I) c RIJ = RIJOLD c ELSE C C COMPUTE CHANGES FROM FOWARD STEP ALONE. C c DENNIJ = ENNIJ - ENNIJ0 c IF(ORBSIJ.AND.SHARE1)THEN c IJ = IJMAT(IJBOND) c DO 300 II=1,NORBSI c DO 290 JJ=1,NORBSJ c DHIJ(II,JJ) = HIJ(II,JJ) - HDIAT(IJ) c DFIJ(II,JJ) = DHIJ(II,JJ) c IJ = IJ + 1 c 290 CONTINUE c 300 CONTINUE c ENDIF c DO 320 II=1,NREPI c DO 310 JJ=1,NREPJ c c DREPIJ(II,JJ) = REPIJ(II,JJ) - REPIJ0(II,JJ) c 310 CONTINUE c 320 CONTINUE c ENDIF C C SINCE ELECTRON-NUCLEAR ATTRACTIONS ARE NOT STORED IN ANY C SORT OF GLOBAL ARRAY, THE CHANGES IN THE DIAGONAL BLOCKS C OF THE 1-ELECTRON MATRIX MUST BE COMPUTED DIRECTLY FROM C THE CHANGES IN THE ELECTRON REPULSIONS. COPY CHANGES TO C FOCK MATRIX. C IF(ORBSI)THEN DO 340 II=1,NORBSI DO 330 JJ=1,II DHII(II,JJ) = -ZNUCJ*DREPIJ(IREPUL(II,JJ),1) DFII(II,JJ) = DHII(II,JJ) 330 CONTINUE 340 CONTINUE ENDIF IF(ORBSJ)THEN DO 360 II=1,NORBSJ DO 350 JJ=1,II DHJJ(II,JJ) = -ZNUCI*DREPIJ(1,IREPUL(II,JJ)) DFJJ(II,JJ) = DHJJ(II,JJ) 350 CONTINUE 360 CONTINUE ENDIF C C ADD ON 2-CENTER 2-ELECTRON TERMS TO THE CHANGES IN THE C OFF-DIAGONAL BLOCK OF THE FOCK MATRIX C IF(ORBSIJ.AND.SHARE1)THEN DO 400 IG=1,NORBSI DO 390 JG=1,NORBSJ DGIJ = 0.0D0 DO 380 IP=1,NORBSI IR = IREPUL(IG,IP) DO 370 JP=1,NORBSJ JR = IREPUL(JG,JP) DGIJ = DGIJ + PIJ(IP,JP)*DREPIJ(IR,JR) 370 CONTINUE 380 CONTINUE DFIJ(IG,JG) = DFIJ(IG,JG) - 0.5D0*DGIJ 390 CONTINUE 400 CONTINUE ENDIF C C ADD ON 2-CENTER 2-ELECTRON TERMS TO CHANGES IN DIAGONAL C BLOCKS OF FOCK MATRIX. C IF(ORBSIJ)THEN DO 440 IG=1,NORBSI DO 430 JG=1,IG JJG = IREPUL(JG,JG) IJG = IREPUL(IG,JG) DGIJ = 0.0D0 DO 420 IP=1,NORBSJ DO 410 JP=1,NORBSJ IJP = IREPUL(IP,JP) DGIJ = DGIJ + DREPIJ(IJG,IJP)*PJJ(IP,JP) 410 CONTINUE 420 CONTINUE DFII(IG,JG) = DFII(IG,JG) + DGIJ 430 CONTINUE 440 CONTINUE DO 480 IG=1,NORBSJ DO 470 JG=1,IG JJG = IREPUL(JG,JG) IJG = IREPUL(IG,JG) DGIJ = 0.0D0 DO 460 IP=1,NORBSI DO 450 JP=1,NORBSI IJP = IREPUL(IP,JP) DGIJ = DGIJ + DREPIJ(IJP,IJG)*PII(IP,JP) 450 CONTINUE 460 CONTINUE DFJJ(IG,JG) = DFJJ(IG,JG) + DGIJ 470 CONTINUE 480 CONTINUE ENDIF C C WE NOW HAVE ALL CHANGES IN THE 1-ELECTRON AND FOCK MATRICES C DUE TO THE I-J INTERACTION. NOW COMPUTE THE CORRESPONDING C CHANGE IN THE TOTAL ENERGY. C DIJ = 0.0D0 DII = 0.0D0 DJJ = 0.0D0 IF(ORBSI)THEN DO 520 II=1,NORBSI DO 510 JJ=1,II DDIJ = (DHII(II,JJ) + DFII(II,JJ))*PII(II,JJ) DIJ = DIJ + DDIJ 510 CONTINUE DII = DII + DDIJ 520 CONTINUE ENDIF IF(ORBSJ)THEN DO 540 II=1,NORBSJ DO 530 JJ=1,II DDIJ = (DHJJ(II,JJ) + DFJJ(II,JJ))*PJJ(II,JJ) DIJ = DIJ + DDIJ 530 CONTINUE DJJ = DJJ + DDIJ 540 CONTINUE ENDIF DELECT = DIJ - 0.5D0*(DII + DJJ) IF(ORBSIJ.AND.SHARE1)THEN DO 560 II=1,NORBSI DO 550 JJ=1,NORBSJ DELECT = DELECT + (DHIJ(II,JJ)+DFIJ(II,JJ))*PIJ(II,JJ) 550 CONTINUE 560 CONTINUE ENDIF C C UPDATE GRADIENT ENTRIES FOR ATOMS I AND J: C DGRAD1 = DELECT + DENNIJ GXYZ(ICART,I) = GXYZ(ICART,I) + DGRAD1 GXYZ(ICART,J) = GXYZ(ICART,J) - DGRAD1 C C RESTORE CARTESIAN COORDINATE OF ATOM I. C c XYZ(ICART,I) = XYZ0 800 CONTINUE C C UPDATE POINTERS FOR FULL AND BOND-BASED PAIRLIST. C IJFULL = IJFULL + 1 IF(SHARE2) IJBOND = IJBOND + 1 900 CONTINUE 1000 CONTINUE C C ADD MOLECULAR MECHANICS TERMS FOR PEPTIDE TORSION BARRIER. C c IF(N2PEP.GT.0)THEN c DCART(1) = DDXYZ c SCALE(1) = 1.0D0 c NSTEP = 1 C C DO FORWARD AND REVERSE STEPS FOR CENTRAL DIFFERENCE. C c IF(CENTRL)THEN c DCART(2) = -DDXYZ c SCALE(1) = 0.5D0 c SCALE(2) = -0.5D0 c NSTEP = 2 c ENDIF c DO 1200 ISTEP=1,NSTEP C C LOOP OVER EACH ATOM IN EACH TORSIONAL ANGLE AND GET THE CHANGE C IN ENERGY DUE TO MOVING THE ATOM BY DCART(ISTEP). C c DO 1180 IPEP=1,N2PEP c EIABC0 = EIABC(IPEP) c DO 1160 I=1,4 c IATM = IABC(I,IPEP) c DO 1140 ICART=1,3 c XYZ0 = XYZ(ICART,IATM) c XYZ(ICART,IATM) = XYZ0 + DCART(ISTEP) c CALL DIHEDR(XYZ,IABC(1,IPEP),IABC(2,IPEP), c . IABC(3,IPEP),IABC(4,IPEP),DIH) c DEPEP = EPEPMX*SIN(DIH)**2 - EIABC0 c GXYZ(ICART,IATM) = GXYZ(ICART,IATM) + DEPEP*SCALE(ISTEP) c XYZ(ICART,IATM) = XYZ0 c 1140 CONTINUE c 1160 CONTINUE c 1180 CONTINUE c 1200 CONTINUE c ENDIF C C CONVERT ENERGY DIFFERENCES TO DERIVATIVES IN KCAL/ANGSTROM. C FACTR = 23.061D0!/Bohr2Ang !/DDXYZ DO 1400 I=1,NATOMS DO 1300 J=1,3 GXYZ(J,I) = GXYZ(J,I)*FACTR 1300 CONTINUE 1400 CONTINUE call etimer(t2) tgrad = t2 - t1 RETURN END