Debugging

test/fortran_main.f90

@@ -6,15 +6,15 @@
 ! unit for rprimd is in Bohr, and for gprimd is in Bohr^-1
 ! 3. gmet : reciprocal space metric (bohr^-2)
 ! 4. ucvol : volume of unit cell (Bohr^3)
-! 5. ngfft, ngfftdg : dimension of FFT grids of the corase and fine grids
+! 5. ngfft, ngfftdg : dimension of FFT grids of the coarse and fine grids
 ! 6. natom, ntypat, typat: #. atoms, #. element types
 ! and typat records the type of each atom
 ! 7. xred : coordinate of each atom, in terms of rprimd (namely, direct coordinate)
 ! 8. filename_list : filename of the PAW xml files for each element
 
 subroutine fortran_main(ecut,ecutpaw,gmet,rprimd,gprimd,ucvol, &
     ngfft,ngfftdg,natom,ntypat,typat,xred,filename_list)
-    use libpaw_mod, only : cplex, pawrhoij
+    use libpaw_mod, only : cplex, pawrhoij, paw_ij
     implicit none
 
     real*8  :: ecut, ecutpaw !a coarse grid, and a fine grid for PAW
@@ -28,18 +28,27 @@ subroutine fortran_main(ecut,ecutpaw,gmet,rprimd,gprimd,ucvol, &
     integer :: nfft
     integer :: nspden,nsppol
     integer :: nrhoijsel,size_rhoij
+    integer :: i,j,k,l,idx
     integer, allocatable :: rhoijselect(:)
     real*8,  allocatable :: rhoijp(:,:)
     real*8,  allocatable :: vtrial(:,:), vxc(:,:)
+    real*8,  allocatable :: vtrial_row_major(:,:), vxc_row_major(:,:)
     real*8,  allocatable :: vloc(:), ncoret(:)
+    real*8,  allocatable :: vloc_col_major(:), ncoret_col_major(:)
+    real*8,  allocatable :: vloc_row_major(:), ncoret_row_major(:)
     real*8,  allocatable :: nhat(:,:), nhatgr(:,:,:)
+    real*8,  allocatable :: nhat_row_major(:,:)
+    real*8,  allocatable :: dij(:,:)
     !integer, allocatable :: typat(:)
     !real*8,  allocatable :: xred(:,:)
 
     integer :: typat(natom)
     real*8  :: xred(3,natom)
 
     character(len=264) :: filename_list(ntypat)
+    real*8  :: epsatm(ntypat)
+    real*8  :: epawdc
+    real*8  :: diff_row_col_majors
 
     write(*,*) '1. Setting up libpaw'
 
@@ -51,14 +60,29 @@ subroutine fortran_main(ecut,ecutpaw,gmet,rprimd,gprimd,ucvol, &
 
     write(*,*) 'Test code for libpaw'
     call prepare_libpaw(ecut,ecutpaw,gmet,rprimd,gprimd,ucvol,ngfft,ngfftdg, &
-        natom,ntypat,typat,xred,ixc,xclevel,filename_list,nspden,nsppol)
+        natom,ntypat,typat,xred,ixc,xclevel,filename_list,nspden,nsppol,epsatm)
 
     nfft = ngfftdg(1) * ngfftdg(2) * ngfftdg(3)
     allocate(ncoret(nfft),vloc(nfft))
+    allocate(ncoret_row_major(nfft),vloc_row_major(nfft))
+    allocate(ncoret_col_major(nfft),vloc_col_major(nfft))
     allocate(vtrial(cplex*nfft,nspden),vxc(cplex*nfft,nspden))
+    allocate(vtrial_row_major(cplex*nfft,nspden),vxc_row_major(cplex*nfft,nspden))
     allocate(nhat(nfft,nspden), nhatgr(nfft,nspden,3))
+    allocate(nhat_row_major(nfft,nspden))
 
     call get_vloc_ncoret(ngfftdg,nfft,natom,ntypat,rprimd,gprimd,gmet,ucvol,xred,ncoret,vloc)
+
+    call col_to_row_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,vloc,vloc_row_major)
+    call col_to_row_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,ncoret,ncoret_row_major)
+
+    open(unit=34,file='vloc_calc.dat')
+    open(unit=35,file='ncoret_calc.dat')
+    do i=1,nfft
+        write(34,*) vloc_row_major(i)
+        write(35,*) ncoret_row_major(i)
+    enddo
+    close(34); close(35)
 
     open(unit=10,file='rhoij')
     do iatom = 1, natom
@@ -74,12 +98,82 @@ subroutine fortran_main(ecut,ecutpaw,gmet,rprimd,gprimd,ucvol, &
     close(10)
 
     call get_nhat(natom,ntypat,xred,ngfft,nfft,nspden,gprimd,rprimd,ucvol,nhat,nhatgr)
+
+    call col_to_row_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,nhat,nhat_row_major)
 
+    open(unit=10,file='nhat_calc.dat')
+    do i=1,nfft
+        write(10,*)  nhat_row_major(i,1)
+    enddo
+    close(10)
     open(unit=10,file='veff')
-    read(10,*) vtrial
-    read(10,*) vxc
+    read(10,*) vtrial_row_major
+    read(10,*) vxc_row_major
     close(10)
 
-    call calculate_dij(natom,ntypat,ixc,xclevel,nfft,nspden,xred,ucvol,gprimd,vtrial,vxc)
+
+    !call row_to_col_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,vloc_row_major,vloc_col_major)
+    !call row_to_col_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,ncoret_row_major,ncoret_col_major)
+    !write(6,*) "test 1 vloc row,col", sum(abs(vloc_col_major-vloc))
+    !write(6,*) "test 2 ncoret row,col", sum(abs(ncoret_col_major-ncoret))
+    !call row_to_col_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,vtrial_row_major,vtrial)
+    !call row_to_col_major_ordering(ngfftdg(1),ngfftdg(2),ngfftdg(3),nfft,vxc_row_major,vxc)
+
+
+    vtrial = vtrial / 2d0
+    vxc = vxc / 2d0
+
+    call calculate_dij(natom,ntypat,ixc,xclevel,nfft,nspden,xred,ucvol,gprimd,vtrial,vxc,epawdc)
+
+    open(unit=15,file='dij_calc')
+    do iatom = 1, natom
+        size_rhoij = size(paw_ij(iatom)%dij)
+        allocate(dij(size_rhoij,nspden))
+        call get_dij(iatom,size_rhoij,nspden,dij)
+        do i=1,size_rhoij
+            write(15,*) dij(i,1)
+        enddo
+        deallocate(dij)
+    enddo
+
+end subroutine
+
+subroutine col_to_row_major_ordering(nx,ny,nz,nn,v_col,v_row)
+    implicit none
+    integer :: nx,ny,nz,nn
+    real*8  :: v_row(nn),v_col(nn)
+    integer :: i,j,k
+    integer :: idx, l
+    if(nx*ny*nz /= nn) stop "nx*ny*nz /= nn in col_to_row"
+    idx=1
+    do i = 1, nx
+        do j = 1, ny
+            do k = 1, nz
+                l = (i - 1) + (j - 1) * nx + (k - 1) * nx * ny + 1
+                v_row(idx) = v_col(l)
+                idx=idx+1
+            end do
+        end do
+    end do
 
 end subroutine
+
+subroutine row_to_col_major_ordering(nx,ny,nz,nn,v_row,v_col)
+    implicit none
+    integer :: nx,ny,nz,nn
+    real*8  :: v_row(nn),v_col(nn)
+    integer :: i,j,k
+    integer :: idx, l
+    if(nx*ny*nz /= nn) stop "nx*ny*nz /= nn in row_to_col"
+    idx=1
+    do i = 1, nx
+        do j = 1, ny
+            do k = 1, nz
+                l = (j - 1) * nx + (i - 1) + (k - 1) * nx * ny + 1
+                v_col(l) = v_row(idx)
+                idx=idx+1
+            end do
+        end do
+    end do
+
+end subroutine

test/get_dij.f90

@@ -15,8 +15,8 @@ subroutine calculate_dij(natom,ntypat,ixc,xclevel,nfft,nspden,xred,ucvol,gprimd,
 
     write(*,*) '4. Generating dij from on-site rhoij, v_ks and v_xc'
 
-    ! Note : in PAW calculations, v_h is calculated with ncomp + ntilde + ncoretilde
-    ! while v_xc is calculated with ntilde + ncoretilde
+    ! Note : in Abinit PAW calculations, v_h is calculated with ncomp + ntilde + ncoretilde
+    ! while v_xc is calculated with ntilde + ncoretilde by default
 
     nucdipmom = 0.0
     qphon = 0.0

test/get_vloc_ncoret.f90

@@ -14,6 +14,7 @@ subroutine get_vloc_ncoret(ngfftdg,nfft,natom,ntypat,rprimd,gprimd,gmet,ucvol,xr
     integer :: nfft, mgrid, qprtrb(3)
     integer :: ngfftdg(3)
     integer :: natom, ntypat
+    integer :: optgr, optn, optn2, optv, optstr, optatm, optdyfr, opteltfr
     real*8  :: rprimd(3,3),gprimd(3,3),gmet(3,3)
     real*8  :: ucvol,xred(3,natom)
 
@@ -27,14 +28,36 @@ subroutine get_vloc_ncoret(ngfftdg,nfft,natom,ntypat,rprimd,gprimd,gmet,ucvol,xr
     rcut = 0d0
     vprtrb = 0d0
 
+    optatm  = 1
+    optdyfr = 0
+    opteltfr = 0
+    optgr   = 0
+    optn    = 1
+    optn2   = 1
+    optstr  = 0
+    optv    = 1
+
     call getph(atindx,natom,ngfftdg(1),ngfftdg(2),ngfftdg(3),ph1d,xred)
 
-    call atm2fft(atindx1,ncoret,vloc,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcutdg, &
-        & mgrid,mqgrid,natom,nattyp,nfft,ngfftdg,ntypat,1,0,0,0,1,1,0,1, &
-        & pawtab,ph1d,qgrid_vl,qprtrb,rcut,dummy5,rprimd,dummy6,dummy7,&
-        & ucvol,1,dummy8,dummy8,dummy8,vprtrb,vlspl,&
-        & ngfftdg(2),fftn2_distrib,ffti2_local,ngfftdg(3),fftn3_distrib,ffti3_local)
-
+    allocate(dummy(3*3*natom*optn*optdyfr))
+    allocate(dummy1(3*3*natom*optn*optdyfr))
+    allocate(dummy2(3*3*natom*optn*optdyfr))
+    allocate(dummy3(3*natom*optn*optgr))
+    allocate(dummy4(3*natom*optn*optgr))
+    allocate(dummy5(2*nfft*optv*max(optgr,optstr,optdyfr,opteltfr)))
+    !allocate(dummy6(6*optn*optstr))
+    !allocate(dummy7(6*optv*optstr))
+    allocate(dummy8(2*nfft*optn*opteltfr))
+    allocate(dummy9(6+3*natom*6))
+
+
+    call atm2fft(atindx1,ncoret,vloc,dummy,dummy2,dummy9,dummy1,gmet,gprimd,&
+                dummy3,dummy4,gsqcutdg, mgrid,mqgrid,natom,nattyp,nfft,ngfftdg,ntypat,&
+                1,0,0,optgr,optn,optn2,optstr,optv, &
+                !1,0,0,0,1,1,0,1, &
+                & pawtab,ph1d,qgrid_vl,qprtrb,rcut,dummy5,rprimd,dummy6,dummy7,ucvol,1,dummy8,dummy8,dummy8,vprtrb,vlspl,&
+                & ngfftdg(2),fftn2_distrib,ffti2_local,ngfftdg(3),fftn3_distrib,ffti3_local)
+
     !to prevent leakage
     if(allocated(dummy)) deallocate(dummy)
     if(allocated(dummy1)) deallocate(dummy1)

test/input

@@ -1,14 +1,14 @@
 ecut 10
-ecutpaw 40
+ecutpaw 10
 gmet 0.01 0 0 0 0.01 0 0 0 0.01
 rprimd 10.0 0 0 0 10.0 0 0 0 10.0
 gprimd 0.1 0 0 0 0.1 0 0 0 0.1
 ucvol 1000
 ngfft 30 30 30
-ngfftdg 60 60 60
+ngfftdg 30 30 30
 natom 5
 ntypat 2
-typat 2 1 1 1 1
+typat 1 2 2 2 2
 xred -0.279789547400000       7.109405980000000E-002  0.000000000000000E+000
  -0.212391634800000      -0.119543389500000       0.000000000000000E+000
  -0.212388153900000       0.166411496700000       0.165096194500000

test/pawfiles

@@ -1,2 +1,2 @@
-   H.xml
    C.xml
+   H.xml

test/reference_data/dij_ref

@@ -0,0 +1,96 @@
+0.322420488884
+-2.98810278859
+25.249948999
+1.23705418153e-06
+-1.53339853181e-05
+-0.0125549406562
+1.03315746154e-06
+-2.88821455395e-05
+2.31334401856e-06
+-0.0126398760676
+-0.000113182074431
+0.00318124441347
+9.20025442006e-05
+-2.53554578503e-05
+-0.015146178977
+-4.32744204436e-06
+5.33712838448e-05
+-0.355126968861
+-1.53571838868e-05
+-0.000611630168176
+6.36559701714
+2.03861552561e-05
+-9.5513499718e-05
+-1.53571838868e-05
+-0.354564829608
+0.000168322750795
+9.4487369401e-05
+6.36215093761
+-0.00225956473644
+0.0106141521312
+-0.000611630168176
+0.000168322750795
+-0.337923230331
+0.00376152088404
+-0.00103563088456
+6.25974595937
+-0.00298448706326
+-0.00328294555739
+-0.00283397563354
+-0.000214613796303
+-0.000342700524797
+0.00228648707883
+1.26356916244e-06
+2.02095706897e-06
+-5.08315523274e-06
+0.00229906653402
+0.000136946874125
+0.000218832273916
+6.37715618417e-06
+-8.3161495404e-08
+0.00234479084662
+-0.00305240728769
+-0.00339123241117
+-0.00300659677946
+9.74553561477e-05
+0.000155598040539
+0.0023143560015
+0.000164780309616
+0.000263077250009
+-3.22776030359e-06
+0.00226445180481
+-7.95460897368e-05
+-0.000127412970014
+-8.75723430872e-07
+-3.69024382727e-06
+0.00232700573936
+-0.00294402227344
+-0.00321839905325
+-0.0027310311822
+0.000100273352205
+0.000160105153581
+0.00230074461143
+-0.000176771497131
+-0.000282255998505
+1.26650666058e-06
+0.00229808434184
+-5.95360944308e-05
+-9.54089114584e-05
+3.42033748632e-07
+2.70462397803e-06
+0.0023452188578
+-0.00303269622931
+-0.00335983580357
+-0.0029565917116
+-2.19519255413e-05
+-3.51106658483e-05
+0.00228379227617
+7.61874869959e-06
+1.2185454108e-05
+-5.68357484387e-06
+0.00229669722719
+-0.000412422381197
+-0.000659038918331
+-3.35862861831e-06
+1.95304204929e-06
+0.00233309379519