← Back to branch summary

~albertog/siesta/psml-chlocal

« back to all changes in this revision

Viewing changes to Src/atom.F

  • Committer: Alberto Garcia
  • Date: 2017-02-02 15:17:57 UTC
  • Revision ID: albertog@icmab.es-20170202151757-wj7xlzlixl9fb82f
New heuristics. Mostly works, with tighter parameters

For Si, a good combination of parameters is:

      psml.chlocal.build.exp_al 1.012
      psml.chlocal.build.rr1  0.001
      psml.chlocal.build.r0  0.004
      psml.chlocal.build.delta 0.005
      psml.chlocal.build.poly_order_rlin 4

although there could be some improvement.

modified:
  Src/Makefile
  Src/atom.F
  Src/hamann.f90

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/atom.F
536
536
 
537
537
             if (use_psml_vlocal) then
538
538
                call get_psml_vlocal_chlocal(basp%psml_handle,is,
539
 
     $                               rofi,drdi,nrval,a,b,Zval,vlocal,
 
539
     $                               rofi,drdi,nrval,Zval,vlocal,
540
540
     $                               chlocal,nchloc)
541
541
 
542
542
                if (write_ion_plot_files) then
2439
2439
          
2440
2440
        end subroutine read_vps
2441
2441
!
2442
 
      subroutine get_psml_vlocal_chlocal(ps,is,rofi,drdi,nrval,a,b,
 
2442
      subroutine get_psml_vlocal_chlocal(ps,is,rofi,drdi,nrval,
2443
2443
     $                                   Zval,vlocal,chlocal,nchloc)
 
2444
 
 
2445
      use m_hamann, only: derivs, dpnint
2444
2446
      
2445
2447
      !! Extracts Vlocal from PSML file and computes Chlocal
2446
2448
      !! by computing its Laplacian
2452
2454
      real(dp), intent(in)          :: rofi(:)
2453
2455
      real(dp), intent(in)          :: drdi(:)
2454
2456
      integer, intent(in)           :: nrval
2455
 
      real(dp), intent(in)          :: a, b  ! r(i) = b*(exp(a(i-1))-1)
2456
2457
      real(dp), intent(in)          :: Zval
2457
2458
      real(dp), intent(out)         :: vlocal(:)
2458
2459
      real(dp), intent(out)         :: chlocal(:)
2459
2460
      integer, intent(out)          :: nchloc
2460
2461
 
2461
 
      integer  :: ir, i
 
2462
      integer  :: ir, i, mmax, n0, npts, nrange
2462
2463
      integer  :: nchloc_charge, nchloc_vlocal
 
2464
      integer  :: poly_order_rlin, poly_order_rofi
2463
2465
      character(len=40) :: vlocal_type
2464
2466
      
2465
 
      real(dp) :: qtot, r, g2, charge
2466
 
      real(dp) :: rcut_chloc
2467
 
      real(dp), dimension(:), allocatable :: vlp, vlpp, vllap
 
2467
      real(dp) :: qtot, r, g2, charge, al, rr1, delta, Z
 
2468
      real(dp) :: rcut_chloc, exp_al, r0
 
2469
      real(dp), dimension(:), allocatable :: vl, vlp, vlpp, vllap
 
2470
      real(dp), dimension(:), allocatable :: rr, ch, rlin, chlocal_rlin
2468
2471
      real(dp), parameter :: pi = 3.141592653589793_dp
2469
2472
 
2470
2473
      real(dp), parameter  :: eps_vlocal=1.0e-5_dp
2472
2475
 
2473
2476
      logical :: use_charge_cutoff
2474
2477
      
2475
 
      external :: derivs, dpnint
2476
 
 
2477
2478
      if (abs(ps_ZPseudo(ps)-Zval) > 1.e-8_dp) then
2478
2479
         call die("Zval mismatch in psml file")
2479
2480
      endif
2493
2494
 
2494
2495
!     Extraction of chlocal
2495
2496
 
2496
 
      allocate(vlp(nrval),vlpp(nrval),vllap(nrval))
 
2497
!     Set up reference grid in the manner of Don Hamann
 
2498
!     so that we can use his 'derivs' routine below
 
2499
 
 
2500
      Z = ps_AtomicNumber(ps)
 
2501
      
 
2502
      ! al = log(1.012_dp)       ! spacing factor
 
2503
      !rr1 = 0.0005_dp / Z       ! first point (r(1))
 
2504
      ! make them larger
 
2505
      rr1 = fdf_get("psml.chlocal.build.rr1",0.001_dp)
 
2506
      exp_al = fdf_get("psml.chlocal.build.exp_al",1.04_dp)
 
2507
      al = log(exp_al)       ! spacing factor
 
2508
      mmax = log(rofi(nrval)/rr1)/al ! range of points
 
2509
 
 
2510
      ! First point to keep
 
2511
      r0 = fdf_get("psml.chlocal.build.r0",0.004_dp)
 
2512
      if (r0 <= rr1) call die("psml: chlocal build: r0 < rr1")
 
2513
      
 
2514
      allocate(rr(mmax),vl(mmax))
 
2515
      do i = 1, mmax
 
2516
         rr(i) = rr1 * exp(al*(i-1))
 
2517
         vl(i) = 2.0_dp * ps_LocalPotential_Value(ps,rr(i))
 
2518
      enddo
 
2519
 
 
2520
      allocate(vlp(mmax),vlpp(mmax),vllap(mmax))
 
2521
      allocate(ch(mmax))
2497
2522
 
2498
2523
      ! Compute laplacian by finite differences
2499
 
      call derivs(nrval,vlocal,a,b,rofi,vlp,vlpp,vllap)
2500
 
      chlocal(:) = - vllap(:) / (8.0_dp*pi) ! Rydberg units
2501
 
 
 
2524
      call derivs(mmax,vl,al,0.0_dp,rr,vlp,vlpp,vllap)
 
2525
      ch(:) = - vllap(:) / (8.0_dp*pi) ! Rydberg units
 
2526
 
 
2527
      call file_out(mmax,rr,ch,
 
2528
     .     trim(global_label)//".Chlocal_rr_psml")
 
2529
 
 
2530
 
 
2531
      ! Now interpolate into linear grid, using
 
2532
      ! the points from r0 onwards
 
2533
      delta = fdf_get("psml.chlocal.build.delta",0.01_dp)
 
2534
      n0 = log(r0/rr1)/al       ! range of points
 
2535
      npts= 2 + rofi(nrval)/delta   ! overkill. will have to reduce range
 
2536
      allocate(rlin(npts))
 
2537
      do i = 1, npts
 
2538
         rlin(i) = delta*(i-1)
 
2539
      enddo
 
2540
 
 
2541
      allocate(chlocal_rlin(npts))
 
2542
      chlocal_rlin(:) = 0.0_dp
 
2543
      ! to avoid out-of-range errors in interp routine
 
2544
      nrange = log(rr(mmax)/rr1)/al - 2  
 
2545
      poly_order_rlin =
 
2546
     $     fdf_get("psml.chlocal.build.poly_order_rlin",7)
 
2547
      call dpnint(rr(n0),ch(n0),mmax-n0+1,rlin,chlocal_rlin,
 
2548
     $            nrange,poly_order=poly_order_rlin)
 
2549
      call file_out(npts,rlin,chlocal_rlin,
 
2550
     .     trim(global_label)//".Chlocal_rlin_psml")
 
2551
 
 
2552
      ! Now interpolate into Siesta's grid
 
2553
      chlocal(:) = 0.0_dp
 
2554
      poly_order_rofi =
 
2555
     $     fdf_get("psml.chlocal.build.poly_order_rofi",7)
 
2556
      call dpnint(rlin,chlocal_rlin,npts,rofi,chlocal,
 
2557
     $     nrval-20,poly_order=poly_order_rofi)
 
2558
 
 
2559
      call file_out(nrval,rofi,chlocal,
 
2560
     .     trim(global_label)//".Chlocal_psml")
 
2561
 
 
2562
      deallocate(rlin, chlocal_rlin)
 
2563
      deallocate(rr, vl, ch, vlp, vlpp, vllap)
 
2564
      
2502
2565
      ! Compute chlocal cutoff
2503
2566
      ! This section tries to keep compatibility with Siesta's behavior
2504
2567
      ! for its own vlocal/chlocal
2584
2647
 
2585
2648
        endif
2586
2649
        
2587
 
        deallocate(vlp, vlpp, vllap)
2588
 
        
2589
2650
      end subroutine get_psml_vlocal_chlocal
2590
2651
 
2591
2652
      subroutine get_psml_kbprojs(is,ps,atm_label,a,b,rofi,nrval,Zval,