TABLE OF CONTENTS
ABINIT/cut3d [ Programs ]
NAME
cut3d
FUNCTION
Main routine for the analysis of the density and potential files, as well as other files with the ABINIT header.
COPYRIGHT
Copyright (C) 1999-2022 ABINIT group (GMR, RC, LSI, XG, NCJ, JFB, MCote, LPizzagalli) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
(main program)
OUTPUT
(main program)
NOTES
natom = number of atoms in the unit cell nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension) ntypat = number of atom types ucvol = unit cell volume (> 0) filrho = name of the density file (binary or netcdf)
SOURCE
31 #if defined HAVE_CONFIG_H 32 #include "config.h" 33 #endif 34 35 #include "abi_common.h" 36 37 program cut3d 38 39 use defs_basis 40 use m_errors 41 use m_build_info 42 use m_xmpi 43 use m_nctk 44 use m_abicore 45 #ifdef HAVE_NETCDF 46 use netcdf 47 #endif 48 #if defined FC_NAG 49 use f90_unix_proc 50 #endif 51 use m_hdr 52 use m_cut3d 53 use m_crystal 54 55 use defs_abitypes, only : MPI_type 56 use m_specialmsg, only : specialmsg_getcount, herald 57 use m_fstrings, only : endswith, sjoin, itoa 58 use m_time, only : timein 59 use m_geometry, only : xred2xcart, metric 60 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 61 use m_fftcore, only : ngfft_seq 62 use m_fft_mesh, only : denpot_project 63 use m_distribfft, only : init_distribfft_seq 64 use m_ioarr, only : fftdatar_write 65 use m_io_tools, only : flush_unit, file_exists, open_file, is_open, get_unit, read_string 66 67 implicit none 68 69 !Local variables------------------------------- 70 character(len=1) :: outputchar,blank=' ' 71 !scalars 72 integer,parameter :: mfiles=10,exchn2n3d0=0 73 integer :: fform0,gridshift1,gridshift2,gridshift3,i1,i2,i3 74 integer :: iatom,ifiles,ii,ii1,ii2,ii3,index,iprompt,ir1,ir2,ir3,ispden,cplex 75 integer :: itask,jfiles,natom,nfiles,nr1,nr2,unt,comm,iomode,nprocs,my_rank 76 integer :: nr3,nr1_stored,nr2_stored,nr3_stored,nrws,nspden,nspden_stored,ntypat,nfft 77 real(dp) :: dotdenpot,maxmz,normz,sumdenpot,ucvol,xm,xnow,xp,ym,ynow,yp,zm,znow,zp,tcpui,twalli 78 character(len=24) :: codename 79 character(len=fnlen) :: filnam,filrho,filrho_tmp 80 character(len=nctk_slen) :: varname 81 type(hdr_type) :: hdr 82 type(abifile_t) :: abifile 83 type(MPI_type) :: mpi_enreg 84 type(crystal_t) :: cryst 85 !arrays 86 integer, allocatable :: isdenpot(:) 87 integer :: ngfft(18) 88 real(dp) :: rprimd(3,3),shift_tau(3),tsec(2) 89 real(dp) :: xcart2(3),gmet(3,3),gprimd(3,3),rmet(3,3) 90 real(dp),allocatable :: grid(:,:,:),grid_full(:,:,:,:),grid_full_stored(:,:,:,:,:),gridtt(:,:,:) !, grid_rot(:,:,:,:) 91 real(dp),allocatable :: tau2(:,:),xcart(:,:),xred(:,:),rhomacu(:,:),gridmz(:,:,:),gridmy(:,:,:),gridmx(:,:,:) 92 character(len=fnlen),allocatable :: filrho_stored(:) 93 character(len=500) :: message 94 95 !****************************************************************** 96 97 !Change communicator for I/O (mandatory!) 98 call abi_io_redirect(new_io_comm=xmpi_world) 99 100 !Initialize MPI 101 call xmpi_init() 102 comm = xmpi_world 103 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 104 ABI_CHECK(nprocs == 1, "cut3d not programmed for parallel execution") 105 106 !Initialize memory profiling if it is activated 107 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0" 108 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally 109 #ifdef HAVE_MEM_PROFILING 110 call abimem_init(0) 111 #endif 112 113 call timein(tcpui,twalli) 114 115 !Default for sequential use 116 !Other values of mpi_enreg are dataset dependent, and should NOT be initialized inside cut3d.F90. 117 call initmpi_seq(mpi_enreg) 118 119 codename='CUT3D '//repeat(' ',18) 120 call herald(codename,abinit_version,std_out) 121 122 !BIG LOOP on files 123 ABI_MALLOC(isdenpot,(mfiles)) 124 isdenpot=0 125 ABI_MALLOC(filrho_stored,(mfiles)) 126 iomode = IO_MODE_FORTRAN 127 128 do ifiles=1,mfiles 129 130 ! Get name of density file 131 write(std_out,*) 132 write(std_out,*) ' What is the name of the 3D function (density, potential or wavef) file ?' 133 if (read_string(filrho, unit=std_in) /= 0) then 134 ABI_ERROR("Fatal error!") 135 end if 136 filrho_tmp=adjustl(filrho) 137 do ii=1,len_trim(filrho_tmp) 138 if(filrho_tmp(ii:ii)==blank)then 139 filrho=trim(filrho_tmp(1:ii-1)) 140 exit 141 end if 142 end do 143 write(std_out,*) ' => Your 3D function file is: ',trim(filrho) 144 write(std_out,*) 145 ! Checking the existence of data file 146 if (nctk_try_fort_or_ncfile(filrho, message) /= 0) then 147 ABI_ERROR(message) 148 end if 149 150 ! Treat the different cases: formatted or unformatted 151 iomode = IO_MODE_FORTRAN; if (endswith(filrho, ".nc")) iomode = IO_MODE_ETSF 152 if (iomode == IO_MODE_FORTRAN) then 153 write(std_out,"(a)") '- Your file contains unformatted binary header + 3D data' 154 else 155 write(std_out,"(a)") '- Your file contains ETSF data' 156 end if 157 158 ! Read the header and extract dimensions. 159 write(std_out,*) 160 call hdr_read_from_fname(hdr, filrho, fform0, comm) 161 ABI_CHECK(fform0 /= 0, "hdr_read returned fform = 0") 162 abifile = abifile_from_fform(fform0) 163 ABI_CHECK(abifile%fform /= 0, "Cannot detect abifile from fform") 164 165 ! Echo part of the header 166 call hdr%echo(fform0, 4) 167 168 nr1=hdr%ngfft(1); nr2=hdr%ngfft(2); nr3=hdr%ngfft(3) 169 natom=hdr%natom 170 nspden=hdr%nspden 171 ntypat=hdr%ntypat 172 rprimd(:,:)=hdr%rprimd(:,:) 173 174 ! Need to know natom in order to allocate xcart 175 ABI_MALLOC(xcart,(3,natom)) 176 ABI_MALLOC(xred,(3,natom)) 177 xred(:,:)=hdr%xred(:,:) 178 call xred2xcart(natom,rprimd,xcart,xred) 179 180 ispden=0 181 if (abifile%class == "density" .or. abifile%class == "potential") then 182 if(nspden/=1)then 183 write(std_out,'(a)' )' ' 184 write(std_out,'(a)' )' * This file contains more than one spin component,' 185 write(std_out,'(a,i3,a)' )' (indeed, nspden=',nspden,' )' 186 write(std_out,'(a)' )' Some of the tasks that you will define later will concern all spin components.' 187 write(std_out,'(a)' )' Others tasks might require you to have chosen among the following:' 188 end if 189 if(nspden==2)then 190 write(std_out,'(a)' )' ispden= 0 ==> Total density' 191 write(std_out,'(a)' )' ispden= 1 ==> spin-up density' 192 write(std_out,'(a)' )' ispden= 2 ==> spin-down density' 193 write(std_out,'(a)' )' ispden= 3 ==> spin-polarization (or magnetization) density' 194 write(std_out,'(a)' )' spin up - spin down difference.' 195 end if 196 if(nspden==4)then 197 write(std_out,'(a)' )' ispden= 0 ==> Total density' 198 write(std_out,'(a)' )' ispden= 1 ==> magnetization in the x direction' 199 write(std_out,'(a)' )' ispden= 2 ==> magnetization in the y direction' 200 write(std_out,'(a)' )' ispden= 3 ==> magnetization in the z direction' 201 write(std_out,'(a)' )' ispden= 4 might be used to plot the magnetization (3D) in the XCrysDen format,' 202 end if 203 if(nspden/=1)then 204 write(std_out,*)' Please define ispden:' 205 read(std_in,*)ispden 206 write(std_out,'(a,i3)' )' You entered ispden=',ispden 207 end if 208 end if 209 210 write(std_out,*) 211 write(std_out,*) '===========================================================' 212 write(std_out,*) 213 214 ! Echo the value of different input parameters 215 write(std_out,*)'ECHO important input variables ...' 216 write(std_out,*) 217 write(std_out,*) ' Dimensional primitive vectors (ABINIT equivalent: rprimd):' 218 write(std_out,'(3es16.6)' ) rprimd(1:3,1) 219 write(std_out,'(3es16.6)' ) rprimd(1:3,2) 220 write(std_out,'(3es16.6)' ) rprimd(1:3,3) 221 222 ! Compute ucvol and test the non-collinearity of rprimd vectors. 223 call metric(gmet,gprimd,dev_null,rmet,rprimd,ucvol) 224 225 write(std_out,'(a,3i5)' ) ' Grid density (ABINIT equivalent: ngfft): ',nr1,nr2,nr3 226 write(std_out,*) ' Number of atoms :',natom 227 write(std_out,*) ' Number of atomic types:',ntypat 228 229 write(std_out,*) 230 write(std_out,*) ' # Atomic positions (cartesian coordinates - Bohr)' 231 do iatom=1,natom 232 write(std_out,'(i4,3es16.6)' )iatom,xcart(1:3,iatom) 233 end do 234 write(std_out,*) 235 236 ! ------------------------------------------------------------------------ 237 ! Branching: either WF file, or DEN/POT file. 238 239 if (abifile%class == "wf_planewave") then 240 write(std_out,*)' This file is a WF file. ' 241 isdenpot(ifiles)=0 242 iprompt = 0 ! this needs to be initialized, as it is used after the loop on files... 243 244 call cut3d_wffile(filrho,hdr%ecut_eff,exchn2n3d0,hdr%istwfk,hdr%kptns,natom,hdr%nband,hdr%nkpt,hdr%npwarr,& 245 & nr1,nr2,nr3,hdr%nspinor,hdr%nsppol,ntypat,rprimd,xcart,hdr%typat,hdr%znucltypat) 246 call hdr%free() 247 248 ! ------------------------------------------------------------------------- 249 ! This is a DEN/POT file 250 else if (abifile%class == "density" .or. abifile%class == "potential") then 251 252 ! This should become a subroutine 253 write(std_out,*)' This file is a Density or Potential file ' 254 isdenpot(ifiles)=1 255 256 ! Read the function on the 3D grid 257 ABI_MALLOC(grid,(nr1,nr2,nr3)) 258 ABI_MALLOC(grid_full,(nr1,nr2,nr3,nspden)) 259 ABI_MALLOC(gridtt,(nr1,nr2,nr3)) 260 ABI_MALLOC(gridmx,(nr1,nr2,nr3)) 261 ABI_MALLOC(gridmy,(nr1,nr2,nr3)) 262 ABI_MALLOC(gridmz,(nr1,nr2,nr3)) 263 264 varname = varname_from_fname(filrho) 265 if (iomode == IO_MODE_ETSF) then 266 call wrtout(std_out, sjoin("- Reading netcdf variable: ", varname)) 267 end if 268 269 call cut3d_rrho(filrho,varname,iomode,grid_full,nr1,nr2,nr3,nspden) 270 271 !ABI_WARNING("Computing (rhor(r) + rho(-r)) / 2") 272 !ABI_MALLOC(grid_rot, (nr1,nr2,nr3,nspden)) 273 !call ngfft_seq(ngfft, [nr1, nr2, nr3]) 274 !ngfft(4:6) = ngfft(1:3) 275 !call denpot_project(1, ngfft, nspden, grid_full, inversion_3d, [zero, zero, zero], grid_rot) 276 !grid_full = grid_rot 277 !ABI_FREE(grid_rot) 278 279 ! Do not forget that the first sub-array of a density file is the total density, 280 ! while the first sub-array of a potential file is the spin-up potential 281 if (abifile%class == "density") then 282 283 ! gridtt= grid --> Total density or potential. 284 ! gridmx= grid --> spin-Up density, or magnetization density in X direction. 285 ! gridmy= grid --> spin-Down density, or magnetization density in Y direction. 286 ! gridmz= grid --> spin-polarization density (Magnetization), 287 ! or magnetization density in Z direction. 288 gridtt(:,:,:)=grid_full(:,:,:,1) 289 if(nspden==2)then 290 gridmx = grid_full(:,:,:,2) 291 gridmy = grid_full(:,:,:,1)-grid_full(:,:,:,2) 292 gridmz = -grid_full(:,:,:,1)+two*grid_full(:,:,:,2) 293 else if(nspden==4)then 294 gridmx = grid_full(:,:,:,2) 295 gridmy = grid_full(:,:,:,3) 296 gridmz = grid_full(:,:,:,4) 297 end if 298 299 if(nspden==1)then 300 grid = grid_full(:,:,:,1) 301 else 302 if(ispden==0)then 303 grid = gridtt 304 else if(ispden==1)then 305 grid = gridmx 306 else if(ispden==2)then 307 grid = gridmy 308 else if(ispden==3)then 309 grid = gridmz 310 ! if(ispden==0)then 311 ! grid(:,:,:)=grid_full(:,:,:,1) 312 ! else if(ispden==1)then 313 ! grid(:,:,:)=grid_full(:,:,:,2) 314 ! else if(ispden==2)then 315 ! grid(:,:,:)=grid_full(:,:,:,1)-grid_full(:,:,:,2) 316 ! else if(ispden==-1)then 317 ! grid(:,:,:)=-grid_full(:,:,:,1)+two*grid_full(:,:,:,2) 318 else if(ispden==4)then 319 write(std_out,*) ' ' 320 else 321 ABI_ERROR(sjoin('bad ispden value = ',itoa(ispden))) 322 end if 323 end if 324 325 else if (abifile%class == "potential") then ! Potential case 326 if(ispden==0)then 327 grid(:,:,:)=grid_full(:,:,:,1) 328 else if(ispden==1 .or. ispden==2)then 329 grid(:,:,:)=grid_full(:,:,:,ispden) 330 else 331 ABI_ERROR(sjoin('bad ispden value = ',itoa(ispden))) 332 end if 333 gridtt = grid 334 end if 335 336 write(std_out,*) 337 write(std_out,*) ' 3D function was read. Ready for further treatment.' 338 write(std_out,*) 339 write(std_out,*) '===========================================================' 340 write(std_out,*) 341 342 ! ------------------------------------------------------------------------ 343 344 ! At this moment all the input is done 345 ! The code knows the geometry of the system, 346 ! and the data file (electron density, potential, etc). 347 ! It will further calculate the electron density by interpolation in 348 ! a point, along a line or in a plane. 349 350 do 351 do 352 write(std_out,*) ' What is your choice ? Type:' 353 write(std_out,*) ' 0 => exit' 354 write(std_out,*) ' 1 => point (interpolation of data for a single point)' 355 write(std_out,*) ' 2 => line (interpolation of data along a line)' 356 write(std_out,*) ' 3 => plane (interpolation of data in a plane)' 357 write(std_out,*) ' 4 => volume (interpolation of data in a volume)' 358 write(std_out,*) ' 5 => 3D formatted data (output the bare 3D data - one column)' 359 write(std_out,*) ' 6 => 3D indexed data (bare 3D data, preceeded by 3D index)' 360 write(std_out,*) ' 7 => 3D Molekel formatted data ' 361 write(std_out,*) ' 8 => 3D data with coordinates (tecplot ASCII format)' 362 write(std_out,*) ' 9 => output .xsf file for XCrysDen' 363 write(std_out,*) ' 11 => compute atomic charge using the Hirshfeld method' 364 write(std_out,*) ' 14 => Gaussian/cube wavefunction module' 365 write(std_out,*) ' 15 => Write data to netcdf file' 366 read(std_in,*) itask 367 write(std_out,'(a,a,i2,a)' ) ch10,' Your choice is ',itask,ch10 368 369 if ((5 <= itask .and. itask <= 9) .or. any(itask == [14, 15]) )then 370 write(std_out,*) ch10,' Enter the name of an output file:' 371 if (read_string(filnam, unit=std_in) /= 0) then 372 ABI_ERROR("Fatal error!") 373 end if 374 write(std_out,*) ' The name of your file is: ',trim(filnam) 375 end if 376 377 select case(itask) 378 379 case(1) ! point calculation 380 call cut3d_pointint(gridtt,gridmx,gridmy,gridmz,nr1,nr2,nr3,nspden,rprimd) 381 exit 382 383 case(2) ! line calculation 384 call cut3d_lineint(gridtt,gridmx,gridmy,gridmz,nr1,nr2,nr3,nspden,rprimd) 385 exit 386 387 case(3) ! plane calculation 388 call cut3d_planeint(gridtt,gridmx,gridmy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart) 389 exit 390 391 case(4) ! volume calculation 392 write(std_out,*) ' Enter volume calculation' 393 call cut3d_volumeint(gridtt,gridmx,gridmy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart) 394 exit 395 396 case(5) 397 ! Rewrite the data on a formatted file, just in one (or four) column(s) 398 if (open_file(filnam,message, newunit=unt, status='unknown', action="write") /= 0) then 399 ABI_ERROR(message) 400 end if 401 402 if(nspden==1)then 403 do i3=1,nr3 404 do i2=1,nr2 405 do i1=1,nr1 406 write(unt,'(4(es22.12))') grid(i1,i2,i3) 407 end do 408 end do 409 end do 410 else 411 do i3=1,nr3 412 do i2=1,nr2 413 do i1=1,nr1 414 write(unt,'(4(es22.12))') gridtt(i1,i2,i3), gridmx(i1,i2,i3), gridmy(i1,i2,i3), gridmz(i1,i2,i3) 415 end do 416 end do 417 end do 418 end if 419 close(unt) 420 exit 421 422 case(6) 423 ! Rewrite the data on a formatted file, 3D index + density 424 if (open_file(filnam,message, newunit=unt, status='unknown', action="write") /= 0) then 425 ABI_ERROR(message) 426 end if 427 428 if(nspden==1)then 429 write(unt,*)' i1 i2 i3 data ' 430 do i3=1,nr3 431 do i2=1,nr2 432 do i1=1,nr1 433 write(unt,'(3i6,4(es24.14))') i1,i2,i3,grid(i1,i2,i3) 434 end do 435 end do 436 end do 437 else 438 if(nspden==2)then 439 write(unt,*)' i1 i2 i3 non-spin-polarized spin up spin down difference ' 440 else if(nspden==4)then 441 write(unt,*)' i1 i2 i3 non-spin-polarized x y z ' 442 end if 443 do i3=1,nr3 444 do i2=1,nr2 445 do i1=1,nr1 446 write(unt,'(3i6,4(es24.14))') i1,i2,i3,gridtt(i1,i2,i3),gridmx(i1,i2,i3),gridmy(i1,i2,i3),gridmz(i1,i2,i3) 447 end do 448 end do 449 end do 450 end if ! nspden 451 close(unt) 452 exit 453 454 case(7) 455 if (open_file(filnam,message, newunit=unt, form='unformatted', action="write") /= 0) then 456 ABI_ERROR(message) 457 end if 458 459 xm=0 ; xp=rprimd(1,1)*Bohr_Ang 460 ym=0 ; yp=rprimd(2,2)*Bohr_Ang 461 zm=0 ; zp=rprimd(3,3)*Bohr_Ang 462 write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms' 463 write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp 464 write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3 465 write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3 466 write(unt) xm,xp,ym,yp,zm,zp,nr1,nr2,nr3 467 ABI_MALLOC(rhomacu,(nr1,nr2)) 468 do i3=1,nr3 469 do i2=1,nr2 470 do i1=1,nr1 471 rhomacu(i1,i2)=grid(i1,i2,i3) 472 end do 473 end do 474 write(unt) rhomacu(:,:) 475 end do 476 close(unt) 477 exit 478 479 case (8) 480 if (open_file(filnam, message, newunit=unt, form='formatted', action="write") /= 0) then 481 ABI_ERROR(message) 482 end if 483 484 write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms' 485 write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp 486 write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3 487 write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3 488 write(unt,'(a)') 'TITLE = " " ' 489 write(unt,'(a)') 'VARIABLES = "X" "Y" "Z" (all three in Angstrom) "DENSITY or POTENTIAL" (atomic units) ' 490 write(unt,'(3(a,i6),a)') 'ZONE I=',nr1, ' J=', nr2, ' K=', nr3, ' F=POINT' 491 do i3=1,nr3 492 do i2=1,nr2 493 do i1=1,nr1 494 xnow = rprimd(1,1)*(i1-1)/nr1 + rprimd(1,2)*(i2-1)/nr2 + rprimd(1,3)*(i3-1)/nr3 495 ynow = rprimd(2,1)*(i1-1)/nr1 + rprimd(2,2)*(i2-1)/nr2 + rprimd(2,3)*(i3-1)/nr3 496 znow = rprimd(3,1)*(i1-1)/nr1 + rprimd(3,2)*(i2-1)/nr2 + rprimd(3,3)*(i3-1)/nr3 497 write(unt,'(4es22.15)') Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow, grid (i1,i2,i3) 498 end do 499 end do 500 end do 501 close(unt) 502 exit 503 504 case (9) 505 if (open_file(filnam, message, newunit=unt, form='formatted', action="write") /= 0) then 506 ABI_ERROR(message) 507 end if 508 xm=0 ; xp=rprimd(1,1)*Bohr_Ang 509 ym=0 ; yp=rprimd(2,2)*Bohr_Ang 510 zm=0 ; zp=rprimd(3,3)*Bohr_Ang 511 write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms' 512 write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp 513 write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1 514 write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1) 515 write(std_out,*) ' znucl = ', hdr%znucltypat, ' type = ', hdr%typat, ' ntypat = ', ntypat 516 517 gridshift1 = 0 518 gridshift2 = 0 519 gridshift3 = 0 520 write(std_out,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?' 521 write(std_out,*) 522 shift_tau(:) = zero 523 read(std_in,"(a)") outputchar 524 if (outputchar == 'y' .or. outputchar == 'Y') then 525 write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,'):' 526 write(std_out,*) 527 read(std_in,*) gridshift1, gridshift2, gridshift3 528 shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1) 529 end if 530 ! 531 ! Generate translated coordinates to match density shift 532 ! 533 ABI_MALLOC(tau2,(3,natom)) 534 do iatom = 1,natom 535 tau2(:,iatom) = xcart(:,iatom) - shift_tau(:) 536 end do 537 ! ################################################################### (LD) 538 ! Option only available for "xcrysden" format as documented at the beginning 539 if (ispden==4) then 540 ! It is necessary to know previously how many atoms will be used. 541 ! in order to plot the necessary magnetization arrows only. 542 write(std_out,*)'Is it possible to decrease the number of arrows in order to improve the' 543 write(std_out,*)'visualization in the screen, and decrease the size of the xcrysden output file.' 544 write(std_out,*)'How many arrows would you like to skip? 0 = take all. 1 = skip every other point...' 545 read (std_in,*) nrws 546 nrws=nrws+1 547 index=natom 548 maxmz=0.0 549 do i1=1,nr1,nrws 550 do i2=1,nr2,nrws 551 do i3=1,nr3,nrws 552 normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2 553 if(normz > maxmz) maxmz=normz 554 end do 555 end do 556 end do 557 if(abs(maxmz)<tol10)then 558 ABI_ERROR('At least, one of the components must differ from zero.') 559 end if 560 do i1=1,nr1,nrws 561 do i2=1,nr2,nrws 562 do i3=1,nr3,nrws 563 normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2 564 if(0.1*maxmz <= normz) index=index+1 565 end do 566 end do 567 end do 568 569 write(unt,'(1X,A)') 'CRYSTAL' 570 write(unt,'(1X,A)') 'PRIMVEC' 571 do i1 = 1,3 572 write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3) 573 end do 574 write(unt,'(1X,A)') 'PRIMCOORD' 575 write(unt,*) index, '1' 576 577 ! write out atom types and positions 578 do iatom = 1,natom 579 write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom) 580 end do 581 582 ! write out magnetization vectors. 583 ! xcrysden consider these as X (dummy) atoms. 584 do i1=1,nr1,nrws 585 do i2=1,nr2,nrws 586 do i3=1,nr3,nrws 587 normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2 588 if(0.1*maxmz <= normz) then 589 xcart2 = matmul (rprimd, (/(i1-one)/nr1, (i2-one)/nr2, (i3-one)/nr3/)) 590 write(unt,'(A,1X,6(ES17.10,2X))')'X',& 591 Bohr_Ang*(xcart2(1)-shift_tau(1)),& 592 Bohr_Ang*(xcart2(2)-shift_tau(2)),& 593 Bohr_Ang*(xcart2(3)-shift_tau(3)),& 594 gridmx(i1,i2,i3),& 595 gridmy(i1,i2,i3),& 596 gridmz(i1,i2,i3) 597 end if 598 end do 599 end do 600 end do 601 else 602 ! ################################################################### (LD) 603 ! 604 ! normal case: output density or potential (scalar field) 605 write(unt,'(1X,A)') 'DIM-GROUP' 606 write(unt,*) '3 1' 607 write(unt,'(1X,A)') 'PRIMVEC' 608 do i1 = 1,3 609 write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3) 610 end do 611 write(unt,'(1X,A)') 'PRIMCOORD' 612 write(unt,*) natom, ' 1' 613 do iatom = 1,natom 614 write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom) 615 end do 616 write(unt,'(1X,A)') 'ATOMS' 617 do iatom = 1,natom 618 write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom) 619 end do 620 ! write(31,'(1X,A)') 'FRAMES' 621 write(unt,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D' 622 write(unt,*) 'datagrids' 623 write(unt,'(1X,A)') 'DATAGRID_3D_DENSITY' 624 write(unt,*) nr1+1,nr2+1,nr3+1 625 write(unt,*) '0.0 0.0 0.0 ' 626 do i1 = 1,3 627 write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3) 628 end do 629 630 index = 0 631 do ir3=gridshift3+1,nr3+1 632 ii3=mod(ir3-1,nr3) + 1 633 do ir2=gridshift2+1,nr2+1 634 ii2=mod(ir2-1,nr2) + 1 635 do ir1=gridshift1+1,nr1+1 636 ii1=mod(ir1-1,nr1) + 1 637 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 638 index = index+1 639 if (mod (index,6) == 0) write (unt,*) 640 end do 641 do ir1=1,gridshift1 642 ii1=mod(ir1-1,nr1) + 1 643 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 644 index = index+1 645 if (mod (index,6) == 0) write (unt,*) 646 end do 647 end do 648 do ir2=1,gridshift2 649 ii2=mod(ir2-1,nr2) + 1 650 do ir1=gridshift1+1,nr1+1 651 ii1=mod(ir1-1,nr1) + 1 652 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 653 index = index+1 654 if (mod (index,6) == 0) write (unt,*) 655 end do 656 do ir1=1,gridshift1 657 ii1=mod(ir1-1,nr1) + 1 658 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 659 index = index+1 660 if (mod (index,6) == 0) write (unt,*) 661 end do 662 end do 663 end do 664 do ir3=1,gridshift3 665 ii3=mod(ir3-1,nr3) + 1 666 do ir2=gridshift2+1,nr2+1 667 ii2=mod(ir2-1,nr2) + 1 668 do ir1=gridshift1+1,nr1+1 669 ii1=mod(ir1-1,nr1) + 1 670 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 671 index = index+1 672 if (mod (index,6) == 0) write (unt,*) 673 end do 674 do ir1=1,gridshift1 675 ii1=mod(ir1-1,nr1) + 1 676 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 677 index = index+1 678 if (mod (index,6) == 0) write (unt,*) 679 end do 680 end do 681 do ir2=1,gridshift2 682 ii2=mod(ir2-1,nr2) + 1 683 do ir1=gridshift1+1,nr1+1 684 ii1=mod(ir1-1,nr1) + 1 685 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 686 index = index+1 687 if (mod (index,6) == 0) write (unt,*) 688 end do 689 do ir1=1,gridshift1 690 ii1=mod(ir1-1,nr1) + 1 691 write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3) 692 index = index+1 693 if (mod (index,6) == 0) write (unt,*) 694 end do 695 end do 696 end do 697 write (unt,*) 698 write(unt,'(1X,A)') 'END_DATAGRID_3D' 699 write(unt,'(1X,A)') 'END_BLOCK_DATAGRID3D' 700 701 end if 702 703 close(unt) 704 exit 705 706 case(11) 707 call cut3d_hirsh(grid,natom,nr1,nr2,nr3,ntypat,rprimd,xcart,hdr%typat,hdr%zionpsp,hdr%znucltypat) 708 exit 709 710 case(14) ! CUBE file format from GAUSSIAN 711 write(std_out,*) 712 write(std_out,*) 'Output a cube file of 3D volumetric data' 713 write(std_out,*) 714 715 ! EXAMPLE FROM THE WEB 716 ! CPMD CUBE FILE. 717 ! OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z 718 ! 3 0.000000 0.000000 0.000000 719 ! 40 0.283459 0.000000 0.000000 720 ! 40 0.000000 0.283459 0.000000 721 ! 40 0.000000 0.000000 0.283459 722 ! 8 0.000000 5.570575 5.669178 5.593517 723 ! 1 0.000000 5.562867 5.669178 7.428055 724 ! 1 0.000000 7.340606 5.669178 5.111259 725 ! -0.25568E-04 0.59213E-05 0.81068E-05 0.10868E-04 0.11313E-04 0.35999E-05 726 727 if (open_file(filnam,message,newunit=unt, status='unknown', form='formatted', action="write") /= 0) then 728 ABI_ERROR(message) 729 end if 730 731 !%% call print_fofr_cube(nr1,nr2,n3,nr1,nr2,nr3,fofr,rprimd,natom,znucl_atom,xcart,unit=unt) 732 write(unt,'(a)') 'ABINIT generated cube file' 733 write(unt,'(a)') 'from cut3d tool' 734 735 write(unt,'(i9,3(1x,f12.6))') natom,0.,0.,0. 736 write(unt,'(i9,3(1x,f12.6))') nr1,(rprimd(ir2,1)/nr1, ir2=1,3) 737 write(unt,'(i9,3(1x,f12.6))') nr2,(rprimd(ir2,2)/nr2, ir2=1,3) 738 write(unt,'(i9,3(1x,f12.6))') nr3,(rprimd(ir2,3)/nr3, ir2=1,3) 739 740 do iatom=1,natom 741 write(unt,'(i9,4(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),0.d0, & 742 & xcart(1,iatom),xcart(2,iatom),xcart(3,iatom) 743 end do 744 745 ! C ordering of the indexes 746 do i1=1,nr1 747 do i2=1,nr2 748 do i3=1,nr3 749 write(unt,'(6(f12.6,2x))') grid(i1,i2,i3) 750 end do 751 end do 752 end do 753 754 close(unt) 755 exit 756 757 case (15) 758 ! Write netcdf file. 759 cryst = hdr%get_crystal() 760 call ngfft_seq(ngfft, [nr1, nr2, nr3]) 761 ngfft(4:6) = ngfft(1:3) 762 nfft = product(ngfft(1:3)) 763 cplex = 1 764 call init_distribfft_seq(mpi_enreg%distribfft, 'c', ngfft(2), ngfft(3), 'all') 765 call init_distribfft_seq(mpi_enreg%distribfft, 'f', ngfft(2), ngfft(3), 'all') 766 767 call fftdatar_write(varname,filnam,IO_MODE_ETSF,hdr,cryst,ngfft,cplex,nfft,nspden,grid_full,mpi_enreg) 768 call cryst%free() 769 770 case(0) 771 write(std_out,*)' Exit requested by user' 772 exit 773 774 case default 775 ABI_ERROR(sjoin("Wrong task:", itoa(itask))) 776 end select 777 end do 778 779 write(std_out,*) ' Task ',itask,' has been done !' 780 write(std_out,*) 781 write(std_out,'(a)') ' More analysis of the 3D file ? ( 0=no ; 1=default=yes ; 2= treat another file - restricted usage)' 782 read(std_in,*) iprompt 783 if(iprompt/=1) then 784 call hdr%free() 785 exit 786 else 787 cycle 788 end if 789 end do 790 791 else 792 ABI_ERROR(sjoin("Don't know how to handle file class ", abifile%class)) 793 end if ! WF file or DEN/POT file 794 795 ! A maximum number of files had been previously specified, but set the actual number of files 796 ! to 1 if one does not read at least one other. 797 if(ifiles==1)then 798 nfiles=1 799 if(iprompt==2)nfiles=mfiles 800 801 ! A data structure for storing the important information should be created ... 802 ! Here, one supposes that the files are compatible ... 803 if(isdenpot(ifiles)==1)then 804 ABI_MALLOC(grid_full_stored,(nr1,nr2,nr3,nspden,nfiles)) 805 nr1_stored=nr1 806 nr2_stored=nr2 807 nr3_stored=nr3 808 nspden_stored=nspden 809 else if(isdenpot(ifiles)/=1 .and. iprompt==2)then 810 ABI_ERROR("in case of storage mode, the first file must be a density/potential file.") 811 end if 812 end if 813 814 if(isdenpot(ifiles)==1) grid_full_stored(:,:,:,:,ifiles)=grid_full(:,:,:,:) 815 if(isdenpot(ifiles)==1) filrho_stored(ifiles)=filrho 816 817 if(allocated(xcart)) then 818 ABI_FREE(xcart) 819 end if 820 if(allocated(xred)) then 821 ABI_FREE(xred) 822 end if 823 if(allocated(grid)) then 824 ABI_FREE(grid) 825 end if 826 if(allocated(grid_full)) then 827 ABI_FREE(grid_full) 828 end if 829 if(allocated(gridtt)) then 830 ABI_FREE(gridtt) 831 end if 832 if(allocated(gridmx)) then 833 ABI_FREE(gridmx) 834 end if 835 if(allocated(gridmy)) then 836 ABI_FREE(gridmy) 837 end if 838 if(allocated(gridmz)) then 839 ABI_FREE(gridmz) 840 end if 841 if(allocated(rhomacu)) then 842 ABI_FREE(rhomacu) 843 end if 844 if(allocated(tau2)) then 845 ABI_FREE(tau2) 846 end if 847 848 if(iprompt/=2) then 849 exit 850 end if 851 852 end do ! End big loop on files 853 854 !Will provide different information on the density and potential files 855 do ifiles=1,nfiles 856 if(isdenpot(ifiles)==1)then 857 write(std_out,*) 858 write(std_out,*) ' Provide some global information about the density and/or potential file(s)' 859 exit 860 end if 861 end do 862 do ifiles=1,nfiles 863 if(isdenpot(ifiles)==1)then 864 write(std_out,*) 865 write(std_out, '(a,i5,3a)' ) '- File number ',ifiles,', with name "',trim(filrho_stored(ifiles)),'"' 866 write(std_out, '(a,i12,a,es14.6)' ) ' Number of grid points =',nr1*nr2*nr3,' ; Volume of real space cell (Bohr^3)=',ucvol 867 do ispden=1,nspden 868 sumdenpot=sum(grid_full_stored(:,:,:,ispden,ifiles)) 869 write(std_out, '(a,i5,3a)' ) ' Spin-component number ',ispden 870 write(std_out, '(a,3es16.6)' ) ' Sum of values, mean, mean times cell volume=',& 871 & sumdenpot,sumdenpot/real(nr1*nr2*nr3),sumdenpot*ucvol/real(nr1*nr2*nr3) 872 end do 873 end if 874 end do 875 876 if(nspden==1)then 877 ! At present, only nspden=1 is correctly implemented, due to specificities of the treatment of the spin-density 878 do ifiles=1,nfiles 879 if(isdenpot(ifiles)==1)then 880 write(std_out,*) 881 write(std_out,'(a)') ' Provide some global joint information about the stored density and potential file(s)' 882 exit 883 end if 884 end do 885 do ifiles=1,nfiles 886 if(isdenpot(ifiles)==1)then 887 do jfiles=ifiles,nfiles 888 if(isdenpot(jfiles)==1)then 889 write(std_out,*) 890 write(std_out, '(a,2i5)' )' File numbers: ',ifiles,jfiles 891 do ispden=1,nspden 892 dotdenpot=zero 893 do ir1=1,nr1 894 do ir2=1,nr2 895 do ir3=1,nr3 896 dotdenpot=dotdenpot+grid_full_stored(ir1,ir2,ir3,ispden,ifiles)*grid_full_stored(ir1,ir2,ir3,ispden,jfiles) 897 end do 898 end do 899 end do 900 write(std_out, '(a,i5,3a)' ) ' Spin-component number ',ispden 901 write(std_out, '(a,3es16.6)' ) ' Dot product of values, mean, mean times cell volume=',& 902 ! write(std_out, '(a,3es20.10)' ) ' Dot product of values, mean, mean times cell volume=',& 903 & dotdenpot,dotdenpot/real(nr1*nr2*nr3),dotdenpot*ucvol/real(nr1*nr2*nr3) 904 end do 905 end if 906 end do 907 end if 908 end do 909 end if 910 911 ABI_FREE(filrho_stored) 912 913 if(allocated(grid_full_stored)) then 914 ABI_FREE(grid_full_stored) 915 end if 916 917 if(allocated(isdenpot)) then 918 ABI_FREE(isdenpot) 919 end if 920 921 call timein(tsec(1),tsec(2)) 922 tsec(1)=tsec(1)-tcpui 923 tsec(2)=tsec(2)-twalli 924 925 write(std_out, '(3a,f13.1,a,f13.1)' )'-',ch10,'- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 926 927 write(std_out,*) 928 write(std_out,*) ' Thank you for using me' 929 write(std_out,*) 930 931 call flush_unit(std_out) 932 call destroy_mpi_enreg(mpi_enreg) 933 call abinit_doctor("__cut3d") 934 call xmpi_end() 935 936 end program cut3d