TABLE OF CONTENTS


ABINIT/out1dm [ Functions ]

[ Top ] [ Functions ]

NAME

 out1dm

FUNCTION

 Output the 1 dimensional mean of potential and density
 on the three reduced axis.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  fnameabo_app_1dm=name of the file in which the data is written, appended with _1DM
  natom=number of atoms in unit cell
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  ntypat=number of types of atoms in unit cell.
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
  rprimd(3,3)=real space dimensional primitive translations (bohr)
  typat(natom)=type integer for each atom in cell
  ucvol=unit cell volume in bohr**3.
  vtrial(nfft,nspden)=INPUT Vtrial(r).
  xred(3,natom)=reduced coordinates of atoms
  znucl(ntypat)=real(dp), nuclear number of atom type

OUTPUT

SIDE EFFECTS

  data written in file fnameabo_app_1dm

NOTES

PARENTS

      outscfcv

CHILDREN

      atomdata_from_znucl,ptabs_fourdp,wrtout,xmpi_sum,xred2xcart

SOURCE

 47 #if defined HAVE_CONFIG_H
 48 #include "config.h"
 49 #endif
 50 
 51 #include "abi_common.h"
 52 
 53 subroutine out1dm(fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,&
 54 &  rhor,rprimd,typat,ucvol,vtrial,xred,znucl)
 55 
 56  use defs_basis
 57  use defs_abitypes
 58  use m_errors
 59  use m_profiling_abi
 60  use m_atomdata
 61 
 62  use m_io_tools, only : open_file
 63  use m_geometry,  only : xred2xcart
 64  use m_mpinfo,   only : ptabs_fourdp
 65  use m_xmpi,     only : xmpi_sum, xmpi_comm_size
 66 
 67 !This section has been created automatically by the script Abilint (TD).
 68 !Do not modify the following lines by hand.
 69 #undef ABI_FUNC
 70 #define ABI_FUNC 'out1dm'
 71  use interfaces_14_hidewrite
 72 !End of the abilint section
 73 
 74  implicit none
 75 
 76 !Arguments ------------------------------------
 77 !scalars
 78  integer,intent(in) :: natom,nfft,nspden,ntypat
 79  real(dp),intent(in) :: ucvol
 80  character(len=fnlen),intent(in) :: fnameabo_app_1dm
 81  type(MPI_type),intent(in) :: mpi_enreg
 82 !arrays
 83  integer,intent(in) :: ngfft(18),typat(natom)
 84  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),vtrial(nfft,nspden)
 85  real(dp),intent(in) :: znucl(ntypat),xred(3,natom)
 86 
 87 !Local variables-------------------------------
 88 !scalars
 89  integer :: ia,ib,idim,ierr,ifft,islice,ispden,na,nb,ndig,nslice,nu,temp_unit
 90  integer :: comm_fft, nproc_fft, me_fft
 91  real(dp) :: global_den,global_pot
 92  character(len=2) :: symbol
 93  character(len=500) :: message
 94  type(atomdata_t) :: atom
 95 !arrays
 96  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
 97  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 98  real(dp),allocatable :: lin_den(:),mean_pot(:),reduced_coord(:),xcart(:,:)
 99  character(len=8),allocatable :: iden(:)
100 
101 ! *************************************************************************
102 
103 !Initialize the file
104  write(message, '(a,a)' ) ' io1dm : about to open file ',fnameabo_app_1dm
105  call wrtout(std_out,message,'COLL')
106  call wrtout(ab_out,message,'COLL')
107 
108  comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
109 
110  if (me_fft == 0) then
111    if (open_file(fnameabo_app_1dm,message,newunit=temp_unit,status='unknown',form='formatted') /= 0) then
112      MSG_ERROR(message)
113    end if
114    rewind(temp_unit)
115  end if
116 
117  write(message, '(a,a)' ) ch10,'# ABINIT package : 1DM file '
118  if (me_fft == 0) write(temp_unit,'(a)') message
119 
120  write(message, '(a,a)' )ch10,'# Primitive vectors of the periodic cell (bohr)'
121  if (me_fft == 0) write(temp_unit,'(a)') message
122  do nu=1,3
123    write(message, '(1x,a,i1,a,3f10.5)' ) '#  R(',nu,')=',rprimd(:,nu)
124    if (me_fft == 0) write(temp_unit,'(a)') message
125  end do
126 
127  write(message, '(a,a)' ) ch10,&
128 & '# Atom list        Reduced coordinates          Cartesian coordinates (bohr)'
129  if (me_fft == 0) write(temp_unit,'(a)') message
130 
131 !Set up a list of character identifiers for all atoms : iden(ia)
132  ABI_ALLOCATE(iden,(natom))
133  do ia=1,natom
134    call atomdata_from_znucl(atom,znucl(typat(ia)))
135    symbol = atom%symbol
136 
137    ndig=int(log10(dble(ia)+0.5_dp))+1
138    if(ndig==1) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')   '
139    if(ndig==2) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')  '
140    if(ndig==3) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,') '
141    if(ndig==4) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')'
142    if(ndig>4)then
143      close(temp_unit)
144      write(message, '(a,i0)' )&
145 &     ' out1dm : cannot handle more than 9999 atoms, while natom=',natom
146      MSG_ERROR(message)
147    end if
148  end do
149 
150 !Compute cartesian coordinates, and print reduced and cartesian coordinates
151  ABI_ALLOCATE(xcart,(3,natom))
152  call xred2xcart(natom,rprimd,xcart,xred)
153  do ia=1,natom
154    write(message, '(a,a,3f10.5,a,3f10.5)' ) &
155 &   '#   ',iden(ia),xred(1:3,ia),'    ',xcart(1:3,ia)
156    if (me_fft == 0) write(temp_unit,'(a)') message
157  end do
158  ABI_DEALLOCATE(iden)
159  ABI_DEALLOCATE(xcart)
160 
161 !Get the distrib associated with this fft_grid
162  call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
163 
164  do idim=1,3
165 
166    nslice=ngfft(idim)
167 
168 !  Dummy initialisation of na, nb and ifft.
169    na=0 ; nb=0; ifft=0
170    select case(idim)
171    case(1)
172      na=ngfft(2) ; nb=ngfft(3)
173    case(2)
174      na=ngfft(1) ; nb=ngfft(3)
175    case(3)
176      na=ngfft(1) ; nb=ngfft(2)
177    end select
178 
179    ABI_ALLOCATE( reduced_coord,(nslice))
180    ABI_ALLOCATE(mean_pot,(nslice))
181    ABI_ALLOCATE(lin_den,(nslice))
182 
183    do ispden=1,nspden
184 
185      if(ispden==1)then
186        write(message, '(a,a,a)' ) ch10,'#===========',&
187 &       '====================================================================='
188        if (me_fft == 0) write(temp_unit,'(a)') message
189      end if
190 
191      select case(idim)
192      case(1)
193        write(message, '(a)' )'# Projection along the first dimension '
194      case(2)
195        write(message, '(a)' )'# Projection along the second dimension '
196      case(3)
197        write(message, '(a)' )'# Projection along the third dimension '
198      end select
199      if (me_fft == 0) write(temp_unit,'(a)') message
200 
201      if(nspden==2)then
202        select case(ispden)
203        case(1)
204          write(message, '(a)' )'# Spin up '
205        case(2)
206          write(message, '(a)' )'# Spin down '
207        end select
208        if (me_fft == 0) write(temp_unit,'(a)') message
209      else if (nspden == 4) then
210        select case(ispden)
211        case(1)
212          write(message, '(a)' )'# Spinor up up '
213        case(2)
214          write(message, '(a)' )'# Spinor down down'
215        case(3)
216          write(message, '(a)' )'# Spinor up down'
217        case(4)
218          write(message, '(a)' )'# Spinor down up'
219        end select
220        if (me_fft == 0) write(temp_unit,'(a)') message
221      end if
222 
223      write(message, '(2a)' ) ch10,&
224 &     '#     Red. coord. Mean KS potential  Linear density  '
225      if (me_fft == 0) write(temp_unit,'(a)') message
226 
227      write(message, '(a)' ) &
228 &     '#                  (Hartree unit)   (electron/red. unit)'
229      if (me_fft == 0) write(temp_unit,'(a)') message
230 
231      global_pot=zero
232      global_den=zero
233      mean_pot(:)=zero
234      lin_den(:)=zero
235      do islice=1,nslice
236        reduced_coord(islice)=(islice-1)/dble(nslice)
237      end do
238      if (idim == 1) then
239        do islice=1,nslice
240          do ib=1,nb
241 ! skip z values not on this processor
242            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
243            do ia=1,na
244              ifft = islice + nslice*( ia    -1 + na    *(ffti3_local(ib)    -1) )
245              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
246              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
247            end do
248          end do
249        end do
250      else if (idim == 2) then
251        do islice=1,nslice
252          do ib=1,nb
253 ! skip z values not on this processor
254            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
255            do ia=1,na
256              ifft = ia     + na    *( islice-1 + nslice*(ffti3_local(ib)    -1) )
257              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
258              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
259            end do
260          end do
261        end do
262      else if (idim == 3) then
263 ! this full z slice is mine in parallel case
264        do islice=1,nslice
265          if (fftn3_distrib(islice)/=mpi_enreg%me_fft) cycle
266          do ib=1,nb
267            do ia=1,na
268              ifft = ia     + na    *( ib    -1 + nb    *(ffti3_local(islice)-1) )
269              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
270              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
271            end do
272          end do
273        end do
274      end if
275      mean_pot(:)=mean_pot(:)/dble(na*nb)
276      lin_den(:)=lin_den(:)/dble(na*nb)*ucvol
277      global_pot=sum(mean_pot)/dble(nslice)
278      global_den=sum(lin_den)/dble(nslice)
279 
280 !    FFT parallelization
281      if(mpi_enreg%nproc_fft>1)then
282        call xmpi_sum(mean_pot  ,mpi_enreg%comm_fft,ierr)
283        call xmpi_sum(global_pot,mpi_enreg%comm_fft,ierr)
284        call xmpi_sum(lin_den   ,mpi_enreg%comm_fft,ierr)
285        call xmpi_sum(global_den,mpi_enreg%comm_fft,ierr)
286      end if
287 
288      do islice=1,ngfft(idim)
289        write(message, '(f10.4,es20.6,es16.6)' )&
290 &       reduced_coord(islice),mean_pot(islice),lin_den(islice)
291        if (me_fft == 0) write(temp_unit,'(a)') message
292      end do
293 
294      write(message, '(a,a,es15.6,es16.6,a)' ) ch10,&
295 &     '# Cell mean       :',global_pot,global_den, ch10
296      if (me_fft == 0) write(temp_unit,'(a)') message
297 
298 
299 !    End of the loop on spins
300    end do
301 
302    ABI_DEALLOCATE(reduced_coord)
303    ABI_DEALLOCATE(mean_pot)
304    ABI_DEALLOCATE(lin_den)
305 
306 !  End of the loops on the three dimensions
307  end do
308 
309  if (me_fft == 0) then
310    write(message, '(a,a,a)' ) ch10,'#===========',&
311 &   '====================================================================='
312    call wrtout(temp_unit,message,'COLL')
313    close(temp_unit)
314  end if
315 
316 end subroutine out1dm