1!===============================================================================
2!
3! Program:
4!
5! (1) paraband            Originally By FHJ      Last Modified Apr/2015 (FHJ)
6!
7! Input is read from file split_spin.inp.
8!
9!===============================================================================
10
11#include "f_defs.h"
12
13program split_spin
14  use global_m
15  use inread_m,         only: pb_params_t, inread
16  use wfn_rho_vxc_io_m, only: read_mf_header, read_binary_gvectors, read_binary_complex_data, &
17                              write_mf_header, write_binary_gvectors, write_binary_complex_data
18  use wfn_io_m,         only: wfn_read
19  use io_utils_m,       only: progress_info, progress_init, progress_step, &
20                              progress_free, print_dealing_with
21  use pseudopot_vkb_m,  only: pseudopot_vkb_t, compare_mf_headers
22
23  implicit none
24
25  type(mf_header_t) :: mf_wfn, mf_vsc, mf_nospin
26  complex(DPC), allocatable :: vsc(:,:), vsc_nospin(:,:), wfn_buf(:,:)
27  integer, allocatable :: isort(:,:)
28  integer, allocatable :: gveck_buf(:,:), gvectmp(:,:)
29  complex(DPC), allocatable :: vkb_buf(:,:)
30  type(pseudopot_vkb_t) :: pseudopot, pseudopot_nospin
31  character(len=128) :: fname, solver_name
32  character(len=32) :: stitle
33  type(pb_params_t) :: params
34  type(progress_info) :: prog_info
35  integer :: is, ikb, ik, ii, ib
36
37
38  !call peinfo_init()
39  write(6,'(1x,a, a)') 'Git version: ', TOSTRING(GIT_COMMIT)
40  if (peinf%npes/=1) then
41    call die('You run this application with a single processor', &
42    only_root_writes=.true.)
43  endif
44
45!------------------------
46! Read input parameters
47  write(6,'(1x,a)') 'Reading parameters from file parabands.inp'
48  call inread(params)
49
50
51!-----------------------
52! Split WFN
53  write(*,'(/1x,a)') 'Reading WFN header'
54  call open_file(7, file=params%fname_wfn_in, status='old', form='unformatted')
55  call read_mf_header(7, mf_wfn, iflavor=SCALARSIZE, sheader='WFN', warn=.false., dont_warn_kgrid=.true.)
56  if (mf_wfn%kp%nspin/=2) then
57    call die('Input WFN file is not spin polarized.')
58  endif
59
60  SAFE_ALLOCATE(mf_wfn%gvec%components, (3, mf_wfn%gvec%ng))
61  call read_binary_gvectors(7, mf_wfn%gvec%ng, mf_wfn%gvec%ng, mf_wfn%gvec%components)
62  SAFE_ALLOCATE(gvectmp, (3,mf_wfn%kp%ngkmax))
63  SAFE_ALLOCATE(wfn_buf, (mf_wfn%kp%ngkmax,mf_wfn%kp%nspin))
64
65  mf_nospin = mf_wfn
66  associate(kp=>mf_nospin%kp)
67    kp%nspin = 1
68    ! By default, assignment of derived types assigns the pointers which are
69    ! members of the types. But if two pointers point to the same memory
70    ! and one of them gets deallocated, all pointers get deallocated.
71    nullify(kp%ifmin)
72    nullify(kp%ifmax)
73    nullify(kp%el)
74    nullify(kp%occ)
75    SAFE_ALLOCATE(kp%ifmin, (kp%nrk,kp%nspin))
76    SAFE_ALLOCATE(kp%ifmax, (kp%nrk,kp%nspin))
77    SAFE_ALLOCATE(kp%el, (kp%mnband,kp%nrk,kp%nspin))
78    SAFE_ALLOCATE(kp%occ, (kp%mnband,kp%nrk,kp%nspin))
79  endassociate
80
81  do is = 1, 2
82    associate(kp=>mf_nospin%kp, kp_wfn=>mf_wfn%kp)
83      kp%ifmin(:,1) = kp_wfn%ifmin(:,is)
84      kp%ifmax(:,1) = kp_wfn%ifmax(:,is)
85      kp%el(:,:,1) = kp_wfn%el(:,:,is)
86      kp%occ(:,:,1) = kp_wfn%occ(:,:,is)
87    endassociate
88    write(fname,'(a,i0)') 'WFN_spin_', is
89    call open_file(10+is, file=trim(fname), status='replace', form='unformatted')
90    call write_mf_header(10+is, mf_nospin)
91    call write_binary_gvectors(10+is, mf_nospin%gvec%ng, mf_nospin%gvec%ng, &
92      mf_nospin%gvec%components)
93  enddo
94
95  call progress_init(prog_info, 'copying WFN', 'state', mf_wfn%kp%nrk*mf_wfn%kp%mnband*2)
96  do ik = 1, mf_wfn%kp%nrk
97    call read_binary_gvectors(7, mf_wfn%kp%ngk(ik), mf_wfn%kp%ngkmax, gvectmp)
98    do is = 1, 2
99      call write_binary_gvectors(10+is, mf_nospin%kp%ngk(ik), mf_nospin%kp%ngkmax, gvectmp)
100    enddo
101    do ib = 1, mf_wfn%kp%mnband
102      call read_binary_complex_data(7, mf_wfn%kp%ngk(ik), mf_wfn%kp%ngkmax, mf_wfn%kp%nspin, &
103        wfn_buf)
104      do is = 1, 2
105         call progress_step(prog_info)
106        call write_binary_complex_data(10+is, mf_nospin%kp%ngk(ik), mf_nospin%kp%ngkmax, mf_nospin%kp%nspin, &
107          wfn_buf(:,is:is))
108      enddo
109    enddo
110  enddo
111  call close_file(7)
112  call close_file(11)
113  call close_file(12)
114  SAFE_DEALLOCATE(gvectmp)
115  SAFE_DEALLOCATE(wfn_buf)
116  call progress_free(prog_info)
117
118
119!------------------------
120! Split self-consistent potential file
121  write(6,'(/1x,a)') 'Reading self-consistent potential from file '//TRUNC(params%fname_vsc)
122  call open_file(10, file=params%fname_vsc, status='old', form='unformatted')
123  call read_mf_header(10, mf_vsc, iflavor=SCALARSIZE, sheader='VXC', warn=.false., dont_warn_kgrid=.true.)
124  if (mf_vsc%kp%nspin/=2) then
125    call die('Input VSC file is not spin polarized.')
126  endif
127  SAFE_ALLOCATE(mf_vsc%gvec%components, (3, mf_vsc%gvec%ng))
128  call read_binary_gvectors(10, mf_vsc%gvec%ng, mf_vsc%gvec%ng, &
129    mf_vsc%gvec%components, bcast=.false.)
130  call compare_mf_headers('WFN', mf_wfn, 'VSC', mf_vsc, .false.)
131  SAFE_ALLOCATE(vsc, (mf_vsc%gvec%ng, mf_vsc%kp%nspin))
132  call read_binary_complex_data(10, mf_vsc%gvec%ng, mf_vsc%gvec%ng, mf_vsc%kp%nspin, vsc)
133  call close_file(10)
134
135  mf_nospin = mf_vsc
136  mf_nospin%kp%nspin = 1
137  SAFE_ALLOCATE(vsc_nospin, (mf_nospin%gvec%ng, mf_nospin%kp%nspin))
138  do is = 1, 2
139    write(*,'(1x,a,i0)') 'Writing spin-unpolarized VSC for spin ', is
140    vsc_nospin(:,1) = vsc(:,is)
141    write(fname,'(a,i0)') 'VSC_spin_', is
142    call open_file(10, file=trim(fname), status='replace', form='unformatted')
143    call write_mf_header(10, mf_nospin)
144    call write_binary_gvectors(10, mf_nospin%gvec%ng, mf_nospin%gvec%ng, &
145      mf_nospin%gvec%components)
146    call write_binary_complex_data(10, mf_nospin%gvec%ng, mf_nospin%gvec%ng, mf_nospin%kp%nspin, vsc_nospin)
147    call close_file(10)
148  enddo
149  SAFE_DEALLOCATE(vsc_nospin)
150
151
152!------------------------
153! Split Kleinman-Bylander projectors file
154  if (params%has_vkb) then
155    call open_file(7, file=params%fname_vkb, status='old', form='unformatted')
156    write(6,'(/1x,a)') 'Reading Kleinman-Bylander projectors from file ' // &
157      TRUNC(params%fname_vkb)
158    call pseudopot%read_header(7, params%kpp, mf_wfn, allocate_vkb_buf=.false.)
159
160    stitle = 'VKB-Complex'
161    pseudopot_nospin = pseudopot
162    pseudopot_nospin%ns = 1
163    SAFE_ALLOCATE(vkb_buf, (pseudopot%ngkmax,1))
164    SAFE_ALLOCATE(gveck_buf, (3,pseudopot%ngkmax))
165    do is = 1, pseudopot%ns
166      ! Write from now on. We should only have pseudopot_nospin%.. below, and no pseudopot%..
167      !write(*,'(1x,a,i0)') 'Writing spin-unpolarized VKB for spin ', is
168
169      write(fname,'(a,i0)') 'VKB_spin_', is
170      call open_file(10+is, file=trim(fname), status='replace', form='unformatted')
171      write(10+is) stitle, mf_wfn%sdate, mf_wfn%stime
172      write(10+is) 1, mf_wfn%gvec%ng, mf_wfn%syms%ntran, mf_wfn%syms%cell_symmetry, &
173        mf_wfn%crys%nat, mf_wfn%gvec%ecutrho, mf_wfn%kp%nrk, &
174        pseudopot_nospin%nsp, pseudopot_nospin%nkb, pseudopot_nospin%nhm, mf_wfn%kp%ngkmax, mf_wfn%kp%ecutwfc
175      write(10+is) mf_wfn%gvec%FFTgrid, mf_wfn%kp%kgrid, mf_wfn%kp%shift
176      write(10+is) mf_wfn%crys%celvol, mf_wfn%crys%alat, mf_wfn%crys%avec, mf_wfn%crys%adot
177      write(10+is) mf_wfn%crys%recvol, mf_wfn%crys%blat, mf_wfn%crys%bvec, mf_wfn%crys%bdot
178      write(10+is) mf_wfn%syms%mtrx(1:3,1:3,1:mf_wfn%syms%ntran)
179      write(10+is) mf_wfn%syms%tnp(1:3,1:mf_wfn%syms%ntran)
180      write(10+is) (mf_wfn%crys%apos(1:3,ii), mf_wfn%crys%atyp(ii), ii=1,mf_wfn%crys%nat)
181      write(10+is) mf_wfn%kp%ngk(1:pseudopot_nospin%nk)
182      write(10+is) mf_wfn%kp%w(1:pseudopot_nospin%nk)
183      write(10+is) mf_wfn%kp%rk(1:3,1:pseudopot_nospin%nk)
184
185      write(10+is) pseudopot_nospin%ityp(:)
186      write(10+is) pseudopot_nospin%nh(:)
187      write(10+is) pseudopot_nospin%deeq(:,:,:,is:is)
188      call write_binary_gvectors(10+is, mf_wfn%gvec%ng, mf_wfn%gvec%ng, mf_wfn%gvec%components)
189    enddo
190
191    call progress_init(prog_info, 'copying VKB', 'projector', &
192      pseudopot_nospin%nk*pseudopot%ns*pseudopot%nkb)
193    do ik = 1, pseudopot%nk
194      call read_binary_gvectors(7, mf_wfn%kp%ngk(ik), pseudopot%ngkmax, gveck_buf(:,:))
195      do is = 1, pseudopot%ns
196        call write_binary_gvectors(10+is, mf_wfn%kp%ngk(ik), pseudopot_nospin%ngkmax, gveck_buf(:,:))
197        do ikb = 1, pseudopot%nkb
198          call progress_step(prog_info)
199          call read_binary_complex_data(7, mf_wfn%kp%ngk(ik), pseudopot%ngkmax, 1, vkb_buf(:,1:1))
200          call write_binary_complex_data(10+is, mf_wfn%kp%ngk(ik), pseudopot_nospin%ngkmax, 1, vkb_buf(:,1:1))
201        enddo
202      enddo
203    enddo
204    call progress_free(prog_info)
205    call close_file(7)
206    call close_file(11)
207    call close_file(12)
208
209    SAFE_DEALLOCATE(gveck_buf)
210    SAFE_DEALLOCATE(vkb_buf)
211  endif ! params%has_vkb
212
213  write(*,'(/1x,a)') 'All done!'
214
215end program split_spin
216