1C> \ingroup wfn1 2C> @{ 3C> 4C> \brief Orthogonalize \f$\Psi_I\f$ on the set of orthogonal states 5C> \f$\Psi_J,\; J=1,I-1\f$ 6C> 7C> Assume we have an orthonormal set of N-electron wavefunctions 8C> \f$\left\{\Psi_J\right\}\f$ where \f$1\le J \lt I\f$. 9C> Assume also that \f$\Psi_I\f$ is normalized. We want to 10C> orthonormalize \f$\Psi_I\f$ onto that set. Because of the way the 11C> superposition principle works for WFN1 style wavefunctions it is 12C> most efficient to use Gramm-Schmidt (even though that is less 13C> accurate than modified Gramm-Schmidt). 14C> 15C> The approach then is to: 16C> 17C> * Calculate \f$\forall_{J<I}, c_{IJ} = \left\langle\Psi_I\right|\left.\Psi_J\right\rangle\f$. 18C> The new wavefunction then is \f$\Psi'_I = A\left(\Psi_I-\sum_{J=1}^{I-1}c_{IJ}\Psi_J\right)\f$ 19C> where \f$A=1/\sqrt{1+\sum_J(c_{IJ}^2-2c_{IJ})}\f$. 20C> 21C> * Calculate the 1-electron density matrix 22C> \f$D = \left|\Psi'_I\right\rangle\left\langle\Psi'_I\right|\f$. 23C> 24C> * Diagonalize \f$D\f$ to obtain new natural orbitals 25C> \f$\left\{^I\phi^{D\sigma}\right\}\f$ and associated occupation 26C> numbers. 27C> 28C> * Update the occupation functions to match the new occupation 29C> numbers. 30C> 31C> In particular the requirement to construct and diagonalize a 32C> density matrix makes the simultaneous orthonormalization on all 33C> lower states attractive (as opposed to a modified Gramm-Schmidt). 34C> 35 subroutine wfn1_orthogonalize(nst,nbf,nea,neb,ov,r_noa,r_nob, 36 & r_ofa,r_ofb,ovla,ovlb,cij, 37 & docca,doccb,da,db) 38 implicit none 39c 40 integer nst !< [Input] The number of states (\f$=I\f$) 41 integer nbf !< [Input] The number of basis functions 42 integer nea !< [Input] The number of \f$\alpha\f$-electrons 43 integer neb !< [Input] The number of \f$\beta\f$-electrons 44c 45 double precision ov(nbf,nbf) !< [Input] The overlap matrix 46 double precision r_noa(nbf,nbf,nst) !< [In/Output] The 47 !< \f$\alpha\f$ natural orbitals 48 double precision r_nob(nbf,nbf,nst) !< [In/Output] The 49 !< \f$\beta\f$ natural orbitals 50 double precision r_ofa(nbf,nbf,nst) !< [In/Output] The 51 !< \f$\alpha\f$ occupation functions 52 double precision r_ofb(nbf,nbf,nst) !< [In/Output] The 53 !< \f$\beta\f$ occupation functions 54 double precision ovla(nea,nea) !< [Scratch] The \f$\alpha\f$ 55 !< overlaps 56 double precision ovlb(neb,neb) !< [Scratch] The \f$\beta\f$ 57 !< overlaps 58 double precision cij(nst) !< [Scratch] The state overlaps 59 double precision docca(nbf) !< [Scratch] The \f$\alpha\f$ 60 !< occupation numbers 61 double precision doccb(nbf) !< [Scratch] The \f$\beta\f$ 62 !< occupation numbers 63 double precision da(nbf,nbf) !< [Scratch] The \f$\alpha\f$ 64 !< density matrix 65 double precision db(nbf,nbf) !< [Scratch] The \f$\beta\f$ 66 !< density matrix 67c 68 integer ii !< Counter over natural orbitals 69 integer ir !< Counter over occupation functions 70 integer ia !< Counter over basis functions 71 integer ib !< Counter over basis functions 72c 73 integer GENEIG !< Parameter for Generalized Eigenvalue Problem 74 parameter(GENEIG = 1) 75c 76c First query LAPACK DSYGV for memory requirements 77c 78 lwork = -1 79 info = 0 80 call dsygv(GENEIG,'V','U',nbf,da,nbf,db,nbf,docca,ww,lwork,info) 81 if (info.ne.0) then 82 call errquit("wfn1_orthogonalize: a: dsygv failed ",info,UERR) 83 endif 84 lwork = ww+0.2d0 85 if (.not.ma_push_get(MT_DBL,lwork,"work dsygv",l_w,k_w)) then 86 call errquit("wfn1_orthogonalize: alloc work failed", 87 + ma_sizeof(MT_DBL,lwork,MT_BYTE),MA_ERR) 88 endif 89c 90c Calculate the overlaps and the normalization factor 91c 92 cij(nst) = 1.0d0 93 do jst = 1, nst-1 94 call wfn1_ovlp_2_states(nea,nbf,r_ofa(1,1,jst),r_noa(1,1,jst), 95 + r_ofa(1,1,nst),r_noa(1,1,nst),ov, 96 + dbl_mb(k_m1),dbl_mb(k_m2),ovla) 97 call wfn1_ovlp_2_states(neb,nbf,r_ofb(1,1,jst),r_nob(1,1,jst), 98 + r_ofb(1,1,nst),r_nob(1,1,nst),ov, 99 + dbl_mb(k_m1),dbl_mb(k_m2),ovlb) 100 Sa_b = wfn1_overlap_bo(nea,ovla,dbl_mb(k_m1),dbl_mb(k_m2)) 101 Sb_b = wfn1_overlap_bo(neb,ovlb,dbl_mb(k_m1),dbl_mb(k_m2)) 102 cij(jst) = Sa_b*Sb_b 103 enddo 104 aa = 1.0d0 105 do jst = 1, nst-1 106 aa = aa + cij(jst)**2-2.0d0*cij(jst) 107 enddo 108c 109c Construct the density matrix of the orthonormalize state 110c 111c Do diagonal terms 112c 113 call wfn1_add_dmat(nbf,nea,neb,1.0d0, 114 + r_noa(1,1,nst),r_nob(1,1,nst), 115 + r_ofa(1,1,nst),r_ofb(1,1,nst), 116 + 0.0d0,da,db,docca,doccb) 117 do jst = 1, nst-1 118 call wfn1_add_dmat(nbf,nea,neb,cij(jst)**2, 119 + r_noa(1,1,jst),r_nob(1,1,jst), 120 + r_ofa(1,1,jst),r_ofb(1,1,jst), 121 + 1.0d0,da,db,docca,doccb) 122 enddo 123c 124c Do cross terms 125c 126 do jst = 1, nst-1 127c 128c Do <Y_J|Y_N> 129c 130 call wfn1_ovlp_2_states(nea,nbf,r_ofa(1,1,jst),r_noa(1,1,jst), 131 + r_ofa(1,1,nst),r_noa(1,1,nst),ov, 132 + dbl_mb(k_m1),dbl_mb(k_m2),ovla) 133 call wfn1_ovlp_2_states(neb,nbf,r_ofb(1,1,jst),r_nob(1,1,jst), 134 + r_ofb(1,1,nst),r_nob(1,1,nst),ov, 135 + dbl_mb(k_m1),dbl_mb(k_m2),ovlb) 136 call wfn1_overlap_ldr(nea,ovla,lpa,rpa,lma,rma,va,ipa,tmp) 137 call wfn1_overlap_ldr(neb,ovlb,lpb,rpb,lmb,rmb,vb,ipb,tmp) 138 call wfn1_LR2d(nbf,nea,lpa,rpa,lma,rma, 139 + r_ofa(1,1,jst),r_ofa(1,1,nst), 140 + r_ofa2(1,1,jst),r_ofa2(1,1,nst),tmp) 141 call wfn1_LR2d(nbf,neb,lpb,rpb,lmb,rmb, 142 + r_ofb(1,1,jst),r_ofb(1,1,nst), 143 + r_ofb2(1,1,jst),r_ofb2(1,1,nst),tmp) 144 call wfn1_add_tdmat(nbf,nea,neb,-cij(jst)*ipa*ipb,va,vb, 145 + r_noa(1,1,jst),r_nob(1,1,jst), 146 + r_ofa2(1,1,jst),r_ofb2(1,1,jst), 147 + r_noa(1,1,nst),r_nob(1,1,nst), 148 + r_ofa2(1,1,nst),r_ofb2(1,1,nst),1.0d0, 149 + da,db,docca,doccb) 150c 151c Do <Y_N|Y_J> 152c 153 call wfn1_ovlp_2_states(nea,nbf,r_ofa(1,1,nst),r_noa(1,1,nst), 154 + r_ofa(1,1,jst),r_noa(1,1,jst),ov, 155 + dbl_mb(k_m1),dbl_mb(k_m2),ovla) 156 call wfn1_ovlp_2_states(neb,nbf,r_ofb(1,1,nst),r_nob(1,1,nst), 157 + r_ofb(1,1,jst),r_nob(1,1,jst),ov, 158 + dbl_mb(k_m1),dbl_mb(k_m2),ovlb) 159 call wfn1_overlap_ldr(nea,ovla,lpa,rpa,lma,rma,va,ipa,tmp) 160 call wfn1_overlap_ldr(neb,ovlb,lpb,rpb,lmb,rmb,vb,ipb,tmp) 161 call wfn1_LR2d(nbf,nea,lpa,rpa,lma,rma, 162 + r_ofa(1,1,nst),r_ofa(1,1,jst), 163 + r_ofa2(1,1,nst),r_ofa2(1,1,jst),tmp) 164 call wfn1_LR2d(nbf,neb,lpb,rpb,lmb,rmb, 165 + r_ofb(1,1,nst),r_ofb(1,1,jst), 166 + r_ofb2(1,1,nst),r_ofb2(1,1,jst),tmp) 167 call wfn1_add_tdmat(nbf,nea,neb,-cij(jst)*ipa*ipb,va,vb, 168 + r_noa(1,1,nst),r_nob(1,1,nst), 169 + r_ofa2(1,1,nst),r_ofb2(1,1,nst), 170 + r_noa(1,1,jst),r_nob(1,1,jst), 171 + r_ofa2(1,1,jst),r_ofb2(1,1,jst),1.0d0, 172 + da,db,docca,doccb) 173c 174 enddo 175c 176c Normalize the new density matrices 177c 178 call dscal(nbf*nbf,1.0d0/aa,da,1) 179 call dscal(nbf*nbf,1.0d0/aa,db,1) 180c 181c Calculate new natural orbitals 182c 183 call dcopy(nbf*nbf,da,1,r_noa(1,1,nst),1) 184 call dcopy(nbf*nbf,ov,1,da,1) 185 call dsygv(GENEIG,'V','U',nbf,r_noa(1,1,nst),nbf,da,nbf,docca, 186 + dbl_mb(k_w),lwork,info) 187 if (info.ne.0) then 188 call errquit("wfn1_orthogonalize: b: dsygv failed ",info,UERR) 189 endif 190 call dcopy(nbf*nbf,db,1,r_nob(1,1,nst),1) 191 call dcopy(nbf*nbf,ov,1,db,1) 192 call dsygv(GENEIG,'V','U',nbf,r_nob(1,1,nst),nbf,db,nbf,doccb, 193 + dbl_mb(k_w),lwork,info) 194 if (info.ne.0) then 195 call errquit("wfn1_orthogonalize: c: dsygv failed ",info,UERR) 196 endif 197c 198c Update the occupation to re-establish the wavefunction 199c 200 call wfn1_match_occ(.false.,nbf,nea,docca,r_ofa(1,1,nst),m1,m2, 201 + m3,m4,v1,.false.) 202 call wfn1_match_occ(.false.,nbf,neb,doccb,r_ofb(1,1,nst),m1,m2, 203 + m3,m4,v1,.false.) 204c 205 end 206C> 207C> @} 208