mopac606_nbo/src/dcart.f

   1       SUBROUTINE DCART (COORD,DXYZ)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       INCLUDE 'SIZES'
   4       DIMENSION COORD(3,*), DXYZ(3,*)
   5       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
   6      1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
   7      2                NCLOSE,NOPEN,NDUMY,FRACT
   8       COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
   9 C***********************************************************************
  10 C
  11 C    DCART CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
  12 C          CARTESIAN COORDINATES. THIS IS DONE BY FINITE DIFFERENCES.
  13 C
  14 C    THE MAIN ARRAYS IN DCART ARE:
  15 C        DXYZ   ON EXIT CONTAINS THE CARTESIAN DERIVATIVES.
  16 C
  17 C***********************************************************************
  18       COMMON /KEYWRD/ KEYWRD
  19       COMMON /EULER / TVEC(3,3), ID
  20       COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
  21       COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
  22       COMMON /DCARTC/ K1L,K2L,K3L,K1U,K2U,K3U
  23       COMMON /NUMCAL/ NUMCAL
  24       CHARACTER*241 KEYWRD
  25       DIMENSION PDI(171),PADI(171),PBDI(171),
  26      1CDI(3,2),NDI(2),LSTOR1(6), LSTOR2(6), ENG(3)
  27       LOGICAL DEBUG, FORCE, MAKEP, ANADER, LARGE
  28       EQUIVALENCE (LSTOR1(1),L1L), (LSTOR2(1), K1L)
  29       SAVE CHNGE, CHNGE2, ANADER, DEBUG, FORCE
  30       DATA ICALCN/0/
  31       DATA CHNGE /1.D-4/
  32       CHNGE2=CHNGE*0.5D0
  33 *
  34 * CHNGE IS A MACHINE-PRECISION DEPENDENT CONSTANT
  35 * CHNGE2=CHNGE/2
  36 *
  37       IF (ICALCN.NE.NUMCAL) THEN
  38          ICALCN=NUMCAL
  39          LARGE = (INDEX(KEYWRD,'LARGE') .NE. 0)
  40          ANADER= (INDEX(KEYWRD,'ANALYT') .NE. 0)
  41          DEBUG = (INDEX(KEYWRD,'DCART') .NE. 0)
  42          FORCE = (INDEX(KEYWRD,'PREC')+INDEX(KEYWRD,'FORCE') .NE. 0)
  43       ENDIF
  44       NCELLS=(L1U-L1L+1)*(L2U-L2L+1)*(L3U-L3L+1)
  45       DO 10 I=1,6
  46          LSTOR2(I)=LSTOR1(I)
  47    10 LSTOR1(I)=0
  48       IOFSET=(NCELLS+1)/2
  49       NUMTOT=NUMAT*NCELLS
  50       DO 20 I=1,NUMTOT
  51          DO 20 J=1,3
  52    20 DXYZ(J,I)=0.D0
  53       IF(ANADER) REWIND 2
  54       DO 130 II=1,NUMAT
  55          III=NCELLS*(II-1)+IOFSET
  56          IM1=II
  57          IF=NFIRST(II)
  58          IM=NMIDLE(II)
  59          IL=NLAST(II)
  60          NDI(2)=NAT(II)
  61          DO 30 I=1,3
  62    30    CDI(I,2)=COORD(I,II)
  63          DO 130 JJ=1,IM1
  64             JJJ=NCELLS*(JJ-1)
  65 C  FORM DIATOMIC MATRICES
  66             JF=NFIRST(JJ)
  67             JM=NMIDLE(JJ)
  68             JL=NLAST(JJ)
  69 C   GET FIRST ATOM
  70             NDI(1)=NAT(JJ)
  71             MAKEP=.TRUE.
  72             DO 120 IK=K1L,K1U
  73                DO 120 JK=K2L,K2U
  74                   DO 120 KL=K3L,K3U
  75                      JJJ=JJJ+1
  76                      KKK=KKK-1
  77                      DO 40 L=1,3
  78    40                CDI(L,1)=COORD(L,JJ)+TVEC(L,1)*IK+TVEC(L,2)*JK+TVEC
  79      1(L,3)*KL
  80                      IF(.NOT. MAKEP) GOTO 90
  81                      MAKEP=.FALSE.
  82                      IJ=0
  83                      DO 50 I=JF,JL
  84                         K=I*(I-1)/2+JF-1
  85                         DO 50 J=JF,I
  86                            IJ=IJ+1
  87                            K=K+1
  88                            PADI(IJ)=PA(K)
  89                            PBDI(IJ)=PB(K)
  90    50                PDI(IJ)=P(K)
  91 C GET SECOND ATOM FIRST ATOM INTERSECTION
  92                      DO 80 I=IF,IL
  93                         L=I*(I-1)/2
  94                         K=L+JF-1
  95                         DO 60 J=JF,JL
  96                            IJ=IJ+1
  97                            K=K+1
  98                            PADI(IJ)=PA(K)
  99                            PBDI(IJ)=PB(K)
 100    60                   PDI(IJ)=P(K)
 101                         K=L+IF-1
 102                         DO 70 L=IF,I
 103                            K=K+1
 104                            IJ=IJ+1
 105                            PADI(IJ)=PA(K)
 106                            PBDI(IJ)=PB(K)
 107    70                   PDI(IJ)=P(K)
 108    80                CONTINUE
 109    90                CONTINUE
 110                      IF(II.EQ.JJ) GOTO  120
 111                      IF(ANADER)THEN
 112                         CALL ANALYT(PDI,PADI,PBDI,CDI,NDI,JF,JL,IF,IL
 113      1,                 ENG)
 114                         DO 100 K=1,3
 115                            DXYZ(K,III)=DXYZ(K,III)-ENG(K)
 116   100                   DXYZ(K,JJJ)=DXYZ(K,JJJ)+ENG(K)
 117                      ELSE
 118                         IF( .NOT. FORCE) THEN
 119                            CDI(1,1)=CDI(1,1)+CHNGE2
 120                            CDI(2,1)=CDI(2,1)+CHNGE2
 121                            CDI(3,1)=CDI(3,1)+CHNGE2
 122                            CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM
 123      1,IL,                 AA,1)
 124                         ENDIF
 125                         DO 110 K=1,3
 126                            IF( FORCE )THEN
 127                               CDI(K,2)=CDI(K,2)-CHNGE2
 128                               CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF
 129      1,IM,IL,                 AA,1)
 130                            ENDIF
 131                            CDI(K,2)=CDI(K,2)+CHNGE
 132                            CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM
 133      1,IL,                 EE,2)
 134                            CDI(K,2)=CDI(K,2)-CHNGE2
 135                            IF( .NOT. FORCE) CDI(K,2)=CDI(K,2)-CHNGE2
 136                            DERIV=(AA-EE)*23.061D0/CHNGE
 137                            DXYZ(K,III)=DXYZ(K,III)-DERIV
 138                            DXYZ(K,JJJ)=DXYZ(K,JJJ)+DERIV
 139   110                   CONTINUE
 140                      ENDIF
 141   120       CONTINUE
 142   130 CONTINUE
 143       IF(NNHCO.NE.0)THEN
 144 C
 145 C   NOW ADD IN MOLECULAR-MECHANICS CORRECTION TO THE H-N-C=O TORSION
 146 C
 147          DEL=1.D-8
 148          DO 160 I=1,NNHCO
 149             DO 150 J=1,4
 150                DO 140 K=1,3
 151                   COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL
 152                   CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4,
 153      1I),ANGLE)
 154                   REFH=HTYPE(ITYPE)*SIN(ANGLE)**2
 155                   COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))+DEL*2.D0
 156                   CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4,
 157      1I),ANGLE)
 158                   COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL
 159                   HEAT=HTYPE(ITYPE)*SIN(ANGLE)**2
 160                   SUM=(REFH-HEAT)/(2.D0*DEL)
 161                   DXYZ(K,NHCO(J,I))=DXYZ(K,NHCO(J,I))-SUM
 162   140          CONTINUE
 163   150       CONTINUE
 164   160    CONTINUE
 165       ENDIF
 166       DO 170 I=1,6
 167   170 LSTOR1(I)=LSTOR2(I)
 168       IF (  .NOT. DEBUG) RETURN
 169       WRITE(6,'(//10X,''CARTESIAN COORDINATE DERIVATIVES'',//3X,
 170      1''NUMBER  ATOM '',5X,''X'',12X,''Y'',12X,''Z'',/)')
 171       IF(NCELLS.EQ.1)THEN
 172          WRITE(6,'(2I6,F13.6,2F13.6)')
 173      1 (I,NAT(I),(DXYZ(J,I),J=1,3),I=1,NUMTOT)
 174       ELSEIF(LARGE)THEN
 175          WRITE(6,'(2I6,F13.6,2F13.6)')
 176      1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I),J=1,3),I=1,NUMTOT)
 177       ELSE
 178          WRITE(6,'(2I6,F13.6,2F13.6)')
 179      1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I)+DXYZ(J,I+1)+DXYZ(J,I+2)
 180      2,J=1,3),I=1,NUMTOT,3)
 181       ENDIF
 182       IF (ANADER) REWIND 2
 183       RETURN
 184       END
 185       SUBROUTINE DHC (P,PA,PB,XI,NAT,IF,IM,IL,JF,JM,JL,DENER,MODE)
 186       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 187       DIMENSION P(*), PA(*), PB(*)
 188       DIMENSION XI(3,*),NFIRST(2),NMIDLE(2),NLAST(2),NAT(*)
 189 C***********************************************************************
 190 C
 191 C  DHC CALCULATES THE ENERGY CONTRIBUTIONS FROM THOSE PAIRS OF ATOMS
 192 C         THAT HAVE BEEN MOVED BY SUBROUTINE DERIV.
 193 C
 194 C***********************************************************************
 195       COMMON /KEYWRD/ KEYWRD
 196      1       /ONELEC/ USS(107),UPP(107),UDD(107)
 197       COMMON /EULER / TVEC(3,3), ID
 198       COMMON /NUMCAL/ NUMCAL
 199       SAVE ICALCN, WLIM, UHF
 200       CHARACTER*241 KEYWRD
 201       LOGICAL UHF, CUTOFF
 202       DIMENSION H(171), SHMAT(9,9), F(171),
 203      1          WJ(100), E1B(10), E2A(10), WK(100), W(100),
 204      2          WJS(100), WKS(100)
 205       DOUBLE PRECISION WJS, WKS
 206       DATA ICALCN /0/
 207       IF( ICALCN.NE.NUMCAL) THEN
 208          ICALCN=NUMCAL
 209          WLIM=4.D0
 210          IF(ID.EQ.0)WLIM=0.D0
 211          UHF=(INDEX(KEYWRD,'UHF') .NE. 0)
 212       ENDIF
 213       NFIRST(1)=1
 214       NMIDLE(1)=IM-IF+1
 215       NLAST(1)=IL-IF+1
 216       NFIRST(2)=NLAST(1)+1
 217       NMIDLE(2)=NFIRST(2)+JM-JF
 218       NLAST(2)=NFIRST(2)+JL-JF
 219       LINEAR=(NLAST(2)*(NLAST(2)+1))/2
 220       DO 10 I=1,LINEAR
 221          F(I)=0.D0
 222    10 H(I)=0.0D00
 223       DO 20 I=1,LINEAR
 224    20 F(I)=H(I)
 225       JA=NFIRST(2)
 226       JB=NLAST(2)
 227       JC=NMIDLE(2)
 228       IA=NFIRST(1)
 229       IB=NLAST(1)
 230       IC=NMIDLE(1)
 231       J=2
 232       I=1
 233       NJ=NAT(2)
 234       NI=NAT(1)
 235       CALL H1ELEC(NI,NJ,XI(1,1),XI(1,2),SHMAT)
 236       IF(NAT(1).EQ.102.OR.NAT(2).EQ.102) THEN
 237          K=(JB*(JB+1))/2
 238          DO 30 J=1,K
 239    30    H(J)=0.D0
 240       ELSE
 241          J1=0
 242          DO 40 J=JA,JB
 243             JJ=J*(J-1)/2
 244             J1=J1+1
 245             I1=0
 246             DO 40 I=IA,IB
 247                JJ=JJ+1
 248                I1=I1+1
 249                H(JJ)=SHMAT(I1,J1)
 250                F(JJ)=SHMAT(I1,J1)
 251    40    CONTINUE
 252       ENDIF
 253       KR=1
 254       IF(ID.EQ.0)THEN
 255          CALL ROTATE (NJ,NI,XI(1,2),XI(1,1),W(KR),KR,E2A,E1B,ENUCLR,100.
 256      1D0)
 257       ELSE
 258          CALL SOLROT (NJ,NI,XI(1,2),XI(1,1),WJ,WK,KR,E2A,E1B,ENUCLR,100.
 259      1D0)
 260       IF(MODE.EQ.1)CUTOFF=(WJ(1).LT.WLIM)
 261          IF(CUTOFF)THEN
 262             DO 50 I=1,KR-1
 263    50       WK(I)=0.D0
 264          ENDIF
 265          DO 60 I=1,KR-1
 266             WJS(I)=WJ(I)
 267             WKS(I)=WK(I)
 268    60    CONTINUE
 269       ENDIF
 270 C
 271 C    * ENUCLR IS SUMMED OVER CORE-CORE REPULSION INTEGRALS.
 272 C
 273       I2=0
 274       DO 70 I1=IA,IC
 275          II=I1*(I1-1)/2+IA-1
 276          DO 70 J1=IA,I1
 277             II=II+1
 278             I2=I2+1
 279             H(II)=H(II)+E1B(I2)
 280    70 F(II)=F(II)+E1B(I2)
 281       DO  80 I1=IC+1,IB
 282          II=(I1*(I1+1))/2
 283          F(II)=F(II)+E1B(1)
 284    80 H(II)=H(II)+E1B(1)
 285       I2=0
 286       DO 90 I1=JA,JC
 287          II=I1*(I1-1)/2+JA-1
 288          DO 90 J1=JA,I1
 289             II=II+1
 290             I2=I2+1
 291             H(II)=H(II)+E2A(I2)
 292    90 F(II)=F(II)+E2A(I2)
 293       DO 100 I1=JC+1,JB
 294          II=(I1*(I1+1))/2
 295          F(II)=F(II)+E2A(1)
 296   100 H(II)=H(II)+E2A(1)
 297       CALL FOCK2(F,P,PA,W, WJS, WKS,2,NFIRST,NMIDLE,NLAST)
 298       EE=HELECT(NLAST(2),PA,H,F)
 299       IF( UHF ) THEN
 300          DO 110 I=1,LINEAR
 301   110    F(I)=H(I)
 302          CALL FOCK2(F,P,PB,W, WJS, WKS,2,NFIRST,NMIDLE,NLAST)
 303          EE=EE+HELECT(NLAST(2),PB,H,F)
 304       ELSE
 305          EE=EE*2.D0
 306       ENDIF
 307       DENER=EE+ENUCLR
 308       RETURN
 309 C
 310       END