1! 2! Copyright (C) 2004 PWSCF group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9! --------------------------------------------------------------- 10subroutine find_coefficients & 11 (ik,psi,energy,r,dx,aenorm,vpot,c,c2,lam,ndm) 12 ! ----------------------------------- 13 ! 14 ! recherche des coefficients du polynome 15 ! 16 use kinds, only : DP 17 USE random_numbers, ONLY : randy 18 implicit none 19 integer, intent(in) :: ndm, lam, ik 20 real(DP), intent(in):: vpot(ndm), psi(ndm), r(ndm), dx, energy 21 real(DP), intent(out):: c2, c(6) 22 ! 23 real(DP ):: amat(6,6), y(6), rc, aenorm 24 integer ipvt(6), i, info, n, ndiv 25 real(DP):: c2o, dc2, dcmin, newvalue, oldvalue, precision 26 real(DP), external:: funz 27 character(len=10) :: prec 28 ! 29 do i = 1,6 30 c(i) = 0.0_dp 31 end do 32 rc = r(ik) 33 34 call fill_matrix(amat,rc,lam) 35 ! 36 ! LU factorization of matrix A = amat 37 ! 38 call DGETRF(6,6,amat,6,ipvt,info) 39 ! 40 ! calculate coefficients y in linear system Ax = y 41 ! 42 call eval_coeff(r,psi,ik,lam,energy,dx,vpot,y) 43 ! 44 ! find value of c2 solving norm-conservation requirement 45 ! This is done by miniming with a random search funz**2 46 ! We are looking for the smallest solution 47 ! 48 c2o = 0.0_dp ! starting point 49 dc2 = 0.1_dp ! starting range 50 dcmin= 1.0e-10_dp ! minimum range 51 n = 0 ! counter of failed attempts 52 ndiv= 200 ! afte ndiv failed attempts reduce range 53 ! 54 precision = 7.e-10_dp! a small number 55 ! 56 oldvalue = funz(amat,ipvt,y,rc,ik,aenorm,c2o,c,c2, & 57 lam,r,dx,ndm)**2 5810 continue 59 c2 = c2o + (0.5_dp - randy())*dc2 60 newvalue = funz(amat,ipvt,y,rc,ik,aenorm,c2,c,c2,lam, & 61 r,dx,ndm)**2 62 if (newvalue < precision) return 63 if (newvalue < oldvalue) then 64 n=0 65 c2o = c2 66 oldvalue=newvalue 67 else 68 n=n+1 69 ! after ndiv failed attempts reduce the size of the interval 70 if (n > ndiv) then 71 n=0 72 dc2=dc2/10.0_dp 73 end if 74 ! if the size of the interval is too small quit 75 if (dc2 < 1.0d-12) then 76 c2=c2o 77 newvalue = funz(amat,ipvt,y,rc,ik, & 78 aenorm,c2,c,c2,lam,r,dx,ndm)**2 79 write(prec,'(e10.4)') newvalue 80 call infomsg('find_coeff','giving up minimization, '//& 81 'the error is still '//prec) 82 return 83 end if 84 end if 85 go to 10 86 87 return 88end subroutine find_coefficients 89 90! 91! --------------------------------------------------------------- 92subroutine fill_matrix(amat,rc,l) 93 ! --------------------------------------------------------------- 94 ! 95 use kinds, only : DP 96 implicit none 97 real(DP) :: amat(6,6),rc 98 integer :: l 99 ! 100 ! this routine fills the matrix with the coefficients taken 101 ! from p,p',p'',p''', p(iv), where p is 102 ! p(r)= c0 + c4 r^4 + c6 r^6 + c8 r^8 + ... 103 integer pr1(6),cr1(6),pr(6),cr(6),i,j 104 ! matrices representing the coefficients (cr) and the powers of r ( pr) 105 data pr1/0,4,6,8,10,12/ 106 data cr1/1,1,1,1,1,1/ 107 108 do i=1,6 109 pr(i) = pr1(i) 110 cr(i) = cr1(i) 111 end do 112 do i = 1,5 113 ! fill matrix row by row 114 do j = 1,6 115 amat(i,j) = cr(j) * rc**(pr(j)*1.0_dp) 116 end do 117 ! derivate polynomial expression 118 do j = 1,6 119 cr(j) = cr(j) * pr(j) 120 pr(j) = max(0,pr(j)-1) 121 end do 122 end do 123 do j = 1,6 124 amat(6,j) = 0.0_dp 125 end do 126 amat(6,2) = 2.0_dp*l +5.0_dp 127 128 return 129end subroutine fill_matrix 130! 131! ---------------------------------------------------------------- 132subroutine eval_coeff (r,psi,ik,l,energy,dx,vpot,y) 133 ! ---------------------------------------------------------------- 134 ! calcule les coefficients dependant de la fct d'onde calculee 135 ! avec tous les electrons. ces coefficients servent a la resolution 136 ! du systeme lineaire. 137 ! en entree : ik,nx,vpot comme dans le programme principal 138 ! en sortie : une matrice colonne y contenant les coefficients 139 ! 140 use kinds, only : DP 141 implicit none 142 integer :: ik,l,lp1 143 real(DP):: r(ik+3),psi(ik+3),vpot(ik+3),energy,dx,y(6) 144 real(DP) p,dpp,d,vae,dvae,ddvae,rc,rc2,rc3 145 real(DP) deriv_7pts,deriv2_7pts 146 external deriv_7pts,deriv2_7pts 147 ! 148 ! evaluer p et p' ( dpp ) 149 rc = r(ik) 150 p = psi(ik) 151 dpp = deriv_7pts(psi,ik,rc,dx) 152 if (p.lt.0.0) then 153 p = -p 154 dpp = -dpp 155 endif 156 d = dpp /p 157 ! evaluer vae, dvae, ddvae 158 vae = vpot(ik) 159 dvae = deriv_7pts(vpot,ik,rc,dx) 160 ddvae = deriv2_7pts(vpot,ik,rc,dx) 161 ! calcul de parametres intervenant dans les calculs successifs 162 lp1 = l + 1 163 rc2 = rc * rc 164 rc3 = rc2* rc 165 ! 166 y(1) = log ( p / rc**lp1 ) 167 y(2) = dpp/p - (lp1 / rc) 168 y(3) = vae - energy + (lp1*lp1)/rc2 - d*d 169 y(4) = dvae - 2.0_dp*(vae - energy + l*lp1/rc2)*d - & 170 2.0_dp*(lp1*lp1)/rc3 & 171 + 2.0_dp*(d*d*d) 172 y(5) = ddvae - 2.0_dp*(dvae - 2.0_dp*l*lp1/rc3)*d & 173 + 6.0_dp*lp1*lp1/(rc3*rc) & 174 - 2.0_dp*(vae - energy + l*lp1/rc2 - 3.0_dp*d*d)* & 175 (vae - energy + l*lp1/rc2 - d*d) 176 return 177end subroutine eval_coeff 178! ---------------------------------------------------------------- 179function deriv2_7pts(f,ik,rc,h) 180 ! ---------------------------------------------------------------- 181 ! 182 ! evaluates the second derivative of function f, the function 183 ! is given numerically on logarithmic mesh r. 184 ! nm = dimension of mesh 185 ! ik = integer : position of the point in which the derivative 186 ! will be evaluated. 187 ! h is distance between x(i) and x(i+1) where 188 ! r(i) = exp(x(i))/znesh & r(j) = exp(x(j))/znesh 189 ! 190 use kinds, only : DP 191 implicit none 192 integer :: a(7),n,nm,i,ik 193 real(DP) :: f(ik+3),rc,h,sum,sum1,deriv_7pts,deriv2_7pts 194 ! coefficients for the formula in abramowitz & stegun p.914 195 ! the formula is used for 7 points. 196 data a/4,-54,540,-980,540,-54,4/ ! these are coefficients 197 198 ! formula for linear mesh 199 sum = 0.0_dp 200 do i=1,7 201 sum = sum + a(i)*f(i-4+ik) 202 end do 203 sum = 2.0_dp*sum/(720.0_dp*h**2) 204 ! transform to logarithmic mesh 205 sum1 = deriv_7pts(f,ik,rc,h) 206 deriv2_7pts = sum/(rc*rc) - sum1 /rc 207 208 return 209end function deriv2_7pts 210! --------------------------------------------------------------- 211function deriv_7pts(f,ik,rc,h) 212 ! --------------------------------------------------------------- 213 ! evaluates the first derivative of function f, the function 214 ! is given numerically on logarithmic mesh r. 215 ! nm = dimension of mesh 216 ! ik = integer : position of the point in which the derivative 217 ! will be evaluated. 218 ! h is distance between x(i) and x(i+1) where 219 ! r(i) = exp(x(i))/znesh & r(j) = exp(x(j))/znesh 220 ! 221 use kinds, only : DP 222 implicit none 223 integer :: a(7),n,ik,i 224 real(DP) :: f(ik+3),rc,h,sum,deriv_7pts 225 ! coefficients for the formula in abramowitz & stegun p.914 226 data a/-12,108,-540,0,540,-108,12/ 227 228 ! formula for linear mesh 229 sum = 0 230 do i=1,7 231 sum = sum + a(i)*f(i-4+ik) 232 end do 233 deriv_7pts = sum/(720.0*h) 234 ! transform to logarithmic mesh 235 deriv_7pts = deriv_7pts /rc 236 237 return 238end function deriv_7pts 239 240! -------------------------------------------------------- 241function funz(amat,ipvt,y,rc,ik,aenorm,x,c,c2, & 242 lam,r,dx,ndm) 243 ! -------------------------------------------------------- 244 ! cette fonction est la fonction de c2 qu'il faut annuler 245 ! pour trouver une valeur de c2 qui verifie la conservation 246 ! de la charge de coeur. 247 ! cette fonction calcule de plus a chaque fois les ci qui 248 ! verifient les equations lineaires avec un c2 donne. 249 ! en entree : une valeur de c2 donnee dans x 250 ! en sortie, la valeur de la fonction pour cette valeur de x: 251 ! c'est la fonction qui correspond a l'equation integrale 252 ! 253 use kinds, only : DP 254 implicit none 255 integer :: ndm,lam,ipvt(6),ik,mesh,i,istart,info 256 real(DP) :: funz, x, c(6), c2, amat(6,6), y(6), rc, aenorm, & 257 r(ndm), dx, chip2, psnorm, f0, f1, f2 258 external chip2 259 260 261 ! resolution du systeme lineaire pour cette valeur de x (=c2) 262 c(1) = y(1) - x*rc**2 ! ajoute les termes en c2 263 c(2) = y(2) - 2*x*rc 264 c(3) = y(3) - 2*x 265 c(4) = y(4) ! pas de coeff en c2 266 c(5) = y(5) 267 c(6) = -x*x 268 ! call dges(amat,6,6,ipvt,c,0) ! resoud le systeme 269 call DGETRS('N',6,1,amat,6,ipvt,c,6,info) ! resoud le systeme 270 ! calcul de la norme de la pseudo-fonction d'onde 271 psnorm = 0.0_DP 272 istart= 2-mod(ik,2) 273 f2 = r(istart) * chip2(c,x,lam,r(istart)) 274 do i = istart,ik-2,2 275 f0 = f2 276 f1 = r(i+1)*chip2(c,x,lam,r(i+1)) 277 f2 = r(i+2)*chip2(c,x,lam,r(i+2)) 278 psnorm = psnorm + f0 + 4*f1 + f2 279 end do 280 psnorm = psnorm * dx/3.0_dp + r(1)**(2*lam+3)/(2*lam+3) 281 funz = log( psnorm / aenorm ) 282 ! 283 return 284end function funz 285 286! -------------------------------------------------------- 287function chip2(c,c2,l,r) 288 ! -------------------------------------------------------- 289 use kinds, only : DP 290 implicit none 291 real(DP):: chip2,c(6),c2,r,r2 292 integer :: l 293 294 r2 = r**2 295 chip2 = r2**(l+1) * exp(2.0_dp* & 296 (((((((c(6)*r2+c(5))*r2+c(4))*r2+c(3))*r2+c(2))*r2+c2)*r2)+c(1))) 297 298 return 299end function chip2 300 301! ---------------------------------------------------------------- 302function der3num(r,f,ik,nm,h) 303 ! ---------------------------------------------------------------- 304 ! calcule la 3eme derivee d'une fonction donnee numeriquement 305 ! sur un maillage logarithmique. 306 ! en entree : r maillage; f fonction ; ik : position de l'eval. 307 ! nm : dimension de f et r ; h comme dx dans 308 ! le programme principal 309 ! attention ! precision pas fantastique ! 310 use kinds, only : DP 311 implicit none 312 integer :: nm, ik,i 313 real(DP):: r(nm),f(nm), y(7), h, deriv_7pts, deriv2_7pts 314 real(DP) :: der3num 315 316 do i=1,7 317 y(i) = deriv2_7pts(f,i+ik-4,r(i+ik-4),h) 318 end do 319 der3num = deriv_7pts(y,4,r(ik),h) 320 321 return 322end function der3num 323! ---------------------------------------------------------------- 324function der4num(r,f,ik,nm,h) 325 ! ---------------------------------------------------------------- 326 ! idem que der3num, mais pour la 4eme derivee 327 ! attention ! precision pas terrible !! 328 use kinds, only : DP 329 implicit none 330 integer:: nm,i,ik 331 real(DP) :: der4num, r(nm),f(nm),y(7),h,deriv2_7pts 332 333 do i=1,7 334 y(i) = deriv2_7pts(f,i+ik-4,r(i+ik-4),h) 335 end do 336 der4num= deriv2_7pts(y,4,r(ik),h) 337 338 return 339end function der4num 340 341! ---------------------------------------------------------------- 342function der3an(l,c,c2,rc) 343 ! ---------------------------------------------------------------- 344 ! 345 ! 3rd derivative of r^(l+1) e^p(r) 346 ! 347 use kinds, only : DP 348 implicit none 349 integer :: l 350 real(DP):: c(6), c2, rc,pr,dexpr,d2expr,d3expr, der3an 351 external pr,dexpr,d2expr,d3expr 352 353 der3an= (l-1)*l*(l+1)*rc**(l-2) *exp(pr(c,c2,rc)) + & 354 3.0_DP * l*(l+1)*rc**(l-1) * dexpr(c,c2,rc) + & 355 3.0_DP * (l+1)*rc**(l ) *d2expr(c,c2,rc) + & 356 rc**(l+1) *d3expr(c,c2,rc) 357 358 return 359end function der3an 360 361! ---------------------------------------------------------------- 362function der4an(l,c,c2,rc) 363 ! ---------------------------------------------------------------- 364 ! 365 ! 4rd derivative of r^(l+1) e^p(r) 366 ! 367 use kinds, only : DP 368 implicit none 369 integer l 370 real(DP):: c(6), c2, rc,pr,dexpr,d2expr,d3expr,d4expr,der4an 371 external pr,dexpr,d2expr,d3expr,d4expr 372 373 der4an = (l-2)*(l-1)*l*(l+1)*rc**(l-3) * exp(pr(c,c2,rc)) + & 374 4.0_DP *(l-1)*l*(l+1)*rc**(l-2) * dexpr(c,c2,rc) + & 375 6.0_DP * l*(l+1)*rc**(l-1) * d2expr(c,c2,rc) + & 376 4.0_DP * (l+1)*rc**(l ) * d3expr(c,c2,rc) + & 377 rc**(l+1) * d4expr(c,c2,rc) 378 379 return 380end function der4an 381! ---------------------------------------------------------------- 382function dexpr(c,c2,rc) 383 ! ---------------------------------------------------------------- 384 ! 385 ! 1st derivative of e^p(r) 386 ! 387 use kinds, only : DP 388 implicit none 389 real(DP) :: c(6), c2, rc, pr, dpr, dexpr 390 external pr, dpr 391 392 dexpr= exp(pr(c,c2,rc)) * dpr(c,c2,rc) 393 394 return 395end function dexpr 396! ---------------------------------------------------------------- 397function d2expr(c,c2,rc) 398 ! ---------------------------------------------------------------- 399 ! 400 ! 2nd derivative of e^p(r) 401 ! 402 use kinds, only : DP 403 implicit none 404 real(DP) :: c(6), c2, rc, pr, dpr, d2pr, d2expr 405 external pr, dpr, d2pr 406 407 d2expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**2+d2pr(c,c2,rc) ) 408 409 return 410end function d2expr 411! ---------------------------------------------------------------- 412function d3expr(c,c2,rc) 413 ! ---------------------------------------------------------------- 414 ! 415 ! 3nd derivative of e^p(r) 416 ! 417 use kinds, only : DP 418 implicit none 419 real(DP):: c(6), c2, rc, pr, dpr, d2pr, d3pr, d3expr 420 external pr, dpr, d2pr, d3pr 421 422 d3expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**3 & 423 + 3*dpr(c,c2,rc)*d2pr(c,c2,rc) & 424 + d3pr(c,c2,rc) ) 425 426 return 427end function d3expr 428! ---------------------------------------------------------------- 429function d4expr(c,c2,rc) 430 ! ---------------------------------------------------------------- 431 ! 432 ! 4th derivative of e^p(r) 433 ! 434 use kinds, only : DP 435 implicit none 436 real(DP):: c(6), c2, rc, pr, dpr, d2pr, d3pr, d4pr, d4expr 437 external pr, dpr, d2pr, d3pr, d4pr 438 439 d4expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**4 & 440 + 6.0_DP* dpr(c,c2,rc)**2*d2pr(c,c2,rc) & 441 + 3.0_DP*d2pr(c,c2,rc)**2 & 442 + 4.0_DP* dpr(c,c2,rc)*d3pr(c,c2,rc) & 443 + d4pr(c,c2,rc) ) 444 445 return 446end function d4expr 447! 448! -------------------------------------------------------- 449function pr(c,c2,x) 450 ! -------------------------------------------------------- 451 ! cette fonction evalue le polynome dont les coefficients 452 ! sont dans c d'apres la forme suivante 453 ! p(x) = c(2)*x^4 + c(3) x^6 + c(4) x^8 + c(5) x^10 454 ! +c(6) x^12 + c2*x^2 + c(1) 455 use kinds, only : DP 456 implicit none 457 real(DP) :: c(6), c2, x, y, pr 458 459 y = x*x 460 pr = (((((c(6)*y+c(5))*y+c(4))*y+c(3))*y+c(2))*y+c2)*y & 461 + c(1) 462 463 return 464end function pr 465! 466function dpr(c,c2,x) 467 ! -------------------------------------------------------- 468 use kinds, only : DP 469 implicit none 470 real(DP) :: c(6),c2,x,dpr 471 472 dpr = 2*c2*x + 4*c(2)*x**3 + 6*c(3)*x**5 + 8*c(4)*x**7 & 473 + 10*c(5)*x**9 + 12*c(6)*x**11 474 475 return 476end function dpr 477! 478function d2pr(c,c2,x) 479 ! -------------------------------------------------------- 480 use kinds, only : DP 481 implicit none 482 real(DP) :: c(6),c2,x,d2pr 483 484 d2pr = 2*c2 + 12*c(2)*x**2 + 30*c(3)*x**4 + 56*c(4)*x**6 & 485 + 90*c(5)*x**8 + 132*c(6)*x**10 486 487 return 488end function d2pr 489! 490function d3pr(c,c2,x) 491 ! -------------------------------------------------------- 492 use kinds, only : DP 493 implicit none 494 real(DP) :: c(6),c2,x,d3pr 495 496 d3pr = 24*c(2)*x + 120*c(3)*x**3 + 336*c(4)*x**5 & 497 + 720*c(5)*x**7 + 1320*c(6)*x**9 498 499 return 500end function d3pr 501! 502function d4pr(c,c2,x) 503 ! -------------------------------------------------------- 504 use kinds, only : DP 505 implicit none 506 real(DP) :: c(6),c2,x,d4pr 507 508 d4pr = 24.0_dp*c(2) + 360.0_dp*c(3)*x**2 + 1680.0_dp*c(4)*x**4 & 509 + 5040.0_dp*c(5)*x**6 + 11880.0_dp*c(6)*x**8 510 511 return 512end function d4pr 513