1!
2!  Dalton, a molecular electronic structure program
3!  Copyright (C) by the authors of Dalton.
4!
5!  This program is free software; you can redistribute it and/or
6!  modify it under the terms of the GNU Lesser General Public
7!  License version 2.1 as published by the Free Software Foundation.
8!
9!  This program is distributed in the hope that it will be useful,
10!  but WITHOUT ANY WARRANTY; without even the implied warranty of
11!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12!  Lesser General Public License for more details.
13!
14!  If a copy of the GNU LGPL v2.1 was not distributed with this
15!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
16!
17!
18c /* deck cclr_fa */
19*=====================================================================*
20       SUBROUTINE CCLR_FA ( LABELA, ISYMTA,  ! inp: label/symmetry A
21     &                      LISTB,  ITAMPB,  ! inp: B resp. amplit.
22     &                      LISTC,  IZETVC,  ! inp: C resp. zeta vec.
23     &                      XINT,            ! inp: integrals of operator
24     &                      WORK,   LWORK   )! work space
25*---------------------------------------------------------------------*
26*
27*    Purpose: transformation of a response vector with a F matrix
28*             where the hamiltonian has been substituted by a
29*             perturbation operator
30*
31*             F^C{A} * t^B = <lambda^C|[[A,t^B],tau_nu]|CC>
32*
33*             JK+OC. CCMM: Allow for input of integrals if
34*             LABELA .eq. 'GIVE INT'
35*
36*    symmetries/variables:
37*
38*           ISYRES : result vector GAMMA1, GAMMA2
39*           ISYCTR : lagrangian multipliers (zeta vector) CTR1, CTR2
40*           ISYMTA : A perturbation
41*           ISYMTB : B response amplitudes
42*
43*    Note: the single and double excitation parts of the result GAMMA2
44*          are returned at the beginning of the work space in
45*          WORK(1)... WORK(NT1AM(ISYRES))
46*          WORK(NT1AM(ISYRES)+1)... WORK(NT1AM(ISYRES)+NT2AM(ISYRES))
47*          (double excitation part will be stored in packed form)
48*
49*     Written by Christof Haettig, October 1996.
50*
51*=====================================================================*
52#if defined (IMPLICIT_NONE)
53      IMPLICIT NONE
54#else
55#  include "implicit.h"
56#endif
57#include "priunit.h"
58#include "iratdef.h"
59#include "cbieri.h"
60#include "mxcent.h"
61#include "eribuf.h"
62#include "maxorb.h"
63#include "distcl.h"
64#include "ccorb.h"
65#include "ccisao.h"
66#include "ccsdsym.h"
67#include "ccsdinp.h"
68
69* local parameters:
70      CHARACTER MSGDBG*(17)
71      PARAMETER (MSGDBG='[debug] CCLR_FA> ')
72#if defined (SYS_CRAY)
73      REAL ONE, TWO, THREE
74#else
75      DOUBLE PRECISION ONE, TWO, THREE
76#endif
77      PARAMETER (ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0)
78      LOGICAL LOCDBG
79      PARAMETER (LOCDBG = .FALSE.)
80      INTEGER KDUM
81      PARAMETER (KDUM = +99 999 999) ! dummy address
82
83      CHARACTER*8 LABELA
84      CHARACTER*10 MODEL
85      CHARACTER LISTB*(*), LISTC*(*)
86      INTEGER ISYRES, ISYCTR, ISYMTA, ISYMTB, LWORK
87      INTEGER ITAMPB, IZETVC
88
89      INTEGER ISYTATB, ISYMJ, ISYMB, ISYMXY, ISYMI, ISYMA, IRREP, ISYM
90      INTEGER KBTAOO, KBTAVV, KPERTA, KT1AMPB, KT2AMPB, KCTR1, KCTR2
91      INTEGER KCTMO, KT1AMP0, KLAMDP0, KLAMDH0, KOFF1, KOFF2, KSCR
92      INTEGER KEND1, KEND2, KEND3, KEND1A, KXBMAT, KYBMAT, IERR
93      INTEGER LEND1, LEND2, LEND3, LEND1A, KGAMMA1, KGAMMA2, KEND0
94      INTEGER IOPT, MAXJ, NIJ, NJI, NAB, NBA, KEMAT1, KEMAT2
95
96#if defined (SYS_CRAY)
97      REAL FREQB
98      REAL SWAP
99      REAL WORK(LWORK), XINT(*)
100#else
101      DOUBLE PRECISION FREQB
102      DOUBLE PRECISION SWAP
103      DOUBLE PRECISION WORK(LWORK), XINT(*)
104#endif
105
106      INTEGER ILSTSYM
107
108*---------------------------------------------------------------------*
109* begin:
110*---------------------------------------------------------------------*
111      IF (CCSDT) THEN
112        WRITE (LUPRI,'(/1x,a)') 'F{A} matrix transformations not '
113     &          //'implemented for triples yet...'
114        CALL QUIT('Triples not implemented for G '//
115     &       'matrix transformations')
116      END IF
117
118      IF ( .not. (CCS .or. CC2 .or. CCSD) ) THEN
119        WRITE (LUPRI,'(/1x,a)') 'CCLR_FA called for a Coupled Cluster '
120     &          //'method not implemented in CCLR_FA...'
121        CALL QUIT('Unknown CC method in CCLR_FA.')
122      END IF
123
124*---------------------------------------------------------------------*
125* set & check symmetries:
126*---------------------------------------------------------------------*
127      ISYMTB  = ILSTSYM(LISTB,ITAMPB)   ! B
128      ISYCTR  = ILSTSYM(LISTC,IZETVC)   ! C
129      ISYTATB = MULD2H(ISYMTA,ISYMTB)   ! A x B
130      ISYRES  = MULD2H(ISYCTR,ISYTATB)  ! A x B x C
131
132      IF (LOCDBG) THEN
133        WRITE (LUPRI,*) 'LISTB,ITAMPB,ISYMTB:',LISTB,ITAMPB,ISYMTB
134        WRITE (LUPRI,*) 'LISTC,IZETVC,ISYCTR:',LISTC,IZETVC,ISYCTR
135      END IF
136
137
138      IF (ISYMOP .NE. 1) THEN
139        WRITE (LUPRI,'(/1x,a)') 'non-total-symmetric MO integrals?!... '
140     &          //'CCLR_G has never been debugged for that!...'
141      END IF
142
143      IF (MULD2H(ISYCTR,ISYTATB) .NE. ISYRES) THEN
144        CALL QUIT('Symmetry mismatch in CCLR_FA.')
145      END IF
146
147*---------------------------------------------------------------------*
148* flush print unit
149*---------------------------------------------------------------------*
150      Call FLSHFO(LUPRI)
151
152      IF (LOCDBG) THEN
153        WRITE (LUPRI,'(/1x,a,i15)') 'work space in CCLR_FA:',LWORK
154      END IF
155*---------------------------------------------------------------------*
156* initialize pointer for work space and allocate memory for
157*  1) single excitation part of the result vector
158*  2) one-index transformed perturbation integrals A^B (occ/occ block)
159*  3) one-index transformed perturbation integrals A^B (vir/vir block)
160*  4) perturbation integrals A
161*  5) singles part of response amplitudes T1^B
162*  6) singles part of zeroth order lagrangian multipliers
163*---------------------------------------------------------------------*
164      KGAMMA1 = 1
165      KBTAOO  = KGAMMA1 + NT1AM(ISYRES)
166      KBTAVV  = KBTAOO  + NMATIJ(ISYTATB)
167      KPERTA  = KBTAVV  + NMATAB(ISYTATB)
168      KT1AMPB = KPERTA  + NT1AM(ISYMTA)
169      KCTR1   = KT1AMPB + NT1AM(ISYMTB)
170      KEND1   = KCTR1   + NT1AM(ISYCTR)
171      LEND1   = LWORK - KEND1
172
173      IF (LEND1 .LT. 0) THEN
174        CALL QUIT('Insufficient work space in CCLR_FA.')
175      END IF
176
177*---------------------------------------------------------------------*
178* initialize single excitation part of result vector GAMMA1:
179*---------------------------------------------------------------------*
180      Call DZERO (WORK(KGAMMA1), NT1AM(ISYRES))
181
182*---------------------------------------------------------------------*
183* for CCS and zeroth-order zeta vector all contributions vanish:
184*---------------------------------------------------------------------*
185      IF (CCS .AND. LISTC(1:2).EQ.'L0') RETURN
186
187*---------------------------------------------------------------------*
188* read singles parts for B response amplitudes and zeta vector:
189*---------------------------------------------------------------------*
190      IOPT = 1
191      CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL,
192     &                  WORK(KT1AMPB),WORK(KDUM)  )
193
194
195      IOPT = 1
196      Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL,
197     &                  WORK(KCTR1),WORK(KDUM))
198
199      IF (LOCDBG) THEN
200        CAll AROUND('response T amplitudes B:')
201        WRITE (LUPRI,*) 'LIST/INDEX:',LISTB,ITAMPB
202        WRITE (LUPRI,*) 'Symmetry:      ',ISYMTB
203        CAll CC_PRP(WORK(KT1AMPB),WORK(KDUM),ISYMTB,1,0)
204        CALL AROUND('CC lagrange multipliers')
205        CALL CC_PRP(WORK(KCTR1), WORK(KDUM),  ISYCTR, 1, 0)
206      END IF
207
208*---------------------------------------------------------------------*
209* read & resort one-electron integrals for operator A:
210*---------------------------------------------------------------------*
211      KCTMO   = KEND1
212      KT1AMP0 = KCTMO   + N2BST(ISYMTA)
213      KLAMDP0 = KT1AMP0 + NT1AM(ISYMOP)
214      KLAMDH0 = KLAMDP0 + NLAMDT
215      KEND1A  = KLAMDH0 + NLAMDT
216      LEND1A  = LWORK - KEND1A
217
218      IF (LEND1A .LT. 0) THEN
219        CALL QUIT('Insufficient work space in CCLR_FA.')
220      END IF
221C
222C     JK+OC, CCSLV
223      IF (LABELA.EQ.'GIVE INT') THEN
224        CALL DCOPY(N2BST(ISYMTA),XINT(1),1,WORK(KCTMO),1)
225      ELSE
226* read the AO integrals:
227       CALL CCPRPAO(LABELA,.TRUE.,WORK(KCTMO),IRREP,ISYM,IERR,
228     &             WORK(KEND1A),LEND1A)
229       IF (IERR.NE.0 .OR. IRREP.NE.ISYMTA) THEN
230        CALL QUIT('CCLR_FA: error while reading operator '//LABELA)
231       END IF
232      END IF
233C
234* get MO coefficients:
235      CALL DZERO(WORK(KT1AMP0),NT1AMX)
236      CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0),
237     &            WORK(KEND1A),LEND1A)
238
239* transform one-electron integrals in place:
240      CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP0),WORK(KLAMDH0),
241     &              WORK(KEND1A),LEND1A,ISYMTA,1,1)
242
243* resort occupied/virtual block to T1 like storage:
244      CALL DZERO(WORK(KPERTA),NT1AM(ISYMTA))
245      DO ISYMJ = 1, NSYM
246        ISYMB = MULD2H(ISYMJ,ISYMTA)
247
248        DO J = 1, NRHF(ISYMJ)
249        DO B = 1, NVIR(ISYMB)
250          KOFF1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) + B
251          KOFF2 = IFCVIR(ISYMJ,ISYMB) + NORB(ISYMJ)*(B-1) + J
252
253          WORK(KPERTA-1+KOFF1) = WORK(KCTMO-1+KOFF2)
254        END DO
255        END DO
256      END DO
257
258      IF (LOCDBG) THEN
259        WRITE (LUPRI,*) MSGDBG,' A integrals in MO basis:'
260        WRITE (LUPRI,*) MSGDBG,' label, symmetry:',LABELA,ISYMTA
261        Call CC_PRP(WORK(KPERTA),WORK(KDUM),ISYMTA,1,0)
262      END IF
263
264*---------------------------------------------------------------------*
265* calculate A perturbation integrals one-index transformed with
266* the B response amplitudes T1^B:
267*---------------------------------------------------------------------*
268* occ/occ block:
269      Call CCG_1ITROO(WORK(KBTAOO), ISYTATB,
270     &                WORK(KPERTA), ISYMTA,
271     &                WORK(KT1AMPB),ISYMTB  )
272
273* vir/vir block:
274      Call CCG_1ITRVV(WORK(KBTAVV), ISYTATB,
275     &                WORK(KPERTA), ISYMTA,
276     &                WORK(KT1AMPB),ISYMTB  )
277
278*=====================================================================*
279*   CCS part:  < Zeta_1 | [tA^B, tau_1] | HF>
280*=====================================================================*
281* do one-index transformation with Zeta vector:
282      IOPT  = 2
283      Call CCG_1ITRVO(WORK(KGAMMA1),ISYRES,WORK(KBTAOO),WORK(KBTAVV),
284     &                ISYTATB,WORK(KCTR1),ISYCTR,IOPT          )
285
286      IF (LOCDBG) THEN
287        WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):'
288        WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB))
289        WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):'
290        WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB))
291        WRITE (LUPRI,*) MSGDBG,
292     *            'contrib. of one-index trans. A to GAMMA:'
293        Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0)
294      END IF
295
296*---------------------------------------------------------------------*
297* end of CCS part
298*---------------------------------------------------------------------*
299
300      If (CCS) Return
301
302*=====================================================================*
303* CC2/CCSD part for the singles: <Zeta_2| [[A,T2B], tau_1] |HF>
304*=====================================================================*
305
306*---------------------------------------------------------------------*
307* memory allocation:
308* 1) double excitation part of response amplitudes T2B (packed)
309* 2) double excitation part of zeta vector (squared)
310* 3) double excitation part of zeta vector (packed)
311* N.B. we account here for the fact, that the packed double excitation
312* part of the result vector will be returned at the beginning of the
313* work space, so we make sure, that there is enough space before
314* the zeta vector to store there later on GAMMA2
315*---------------------------------------------------------------------*
316       KT2AMPB = KEND1
317       KCTR2   = KT2AMPB + MAX( NT2AM(ISYMTB), NT2AM(ISYRES) )
318       KEND2   = KCTR2 + NT2SQ(ISYCTR)
319       LEND2   = LWORK - KEND2
320
321       IF (LEND2 .LT. NT2AM(ISYCTR) ) THEN
322         CALL QUIT('Insufficient work space in CCLR_FA.')
323       END IF
324
325*---------------------------------------------------------------------*
326* read response amplitudes T2B and scale the diagonal:
327*---------------------------------------------------------------------*
328      IOPT = 2
329      CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL,
330     &                  WORK(KDUM),WORK(KT2AMPB)  )
331
332      CAll CCLR_DIASCL(WORK(KT2AMPB),TWO,ISYMTB)
333
334      IF (LOCDBG) THEN
335        WRITE (LUPRI,*) MSGDBG, 'B response amplitudes:'
336        Call CC_PRP(WORK(KT1AMPB),WORK(KT2AMPB),ISYMTB,1,1)
337      END IF
338
339*---------------------------------------------------------------------*
340* read packed lagrangian multipliers and square them up:
341*---------------------------------------------------------------------*
342      IOPT = 2
343      Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL,
344     &                  WORK(KDUM),WORK(KEND2))
345
346      CALL CC_T2SQ (WORK(KEND2), WORK(KCTR2), ISYCTR)
347
348*---------------------------------------------------------------------*
349* calculate X^B and Y^B intermediates:
350*---------------------------------------------------------------------*
351      ISYMXY  = MULD2H(ISYCTR,ISYMTB)
352
353      KXBMAT  = KEND2
354      KYBMAT  = KXBMAT  + NMATIJ(ISYMXY)
355      KSCR    = KYBMAT  + NMATAB(ISYMXY)
356      KEND3   = KSCR    + NT1AM(ISYRES)
357      LEND3   = LWORK - KEND3
358
359      If (LEND3 .LT. 0) THEN
360        CALL QUIT('Insufficient work space in CCLR_FA.')
361      END IF
362
363* calculate X^C & Y^C intermediate:
364      Call CC_XI(WORK(KXBMAT),WORK(KCTR2), ISYCTR,
365     &           WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3)
366
367      Call CC_YI(WORK(KYBMAT),WORK(KCTR2), ISYCTR,
368     &           WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3)
369
370      IF (LOCDBG) THEN
371        WRITE (LUPRI,*) MSGDBG,'response X intermediate:'
372        WRITE (LUPRI,'(5f12.6)') (WORK(KXBMAT-1+I),I=1,NMATIJ(ISYMXY))
373        WRITE (LUPRI,*) MSGDBG,'response Y intermediate:'
374        WRITE (LUPRI,'(2f12.6)') (WORK(KYBMAT-1+I),I=1,NMATAB(ISYMXY))
375      END IF
376
377* calculate XY^B:  XY^B_ij = X^B_ji,  XY^B_bd = -Y^B_bd
378* 1.) transpose X^B intermediate
379      DO ISYMI = 1, NSYM
380        ISYMJ = MULD2H(ISYMI,ISYMXY)
381        IF (ISYMJ .LE. ISYMI) THEN
382          DO I = 1, NRHF(ISYMI)
383            MAXJ =  NRHF(ISYMJ)
384            IF (ISYMJ .EQ. ISYMI) MAXJ = I-1
385          DO J = 1, MAXJ
386            NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I
387            NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J
388            SWAP = WORK(KXBMAT-1+NIJ)
389            WORK(KXBMAT-1+NIJ) = WORK(KXBMAT-1+NJI)
390            WORK(KXBMAT-1+NJI) = SWAP
391          END DO
392          END DO
393        END IF
394      END DO
395
396* 2.) multiply Y^B intermediate with -1:
397      Call DSCAL(NMATAB(ISYMXY), -ONE, WORK(KYBMAT), 1)
398
399
400* do one-index transformation of XY^B with A integrals:
401      IOPT  = 2
402      Call CCG_1ITRVO(WORK(KSCR),ISYRES,
403     &                WORK(KXBMAT),WORK(KYBMAT),ISYMXY,
404     &                WORK(KPERTA),ISYMTA,      IOPT    )
405
406* add contribution to GAMMA1:
407      Call DAXPY (NT1AM(ISYRES), ONE, WORK(KSCR),1, WORK(KGAMMA1), 1)
408
409      IF (LOCDBG) THEN
410        WRITE (LUPRI,*) MSGDBG,'A integrals:'
411        WRITE (LUPRI,'(5f12.6)') (WORK(KPERTA-1+I),I=1,NT1AM(ISYMTA))
412        WRITE (LUPRI,*) MSGDBG, 'CC2/CCSD contribution to singles part:'
413        Call CC_PRP(WORK(KSCR),WORK(KDUM),ISYRES,1,0)
414        WRITE (LUPRI,*) MSGDBG, 'GAMMA1 now:'
415        Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0)
416      END IF
417
418*---------------------------------------------------------------------*
419
420*=====================================================================*
421* CC2/CCSD part for the doubles: <Zeta_2| [[A,T1B], tau_2] |HF>
422*=====================================================================*
423
424*---------------------------------------------------------------------*
425* reorganize work space, so that the result vector GAMMA2 can be
426* stored at the early beginning of the work space
427*---------------------------------------------------------------------*
428      KGAMMA1 = 1
429      KGAMMA2 = KGAMMA1 + NT1AM(ISYRES)
430      KEND0   = KGAMMA2 + NT2AM(ISYRES)
431
432      IF (KEND0 .GT. KCTR2) THEN
433        CALL QUIT('memory organization mixed up in CCLR_FA.')
434      END IF
435
436      KEMAT1 = KCTR2  + NT2SQ(ISYCTR)
437      KEMAT2 = KEMAT1 + NMATAB(ISYTATB)
438      KEND3  = KEMAT2 + NMATIJ(ISYTATB)
439      LEND3  = LWORK - KEND3
440
441      IF ( LEND3 .LT. 0 ) THEN
442        CALL QUIT('Insufficient work space in CCLR_FA.')
443      END IF
444
445*---------------------------------------------------------------------*
446* transpose tA^B(a b) --> EMAT1(b a)
447*---------------------------------------------------------------------*
448      DO ISYMA = 1, NSYM
449        ISYMB = MULD2H(ISYMA,ISYTATB)
450        DO A = 1, NVIR(ISYMA)
451        DO B = 1, NVIR(ISYMB)
452          NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B-1) + A
453          NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A-1) + B
454
455          WORK(KEMAT1 - 1 + NBA) = WORK(KBTAVV - 1 + NAB)
456        END DO
457        END DO
458      END DO
459
460
461*---------------------------------------------------------------------*
462* transpose tA^B(i j) --> EMAT2(j i)
463*---------------------------------------------------------------------*
464      DO ISYMI = 1, NSYM
465        ISYMJ = MULD2H(ISYMI,ISYTATB)
466        DO I = 1, NRHF(ISYMI)
467        DO J = 1, NRHF(ISYMJ)
468          NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I
469          NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J
470
471          WORK(KEMAT2 - 1 + NJI) = WORK(KBTAOO - 1 + NIJ)
472        END DO
473        END DO
474      END DO
475
476      IF (LOCDBG) THEN
477        WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):'
478        WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB))
479        WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):'
480        WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB))
481        WRITE (LUPRI,*) MSGDBG,'EMAT2:'
482        WRITE (LUPRI,'(5f12.6)') (WORK(KEMAT2-1+I),I=1,NMATIJ(ISYTATB))
483        WRITE (LUPRI,*) MSGDBG,'EMAT1:'
484        WRITE (LUPRI,'(2f12.6)') (WORK(KEMAT1-1+I),I=1,NMATAB(ISYTATB))
485      END IF
486
487*---------------------------------------------------------------------*
488* combine EMAT1/EMAT2 with lagrangian multipliers:
489* (note: this overwrites the intermedites stored at the beginning
490*        of the work space...)
491*---------------------------------------------------------------------*
492* initialize GAMMA2:
493      CALL DZERO(WORK(KGAMMA2),NT2AM(ISYRES))
494
495* do the caculation:
496      Call CCRHS_E(WORK(KGAMMA2),WORK(KCTR2),WORK(KEMAT1),
497     &             WORK(KEMAT2), WORK(KEND3), LEND3, ISYCTR, ISYTATB)
498
499*---------------------------------------------------------------------*
500      IF (LOCDBG) THEN
501        WRITE (LUPRI,*) MSGDBG, 'GAMMA:'
502        Call CC_PRP(WORK(KGAMMA1),WORK(KGAMMA2),ISYRES,1,1)
503      END IF
504C
505      RETURN
506      END
507*=====================================================================*
508*                  END OF SUBROUTINE CCLR_FA                          *
509*=====================================================================*
510