1C  /* Deck cc_freeze_tripletexci */
2      SUBROUTINE CC_freeze_tripletexci(CAM1,CAMP,CAMM,ISYMTR,
3     &           MAXCORE, MAXION,
4     &           NRHFCORE,IRHFCORE,NVIRION,IVIRION,
5     &           LBOTH)
6C
7C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
8C     Purpose: Project out specific excitations
9C              from a trial vector (by zeroing
10C              specific elements)
11C     Ex1: zero all ai and aibj elements where i and j
12C     are valence orbitals (CORE-VALENCE SEPARATION)
13C     Ex2: zero all a an b elements that do not correspond
14C     to a specific virtual orbitals
15C     Sonia and Eirik 2016
16! Based on cc_pram3()
17! CAM is the vector analyzed, of symmetry ISYMTR
18! LBOTH checks if both CVS and IONISATION are requested
19! Control is passed via argument list, not via common block
20C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
21C
22      Implicit none
23C
24#include "priunit.h"
25#include "ccorb.h"
26#include "ccsdsym.h"
27#include "ccsdinp.h"
28C
29      Double precision CAM1(*), CAMP(*), CAMM(*)
30      Integer MAXCORE, NRHFCORE(8),IRHFCORE(MAXCORE,8)
31      Integer MAXION,NVIRION(8),IVIRION(MAXION,8)
32      integer ISYMTR,ISYMAI,ISYMI,ISYMA,ISYMJ,ISYMB,ISYMBJ
33      Double precision TWO, THR1, THR2, zero, thprt
34      PARAMETER (TWO = 2.0D0,zero=0.0d0, THPRT = 1.0D-9)
35      Logical LOCDBG, ikeep, LBOTH
36      Parameter (Locdbg = .false.)
37      Integer AA, II, MA, MI, JJ, BB, NBJ, NAI, MJ, MB
38      Integer KAIBJ, NAIBJ, INDEX
39      double precision c1nosq, C2MNOSQ, C2PNOSQ, cnosq
40      double precision pt1, ptm, ptp, sumofp, ddot
41      integer nl, n1,n2
42
43C
44      INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J
45C
46      CALL QENTER('CC_FREEZE_TRIPLETEXCI')
47C
48      CALL FLSHFO(LUPRI)
49C
50
51C---------------------------------------
52C     Loop through single excitation part.
53C---------------------------------------
54C
55!      WRITE(LUPRI,'(//A)')
56!     *     ' +=============================================='
57!     *    //'===============================+'
58!      WRITE(LUPRI,'(1X,A)')
59!     *     '| symmetry|  orbital index  |   Excitation Numbers'
60!     *     //'             |   Amplitude  |'
61!      WRITE(LUPRI,'(1X,A)')
62!     *     '|  Index  |   a   b   i   j |      NAI      NBJ |'
63!     *     //'     NAIBJ    |              |'
64!      WRITE(LUPRI,'(A)')
65!     *     ' +=============================================='
66!     *    //'===============================+'
67C
68      ISYMAI = MULD2H(ISYMTR,ISYMOP)
69C
70      N1 = 0
71      DO 100 ISYMA = 1,NSYM
72         ISYMI = MULD2H(ISYMAI,ISYMA)
73         DO 110 I = 1,NRHF(ISYMI)
74            MI = IORB(ISYMI) + I
75            DO 120 A=1,NVIR(ISYMA)
76               NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A
77               MA = IORB(ISYMA) + NRHF(ISYMA) +  A
78               ikeep = .false.
79               IF (LBOTH) THEN
80                 do ii = 1, NRHFCORE(ISYMI)
81                  IF (I==IRHFCORE(II,ISYMI)) THEN
82                     do aa = 1, NVIRION(ISYMA)
83                        IF (A==IVIRION(AA,ISYMA)) THEN
84                           ikeep = .true.
85                           go to 333
86                        END IF
87                     end do
88                  END IF
89                 end do
90               ELSE
91                 do ii = 1, NRHFCORE(ISYMI)
92                   IF (I==IRHFCORE(II,ISYMI)) THEN
93                      ikeep = .true.
94                      go to 333
95                   end if
96                 end do
97                 do aa = 1, NVIRION(ISYMA)
98                   IF (A==IVIRION(AA,ISYMA)) THEN
99                    ikeep = .true.
100                    go to 333
101                   END IF
102                 end do
103               end if
104  333          continue
105               if (.not.ikeep) CAM1(NAI) = zero
106C
107  120       CONTINUE
108  110    CONTINUE
109  100 CONTINUE
110C
111C
112      CALL FLSHFO(LUPRI)
113C--------------------------------------------
114C     Loop through Doublee excitation vector.
115C     If not ccs or ccp2
116C--------------------------------------------
117C
118      IF (.NOT. ( CCS .OR. CCP2 )) THEN
119C
120!      WRITE(LUPRI,'(A)')
121!     *     ' +----------------------------------------------'
122!     *    //'-------------------------------+'
123C
124
125      N2 = 0
126C
127      ikeep = .false.
128      DO 200 ISYMAI = 1,NSYM
129         ISYMBJ = MULD2H(ISYMAI,ISYMTR)
130         IF (ISYMAI.LT.ISYMBJ) CYCLE
131         DO 210 ISYMJ = 1,NSYM
132            ISYMB = MULD2H(ISYMJ,ISYMBJ)
133            DO 220 ISYMI = 1,NSYM
134               ISYMA = MULD2H(ISYMI,ISYMAI)
135               DO 230 J = 1,NRHF(ISYMJ)
136                  MJ = IORB(ISYMJ) + J
137                  DO 240 B = 1,NVIR(ISYMB)
138                     NBJ = IT1AM(ISYMB,ISYMJ)
139     *                   + NVIR(ISYMB)*(J - 1) + B
140                     MB = IORB(ISYMB) + NRHF(ISYMB) + B
141                     DO 250 I = 1,NRHF(ISYMI)
142                        MI = IORB(ISYMI) + I
143                        DO 260 A = 1,NVIR(ISYMA)
144                           NAI = IT1AM(ISYMA,ISYMI)
145     *                         + NVIR(ISYMA)*(I - 1) + A
146                           MA = IORB(ISYMA) + NRHF(ISYMA) +  A
147                           IF ((ISYMAI.EQ.ISYMBJ).AND.
148     *                         (NAI .LT. NBJ))
149     *                          GOTO 260
150                           IF (ISYMAI.EQ.ISYMBJ) THEN
151                              NAIBJ = IT2AM(ISYMAI,ISYMBJ)
152     *                             + INDEX(NAI,NBJ)
153                           ELSE
154                               NAIBJ = IT2AM(ISYMAI,ISYMBJ)
155     *                           + NT1AM(ISYMBJ)*(NAI-1)+NBJ
156                           ENDIF
157!***
158                           ikeep = .false.
159                           if (LBOTH) then
160                              do ii = 1, nrhfcore(isymi)
161                                 if (i==IRHFCORE(II,ISYMI)) then
162                                   do aa = 1, NVIRION(ISYMA)
163                                     IF (A.EQ.IVIRION(AA,ISYMA)) THEN
164                                        ikeep = .true.
165                                        go to 444
166                                     END IF
167                                   end do
168                                   do bb = 1, NVIRION(ISYMB)
169                                     IF (B==IVIRION(BB,ISYMB)) THEN
170                                        ikeep = .true.
171                                        go to 444
172                                     END IF
173                                   end do
174                                 end if
175                              end do
176                              do jj = 1, nrhfcore(isymj)
177                                 if (j==IRHFCORE(jj,ISYMJ)) then
178                                   do aa = 1, NVIRION(ISYMA)
179                                     IF (A.EQ.IVIRION(AA,ISYMA)) THEN
180                                        ikeep = .true.
181                                        go to 444
182                                     END IF
183                                   end do
184                                   do bb = 1, NVIRION(ISYMB)
185                                     IF (B==IVIRION(BB,ISYMB)) THEN
186                                        ikeep = .true.
187                                        go to 444
188                                     END IF
189                                   end do
190                                 end if
191                              end do
192                           else
193                             do ii = 1, nrhfcore(isymi)
194                              if (i==IRHFCORE(II,ISYMI)) then
195                                 ikeep = .true.
196                                 go to 444
197                              end if
198                             end do
199                             do aa = 1, nvirion(isyma)
200                              if (a==IVIRION(aa,ISYMA)) then
201                                 ikeep = .true.
202                                 go to 444
203                              end if
204                             end do
205                             do jj = 1, nrhfcore(isymj)
206                              if (j==IRHFCORE(JJ,ISYMJ)) then
207                                 ikeep = .true.
208                                 go to 444
209                              end if
210                             end do
211                             do bb = 1, nvirion(isymb)
212                              if (b==IVIRION(bb,ISYMB)) then
213                                 ikeep = .true.
214                                 go to 444
215                              end if
216                             end do
217                           end if
218!========================================
219  444                      continue
220                           if (.not.ikeep) CAMP(NAIBJ) = zero
221                           if (.not.ikeep) CAMM(NAIBJ) = zero
222  260                   CONTINUE
223  250                CONTINUE
224  240             CONTINUE
225  230          CONTINUE
226  220       CONTINUE
227  210    CONTINUE
228  200 CONTINUE
229C
230C
231      ENDIF
232C
233      CALL QEXIT('CC_FREEZE_TRIPLETEXCI')
234C
235      RETURN
236      END
237
238C  /* Deck cc_freeze_tripletstart */
239      SUBROUTINE CC_freeze_tripletstart(CAM1,CAMP,CAMM,ISYMTR,
240     &           MAXCORE, MAXION,
241     &           NRHFCORE,IRHFCORE,NVIRION,IVIRION,
242     &           LBOTH)
243C
244C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
245C     Filter dia for start vectors
246C     Purpose: allow automated selection of start vectors
247C              for core excitations and ionizations
248C              Achieved by setting all elements of the
249C              diagonal of the Fock matrix non refering to
250C              the selected core/diffuse orbital to a huge
251C              number so that it will be discarded by the
252C              FNDM3 routine later on (the routine picks up
253C              the indices of the lowest energy eigenvalue).
254C    In other words, we push the valence eigenvalues above the core
255C              specific elements)
256C
257C     Sonia and Eirik 2016
258C
259! Based on cc_pram3()
260! CAM is the vector analyzed, of symmetry ISYMTR
261! LBOTH checks if both CVS and IONISATION are requested
262! Control is passed via argument list, not via common block
263C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
264C
265      Implicit none
266C
267#include "priunit.h"
268#include "ccorb.h"
269#include "ccsdsym.h"
270#include "ccsdinp.h"
271C
272      Double precision CAM1(*), CAMP(*), CAMM(*)
273      Integer MAXCORE, NRHFCORE(8),IRHFCORE(MAXCORE,8)
274      Integer MAXION,NVIRION(8),IVIRION(MAXION,8)
275      integer ISYMTR,ISYMAI,ISYMI,ISYMA,ISYMJ,ISYMB,ISYMBJ
276      Double precision TWO, THR1, THR2, crazy, thprt
277      PARAMETER (TWO = 2.0D0,crazy=1.0d9, THPRT = 1.0D-9)
278      Logical LOCDBG, ikeep, LBOTH
279      Parameter (Locdbg = .false.)
280      Integer AA, II, MA, MI, JJ, BB, NBJ, NAI, MJ, MB
281      Integer KAIBJ, NAIBJ, INDEX
282      double precision c1nosq, C2MNOSQ, C2PNOSQ, cnosq
283      double precision pt1, ptm, ptp, sumofp, ddot
284      integer nl, n1,n2
285
286C
287      INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J
288C
289      CALL QENTER('CC_FREEZE_TRIPLETSTART')
290C
291C      CALL FLSHFO(LUPRI)
292C
293      ISYMAI = MULD2H(ISYMTR,ISYMOP)
294C
295      N1 = 0
296      DO 100 ISYMA = 1,NSYM
297         ISYMI = MULD2H(ISYMAI,ISYMA)
298         DO 110 I = 1,NRHF(ISYMI)
299            MI = IORB(ISYMI) + I
300            DO 120 A=1,NVIR(ISYMA)
301               NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A
302               MA = IORB(ISYMA) + NRHF(ISYMA) +  A
303               ikeep = .false.
304               IF (LBOTH) THEN
305                 do ii = 1, NRHFCORE(ISYMI)
306                  IF (I==IRHFCORE(II,ISYMI)) THEN
307                     do aa = 1, NVIRION(ISYMA)
308                        IF (A==IVIRION(AA,ISYMA)) THEN
309                           ikeep = .true.
310                           go to 333
311                        END IF
312                     end do
313                  END IF
314                 end do
315               ELSE
316                 do ii = 1, NRHFCORE(ISYMI)
317                   IF (I==IRHFCORE(II,ISYMI)) THEN
318                      ikeep = .true.
319                      go to 333
320                   end if
321                 end do
322                 do aa = 1, NVIRION(ISYMA)
323                   IF (A==IVIRION(AA,ISYMA)) THEN
324                    ikeep = .true.
325                    go to 333
326                   END IF
327                 end do
328               end if
329  333          continue
330               if (.not.ikeep) CAM1(NAI) = crazy
331C
332  120       CONTINUE
333  110    CONTINUE
334  100 CONTINUE
335C
336C
337      CALL FLSHFO(LUPRI)
338C
339C--------------------------------------------
340C     Loop through Doublee excitation vector.
341C     If not ccs or ccp2
342C--------------------------------------------
343C
344      IF (.NOT. ( CCS .OR. CCP2 )) THEN
345C
346      N2 = 0
347C
348      ikeep = .false.
349      DO 200 ISYMAI = 1,NSYM
350         ISYMBJ = MULD2H(ISYMAI,ISYMTR)
351         IF (ISYMAI.LT.ISYMBJ) CYCLE
352         DO 210 ISYMJ = 1,NSYM
353            ISYMB = MULD2H(ISYMJ,ISYMBJ)
354            DO 220 ISYMI = 1,NSYM
355               ISYMA = MULD2H(ISYMI,ISYMAI)
356               DO 230 J = 1,NRHF(ISYMJ)
357                  MJ = IORB(ISYMJ) + J
358                  DO 240 B = 1,NVIR(ISYMB)
359                     NBJ = IT1AM(ISYMB,ISYMJ)
360     *                   + NVIR(ISYMB)*(J - 1) + B
361                     MB = IORB(ISYMB) + NRHF(ISYMB) + B
362                     DO 250 I = 1,NRHF(ISYMI)
363                        MI = IORB(ISYMI) + I
364                        DO 260 A = 1,NVIR(ISYMA)
365                           NAI = IT1AM(ISYMA,ISYMI)
366     *                         + NVIR(ISYMA)*(I - 1) + A
367                           MA = IORB(ISYMA) + NRHF(ISYMA) +  A
368                           IF ((ISYMAI.EQ.ISYMBJ).AND.
369     *                         (NAI .LT. NBJ))
370     *                          GOTO 260
371                           IF (ISYMAI.EQ.ISYMBJ) THEN
372                              NAIBJ = IT2AM(ISYMAI,ISYMBJ)
373     *                             + INDEX(NAI,NBJ)
374                           ELSE
375                               NAIBJ = IT2AM(ISYMAI,ISYMBJ)
376     *                           + NT1AM(ISYMBJ)*(NAI-1)+NBJ
377                           ENDIF
378                           ikeep = .false.
379                           if (LBOTH) then
380                              do ii = 1, nrhfcore(isymi)
381                                 if (i==IRHFCORE(II,ISYMI)) then
382                                   do aa = 1, NVIRION(ISYMA)
383                                     IF (A.EQ.IVIRION(AA,ISYMA)) THEN
384                                        ikeep = .true.
385                                        go to 444
386                                     END IF
387                                   end do
388                                   do bb = 1, NVIRION(ISYMB)
389                                     IF (B==IVIRION(BB,ISYMB)) THEN
390                                        ikeep = .true.
391                                        go to 444
392                                     END IF
393                                   end do
394                                 end if
395                              end do
396                              do jj = 1, nrhfcore(isymj)
397                                 if (j==IRHFCORE(jj,ISYMJ)) then
398                                   do aa = 1, NVIRION(ISYMA)
399                                     IF (A.EQ.IVIRION(AA,ISYMA)) THEN
400                                        ikeep = .true.
401                                        go to 444
402                                     END IF
403                                   end do
404                                   do bb = 1, NVIRION(ISYMB)
405                                     IF (B==IVIRION(BB,ISYMB)) THEN
406                                        ikeep = .true.
407                                        go to 444
408                                     END IF
409                                   end do
410                                 end if
411                              end do
412                           else
413                             do ii = 1, nrhfcore(isymi)
414                              if (i==IRHFCORE(II,ISYMI)) then
415                                 ikeep = .true.
416                                 go to 444
417                              end if
418                             end do
419                             do aa = 1, nvirion(isyma)
420                              if (a==IVIRION(aa,ISYMA)) then
421                                 ikeep = .true.
422                                 go to 444
423                              end if
424                             end do
425                             do jj = 1, nrhfcore(isymj)
426                              if (j==IRHFCORE(JJ,ISYMJ)) then
427                                 ikeep = .true.
428                                 go to 444
429                              end if
430                             end do
431                             do bb = 1, nvirion(isymb)
432                              if (b==IVIRION(bb,ISYMB)) then
433                                 ikeep = .true.
434                                 go to 444
435                              end if
436                             end do
437                           end if
438!========================================
439  444                      continue
440                           if (.not.ikeep) CAMP(NAIBJ) = crazy
441                           if (.not.ikeep) CAMM(NAIBJ) = crazy
442  260                   CONTINUE
443  250                CONTINUE
444  240             CONTINUE
445  230          CONTINUE
446  220       CONTINUE
447  210    CONTINUE
448  200 CONTINUE
449C
450C
451      ENDIF
452C
453C
454      CALL QEXIT('CC_FREEZE_TRIPLETSTART')
455C
456      RETURN
457      END
458