1C> \ingroup nwpwxc 2C> @{ 3C> 4C> \file nwpwxc_c_opt.F 5C> The OPT correlation functional 6C> 7C> @} 8C> 9C> \ingroup nwpwxc_priv 10C> @{ 11C> 12C> \brief The OPT correlation functional 13C> 14C> The OPTimized correlation functional [1]. 15C> 16C> [1] N.C. Handy, A.J. Cohen, "Dynamic correlation", 17C> Mol. Phys. <b>99</b>, 607-615 (2001), DOI: 18C> <A HREF="https://doi.org/10.1080/00268970010023435"> 19C> 10.1080/00268970010023435</A>. 20c 21C$Id$ 22c 23 subroutine nwpwxc_c_opt(tol_rho,ipol,nq,wght,rho,rgamma,func, 24 & Amat, Cmat) 25 implicit none 26c 27#include "nwpwxc_param.fh" 28c 29 double precision tol_rho !< [Input] The lower limit on the density 30 integer nq !< [Input] The number of points 31 integer ipol !< [Input] The number of spin channels 32 double precision wght !< [Input] The weight of the functional 33c 34c Charge Density 35c 36 double precision rho(nq,*) !< [Input] The density 37c 38c Charge Density Gradient 39c 40 double precision rgamma(nq,*) !< [Input] The norm of the density gradients 41c 42c Sampling Matrices for the XC Potential & Energy 43c 44 double precision func(nq) !< [Output] The value of the functional 45 double precision Amat(nq,*) !< [Output] The derivative wrt rho 46 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 47c 48c References: 49c 50c Handy NC, Cohen AJ, Mol Phys 99 (7); 607-615 2001 51c 52 integer l_rho,k_rho,l_delrho,k_delrho 53 double precision c1,c2 54! parameter (c1=1.1015d0,c2=0.6625d0) 55 parameter (c1=1.d0,c2=0d0) 56c 57 integer iq ! counter over grid points 58 integer ii ! offset 59 integer iqt ! upper limit 60 integer num ! the number of grid points in current batch 61 integer maxp ! the maximum number of points in a batch 62 parameter (maxp = 16) 63 double precision funcl(maxp) 64 double precision rhol(maxp,2) 65 double precision rgammal(maxp,3) 66 double precision Amatl(maxp,2) 67 double precision Cmatl(maxp,3) 68 double precision fac 69c 70c*************************************************************************** 71c 72c 73 if (ipol.eq.1) then 74c 75c ======> SPIN-RESTRICTED <====== 76c 77 do iq = 1, nq, maxp 78 iqt = min(nq,iq+maxp-1) 79 num = iqt-iq+1 80c 81c c1*Ec[a,b] 82c 83 call dfill(maxp,0.0d0,funcl,1) 84 call dfill(maxp*2,0.0d0,rhol,1) 85 call dfill(maxp*3,0.0d0,rgammal,1) 86 call dfill(maxp*2,0.0d0,Amatl,1) 87 call dfill(maxp*3,0.0d0,Cmatl,1) 88 do ii = 0, num - 1 89 rhol(ii+1,R_T) = rho(iq+ii,R_T) 90 rgammal(ii+1,G_TT) = rgamma(iq+ii,G_TT) 91 enddo 92 fac = c1*wght 93 call nwpwxc_c_pw91lda(tol_rho,1,maxp,fac,rhol,funcl,Amatl) 94 call nwpwxc_c_p91(tol_rho,1,maxp,fac,rhol,rgammal,funcl, 95 + Amatl,Cmatl) 96 do ii = 0, num - 1 97 func(iq+ii) = func(iq+ii) + funcl(ii+1) 98 Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA) 99 Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA) 100 + + Cmatl(ii+1,D1_GAA) 101 Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB) 102 + + Cmatl(ii+1,D1_GAB) 103 enddo 104c 105c {(c2-c1)*Ec[a,0]}(a) + {(c2-c1)*Ec[0,b]}(a) 106c 107 fac = (c2-c1)*wght 108 call dfill(maxp,0.0d0,funcl,1) 109 call dfill(maxp*2,0.0d0,Amatl,1) 110 call dfill(maxp*3,0.0d0,Cmatl,1) 111 do ii = 0, num - 1 112 rhol(ii+1,R_A) = 0.5d0*rhol(ii+1,R_T) 113 rgammal(ii+1,G_AA) = 0.25d0*rgammal(ii+1,G_TT) 114 enddo 115 call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl) 116 call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl, 117 + Amatl,Cmatl) 118 do ii = 0, num - 1 119 func(iq+ii) = func(iq+ii) + 2*funcl(ii+1) 120 Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) 121 + + Amatl(ii+1,D1_RA) + Amatl(ii+1,D1_RB) 122 Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA) 123 + + (2.0d0)*Cmatl(ii+1,D1_GAA) 124 Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB) 125 + + (2.0d0)*Cmatl(ii+1,D1_GAB) 126 enddo 127 enddo 128c 129 else 130c 131c ======> SPIN-UNRESTRICTED <====== 132c 133 do iq = 1, nq, maxp 134 iqt = min(nq,iq+maxp-1) 135 num = iqt-iq+1 136c 137c c1*Ec[a,b] 138c 139 call dfill(maxp,0.0d0,funcl,1) 140 call dfill(maxp*2,0.0d0,rhol,1) 141 call dfill(maxp*3,0.0d0,rgammal,1) 142 call dfill(maxp*2,0.0d0,Amatl,1) 143 call dfill(maxp*3,0.0d0,Cmatl,1) 144 do ii = 0, num - 1 145 rhol(ii+1,R_A) = rho(iq+ii,R_A) 146 rhol(ii+1,R_B) = rho(iq+ii,R_B) 147 rgammal(ii+1,G_AA) = rgamma(iq+ii,G_AA) 148 rgammal(ii+1,G_AB) = rgamma(iq+ii,G_AB) 149 rgammal(ii+1,G_BB) = rgamma(iq+ii,G_BB) 150 enddo 151 fac = c1*wght 152 call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl) 153 call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl, 154 + Amatl,Cmatl) 155 do ii = 0, num - 1 156 func(iq+ii) = func(iq+ii) + funcl(ii+1) 157 Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA) 158 Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB) 159 Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA) 160 + + Cmatl(ii+1,D1_GAA) 161 Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB) 162 + + Cmatl(ii+1,D1_GAB) 163 Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB) 164 + + Cmatl(ii+1,D1_GBB) 165 enddo 166c 167c (c2-c1)*Ec[a,0] 168c 169 fac = (c2-c1)*wght 170 call dfill(maxp ,0.0d0,funcl,1) 171 call dfill(maxp*2,0.0d0,Amatl,1) 172 call dfill(maxp*3,0.0d0,Cmatl,1) 173 do ii = 0, num - 1 174 rhol(ii+1,R_A) = rho(iq+ii,R_A) 175 rhol(ii+1,R_B) = 0.0d0 176 rgammal(ii+1,G_AA) = rgamma(iq+ii,G_AA) 177 rgammal(ii+1,G_AB) = 0.0d0 178 rgammal(ii+1,G_BB) = 0.0d0 179 enddo 180 call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl) 181 call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl, 182 + Amatl,Cmatl) 183 do ii = 0, num - 1 184 func(iq+ii) = func(iq+ii) + funcl(ii+1) 185 Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA) 186 Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB) 187 Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA) 188 + + Cmatl(ii+1,D1_GAA) 189 Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB) 190 + + Cmatl(ii+1,D1_GAB) 191 Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB) 192 + + Cmatl(ii+1,D1_GBB) 193 enddo 194c 195c (c2-c1)*Ec[0,b] 196c 197 fac = (c2-c1)*wght 198 call dfill(maxp ,0.0d0,funcl,1) 199 call dfill(maxp*2,0.0d0,Amatl,1) 200 call dfill(maxp*3,0.0d0,Cmatl,1) 201 do ii = 0, num - 1 202 rhol(ii+1,R_A) = 0.0d0 203 rhol(ii+1,R_B) = rho(iq+ii,R_B) 204 rgammal(ii+1,G_AA) = 0.0d0 205 rgammal(ii+1,G_AB) = 0.0d0 206 rgammal(ii+1,G_BB) = rgamma(iq+ii,G_BB) 207 enddo 208 call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl) 209 call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl, 210 + Amatl,Cmatl) 211 do ii = 0, num - 1 212 func(iq+ii) = func(iq+ii) + funcl(ii+1) 213 Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA) 214 Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB) 215 Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA) 216 + + Cmatl(ii+1,D1_GAA) 217 Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB) 218 + + Cmatl(ii+1,D1_GAB) 219 Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB) 220 + + Cmatl(ii+1,D1_GBB) 221 enddo 222 enddo 223c 224 endif 225c 226 return 227 end 228C> 229C> @} 230