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