TABLE OF CONTENTS
ABINIT/vgh_rho [ Functions ]
NAME
vgh_rho
FUNCTION
The general procedure to obtain the value, the gradient and the hessian of the density of electrons in the point vv (in cart.coord).
COPYRIGHT
Copyright (C) 2002-2017 ABINIT group (PCasek,FF,XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt . WARNING This file does not follow the ABINIT coding rules (yet)
INPUTS
vv(3)=position
chs 1 only valence density
2 only core density
0 total density
-2 iat, ipos are nulify and ignored
-1 iat,ipos = index of atom if vv is inside
the "core sphere (rminl)", 0 otherwise
OUTPUT
rho,grho(3),hrho(3,3) - density, gradient of density, hessian of density
(cart. coord)
iat, ipos - index of the nearest atom (except chs < 0 see above)
rdmin - the distance to the nearest atom
SIDE EFFECTS
This routine also works on the data contained in the defs_aimprom and defs_aimfields modules
PARENTS
addout,aim_follow,cpdrv,critic,critics,integrho,onestep,plint,rsurf
CHILDREN
bschg1,bschg2
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 53 subroutine vgh_rho(vv,rho,grho,hrho,rdmin,iat,ipos,chs) 54 55 use m_profiling_abi 56 57 use defs_basis 58 use defs_parameters 59 use defs_aimprom 60 use defs_aimfields 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'vgh_rho' 66 use interfaces_63_bader, except_this_one => vgh_rho 67 !End of the abilint section 68 69 implicit none 70 71 !Arguments ------------------------------------ 72 !scalars 73 integer,intent(in) :: chs 74 integer,intent(inout) :: iat,ipos 75 real(dp),intent(out) :: rdmin,rho 76 !arrays 77 real(dp),intent(in) :: vv(3) 78 real(dp),intent(out) :: grho(3),hrho(3,3) 79 80 !Local variables ------------------------------ 81 !scalars 82 integer :: ii,inmax,inmin,inx,jj,kk,ll,nn,oii,omm,onn 83 integer :: selct 84 ! real(dp),save :: cumul_cpu=0.0_dp,cumul_cpu_old=0.0_dp 85 real(dp),save :: tcpui,tcpuo,twalli 86 real(dp),save :: twallo 87 real(dp) :: aa,bb,cc,cgrad1_rr_inv,coeff,dd,rr,rr2,rr_inv 88 real(dp) :: rrad2_nn,rrad_nn,ss,uu,uu_inv,val,vt1,vt2,vt3,vw1,vw2 89 ! real(dp) :: ss_inv 90 real(dp) :: vw3 91 !arrays 92 integer :: indx(3),inii(4,3) 93 real(dp) :: cgrad(3),ches(3,3),cof(2,3),ddstar(6),ddu(2),grd(4) 94 real(dp) :: hh(4,2),hrh(2),lder(4),pom2sq(2,3),pomsq(2,3) 95 real(dp) :: rhstar(6),sqder(6,4),sqvlr(6,4),trsf(3,3),xx(3) 96 real(dp),pointer :: ptddx(:,:,:),ptddy(:,:,:),ptddz(:,:,:),ptrho(:,:,:) 97 98 !************************************************************************ 99 tcpui=0.0_dp 100 tcpuo=0.0_dp 101 twalli=0.0_dp 102 twallo=0.0_dp 103 104 nullify(ptddx,ptddy,ptddz,ptrho) 105 106 selct=chs 107 108 if (selct/=2) then 109 110 ! call timein(tcpui,twalli) 111 112 ! TRANSFORMATION TO THE REDUCED COORD. 113 114 xx(:)=vv(:) 115 call bschg1(xx,-1) 116 117 ! call timein(tcpuo,twallo) 118 ! cumul_cpu=cumul_cpu+(tcpuo-tcpui) 119 120 ! REDUCTION TO THE PRIMITIVE CELL 121 122 do ii=1,3 123 if (xx(ii) >= one-tol12 ) then 124 xx(ii)=xx(ii)-aint(xx(ii)) 125 elseif (xx(ii) < -tol12 ) then 126 xx(ii)=xx(ii)-floor(xx(ii)) 127 end if 128 end do 129 130 131 ! DETERMINATION OF THE INDEX IN THE GRID 132 133 do ii=1,3 134 indx(ii)=aint(xx(ii)*ngfft(ii)) 135 bb=(xx(ii)-indx(ii)*dix(ii))*ngfft(ii) 136 if (indx(ii)==ngfft(ii)) then 137 indx(ii)=1 138 xx(ii)=0._dp 139 else 140 indx(ii)=indx(ii)+1 141 end if 142 143 ! Explicit handling to avoid numeric problems 144 145 if (bb > 1._dp+tol12 ) then 146 cof(1,ii)=0._dp 147 cof(2,ii)=1._dp 148 elseif (bb < -tol12 ) then 149 cof(1,ii)=1._dp 150 cof(2,ii)=0._dp 151 else 152 cof(1,ii)=1._dp-bb 153 cof(2,ii)=bb 154 end if 155 end do 156 157 ! 3D INTERPOLATION OF THE VALENCE DENSITY 158 159 ! determination of the values of density and of its second derivative 160 ! at the "star" = constructed at vv with primitive directions 161 ! To interpolation the values at the faces of the grid cell are needed 162 163 rhstar(:)=0._dp 164 sqder(:,:)=0._dp 165 sqvlr(:,:)=0._dp 166 ddstar(:)=0._dp 167 pomsq(:,:)=0._dp 168 pom2sq(:,:)=0._dp 169 170 oii=1; onn=1; omm=1 171 if (indx(1)==ngfft(1)) oii=1-ngfft(1) 172 if (indx(2)==ngfft(2)) onn=1-ngfft(2) 173 if (indx(3)==ngfft(3)) omm=1-ngfft(3) 174 175 ! the values in the corners of the grid cell 176 177 ptddx=>ddx(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 178 ptddy=>ddy(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 179 ptddz=>ddz(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 180 ptrho=>dvl(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 181 182 ! the coefficients for spline interpolation of density and its derivation 183 do ii=1,3 184 do jj=1,2 185 pomsq(jj,ii)=(cof(jj,ii)*cof(jj,ii)*cof(jj,ii)-cof(jj,ii))/6._dp*dix(ii)*dix(ii) 186 pom2sq(jj,ii)=(3._dp*cof(jj,ii)*cof(jj,ii)-1._dp)/6._dp*dix(ii) 187 if (jj==1) pom2sq(jj,ii)=-pom2sq(jj,ii) 188 end do 189 end do 190 191 192 do ii=1,2 193 do jj=1,2 194 do kk=1,2 195 ddstar(ii)=ddstar(ii)+cof(jj,2)*cof(kk,3)*ptddx(ii,jj,kk) 196 ddstar(ii+2)=ddstar(ii+2)+cof(jj,3)*cof(kk,1)*ptddy(kk,ii,jj) 197 ddstar(ii+4)=ddstar(ii+4)+cof(jj,1)*cof(kk,2)*ptddz(jj,kk,ii) 198 sqder(ii,jj)=sqder(ii,jj)+cof(kk,2)*ptddz(ii,kk,jj) 199 sqder(ii,jj+2)=sqder(ii,jj+2)+cof(kk,3)*ptddy(ii,jj,kk) 200 sqder(ii+2,jj)=sqder(ii+2,jj)+cof(kk,3)*ptddx(jj,ii,kk) 201 sqder(ii+2,jj+2)=sqder(ii+2,jj+2)+cof(kk,1)*ptddz(kk,ii,jj) 202 sqder(ii+4,jj)=sqder(ii+4,jj)+cof(kk,1)*ptddy(kk,jj,ii) 203 sqder(ii+4,jj+2)=sqder(ii+4,jj+2)+cof(kk,2)*ptddx(jj,kk,ii) 204 sqvlr(ii,jj)=sqvlr(ii,jj)+cof(kk,2)*ptrho(ii,kk,jj)+pomsq(kk,2)*ptddy(ii,kk,jj) 205 sqvlr(ii,jj+2)=sqvlr(ii,jj+2)+cof(kk,3)*ptrho(ii,jj,kk)+pomsq(kk,3)*ptddz(ii,jj,kk) 206 sqvlr(ii+2,jj+2)=sqvlr(ii+2,jj+2)+cof(kk,1)*ptrho(kk,ii,jj)+pomsq(kk,1)*ptddx(kk,ii,jj) 207 end do 208 end do 209 end do 210 211 do ii=1,2 212 do jj=1,2 213 sqvlr(ii+2,jj)=sqvlr(jj,ii+2) 214 sqvlr(ii+4,jj)=sqvlr(jj+2,ii+2) 215 sqvlr(ii+4,jj+2)=sqvlr(jj,ii) 216 end do 217 end do 218 219 do ii=1,2 220 do jj=1,2 221 rhstar(ii)=rhstar(ii)+cof(jj,3)*sqvlr(ii,jj)+pomsq(jj,3)*sqder(ii,jj)+& 222 & cof(jj,2)*sqvlr(ii,jj+2)+pomsq(jj,2)*sqder(ii,jj+2) 223 rhstar(ii+2)=rhstar(ii+2)+cof(jj,1)*sqvlr(ii+2,jj)+pomsq(jj,1)*sqder(ii+2,jj)+& 224 & cof(jj,3)*sqvlr(ii+2,jj+2)+pomsq(jj,3)*sqder(ii+2,jj+2) 225 rhstar(ii+4)=rhstar(ii+4)+cof(jj,2)*sqvlr(ii+4,jj)+pomsq(jj,2)*sqder(ii+4,jj)+& 226 & cof(jj,1)*sqvlr(ii+4,jj+2)+pomsq(jj,1)*sqder(ii+4,jj+2) 227 end do 228 end do 229 rhstar(:)=rhstar(:)/2._dp 230 231 rho=0._dp 232 grho(:)=0._dp 233 hrho(:,:)=0._dp 234 kk=1; nn=1 235 do ii=1,5,2 236 do jj=1,2 237 nn=-nn 238 rho=rho+cof(jj,kk)*rhstar(ii+jj-1)+pomsq(jj,kk)*ddstar(ii+jj-1) 239 grho(kk)=grho(kk)+pom2sq(jj,kk)*ddstar(ii+jj-1) 240 hrho(kk,kk)=hrho(kk,kk)+cof(jj,kk)*ddstar(ii+jj-1) 241 grho(kk)=grho(kk)+nn*rhstar(ii+jj-1)/dix(kk) 242 end do 243 kk=kk+1 244 end do 245 rho=rho/3._dp 246 247 ! Off-diagonal elements of the hessian 248 249 ! for the speed reasons the polynomial interpolation 250 ! for second derivation fields is used in this case 251 ! but the last step is always done by spline interpolation. 252 253 254 do ii=1,3 255 do jj=-1,2 256 inii(jj+2,ii)=indx(ii)+jj 257 if (inii(jj+2,ii) < 1) inii(jj+2,ii)=inii(jj+2,ii)+ngfft(ii) 258 if (inii(jj+2,ii) > ngfft(ii)) inii(jj+2,ii)=inii(jj+2,ii)-ngfft(ii) 259 end do 260 end do 261 262 ! Not very nice 263 264 do ii=1,3 265 select case (ii) 266 case (1) 267 do jj=1,4 268 ddu(1)=cof(1,2)*ddz(inii(jj,1),inii(2,2),inii(2,3))+cof(2,2)*ddz(inii(jj,1),inii(3,2),inii(2,3)) 269 ddu(2)=cof(1,2)*ddz(inii(jj,1),inii(2,2),inii(3,3))+cof(2,2)*ddz(inii(jj,1),inii(3,2),inii(3,3)) 270 hrh(1)=cof(1,2)*dvl(inii(jj,1),inii(2,2),inii(2,3))+cof(2,2)*dvl(inii(jj,1),inii(3,2),inii(2,3))+& 271 & pomsq(1,2)*ddy(inii(jj,1),inii(2,2),inii(2,3))+pomsq(2,2)*ddy(inii(jj,1),inii(3,2),inii(2,3)) 272 hrh(2)=cof(1,2)*dvl(inii(jj,1),inii(2,2),inii(3,3))+cof(2,2)*dvl(inii(jj,1),inii(3,2),inii(3,3))+& 273 & pomsq(1,2)*ddy(inii(jj,1),inii(2,2),inii(3,3))+pomsq(2,2)*ddy(inii(jj,1),inii(3,2),inii(3,3)) 274 hh(jj,2)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 275 276 ddu(1)=cof(1,3)*ddy(inii(jj,1),inii(2,2),inii(2,3))+cof(2,3)*ddy(inii(jj,1),inii(2,2),inii(3,3)) 277 ddu(2)=cof(1,3)*ddy(inii(jj,1),inii(3,2),inii(2,3))+cof(2,3)*ddy(inii(jj,1),inii(3,2),inii(3,3)) 278 hrh(1)=cof(1,3)*dvl(inii(jj,1),inii(2,2),inii(2,3))+cof(2,3)*dvl(inii(jj,1),inii(2,2),inii(3,3))+& 279 & pomsq(1,3)*ddz(inii(jj,1),inii(2,2),inii(2,3))+pomsq(2,3)*ddz(inii(jj,1),inii(2,2),inii(3,3)) 280 hrh(2)=cof(1,3)*dvl(inii(jj,1),inii(3,2),inii(2,3))+cof(2,3)*dvl(inii(jj,1),inii(3,2),inii(3,3))+& 281 & pomsq(1,3)*ddz(inii(jj,1),inii(3,2),inii(2,3))+pomsq(2,3)*ddz(inii(jj,1),inii(3,2),inii(3,3)) 282 hh(jj,1)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 283 end do 284 case (2) 285 do jj=1,4 286 ddu(1)=cof(1,3)*ddx(inii(2,1),inii(jj,2),inii(2,3))+cof(2,3)*ddx(inii(2,1),inii(jj,2),inii(3,3)) 287 ddu(2)=cof(1,3)*ddx(inii(3,1),inii(jj,2),inii(2,3))+cof(2,3)*ddx(inii(3,1),inii(jj,2),inii(3,3)) 288 hrh(1)=cof(1,3)*dvl(inii(2,1),inii(jj,2),inii(2,3))+cof(2,3)*dvl(inii(2,1),inii(jj,2),inii(3,3))+& 289 & pomsq(1,3)*ddz(inii(2,1),inii(jj,2),inii(2,3))+pomsq(2,3)*ddz(inii(2,1),inii(jj,2),inii(3,3)) 290 hrh(2)=cof(1,3)*dvl(inii(3,1),inii(jj,2),inii(2,3))+cof(2,3)*dvl(inii(3,1),inii(jj,2),inii(3,3))+& 291 & pomsq(1,3)*ddz(inii(3,1),inii(jj,2),inii(2,3))+pomsq(2,3)*ddz(inii(3,1),inii(jj,2),inii(3,3)) 292 hh(jj,2)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 293 294 ddu(1)=cof(1,1)*ddz(inii(2,1),inii(jj,2),inii(2,3))+cof(2,1)*ddz(inii(3,1),inii(jj,2),inii(2,3)) 295 ddu(2)=cof(1,1)*ddz(inii(2,1),inii(jj,2),inii(3,3))+cof(2,1)*ddz(inii(3,1),inii(jj,2),inii(3,3)) 296 hrh(1)=cof(1,1)*dvl(inii(2,1),inii(jj,2),inii(2,3))+cof(2,1)*dvl(inii(3,1),inii(jj,2),inii(2,3))+& 297 & pomsq(1,1)*ddx(inii(2,1),inii(jj,2),inii(2,3))+pomsq(2,1)*ddx(inii(3,1),inii(jj,2),inii(2,3)) 298 hrh(2)=cof(1,1)*dvl(inii(2,1),inii(jj,2),inii(3,3))+cof(2,1)*dvl(inii(3,1),inii(jj,2),inii(3,3))+& 299 & pomsq(1,1)*ddx(inii(2,1),inii(jj,2),inii(3,3))+pomsq(2,1)*ddx(inii(3,1),inii(jj,2),inii(3,3)) 300 hh(jj,1)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 301 end do 302 case (3) 303 do jj=1,4 304 ddu(1)=cof(1,1)*ddy(inii(2,1),inii(2,2),inii(jj,3))+cof(2,1)*ddy(inii(3,1),inii(2,2),inii(jj,3)) 305 ddu(2)=cof(1,1)*ddy(inii(2,1),inii(3,2),inii(jj,3))+cof(2,1)*ddy(inii(3,1),inii(3,2),inii(jj,3)) 306 hrh(1)=cof(1,1)*dvl(inii(2,1),inii(2,2),inii(jj,3))+cof(2,1)*dvl(inii(3,1),inii(2,2),inii(jj,3))+& 307 & pomsq(1,1)*ddx(inii(2,1),inii(2,2),inii(jj,3))+pomsq(2,1)*ddx(inii(3,1),inii(2,2),inii(jj,3)) 308 hrh(2)=cof(1,1)*dvl(inii(2,1),inii(3,2),inii(jj,3))+cof(2,1)*dvl(inii(3,1),inii(3,2),inii(jj,3))+& 309 & pomsq(1,1)*ddx(inii(2,1),inii(3,2),inii(jj,3))+pomsq(2,1)*ddx(inii(3,1),inii(3,2),inii(jj,3)) 310 hh(jj,2)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 311 312 ddu(1)=cof(1,2)*ddx(inii(2,1),inii(2,2),inii(jj,3))+cof(2,2)*ddx(inii(2,1),inii(3,2),inii(jj,3)) 313 ddu(2)=cof(1,2)*ddx(inii(3,1),inii(2,2),inii(jj,3))+cof(2,2)*ddx(inii(3,1),inii(3,2),inii(jj,3)) 314 hrh(1)=cof(1,2)*dvl(inii(2,1),inii(2,2),inii(jj,3))+cof(2,2)*dvl(inii(2,1),inii(3,2),inii(jj,3))+& 315 & pomsq(1,2)*ddy(inii(2,1),inii(2,2),inii(jj,3))+pomsq(2,2)*ddy(inii(2,1),inii(3,2),inii(jj,3)) 316 hrh(2)=cof(1,2)*dvl(inii(3,1),inii(2,2),inii(jj,3))+cof(2,2)*dvl(inii(3,1),inii(3,2),inii(jj,3))+& 317 & pomsq(1,2)*ddy(inii(3,1),inii(2,2),inii(jj,3))+pomsq(2,2)*ddy(inii(3,1),inii(3,2),inii(jj,3)) 318 hh(jj,1)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 319 end do 320 end select 321 do jj=-2,1 322 grd(jj+3)=(indx(ii)+jj)*dix(ii) 323 end do 324 325 ! write(std_out,'("hh: ",/,4F16.8,/,4F16.8)') ((hh(kk,jj),kk=1,4),jj=1,2) 326 ! write(std_out,'("grad: ",3F16.8)') (grho(kk),kk=1,3) 327 ! write(std_out,'("dix: ",3F16.8)') (dix(kk),kk=1,3) 328 ! write(std_out,'("grd: ",4F16.8)') (grd(kk),kk=1,4) 329 ! write(std_out,'("inii: ",4I4)') (inii(kk,ii),kk=1,4) 330 331 do jj=1,2 332 333 ! polynomial interpolation 334 335 do kk=1,3 336 do ll=4,kk+1,-1 337 hh(ll,jj)=(hh(ll,jj)-hh(ll-1,jj))/(grd(ll)-grd(ll-1)) 338 end do 339 end do 340 lder(4)=hh(4,jj) 341 do kk=3,1,-1 342 lder(kk)=hh(kk,jj)+(xx(ii)-grd(kk))*lder(kk+1) 343 end do 344 do kk=1,2 345 do ll=3,kk+1,-1 346 lder(ll)=lder(ll)+(xx(ii)-grd(ll-kk))*lder(ll+1) 347 end do 348 end do 349 nn=ii+jj 350 if (nn > 3) nn=nn-3 351 hrho(ii,nn)=hrho(ii,nn)+lder(2) 352 hrho(nn,ii)=hrho(nn,ii)+lder(2) 353 end do 354 end do 355 356 ! averaging of the mixed derivations obtained in different order 357 358 do ii=1,3 359 do jj=1,3 360 if (ii /= jj) hrho(ii,jj)=hrho(ii,jj)/2._dp 361 end do 362 end do 363 364 365 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 366 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 367 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 368 ! stop 369 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 370 ! write(std_out,'(":GRAD pred tr ",3F16.8)') grho 371 ! write(std_out,'(":HESSIAN pred tr",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 372 373 374 ! Transformation back to Cart. coordonnes 375 376 call bschg1(grho,2) 377 call bschg2(hrho,2) 378 379 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 380 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 381 ! stop 382 383 nullify(ptddx,ptddy,ptddz,ptrho) 384 385 if (selct==1) return 386 387 end if 388 389 !write(51,'(":GRADv ",3F16.8)') grho 390 !write(52,'(":LAPv ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 391 !write(52,'(":HESNv ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 392 393 !INTERPOLATION OF THE CORE DENSITY 394 395 if (selct/=1) then 396 397 if (selct==2) then 398 grho(:)=0._dp 399 hrho(:,:)=0._dp 400 rho=0._dp 401 end if 402 403 ! SEARCH OF THE NEIGHBOUR ATOMS 404 405 if (selct /= -2) then 406 iat=0 407 ipos=0 408 end if 409 rdmin=20._dp 410 411 do jj=1,natom 412 nn=typat(jj) 413 rrad_nn=rrad(corlim(nn),nn) 414 rrad2_nn=rrad_nn*rrad_nn 415 vw1=vv(1)-xatm(1,jj) 416 vw2=vv(2)-xatm(2,jj) 417 vw3=vv(3)-xatm(3,jj) 418 419 do kk=1,nnpos 420 421 vt1=vw1-atp(1,kk) 422 vt2=vw2-atp(2,kk) 423 vt3=vw3-atp(3,kk) 424 rr2=vt1*vt1+vt2*vt2+vt3*vt3 425 426 ! rr=vnorm(vt,0) 427 428 ! Only contribution > rhocormin (adhoc.f90) are considered 429 430 if (rr2 < rrad2_nn .and.(.not.((selct==-2).and.(iat==jj).and.(ipos==kk)))) then 431 ! if (rr /= 0.0_dp) then ! XG020629 : never test a real number against zero (not portable) 432 if (rr2 > 1.0d-28) then ! SEARCH INDEX 433 434 rr=sqrt(rr2) 435 rr_inv=1.0_dp/rr 436 437 if (rr < rrad(1,nn)) then 438 inx=-1 439 elseif (rr > rrad(ndat(nn),nn)) then 440 inx=ndat(nn) 441 else 442 ! Find the index of the radius by bissection 443 inmin=1 444 inmax=ndat(nn) 445 inx=1 446 do 447 if(inmax-inmin==1)exit 448 inx=(inmin+inmax)/2 449 if(rr>=rrad(inx,nn))then 450 inmin=inx 451 else 452 inmax=inx 453 end if 454 end do 455 inx=inmin 456 457 ! XG020629 : old coding, slower 458 ! inx=0 459 ! do while (rr >= rrad(inx+1,nn)) 460 ! inx=inx+1 461 ! end do 462 463 end if 464 465 ! Transformation matrix radial -> cart. coord 466 ss=sqrt(vt1*vt1+vt2*vt2) 467 ! if (ss /=0._dp) then ! XG020629 : never test a real number against zero (not portable) 468 if (ss*ss > 1.0d-28) then ! ss non-zero 469 ! XG020629 : very strange : only trsf(:,1) is needed in what follows ? ! 470 ! ss_inv=1.0_dp/ss 471 trsf(1,1)=vt1*rr_inv 472 ! trsf(1,2)=-vt2*ss_inv 473 ! trsf(1,3)=vt3*vt1*rr_inv*ss_inv 474 trsf(2,1)=vt2*rr_inv 475 ! trsf(2,2)=vt1*ss_inv 476 ! trsf(2,3)=vt3*vt2*rr_inv*ss_inv 477 trsf(3,1)=vt3*rr_inv 478 ! trsf(3,2)=0._dp 479 ! trsf(3,3)=-ss*rr_inv 480 ! XG020629 Not needed 481 ! do ii=1,3 482 ! do ll=1,3 483 ! ches(ii,ll)=0._dp 484 ! end do 485 ! cgrad(ii)=0._dp 486 ! end do 487 else ! ss zero 488 do ii=1,3 489 do ll=1,3 490 trsf(ii,ll)=0._dp 491 end do 492 trsf(ii,4-ii)=1._dp 493 end do 494 end if ! ss zero or non-zero 495 496 if (inx == -1) then ! LEFT EXTRAPOLATION y=a*x^2+b (a<0) 497 val=sp2(1,nn)*0.5_dp*rr*rr/rrad(1,nn)+crho(1,nn)-sp2(1,nn)*rrad(1,nn)*0.5_dp 498 cgrad(1)=sp2(1,nn)*rr/rrad(1,nn) 499 ches(1,1)=sp2(1,nn)/rrad(1,nn) 500 elseif (inx == ndat(nn) ) then ! RIGHT EXTRAPOLATION y=a*exp(b*x) 501 val=rrad(ndat(nn),nn)*exp(sp2(ndat(nn),nn)*(rr-rrad(ndat(nn),nn))/crho(ndat(nn),nn)) 502 cgrad(1)=val*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 503 ches(1,1)=cgrad(1)*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 504 else ! INTERPOLATION 505 uu=rrad(inx+1,nn)-rrad(inx,nn) 506 uu_inv=1.0_dp/uu 507 aa=(rrad(inx+1,nn)-rr)*uu_inv 508 bb=(rr-rrad(inx,nn))*uu_inv 509 cc=(aa*aa*aa-aa)*uu*uu*0.16666666666666666_dp 510 dd=(bb*bb*bb-bb)*uu*uu*0.16666666666666666_dp 511 val=aa*crho(inx,nn)+bb*crho(inx+1,nn)+cc*sp3(inx,nn)+dd*sp3(inx+1,nn) 512 cgrad(1)=(crho(inx+1,nn)-crho(inx,nn))*uu_inv& 513 & -(3._dp*aa*aa-1._dp)*uu*0.16666666666666666_dp*sp3(inx,nn)+& 514 & (3._dp*bb*bb-1._dp)*uu*0.16666666666666666_dp*sp3(inx+1,nn) 515 ches(1,1)=aa*sp3(inx,nn)+bb*sp3(inx+1,nn) 516 517 end if ! TRANSFORMATION TO CARTEZ. COORD. 518 519 cgrad1_rr_inv=cgrad(1)*rr_inv 520 coeff=(ches(1,1)-cgrad1_rr_inv)*rr_inv*rr_inv 521 cgrad(3)=trsf(3,1)*cgrad(1) 522 cgrad(2)=trsf(2,1)*cgrad(1) 523 cgrad(1)=trsf(1,1)*cgrad(1) 524 ches(1,1)=coeff*vt1*vt1+cgrad1_rr_inv 525 ches(2,2)=coeff*vt2*vt2+cgrad1_rr_inv 526 ches(3,3)=coeff*vt3*vt3+cgrad1_rr_inv 527 ches(1,2)=coeff*vt1*vt2 ; ches(2,1)=coeff*vt1*vt2 528 ches(1,3)=coeff*vt1*vt3 ; ches(3,1)=coeff*vt1*vt3 529 ches(2,3)=coeff*vt2*vt3 ; ches(3,2)=coeff*vt2*vt3 530 531 else ! case rr==0 532 533 val=crho(1,nn)-sp2(1,nn)*rrad(1,nn)/2._dp 534 do ii=1,3 535 do ll=1,3 536 ches(ii,ll)=0._dp 537 end do 538 cgrad(ii)=0._dp 539 ches(ii,ii)=sp2(1,nn)/rrad(1,nn) 540 end do 541 542 end if ! rr>0 or rr==0 543 544 do ii=1,3 545 do ll=1,3 546 hrho(ii,ll)=hrho(ii,ll)+ches(ii,ll) 547 end do 548 grho(ii)=grho(ii)+cgrad(ii) 549 end do 550 rho=rho+val 551 552 end if ! rr2< rrad_nn*rrad_nn 553 554 if (selct==-1) then 555 if (rr2 < rminl(jj)*rminl(jj) ) then 556 iat=jj 557 ipos=kk 558 rdmin=sqrt(rr2) 559 end if 560 elseif (selct==-2) then 561 cycle 562 else 563 if (rr2 < rdmin*rdmin) then 564 iat=jj 565 ipos=kk 566 rdmin=sqrt(rr2) 567 end if 568 end if 569 570 end do 571 end do 572 573 end if 574 575 !write(51,'(":GRADt ",3F16.8)') grho 576 !write(52,'(":LAPt ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 577 !write(52,'(":HESNt ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 578 579 !if(abs(cumul_cpu-cumul_cpu_old)>0.499)then 580 !write(std_out,'(a,f7.1)' )' vgh_rho : cumul_cpu=',cumul_cpu 581 !cumul_cpu_old=cumul_cpu 582 !end if 583 584 end subroutine vgh_rho