mopac606_nbo/src/fmat.f

   1       SUBROUTINE FMAT(FMATRX, NREAL, TSCF, TDER, DELDIP, HEAT)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       INCLUDE 'SIZES'
   4       DIMENSION FMATRX(*), DELDIP(3,*)
   5 ***********************************************************************
   6 *
   7 *  VALUE CALCULATES THE SECOND-ORDER OF THE ENERGY WITH
   8 *        RESPECT TO THE CARTESIAN COORDINATES I AND J AND PLACES IT
   9 *        IN FMATRX
  10 *
  11 *  ON INPUT NATOMS  = NUMBER OF ATOMS IN THE SYSTEM.
  12 *           XPARAM  = INTERNAL COORDINATES OF MOLECULE STORED LINEARLY
  13 *
  14 *  VARIABLES USED
  15 *           COORDL  = ARRAY OF CARTESIAN COORDINATES, STORED LINEARLY.
  16 *           I       = INDEX OF CARTESIAN COORDINATE.
  17 *           J       = INDEX OF CARTESIAN COORDINATE.
  18 *
  19 *  ON OUTPUT FMATRX = SECOND DERIVATIVE OF THE ENERGY WITH RESPECT TO
  20 *                    CARTESIAN COORDINATES I AND J.
  21 ***********************************************************************
  22       COMMON /KEYWRD/ KEYWRD
  23       COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
  24      1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
  25       COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY, DUMY(MAXPAR)
  26       COMMON /DENSTY/ P(MPACK),PDUMY(2,MPACK)
  27       COMMON /TIMDMP/ TLEFT, TDUMP
  28       COMMON /ATMASS/ ATMASS(NUMATM)
  29       COMMON /TIME  / TIME0
  30       COMMON /CORE  / CORE(107)
  31       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
  32      1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
  33      2                NCLOSE,NOPEN,NDUMY,FRACT
  34       COMMON /COORD / COORD(3,NUMATM)
  35       COMMON /NLLCOM/ EVECS(MAXPAR*MAXPAR),BMAT(MAXPAR,MAXPAR*2)
  36       DIMENSION GRAD(MAXPAR),
  37      1GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3), G2OLD(MAXPAR)
  38      2, EIGS(MAXPAR), G2RAD(MAXPAR),
  39      3 FCONST(MAXPAR)
  40       CHARACTER*241 KEYWRD
  41       SAVE  FACT
  42       LOGICAL DEBUG, RESTRT, PRNT, RESFIL, PRECIS, BIG, LOG
  43       EQUIVALENCE (COORD(1,1),COORDL(1))
  44       DATA FACT/6.95125D-3/
  45 C
  46 C    FACT IS THE CONVERSION FACTOR FROM KCAL/MOLE TO ERGS
  47 C
  48 C SET UP CONSTANTS AND FLAGS
  49       NA(1)=99
  50 C
  51 C  SET UP THE VARIABLES IN XPARAM ANDLOC,THESE ARE IN CARTESIAN COORDINA
  52 C
  53       NUMAT=0
  54 C$DOIT ASIS
  55       DO 10 I=1,NATOMS
  56          IF(LABELS(I).NE.99.AND.LABELS(I).NE.107) THEN
  57             NUMAT=NUMAT+1
  58             LABELS(NUMAT)=LABELS(I)
  59          ENDIF
  60    10 CONTINUE
  61       NATOMS=NUMAT
  62 C
  63 C   THIS IS A QUICK, IF CLUMSY, WAY TO CALCULATE NUMAT, AND TO REMOVE
  64 C   THE DUMMY ATOMS FROM THE ARRAY LABELS.
  65 C
  66       NVAR=NUMAT*3
  67       DO 20 I=1,NUMAT
  68          LOC(1,(I-1)*3+1)=I
  69          LOC(2,(I-1)*3+1)=1
  70 C
  71          LOC(1,(I-1)*3+2)=I
  72          LOC(2,(I-1)*3+2)=2
  73 C
  74          LOC(1,(I-1)*3+3)=I
  75          LOC(2,(I-1)*3+3)=3
  76    20 CONTINUE
  77       LIN=(NVAR*(NVAR+1))/2
  78       DO 30 I=1,LIN
  79    30 FMATRX(I)=0.D0
  80       PRNT   =(INDEX(KEYWRD,'IRC=') .EQ. 0)
  81       LOG    =(INDEX(KEYWRD,'NOLOG') .EQ. 0)
  82       PRECIS =(INDEX(KEYWRD,'PREC') .NE. 0)
  83       RESTRT =(INDEX(KEYWRD,'RESTART') .NE. 0)
  84       IF(INDEX(KEYWRD,'NLLSQ') .NE. 0) RESTRT=.FALSE.
  85       DEBUG =(INDEX(KEYWRD,'FMAT') .NE. 0)
  86       BIG    =(INDEX(KEYWRD,'LARGE') .NE. 0 .AND. DEBUG)
  87       IF(PRNT)WRITE(6,'(//4X,''FIRST DERIVATIVES WILL BE USED IN THE''
  88      1,'' CALCULATION OF SECOND DERIVATIVES'')')
  89       TLAST=TLEFT
  90       RESFIL=.FALSE.
  91       IF(RESTRT) THEN
  92          ISTART = 0
  93          I=0
  94          CALL FORSAV(TOTIME,DELDIP,ISTART,FMATRX, COORD, NVAR,HEAT,
  95      1                EVECS,JSTART,FCONST)
  96          KOUNTF=(ISTART*(ISTART+1))/2
  97          ISTART=ISTART+1
  98          JSTART=JSTART+1
  99          TIME2 = SECOND()
 100       ELSE
 101          KOUNTF=0
 102          TOTIME=0.D0
 103          IF (TSCF.GT.0.D0)TLEFT=TLEFT-TSCF-TDER
 104          ISTART=1
 105       ENDIF
 106 C CALCULATE FMATRX
 107       IF(ISTART.GT.1) THEN
 108          ESTIME=(NVAR-ISTART+1)*TOTIME/(ISTART-1.D0)
 109       ELSE
 110          ESTIME=NVAR*(TSCF+TDER)*2.D0
 111          IF (PRECIS) ESTIME=ESTIME*2.D0
 112       ENDIF
 113       IF(TSCF.GT.0)
 114      1WRITE(6,'(/10X,''ESTIMATED TIME TO COMPLETE CALCULATION =''
 115      2,F9.2,'' SECONDS'')')ESTIME
 116       IF(RESTRT) THEN
 117          IF(ISTART.LE.NVAR)
 118      1    WRITE(6,'(/10X,''STARTING AGAIN AT LINE'',18X,I4)')ISTART
 119          WRITE(6,'(/10X,''TIME USED UP TO RESTART ='',F22.2)')TOTIME
 120       ENDIF
 121       LU=KOUNTF
 122       NUMAT=NVAR/3
 123       DO 50 I=1,NVAR
 124    50 EIGS(I)=0.D0
 125       DO 120 I=ISTART,NVAR
 126          TIME2 = SECOND()
 127          DELTA=1.D0/120.D0
 128          IF(PRECIS)THEN
 129 C
 130 C   DETERMINE A GOOD STEP SIZE
 131 C
 132             G2OLD(1)=100.D0
 133             COORDL(I)=COORDL(I)+DELTA
 134             CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
 135             COORDL(I)=COORDL(I)-DELTA
 136             DELTA=DELTA*10.D0/SQRT(DOT(G2OLD,G2OLD,NVAR))
 137 C
 138 C   CONSTRAIN DELTA TO A 'REASONABLE' VALUE
 139 C
 140             DELTA=MIN(0.05D0,MAX(0.005D0,DELTA))
 141             IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
 142      1A
 143             G2OLD(1)=100.D0
 144             COORDL(I)=COORDL(I)+DELTA
 145             CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
 146             IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +1.0*DELTA',
 147      1SQRT(DOT(G2OLD,G2OLD,NVAR))
 148             COORDL(I)=COORDL(I)-DELTA*2.D0
 149             G2RAD(1)=100.D0
 150             CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., G2RAD, .TRUE.)
 151             COORDL(I)=COORDL(I)+DELTA
 152             IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -1.0*DELTA',
 153      1SQRT(DOT(G2RAD,G2RAD,NVAR))
 154          ELSE
 155             IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
 156      1A
 157          ENDIF
 158          COORDL(I)=COORDL(I)+0.5D0*DELTA
 159          GROLD(1)=100.D0
 160          CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.)
 161          IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +0.5*DELTA',
 162      1SQRT(DOT(GROLD,GROLD,NVAR))
 163          CALL CHRGE(P,Q)
 164          DO 60 II=1,NUMAT
 165    60    Q(II)=CORE(LABELS(II))-Q(II)
 166          SUM = DIPOLE(P,Q,COORDL,DELDIP(1,I),0)
 167          COORDL(I)=COORDL(I)-DELTA
 168          GRAD(1)=100.D0
 169          CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .TRUE.)
 170          IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -0.5*DELTA',
 171      1SQRT(DOT(GRAD,GRAD,NVAR))
 172          CALL CHRGE(P,Q)
 173          DO 70 II=1,NUMAT
 174    70    Q(II)=CORE(LABELS(II))-Q(II)
 175          SUM = DIPOLE(P,Q,COORDL,DEL2,0)
 176          COORDL(I)=COORDL(I)+DELTA*0.5D0
 177          DELDIP(1,I)=(DELDIP(1,I)-DEL2(1))*0.5D0/DELTA
 178          DELDIP(2,I)=(DELDIP(2,I)-DEL2(2))*0.5D0/DELTA
 179          DELDIP(3,I)=(DELDIP(3,I)-DEL2(3))*0.5D0/DELTA
 180          LL=LU+1
 181          LU=LL+I-1
 182          L=0
 183          IF(PRECIS)THEN
 184             DO 80 KOUNTF=LL,LU
 185                L=L+1
 186 C
 187 C       G2OLD = X + 1.0*DELTA
 188 C       GROLD = X + 0.5*DELTA
 189 C       GRAD  = X - 0.5*DELTA
 190 C       G2RAD = X - 1.0*DELTA
 191 C
 192                DUMY(L)= (8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
 193      1          /DELTA*FACT/24.D0
 194                EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
 195      1          /DELTA**3*FACT/56.D0
 196 C
 197 C  CORRECT FOR 4'TH ORDER CONTAMINATION
 198 C
 199 C#             CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
 200 C#             DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
 201                FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
 202    80       CONTINUE
 203             L=L-1
 204             DO 90 K=I,NVAR
 205                L=L+1
 206                KK=(K*(K-1))/2+I
 207                DUMY(L)=(8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
 208      1          /DELTA*FACT/24.D0
 209                EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
 210      1          /DELTA**3*FACT/56.D0
 211 C
 212 C  CORRECT FOR 4'TH ORDER CONTAMINATION
 213 C
 214 C#             CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
 215 C#             DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
 216                FMATRX(KK)=FMATRX(KK)+DUMY(L)
 217    90       CONTINUE
 218          ELSE
 219             DO 100 KOUNTF=LL,LU
 220                L=L+1
 221                DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
 222                FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
 223   100       CONTINUE
 224             L=L-1
 225             DO 110 K=I,NVAR
 226                L=L+1
 227                KK=(K*(K-1))/2+I
 228                DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
 229                FMATRX(KK)=FMATRX(KK)+DUMY(L)
 230   110       CONTINUE
 231          ENDIF
 232          IF(BIG)THEN
 233             WRITE(6,'(A)')' CONTRIBUTIONS TO F-MATRIX'
 234             WRITE(6,'(A)')' ELEMENT  +1.0*DELTA  +0.5*DELTA  -0.5*DEL'
 235      1//'TA  -1.0*DELTA   2''ND ORDER 4TH ORDER'
 236             WRITE(6,'(I7,6F12.6)')(L,G2OLD(L),GROLD(L),GRAD(L),G2RAD(L),
 237      1DUMY(L),EIGS(L),L=1,NVAR)
 238          ENDIF
 239          TIME3 = SECOND()
 240          TSTEP=TIME3-TIME2
 241          TLEFT= MAX(0.1D0,TLEFT-TSTEP)
 242          IF(TSTEP.GT.1.D6)TSTEP=TSTEP-1.D6
 243          TOTIME= TOTIME+TSTEP
 244          IF(RESFIL)THEN
 245             WRITE(6,'('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL ='
 246      1',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
 247             WRITE(0,'('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL ='
 248      1',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
 249          IF(LOG)WRITE(11,'('' STEP:'',I4,'' RESTART FILE WRITTEN, '',
 250      +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
 251             RESFIL=.FALSE.
 252          ELSE
 253             WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL =
 254      1'',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
 255             WRITE(0,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL =
 256      1'',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
 257         IF(LOG) WRITE(11,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, '',
 258      +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
 259          ENDIF
 260          ESTIM = TOTIME/I
 261          IF(TLAST-TLEFT.GT.TDUMP)THEN
 262             TLAST=TLEFT
 263             RESFIL=.TRUE.
 264             JSTART=1
 265             II=I
 266             CALL FORSAV(TOTIME,DELDIP,II,FMATRX, COORD,NVAR,HEAT,
 267      1                EVECS,JSTART,FCONST)
 268          ENDIF
 269          IF(I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN
 270             WRITE(6,'(//10X,''- - - - - - - TIME UP - - - - - - -'',//)'
 271      1)
 272             WRITE(6,'(/10X,'' POINT REACHED ='',I4)')I
 273             WRITE(6,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')')
 274             WRITE(6,'(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2,
 275      1'' SECONDS'')')ESTIM
 276             JSTART=1
 277             II=I
 278             CALL FORSAV(TOTIME,DELDIP,II,FMATRX, COORD,NVAR,HEAT,
 279      1                EVECS,JSTART,FCONST)
 280             WRITE(6,'(//10X,''FORCE MATRIX WRITTEN TO DISK'')')
 281             NREAL=-1
 282             RETURN
 283          ENDIF
 284   120 CONTINUE
 285       DO 130 I=1,NATOMS
 286          IF(ATMASS(I).LT.1.D-20.AND.LABELS(I).LT.99)THEN
 287             CALL FORSAV(TOTIME,DELDIP,NVAR,FMATRX, COORD,NVAR,HEAT,
 288      1                EVECS,ILOOP,FCONST)
 289             WRITE(6,'(A)')' AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART
 290      1'
 291             WRITE(6,'(A)')' FILE HAS BEEN WRITTEN AND THE JOB STOPPED'
 292             STOP
 293          ENDIF
 294   130 CONTINUE
 295       IF(ISTART.LE.NVAR .AND. INDEX(KEYWRD,'ISOT') .NE. 0)
 296      1CALL FORSAV(TOTIME,DELDIP,NVAR,FMATRX, COORD,NVAR,HEAT,
 297      2                EVECS,ILOOP,FCONST)
 298       RETURN
 299       END