← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/tce/mrcc/tce_mrcc_heff.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tce/mrcc/tce_mrcc_heff.F
 
1
c
 
2
c ==============================================
 
3
c     Create effective Hamiltonian
 
4
c ==============================================
 
5
c
 
6
       subroutine tce_heff(d_em,k_e_offsetm,k_r1_offsetm,
 
7
     1 k_r2_offsetm,k_r3_offsetm,k_r4_offsetm,d_r1m,d_r2m,d_r3m,d_r4m,
 
8
     2 needt1,needt2,needr3act,needr4act,rtdb)
 
9
        implicit none
 
10
#include "tce.fh"
 
11
#include "mafdecls.fh"
 
12
#include "stdio.fh"
 
13
#include "rtdb.fh"
 
14
#include "errquit.fh"
 
15
#include "sym.fh"
 
16
#include "tce_mrcc.fh"
 
17
#include "global.fh"
 
18
#include "tce_main.fh"
 
19
 
 
20
       integer d_em(maxref)
 
21
       integer k_e_offsetm(maxref)
 
22
       integer iref,i
 
23
       double precision corr
 
24
       logical nodezero
 
25
       integer k_r1_offsetm(maxref)
 
26
       integer k_r2_offsetm(maxref)
 
27
       integer k_r3_offsetm(maxref)
 
28
       integer k_r4_offsetm(maxref)
 
29
       integer d_r1m(maxref),d_r2m(maxref)
 
30
       integer d_r3m(maxref),d_r4m(maxref)
 
31
       logical needt1,needt2,needr3act
 
32
       logical needr4act
 
33
       integer rtdb
 
34
 
 
35
       nodezero = (ga_nodeid().eq.0)
 
36
 
 
37
       if(lusesub)call ga_zero(g_heff)
 
38
 
 
39
       do i=1,nref*nref
 
40
         dbl_mb(k_heff+i-1) = 0.0d0
 
41
       enddo
 
42
 
 
43
       do iref=1,nref
 
44
 
 
45
        if((int_mb(k_refafi+iref-1).eq.int_mb(k_innodes
 
46
     1 +ga_nnodes()+ga_nodeid())).or.(.not.lusesub)) then
 
47
 
 
48
          call get_block(d_em(iref),corr,1,0)
 
49
          dbl_mb(k_heff+iref-1+(iref-1)*nref) = corr+duhfens(iref)
 
50
 
 
51
        if(lusesub) then
 
52
 
 
53
c          write(6,*)ga_nodeid(),corr+duhfens(iref),iref
 
54
          call ga_put(g_heff,nref*(iref-1)+iref,nref*(iref-1)+iref,1,1,
 
55
     1    corr+duhfens(iref),1)
 
56
        endif
 
57
 
 
58
        endif
 
59
 
 
60
       enddo
 
61
 
 
62
       call tce_heff_offdiagonal(k_r1_offsetm,k_r2_offsetm,
 
63
     1 k_r3_offsetm,k_r4_offsetm,d_r1m,d_r2m,d_r3m,d_r4m,needt1,
 
64
     2 needt2,needr3act,needr4act,rtdb) 
 
65
 
 
66
c       if(nodezero) then
 
67
c         call ma_print(dbl_mb(k_heff),nref,nref,'Heff')
 
68
c       endif
 
69
 
 
70
       return
 
71
       end
 
72
c
 
73
c ==============================================
 
74
c     Add offdiagonal elements of Heff
 
75
c ==============================================
 
76
c
 
77
       subroutine tce_heff_offdiagonal(k_r1_offsetm,
 
78
     1 k_r2_offsetm,k_r3_offsetm,k_r4_offsetm,d_r1m,d_r2m,d_r3m,
 
79
     2 d_r4m,needt1,needt2,needr3act,needr4act,rtdb)
 
80
        implicit none
 
81
#include "tce.fh"
 
82
#include "mafdecls.fh"
 
83
#include "stdio.fh"
 
84
#include "rtdb.fh"
 
85
#include "errquit.fh"
 
86
#include "sym.fh"
 
87
#include "tce_mrcc.fh"
 
88
#include "global.fh"
 
89
#include "tce_main.fh"
 
90
 
 
91
       integer rtdb
 
92
       logical nodezero
 
93
       integer k_r1_offsetm(maxref)
 
94
       integer k_r2_offsetm(maxref)
 
95
       integer k_r3_offsetm(maxref)
 
96
       integer k_r4_offsetm(maxref)
 
97
       integer d_r1m(maxref),d_r2m(maxref)
 
98
       integer d_r3m(maxref)
 
99
       integer d_r4m(maxref)
 
100
       integer iref
 
101
       integer i,j,p5b,h6b,k
 
102
       integer size,l,m,n,o
 
103
       integer l_r1,k_r1,l_r2,k_r2
 
104
       integer l_r3,k_r3
 
105
       integer l_r4,k_r4
 
106
       integer p1b,p2b,h3b,h4b
 
107
       integer p3b,p7b,p8b
 
108
       INTEGER p4b
 
109
       INTEGER p6b
 
110
       INTEGER h1b
 
111
       INTEGER h2b
 
112
       integer h1,h2,h3,p4,p5,p6
 
113
       integer h4,p7,p8
 
114
       integer orbindex(8)
 
115
       integer orbspin(8)
 
116
       integer ioccnew(maxorb,2)
 
117
       integer iu
 
118
       logical isfound
 
119
       logical needt1,needt2,needr3act
 
120
       logical needr4act
 
121
c       logical lusescffv
 
122
       integer is,iocc0(maxorb,2)
 
123
       integer i1,i2,dist
 
124
       double precision dsmult
 
125
c       logical limprovet
 
126
       EXTERNAL NXTASKsub
 
127
       EXTERNAL NXTASK
 
128
       INTEGER NXTASKsub
 
129
       INTEGER NXTASK
 
130
       INTEGER nxt
 
131
       INTEGER nprocs
 
132
       INTEGER count
 
133
 
 
134
       nodezero = (ga_nodeid().eq.0)
 
135
       isfound = .false.
 
136
 
 
137
c       if (.not.rtdb_get(rtdb,'mrcc:usescffermiv',mt_log,1,lusescffv))
 
138
c     1 lusescffv = .false.
 
139
c       if (.not.rtdb_get(rtdb,'mrcc:improvetiling',mt_log,1,limprovet))
 
140
c     1 limprovet = .false.
 
141
 
 
142
      do is=1,2
 
143
        do i=1,nmo(is)
 
144
          if(i.le.nocc(is)) then
 
145
            iocc0(i,is) = 1
 
146
          else
 
147
            iocc0(i,is) = 0
 
148
          end if
 
149
        end do
 
150
      end do
 
151
 
 
152
       do iref=1,nref
 
153
 
 
154
        if((int_mb(k_refafi+iref-1).eq.int_mb(k_innodes
 
155
     1 +ga_nnodes()+ga_nodeid())).or.(.not.lusesub)) then
 
156
 
 
157
         k_sym = k_symm(iref)
 
158
         k_offset = k_offsetm(iref)
 
159
         k_range = k_rangem(iref)
 
160
         k_spin = k_spinm(iref)
 
161
         k_movecs_sorted = k_movecs_sortedm(iref)
 
162
         k_active = k_active_tmpm(iref)
 
163
 
 
164
         noa = nblcks(1,iref)
 
165
         nob = nblcks(2,iref)
 
166
         nva = nblcks(3,iref)
 
167
         nvb = nblcks(4,iref)
 
168
 
 
169
         noab = noa+nob
 
170
         nvab = nva+nvb
 
171
c
 
172
c ---------------
 
173
c    R1 active
 
174
c ---------------
 
175
c
 
176
         do p5b = noab+1,noab+nvab
 
177
         do h6b = 1,noab
 
178
 
 
179
      k = 0
 
180
 
 
181
      if (int_mb(k_spin+p5b-1) .eq. int_mb(k_spin+h6b-1)) then
 
182
      if (ieor(int_mb(k_sym+p5b-1),int_mb(k_sym+h6b-1)).eq.irrep_t)then
 
183
      if ((.not.restricted).or.(int_mb(k_spin+p5b-1)+int_mb(k_spin+h6b-1
 
184
     &).ne.4)) then
 
185
      if(log_mb(k_isactive(iref)+p5b-1).and.
 
186
     &log_mb(k_isactive(iref)+h6b-1)) then
 
187
 
 
188
        size = int_mb(k_range+h6b-1) * int_mb(k_range+p5b-1)
 
189
 
 
190
        if (.not.ma_push_get(mt_dbl,size,'r1mi',l_r1,k_r1))
 
191
     1   call errquit('tce_mrcc_iface_r1: MA problem',0,MA_ERR)
 
192
 
 
193
        call get_hash_block(d_r1m(iref),dbl_mb(k_r1),size,
 
194
     1   int_mb(k_r1_offsetm(iref)),h6b-1+noab*(p5b-noab-1))
 
195
 
 
196
        k=0
 
197
        do i=1,int_mb(k_range+p5b-1)
 
198
        do j=1,int_mb(k_range+h6b-1)
 
199
          k = k + 1
 
200
 
 
201
        orbspin(1) = int_mb(k_spin+p5b-1)-1
 
202
        orbspin(2) = int_mb(k_spin+h6b-1)-1
 
203
 
 
204
        orbindex(1) = (1 - orbspin(1)+
 
205
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p5b-1)+i-1))/2
 
206
        orbindex(2) = (1 - orbspin(2)+
 
207
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h6b-1)+j-1))/2
 
208
 
 
209
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
210
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
211
 
 
212
cjb ===========================================================================
 
213
 
 
214
        if(isactive(orbindex(1),orbspin(1)+1).and.
 
215
     1 isactive(orbindex(2),orbspin(2)+1).or.(.not.limprovet)) then
 
216
 
 
217
ccc =====
 
218
c       if(nodezero)write(6,"('ACTIVITY: ',2L2)")
 
219
c     1 isactive(orbindex(1),orbspin(1)+1),
 
220
c     1 isactive(orbindex(2),orbspin(2)+1)
 
221
c       if(nodezero)write(6,"('DEBUG: ',5I4)")
 
222
c     1 orbindex(1),orbspin(1),
 
223
c     1 orbindex(2),orbspin(2),iref
 
224
ccc =====
 
225
 
 
226
       do iu=1,2
 
227
        do n=1,nmo(iu)
 
228
          ioccnew(n,iu) = iocc(n,iref,iu)
 
229
        enddo
 
230
       enddo
 
231
 
 
232
      if(iocc(orbindex(1),iref,orbspin(1)+1).eq.
 
233
     1 iocc(orbindex(2),iref,orbspin(2)+1)) then
 
234
          goto 200
 
235
       endif
 
236
 
 
237
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(2),iref,
 
238
     1 orbspin(2)+1)
 
239
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(1),iref,
 
240
     1 orbspin(1)+1)
 
241
 
 
242
       do n=1,nref
 
243
       isfound = .true.
 
244
        do iu=1,2
 
245
         do o=1,nmo(iu)
 
246
          isfound = isfound.and.(iocc(o,n,iu).eq.ioccnew(o,iu))
 
247
         enddo
 
248
        enddo
 
249
       if(isfound) then
 
250
c          write(LuOut,"('Internal amplitude',I4,'->',I4)")iref,n
 
251
c        write(LuOut,"('1Internal amplitude',I4,'->',I4,2F16.12)")iref,n,
 
252
c     1 dbl_mb(k_r1+m-1)
 
253
            if(iref.ne.n) then
 
254
 
 
255
             dist = 0
 
256
       do iu=1,2
 
257
        do i1=1,nmo(iu)
 
258
          ioccnew(i1,iu) = iocc(i1,iref,iu)
 
259
        enddo
 
260
       enddo
 
261
 
 
262
             i2 = 0
 
263
             do i1=min(orbindex(1),orbindex(2)),
 
264
     1 max(orbindex(1),orbindex(2))
 
265
               i2 = i2 + 1
 
266
               if(i2.lt.abs(orbindex(1)-orbindex(2))) then
 
267
                  if(iocc(i1+1,iref,orbspin(1)+1).eq.1) then
 
268
                     dist=dist+1
 
269
                  endif
 
270
               endif
 
271
             enddo
 
272
 
 
273
         dsmult = 1.0d0
 
274
 
 
275
         if(mod(dist,2).eq.1) dsmult = -dsmult
 
276
 
 
277
         dbl_mb(k_heff+n-1+(iref-1)*nref) = dbl_mb(k_r1+k-1)*dsmult
 
278
c     1 dbl_mb(k_heff+n-1+(iref-1)*nref)
 
279
 
 
280
           if(lusesub) then
 
281
         call ga_put(g_heff,nref*(iref-1)+n,nref*(iref-1)+n,1,1,
 
282
     1   dbl_mb(k_r1+k-1)*dsmult,1)
 
283
           endif
 
284
 
 
285
            endif
 
286
          goto 200
 
287
       endif
 
288
       enddo
 
289
 
 
290
 200   continue
 
291
 
 
292
        if((.not.isfound).and.(abs(dbl_mb(k_r1+k-1)).gt.1.0d-5)) then
 
293
         if(nodezero) then
 
294
           write(LuOut,"('DEBUG: ',4F16.12)")dbl_mb(k_r1+k-1)
 
295
           write(LuOut,"('YOU ARE USING INCOMPLETE MODEL SPACE!')")
 
296
         endif
 
297
c         call errquit('YOU ARE USING INCOMPLETE MODEL SPACE!',1,MA_ERR)
 
298
        endif
 
299
 
 
300
        endif ! active
 
301
 
 
302
        enddo
 
303
        enddo
 
304
 
 
305
        if (.not.ma_pop_stack(l_r1))
 
306
     1   call errquit('tce_mrcc_iface_r1: MA problem',1,MA_ERR)
 
307
 
 
308
      endif
 
309
      endif
 
310
      endif
 
311
      endif
 
312
      enddo
 
313
      enddo
 
314
c
 
315
c ---------------
 
316
c    R2 active
 
317
c ---------------
 
318
c
 
319
      nxt = 0
 
320
      count = 0
 
321
      if(limprovet) then
 
322
       if(lusesub) then
 
323
         nprocs=GA_pgroup_NNODES(mypgid)
 
324
         nxt=NXTASKsub(nprocs,1,mypgid)
 
325
       else
 
326
         nprocs=GA_NNODES()
 
327
         nxt=NXTASK(nprocs,1)
 
328
       endif
 
329
      count = 0
 
330
      endif
 
331
 
 
332
      DO p1b = noab+1,noab+nvab
 
333
      DO p2b = p1b,noab+nvab
 
334
      DO h3b = 1,noab
 
335
      DO h4b = h3b,noab
 
336
 
 
337
      IF ((nxt.eq.count).or.(.not.limprovet)) THEN
 
338
 
 
339
      IF (int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1) .eq. int_mb(k_spin+h
 
340
     &3b-1)+int_mb(k_spin+h4b-1)) THEN
 
341
 
 
342
      IF (ieor(int_mb(k_sym+p1b-1),ieor(int_mb(k_sym+p2b-1),ieor(int_mb(
 
343
     &k_sym+h3b-1),int_mb(k_sym+h4b-1)))) .eq. irrep_t) THEN
 
344
 
 
345
      IF ((.not.restricted).or.(int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1
 
346
     &)+int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1).ne.8)) THEN
 
347
 
 
348
      if(log_mb(k_isactive(iref)+p1b-1).and.
 
349
     1 log_mb(k_isactive(iref)+p2b-1).and.
 
350
     2 log_mb(k_isactive(iref)+h3b-1).and.
 
351
     3 log_mb(k_isactive(iref)+h4b-1)) then
 
352
 
 
353
      size = int_mb(k_range+p1b-1) * int_mb(k_range+p2b-1) * int_
 
354
     &mb(k_range+h3b-1) * int_mb(k_range+h4b-1)
 
355
 
 
356
        if (.not.ma_push_get(mt_dbl,size,'r2mi',l_r2,k_r2))
 
357
     1   call errquit('tce_mrcc_iface_r2: MA problem',0,MA_ERR)
 
358
 
 
359
        call get_hash_block(d_r2m(iref),dbl_mb(k_r2),size,
 
360
     1   int_mb(k_r2_offsetm(iref)),h4b-1+noab*(h3b-1+noab*(p2b-
 
361
     &noab-1+nvab*(p1b - noab - 1))))
 
362
 
 
363
c         write(LuOut,"(I4,L3,L3,L3,L3)")
 
364
c     1 iref,log_mb(k_isactive(iref)+p1b-1),
 
365
c     1 log_mb(k_isactive(iref)+p2b-1),log_mb(k_isactive(iref)+h3b-1),
 
366
c     1 log_mb(k_isactive(iref)+h4b-1)
 
367
       m = 0
 
368
 
 
369
        do i=1,int_mb(k_range+p1b-1)
 
370
        do j=1,int_mb(k_range+p2b-1)
 
371
        do k=1,int_mb(k_range+h3b-1)
 
372
        do l=1,int_mb(k_range+h4b-1)
 
373
         m = m + 1
 
374
c         write(LuOut,"(I4,'(',I4,I4,I4,I4,'):',2F16.12)")
 
375
c     1 iref,i,j,k,l,dbl_mb(k_r2+m-1)
 
376
c         write(LuOut,*)int_mb(k_spin+p1b-1)
 
377
 
 
378
        orbspin(1) = int_mb(k_spin+p1b-1)-1
 
379
        orbspin(2) = int_mb(k_spin+p2b-1)-1
 
380
        orbspin(3) = int_mb(k_spin+h3b-1)-1
 
381
        orbspin(4) = int_mb(k_spin+h4b-1)-1
 
382
 
 
383
        orbindex(1) = (1 - orbspin(1)+ 
 
384
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p1b-1)+i-1))/2
 
385
        orbindex(2) = (1 - orbspin(2)+
 
386
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p2b-1)+j-1))/2
 
387
        orbindex(3) = (1 - orbspin(3)+
 
388
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h3b-1)+k-1))/2
 
389
        orbindex(4) = (1 - orbspin(4)+
 
390
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h4b-1)+l-1))/2
 
391
 
 
392
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
393
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
394
        orbindex(3) = moindexes(orbindex(3),orbspin(3)+1,iref)
 
395
        orbindex(4) = moindexes(orbindex(4),orbspin(4)+1,iref)
 
396
 
 
397
cjb ===========================================================================
 
398
 
 
399
        if(isactive(orbindex(1),orbspin(1)+1).and.
 
400
     1 isactive(orbindex(2),orbspin(2)+1).and.
 
401
     2 isactive(orbindex(3),orbspin(3)+1).and.
 
402
     3 isactive(orbindex(4),orbspin(4)+1).or.(.not.limprovet)) then
 
403
 
 
404
c        write(LuOut,"('Real indexes: [',I4,I4,I4,I4,']',
 
405
c     1 '[',I2,I2,I2,I2,']')")
 
406
c     1 orbindex(1),orbindex(2),orbindex(3),orbindex(4),
 
407
c     1 orbspin(1),orbspin(2),orbspin(3),orbspin(4)
 
408
 
 
409
       do iu=1,2
 
410
        do n=1,nmo(iu)
 
411
          ioccnew(n,iu) = iocc(n,iref,iu) 
 
412
        enddo
 
413
       enddo 
 
414
        
 
415
       if(((orbindex(1).eq.orbindex(2)).and.(orbspin(1).eq.orbspin(2)))
 
416
     1 .or.((orbindex(3).eq.orbindex(4)).and.(orbspin(3).eq.
 
417
     2 orbspin(4)))) then
 
418
         goto 100
 
419
       endif
 
420
 
 
421
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(3),iref,
 
422
     1 orbspin(3)+1)
 
423
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(4),iref,
 
424
     1 orbspin(4)+1)
 
425
         ioccnew(orbindex(3),orbspin(3)+1) = iocc(orbindex(1),iref,
 
426
     1 orbspin(1)+1)
 
427
         ioccnew(orbindex(4),orbspin(4)+1) = iocc(orbindex(2),iref,
 
428
     1 orbspin(2)+1)
 
429
 
 
430
       do n=1,nref
 
431
       isfound = .true.
 
432
        do iu=1,2
 
433
         do o=1,nmo(iu)
 
434
          isfound = isfound.and.(iocc(o,n,iu).eq.ioccnew(o,iu))
 
435
         enddo
 
436
        enddo
 
437
       if(isfound) then
 
438
c        write(LuOut,"('2Internal amplitude',I4,'->',I4,2F16.12)")iref,n,
 
439
c     1 dbl_mb(k_r2+m-1)
 
440
            if(iref.ne.n) then
 
441
 
 
442
             dist = 0
 
443
       do iu=1,2
 
444
        do i1=1,nmo(iu)
 
445
          ioccnew(i1,iu) = iocc(i1,iref,iu)
 
446
        enddo
 
447
       enddo
 
448
 
 
449
             i2 = 0
 
450
             do i1=min(orbindex(1),orbindex(3)),
 
451
     1 max(orbindex(1),orbindex(3))
 
452
               i2 = i2 + 1
 
453
               if(i2.lt.abs(orbindex(1)-orbindex(3))) then
 
454
                  if(iocc(i1+1,iref,orbspin(1)+1).eq.1) then
 
455
                     dist=dist+1
 
456
                  endif
 
457
               endif
 
458
             enddo
 
459
 
 
460
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(3),iref,
 
461
     1 orbspin(3)+1)
 
462
         ioccnew(orbindex(3),orbspin(3)+1) = iocc(orbindex(1),iref,
 
463
     1 orbspin(1)+1)
 
464
 
 
465
             i2 = 0
 
466
             do i1=min(orbindex(2),orbindex(4)),
 
467
     1 max(orbindex(2),orbindex(4))
 
468
               i2 = i2 + 1
 
469
               if(i2.lt.abs(orbindex(2)-orbindex(4))) then
 
470
                  if(ioccnew(i1+1,orbspin(2)+1).eq.1) then
 
471
                     dist=dist+1
 
472
                  endif
 
473
               endif
 
474
             enddo
 
475
 
 
476
            dsmult = 1.0d0
 
477
         if(mod(dist,2).eq.1) dsmult = -dsmult
 
478
 
 
479
 
 
480
         dbl_mb(k_heff+n-1+(iref-1)*nref) = dbl_mb(k_r2+m-1)*dsmult
 
481
c     1 dbl_mb(k_heff+n-1+(iref-1)*nref)
 
482
 
 
483
           if(lusesub) then
 
484
         call ga_put(g_heff,nref*(iref-1)+n,nref*(iref-1)+n,1,1,
 
485
     1   dbl_mb(k_r2+m-1)*dsmult,1)
 
486
           endif
 
487
 
 
488
 
 
489
            endif
 
490
          goto 100
 
491
       endif
 
492
       enddo
 
493
 
 
494
 100    continue
 
495
 
 
496
        if((.not.isfound).and.(abs(dbl_mb(k_r2+m-1)).gt.1.0d-5)) then
 
497
         if(nodezero) then
 
498
           write(LuOut,"('DEBUG: ',4F16.12)")dbl_mb(k_r2+m-1)
 
499
           write(LuOut,"('YOU ARE USING INCOMPLETE MODEL SPACE!')")
 
500
         endif
 
501
c         call errquit('YOU ARE USING INCOMPLETE MODEL SPACE!',2,MA_ERR)
 
502
        endif
 
503
 
 
504
        endif ! act
 
505
 
 
506
        enddo
 
507
        enddo
 
508
        enddo
 
509
        enddo
 
510
 
 
511
        if (.not.ma_pop_stack(l_r2))
 
512
     1   call errquit('tce_mrcc_iface_r2: MA problem',1,MA_ERR)
 
513
      END IF
 
514
      END IF
 
515
      END IF
 
516
      endif
 
517
      if(limprovet) then
 
518
       if(lusesub) then
 
519
        nxt=NXTASKsub(nprocs,1,mypgid)
 
520
       else
 
521
        nxt=NXTASK(nprocs,1)
 
522
       endif
 
523
      endif
 
524
      endif
 
525
       if(limprovet)count = count + 1
 
526
      END DO
 
527
      END DO
 
528
      END DO
 
529
      END DO
 
530
 
 
531
      if(limprovet) then
 
532
      if(lusesub) then
 
533
      nxt=NXTASKsub(-nprocs,1,mypgid)
 
534
      call GA_Pgroup_SYNC(mypgid)
 
535
      else
 
536
      nxt=NXTASKsub(-nprocs,1)
 
537
      call GA_SYNC()
 
538
      endif
 
539
      endif
 
540
c
 
541
c ---------------
 
542
c    R3 active
 
543
c ---------------
 
544
c
 
545
      if(needr3act) then
 
546
      DO p4b = noab+1,noab+nvab
 
547
      DO p5b = p4b,noab+nvab
 
548
      DO p6b = p5b,noab+nvab
 
549
      DO h1b = 1,noab
 
550
      DO h2b = h1b,noab
 
551
      DO h3b = h2b,noab
 
552
      IF ((.not.restricted).or.(int_mb(k_spin+p4b-1)+int_mb(k_spin+p5b-1
 
553
     &)+int_mb(k_spin+p6b-1)+int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-1)+i
 
554
     &nt_mb(k_spin+h3b-1).ne.12)) THEN
 
555
      IF (int_mb(k_spin+p4b-1)+int_mb(k_spin+p5b-1)+int_mb(k_spin+p6b-1)
 
556
     & .eq. int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-1)+int_mb(k_spin+h3b-
 
557
     &1)) THEN
 
558
      IF (ieor(int_mb(k_sym+p4b-1),ieor(int_mb(k_sym+p5b-1),ieor(int_mb(
 
559
     &k_sym+p6b-1),ieor(int_mb(k_sym+h1b-1),ieor(int_mb(k_sym+h2b-1),int
 
560
     &_mb(k_sym+h3b-1)))))) .eq. ieor(irrep_v,irrep_t)) THEN
 
561
      IF ((log_mb(k_active+p4b-1).eqv..true.).and.(log_mb(k_active+p5b-1
 
562
     &).eqv..true.).and.(log_mb(k_active+p6b-1).eqv..true.).and.(log_mb(
 
563
     &k_active+h1b-1).eqv..true.).and.(log_mb(k_active+h2b-1).eqv..true.
 
564
     &).and.(log_mb(k_active+h3b-1).eqv..true.)) THEN
 
565
 
 
566
      size = int_mb(k_range+p4b-1) * int_mb(k_range+p5b-1) * int_
 
567
     &mb(k_range+p6b-1) * int_mb(k_range+h1b-1) * int_mb(k_range+h2b-1)
 
568
     &* int_mb(k_range+h3b-1)
 
569
 
 
570
      if (.not.ma_push_get(mt_dbl,size,'r3mi',l_r3,k_r3))
 
571
     1   call errquit('tce_mrcc_iface_r3: MA problem',0,MA_ERR)
 
572
 
 
573
      call get_hash_block(d_r3m(iref),dbl_mb(k_r3),size,
 
574
     1   int_mb(k_r3_offsetm(iref)),h3b - 1 + noab * (h2b - 1 + noab *
 
575
     1 (h1b - 1 + noab * (p6b - noab - 1 + nvab * (p5b - noab - 1 + 
 
576
     1 nvab * (p4b - noab - 1))))))
 
577
 
 
578
       m = 0
 
579
 
 
580
        orbspin(1) = int_mb(k_spin+p4b-1)-1
 
581
        orbspin(2) = int_mb(k_spin+p5b-1)-1
 
582
        orbspin(3) = int_mb(k_spin+p6b-1)-1
 
583
        orbspin(4) = int_mb(k_spin+h1b-1)-1
 
584
        orbspin(5) = int_mb(k_spin+h2b-1)-1
 
585
        orbspin(6) = int_mb(k_spin+h3b-1)-1
 
586
 
 
587
        do p4=1,int_mb(k_range+p4b-1)
 
588
        do p5=1,int_mb(k_range+p5b-1)
 
589
        do p6=1,int_mb(k_range+p6b-1)
 
590
        do h1=1,int_mb(k_range+h1b-1)
 
591
        do h2=1,int_mb(k_range+h2b-1)
 
592
        do h3=1,int_mb(k_range+h3b-1)
 
593
 
 
594
        m = m + 1
 
595
 
 
596
        orbindex(1) = (1 - orbspin(1)+
 
597
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p4b-1)+p4-1))/2
 
598
        orbindex(2) = (1 - orbspin(2)+
 
599
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p5b-1)+p5-1))/2
 
600
        orbindex(3) = (1 - orbspin(3)+
 
601
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p6b-1)+p6-1))/2
 
602
        orbindex(4) = (1 - orbspin(4)+
 
603
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h1b-1)+h1-1))/2
 
604
        orbindex(5) = (1 - orbspin(5)+
 
605
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h2b-1)+h2-1))/2
 
606
        orbindex(6) = (1 - orbspin(6)+
 
607
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h3b-1)+h3-1))/2
 
608
 
 
609
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
610
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
611
        orbindex(3) = moindexes(orbindex(3),orbspin(3)+1,iref)
 
612
        orbindex(4) = moindexes(orbindex(4),orbspin(4)+1,iref)
 
613
        orbindex(5) = moindexes(orbindex(5),orbspin(5)+1,iref)
 
614
        orbindex(6) = moindexes(orbindex(6),orbspin(6)+1,iref)
 
615
 
 
616
c        write(LuOut,"('Real indexes: [',I4,I4,I4,I4,I4,I4,']')")
 
617
c     1 orbindex(1),orbindex(2),orbindex(3),orbindex(4),orbindex(5),
 
618
c     1 orbindex(6)
 
619
c        write(LuOut,"('Spin indexes : [',I4,I4,I4,I4,I4,I4,']')")
 
620
c     1 orbspin(1),orbspin(2),orbspin(3),orbspin(4),orbspin(5),orbspin(6)
 
621
 
 
622
       do iu=1,2
 
623
        do n=1,nbf
 
624
          ioccnew(n,iu) = iocc(n,iref,iu)
 
625
        enddo
 
626
       enddo
 
627
 
 
628
      if((iocc(orbindex(1),iref,orbspin(1)+1).eq.
 
629
     1 iocc(orbindex(4),iref,orbspin(4)+1)).or.
 
630
     2 (iocc(orbindex(2),iref,orbspin(2)+1).eq.
 
631
     3 iocc(orbindex(5),iref,orbspin(5)+1)).or.
 
632
     4 (iocc(orbindex(3),iref,orbspin(3)+1).eq.
 
633
     1 iocc(orbindex(6),iref,orbspin(6)+1))) then
 
634
          goto 300
 
635
       endif
 
636
 
 
637
      if((orbspin(1).ne.orbspin(4)).or.
 
638
     1   (orbspin(2).ne.orbspin(5)).or.
 
639
     2   (orbspin(3).ne.orbspin(6))) then
 
640
          goto 300
 
641
       endif
 
642
 
 
643
       if(
 
644
     1 ((orbindex(1).eq.orbindex(2)).and.(orbspin(1).eq.orbspin(2))).or.
 
645
     1 ((orbindex(1).eq.orbindex(3)).and.(orbspin(1).eq.orbspin(3))).or.
 
646
     1 ((orbindex(2).eq.orbindex(3)).and.(orbspin(2).eq.orbspin(3))).or.
 
647
     1 ((orbindex(4).eq.orbindex(5)).and.(orbspin(4).eq.orbspin(5))).or.
 
648
     1 ((orbindex(4).eq.orbindex(6)).and.(orbspin(4).eq.orbspin(6))).or.
 
649
     1 ((orbindex(5).eq.orbindex(6)).and.(orbspin(5).eq.orbspin(6)))
 
650
     1 ) then
 
651
          goto 300
 
652
       endif
 
653
 
 
654
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(4),iref,
 
655
     1 orbspin(4)+1)
 
656
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(5),iref,
 
657
     1 orbspin(5)+1)
 
658
         ioccnew(orbindex(3),orbspin(3)+1) = iocc(orbindex(6),iref,
 
659
     1 orbspin(6)+1)
 
660
         ioccnew(orbindex(4),orbspin(4)+1) = iocc(orbindex(1),iref,
 
661
     1 orbspin(1)+1)
 
662
         ioccnew(orbindex(5),orbspin(5)+1) = iocc(orbindex(2),iref,
 
663
     1 orbspin(2)+1)
 
664
         ioccnew(orbindex(6),orbspin(6)+1) = iocc(orbindex(3),iref,
 
665
     1 orbspin(3)+1)
 
666
 
 
667
 
 
668
       do n=1,nref
 
669
       isfound = .true.
 
670
        do iu=1,2
 
671
         do o=1,nbf
 
672
          isfound = isfound.and.(iocc(o,n,iu).eq.ioccnew(o,iu))
 
673
         enddo
 
674
        enddo
 
675
       if(isfound) then
 
676
c         write(LuOut,"('Internal amplitude',I4,'->',I4,2F16.12)")iref,n,
 
677
c     1 dbl_mb(k_r3+m-1)
 
678
            if(iref.ne.n) then
 
679
cJB START
 
680
             dist = 0
 
681
       do iu=1,2
 
682
        do i1=1,nmo(iu)
 
683
          ioccnew(i1,iu) = iocc(i1,iref,iu)
 
684
        enddo
 
685
       enddo
 
686
 
 
687
             i2 = 0
 
688
             do i1=min(orbindex(1),orbindex(4)),
 
689
     1 max(orbindex(1),orbindex(4))
 
690
               i2 = i2 + 1
 
691
               if(i2.lt.abs(orbindex(1)-orbindex(4))) then
 
692
                  if(iocc(i1+1,iref,orbspin(1)+1).eq.1) then
 
693
                     dist=dist+1
 
694
                  endif
 
695
               endif
 
696
             enddo
 
697
 
 
698
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(4),iref,
 
699
     1 orbspin(4)+1)
 
700
         ioccnew(orbindex(4),orbspin(4)+1) = iocc(orbindex(1),iref,
 
701
     1 orbspin(1)+1)
 
702
 
 
703
             i2 = 0
 
704
             do i1=min(orbindex(2),orbindex(5)),
 
705
     1 max(orbindex(2),orbindex(5))
 
706
               i2 = i2 + 1
 
707
               if(i2.lt.abs(orbindex(2)-orbindex(5))) then
 
708
                  if(ioccnew(i1+1,orbspin(5)+1).eq.1) then
 
709
                     dist=dist+1
 
710
                  endif
 
711
               endif
 
712
             enddo
 
713
 
 
714
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(5),iref,
 
715
     1 orbspin(5)+1)
 
716
         ioccnew(orbindex(5),orbspin(5)+1) = iocc(orbindex(2),iref,
 
717
     1 orbspin(2)+1)
 
718
 
 
719
             i2 = 0
 
720
             do i1=min(orbindex(3),orbindex(6)),
 
721
     1 max(orbindex(3),orbindex(6))
 
722
               i2 = i2 + 1
 
723
               if(i2.lt.abs(orbindex(3)-orbindex(6))) then
 
724
                  if(ioccnew(i1+1,orbspin(6)+1).eq.1) then
 
725
                     dist=dist+1
 
726
                  endif
 
727
               endif
 
728
             enddo
 
729
 
 
730
            dsmult = 1.0d0
 
731
         if(mod(dist,2).eq.1) dsmult = -dsmult
 
732
 
 
733
cJB END
 
734
              dbl_mb(k_heff+n-1+(iref-1)*nref) = 
 
735
     1 dbl_mb(k_r3+m-1)*dsmult
 
736
 
 
737
            endif
 
738
          goto 300
 
739
       endif
 
740
       enddo
 
741
 
 
742
 300   continue
 
743
 
 
744
        enddo
 
745
        enddo
 
746
        enddo
 
747
        enddo
 
748
        enddo
 
749
        enddo
 
750
 
 
751
        if (.not.ma_pop_stack(l_r3))
 
752
     1   call errquit('tce_mrcc_iface_r3: MA problem',1,MA_ERR)
 
753
 
 
754
      END IF
 
755
      END IF
 
756
      END IF
 
757
      END IF
 
758
      END DO
 
759
      END DO
 
760
      END DO
 
761
      END DO
 
762
      END DO
 
763
      END DO
 
764
      endif
 
765
c
 
766
c ---------------
 
767
c    R4 active
 
768
c ---------------
 
769
c
 
770
      if(needr4act) then
 
771
c      DO p5b = noab+1,noab+nvab
 
772
c      DO p6b = p5b,noab+nvab
 
773
c      DO p7b = noab+1,noab+nvab
 
774
c      DO p8b = p7b,noab+nvab
 
775
c      DO h1b = 1,noab
 
776
c      DO h2b = 1,noab
 
777
c      DO h3b = h2b,noab
 
778
c      DO h4b = h3b,noab
 
779
 
 
780
      DO p5b = noab+1,noab+nvab
 
781
      DO p6b = p5b,noab+nvab
 
782
      DO p7b = p6b,noab+nvab
 
783
      DO p8b = p7b,noab+nvab
 
784
      DO h1b = 1,noab
 
785
      DO h2b = h1b,noab
 
786
      DO h3b = h2b,noab
 
787
      DO h4b = h3b,noab
 
788
 
 
789
      IF ((.not.restricted).or.(int_mb(k_spin+p5b-1)+int_mb(k_spin+p6b-1
 
790
     &)+int_mb(k_spin+p7b-1)+int_mb(k_spin+p8b-1)+int_mb(k_spin+h1b-1)+i
 
791
     &nt_mb(k_spin+h2b-1)+int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1).ne.1
 
792
     &6)) THEN
 
793
      IF (int_mb(k_spin+p5b-1)+int_mb(k_spin+p6b-1)+int_mb(k_spin+p7b-1)
 
794
     &+int_mb(k_spin+p8b-1) .eq. int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-
 
795
     &1)+int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1)) THEN
 
796
      IF (ieor(int_mb(k_sym+p5b-1),ieor(int_mb(k_sym+p6b-1),ieor(int_mb(
 
797
     &k_sym+p7b-1),ieor(int_mb(k_sym+p8b-1),ieor(int_mb(k_sym+h1b-1),ieo
 
798
     &r(int_mb(k_sym+h2b-1),ieor(int_mb(k_sym+h3b-1),int_mb(k_sym+h4b-1)
 
799
     &))))))) .eq. ieor(irrep_v,ieor(irrep_t,irrep_t))) THEN
 
800
      IF ((log_mb(k_active+p5b-1).eqv..true.).and.(log_mb(k_active+p6b-1
 
801
     &).eqv..true.).and.(log_mb(k_active+p7b-1).eqv..true.).and.(log_mb(
 
802
     &k_active+p8b-1).eqv..true.).and.(log_mb(k_active+h1b-1).eqv..true.
 
803
     &).and.(log_mb(k_active+h2b-1).eqv..true.).and.(log_mb(k_active+h3b
 
804
     &-1).eqv..true.).and.(log_mb(k_active+h4b-1).eqv..true.)) THEN
 
805
 
 
806
      size = int_mb(k_range+p5b-1) * int_mb(k_range+p6b-1) * int_
 
807
     &mb(k_range+p7b-1) * int_mb(k_range+p8b-1) * int_mb(k_range+h1b-1)
 
808
     &* int_mb(k_range+h2b-1) * int_mb(k_range+h3b-1) * int_mb(k_range+h
 
809
     &4b-1)
 
810
 
 
811
      if (.not.ma_push_get(mt_dbl,size,'r4mi',l_r4,k_r4))
 
812
     1   call errquit('tce_mrcc_iface_r4: MA problem',0,MA_ERR)
 
813
 
 
814
      call get_hash_block(d_r4m(iref),dbl_mb(k_r4),size,
 
815
     1   int_mb(k_r4_offsetm(iref)),(h4b - 1 + noab * (h3b - 1 + noab *
 
816
     1(h2b - 1 + noab * (h1b - 1 + noab * (p8b - noab - 1 + nvab * (p7b
 
817
     1 - noab - 1 + nvab * (p6b - noab - 1 + nvab * (p5b - noab - 1))))
 
818
     1)))))
 
819
 
 
820
       m = 0
 
821
 
 
822
        orbspin(1) = int_mb(k_spin+p5b-1)-1
 
823
        orbspin(2) = int_mb(k_spin+p6b-1)-1
 
824
        orbspin(3) = int_mb(k_spin+p7b-1)-1
 
825
        orbspin(4) = int_mb(k_spin+p8b-1)-1
 
826
        orbspin(5) = int_mb(k_spin+h1b-1)-1
 
827
        orbspin(6) = int_mb(k_spin+h2b-1)-1
 
828
        orbspin(7) = int_mb(k_spin+h3b-1)-1
 
829
        orbspin(8) = int_mb(k_spin+h4b-1)-1
 
830
 
 
831
        do p5=1,int_mb(k_range+p5b-1)
 
832
        do p6=1,int_mb(k_range+p6b-1)
 
833
        do p7=1,int_mb(k_range+p7b-1)
 
834
        do p8=1,int_mb(k_range+p8b-1)
 
835
        do h1=1,int_mb(k_range+h1b-1)
 
836
        do h2=1,int_mb(k_range+h2b-1)
 
837
        do h3=1,int_mb(k_range+h3b-1)
 
838
        do h4=1,int_mb(k_range+h4b-1)
 
839
 
 
840
        m = m + 1
 
841
 
 
842
        orbindex(1) = (1 - orbspin(1)+
 
843
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p5b-1)+p5-1))/2
 
844
        orbindex(2) = (1 - orbspin(2)+
 
845
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p6b-1)+p6-1))/2
 
846
        orbindex(3) = (1 - orbspin(3)+
 
847
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p7b-1)+p7-1))/2
 
848
        orbindex(4) = (1 - orbspin(4)+
 
849
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p8b-1)+p8-1))/2
 
850
        orbindex(5) = (1 - orbspin(5)+
 
851
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h1b-1)+h1-1))/2
 
852
        orbindex(6) = (1 - orbspin(6)+
 
853
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h2b-1)+h2-1))/2
 
854
        orbindex(7) = (1 - orbspin(7)+
 
855
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h3b-1)+h3-1))/2
 
856
        orbindex(8) = (1 - orbspin(8)+
 
857
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h4b-1)+h4-1))/2
 
858
 
 
859
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
860
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
861
        orbindex(3) = moindexes(orbindex(3),orbspin(3)+1,iref)
 
862
        orbindex(4) = moindexes(orbindex(4),orbspin(4)+1,iref)
 
863
        orbindex(5) = moindexes(orbindex(5),orbspin(5)+1,iref)
 
864
        orbindex(6) = moindexes(orbindex(6),orbspin(6)+1,iref)
 
865
        orbindex(7) = moindexes(orbindex(7),orbspin(7)+1,iref)
 
866
        orbindex(8) = moindexes(orbindex(8),orbspin(8)+1,iref)
 
867
 
 
868
       do iu=1,2
 
869
        do n=1,nbf
 
870
          ioccnew(n,iu) = iocc(n,iref,iu)
 
871
        enddo
 
872
       enddo
 
873
 
 
874
      if((iocc(orbindex(1),iref,orbspin(1)+1).eq.
 
875
     1 iocc(orbindex(5),iref,orbspin(5)+1)).or.
 
876
     2 (iocc(orbindex(2),iref,orbspin(2)+1).eq.
 
877
     3 iocc(orbindex(6),iref,orbspin(6)+1)).or.
 
878
     4 (iocc(orbindex(3),iref,orbspin(3)+1).eq.
 
879
     1 iocc(orbindex(7),iref,orbspin(7)+1)).or.
 
880
     2 (iocc(orbindex(4),iref,orbspin(4)+1).eq.
 
881
     3 iocc(orbindex(8),iref,orbspin(8)+1))) then
 
882
          goto 400
 
883
       endif
 
884
 
 
885
      if((orbspin(1).ne.orbspin(5)).or.
 
886
     1   (orbspin(2).ne.orbspin(6)).or.
 
887
     2   (orbspin(3).ne.orbspin(7)).or.
 
888
     3   (orbspin(4).ne.orbspin(8))) then
 
889
          goto 400
 
890
       endif
 
891
 
 
892
      if(
 
893
     1 ((orbindex(1).eq.orbindex(2)).and.(orbspin(1).eq.orbspin(2))).or.
 
894
     1 ((orbindex(1).eq.orbindex(3)).and.(orbspin(1).eq.orbspin(3))).or.
 
895
     1 ((orbindex(1).eq.orbindex(4)).and.(orbspin(1).eq.orbspin(4))).or.
 
896
     1 ((orbindex(2).eq.orbindex(3)).and.(orbspin(2).eq.orbspin(3))).or.
 
897
     1 ((orbindex(2).eq.orbindex(4)).and.(orbspin(2).eq.orbspin(4))).or.
 
898
     1 ((orbindex(3).eq.orbindex(4)).and.(orbspin(3).eq.orbspin(4))).or.
 
899
     1 ((orbindex(5).eq.orbindex(6)).and.(orbspin(5).eq.orbspin(6))).or.
 
900
     1 ((orbindex(5).eq.orbindex(7)).and.(orbspin(5).eq.orbspin(7))).or.
 
901
     1 ((orbindex(5).eq.orbindex(8)).and.(orbspin(5).eq.orbspin(8))).or.
 
902
     1 ((orbindex(6).eq.orbindex(7)).and.(orbspin(6).eq.orbspin(7))).or.
 
903
     1 ((orbindex(6).eq.orbindex(8)).and.(orbspin(6).eq.orbspin(8))).or.
 
904
     1 ((orbindex(7).eq.orbindex(8)).and.(orbspin(7).eq.orbspin(8)))
 
905
     1 ) then
 
906
          goto 400
 
907
       endif
 
908
 
 
909
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(5),iref,
 
910
     1 orbspin(5)+1)
 
911
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(6),iref,
 
912
     1 orbspin(6)+1)
 
913
         ioccnew(orbindex(3),orbspin(3)+1) = iocc(orbindex(7),iref,
 
914
     1 orbspin(7)+1)
 
915
         ioccnew(orbindex(4),orbspin(4)+1) = iocc(orbindex(8),iref,
 
916
     1 orbspin(8)+1)
 
917
         ioccnew(orbindex(5),orbspin(5)+1) = iocc(orbindex(1),iref,
 
918
     1 orbspin(1)+1)
 
919
         ioccnew(orbindex(6),orbspin(6)+1) = iocc(orbindex(2),iref,
 
920
     1 orbspin(2)+1)
 
921
         ioccnew(orbindex(7),orbspin(7)+1) = iocc(orbindex(3),iref,
 
922
     1 orbspin(3)+1)
 
923
         ioccnew(orbindex(8),orbspin(8)+1) = iocc(orbindex(4),iref,
 
924
     1 orbspin(4)+1)
 
925
 
 
926
       do n=1,nref
 
927
       isfound = .true.
 
928
        do iu=1,2
 
929
         do o=1,nbf
 
930
          isfound = isfound.and.(iocc(o,n,iu).eq.ioccnew(o,iu))
 
931
         enddo
 
932
        enddo
 
933
       if(isfound) then
 
934
c         write(LuOut,"('Internal amplitude',I4,'->',I4,2F16.12)")iref,n,
 
935
c     1 dbl_mb(k_r4+m-1)
 
936
            if(iref.ne.n) then
 
937
c ===============
 
938
             dist = 0
 
939
       do iu=1,2
 
940
        do i1=1,nmo(iu)
 
941
          ioccnew(i1,iu) = iocc(i1,iref,iu)
 
942
        enddo
 
943
       enddo
 
944
 
 
945
             i2 = 0
 
946
             do i1=min(orbindex(1),orbindex(5)),
 
947
     1 max(orbindex(1),orbindex(5))
 
948
               i2 = i2 + 1
 
949
               if(i2.lt.abs(orbindex(1)-orbindex(5))) then
 
950
                  if(iocc(i1+1,iref,orbspin(1)+1).eq.1) then
 
951
                     dist=dist+1
 
952
                  endif
 
953
               endif
 
954
             enddo
 
955
 
 
956
         ioccnew(orbindex(1),orbspin(1)+1) = iocc(orbindex(5),iref,
 
957
     1 orbspin(5)+1)
 
958
         ioccnew(orbindex(5),orbspin(5)+1) = iocc(orbindex(1),iref,
 
959
     1 orbspin(1)+1)
 
960
 
 
961
             i2 = 0
 
962
             do i1=min(orbindex(2),orbindex(6)),
 
963
     1 max(orbindex(2),orbindex(6))
 
964
               i2 = i2 + 1
 
965
               if(i2.lt.abs(orbindex(2)-orbindex(6))) then
 
966
                  if(iocc(i1+1,iref,orbspin(2)+1).eq.1) then
 
967
                     dist=dist+1
 
968
                  endif
 
969
               endif
 
970
             enddo
 
971
 
 
972
         ioccnew(orbindex(2),orbspin(2)+1) = iocc(orbindex(6),iref,
 
973
     1 orbspin(6)+1)
 
974
         ioccnew(orbindex(6),orbspin(6)+1) = iocc(orbindex(2),iref,
 
975
     1 orbspin(2)+1)
 
976
 
 
977
             i2 = 0
 
978
             do i1=min(orbindex(3),orbindex(7)),
 
979
     1 max(orbindex(3),orbindex(7))
 
980
               i2 = i2 + 1
 
981
               if(i2.lt.abs(orbindex(3)-orbindex(7))) then
 
982
                  if(iocc(i1+1,iref,orbspin(3)+1).eq.1) then
 
983
                     dist=dist+1
 
984
                  endif
 
985
               endif
 
986
             enddo
 
987
 
 
988
         ioccnew(orbindex(3),orbspin(3)+1) = iocc(orbindex(7),iref,
 
989
     1 orbspin(7)+1)
 
990
         ioccnew(orbindex(7),orbspin(7)+1) = iocc(orbindex(3),iref,
 
991
     1 orbspin(3)+1)
 
992
 
 
993
             i2 = 0
 
994
             do i1=min(orbindex(4),orbindex(8)),
 
995
     1 max(orbindex(4),orbindex(8))
 
996
               i2 = i2 + 1
 
997
               if(i2.lt.abs(orbindex(4)-orbindex(8))) then
 
998
                  if(iocc(i1+1,iref,orbspin(4)+1).eq.1) then
 
999
                     dist=dist+1
 
1000
                  endif
 
1001
               endif
 
1002
             enddo
 
1003
 
 
1004
            dsmult = 1.0d0
 
1005
         if(mod(dist,2).eq.1) dsmult = -dsmult
 
1006
 
 
1007
c         if(nodezero)then
 
1008
c          write(6,"('T4 iref/n:',2I4,f4.1)")iref,n,dsmult
 
1009
c         endif
 
1010
c ===============
 
1011
              dbl_mb(k_heff+n-1+(iref-1)*nref) =
 
1012
     1 dbl_mb(k_r4+m-1)*dsmult
 
1013
            endif
 
1014
          goto 400
 
1015
       endif
 
1016
       enddo
 
1017
 
 
1018
 400  continue
 
1019
 
 
1020
        enddo
 
1021
        enddo
 
1022
        enddo
 
1023
        enddo
 
1024
        enddo
 
1025
        enddo
 
1026
        enddo
 
1027
        enddo
 
1028
 
 
1029
        if (.not.ma_pop_stack(l_r4))
 
1030
     1   call errquit('tce_mrcc_iface_r4: MA problem',1,MA_ERR)
 
1031
 
 
1032
      END IF
 
1033
      END IF
 
1034
      END IF
 
1035
      END IF
 
1036
      END DO
 
1037
      END DO
 
1038
      END DO
 
1039
      END DO
 
1040
      END DO
 
1041
      END DO
 
1042
      END DO
 
1043
      END DO
 
1044
 
 
1045
      endif
 
1046
 
 
1047
       endif !sub
 
1048
 
 
1049
       enddo !iref
 
1050
 
 
1051
       return
 
1052
       end
 
1053
c
 
1054
c ==============================================
 
1055
c     Diagonalize effective Hamiltonian
 
1056
c ==============================================
 
1057
c
 
1058
        subroutine tce_diagonalize_heff(rtdb)
 
1059
        implicit none
 
1060
#include "global.fh"
 
1061
#include "mafdecls.fh"
 
1062
#include "sym.fh"
 
1063
#include "util.fh"
 
1064
#include "stdio.fh"
 
1065
#include "errquit.fh"
 
1066
#include "tce.fh"
 
1067
#include "tce_mrcc.fh"
 
1068
#include "tce_main.fh"
 
1069
#include "rtdb.fh"
 
1070
#include "tcgmsg.fh"
 
1071
#include "msgtypesf.h"
 
1072
#include "msgids.fh"
 
1073
 
 
1074
c        integer nref
 
1075
        double precision heff(nref,nref)
 
1076
        double precision vl(nref,nref)
 
1077
        double precision vr(nref,nref)
 
1078
        double precision er(nref)
 
1079
        double precision ei(nref)
 
1080
        double precision work(4*nref)
 
1081
        integer info
 
1082
        integer i,j,l
 
1083
        double precision x
 
1084
        integer k
 
1085
        logical nodezero
 
1086
        integer rtdb,itarget
 
1087
c        integer nrootmuc
 
1088
        double precision dvalue,dsum
 
1089
        double precision dfin
 
1090
        character*4 ds
 
1091
        integer l_buff,k_buff
 
1092
c        integer iignore
 
1093
        double precision isum
 
1094
        integer ddblsize,inntsize
 
1095
 
 
1096
        nodezero = (ga_nodeid().eq.0)  
 
1097
        ddblsize=MA_sizeof(MT_DBL,1,MT_BYTE)
 
1098
        inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1099
 
 
1100
        call ga_sync()
 
1101
 
 
1102
        if(lusesub) then
 
1103
        do i=1,nref*nref
 
1104
          dbl_mb(k_heff+i-1) = 0.0d0
 
1105
        enddo
 
1106
         call ga_sync()
 
1107
         call ga_get(g_heff,1,nref*nref,1,1,
 
1108
     1   dbl_mb(k_heff),1)
 
1109
c         call ga_print(g_heff)
 
1110
        endif
 
1111
 
 
1112
        do i=1,nref
 
1113
         do j=1,nref
 
1114
          x = dbl_mb(k_heff+j-1+(i-1)*nref)
 
1115
          heff(j,i) = x
 
1116
         end do
 
1117
        end do
 
1118
 
 
1119
        do i=1,nref
 
1120
          er(i) = 0.0d0
 
1121
          ei(i) = 0.0d0
 
1122
          do j=1,nref
 
1123
            vl(i,j)=0.0d0
 
1124
            vr(i,j)=0.0d0
 
1125
          enddo
 
1126
        enddo
 
1127
 
 
1128
       if(nodezero.and.(nref.lt.21)) then
 
1129
c           call ma_print(dbl_mb(k_heff),nref,nref,'Heff')
 
1130
 
 
1131
        write(LuOut,"(/,'Heff',/,
 
1132
     1 '=============================================')")
 
1133
        do i=1,nref
 
1134
         write(LuOut,"(i5,i5,100F14.8)")ga_nodeid(),i,
 
1135
     1 (dbl_mb(k_heff+(j-1)*nref+i-1),j=1,nref)
 
1136
        enddo
 
1137
 
 
1138
       endif
 
1139
c       call ga_sync()
 
1140
c       call util_flush(LuOut)
 
1141
c       if((ga_nodeid().eq.5).and.(nref.lt.21)) then
 
1142
c           call ma_print(dbl_mb(k_heff),nref,nref,'Heff')
 
1143
c
 
1144
c        write(LuOut,"(/,'Heff 5',/,
 
1145
c     1 '=============================================')")
 
1146
c        do i=1,nref
 
1147
c         write(LuOut,"(i5,i5,100F14.8)")ga_nodeid(),i,
 
1148
c     1 (dbl_mb(k_heff+(j-1)*nref+i-1),j=1,nref)
 
1149
c        enddo
 
1150
 
 
1151
c       endif
 
1152
c       call util_flush(LuOut)
 
1153
c       call ga_sync()
 
1154
 
 
1155
c         call util_flush(LuOut)
 
1156
c        write(6,*)'BEFORE',ga_nodeid()
 
1157
c         call util_flush(LuOut)
 
1158
        call ga_sync()
 
1159
c       if(nodezero)write(6,*)'TEST 3'
 
1160
        if(nodezero) then
 
1161
c        call DGEEV('V','V',nref,heff,nref,er,ei,vl,nref,vr,
 
1162
        call util_dgeev('V','V',nref,heff,nref,er,ei,vl,nref,vr,
 
1163
     $                  nref,work,4*nref,info)
 
1164
        if(info .ne. 0) call errquit('Heff diagonalization',0,CALC_ERR)
 
1165
        call amp_stabilization(vl,vr,nref)
 
1166
        endif
 
1167
c      if(nodezero)write(6,*)'TEST 4'
 
1168
        call ga_sync()
 
1169
 
 
1170
c         call util_flush(LuOut)
 
1171
c        write(6,*)'AFTER',ga_nodeid()
 
1172
c         call util_flush(LuOut)
 
1173
         
 
1174
c       if(nodezero.and..not.lconverged) then
 
1175
c           call ma_print(dbl_mb(k_heff),nref,nref,'Heff')
 
1176
c       endif
 
1177
 
 
1178
c       if(lconverged.and.nodezero) then
 
1179
       if(nodezero) then
 
1180
        if(nref.lt.21) then
 
1181
 
 
1182
        write(LuOut,"(/,'Eigenvalues (real and imaginary)',/,
 
1183
     1 '=============================================')")
 
1184
        do i=1,nref
 
1185
           write(LuOut,"(F18.12,100F14.8)")er(i),ei(i)
 
1186
        enddo
 
1187
 
 
1188
        write(LuOut,"(/,'Left eigenvectors',/,
 
1189
     1 '=============================================')")
 
1190
        do i=1,nref
 
1191
           write(LuOut,"(i5,100F14.8)")i,(vl(i,j),j=1,nref)
 
1192
        enddo
 
1193
 
 
1194
c       write(LuOut,"(/,'Left eigenvectors - squares',/,
 
1195
c     1 '=============================================')")
 
1196
c        do i=1,nref
 
1197
c           write(LuOut,"(i5,35f18.12)")i,((vl(i,j)*vl(i,j)),j=1,nref)
 
1198
c        enddo
 
1199
        
 
1200
c        endif
 
1201
c        endif      
 
1202
 
 
1203
        write(LuOut,"(/,'Right eigenvectors',/,
 
1204
     1 '=============================================')")
 
1205
        do i=1,nref
 
1206
           write(LuOut,"('VR',i5,100F14.8)")i,(vr(i,j),j=1,nref)
 
1207
        enddo
 
1208
 
 
1209
        endif
 
1210
        endif
 
1211
 
 
1212
        do i=1,nref
 
1213
         do j=1,nref
 
1214
          dbl_mb(k_sqc+(i-1)*nref+j-1)=vr(i,j)
 
1215
          dbl_mb(k_sqcl+(i-1)*nref+j-1)=vl(i,j)
 
1216
         enddo
 
1217
        enddo
 
1218
 
 
1219
       call ga_brdcst(Msg_Vec_EVal+21,dbl_mb(k_sqc),
 
1220
     1 ddblsize*nref*nref, 0)
 
1221
       call ga_sync()
 
1222
       call ga_brdcst(Msg_Vec_EVal+20,dbl_mb(k_sqcl),
 
1223
     1 ddblsize*nref*nref, 0)
 
1224
       call ga_sync()
 
1225
 
 
1226
       if (.not.rtdb_get(rtdb,'mrcc:zignore',mt_int,1,iignore))
 
1227
     1 iignore = 0
 
1228
 
 
1229
        if(.not.nodezero) then
 
1230
        do i=1,nref
 
1231
         do j=1,nref
 
1232
          vr(i,j) = dbl_mb(k_sqc+(i-1)*nref+j-1)
 
1233
          vl(i,j) = dbl_mb(k_sqcl+(i-1)*nref+j-1)
 
1234
         enddo
 
1235
        enddo
 
1236
        endif
 
1237
 
 
1238
c--------
 
1239
 
 
1240
        epsilon = 0.0d0
 
1241
        do i=1,nref
 
1242
           isum = 0.0d0
 
1243
           do j=1,nref
 
1244
             isum = isum + vr(j,i)
 
1245
           enddo
 
1246
           if((iignore.eq.0).or.(abs(isum).gt.1.0d-5)) then
 
1247
             if(epsilon.gt.min(epsilon,er(i)))mkroot=i
 
1248
             epsilon = min(epsilon,er(i))
 
1249
           endif
 
1250
        enddo
 
1251
 
 
1252
c--------
 
1253
 
 
1254
       if (.not.rtdb_get(rtdb,'mrcc:rootmuc',mt_int,1,nrootmuc))
 
1255
     1 nrootmuc = 0
 
1256
c ------
 
1257
      if(nrootmuc.gt.0) then ! 1
 
1258
        dfin = 0.0d0
 
1259
 
 
1260
        do j=1,nref
 
1261
        dsum = 0.0d0
 
1262
        isum = 0.0d0
 
1263
        do i=1,nref
 
1264
           isum = isum + vr(i,j)
 
1265
        enddo
 
1266
        if((iignore.eq.0).or.(abs(isum).gt.1.0d-5)) then
 
1267
        do i=1,nref
 
1268
          write(ds,"(I3.3)")i
 
1269
       if (.not.rtdb_get(rtdb,'mrcc:rootmuc'//ds,mt_dbl,1,dvalue))
 
1270
     1 dvalue = 0.0d0
 
1271
          if(dvalue.lt.0.0d0) then
 
1272
            if(abs(vr(i,j)).lt.abs(dvalue))
 
1273
     1     dsum = dsum + abs(dvalue)*abs(vr(i,j))
 
1274
          else
 
1275
           dsum = dsum + abs(dvalue)*abs(vr(i,j))
 
1276
          endif
 
1277
        enddo
 
1278
c            write(6,"('SUM:',I4,4F16.12)")j,
 
1279
c     1 abs(dsum)
 
1280
          if(dfin.lt.abs(dsum)) then
 
1281
c            write(6,"('I am watching reference #',I4,4F16.12)")j
 
1282
            mkroot=j
 
1283
            epsilon = er(j)
 
1284
            dfin = abs(dsum)
 
1285
          endif
 
1286
        endif
 
1287
        enddo
 
1288
 
 
1289
      else ! 1
 
1290
 
 
1291
      if (rtdb_get(rtdb,'bwcc:targetroot',mt_int,1,itarget)) then
 
1292
c       mkroot = itarget
 
1293
       do i=1,nref
 
1294
 
 
1295
           isum = 0.0d0
 
1296
           do j=1,nref
 
1297
c             if(abs(vr(j,i)).lt.1.0d-8) 
 
1298
              isum = isum + vr(j,i)
 
1299
           enddo
 
1300
           if((iignore.eq.0).or.(abs(isum).gt.1.0d-5)) then
 
1301
 
 
1302
        dvalue = er(i)
 
1303
        k = 0     
 
1304
        do j=1,nref
 
1305
        if(j.ne.i) then
 
1306
 
 
1307
           isum = 0.0d0
 
1308
           do l=1,nref
 
1309
c             if(abs(vr(l,j)).lt.1.0d-8) 
 
1310
             isum = isum + vr(l,j)
 
1311
           enddo
 
1312
           if((iignore.eq.0).or.(abs(isum).gt.1.0d-5)) then
 
1313
 
 
1314
 
 
1315
         if(er(j).lt.er(i))k=k+1
 
1316
 
 
1317
           endif
 
1318
 
 
1319
        endif
 
1320
        enddo
 
1321
        if(k.eq.(itarget-1)) then
 
1322
         itarget = i
 
1323
         goto 63454
 
1324
        endif
 
1325
 
 
1326
            endif
 
1327
 
 
1328
       enddo
 
1329
 
 
1330
63454  continue
 
1331
 
 
1332
         mkroot = itarget
 
1333
         epsilon = er(itarget)
 
1334
 
 
1335
      endif
 
1336
      endif
 
1337
 
 
1338
        if (nodezero.and.(nref.lt.21)) then
 
1339
           write(6,"('Target root: ',I4)")mkroot
 
1340
        end if
 
1341
 
 
1342
ckbn    Introduce some checks before proceeding further.
 
1343
 
 
1344
ckbn    check of mkroot, crash if mkroot gt nref
 
1345
        if(nodezero) then 
 
1346
         if(mkroot.gt.nref) call errquit
 
1347
     +    ('root followed is greater than total number of references',0,
 
1348
     +    CALC_ERR)
 
1349
        endif
 
1350
 
 
1351
 
 
1352
ckbn check whether there is imaginary eigen values
 
1353
        do i=1,nref
 
1354
c         write(LuOut,'(A,F17.10)')'Imaginary eigenvalue',ei(i)
 
1355
c         if (nodezero) call util_flush(LuOut)
 
1356
         if(abs(ei(i)).ge.toleiimag) then
 
1357
          write(LuOut,'(A,F15.10)')
 
1358
     +            'Warning: complex Heff eigenvalue detected',ei(i)
 
1359
          if(i.eq.mkroot) then
 
1360
           write(LuOut,*) "ignorecomplex1 ", ignorecomplex
 
1361
c          if (rtdb_get(rtdb,'mrcc:ignorecomplex',mt_log,1,
 
1362
c     +     ignorecomplex)) ignorecomplex = .true.
 
1363
c           call errquit('Complex root found and no ignorecomplex option',
 
1364
c     +      0,RTDB_ERR)
 
1365
c          else
 
1366
c           if(nodezero) write(*,*) "ignorecomplex1. ", ignorecomplex
 
1367
c           ignorecomplex = .true.
 
1368
c          endif
 
1369
c           if(nodezero) write(*,*) "ignorecomplex2 ", ignorecomplex
 
1370
           if(ignorecomplex) then
 
1371
            write(LuOut,'(A,F15.10)')
 
1372
     +       'Warning: Proceeding with complex Heff eigenvalue ',ei(i)
 
1373
           else
 
1374
            call errquit('Warning:complex Heff eigenvalue detected',0,
 
1375
     +                    UNKNOWN_ERR)
 
1376
           endif
 
1377
          endif
 
1378
         endif
 
1379
        enddo
 
1380
 
 
1381
 
 
1382
 
 
1383
 
 
1384
        if(nref.gt.20) then
 
1385
         if(nodezero)write(LuOut,"(/,'Right eigenvector for target',I7,
 
1386
     1 /)")mkroot
 
1387
c      if (.not.ma_push_get(mt_dbl,nref,'buff',l_buff,k_buff))
 
1388
c     1   call errquit('tce_mrcc_iface_buff: MA problem',0,MA_ERR)
 
1389
         do i=1,nref
 
1390
          if(abs(vr(i,mkroot)).gt.0.05d0) then
 
1391
            if(nodezero)write(LuOut,"(I7,'  ',F11.8)")i,vr(i,mkroot)
 
1392
          endif
 
1393
         enddo
 
1394
c        if (.not.ma_pop_stack(l_buff))
 
1395
c     1   call errquit('tce_mrcc_iface_buff: MA problem',1,MA_ERR)
 
1396
 
 
1397
        endif
 
1398
 
 
1399
        if (nodezero) call util_flush(LuOut)
 
1400
 
 
1401
cjb broadcasts
 
1402
 
 
1403
        call ga_brdcst(Msg_Vec_EVal+MSGINT+30,mkroot,inntsize, 0)
 
1404
        call ga_brdcst(Msg_Vec_EVal+MSGDBL+31,epsilon,ddblsize, 0)
 
1405
        call ga_sync()
 
1406
 
 
1407
c        write(LuOut,"(/,'Epsilon: ',2F16.12,/)") epsilon
 
1408
 
 
1409
        return
 
1410
        end
 
1411
 
 
1412
c
 
1413
c ==============================================
 
1414
c     Clean internal amplitudes
 
1415
c ==============================================
 
1416
c
 
1417
        subroutine tce_internal_t_zero(d_t1m,d_t2m,k_t1_offsetm,
 
1418
     1 k_t2_offsetm,lneedt3,d_t3m,k_t3_offsetm,rtdb)
 
1419
c     1 k_t2_offsetm,nref,lneedt3,d_t3m,k_t3_offsetm,rtdb)
 
1420
        implicit none
 
1421
#include "global.fh"
 
1422
#include "mafdecls.fh"
 
1423
#include "sym.fh"
 
1424
#include "util.fh"
 
1425
#include "stdio.fh"
 
1426
#include "errquit.fh"
 
1427
#include "tce.fh"
 
1428
#include "tce_mrcc.fh"
 
1429
#include "tce_main.fh"
 
1430
#include "rtdb.fh"
 
1431
 
 
1432
        integer d_t1m(maxref),d_t2m(maxref)
 
1433
        integer d_t3m(maxref)
 
1434
        integer k_t1_offsetm(maxref),k_t2_offsetm(maxref)
 
1435
        integer k_t3_offsetm(maxref)
 
1436
        integer size,p5b,h6b
 
1437
c        integer nref
 
1438
        integer l_t1,k_t1
 
1439
        integer l_t2,k_t2
 
1440
        integer l_t3,k_t3
 
1441
        integer i,j,k,l,m
 
1442
        integer iref,rtdb
 
1443
        integer p1b,p2b,h3b,h4b
 
1444
        integer p3b,p4b,h5b,p6b
 
1445
        integer h1b,h2b
 
1446
c        logical lneedt3,limprovet
 
1447
        logical lneedt3
 
1448
        integer orbindex(6),orbspin(6)
 
1449
        EXTERNAL NXTASKsub
 
1450
        EXTERNAL NXTASK
 
1451
        INTEGER NXTASKsub
 
1452
        INTEGER NXTASK
 
1453
        INTEGER nxt
 
1454
        INTEGER nprocs
 
1455
        INTEGER count
 
1456
 
 
1457
c       if (.not.rtdb_get(rtdb,'mrcc:improvetiling',mt_log,1,limprovet))
 
1458
c     1 limprovet = .false.
 
1459
 
 
1460
      do iref=1,nref
 
1461
 
 
1462
        if((int_mb(k_refafi+iref-1).eq.int_mb(k_innodes
 
1463
     1 +ga_nnodes()+ga_nodeid())).or.(.not.lusesub)) then
 
1464
 
 
1465
         k_sym = k_symm(iref)
 
1466
         k_offset = k_offsetm(iref)
 
1467
         k_range = k_rangem(iref)
 
1468
         k_spin = k_spinm(iref)
 
1469
         k_movecs_sorted = k_movecs_sortedm(iref)
 
1470
 
 
1471
         noa = nblcks(1,iref)
 
1472
         nob = nblcks(2,iref)
 
1473
         nva = nblcks(3,iref)
 
1474
         nvb = nblcks(4,iref)
 
1475
 
 
1476
         noab = noa+nob
 
1477
         nvab = nva+nvb
 
1478
 
 
1479
         do p5b = noab+1,noab+nvab
 
1480
         do h6b = 1,noab
 
1481
 
 
1482
      if (int_mb(k_spin+p5b-1) .eq. int_mb(k_spin+h6b-1)) then
 
1483
      if (ieor(int_mb(k_sym+p5b-1),int_mb(k_sym+h6b-1)).eq.irrep_t)then
 
1484
      if ((.not.restricted).or.(int_mb(k_spin+p5b-1)+int_mb(k_spin+h6b-1
 
1485
     &).ne.4)) then
 
1486
      if(log_mb(k_isactive(iref)+p5b-1).and.
 
1487
     &log_mb(k_isactive(iref)+h6b-1)) then
 
1488
 
 
1489
        size = int_mb(k_range+h6b-1) * int_mb(k_range+p5b-1)
 
1490
        if (.not.ma_push_get(mt_dbl,size,'t1',l_t1,k_t1))
 
1491
     1   call errquit('tce_mrcc_iface_t1: MA problem',0,MA_ERR)
 
1492
 
 
1493
c         call dfill(size,0.0d0,dbl_mb(k_t1),1)
 
1494
         do i=1,size
 
1495
          dbl_mb(k_t1+i-1)=0.0d0
 
1496
         enddo
 
1497
 
 
1498
cjb ============================
 
1499
 
 
1500
         if(limprovet) then
 
1501
 
 
1502
        call get_hash_block(d_t1m(iref),dbl_mb(k_t1),size,
 
1503
     1   int_mb(k_t1_offsetm(iref)),h6b-1+noab*(p5b-noab-1))
 
1504
 
 
1505
        k=0
 
1506
        do i=1,int_mb(k_range+p5b-1)
 
1507
        do j=1,int_mb(k_range+h6b-1)
 
1508
          k = k + 1
 
1509
 
 
1510
        orbspin(1) = int_mb(k_spin+p5b-1)-1
 
1511
        orbspin(2) = int_mb(k_spin+h6b-1)-1
 
1512
 
 
1513
        orbindex(1) = (1 - orbspin(1)+
 
1514
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p5b-1)+i-1))/2
 
1515
        orbindex(2) = (1 - orbspin(2)+
 
1516
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h6b-1)+j-1))/2
 
1517
 
 
1518
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
1519
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
1520
 
 
1521
        if(isactive(orbindex(1),orbspin(1)+1).and.
 
1522
     1 isactive(orbindex(2),orbspin(2)+1).or.(.not.limprovet)) then
 
1523
 
 
1524
          dbl_mb(k_t1+k-1) = 0.0d0
 
1525
 
 
1526
          endif
 
1527
          enddo
 
1528
          enddo
 
1529
 
 
1530
          endif ! limprovet
 
1531
 
 
1532
c =============================
 
1533
 
 
1534
c         do i=1,size
 
1535
c           dbl_mb(k_t1+i-1) = 0.0d0
 
1536
c         enddo
 
1537
 
 
1538
         call put_hash_block(d_t1m(iref),dbl_mb(k_t1),size,
 
1539
     1   int_mb(k_t1_offsetm(iref)),((p5b-noab-1)*noab+h6b-1))
 
1540
 
 
1541
        if (.not.ma_pop_stack(l_t1))
 
1542
     1   call errquit('tce_mrcc_iface_t1: MA problem',1,MA_ERR)
 
1543
 
 
1544
      endif
 
1545
      endif
 
1546
      endif
 
1547
      endif
 
1548
      enddo
 
1549
      enddo
 
1550
 
 
1551
cjb new in parallel
 
1552
 
 
1553
      nxt = 0
 
1554
      if(limprovet) then
 
1555
       if(lusesub) then
 
1556
         nprocs=GA_pgroup_NNODES(mypgid)
 
1557
         nxt=NXTASKsub(nprocs,1,mypgid)
 
1558
       else
 
1559
         nprocs=GA_NNODES()
 
1560
         nxt=NXTASK(nprocs,1)
 
1561
       endif
 
1562
      count = 0
 
1563
      endif
 
1564
 
 
1565
      DO p1b = noab+1,noab+nvab
 
1566
      DO p2b = p1b,noab+nvab
 
1567
      DO h3b = 1,noab
 
1568
      DO h4b = h3b,noab
 
1569
 
 
1570
      IF ((nxt.eq.count).or.(.not.limprovet)) THEN
 
1571
 
 
1572
      IF (int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1) .eq. int_mb(k_spin+h
 
1573
     &3b-1)+int_mb(k_spin+h4b-1)) THEN
 
1574
      IF (ieor(int_mb(k_sym+p1b-1),ieor(int_mb(k_sym+p2b-1),ieor(int_mb(
 
1575
     &k_sym+h3b-1),int_mb(k_sym+h4b-1)))) .eq. irrep_t) THEN
 
1576
      IF ((.not.restricted).or.(int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1
 
1577
     &)+int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1).ne.8)) THEN
 
1578
      if(log_mb(k_isactive(iref)+p1b-1).and.
 
1579
     1 log_mb(k_isactive(iref)+p2b-1).and.
 
1580
     2 log_mb(k_isactive(iref)+h3b-1).and.
 
1581
     3 log_mb(k_isactive(iref)+h4b-1)) then
 
1582
 
 
1583
      size = int_mb(k_range+p1b-1) * int_mb(k_range+p2b-1) * int_
 
1584
     &mb(k_range+h3b-1) * int_mb(k_range+h4b-1)
 
1585
 
 
1586
        if (.not.ma_push_get(mt_dbl,size,'t2',l_t2,k_t2))
 
1587
     1   call errquit('tce_mrcc_iface_t2: MA problem',0,MA_ERR)
 
1588
 
 
1589
c         call dfill(size,0.0d0,dbl_mb(k_t2),1)
 
1590
         do i=1,size
 
1591
          dbl_mb(k_t2+i-1)=0.0d0
 
1592
         enddo
 
1593
 
 
1594
c ===============================================================
 
1595
 
 
1596
         if(limprovet) then
 
1597
 
 
1598
        call get_hash_block(d_t2m(iref),dbl_mb(k_t2),size,
 
1599
     1   int_mb(k_t2_offsetm(iref)),h4b-1+noab*(h3b-1+noab*(p2b-
 
1600
     &noab-1+nvab*(p1b - noab - 1))))
 
1601
 
 
1602
       m = 0
 
1603
        do i=1,int_mb(k_range+p1b-1)
 
1604
        do j=1,int_mb(k_range+p2b-1)
 
1605
        do k=1,int_mb(k_range+h3b-1)
 
1606
        do l=1,int_mb(k_range+h4b-1)
 
1607
         m = m + 1
 
1608
 
 
1609
        orbspin(1) = int_mb(k_spin+p1b-1)-1
 
1610
        orbspin(2) = int_mb(k_spin+p2b-1)-1
 
1611
        orbspin(3) = int_mb(k_spin+h3b-1)-1
 
1612
        orbspin(4) = int_mb(k_spin+h4b-1)-1
 
1613
 
 
1614
        orbindex(1) = (1 - orbspin(1)+
 
1615
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p1b-1)+i-1))/2
 
1616
        orbindex(2) = (1 - orbspin(2)+
 
1617
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+p2b-1)+j-1))/2
 
1618
        orbindex(3) = (1 - orbspin(3)+
 
1619
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h3b-1)+k-1))/2
 
1620
        orbindex(4) = (1 - orbspin(4)+
 
1621
     1 int_mb(k_mo_indexm(iref)+int_mb(k_offset+h4b-1)+l-1))/2
 
1622
 
 
1623
        orbindex(1) = moindexes(orbindex(1),orbspin(1)+1,iref)
 
1624
        orbindex(2) = moindexes(orbindex(2),orbspin(2)+1,iref)
 
1625
        orbindex(3) = moindexes(orbindex(3),orbspin(3)+1,iref)
 
1626
        orbindex(4) = moindexes(orbindex(4),orbspin(4)+1,iref)
 
1627
 
 
1628
        if(isactive(orbindex(1),orbspin(1)+1).and.
 
1629
     1 isactive(orbindex(2),orbspin(2)+1).and.
 
1630
     2 isactive(orbindex(3),orbspin(3)+1).and.
 
1631
     3 isactive(orbindex(4),orbspin(4)+1).or.(.not.limprovet)) then
 
1632
 
 
1633
         dbl_mb(k_t2+m-1)=0.0d0
 
1634
 
 
1635
          endif
 
1636
         enddo
 
1637
         enddo
 
1638
         enddo
 
1639
         enddo
 
1640
 
 
1641
         endif
 
1642
 
 
1643
c ===============================================================
 
1644
 
 
1645
c        do i=1,size
 
1646
c           write(LuOut,*)dbl_mb(k_t2+i-1),'->','0.00000000'
 
1647
c           dbl_mb(k_t2+i-1) = 0.0d0
 
1648
c        enddo
 
1649
 
 
1650
        call put_hash_block(d_t2m(iref),dbl_mb(k_t2),size,
 
1651
     1   int_mb(k_t2_offsetm(iref)),((((p1b-noab-1)*nvab+p2b-noab-1)
 
1652
     2   *noab+h3b-1)*noab+h4b-1))
 
1653
 
 
1654
        if (.not.ma_pop_stack(l_t2))
 
1655
     1   call errquit('tce_mrcc_iface_t2: MA problem',1,MA_ERR)
 
1656
 
 
1657
      END IF
 
1658
      END IF
 
1659
      END IF
 
1660
      endif
 
1661
 
 
1662
      if(limprovet) then
 
1663
       if(lusesub) then
 
1664
        nxt=NXTASKsub(nprocs,1,mypgid)
 
1665
       else
 
1666
        nxt=NXTASK(nprocs,1)
 
1667
       endif
 
1668
      endif
 
1669
 
 
1670
      endif
 
1671
 
 
1672
       if(limprovet)count = count + 1
 
1673
 
 
1674
      END DO
 
1675
      END DO
 
1676
      END DO
 
1677
      END DO
 
1678
 
 
1679
      if(limprovet) then
 
1680
      if(lusesub) then
 
1681
      nxt=NXTASKsub(-nprocs,1,mypgid)
 
1682
      call GA_Pgroup_SYNC(mypgid)
 
1683
      else
 
1684
      nxt=NXTASKsub(-nprocs,1)
 
1685
      call GA_SYNC()
 
1686
      endif
 
1687
      endif
 
1688
 
 
1689
      if(lneedt3) then
 
1690
 
 
1691
      DO p2b = noab+1,noab+nvab
 
1692
      DO p3b = p2b,noab+nvab
 
1693
      DO p4b = p3b,noab+nvab
 
1694
      DO h1b = 1,noab
 
1695
      DO h5b = h1b,noab
 
1696
      DO h6b = h5b,noab
 
1697
      IF (int_mb(k_spin+p2b-1)+int_mb(k_spin+p3b-1)+int_mb(k_spin+p4b-1)
 
1698
     & .eq. int_mb(k_spin+h1b-1)+int_mb(k_spin+h5b-1)+int_mb(k_spin+h6b-
 
1699
     &1)) THEN
 
1700
      IF (ieor(int_mb(k_sym+p2b-1),ieor(int_mb(k_sym+p3b-1),ieor(int_mb(
 
1701
     &k_sym+p4b-1),ieor(int_mb(k_sym+h1b-1),ieor(int_mb(k_sym+h5b-1),int
 
1702
     &_mb(k_sym+h6b-1)))))) .eq. irrep_t) THEN
 
1703
      IF ((.not.restricted).or.(int_mb(k_spin+p2b-1)+int_mb(k_spin+p3b-1
 
1704
     &)+int_mb(k_spin+p4b-1)+int_mb(k_spin+h1b-1)+int_mb(k_spin+h5b-1)+i
 
1705
     &nt_mb(k_spin+h6b-1).ne.12)) THEN
 
1706
      IF ((log_mb(k_isactive(iref)+p2b-1).eqv..true.).and.
 
1707
     1 (log_mb(k_isactive(iref)+p3b-1).eqv..true.).and.
 
1708
     2 (log_mb(k_isactive(iref)+p4b-1).eqv..true.).and.
 
1709
     3 (log_mb(k_isactive(iref)+h1b-1).eqv..true.).and.
 
1710
     4 (log_mb(k_isactive(iref)+h5b-1).eqv..true.).and.
 
1711
     5 (log_mb(k_isactive(iref)+h6b-1).eqv..true.)) THEN
 
1712
 
 
1713
      size = int_mb(k_range+p2b-1) * int_mb(k_range+p3b-1) * int_
 
1714
     &mb(k_range+p4b-1) * int_mb(k_range+h1b-1) * int_mb(k_range+h5b-1)
 
1715
     &* int_mb(k_range+h6b-1)
 
1716
 
 
1717
        if (.not.ma_push_get(mt_dbl,size,'t3',l_t3,k_t3))
 
1718
     1   call errquit('tce_mrcc_iface_t3: MA problem',0,MA_ERR)
 
1719
 
 
1720
        do i=1,size
 
1721
           dbl_mb(k_t3+i-1) = 0.0d0
 
1722
        enddo
 
1723
        call put_hash_block(d_t3m(iref),dbl_mb(k_t3),size,
 
1724
     1   int_mb(k_t3_offsetm(iref)),(h6b - 1 + noab * 
 
1725
     2 (h5b - 1 + noab * (h1b
 
1726
     &- 1 + noab * (p4b - noab - 1 + nvab * (p3b - noab - 1 + nvab * (p2
 
1727
     &b - noab - 1)))))))
 
1728
 
 
1729
        if (.not.ma_pop_stack(l_t3))
 
1730
     1   call errquit('tce_mrcc_iface_t3: MA problem',1,MA_ERR)
 
1731
 
 
1732
      END IF
 
1733
      END IF
 
1734
      END IF
 
1735
      END IF
 
1736
      END DO
 
1737
      END DO
 
1738
      END DO
 
1739
      END DO
 
1740
      END DO
 
1741
      END DO
 
1742
 
 
1743
      endif !needt3
 
1744
 
 
1745
c      write(6,"('CPU BEFORE',I4,F16.12)")ga_nodeid(),util_cpusec()
 
1746
c      call ga_sync()
 
1747
      if(lusesub) call ga_pgroup_sync(
 
1748
     1 int_mb(k_innodes+ga_nnodes()+ga_nodeid()))
 
1749
c      write(6,"('CPU AFTER',I4,F16.12)")ga_nodeid(),util_cpusec()
 
1750
 
 
1751
      endif !sub
 
1752
 
 
1753
      enddo !iref
 
1754
 
 
1755
        return
 
1756
        end
 
1757
 
 
1758