1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
SUBROUTINE DERITR(ERRFN,GEO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
DIMENSION GEO(3,NUMATM), ERRFN(MAXPAR)
************************************************************************
*
* DERITR CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
* INTERNAL COORDINATES. THIS IS DONE BY FINITE DIFFERENCES
* USING FULL SCF CALCULATIONS.
*
* THIS IS VERY TIME-CONSUMING, AND SHOULD ONLY BE USED WHEN
* NO OTHER DERIVATIVE CALCULATION WILL DO.
*
* THE MAIN ARRAYS IN DERIV ARE:
* LOC INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM
* INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE
* DERIVATIVE CALCULATION.
* GEO ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES.
*
************************************************************************
COMMON / EULER/ TVEC(3,3), ID
COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, DUMMY(MAXPAR)
COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
2 NCLOSE,NOPEN,NDUMY,FRACT
COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
1NA(NUMATM),NB(NUMATM),NC(NUMATM)
COMMON /GEOSYM/ NDEP, IDUMYS(MAXPAR,3)
COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
COMMON /ENUCLR/ ENUCLR
COMMON /NUMCAL/ NUMCAL
COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
COMMON /WMATRX/ WJ(N2ELEC), WK(N2ELEC)
COMMON /HMATRX/ H(MPACK)
COMMON /KEYWRD / KEYWRD
CHARACTER*241 KEYWRD
DIMENSION CHANGE(3), COORD(3,NUMATM)
1, XDERIV(3), XPARAM(MAXPAR), XJUC(3), W(N2ELEC)
DOUBLE PRECISION WJ, WK
SAVE IDELTA, XDERIV
LOGICAL DEBUG
EQUIVALENCE (W,WJ)
DATA ICALCN /0/
IF(ICALCN.NE.NUMCAL) THEN
DEBUG = (INDEX(KEYWRD,'DERITR') .NE. 0)
ICALCN=NUMCAL
*
* IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER
*
IDELTA=-3
CHANGE(1)= 10.D0**IDELTA
CHANGE(2)= 10.D0**IDELTA
CHANGE(3)= 10.D0**IDELTA
C
C CHANGE(I) IS THE STEP SIZE USED IN CALCULATING THE DERIVATIVES.
C BECAUSE FULL SCF CALCULATIONS ARE BEING DONE QUITE LARGE STEPS
C ARE NEEDED. ON THE OTHER HAND, THE STEP CANNOT BE VERY LARGE,
C AS THE SECOND DERIVITIVE IN FLEPO IS CALCULATED FROM THE
C DIFFERENCES OF TWO FIRST DERIVATIVES. CHANGE(1) IS FOR CHANGE IN
C BOND LENGTH, (2) FOR ANGLE, AND (3) FOR DIHEDRAL.
C
XDERIV(1)= 0.5D0/CHANGE(1)
XDERIV(2)= 0.5D0/CHANGE(2)
XDERIV(3)= 0.5D0/CHANGE(3)
ENDIF
DO 10 I=1,NVAR
10 XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
IF(NDEP.NE.0) CALL SYMTRY
CALL GMETRY(GEO,COORD)
C
C ESTABLISH THE ENERGY AT THE CURRENT POINT
C
CALL HCORE(COORD,H,W, WJ, WK, ENUCLR)
IF(NORBS*NELECS.GT.0)THEN
CALL ITER(H, W, WJ, WK, AA,.TRUE.,.FALSE.)
ELSE
AA=0.D0
ENDIF
LINEAR=(NORBS*(NORBS+1))/2
C
C RESTORE THE DENSITY MATRIX (WHY?)
C
DO 20 I=1,LINEAR
20 P(I)=PA(I)*2.D0
AA=(AA+ENUCLR)
IJ=0
DO 60 II=1,NUMAT
DO 50 IL=L1L,L1U
DO 50 JL=L2L,L2U
DO 50 KL=L3L,L3U
DO 30 LL=1,3
30 XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*JL+TVEC
1(LL,3)*KL
IJ=IJ+1
50 CONTINUE
60 CONTINUE
DO 90 I=1,NVAR
K=LOC(1,I)
L=LOC(2,I)
XSTORE=XPARAM(I)
DO 70 J=1,NVAR
70 GEO(LOC(2,J),LOC(1,J))=XPARAM(J)
GEO(L,K)=XSTORE-CHANGE(L)
IF(NDEP.NE.0) CALL SYMTRY
CALL GMETRY(GEO,COORD)
C
C IF NEEDED, CALCULATE "EXACT" DERIVITIVES.
C
CALL HCORE(COORD,H,W, WJ, WK,ENUCLR)
IF(NORBS*NELECS.GT.0)THEN
CALL ITER(H,W, WJ, WK,EE,.TRUE.,.FALSE.)
ELSE
EE=0.D0
ENDIF
DO 80 II=1,LINEAR
80 P(II)=PA(II)*2.D0
EE=(EE+ENUCLR)
ERRFN(I)=(AA-EE)*23.061D0*XDERIV(L)*2.D0
90 CONTINUE
IF(DEBUG)THEN
WRITE(6,'('' ERROR FUNCTION'')')
WRITE(6,'(10F8.3)')(ERRFN(I),I=1,NVAR)
ENDIF
RETURN
END
|
