1!
2!    Argonne Leadership Computing Facility
3!    BlueGene/P version
4!    Written by Vitali Morozov <morozov@anl.gov>
5!    Updated 20091022
6!    -O5 -qipa=noinline -qsmp=omp
7!
8!    modified from the original (use of einv)
9!
10      subroutine ccsd_tengy_bgp(f1n,f1t,f2n,f2t,f3n,f3t,f4n,f4t,
11     1                          dintc,dintx,t1v,einv,
12     2                          emp4,emp5,nvir)
13      implicit none
14      integer nvir,b,c
15      double precision emp4,emp5,denom
16      double precision two,three,four
17      double precision f1n(nvir,nvir),f1t(nvir,nvir)
18      double precision f2n(nvir,nvir),f2t(nvir,nvir)
19      double precision f3n(nvir,nvir),f3t(nvir,nvir)
20      double precision f4n(nvir,nvir),f4t(nvir,nvir)
21      double precision dintc(nvir),dintx(nvir),t1v(nvir)
22      double precision einv(nvir,nvir)
23      double precision e1, e2, e3, e4
24      double precision t0, t1, t2, t3, t4, t5, t6, t7, t8, t9
25      double precision z0, z1, z3, z5, z9, s0
26c
27!$omp parallel do
28!$omp& private(b,c,e1,e2,e3,e4,z0,z1,z3,z5,z9,s0,
29!$omp&         t0,t1,t2,t3,t4,t5,t6,t7,t8,t9)
30!$omp& shared(f1n,f2n,f3n,f4n,f1t,f2t,f3t,f4t,dintx,dintc,t1v,einv)
31!$omp& reduction(+:emp5,emp4) schedule(static)
32      do b = 1, nvir
33         do c = 1, nvir
34
35            s0 = einv(c,b)
36
37            e1 = f1n(c,b) + f4n(c,b) - 2d0*( f2n(c,b) + f3n(c,b) )
38            e2 = e1 + 3d0*f1n(c,b)
39
40            e3 = f1t(c,b) -2d0 * f2t(c,b)
41            e4 = e3 + f2n(c,b) - 2d0*( f3t(c,b) + f4n(c,b) + f1n(c,b)
42     1                                - 2d0*( f3n(c,b) + f4t(c,b) ) )
43
44            emp5 = emp5 + ( t1v(c)*( e1*dintc(b) + e4*dintx(b) )
45     1                    + t1v(b)*( e2*dintx(c) + e3*dintc(c) ) ) * s0
46
47            t1 = f1t(c,b) + f3n(c,b)
48            t2 = f1n(c,b) + f2t(c,b) + f4n(c,b)
49            t9 = f1t(c,b) + f4t(c,b)
50            t0 = f2t(c,b) + f3t(c,b)
51            t3 = t9 - 2d0*t0
52            t4 = t1 + f2n(c,b)
53            t5 = 2d0*t9 - t0
54            t6 = f3n(c,b) + 2d0*f4t(c,b)
55            t7 = f3n(c,b) * f4t(c,b)
56            t8 = t1*t3 - f1n(c,b)*t5 +3d0*( f1n(c,b)*f1n(c,b) + t7 )
57
58            z0 = f2t(b,c) + f3t(b,c)
59            z9 = f1t(b,c) + f4t(b,c)
60            z1 = f1t(b,c) + f3n(b,c)
61            z3 = z9 - 2d0*z0
62            z5 = 2d0*z9 - z0
63
64            emp4 = emp4 + ( t8 + t2*z3 - t4*z5
65     1           + 3d0*( t6*f1n(b,c) + f2t(c,b)*f2n(b,c) ) ) * s0
66         enddo
67      enddo
68c
69      return
70      end
71