1C> \ingroup wfn1_nxt
2C> @{
3C>
4C> \brief Calculate the occupation function penalty energy and
5C> gradient contributions
6C>
7C> Compute the occupation function penalty energy contribution as
8C> \f{eqnarray*}{
9C>   L_1 &=& \sum_\sigma\sum_{pq}\lambda^{O\sigma}_{pq}\left(
10C>           \sum_i O^\sigma_{ip}O^\sigma_{iq} - \delta_{pq}\right)^2
11C> \f}
12C> and the corresponding gradient contributions to the occupation
13C> function and penalty coefficient components. See wfn1_next_step.dox
14C> for details.
15C>
16      subroutine wfn1_nxt_Po_eg(nbf,oa,ob,poa,pob,
17     +           doa,dob,dpoa,dpob,EPo,ta,tb)
18      implicit none
19c
20      integer nbf !< [Input] The number of basis functions
21c
22      double precision oa(nbf,nbf)  !< [Input] The alpha occupation
23                                    !< functions
24      double precision ob(nbf,nbf)  !< [Input] The beta occupation
25                                    !< functions
26      double precision poa(nbf,nbf) !< [Input] The alpha occupation
27                                    !< functions penalty coefficients
28      double precision pob(nbf,nbf) !< [Input] The beta occupation
29                                    !< functions penalty coefficients
30c
31      double precision doa(nbf,nbf)  !< [In/Output] The alpha occupation
32                                     !< function gradient
33      double precision dob(nbf,nbf)  !< [In/Output] The beta occupation
34                                     !< function gradient
35      double precision dpoa(nbf,nbf) !< [In/Output] The alpha occupation
36                                     !< function penalty coefficient
37                                     !< gradient
38      double precision dpob(nbf,nbf) !< [In/Output] The beta occupation
39                                     !< function penalty coefficient
40                                     !< gradient
41c
42      double precision EPo !< [Output] The occupation function
43                           !< penalty energy
44c
45      double precision ta(nbf,nbf) !< [Scratch]
46      double precision tb(nbf,nbf) !< [Scratch]
47c
48      double precision tt
49c
50      double precision ddot
51      external         ddot
52c
53      integer ip, iq
54      integer ii
55c
56c     Compute the penalty function contributions
57c
58      EPo = 0.0d0
59      do ip = 1, nbf
60        do iq = 1, nbf
61          tt = ddot(nbf,oa(1,ip),1,oa(1,iq),1)
62          if (ip.eq.iq) tt = tt - 1.0d0
63          ta(ip,iq) = tt
64          EPo = EPo + poa(ip,iq)*(tt**2)
65          dpoa(ip,iq) = dpoa(ip,iq) + tt**2
66        enddo
67      enddo
68      do ip = 1, nbf
69        do iq = 1, nbf
70          tt = ddot(nbf,ob(1,ip),1,ob(1,iq),1)
71          if (ip.eq.iq) tt = tt - 1.0d0
72          tb(ip,iq) = tt
73          EPo = EPo + pob(ip,iq)*(tt**2)
74          dpob(ip,iq) = dpob(ip,iq) + tt**2
75        enddo
76      enddo
77c
78c     Compute derivatives wrt occupation function coefficients
79c
80      do iq = 1, nbf
81        do ip = 1, nbf
82          ta(ip,iq) = 2*ta(ip,iq)*poa(ip,iq)
83          tb(ip,iq) = 2*tb(ip,iq)*pob(ip,iq)
84        enddo
85      enddo
86      do ip = 1, nbf
87        do ii = 1, nbf
88          doa(ii,ip) = doa(ii,ip) + ddot(nbf,oa(ii,1),nbf,ta(1,ip),1)
89          doa(ii,ip) = doa(ii,ip) + ddot(nbf,oa(ii,1),nbf,ta(ip,1),nbf)
90        enddo
91      enddo
92      do ip = 1, nbf
93        do ii = 1, nbf
94          dob(ii,ip) = dob(ii,ip) + ddot(nbf,ob(ii,1),nbf,tb(1,ip),1)
95          dob(ii,ip) = dob(ii,ip) + ddot(nbf,ob(ii,1),nbf,tb(ip,1),nbf)
96        enddo
97      enddo
98c
99      end
100C>
101C> @}
102