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