1! 2! Copyright (C) 2020 Quantum ESPRESSO Foundation 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 9MODULE cpmd_module 10 ! 11 ! Read and convert a pseudopotential written in the CPMD format, TYPE: 12 ! NORMCONSERVING [ NUMERIC, LOGARITHMIC, CAR ], single radial grid 13 ! to unified pseudopotential format (v.2) 14 ! 15 USE upf_kinds, ONLY: dp 16 ! 17 IMPLICIT NONE 18 PRIVATE 19 PUBLIC :: read_cpmd 20 ! 21 CHARACTER (len=80) title 22 ! 23 ! All variables read from CPMD file format 24 ! 25 INTEGER :: ixc, pstype = 0 26 REAL(dp) :: alphaxc 27 REAL(dp) :: z, zv 28 ! Grid variables 29 INTEGER :: mesh 30 REAL(dp) :: amesh, rmax, xmin 31 REAL(dp), ALLOCATABLE :: r(:) 32 ! PP variables 33 INTEGER, PARAMETER :: lmaxx=3 34 INTEGER ::lmax, nwfc=0 35 ! Car PP variables 36 REAL(dp) :: alphaloc, alpha(0:lmaxx), a(0:lmaxx), b(0:lmaxx) 37 ! Goedecker PP variables 38 INTEGER, PARAMETER :: ncmax=4, nlmax=3 39 INTEGER :: nc, nl(0:lmaxx) 40 REAL(dp) :: rc, rl(0:lmaxx), c(ncmax), h(0:lmaxx, nlmax*(nlmax+1)/2 ) 41 ! Numeric PP variables 42 REAL(dp), ALLOCATABLE :: vnl(:,:) 43 REAL(dp), ALLOCATABLE :: chi(:,:) 44 ! Core correction variables 45 LOGICAL :: nlcc=.false. 46 REAL(dp), ALLOCATABLE :: rho_atc(:) 47 ! Variables used for reading and for checks 48 INTEGER :: maxinfo, info_lines 49 PARAMETER (maxinfo = 100) 50 CHARACTER (len=80), ALLOCATABLE :: info_sect(:) 51 !------------------------------ 52 ! 53CONTAINS 54 ! 55 ! ---------------------------------------------------------- 56SUBROUTINE read_cpmd(iunps, upf) 57 ! ---------------------------------------------------------- 58 ! 59 USE pseudo_types, ONLY : pseudo_upf 60 USE upf_utils, ONLY : matches 61 ! 62 IMPLICIT NONE 63 INTEGER, INTENT(in) :: iunps 64 TYPE (pseudo_upf) :: upf 65 ! 66 INTEGER :: found = 0, closed = 0, unknown = 0 67 INTEGER :: i, l, dum, ios 68 CHARACTER (len=256) line 69 CHARACTER (len=4) token 70 REAL (dp) :: amesh_, vnl0(0:3) 71 LOGICAL :: grid_read = .false., wfc_read=.false. 72 REAL(dp), EXTERNAL :: upf_erf 73 ! 74 info_lines = 0 7510 READ (iunps,'(A)',end=20,err=20) line 76 IF (matches ("&ATOM", trim(line)) ) THEN 77 found = found + 1 78 ! Z 79 READ (iunps,'(a)',end=200,err=200) line 80 l = len_trim(line) 81 i = INDEX(line,'=') 82 READ (line(i+1:l),*) z 83 ! Zv 84 READ (iunps,'(a)',end=200,err=200) line 85 l = len_trim(line) 86 i = INDEX(line,'=') 87 READ (line(i+1:l),*) zv 88 ! XC 89 READ (iunps,'(a)',end=200,err=200) line 90 l = len_trim(line) 91 i = INDEX(line,'=') 92 READ (line(i+1:l),*) ixc, alphaxc 93 ! TYPE 94 READ (iunps,'(a)',end=200,err=200) line 95 IF ( matches("NORMCONSERVING",line) ) THEN 96 IF ( matches("NUMERIC",line) .or. matches("LOGARITHMIC",line) ) THEN 97 pstype = 1 98 ELSEIF ( matches("CAR",line) ) THEN 99 pstype = 2 100 ELSEIF ( matches("GOEDECKER",line) ) THEN 101 pstype = 3 102 ENDIF 103 ENDIF 104 IF (pstype == 0 ) CALL upf_error('read_cpmd','unknown type: '//line,1) 105 ELSEIF (matches ("&INFO", trim(line)) ) THEN 106 found = found + 1 107 ! read (iunps,'(a)') title 108 ! store info section for later perusal 109 ALLOCATE (info_sect(maxinfo)) 110 DO i=1,maxinfo 111 READ (iunps,'(a)',end=20,err=20) title 112 IF (matches ("&END", trim(title)) ) THEN 113 closed = closed + 1 114 GOTO 10 115 ELSE 116 info_sect(i) = trim(title) 117 info_lines = i 118 ENDIF 119 ENDDO 120 ELSEIF (matches ("&POTENTIAL", trim(line)) ) THEN 121 found = found + 1 122 READ (iunps,'(a)') line 123 IF ( pstype == 1 ) THEN 124 ! 125 ! NORMCONSERVING NUMERIC 126 ! 127 READ (line,*,iostat=ios) mesh, amesh_ 128 IF ( ios /= 0) THEN 129 READ (line,*,iostat=ios) mesh 130 amesh_ = -1.0d0 131 ENDIF 132 IF ( .not. grid_read ) ALLOCATE (r(mesh)) 133 ! 134 ! determine the number of angular momenta 135 ! 136 READ (iunps, '(a)') line 137 ios = 1 138 lmax= 4 139 DO WHILE (ios /= 0) 140 lmax = lmax - 1 141 READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,lmax) 142 ENDDO 143 ALLOCATE (vnl(mesh,0:lmax)) 144 vnl(1,0:lmax) = vnl0(0:lmax) 145 DO i=2,mesh 146 READ(iunps, *) r(i),(vnl(i,l),l=0,lmax) 147 ENDDO 148 IF ( .not.grid_read ) THEN 149 CALL check_radial_grid ( amesh_, mesh, r, amesh ) 150 grid_read = .true. 151 ENDIF 152 ELSEIF ( pstype == 2 ) THEN 153 ! 154 ! NORMCONSERVING CAR 155 ! 156 READ(iunps, *) alphaloc 157 ! convert r_c's written in file to alpha's: alpha = 1/r_c^2 158 alphaloc = 1.d0/alphaloc**2 159 DO lmax=-1,2 160 READ(iunps, '(A)') line 161 IF (matches ("&END", trim(line)) ) THEN 162 closed = closed + 1 163 exit 164 ENDIF 165 READ(line, *) alpha(lmax+1), a(lmax+1), b(lmax+1) 166 alpha(lmax+1) = 1.d0/alpha(lmax+1)**2 167 ENDDO 168 ELSEIF ( pstype == 3 ) THEN 169 ! 170 ! NORMCONSERVING GOEDECKER 171 ! 172 c(:) = 0.d0 173 rl(:) = 0.d0 174 nl(:) = 0 175 h(:,:) = 0.d0 176 READ(iunps, *) lmax 177 lmax = lmax - 1 178 IF ( lmax > lmaxx ) & 179 CALL upf_error('read_cpmd','incorrect parameter read',1) 180 READ(iunps, *) rc 181 READ(iunps, '(A)') line 182 READ(line, *) nc 183 IF ( nc > ncmax ) & 184 CALL upf_error('read_cpmd','incorrect parameter read',2) 185 ! I am not sure if it is possible to use nc in the same line 186 ! where it is read. Just in case, better to read twice 187 READ(line, *) dum, (c(i), i=1,nc) 188 DO l=0,lmax+1 189 READ(iunps, '(A)') line 190 IF ( matches ("&END", trim(line)) ) THEN 191 closed = closed + 1 192 exit 193 ENDIF 194 READ(line, *) rl(l), nl(l) 195 IF ( nl(l) > nlmax ) & 196 CALL upf_error('read_cpmd','incorrect parameter read',3) 197 IF ( nl(l) > 0 ) & 198 READ(line, *) rl(l), dum, ( h(l,i), i=1,nl(l)*(nl(l)+1)/2) 199 ENDDO 200 ENDIF 201 ELSEIF (matches ("&WAVEFUNCTION", trim(line)) ) THEN 202 wfc_read=.true. 203 found = found + 1 204 READ (iunps,'(A)') line 205 READ (line,*,iostat=ios) mesh, amesh_ 206 IF ( ios /= 0) THEN 207 READ (line,*,iostat=ios) mesh 208 amesh_ = -1.0d0 209 ENDIF 210 IF ( .not. grid_read ) ALLOCATE(r(mesh)) 211 ! find number of atomic wavefunctions 212 READ (iunps,'(A)') line 213 DO nwfc = lmax+1,1,-1 214 READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,nwfc-1) 215 IF ( ios == 0 ) exit 216 ENDDO 217 IF ( ios /= 0 ) & 218 CALL upf_error('read_cpmd','at least one atomic wvfct should be present',1) 219 ALLOCATE(chi(mesh,nwfc)) 220 chi(1,1:nwfc) = vnl0(0:nwfc-1) 221 DO i=2,mesh 222 READ(iunps, *) r(i),(chi(i,l),l=1,nwfc) 223 ENDDO 224 IF ( .not.grid_read ) THEN 225 CALL check_radial_grid ( amesh_, mesh, r, amesh ) 226 grid_read = .true. 227 ENDIF 228 ELSEIF (matches ("&NLCC", trim(line)) ) THEN 229 found = found + 1 230 READ (iunps, '(a)') line 231 READ(iunps, *) mesh 232 nlcc = ( mesh > 0 ) 233 IF (nlcc) THEN 234 IF ( .not. matches ("NUMERIC", trim(line)) ) & 235 CALL upf_error('read_cpmd',' only NUMERIC core-correction supported',1) 236 ALLOCATE (rho_atc(mesh)) 237 READ(iunps, * ) (r(i), rho_atc(i), i=1,mesh) 238 ENDIF 239 ELSEIF (matches ("&ATDENS", trim(line)) ) THEN 240 ! skip over &ATDENS section, add others here, if there are more. 241 DO WHILE(.not. matches("&END", trim(line))) 242 READ (iunps,'(a)') line 243 ENDDO 244 ELSEIF (matches ("&END", trim(line)) ) THEN 245 closed = closed + 1 246 ELSE 247 PRINT*, 'line ignored: ', line 248 unknown = unknown + 1 249 ENDIF 250 GOTO 10 251 25220 CONTINUE 253 254 IF ( pstype /= 3 ) THEN 255 IF (nlcc .and. found /= 5 .or. .not.nlcc .and. found /= 4) & 256 CALL upf_error('read_cpmd','some &FIELD card missing',found) 257 ELSE 258 IF (found /= 3) & 259 CALL upf_error('read_cpmd','some &FIELD card missing',found) 260 ENDIF 261 IF (closed /= found) & 262 CALL upf_error('read_cpmd','some &END card missing',closed) 263 IF (unknown /= 0 ) PRINT '("WARNING: ",i3," cards not read")', unknown 264 ! 265 IF ( .not. grid_read ) THEN 266 xmin = -7.0d0 267 amesh=0.0125d0 268 rmax =100.0d0 269 PRINT '("A radial grid must be provided. We use the following one:")' 270 PRINT '("r_i = e^{xmin+(i-1)*dx}/Z, i=1,mesh, with parameters:")' 271 PRINT '("Z=",f6.2,", xmin=",f6.2," dx=",f8.4," rmax=",f6.1)', & 272 z, xmin, amesh, rmax 273 mesh = 1 + (log(z*rmax)-xmin)/amesh 274 mesh = (mesh/2)*2+1 ! mesh is odd (for historical reasons?) 275 ALLOCATE (r(mesh)) 276 DO i=1, mesh 277 r(i) = exp (xmin+(i-1)*amesh)/z 278 ENDDO 279 PRINT '(I4," grid points, rmax=",f8.4)', mesh, r(mesh) 280 grid_read = .true. 281 ENDIF 282 rmax = r(mesh) 283 xmin = log(z*r(1)) 284 ! 285 IF ( .not. wfc_read ) PRINT '("Notice: atomic wfcs not found")' 286 ! 287 IF ( pstype == 2 ) THEN 288 ALLOCATE (vnl(mesh,0:lmax)) 289 DO l=0, lmax 290 DO i=1, mesh 291 vnl(i,l) = ( a(l) + b(l)*r(i)**2 ) * exp (-alpha(l)*r(i)**2) - & 292 zv * upf_erf (sqrt(alphaloc)*r(i))/r(i) 293 ENDDO 294 ENDDO 295 ENDIF 296 297 CALL convert_cpmd(upf) 298 RETURN 299200 CALL upf_error('read_cpmd','error in reading file',1) 300 301END SUBROUTINE read_cpmd 302 303! ---------------------------------------------------------- 304SUBROUTINE convert_cpmd(upf) 305 ! ---------------------------------------------------------- 306 ! 307 USE pseudo_types, ONLY : pseudo_upf 308 USE upf_const, ONLY : e2 309 ! 310 IMPLICIT NONE 311 ! 312 TYPE(pseudo_upf) :: upf 313 ! 314 REAL(dp), ALLOCATABLE :: aux(:) 315 REAL(dp) :: x, x2, vll, rcloc, fac 316 REAL(dp), EXTERNAL :: upf_erf 317 CHARACTER (len=20):: dft 318 CHARACTER (len=2):: label 319 CHARACTER (len=1):: spdf(0:3) = ['S','P','D','F'] 320 CHARACTER (len=2), EXTERNAL :: atom_name 321 INTEGER :: lloc, my_lmax, l, i, j, ij, ir, iv, jv 322 ! 323 ! NOTE: many CPMD pseudopotentials created with the 'Hamann' code 324 ! from Juerg Hutter's homepage have additional (bogus) entries for 325 ! pseudo-potential and wavefunction. In the 'report' they have 326 ! the same rc and energy eigenvalue than the previous angular momentum. 327 ! we need to be able to ignore that part or the resulting UPF file 328 ! will be useless. so we first print the info section and ask 329 ! for the LMAX to really use. AK 2005/03/30. 330 ! 331 DO i=1,info_lines 332 PRINT '(A)', info_sect(i) 333 ENDDO 334 IF ( pstype == 3 ) THEN 335 ! not actually used, except by write_upf to write a meaningful message 336 lloc = -3 337 rcloc=0.0 338 ELSE 339 WRITE(*,'("max L to use ( <= ",I1," ) > ")', advance="NO") lmax 340 READ (5,*) my_lmax 341 IF ((my_lmax <= lmax) .and. (my_lmax >= 0)) lmax = my_lmax 342 WRITE(*,'("local L ( <= ",I1," ), Rc for local pot (au) > ")', advance="NO") lmax 343 READ (5,*) lloc, rcloc 344 ENDIF 345 ! 346 IF ( pstype == 3 ) THEN 347 upf%generated= "Generated in analytical, separable form" 348 upf%author = "Goedecker/Hartwigsen/Hutter/Teter" 349 upf%date = "Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)" 350 ELSE 351 upf%generated= "Generated using unknown code" 352 upf%author = "unknown" 353 upf%date = "unknown" 354 ENDIF 355 upf%nv = "2.0.1" 356 upf%comment = "Info: automatically converted from CPMD format" 357 upf%psd = atom_name ( nint(z) ) 358 ! reasonable assumption 359 IF (z > 18) THEN 360 upf%rel = 'no' 361 ELSE 362 upf%rel = 'scalar' 363 ENDIF 364 IF ( pstype == 3 ) THEN 365 upf%typ = 'NC' 366 ELSE 367 upf%typ = 'SL' 368 ENDIF 369 upf%tvanp = .false. 370 upf%tpawp = .false. 371 upf%tcoulombp=.false. 372 upf%nlcc = nlcc 373 ! 374 IF (ixc==900) THEN 375 PRINT '("Pade approx. not implemented! assuming Perdew-Zunger LDA")' 376 upf%dft='SLA-PZ-NOGX-NOGC' 377 ELSEIF (ixc==1100) THEN 378 upf%dft='SLA-PZ-NOGX-NOGC' 379 ELSEIF (ixc==1111) THEN 380 upf%dft='SLA-PZ-B86-P88' 381 ELSEIF (ixc==1134 .or. ixc==1434) THEN 382 upf%dft='SLA-PW-PBX-PBC' 383 ELSEIF (ixc==1134) THEN 384 upf%dft='revPBE' 385 ELSEIF (ixc==1197) THEN 386 upf%dft='PBESOL' 387 ELSEIF (ixc==1312) THEN 388 upf%dft='BLYP' 389 ELSEIF (ixc==362) THEN 390 upf%dft='OLYP' 391 ELSEIF (ixc==1372) THEN 392 upf%dft='XLYP' 393 ELSEIF (ixc==55) THEN 394 upf%dft='HCTH' 395 ELSE 396 WRITE(*,'("Unknown DFT ixc=",i4,". Please provide a DFT name > ")', advance="NO") ixc 397 READ *, upf%dft 398 ENDIF 399 PRINT '("Assuming DFT: ",A," . Please check this is what you want")', & 400 trim(upf%dft) 401 ! 402 upf%zp = zv 403 upf%etotps =0.0d0 404 upf%ecutrho=0.0d0 405 upf%ecutwfc=0.0d0 406 IF ( lmax == lloc) THEN 407 upf%lmax = lmax-1 408 ELSE 409 upf%lmax = lmax 410 ENDIF 411 upf%lloc = lloc 412 upf%lmax_rho = 0 413 upf%nwfc = nwfc 414 ! 415 ALLOCATE( upf%els(upf%nwfc) ) 416 ALLOCATE( upf%oc(upf%nwfc) ) 417 ALLOCATE( upf%epseu(upf%nwfc) ) 418 ALLOCATE( upf%lchi(upf%nwfc) ) 419 ALLOCATE( upf%nchi(upf%nwfc) ) 420 ALLOCATE( upf%rcut_chi (upf%nwfc) ) 421 ALLOCATE( upf%rcutus_chi(upf%nwfc) ) 422 423 DO i=1, upf%nwfc 42410 WRITE(*,'("Wavefunction # ",i1,": label (e.g. 4s), occupancy > ")', advance="NO") i 425 READ (5,*) label, upf%oc(i) 426 READ (label(1:1),*, err=10) l 427 upf%els(i) = label 428 upf%nchi(i) = l 429 IF ( label(2:2) == 's' .or. label(2:2) == 'S') THEN 430 l=0 431 ELSEIF ( label(2:2) == 'p' .or. label(2:2) == 'P') THEN 432 l=1 433 ELSEIF ( label(2:2) == 'd' .or. label(2:2) == 'D') THEN 434 l=2 435 ELSEIF ( label(2:2) == 'f' .or. label(2:2) == 'F') THEN 436 l=3 437 ELSE 438 l=i-1 439 ENDIF 440 upf%lchi(i) = l 441 upf%rcut_chi(i) = 0.0d0 442 upf%rcutus_chi(i)= 0.0d0 443 upf%epseu(i) = 0.0d0 444 ENDDO 445 446 upf%mesh = mesh 447 upf%dx = amesh 448 upf%rmax = rmax 449 upf%xmin = xmin 450 upf%zmesh= z 451 ALLOCATE(upf%rab(upf%mesh)) 452 ALLOCATE(upf%r(upf%mesh)) 453 upf%r(:) = r(1:upf%mesh) 454 upf%rab(:)=upf%r(:)*amesh 455 456 ALLOCATE (upf%rho_atc(upf%mesh)) 457 IF (upf%nlcc) upf%rho_atc(:) = rho_atc(1:upf%mesh) 458 459 upf%rcloc = rcloc 460 ALLOCATE (upf%vloc(upf%mesh)) 461 ! 462 ! the factor e2=2 converts from Hartree to Rydberg 463 ! 464 IF ( upf%typ == "SL" ) THEN 465 upf%vloc(:) = vnl(1:upf%mesh,upf%lloc)*e2 466 ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,1)) 467 upf%vnl(:,:,1) = vnl(1:upf%mesh,0:upf%lmax)*e2 468 upf%nbeta= lmax 469 ELSE 470 DO i=1,upf%mesh 471 x = upf%r(i)/rc 472 x2=x**2 473 upf%vloc(i) = e2 * ( -upf%zp*upf_erf(x/sqrt(2.d0))/upf%r(i) + & 474 exp ( -0.5d0*x2 ) * (c(1) + x2*( c(2) + x2*( c(3) + x2*c(4) ) ) ) ) 475 ENDDO 476 upf%nbeta=0 477 DO l=0,upf%lmax 478 upf%nbeta = upf%nbeta + nl(l) 479 ENDDO 480 ENDIF 481 482 IF (upf%nbeta > 0) THEN 483 ALLOCATE(upf%els_beta(upf%nbeta) ) 484 ALLOCATE(upf%lll(upf%nbeta)) 485 ALLOCATE(upf%kbeta(upf%nbeta)) 486 IF ( pstype == 3 ) THEN 487 upf%kbeta(:) = upf%mesh 488 ELSE 489 iv=0 ! counter on beta functions 490 DO i=1,upf%nwfc 491 l=upf%lchi(i) 492 IF (l/=lloc) THEN 493 iv=iv+1 494 upf%kbeta(iv)=upf%mesh 495 DO ir = upf%mesh,1,-1 496 IF ( abs ( vnl(ir,l) - vnl(ir,lloc) ) > 1.0E-6 ) THEN 497 ! include points up to the last with nonzero value 498 upf%kbeta(iv)=ir+1 499 exit 500 ENDIF 501 ENDDO 502 ENDIF 503 ENDDO 504 ! the number of points used in the evaluation of integrals 505 ! should be even (for simpson integration) 506 DO i=1,upf%nbeta 507 IF ( mod (upf%kbeta(i),2) == 0 ) upf%kbeta(i)=upf%kbeta(i)+1 508 upf%kbeta(i)=MIN(upf%mesh,upf%kbeta(i)) 509 ENDDO 510 upf%kkbeta = maxval(upf%kbeta(:)) 511 ENDIF 512 ALLOCATE(upf%beta(upf%mesh,upf%nbeta)) 513 ALLOCATE(upf%dion(upf%nbeta,upf%nbeta)) 514 upf%beta(:,:) =0.d0 515 upf%dion(:,:) =0.d0 516 ALLOCATE(upf%rcut (upf%nbeta)) 517 ALLOCATE(upf%rcutus(upf%nbeta)) 518 519 IF ( pstype == 3 ) THEN 520 iv=0 ! counter on beta functions 521 DO l=0,upf%lmax 522 ij = 0 523 DO i=1, nl(l) 524 iv = iv+1 525 upf%lll(iv)=l 526 WRITE (upf%els_beta(iv), '(I1,A1)' ) i, spdf(l) 527 DO j=i, nl(l) 528 jv = iv+j-i 529 ij=ij+1 530 upf%dion(iv,jv) = h(l,ij)/e2 531 IF ( j > i ) upf%dion(jv,iv) = upf%dion(iv,jv) 532 ENDDO 533 fac= sqrt(2d0*rl(l)) / ( rl(l)**(l+2*i) * sqrt(mygamma(l+2*i)) ) 534 DO ir=1,upf%mesh 535 x2 = (upf%r(ir)/rl(l))**2 536 upf%beta(ir,iv) = upf%r(ir)**(l+2*(i-1)) * & 537 exp ( -0.5d0*x2 ) * fac * e2 538 ! ...remember: the beta functions in the UPF format 539 ! ...have to be multiplied by a factor r !!! 540 upf%beta(ir,iv) = upf%beta(ir,iv)*upf%r(ir) 541 ! 542 ENDDO 543 ! look for index kbeta such that v(i)=0 if i>kbeta 544 DO ir=upf%mesh,1,-1 545 IF ( abs(upf%beta(ir,iv)) > 1.D-12 ) exit 546 ENDDO 547 IF ( ir < 2 ) THEN 548 CALL upf_error('cpmd2upf','zero beta function?!?',iv) 549 ELSEIF ( mod(ir,2) /= 0 ) THEN 550 ! even index 551 upf%kbeta(iv) = ir 552 ELSEIF ( ir < upf%mesh .and. mod(ir,2) == 0 ) THEN 553 ! odd index 554 upf%kbeta(iv) = ir+1 555 ELSE 556 upf%kbeta(iv) = upf%mesh 557 ENDIF 558 ! not really the same thing as rc in PP generation 559 upf%rcut (iv) = upf%r(upf%kbeta(iv)) 560 upf%rcutus(iv) = 0.0 561 ENDDO 562 ENDDO 563 upf%kkbeta = maxval(upf%kbeta(:)) 564 ELSE 565 ALLOCATE(aux(upf%kkbeta)) 566 iv=0 ! counter on beta functions 567 DO i=1,upf%nwfc 568 l=upf%lchi(i) 569 IF (l/=lloc) THEN 570 iv=iv+1 571 upf%lll(iv)=l 572 upf%els_beta(iv)=upf%els(i) 573 DO ir=1,upf%kbeta(iv) 574 ! the factor e2 converts from Hartree to Rydberg 575 upf%beta(ir,iv) = e2 * chi(ir,l+1) * & 576 ( vnl(ir,l) - vnl(ir,lloc) ) 577 aux(ir) = chi(ir,l+1) * upf%beta(ir,iv) 578 ENDDO 579 upf%rcut (iv) = upf%r(upf%kbeta(iv)) 580 upf%rcutus(iv) = 0.0 581 CALL simpson(upf%kbeta(iv),aux,upf%rab,vll) 582 upf%dion(iv,iv) = 1.0d0/vll 583 ENDIF 584 ENDDO 585 DEALLOCATE(aux) 586 ENDIF 587 ELSE 588 ! prevents funny errors when writing file 589 ALLOCATE(upf%dion(upf%nbeta,upf%nbeta)) 590 ENDIF 591 592 ALLOCATE (upf%chi(upf%mesh,upf%nwfc)) 593 upf%chi(:,:) = chi(1:upf%mesh,1:upf%nwfc) 594 595 ALLOCATE (upf%rho_at(upf%mesh)) 596 upf%rho_at(:) = 0.d0 597 DO i=1,upf%nwfc 598 upf%rho_at(:) = upf%rho_at(:) + upf%oc(i) * upf%chi(:,i) ** 2 599 ENDDO 600 601 ! ---------------------------------------------------------- 602 WRITE (6,'(a)') 'Pseudopotential successfully converted' 603 ! ---------------------------------------------------------- 604 605 RETURN 606END SUBROUTINE convert_cpmd 607! 608! ------------------------------------------------------------------ 609SUBROUTINE check_radial_grid ( amesh_, mesh, r, amesh ) 610! ------------------------------------------------------------------ 611! 612 IMPLICIT NONE 613 INTEGER, INTENT (in) :: mesh 614 REAL(dp), INTENT (in) :: amesh_, r(mesh) 615 REAL(dp), INTENT (out) :: amesh 616 INTEGER :: i 617 ! 618 ! get amesh if not available directly, check its value otherwise 619 PRINT "('Radial grid r(i) has ',i4,' points')", mesh 620 PRINT "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')" 621 IF (amesh_ < 0.0d0) THEN 622 amesh = log (r(mesh)/r(1))/(mesh-1) 623 PRINT "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh 624 ELSE 625 ! not clear whether the value of amesh read from file 626 ! matches the above definition, or if it is exp(amesh) ... 627 amesh = log (r(mesh)/r(1))/(mesh-1) 628 IF ( abs ( amesh - amesh_ ) > 1.0d-5 ) THEN 629 IF ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) THEN 630 amesh = log(amesh_) 631 PRINT "('amesh = log (value read from file) = ',f10.6)",amesh 632 ELSE 633 CALL upf_error ('cpmd2upf', 'unknown real-space grid',2) 634 ENDIF 635 ELSE 636 amesh = amesh_ 637 PRINT "('amesh = value read from file = ',f10.6)",amesh 638 ENDIF 639 ENDIF 640 ! check if the grid is what we expect 641 DO i=2,mesh 642 IF ( abs(r(i) - exp((i-1)*amesh)*r(1)) > 1.0d-5) THEN 643 PRINT "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",& 644 i, r(i), exp((i-1)*amesh)*r(1) 645 CALL upf_error ('cpmd2upf', 'unknown real-space grid',1) 646 ENDIF 647 ENDDO 648 RETURN 649 END SUBROUTINE check_radial_grid 650! ------------------------------------------------------------------ 651FUNCTION mygamma ( n ) 652 !------------------------------------------------------------------ 653 ! 654 ! mygamma(n) = \Gamma(n-1/2) = sqrt(pi)*(2n-3)!!/2**(n-1) 655 ! 656 USE upf_const, ONLY : pi 657 IMPLICIT NONE 658 REAL(dp) :: mygamma 659 INTEGER, INTENT(in) :: n 660 INTEGER :: i 661 ! 662 IF ( n < 2 ) CALL upf_error('mygamma','unexpected input argument',1) 663 mygamma = sqrt(pi) / 2.d0**(n-1) 664 ! next factor is (2n-3)!! 665 do i = 2*n-3, 1, -2 666 mygamma = i*mygamma 667 end do 668 ! 669END FUNCTION mygamma 670 671END MODULE cpmd_module 672