← Back to branch summary

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

« back to all changes in this revision

Viewing changes to .pc/09_backported_6.1.1_fixes.patch/src/mp2_grad/mp2_back_transform.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:
.pc/09_backported_6.1.1_fixes.patch/src/mp2_grad/mp2_back_transform.F
1
 
      subroutine mp2_nonsep_uhf( rtdb, geom,
2
 
     $     basis, oskel,
3
 
     $     nbf, 
4
 
     $     nir, 
5
 
     $     noa, nva, nva_lo, nva_hi, num_va, num_oa, 
6
 
     $     nob, nvb, nvb_lo, nvb_hi, num_vb, num_ob,
7
 
     $     sym_lo_oa, sym_hi_oa, sym_lo_va, sym_hi_va,
8
 
     $     sym_lo_ob, sym_hi_ob, sym_lo_vb, sym_hi_vb,
9
 
     $     oseg_lo, oseg_hi, 
10
 
     $     irs_a, irs_b,
11
 
     $     c_a, c_b,            ! Better left in a GA to conserve memory
12
 
     $     nva_lo_local, nva_hi_local, 
13
 
     $     tunita, tunitb, grad,
14
 
     P     p_file_size)
15
 
*
16
 
* $Id: mp2_back_transform.F 21424 2011-11-07 19:21:05Z edo $
17
 
*
18
 
      implicit none
19
 
#include "errquit.fh"
20
 
#include "mafdecls.fh"
21
 
#include "global.fh"
22
 
#include "util.fh"
23
 
#include "bas.fh"
24
 
#include "eaf.fh"
25
 
#include "cmp2ps.fh"
26
 
#include "geom.fh"
27
 
#include "rtdb.fh"
28
 
      integer basis, rtdb, geom
29
 
      logical oskel
30
 
      integer nbf               ! No. of basis functions
31
 
      integer nir               ! No. of irreducible representations
32
 
      integer noa, nob          ! No. of occupied orbitals
33
 
      integer nva, nva_lo, nva_hi ! Number and ranges of virtual orbitals
34
 
      integer nvb, nvb_lo, nvb_hi
35
 
      integer num_va(0:nir-1), num_vb(0:nir-1) ! No. of vir of each symmetry
36
 
      integer num_oa(0:nir-1), num_ob(0:nir-1) ! No. of occ of each symmetry
37
 
      integer sym_lo_oa(0:nir-1), sym_hi_oa(0:nir-1),
38
 
     $     sym_lo_va(0:nir-1), sym_hi_va(0:nir-1) ! Ranges of each symmetry
39
 
      integer sym_lo_ob(0:nir-1), sym_hi_ob(0:nir-1),
40
 
     $     sym_lo_vb(0:nir-1), sym_hi_vb(0:nir-1)
41
 
      integer oseg_lo, oseg_hi  ! Range of occupied for this pass
42
 
      integer irs_a(nbf), irs_b(nbf) ! Orbital symmetries
43
 
      double precision c_a(nbf,*), c_b(nbf,*) ! MO coefficients
44
 
      integer nva_lo_local, nva_hi_local ! Range of virtuals on this node
45
 
      integer tunita, tunitb    ! Unit no.s for pure and mixed spin T
46
 
      double precision grad(3,*)
47
 
c
48
 
c     Allocate memory for back transformation routine
49
 
c
50
 
      integer g_buf
51
 
      integer nshpair, nshpairlocal, nbfpair, nbfpairlocal
52
 
      integer l_shpairs, k_shpairs, l_shpairslocal, k_shpairslocal
53
 
      integer l_shdim, k_shdim, l_shlo, k_shlo
54
 
      integer l_t, k_t, l_tmp, k_tmp, l_iauv, k_iauv, l_map, k_map
55
 
      integer l_c_t, k_c_t
56
 
      integer l_act, k_act, l_actsh, k_actsh
57
 
      integer twopdmunit
58
 
      integer junk, ninseg, ierr, i, j, ish, ishlo, ishhi, shmax, tdim
59
 
      integer nsh, natoms, nactive, nblock
60
 
c
61
 
      integer k_scr, l_scr, k_lab, l_lab, k_eri, l_eri, leneri, lenscr
62
 
      integer k_bftosh, l_bftosh, k_bftoce, l_bftoce, k_pdm, l_pdm
63
 
      integer eaftype,eaf_size_in_mb,inntsize
64
 
      double precision p_file_size
65
 
c
66
 
      character*(nw_max_path_len) fname
67
 
c
68
 
      double precision tol2e
69
 
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
70
 
c
71
 
      if (.not. rtdb_get(rtdb, 'mp2:backtol', mt_dbl, 1, tol2e))
72
 
     $     tol2e = 1d-9
73
 
c
74
 
      if (.not. bas_numcont(basis, nsh)) call errquit
75
 
     $     ('mp2: backt bad basis handle ', basis, BASIS_ERR)
76
 
      ninseg = oseg_hi - oseg_lo + 1
77
 
c
78
 
c     Make mapping of shell pairs to processors
79
 
c
80
 
      call int_init(rtdb,1,basis)     ! intd_init overwrite mem estimates
81
 
      call schwarz_init(geom,basis)
82
 
      call int_terminate()
83
 
c
84
 
      call mp2_backt_info(basis, tol2e, oskel,
85
 
     $     .false., nshpair, nshpairlocal,
86
 
     $     nbfpair, nbfpairlocal, shmax, junk, junk, junk, junk, junk)
87
 
      if (.not. ma_push_get(mt_int, nshpair*2, 'mp2: shpairs',
88
 
     $     l_shpairs, k_shpairs)) call errquit
89
 
     $     ('mp2: insufficient memory : shpairs ', 2*nshpair, MA_ERR)
90
 
      if (.not. ma_push_get(mt_int, nshpairlocal*2, 'mp2:shpairslocal',
91
 
     $     l_shpairslocal, k_shpairslocal)) call errquit
92
 
     $     ('mp2: insufficient memory : shpairs ', 2*nshpairlocal,
93
 
     &       MA_ERR)
94
 
      if (.not. ma_push_get(mt_int, ga_nnodes(), 'mp2:map',
95
 
     $     l_map, k_map)) call errquit
96
 
     $     ('mp2: insufficient memory : map ', ga_nnodes(), MA_ERR)
97
 
      if (.not. ma_push_get(mt_int, nsh, 'mp2:shlo',
98
 
     $     l_shlo, k_shlo)) call errquit
99
 
     $     ('mp2: insufficient memory : shlo ', nsh, MA_ERR)
100
 
      if (.not. ma_push_get(mt_int, nsh, 'mp2:shdim',
101
 
     $     l_shdim, k_shdim)) call errquit
102
 
     $     ('mp2: insufficient memory : shdim ', nsh, MA_ERR)
103
 
      call mp2_backt_info(basis, tol2e, oskel,
104
 
     $     .true., nshpair, nshpairlocal,
105
 
     $     nbfpair, nbfpairlocal, shmax,
106
 
     $     int_mb(k_shdim), int_mb(k_shlo),
107
 
     $     int_mb(k_shpairs), int_mb(k_shpairslocal), int_mb(k_map))
108
 
c
109
 
c     Allocate remaining memory and open file to hold the 
110
 
c     3-parts transformed density matrix.
111
 
c
112
 
      tdim = max(nbf*nbf, nbf*shmax*shmax)
113
 
c         
114
 
      if (.not. ma_push_get(mt_dbl, tdim, 'mp2: backt t',
115
 
     $     l_t, k_t)) call errquit
116
 
     $     ('mp2: insufficient memory : t ', tdim, MA_ERR)
117
 
      if (.not. ma_push_get(mt_dbl, tdim, 'mp2: backt tmp',
118
 
     $     l_tmp, k_tmp)) call errquit
119
 
     $     ('mp2: insufficient memory : tmp ', tdim, MA_ERR)
120
 
      if (.not. ma_push_get(mt_dbl, nbf*nbf, 'mp2: backt iauv',
121
 
     $     l_iauv, k_iauv)) call errquit
122
 
     $     ('mp2: insufficient memory : iauv ', nbf*nbf, MA_ERR)
123
 
c
124
 
c     Right now we have an N**3 array in GA.  This is easily
125
 
c     replaced with an N*B array with B=P*S*S (P procs, S=max shell size)
126
 
c     by blocking the transformation.
127
 
c
128
 
c     For small problems some nodes may not have data in which case
129
 
c     map(proc) will be set to a value larger then nbfpair (which is the 
130
 
c     actual maximum dimension) ... discard these before creating
131
 
c     the array.
132
 
c
133
 
      nblock = 0
134
 
      do i = 1, ga_nnodes()
135
 
         if (int_mb(k_map+i-1) .le. nbfpair) nblock = nblock + 1
136
 
      end do
137
 
c
138
 
      if (.not. ga_create_irreg(mt_dbl, nva, nbfpair,
139
 
     $     'mp2: backt', 1, 1, int_mb(k_map), nblock, g_buf))
140
 
     $     call errquit('mp2: backt: ga_create failed', nva*nbfpair,
141
 
     &       GA_ERR)
142
 
c     
143
 
      call util_file_name('2pdm', .true., .true., fname)
144
 
#ifdef EAFHACK
145
 
         if (.not. rtdb_get(rtdb, 'mp2:eaf_size_in_mb',
146
 
     $        MT_INT, 1, eaf_size_in_mb)) eaf_size_in_mb = 0
147
 
             if(eaf_size_in_mb.ne.0) then
148
 
                eaftype=1000000+eaf_size_in_mb*1024*1024
149
 
             else
150
 
                eaftype=1000000+inntsize*p_file_size*11/10
151
 
             endif
152
 
             if(ga_nodeid().eq.0) then
153
 
                write(6,*)' EAF nonsep ',
154
 
     A            eaf_size_in_mb, eaftype
155
 
                call util_flush(6)
156
 
             endif
157
 
#else
158
 
             eaftype=eaf_rw
159
 
#endif
160
 
      if (eaf_open(fname,eaftype,twopdmunit) .ne. 0)
161
 
     $     call errquit('mp2: backt: failed to open file',0, DISK_ERR)
162
 
c     
163
 
c     Now do the back transformation
164
 
c
165
 
      call ga_sync()
166
 
      call pstat_on(ps_backt)
167
 
      call mp2_back_transform_uhf(
168
 
     $     nbf, 
169
 
     $     nir, 
170
 
     $     nva, nva_lo, nva_hi, num_va, num_oa, 
171
 
     $     nvb, nvb_lo, nvb_hi, num_vb, num_ob,
172
 
     $     sym_lo_oa, sym_hi_oa, sym_lo_va, sym_hi_va,
173
 
     $     sym_lo_ob, sym_hi_ob, sym_lo_vb, sym_hi_vb,
174
 
     $     oseg_lo, oseg_hi, 
175
 
     $     irs_a, irs_b,
176
 
     $     c_a, c_b, 
177
 
     $     nva_lo_local, nva_hi_local, 
178
 
     $     tunita, tunitb, twopdmunit,
179
 
     $     nshpair, nbfpair, int_mb(k_shpairs),
180
 
     $     nshpairlocal, nbfpairlocal, int_mb(k_shpairslocal),
181
 
     $     int_mb(k_shdim), int_mb(k_shlo),
182
 
     $     dbl_mb(k_t), dbl_mb(k_iauv), dbl_mb(k_tmp),
183
 
     $     int_mb(k_map), g_buf)
184
 
      call pstat_off(ps_backt)
185
 
      call ga_sync()
186
 
c
187
 
c     Free up some memory and then contract with gradient integrals.
188
 
c
189
 
      if (.not. ga_destroy(g_buf)) call errquit
190
 
     $     ('mp2: backt: ga destroy failed', 0, GA_ERR)
191
 
      if (.not. ma_chop_stack(l_t)) call errquit
192
 
     $     ('mp2:backt: first chop stack failed', 0, MA_ERR)
193
 
c
194
 
      tdim = shmax*shmax*nbf*ninseg
195
 
      if (.not. ma_push_get(mt_dbl, tdim, 'mp2:back tt',
196
 
     $     l_t, k_t)) call errquit
197
 
     $     ('mp2: backt: failed ma for nonsep test', tdim, MA_ERR)
198
 
      tdim = shmax*shmax*nbf
199
 
      if (.not. ma_push_get(mt_dbl, tdim*ninseg, 'mp2:back tbuf',
200
 
     $     l_tmp, k_tmp)) call errquit
201
 
     $     ('mp2: backt: failed ma for nonsep test', tdim, MA_ERR)
202
 
      if (.not. ma_push_get(mt_dbl, ninseg*nbf, 'mp2:back c_t',
203
 
     $     l_c_t, k_c_t)) call errquit
204
 
     $     ('mp2: backt: failed ma for nonsep test', ninseg*nbf, MA_ERR)
205
 
      do i = oseg_lo, oseg_hi
206
 
         call dcopy(nbf, c_a(1,i), 1, dbl_mb(k_c_t+i-oseg_lo), ninseg)
207
 
      end do
208
 
c
209
 
c     Determine list of active centers
210
 
c
211
 
      if (.not. geom_ncent(geom, natoms)) call errquit
212
 
     $     ('mp2_backt: geom ?',0, GEOM_ERR)
213
 
      if (.not. ma_push_get(mt_log, natoms, 'mp2:back act',
214
 
     $     l_act, k_act)) call errquit('mp2:back ma ', natoms, MA_ERR)
215
 
      call grad_active_atoms(rtdb, natoms, log_mb(k_act), nactive)
216
 
c
217
 
c     Turn this into a list of active shells
218
 
c
219
 
      if (.not. ma_push_get(mt_log, nsh, 'mp2:back actsh',
220
 
     $     l_actsh, k_actsh)) call errquit('mp2:back ma ', nsh, MA_ERR)
221
 
      do i = 1, natoms
222
 
         if (.not. bas_ce2cnr(basis, i, ishlo, ishhi))
223
 
     $        call errquit('mp2:backt basis?',0, BASIS_ERR)
224
 
         do ish = ishlo,ishhi
225
 
            log_mb(k_actsh+ish-1) = log_mb(k_act+i-1)
226
 
         end do
227
 
      end do
228
 
c
229
 
      call intd_init(rtdb, 1, basis)
230
 
      call intb_mem_2e4c(leneri, lenscr) ! blocking algorithm
231
 
      leneri = leneri/12
232
 
      leneri = max(leneri,15**4) ! 1 G quartet = 39 D quartets ... SEE MP2_MEMORY
233
 
      if (.not. ma_push_get(mt_dbl,12*leneri,'deriv buffer',
234
 
     $     l_eri,k_eri)) call errquit
235
 
     $     ('mp2:backt could not allocate buffer',12*leneri, MA_ERR)
236
 
      if (.not. ma_push_get(mt_dbl,shmax**4,'pdm buffer',
237
 
     $     l_pdm,k_pdm)) call errquit
238
 
     $     ('mp2:backt could not allocate pdm',shmax*4, MA_ERR)
239
 
      if (.not. ma_push_get(mt_int,4*leneri,'deriv labels',l_lab,k_lab))
240
 
     $     call errquit('mp2:backt could not allocate labels',leneri,
241
 
     &       MA_ERR)
242
 
      if (.not. ma_push_get(mt_dbl,lenscr,'deriv scratch',
243
 
     $     l_scr,k_scr))
244
 
     $     call errquit('mp2:backt could not allocate scr',lenscr,
245
 
     &       MA_ERR)
246
 
      if (.not. ma_push_get(mt_int,nbf,'shmap',l_bftosh,k_bftosh))
247
 
     $     call errquit('mp2:backt could not allocate bftosh',nbf,
248
 
     &       MA_ERR)
249
 
      if (.not. ma_push_get(mt_int,nbf,'cemap',l_bftoce,k_bftoce))
250
 
     $     call errquit('mp2:backt could not allocate bftoce',nbf,
251
 
     &       MA_ERR)
252
 
c      
253
 
      call ga_sync()
254
 
      call pstat_on(ps_nonsep)
255
 
      call mp2_nonsep(rtdb, basis,
256
 
     $     twopdmunit,
257
 
     $     dbl_mb(k_c_t),
258
 
     $     nsh, nbf, ninseg,
259
 
     $     nshpairlocal, int_mb(k_shpairslocal),
260
 
     $     int_mb(k_shdim), int_mb(k_shlo),
261
 
     $     dbl_mb(k_t), dbl_mb(k_tmp), grad,
262
 
     $     tol2e, oskel,
263
 
     $     log_mb(k_actsh), leneri, dbl_mb(k_eri), int_mb(k_lab),
264
 
     $     lenscr, dbl_mb(k_scr), int_mb(k_bftosh), int_mb(k_bftoce),
265
 
     $     dbl_mb(k_pdm))
266
 
      call pstat_off(ps_nonsep)
267
 
      call ga_sync()
268
 
      call intd_terminate()
269
 
c
270
 
c     Zero out gradients on inactive atoms
271
 
c
272
 
      do i = 1, natoms
273
 
         if (.not. log_mb(k_act+i-1)) then
274
 
            do j = 1, 3
275
 
               grad(j,i) = 0.0d0
276
 
            end do
277
 
         end if
278
 
      end do
279
 
c
280
 
      call schwarz_tidy()
281
 
c
282
 
c     Done.
283
 
c     
284
 
      if (.not. ma_chop_stack(l_shpairs)) call errquit
285
 
     $     ('mp2: backt: failed chopping stack',0, MA_ERR)
286
 
c
287
 
      if (util_print('iostats', print_high) .and. 
288
 
     $     ga_nodeid().eq.0) call eaf_print_stats(twopdmunit)
289
 
      if (eaf_close(twopdmunit) .ne. 0)
290
 
     $     call errquit('mp2: backt: closing 2pdm',ierr, DISK_ERR)
291
 
      call util_file_unlink(fname)
292
 
c
293
 
      end
294
 
      subroutine mp2_back_transform_uhf(
295
 
     $     nbf, 
296
 
     $     nir, 
297
 
     $     nva, nva_lo, nva_hi, num_va, num_oa, 
298
 
     $     nvb, nvb_lo, nvb_hi, num_vb, num_ob,
299
 
     $     sym_lo_oa, sym_hi_oa, sym_lo_va, sym_hi_va,
300
 
     $     sym_lo_ob, sym_hi_ob, sym_lo_vb, sym_hi_vb,
301
 
     $     oseg_lo, oseg_hi, 
302
 
     $     irs_a, irs_b,
303
 
     $     c_a, c_b, 
304
 
     $     nva_lo_local, nva_hi_local, 
305
 
     $     tunita, tunitb, twopdmunit,
306
 
     $     nshpair, nbfpair, shpairs,
307
 
     $     nshpairlocal, nbfpairlocal, shpairslocal,
308
 
     $     shdim, shlo,
309
 
     $     t, ia_uv, tmp,
310
 
     $     map, g_buf)
311
 
      implicit none
312
 
#include "errquit.fh"
313
 
#include "global.fh"
314
 
#include "mafdecls.fh"
315
 
#include "eaf.fh"
316
 
#include "util.fh"
317
 
      integer nbf               ! No. of basis functions
318
 
      integer nir               ! No. of irreducible representations
319
 
      integer nva, nva_lo, nva_hi ! Number and ranges of virtual orbitals
320
 
      integer nvb, nvb_lo, nvb_hi
321
 
      integer num_va(0:nir-1), num_vb(0:nir-1) !  No. of vir of each symmetry
322
 
      integer num_oa(0:nir-1), num_ob(0:nir-1) !  No. of occ of each symmetry
323
 
      integer sym_lo_oa(0:nir-1), sym_hi_oa(0:nir-1),
324
 
     $     sym_lo_va(0:nir-1), sym_hi_va(0:nir-1) ! Ranges of each symmetry
325
 
      integer sym_lo_ob(0:nir-1), sym_hi_ob(0:nir-1),
326
 
     $     sym_lo_vb(0:nir-1), sym_hi_vb(0:nir-1)
327
 
      integer oseg_lo, oseg_hi  ! Range of occupied for this pass
328
 
      integer irs_a(nbf), irs_b(nbf) ! Orbital symmetries
329
 
      double precision c_a(nbf,*), c_b(nbf,*) ! MO coefficients
330
 
      integer nva_lo_local, nva_hi_local ! Range of virtuals on this node
331
 
      integer tunita, tunitb    ! Unit no.s for pure and mixed spin T
332
 
      integer nshpair           ! total no. of non-zero shell pairs
333
 
      integer nshpairlocal      ! no. of pairs assigned to this processor
334
 
      integer nbfpair           ! dimension of bf pairs
335
 
      integer nbfpairlocal      ! dimension of local bf pairs
336
 
      integer shpairs(2,nshpair) ! (u>=v) shell pair indices
337
 
      integer shpairslocal(2,*) ! local shell pair indices
338
 
      integer shdim(*)          ! no. of bf in shell
339
 
      integer shlo(*)           ! first bf in shell
340
 
      double precision t(*)     ! Scratch max(nbf*nbf,nbf*S*S)
341
 
      double precision tmp(*)   ! Scratch max(nbf*nbf,nbf*S*S)
342
 
      double precision ia_uv(nbf,nbf) ! Scratch
343
 
      integer twopdmunit        ! Unit no. for part transformed two-PDM
344
 
      integer map(0:*)          ! First index of g_buf on each node
345
 
      integer g_buf             ! Global array handle for transpose buffer
346
 
c     
347
 
*     double precision tsum
348
 
 
349
 
      double precision zero, one, scale
350
 
      double precision tunitptra, tunitptrb ! Pointers into files
351
 
      double precision twopdmunitptr
352
 
      integer count, tcount, ind,  ptr
353
 
      integer symi, symj, syma, symb, symia
354
 
      integer i, a, u, v, ush, vsh, udim, vdim, ishpair
355
 
c
356
 
#include "bitops.fh"
357
 
c     
358
 
      zero=0.0d0
359
 
      one=1.0d0
360
 
c     
361
 
      tunitptra=1
362
 
      tunitptrb=1
363
 
*     tsum = 0
364
 
c     
365
 
      do i=oseg_lo,oseg_hi 
366
 
         call ga_fill_patch(g_buf,1,nva,1,nbfpair,99.0d0) ! For debug
367
 
         symi=irs_a(i)
368
 
         do a=nva_lo_local,nva_hi_local
369
 
            syma=irs_a(a)
370
 
            symia=ieor(syma,symi)
371
 
            call dfill((nbf*nbf),zero,ia_uv,1) ! Will add pure and mixed in
372
 
c     
373
 
            call mp2_read_tijab(nva_lo, nva_hi, irs_a, symia,
374
 
     $           num_oa, sym_hi_oa, sym_lo_oa, tunita, tunitptra, t)
375
 
c     
376
 
            tcount=1
377
 
            do symb=0,nir-1
378
 
               symj=ieor(symia,symb)
379
 
               if(num_va(symb).gt.0.and.num_oa(symj).gt.0) then
380
 
                  call dgemm('n','t', !             t(j,b)Cbv -> t(j,v)
381
 
     $                 num_oa(symj), nbf, num_va(symb),
382
 
     $                 one, t(tcount), num_oa(symj),
383
 
     $                 c_a(1,sym_lo_va(symb)), nbf,
384
 
     $                 zero, tmp, num_oa(symj))
385
 
                  call dgemm('n','n', !             Cuj t(j,v) -> t(u,v)
386
 
     $                 nbf, nbf, num_oa(symj),
387
 
     $                 one, c_a(1,sym_lo_oa(symj)), nbf,
388
 
     $                 tmp, num_oa(symj),
389
 
     $                 one, ia_uv, nbf)
390
 
                  tcount=tcount+num_oa(symj)*num_va(symb)
391
 
               end if
392
 
            end do
393
 
c     
394
 
            call mp2_read_tijab(nvb_lo, nvb_hi, irs_b, symia,
395
 
     $           num_ob, sym_hi_ob, sym_lo_ob, tunitb, tunitptrb, t)
396
 
c     
397
 
            tcount=1
398
 
            do symb=0,nir-1
399
 
               symj=ieor(symia,symb)
400
 
               if(num_vb(symb).gt.0.and.num_ob(symj).gt.0) then
401
 
                  call dgemm('n','t', !             t(j,b)Cbv -> t(j,v)
402
 
     $                 num_ob(symj), nbf, num_vb(symb),
403
 
     $                 one, t(tcount), num_ob(symj),
404
 
     $                 c_b(1,sym_lo_vb(symb)), nbf,
405
 
     $                 zero, tmp, num_ob(symj))
406
 
                  call dgemm('n','n', !             Cuj t(j,v) -> t(u,v)
407
 
     $                 nbf, nbf, num_ob(symj),
408
 
     $                 one, c_b(1,sym_lo_ob(symj)), nbf,
409
 
     $                 tmp, num_ob(symj),
410
 
     $                 one, ia_uv, nbf)
411
 
                  tcount=tcount+num_ob(symj)*num_vb(symb)
412
 
               end if
413
 
            end do
414
 
c     
415
 
c     Pack down into sparse, symmetry unique list and symmetrize over uv
416
 
c     Note we put in the whole square of diagonal shell blocks and divide
417
 
c     it by two.
418
 
c     
419
 
            ind = 0
420
 
            do ishpair = 1, nshpair
421
 
               ush = shpairs(1,ishpair)
422
 
               vsh = shpairs(2,ishpair)
423
 
               scale = 1.0d0
424
 
               if (ush.eq.vsh) scale = 0.5d0
425
 
               do u = shlo(ush), shlo(ush)+shdim(ush)-1
426
 
                  do v = shlo(vsh), shlo(vsh)+shdim(vsh)-1
427
 
                     ind = ind + 1
428
 
                     tmp(ind) = (ia_uv(u,v)+ia_uv(v,u))*scale
429
 
                  end do
430
 
               end do
431
 
            end do
432
 
            if (ind .ne. nbfpair) call errquit('mp2bt: ind?', ind,
433
 
     &       UNKNOWN_ERR)
434
 
c     
435
 
c     shpairs(1:2,1:nshpair)    (u>=v) shell pair indices
436
 
c     shdim(1:nsh)              dimension of this sh (s=1,p=2,...)
437
 
c     shlo(1:nsh)               first basis function in shell
438
 
c     nbfpair                   total no. of non-zero bf pairs
439
 
c     nshpair                   total no. of non-zero shell pairs
440
 
c     nshpairlocal              no. of shellpairs assigned to this processor
441
 
c     nbfpairlocal              no. of bfpairs assigned to this processor
442
 
c     shpairslocal(1:2,1:nshpairlocal) (u>=v) local shell pair indices
443
 
c     
444
 
            call ga_put(g_buf, a-nva_lo+1, a-nva_lo+1,
445
 
     $           1, nbfpair, tmp, 1)
446
 
c     
447
 
         end do                 ! End of a
448
 
         call ga_sync
449
 
c     
450
 
c     Now have locally (a,1:nbfpairlocal).  Loop thru local shell pairs
451
 
c     transform the a index and write to disk.
452
 
c     
453
 
         ind = map(ga_nodeid())
454
 
         ptr = 0
455
 
         do ishpair = 1, nshpairlocal
456
 
            ush = shpairslocal(1,ishpair)
457
 
            vsh = shpairslocal(2,ishpair)
458
 
            udim = shdim(ush)
459
 
            vdim = shdim(vsh)
460
 
c     
461
 
            call ga_get(g_buf,1,nva,ind,ind+udim*vdim-1,tmp,nva)
462
 
            ind = ind + udim*vdim
463
 
            call dgemm('n','n', ! C(s,a)*F(a,uv) -> F(s,uv)
464
 
     $           nbf, udim*vdim, nva,
465
 
     $           one, c_a(1,nva_lo), nbf,
466
 
     $           tmp, nva,
467
 
     $           zero, t, nbf)
468
 
c     
469
 
c     On disk is D(1:nbf,1:vdim(v),1:udim(u),oseg_lo:oseg_hi,1:nshpairlocal)
470
 
c     with indices D(s,vf,uf,i,ishpair)
471
 
c     
472
 
c     Note - this presumes that the routine reading the data back in will
473
 
c     be using an algorithm with o*n*s*s memory requirement (n=#bf, s=#bf 
474
 
c     in shell, o=occupied).
475
 
c     
476
 
            count = nbf*udim*vdim ! Amount of data to write
477
 
            twopdmunitptr = 8.0d0*(ptr + (i-oseg_lo)*count)
478
 
c     
479
 
            if (eaf_write(twopdmunit, twopdmunitptr, t, count*8).ne.0)
480
 
     $           call errquit('mp2_bt: write of two pdm?', 0,
481
 
     &       DISK_ERR)
482
 
c     
483
 
            ptr = ptr + count*(oseg_hi-oseg_lo+1)
484
 
         end do
485
 
c     
486
 
         call ga_sync
487
 
c     
488
 
      end do                    ! End of i
489
 
c     
490
 
      end
491
 
      subroutine mp2_backt_info(
492
 
     $     basis,
493
 
     $     tol2e, oskel,
494
 
     $     omakearrays,
495
 
     $     nshpair, nshpairlocal, 
496
 
     $     nbfpair, nbfpairlocal, shmax,
497
 
     $     shdim, shlo, shpairs, shpairslocal, map)
498
 
      implicit none
499
 
#include "errquit.fh"
500
 
#include "bas.fh"
501
 
#include "schwarz.fh"
502
 
#include "global.fh"
503
 
#include "util.fh"
504
 
#include "sym.fh"
505
 
      integer basis             ! [input] basis set handle
506
 
      double precision tol2e    ! [input] screening threshold
507
 
      logical oskel             ! [input] if true use skeleton symm
508
 
      logical omakearrays       ! [input] if true then make the arrays
509
 
      integer nshpair           ! [output] no. of non-zero shell pairs
510
 
      integer nshpairlocal      ! [output] no. of local shell pairs
511
 
      integer shmax             ! [output] Max AO shell dimension
512
 
      integer nbfpair           ! [output] sum of all non-zero pair dims
513
 
      integer nbfpairlocal      ! [output] sum of local non-zero pair dims
514
 
      integer shdim(*)          ! [output] dimension of each shell
515
 
      integer shlo(*)           ! [output] first bf in each shell
516
 
      integer shpairs(2,*)      ! [output] 1->u, 2->v u>=v shells in pair
517
 
      integer shpairslocal(2,*) ! [output] ditto but only for local pairs
518
 
      integer map(0:*)          ! [output] Map for g_create_irreg for g_buf
519
 
c     
520
 
c     Form a list of interacting shell pairs (u>=v) and assign them
521
 
c     to processors.  Ideally do this so as to optimize the efficiency
522
 
c     of derivative integral evaluation and to provide good load balance 
523
 
c     for the back transformation.  Right now we just do round-robin.
524
 
c     
525
 
c     If (omakearrays) then 
526
 
c     actually make the arrays
527
 
c     else
528
 
c     just return the scalar results
529
 
c     end if
530
 
c     
531
 
c     The Schwarz package must be initialized before entry
532
 
c     
533
 
      integer u, v, nbf, nsh, lo, hi, udim, vdim, me, nproc, owner
534
 
      integer count, tmp, ishpair, first, last, iproc, uvdim
535
 
      logical odoit
536
 
      double precision q2, sss
537
 
c     
538
 
      integer npairblock        ! Crude attempt to load balance the pairs
539
 
      parameter (npairblock = 7)
540
 
      integer pairblock(npairblock), ipairb, npairmin, npairmax
541
 
      integer logmin, logmax
542
 
      intrinsic log10
543
 
      data pairblock /1,3,6,9,18,36,1000000/ ! Selected from s*s,s*p,p*p ...
544
 
c
545
 
      if (.not. bas_numbf(basis, nbf))
546
 
     $     call errquit('mp2_backt_info: bad basis handle', 0,
547
 
     &       BASIS_ERR)
548
 
      if (.not. bas_numcont(basis, nsh))
549
 
     $     call errquit('mp2_backt_info: bad basis handle', 0,
550
 
     &       BASIS_ERR)
551
 
c     
552
 
      me = ga_nodeid()
553
 
      nproc = ga_nnodes()
554
 
      shmax = 0
555
 
c     
556
 
      if (omakearrays) then
557
 
c     
558
 
c     First compute offset to first shell pair on each processor
559
 
c     
560
 
         nshpair = 0
561
 
         call ifill(nproc, 0, map, 1)
562
 
         logmax = 1000000
563
 
         do logmin = -1,-13,-4  ! NOTE THAT SCREENING IS IMPLIED IN THIS LOOP
564
 
            npairmin = 0
565
 
            do ipairb = 1, npairblock ! Loop thru types of pairs
566
 
               npairmax = pairblock(ipairb)
567
 
               do u = nsh, 1, -1
568
 
                  if (.not. bas_cn2bfr(basis, u, lo, hi))
569
 
     $                 call errquit('mp2_backt_info: bas range', u,
570
 
     &       BASIS_ERR)
571
 
                  udim = hi - lo + 1
572
 
                  do v = 1, u
573
 
                     if (.not. bas_cn2bfr(basis, v, lo, hi))
574
 
     $                    call errquit('mp2_backt_info:bas range',v,
575
 
     &       BASIS_ERR)
576
 
                     vdim = hi - lo + 1
577
 
                     uvdim = udim*vdim
578
 
                     sss = schwarz_shell(u,v)*schwarz_max() + 1d-100
579
 
                     if (uvdim.gt.npairmin .and. uvdim.le.npairmax .and.
580
 
     $                    log10(sss).lt.logmax .and. 
581
 
     $                    log10(sss).ge.logmin) then
582
 
c     
583
 
                        odoit = sss .gt. tol2e*0.1d0
584
 
                        if (odoit .and. oskel) 
585
 
     $                       odoit = sym_shell_pair(basis,u,v,q2)
586
 
                        if (odoit) then
587
 
                           owner = mod(nshpair,nproc)
588
 
                           map(owner) = map(owner) + 1
589
 
                           nshpair = nshpair + 1
590
 
                        end if
591
 
c     
592
 
                     end if
593
 
                  end do
594
 
               end do
595
 
               npairmin = npairmax
596
 
            end do
597
 
            logmax = logmin
598
 
         end do
599
 
         count = 1
600
 
         do u = 0, nproc-1
601
 
            tmp = map(u)
602
 
            map(u) = count
603
 
            count = count + tmp
604
 
         end do
605
 
      end if
606
 
c     
607
 
c     Make list of all shell pairs ordered by their assigned processor,
608
 
c     the list of local shell pairs, and also count the other information.
609
 
c     
610
 
      nshpairlocal = 0
611
 
      nbfpair = 0
612
 
      nbfpairlocal = 0
613
 
      nshpair = 0
614
 
      logmax = 1000000
615
 
      do logmin = -1,-13,-4  ! NOTE THAT SCREENING IS IMPLIED IN THIS LOOP
616
 
         npairmin = 0
617
 
         do ipairb = 1, npairblock ! Loop thru types of pairs
618
 
            npairmax = pairblock(ipairb)
619
 
            do u = nsh, 1, -1
620
 
               if (.not. bas_cn2bfr(basis, u, lo, hi))
621
 
     $              call errquit('mp2_backt_info: bas range', u,
622
 
     &       BASIS_ERR)
623
 
               shmax = max(shmax,hi-lo+1)
624
 
               if (omakearrays) then
625
 
                  shdim(u) = hi - lo + 1
626
 
                  shlo(u) = lo
627
 
               end if
628
 
               udim = hi - lo + 1
629
 
               do v = 1, u
630
 
                  if (.not. bas_cn2bfr(basis, v, lo, hi))
631
 
     $                 call errquit('mp2_backt_info: bas range', v,
632
 
     &       BASIS_ERR)
633
 
                  vdim = hi - lo + 1
634
 
                  uvdim = udim*vdim
635
 
                  sss = schwarz_shell(u,v)*schwarz_max() + 1d-100
636
 
                  if (uvdim.gt.npairmin .and. uvdim.le.npairmax .and.
637
 
     $                 log10(sss).lt.logmax .and. 
638
 
     $                 log10(sss).ge.logmin) then
639
 
c     
640
 
                     odoit = sss.gt.tol2e*0.1d0
641
 
                     if (odoit .and. oskel) 
642
 
     $                    odoit = sym_shell_pair(basis,u,v,q2)
643
 
                     if (odoit) then
644
 
                        owner = mod(nshpair,nproc)
645
 
                        if (owner.eq.me) then
646
 
                           nshpairlocal = nshpairlocal + 1
647
 
                           nbfpairlocal = nbfpairlocal + udim*vdim
648
 
                           if (omakearrays) then
649
 
                              shpairslocal(1,nshpairlocal) = u
650
 
                              shpairslocal(2,nshpairlocal) = v
651
 
                           end if
652
 
                        end if
653
 
                        nshpair = nshpair + 1
654
 
                        nbfpair = nbfpair + udim*vdim
655
 
                        if (omakearrays) then
656
 
                           shpairs(1,map(owner)) = u
657
 
                           shpairs(2,map(owner)) = v
658
 
                           map(owner) = map(owner) + 1
659
 
                        end if
660
 
                     end if
661
 
                  end if
662
 
               end do
663
 
            end do
664
 
            npairmin = npairmax
665
 
         end do
666
 
         logmax = logmin
667
 
      end do
668
 
c     
669
 
c     Now replace map with the map creating the global array for transposing
670
 
c     
671
 
      if (omakearrays) then
672
 
         count = 1
673
 
         first = 1
674
 
         do iproc = 0, nproc-1
675
 
            last = map(iproc) - 1
676
 
            map(iproc) = count
677
 
            do ishpair = first, last
678
 
               u = shpairs(1,ishpair)
679
 
               v = shpairs(2,ishpair)
680
 
               count = count + shdim(u)*shdim(v)
681
 
            end do
682
 
            first = last + 1
683
 
         end do
684
 
      end if
685
 
c     
686
 
      if (util_print('mp2_backt', print_debug)) then
687
 
         write(6,*) me, ' nshpair ', nshpair
688
 
         write(6,*) me, ' nshpairlocal ', nshpairlocal
689
 
         write(6,*) me, ' nbfpair ', nbfpair
690
 
         write(6,*) me, ' nbfpairlocal ', nbfpairlocal
691
 
         if (omakearrays) then
692
 
            do u = 1, nsh
693
 
               write(6,*) me, ' u shlo shdim ', u, shlo(u), shdim(u)
694
 
            end do
695
 
            do u = 1, nshpairlocal
696
 
               write(6,*) me, ' pair shpairslocal ', u,
697
 
     $              shpairslocal(1,u), shpairslocal(2,u)
698
 
            end do
699
 
            do u = 0, nproc-1
700
 
               write(6,*) me, ' proc map ', u, map(u)
701
 
            end do
702
 
         end if
703
 
      end if
704
 
c     
705
 
      end
706
 
      subroutine mp2_nonsep(
707
 
     $     rtdb,
708
 
     $     basis,
709
 
     $     twopdmunit,
710
 
     $     c_t,
711
 
     $     nsh, nbf, ninseg,
712
 
     $     nshpairlocal, shpairslocal, shdim, shlo,
713
 
     $     t, tbuf, grad,
714
 
     $     tol2e, oskel, oactive, leneri, eri, labels, lenscr, 
715
 
     $     scratch, shmap, cemap, pdm)
716
 
      implicit none
717
 
#include "errquit.fh"
718
 
#include "schwarz.fh"
719
 
#include "mafdecls.fh"
720
 
#include "util.fh"
721
 
#include "bas.fh"
722
 
#include "geom.fh"
723
 
#include "global.fh"
724
 
#include "eaf.fh"
725
 
#include "rtdb.fh"
726
 
#include "sym.fh"
727
 
#include "stdio.fh"
728
 
c     
729
 
      integer rtdb              ! [input]
730
 
      integer basis             ! [input]
731
 
      integer twopdmunit        ! [input]
732
 
      integer nsh, nbf, ninseg  ! [input] 
733
 
      double precision c_t(ninseg,nbf) ! [input] Transposed MOs in segment
734
 
      integer nshpairlocal, shpairslocal(2,*), shdim(nsh), shlo(nsh) ! [input]
735
 
      double precision t(*)     ! [scratch] ninseg*nbf*S*S
736
 
      double precision tbuf(*)  ! [scratch] nbf*S*S
737
 
      double precision grad(3,*) ! [input/output]
738
 
      double precision tol2e    ! [input]
739
 
      logical oskel             ! [input]
740
 
      logical oactive(*)        ! [input] oactive(ish)=true if shell is active
741
 
      integer leneri, lenscr    ! [input]
742
 
      integer labels(leneri,4), nq, nint ! [scratch]
743
 
      double precision eri(3,4,leneri), scratch(lenscr) ! [scratch]
744
 
      integer shmap(*), cemap(*) ! [scratch]
745
 
      double precision pdm(*)   ! [scratch]
746
 
c     
747
 
      double precision fileptr, energy, scale, psum, q4, block_eff
748
 
      integer i, count, ishpair
749
 
      integer ush, vsh, xsh, ysh, u, v, x, y
750
 
      integer udim, vdim, xdim, ydim, xlo, ylo
751
 
      integer ush_cur, vsh_cur, xsh_cur, ysh_cur
752
 
      integer ush_prev, vsh_prev, xsh_prev, ysh_prev
753
 
      integer xcent, ycent, ucent, vcent, ijkl, xsh_start, ysh_start
754
 
      logical omore
755
 
c
756
 
      integer maxq
757
 
      parameter (maxq=10000)
758
 
      integer sh_list(maxq,4)
759
 
      double precision q4_list(maxq)
760
 
c
761
 
      logical intbd_init4c, intbd_2e4c
762
 
      external intbd_init4c, intbd_2e4c
763
 
c     
764
 
      integer nat, geom,
765
 
     $     ulo_cur, uhi_cur, vlo_cur, vhi_cur, xlo_cur, xhi_cur, 
766
 
     $     ylo_cur, yhi_cur, udim_cur, vdim_cur, xdim_cur, ydim_cur
767
 
      logical status, odebug, odoit, sym_shell_quartet, oenergy
768
 
      external sym_shell_quartet
769
 
c     
770
 
      odebug = util_print('mp2_backt', print_debug)
771
 
      oenergy = util_print('backtenergy', print_debug)
772
 
c     
773
 
      energy = 0.0d0
774
 
c     
775
 
      status = bas_geom ( basis, geom )
776
 
      status = geom_ncent ( geom, nat )
777
 
c     
778
 
      do u = 1, nbf
779
 
         if (.not. bas_bf2ce(basis,u,ucent))
780
 
     $        call errquit('mp2g: bad something?',0, BASIS_ERR)
781
 
         if (.not. bas_bf2cn(basis,u,ush))
782
 
     $        call errquit('mp2g: bad something?',0, BASIS_ERR)
783
 
         cemap(u) = ucent
784
 
         shmap(u) = ush
785
 
      enddo
786
 
c     
787
 
      if (odebug) then
788
 
         write(6,*) ' Transposed occupied MOS in segment '
789
 
         call output(c_t, 1, ninseg, 1, nbf, ninseg, nbf, 1)
790
 
      end if
791
 
c     
792
 
      fileptr = 0.0d0
793
 
      do ishpair = 1, nshpairlocal
794
 
         ush = shpairslocal(1,ishpair)
795
 
         vsh = shpairslocal(2,ishpair)
796
 
         udim = shdim(ush)
797
 
         vdim = shdim(vsh)
798
 
         count = nbf*udim*vdim
799
 
c     
800
 
c     Part transformed density is stored as 
801
 
c     t(1:nbf,1:vdim,1:udim,1:ninseg,1:nshpairlocal) -> 
802
 
c     t(x,v,u,i,ush,vsh)
803
 
c     
804
 
c     Read this in and transpose to in core structure t(i,x,v,u)
805
 
c     
806
 
        if (eaf_read(twopdmunit, fileptr, tbuf, ninseg*count*8).ne.0)
807
 
     $       call errquit('mp2: ao test: failed reading density', 0,
808
 
     &     DISK_ERR)
809
 
         do i = 1, ninseg
810
 
            call dcopy(count, tbuf(count*(i-1)+1), 1, t(i), ninseg)
811
 
         end do
812
 
            fileptr = fileptr + ninseg*count*8.0d0
813
 
c     
814
 
         xsh_start = 1          ! For braindead multipassing
815
 
         ysh_start = 1
816
 
c
817
 
 333     nq = 0
818
 
         do xsh = xsh_start, nsh
819
 
            do ysh = ysh_start, xsh
820
 
               odoit = schwarz_shell(xsh,ysh)*schwarz_shell(ush,vsh)
821
 
     $              .gt. tol2e*0.1d0
822
 
               odoit = odoit .and. 
823
 
     $              (oactive(ush).or.oactive(vsh).or.
824
 
     $              oactive(xsh).or.oactive(ysh))
825
 
               q4 = 1.0d0
826
 
               if (odoit .and. oskel) odoit = 
827
 
     $              sym_shell_quartet(basis, ush, vsh, xsh, ysh, q4)
828
 
c     
829
 
               odoit = .true.
830
 
               if (odoit) then
831
 
c     
832
 
                  scale = q4
833
 
                  if (xsh .eq. ysh) scale = scale*0.5d0
834
 
                  xdim = shdim(xsh)
835
 
                  ydim = shdim(ysh)
836
 
                  xlo  = shlo(xsh)
837
 
                  ylo  = shlo(ysh)
838
 
 
839
 
                  call mp2_make_two_particle_density(psum,
840
 
     $                 udim, vdim, xdim, ydim, xlo, ylo, 
841
 
     $                 ninseg, nbf, 
842
 
     $                 scale, c_t, t, pdm)
843
 
c     
844
 
c     use density to screen evaluation of gradient integrals
845
 
                     if (oenergy) then
846
 
                        call int_2e4c(basis, ush, vsh, basis, xsh, ysh, 
847
 
     $                       lenscr, scratch, leneri, eri)
848
 
                        energy = energy + 
849
 
     $                       ddot(udim*vdim*xdim*ydim,eri,1,pdm,1)
850
 
                     endif
851
 
c
852
 
                  if (schwarz_shell(xsh,ysh)*schwarz_shell(ush,vsh)*psum
853
 
     $                 .gt. tol2e) then
854
 
c
855
 
                     nq = nq + 1
856
 
                     sh_list(nq,1) = ysh
857
 
                     sh_list(nq,2) = xsh
858
 
                     sh_list(nq,3) = vsh
859
 
                     sh_list(nq,4) = ush
860
 
                     q4_list(nq)   = scale
861
 
                     if (nq .eq. maxq) goto 666 ! UGLY 
862
 
                  endif         ! schwarz and density
863
 
               endif            ! schwarz
864
 
            enddo               ! ysh
865
 
            ysh_start = 1
866
 
         enddo                  ! xsh
867
 
c
868
 
 666     xsh_start = xsh
869
 
         ysh_start = ysh+1
870
 
c     
871
 
         omore = .false.
872
 
c     
873
 
         if (.not. intbd_init4c(
874
 
     $        basis, sh_list(1,1), sh_list(1,2),
875
 
     $        basis, sh_list(1,3), sh_list(1,4),
876
 
     $        nq, q4_list, .true., lenscr, scratch, leneri, 
877
 
     $        block_eff)) call errquit('mp2:backt: txs init?',nq,
878
 
     &       INT_ERR)
879
 
c     
880
 
 1000    omore = intbd_2e4c(
881
 
     $        basis, sh_list(1,1), sh_list(1,2),
882
 
     $        basis, sh_list(1,3), sh_list(1,4),
883
 
     $        nq, q4_list, .true., tol2e, .false.,
884
 
     $        labels(1,1),labels(1,2), 
885
 
     $        labels(1,3), labels(1,4), 
886
 
     $        eri, leneri, nint, lenscr, scratch)
887
 
 
888
 
         if (nint .gt. leneri) call errquit('mp2_nonsep: nint', nint,
889
 
     &       INT_ERR)
890
 
c     
891
 
         ush_prev = -1
892
 
         vsh_prev = -1
893
 
         ysh_prev = -1
894
 
         xsh_prev = -1
895
 
         do ijkl = 1, nint
896
 
            u = labels(ijkl,4)
897
 
            v = labels(ijkl,3)
898
 
            x = labels(ijkl,2)
899
 
            y = labels(ijkl,1)
900
 
            ush_cur = shmap(u)
901
 
            vsh_cur = shmap(v)
902
 
            xsh_cur = shmap(x)
903
 
            ysh_cur = shmap(y)
904
 
c
905
 
            if ( ush_cur.ne.ush_prev .or.
906
 
     $           vsh_cur.ne.vsh_prev .or.
907
 
     $           xsh_cur.ne.xsh_prev .or.
908
 
     $           ysh_cur.ne.ysh_prev ) then
909
 
c     
910
 
               ush_prev = ush_cur
911
 
               vsh_prev = vsh_cur
912
 
               xsh_prev = xsh_cur
913
 
               ysh_prev = ysh_cur
914
 
c     
915
 
               ulo_cur  = shlo(ush_cur)
916
 
               vlo_cur  = shlo(vsh_cur)
917
 
               xlo_cur  = shlo(xsh_cur)
918
 
               ylo_cur  = shlo(ysh_cur)
919
 
c     
920
 
               udim_cur = shdim(ush_cur)
921
 
               vdim_cur = shdim(vsh_cur)
922
 
               xdim_cur = shdim(xsh_cur)
923
 
               ydim_cur = shdim(ysh_cur)
924
 
c     
925
 
               uhi_cur  = ulo_cur + udim_cur - 1
926
 
               vhi_cur  = vlo_cur + vdim_cur - 1
927
 
               xhi_cur  = xlo_cur + xdim_cur - 1
928
 
               yhi_cur  = ylo_cur + ydim_cur - 1
929
 
c     
930
 
               ucent = cemap(u)
931
 
               vcent = cemap(v)
932
 
               xcent = cemap(x)
933
 
               ycent = cemap(y)
934
 
c     
935
 
               call mp2_make_two_particle_density( psum,
936
 
     $              udim_cur, vdim_cur, xdim_cur, ydim_cur, 
937
 
     $              xlo_cur, ylo_cur, 
938
 
     $              ninseg, nbf, 
939
 
     $              1.0d0, c_t, t, pdm)
940
 
            endif
941
 
c     
942
 
            call make_mp2grad(ucent,vcent,xcent,ycent, u,v,x,y,
943
 
     $           ulo_cur, uhi_cur, vlo_cur, vhi_cur,
944
 
     $           xlo_cur, xhi_cur, ylo_cur, yhi_cur,
945
 
     $           pdm, eri(1,1,ijkl), grad )
946
 
c     
947
 
         enddo
948
 
c     
949
 
         if (omore) goto 1000   ! Texas split the request
950
 
c
951
 
         if (xsh_start .le. nsh) goto 333
952
 
c     
953
 
      end do
954
 
c     
955
 
      if (oenergy) then
956
 
         call ga_dgop(1,energy,1,'+')
957
 
         if (ga_nodeid().eq.0) 
958
 
     +     write(LuOut,*) ' The energy from the two PDM is ', energy
959
 
      end if
960
 
c     
961
 
      end
962
 
      subroutine mp2_make_two_particle_density(
963
 
     $     psum, udim, vdim, xdim, ydim, xlo, ylo, ninseg, nbf, 
964
 
     $     scale, c_t, t, density)
965
 
      implicit none
966
 
c     
967
 
      double precision psum     ! [output] Norm of density
968
 
      integer udim, vdim, xdim, ydim ! [input] Shell dimensions
969
 
      integer xlo, ylo          ! [input] Start of shells
970
 
      integer ninseg            ! [input] No. of i in the batch
971
 
      integer nbf               ! [input]
972
 
      double precision scale    ! [input] Q4 and other factors
973
 
      double precision c_t(ninseg,nbf) ! [input] Transposed MO vectors in batch
974
 
      double precision t(ninseg,nbf,vdim,udim) ! [input] Partially transformed amplitudes
975
 
      double precision density(ydim,xdim,vdim,udim) ! [output] 
976
 
c
977
 
      integer u, v, x, y, i
978
 
      double precision p
979
 
c
980
 
      psum = 0.0d0
981
 
      do u = 1, udim
982
 
         do v = 1, vdim
983
 
            do x = 1, xdim
984
 
               do y = 1, ydim
985
 
                  p = 0.0d0
986
 
                  do i = 1, ninseg
987
 
                     p = p + t(i,x+xlo-1,v,u)*c_t(i,y+ylo-1)
988
 
     $                     + t(i,y+ylo-1,v,u)*c_t(i,x+xlo-1)
989
 
                  end do
990
 
                  p = p*scale
991
 
                  psum = psum + p*p
992
 
                  density(y,x,v,u) = p
993
 
               end do
994
 
            end do
995
 
         end do
996
 
      end do
997
 
c
998
 
      psum = sqrt(psum)
999
 
c
1000
 
      end
1001
 
      subroutine make_mp2grad(ucent,vcent,xcent,ycent, u,v,x,y,
1002
 
     $                        ulo, uhi, vlo, vhi, xlo, xhi, ylo, yhi,
1003
 
     $                        density, eri, grad )
1004
 
      implicit none
1005
 
      integer ucent,vcent,xcent,ycent, u,v,x,y 
1006
 
      integer ulo, uhi, vlo, vhi, xlo, xhi, ylo, yhi 
1007
 
      integer icart
1008
 
      double precision density(ylo:yhi,xlo:xhi,vlo:vhi,ulo:uhi)
1009
 
      double precision eri(3,4) , grad(3,*)
1010
 
c
1011
 
c ycent - eri(*,1)
1012
 
c xcent - eri(*,2)
1013
 
c vcent - eri(*,3)
1014
 
c ucent - eri(*,4)
1015
 
c
1016
 
      do icart=1,3
1017
 
         grad(icart,ucent)=grad(icart,ucent)+
1018
 
     $                     density(y,x,v,u)*eri(icart,4)
1019
 
         grad(icart,vcent)=grad(icart,vcent)+
1020
 
     $                     density(y,x,v,u)*eri(icart,3)
1021
 
         grad(icart,xcent)=grad(icart,xcent)+
1022
 
     $                     density(y,x,v,u)*eri(icart,2)
1023
 
         grad(icart,ycent)=grad(icart,ycent)+
1024
 
     $                     density(y,x,v,u)*eri(icart,1)
1025
 
      enddo
1026
 
c
1027
 
      end
1028
 
c--------------------------------------------------