TABLE OF CONTENTS
ABINIT/printmagvtk [ Functions ]
NAME
printmagvtk
FUNCTION
Auxiliary routine for printing out magnetization density in VTK format. Output file name is DEN.vtk
COPYRIGHT
Copyright (C) 2017-2018 ABINIT group (SPr) 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
mpi_enreg = information about adopted parallelization strategy
nspden = number of components of density matrix (possible values re 1,2, or 4)
nspden: 1 -> rho
nspden: 2 -> rho_up,rho_dwn
nspden: 4 -> rho,mx,my,mz
nfft = number of fft points per FFT processor
ngfft = full information about FFT mesh
rhor = density array stored in the memory of current FFT processor
rprimd = array of lattice vectors
OUTPUT
SIDE EFFECTS
NOTES
At the moment this routine is mainly used for development and debugging
of gs and dfpt calculations with non-collinear spins. If needed, can be used
to print final density in vtk format.
IMPORTANT: implementation is thoroughly checked only for npspinor = 1,
for other case might need to change the part gathering
the FFT mesh info
PARENTS
CHILDREN
xmpi_gather
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 53 subroutine printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,fname) 54 55 use defs_basis 56 use m_errors 57 use m_profiling_abi 58 use defs_abitypes 59 use m_xmpi 60 use m_io_tools, only : open_file 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'printmagvtk' 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 type(MPI_type),intent(in) :: mpi_enreg 73 integer,intent(in) :: nfft,nspden,cplex 74 !arrays 75 integer,intent(in) :: ngfft(18) 76 real(dp),intent(in) :: rhor(cplex*nfft,nspden),rprimd(3,3) 77 character(len=*),intent(in) :: fname 78 79 !Local variables------------------------------- 80 !scalars 81 integer :: denvtk,denxyz,denxyz_im,nfields 82 integer :: nx,ny,nz,nfft_tot 83 integer :: ii,jj,kk,ind,ispden 84 integer :: mpi_comm,mpi_head,mpi_rank,ierr 85 real :: rx,ry,rz 86 integer :: nproc_fft,ir 87 character(len=500) :: msg 88 character(len=10) :: outformat 89 character(len=50) :: fname_vtk 90 character(len=50) :: fname_xyz 91 character(len=50) :: fname_xyz_re 92 character(len=50) :: fname_xyz_im 93 !arrays 94 real(dp),allocatable :: rhorfull(:,:) 95 96 97 ! ************************************************************************* 98 99 DBG_ENTER("COLL") 100 101 ! if (option/=1 .and. option/=2 ) then 102 ! write(msg,'(3a,i0)')& 103 !& 'The argument option should be 1 or 2,',ch10,& 104 !& 'however, option=',option 105 ! MSG_BUG(msg) 106 ! end if 107 ! 108 ! if (sizein<1) then 109 ! write(msg,'(3a,i0)')& 110 !& 'The argument sizein should be a positive number,',ch10,& 111 !& 'however, sizein=',sizein 112 ! MSG_ERROR(msg) 113 ! end if 114 115 DBG_EXIT("COLL") 116 117 fname_vtk=adjustl(adjustr(fname)//".vtk") 118 fname_xyz=adjustl(adjustr(fname)//".xyz") 119 fname_xyz_re=adjustl(adjustr(fname)//"_re.xyz") 120 fname_xyz_im=adjustl(adjustr(fname)//"_im.xyz") 121 !write(std_out,*) ' Writing out .vtk file: ',fname_vtk 122 !write(std_out,*) ' Writing out .xyz file: ',fname_xyz 123 124 !if 1 or two component density then write out either 1 or 2 scalar density fields 125 !if 4, then write one scalar field (density) and one vector field (magnetization density) 126 if(nspden/=4)then 127 nfields=nspden 128 else 129 nfields=2 130 end if 131 132 nfields=nfields*cplex 133 134 ! FFT mesh specifications: full grid 135 nx=ngfft(1) ! number of points along 1st lattice vector 136 ny=ngfft(2) ! number of points along 2nd lattice vector 137 nz=ngfft(3) ! number of points along 3rd lattice vector 138 nfft_tot=nx*ny*nz ! total number of fft mesh points (can be different from nfft in case of distributed memory of nproc_fft processors) 139 140 141 ! Gather information about memory distribution 142 mpi_head=0 143 mpi_comm = mpi_enreg%comm_fft 144 mpi_rank = xmpi_comm_rank(mpi_comm) 145 nproc_fft=ngfft(10) 146 147 ! Create array to host full FFT mesh 148 if(mpi_rank==mpi_head)then 149 ABI_ALLOCATE(rhorfull,(cplex*nfft_tot,nspden)) 150 end if 151 152 ! Fill in the full mesh 153 if(nproc_fft==1)then 154 rhorfull=rhor 155 else 156 do ir=1,nspden 157 call xmpi_gather(rhor(:,ir),cplex*nfft,rhorfull(:,ir),cplex*nfft,mpi_head,mpi_comm,ierr) 158 end do 159 end if 160 161 if(mpi_rank==mpi_head)then 162 163 ! Open the output vtk file 164 if (open_file(fname_vtk,msg,newunit=denvtk,status='replace',form='formatted') /=0) then 165 MSG_WARNING(msg) 166 RETURN 167 end if 168 169 if(cplex==1) then 170 if (open_file(fname_xyz,msg,newunit=denxyz,status='replace',form='formatted') /=0) then 171 MSG_WARNING(msg) 172 RETURN 173 end if 174 else if (cplex==2) then 175 if (open_file(fname_xyz_re,msg,newunit=denxyz,status='replace',form='formatted') /=0) then 176 MSG_WARNING(msg) 177 RETURN 178 end if 179 if (open_file(fname_xyz_im,msg,newunit=denxyz_im,status='replace',form='formatted') /=0) then 180 MSG_WARNING(msg) 181 RETURN 182 end if 183 end if 184 185 ! Write the header of the output vtk file 186 write(denvtk,"(a)") '# vtk DataFile Version 2.0' 187 write(denvtk,"(a)") 'Electron density components' 188 write(denvtk,"(a)") 'ASCII' 189 write(denvtk,"(a)") 'DATASET STRUCTURED_GRID' 190 write(denvtk,"(a,3i6)") 'DIMENSIONS ', nx,ny,nz 191 write(denvtk,"(a,i18,a)") 'POINTS ',nfft_tot,' double' 192 193 if (nspden==1) then 194 outformat="(4e16.8)" 195 else if (nspden==2) then 196 outformat="(5e16.8)" 197 else 198 outformat="(7e16.8)" 199 end if 200 201 ! Write out information about grid points 202 do kk=0,nz-1 203 do jj=0,ny-1 204 do ii=0,nx-1 205 206 rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 207 ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 208 rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 209 write(denvtk,'(3f16.8)') rx,ry,rz !coordinates of the grid point 210 ind=1+ii+nx*(jj+ny*kk) 211 if (cplex==1) then 212 write(denxyz,outformat) rx,ry,rz,(rhorfull(ind,ispden),ispden=1,nspden) 213 else 214 write(denxyz,outformat) rx,ry,rz,(rhorfull(2*ind-1,ispden),ispden=1,nspden) 215 write(denxyz_im,outformat) rx,ry,rz,(rhorfull(2*ind ,ispden),ispden=1,nspden) 216 end if 217 end do 218 end do 219 end do 220 221 if(cplex==1) then 222 close(denxyz) 223 else 224 close(denxyz) 225 close(denxyz_im) 226 end if 227 228 ! Write out information about field defined on the FFT mesh 229 write(denvtk,"(a,i18)") 'POINT_DATA ',nfft_tot 230 write(denvtk,"(a,i6)") 'FIELD Densities ',nfields 231 232 233 ! Write out different fields depending on the number of density matrix components 234 if(nspden==1)then 235 236 !single component, so just write out the density 237 if(cplex==1) then 238 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 239 do kk=0,nz-1 240 do jj=0,ny-1 241 do ii=0,nx-1 242 ind=1+ii+nx*(jj+ny*kk) 243 write(denvtk,'(f16.8)') rhorfull(ind,1) 244 end do 245 end do 246 end do 247 else 248 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 249 do kk=0,nz-1 250 do jj=0,ny-1 251 do ii=0,nx-1 252 ind=2*(1+ii+nx*(jj+ny*kk))-1 253 write(denvtk,'(f16.8)') rhorfull(ind,1) 254 end do 255 end do 256 end do 257 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 258 do kk=0,nz-1 259 do jj=0,ny-1 260 do ii=0,nx-1 261 ind=2*(1+ii+nx*(jj+ny*kk)) 262 write(denvtk,'(f16.8)') rhorfull(ind,1) 263 end do 264 end do 265 end do 266 end if 267 268 else if(nspden==2)then 269 270 !two component, write the density for spin_up and spin_down channels 271 if(cplex==1) then 272 273 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 274 do kk=0,nz-1 275 do jj=0,ny-1 276 do ii=0,nx-1 277 ind=1+ii+nx*(jj+ny*kk) 278 write(denvtk,'(f16.8)') rhorfull(ind,1) 279 end do 280 end do 281 end do 282 write(denvtk,"(a,i18,a)") 'mag 1 ',nfft_tot,' double' 283 do kk=0,nz-1 284 do jj=0,ny-1 285 do ii=0,nx-1 286 ind=1+ii+nx*(jj+ny*kk) 287 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 288 end do 289 end do 290 end do 291 292 else 293 294 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 295 do kk=0,nz-1 296 do jj=0,ny-1 297 do ii=0,nx-1 298 ind=2*(1+ii+nx*(jj+ny*kk))-1 299 write(denvtk,'(f16.8)') rhorfull(ind,1) 300 end do 301 end do 302 end do 303 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 304 do kk=0,nz-1 305 do jj=0,ny-1 306 do ii=0,nx-1 307 ind=2*(1+ii+nx*(jj+ny*kk)) 308 write(denvtk,'(f16.8)') rhorfull(ind,1) 309 end do 310 end do 311 end do 312 write(denvtk,"(a,i18,a)") 'Re_mag 1 ',nfft_tot,' double' 313 do kk=0,nz-1 314 do jj=0,ny-1 315 do ii=0,nx-1 316 ind=2*(1+ii+nx*(jj+ny*kk))-1 317 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 318 end do 319 end do 320 end do 321 write(denvtk,"(a,i18,a)") 'Im_mag 1 ',nfft_tot,' double' 322 do kk=0,nz-1 323 do jj=0,ny-1 324 do ii=0,nx-1 325 ind=2*(1+ii+nx*(jj+ny*kk)) 326 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 327 end do 328 end do 329 end do 330 331 end if 332 333 else !here is the last option: nspden==4 334 335 if(cplex==1) then 336 337 !four component, write the density (scalar field) and magnetization density (vector field) 338 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 339 do kk=0,nz-1 340 do jj=0,ny-1 341 do ii=0,nx-1 342 ind=1+ii+nx*(jj+ny*kk) 343 write(denvtk,'(f16.8)') rhorfull(ind,1) 344 end do 345 end do 346 end do 347 write(denvtk,"(a,i18,a)") 'mag 3 ',nfft_tot,' double' 348 do kk=0,nz-1 349 do jj=0,ny-1 350 do ii=0,nx-1 351 ind=1+ii+nx*(jj+ny*kk) 352 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 353 end do 354 end do 355 end do 356 357 else 358 359 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 360 do kk=0,nz-1 361 do jj=0,ny-1 362 do ii=0,nx-1 363 ind=2*(1+ii+nx*(jj+ny*kk))-1 364 write(denvtk,'(f16.8)') rhorfull(ind,1) 365 end do 366 end do 367 end do 368 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 369 do kk=0,nz-1 370 do jj=0,ny-1 371 do ii=0,nx-1 372 ind=2*(1+ii+nx*(jj+ny*kk)) 373 write(denvtk,'(f16.8)') rhorfull(ind,1) 374 end do 375 end do 376 end do 377 write(denvtk,"(a,i18,a)") 'Re_mag 3 ',nfft_tot,' double' 378 do kk=0,nz-1 379 do jj=0,ny-1 380 do ii=0,nx-1 381 ind=2*(1+ii+nx*(jj+ny*kk))-1 382 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 383 end do 384 end do 385 end do 386 write(denvtk,"(a,i18,a)") 'Im_mag 3 ',nfft_tot,' double' 387 do kk=0,nz-1 388 do jj=0,ny-1 389 do ii=0,nx-1 390 ind=2*(1+ii+nx*(jj+ny*kk)) 391 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 392 end do 393 end do 394 end do 395 396 end if 397 398 end if ! nspden options condition 399 400 close (denvtk) 401 402 !clean up the gathered FFT mesh 403 ABI_DEALLOCATE(rhorfull) 404 405 end if 406 407 408 end subroutine printmagvtk