1      SUBROUTINE COMPFG(XPARAM,INT,ESCF,FULSCF,GRAD,LGRAD)
2      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3      INCLUDE 'SIZES'
4      DIMENSION XPARAM(MAXPAR),GRAD(MAXPAR)
5      LOGICAL LGRAD, FULSCF, LIMSCF
6      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY,DUMY(MAXPAR)
7      COMMON /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR)
8      COMMON /GEOM  / GEO(3,NUMATM), XCOORD(3,NUMATM)
9      COMMON /ATHEAT/ ATHEAT
10      COMMON /WMATRX/ WJ(N2ELEC), WK(N2ELEC)
11      COMMON /ENUCLR/ ENUCLR
12      COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
13      COMMON /ELECT / ELECT
14      PARAMETER (MDUMY=MAXPAR**2-MPACK)
15      COMMON /SCRACH/ RXYZ(MPACK), XDUMY(MDUMY)
16      COMMON /HMATRX/ H(MPACK)
17      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
18     1                NA(NUMATM), NB(NUMATM), NC(NUMATM)
19      COMMON /ERRFN  / ERRFN(MAXPAR), AICORR(MAXPAR)
20      COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CBETA(MORB2),EIGB(MAXORB)
21      COMMON /LAST  / LAST
22      COMMON /NUMCAL/ NUMCAL
23      COMMON /SCFTYP/ EMIN, LIMSCF
24      COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
25     1       /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
26     2                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
27     3                NCLOSE,NOPEN,NDUMY,FRACT
28C COSMO change A. Klamt
29      LOGICAL ISEPS, USEPS , UPDA
30      COMMON /ISEPS/  ISEPS, USEPS, UPDA
31C end of COSMO change
32C***********************************************************************
33C
34C   COMPFG CALCULATES (A) THE HEAT OF FORMATION OF THE SYSTEM, AND
35C                     (B) THE GRADIENTS, IF LGRAD IS .TRUE.
36C
37C   ON INPUT  XPARAM = ARRAY OF PARAMETERS TO BE USED IN INTERNAL COORDS
38C             LGRAD  = .TRUE. IF GRADIENTS ARE NEEDED, .FALSE. OTHERWISE
39C             INT    = .TRUE. IF HEAT OF FORMATION IS TO BE CALCULATED
40C             FULSCF = .TRUE. IF FULL SCF TO BE DONE, .FALSE. OTHERWISE.
41C
42C   ON OUTPUT ESCF  = HEAT OF FORMATION.
43C             GRAD   = ARRAY OF GRADIENTS, IF LGRAD = .TRUE.
44C
45C***********************************************************************
46      COMMON /KEYWRD/KEYWRD
47      CHARACTER*241 KEYWRD
48      LOGICAL DEBUG, INT, PRINT, ANALYT, LARGE, USEDCI,
49     1FORCE, TIMES, AIDER
50      DIMENSION COORD(3,NUMATM), W(N2ELEC), DEGREE(3), XPAREF(MAXPAR)
51     1,DELTAP(NMECI**2) ,DELTA(NMECI*MAXORB)
52      SAVE DEGREE, PRINT, DEBUG
53      EQUIVALENCE (W,WJ)
54      DATA ICALCN /0/
55C                 MNDO     AM1      PM3      MINDO/
56      IF (ICALCN.NE.NUMCAL) THEN
57         ICALCN=NUMCAL
58         HTYPE(1)=6.1737D0
59         HTYPE(2)=3.3191D0
60         HTYPE(3)=7.1853D0
61         HTYPE(4)=1.7712D0
62         LTYPE=0
63         DO 30 I=1,NUMAT
64            IF(NAT(I).LT.99)THEN
65               DO 10 J=1,LTYPE
66   10          IF(NAT(I).EQ.MTYPE(J)) GOTO 20
67               LTYPE=LTYPE+1
68               MTYPE(LTYPE)=NAT(I)
69               NZTYPE(NAT(I))=LTYPE
70C
71C       LTYPE = NUMBER OF TYPES OF REAL ATOM PRESENT
72C       MTYPE = TYPES OF REAL ATOMS PRESENT
73               J=LTYPE
74   20          CONTINUE
75            ENDIF
76   30    CONTINUE
77         AIDER=(INDEX(KEYWRD,'AIDER').NE.0)
78         TIMES=(INDEX(KEYWRD,'TIMES').NE.0)
79         ANALYT=(INDEX(KEYWRD,'ANALYT').NE.0)
80         IF(INT.AND.ANALYT)CALL SETUPG
81         DEGREE(1)=1.D0
82         IF(INDEX(KEYWRD,' XYZ').NE.0)THEN
83            DEGREE(2)=1.D0
84         ELSE
85            DEGREE(2)=180.D0/3.141592652589D0
86         ENDIF
87         DEGREE(3)=DEGREE(2)
88         USEDCI=(NCLOSE.NE.NOPEN.AND.FRACT.NE.2.D0.AND.FRACT.NE.0.D0
89     1         .OR.(INDEX(KEYWRD,'C.I.').NE.0))
90         FORCE=(INDEX(KEYWRD,'FORCE').NE.0)
91         LARGE=(INDEX(KEYWRD,'LARGE') .NE. 0)
92         PRINT=(INDEX(KEYWRD,'COMPFG') .NE. 0)
93         DEBUG=(INDEX(KEYWRD,'DEBUG') .NE. 0 .AND. PRINT)
94         EMIN=0.D0
95         DO 40 I=1,NVAR
96   40    XPAREF(I)=XPARAM(I)
97      ENDIF
98C
99C SET UP COORDINATES FOR CURRENT CALCULATION
100C
101C       PLACE THE NEW VALUES OF THE VARIABLES IN THE ARRAY GEO.
102C       MAKE CHANGES IN THE GEOMETRY.
103      DO 50 I=1,NVAR
104         K=LOC(1,I)
105         L=LOC(2,I)
106   50 GEO(L,K)=XPARAM(I)
107C      IMPOSE THE SYMMETRY CONDITIONS + COMPUTE THE DEPENDENT-PARAMETERS
108      IF(NDEP.NE.0) CALL SYMTRY
109C      NOW COMPUTE THE ATOMIC COORDINATES.
110      IF( DEBUG ) THEN
111         IF( LARGE ) THEN
112            K=NATOMS
113         ELSE
114            K=MIN(5,NATOMS)
115         ENDIF
116         WRITE(6,FMT='('' INTERNAL COORDS'',/100(/,3F12.6))')
117     1            ((GEO(J,I)*DEGREE(J),J=1,3),I=1,K)
118      END IF
119      CALL GMETRY(GEO,COORD)
120      IF( DEBUG ) THEN
121         IF( LARGE ) THEN
122            K=NUMAT
123         ELSE
124            K=MIN(5,NUMAT)
125         ENDIF
126         WRITE(6,FMT='('' CARTESIAN COORDS'',/100(/,3F16.9))')
127     1            ((COORD(J,I),J=1,3),I=1,K)
128      ENDIF
129      IF(INT.AND.ANALYT)REWIND 2
130C COSMO change A. Klamt
131*     IF (.NOT. USEPS) THEN
132      IF (.NOT. ISEPS) THEN
133C end of COSMO change
134      IF(TIMES)CALL TIMER('BEFORE HCORE')
135      IF(INT)CALL HCORE(COORD, H, W, WJ, WK, ENUCLR)
136      IF(TIMES)CALL TIMER('AFTER HCORE')
137C
138C COMPUTE THE HEAT OF FORMATION.
139C
140      IF(NORBS.GT.0.AND.NELECS.GT.0) THEN
141         IF(TIMES)CALL TIMER('BEFORE ITER')
142         IF(INT) CALL ITER(H, W, WJ, WK, ELECT, FULSCF,.TRUE.)
143         IF(TIMES)CALL TIMER('AFTER ITER')
144      ELSE
145         ELECT=0.D0
146      ENDIF
147      ESCF=(ELECT+ENUCLR)*23.061D0+ATHEAT
148      IF(ESCF.LT.EMIN.OR.EMIN.EQ.0.D0)EMIN=ESCF
149      DO 60 I=1,NNHCO
150         CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4,I),ANGLE)
151         ESCF=ESCF+HTYPE(ITYPE)*SIN(ANGLE)**2
152   60 CONTINUE
153C COSMO change A. Klamt 18.7.91
154      ENDIF
155      IF (ISEPS) THEN
156        INDEPS=INDEX(KEYWRD,'EPS=')
157        CALL INITSV (INDEPS)
158C The following routine constructs the dielectric screening surface
159           CALL CONSTS (COORD)
160C The following routine constructs dielectric response matrix CCMAT
161        CALL BTOC (COORD)
162C A. Klamt 18.7.91
163        USEPS = .TRUE.
164        IF(TIMES) CALL TIMER('BEFORE HCORE')
165        IF(INT) CALL HCORE(COORD, H, W, WJ, WK, ENUCLR)
166        IF(TIMES) CALL TIMER('AFTER HCORE')
167C
168C COMPUTE THE HEAT OF FORMATION.
169C
170        IF(NORBS.GT.0.AND.NELECS.GT.0) THEN
171          IF(TIMES) CALL TIMER('BEFORE ITER')
172          IF(INT) CALL ITER(H, W, WJ, WK, ELECT, FULSCF,.TRUE.)
173          IF(TIMES) CALL TIMER('AFTER ITER')
174        ELSE
175          ELECT=0.D0
176        ENDIF
177        ESCF=(ELECT+ENUCLR)*23.061D0+ATHEAT
178        IF(ESCF.LT.EMIN.OR.EMIN.EQ.0.D0) EMIN=ESCF
179        DO 61 I=1,NNHCO
180         CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4,I),ANGLE)
181         ESCF=ESCF+HTYPE(ITYPE)*SIN(ANGLE)**2
182   61   CONTINUE
183      ENDIF
184C end of COSMO change
185C
186C FIND DERIVATIVES IF DESIRED
187C
188      IF(LGRAD) THEN
189         IF(TIMES)CALL TIMER('BEFORE DERIV')
190         CALL DERIV(GEO,GRAD)
191         IF(TIMES)CALL TIMER('AFTER DERIV')
192      ENDIF
193      IF(AIDER)THEN
194C
195C  ADD IN AB INITIO CORRECTION
196C
197         DO 70 I=1,NVAR
198   70    ESCF=ESCF+(XPARAM(I)-XPAREF(I))*AICORR(I)
199      ENDIF
200      IF(INT.AND.PRINT)
201     1WRITE(6,'(/10X,'' HEAT OF FORMATION'',G30.17)')ESCF
202      IF(PRINT.AND.LGRAD)
203     1   WRITE(6,FMT='('' GRADIENT       '',8F8.2,(/10F8.2))')
204     2                (GRAD(I),I=1,NVAR)
205C
206C REFORM DENSITY MATRIX, IF A C.I. DONE AND EITHER THE LAST SCF OR A
207C FORCE CALCULATION
208C
209      IF(USEDCI.AND. (LAST.EQ.1 .OR. FORCE))
210     1CALL MECIP(C,NORBS,DELTAP,DELTA)
211      RETURN
212      END
213