TABLE OF CONTENTS
- ABINIT/m_cut3d
- m_cut3d/cut3d_hirsh
- m_cut3d/cut3d_lineint
- m_cut3d/cut3d_planeint
- m_cut3d/cut3d_pointint
- m_cut3d/cut3d_rrho
- m_cut3d/cut3d_volumeint
- m_cut3d/cut3d_wffile
- m_cut3d/normalize
- m_cut3d/reduce
- m_cut3d/vdot
ABINIT/m_cut3d [ Modules ]
NAME
m_cut3d
FUNCTION
This module the predures used by cut3d
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (XG,MVerstraete,GMR,RC,LSI,JFB,MCote,MB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_cut3d 24 25 use defs_basis 26 use defs_abitypes 27 use m_abicore 28 use m_errors 29 use m_splines 30 use m_hdr 31 #ifdef HAVE_NETCDF 32 use netcdf 33 #endif 34 use m_nctk 35 use m_wfk 36 use m_xmpi 37 use m_sort 38 39 use m_io_tools, only : get_unit, iomode_from_fname, open_file, file_exists, read_string 40 use m_numeric_tools, only : interpol3d 41 use m_symtk, only : matr3inv 42 use m_fstrings, only : int2char10, sjoin, itoa 43 use m_geometry, only : xcart2xred, metric 44 use m_special_funcs, only : jlspline_t, jlspline_new, jlspline_free, jlspline_integral 45 use m_pptools, only : print_fofr_ri, print_fofr_xyzri , print_fofr_cube 46 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 47 use m_cgtools, only : cg_getspin 48 use m_gsphere, only : getkpgnorm 49 use m_epjdos, only : recip_ylm, dens_in_sph 50 use m_dens, only : dens_hirsh 51 use m_kg, only : kpgio, ph1d3d, getph 52 use m_fftcore, only : sphereboundary 53 use m_initylmg, only : initylmg 54 use m_fft, only : fourwf 55 56 implicit none 57 58 private 59 60 public :: cut3d_hirsh 61 public :: cut3d_rrho 62 public :: cut3d_volumeint 63 public :: cut3d_planeint 64 public :: cut3d_lineint 65 public :: cut3d_pointint 66 public :: cut3d_wffile 67 68 CONTAINS !===========================================================
m_cut3d/cut3d_hirsh [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_hirsh
FUNCTION
Compute the Hirshfeld charges
INPUTS
grid_den(nrx,nry,nrz)= density on the grid natom = number of atoms in the unit cell nrx,nry,nrz= number of points in the grid for the three directions ntypat=number of types of atoms in unit cell. rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(natom)=type of each atom xcart(3,natom) = different positions of the atoms in the unit cell zion=(ntypat)gives the ionic charge for each type of atom znucl(ntypat)=gives the nuclear number for each type of atom
OUTPUT
write the Hirshfeld charge decomposition
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
104 subroutine cut3d_hirsh(grid_den,natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,znucl) 105 106 107 !This section has been created automatically by the script Abilint (TD). 108 !Do not modify the following lines by hand. 109 #undef ABI_FUNC 110 #define ABI_FUNC 'cut3d_hirsh' 111 !End of the abilint section 112 113 implicit none 114 115 !Arguments ------------------------------------ 116 !scalars 117 integer,intent(in) :: natom,nrx,nry,nrz,ntypat 118 !arrays 119 integer,intent(in) :: typat(natom) 120 real(dp),intent(in) :: grid_den(nrx,nry,nrz),rprimd(3,3),zion(ntypat) 121 real(dp),intent(in) :: znucl(ntypat) 122 real(dp),intent(in) :: xcart(3,natom) 123 124 !Local variables ------------------------- 125 !scalars 126 integer,parameter :: prtcharge1=1 127 integer :: ierr,ipoint,itypat,mpoint,temp_unit 128 real(dp) :: minimal_den 129 real(dp) :: param1,param2,xx,yy 130 character(len=fnlen) :: file_allelectron 131 character(len=500) :: msg 132 !arrays 133 integer,allocatable :: npoint(:) 134 real(dp),allocatable :: aeden(:,:),hcharge(:),hden(:),hweight(:),radii(:,:) 135 136 ! ********************************************************************* 137 138 !1. Read the 1D all-electron atomic files 139 !Store the radii in radii(:,itypat), and the all-electron 140 !densities in aeden(:,itypat). The number of the last 141 !point with significant density is stored in npoint(itypat) 142 143 minimal_den=tol6 144 mpoint=4000 145 ABI_ALLOCATE(npoint,(ntypat)) 146 ABI_ALLOCATE(radii,(4000,ntypat)) 147 ABI_ALLOCATE(aeden,(4000,ntypat)) 148 do itypat=1,ntypat 149 write(std_out,'(a)' )' Please, give the filename of the all-electron density file' 150 write(std_out,'(a,es16.6)' )' for the first type of atom, with atomic number=',znucl(itypat) 151 if (read_string(file_allelectron, unit=std_in) /= 0) then 152 MSG_ERROR("Fatal error!") 153 end if 154 write(std_out,*)' The name you entered is : ',trim(file_allelectron),ch10 155 ierr = open_file(file_allelectron,msg,newunit=temp_unit,form='formatted',status='old') 156 if (ierr/=0) then 157 MSG_ERROR(msg) 158 else 159 read(temp_unit, *) param1, param2 160 do ipoint=1,mpoint 161 ! Either the file is finished 162 read(temp_unit, *, end=888) xx,yy 163 radii(ipoint,itypat)=xx 164 aeden(ipoint,itypat)=yy 165 ! Or the density is lower than the minimal significant value 166 if(yy<minimal_den)exit 167 end do 168 888 continue 169 npoint(itypat)=ipoint-1 170 if(ipoint==mpoint)then 171 write(std_out,*)' hirsh : mpoint is too low, increase its value to match ipoint.' 172 end if 173 end if 174 close(temp_unit) 175 end do 176 177 ABI_MALLOC(hden,(natom)) 178 ABI_MALLOC(hcharge,(natom)) 179 ABI_MALLOC(hweight,(natom)) 180 181 call dens_hirsh(mpoint,radii,aeden,npoint,minimal_den,grid_den, & 182 natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,prtcharge1,hcharge,hden,hweight) 183 184 ABI_FREE(hweight) 185 ABI_FREE(aeden) 186 ABI_FREE(hcharge) 187 ABI_FREE(hden) 188 ABI_FREE(npoint) 189 ABI_FREE(radii) 190 191 end subroutine cut3d_hirsh
m_cut3d/cut3d_lineint [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_lineint
FUNCTION
Computes the values along a line defined by two points
INPUTS
gridtt(nr1,nr2,nr3)=Total density gridux(nr1,nr2,nr3)=spin-Up density, or magnetization density in X direction griddy(nr1,nr2,nr3)=spin-Down density, or magnetization density in Y direction gridmz(nr1,nr2,nr3)=spin-polarization density or magnetization density in Z direction nr1=grid size along x nr2=grid size along y nr3=grid size along z nspden=number of spin-density components rprimd(3,3)=orientation of the unit cell in 3D
OUTPUT
only writing
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
227 subroutine cut3d_lineint(gridtt,gridux,griddy,gridmz,nr1,nr2,nr3,nspden,rprimd) 228 229 230 !This section has been created automatically by the script Abilint (TD). 231 !Do not modify the following lines by hand. 232 #undef ABI_FUNC 233 #define ABI_FUNC 'cut3d_lineint' 234 !End of the abilint section 235 236 implicit none 237 238 !Arguments------------------------------------------------------------- 239 !scalars 240 integer,intent(in) :: nr1,nr2,nr3,nspden 241 !arrays 242 real(dp),intent(in) :: griddy(nr1,nr2,nr3) 243 real(dp),intent(in) :: gridmz(nr1,nr2,nr3),gridtt(nr1,nr2,nr3) 244 real(dp),intent(in) :: gridux(nr1,nr2,nr3),rprimd(3,3) 245 246 !Local variables-------------------------------------------------------- 247 !scalars 248 integer :: inpopt,inpopt2,k2,nresol,okline,unt 249 real(dp) :: denvaldy,denvalmz,denvaltt,denvalux,dx,dy,dz,length 250 character(len=fnlen) :: filnam 251 character(len=500) :: msg 252 !arrays 253 real(dp) :: cent(3),r1(3),r2(3),rcart(3),rr(3),x1(3),x2(3) 254 255 ! ********************************************************************* 256 257 okline=0 258 do while (okline==0) 259 write(std_out,*) ' Type 1) for a line between two cartesian-defined points' 260 write(std_out,*) ' or 2) for a line between two crystallographic-defined points ' 261 write(std_out,*) ' or 3) for a line defined by its direction in cartesion coordinates' 262 write(std_out,*) ' or 4) for a line defined by its direction in crystallographic coordinates' 263 read(std_in,*) inpopt 264 write(std_out,*) ' You typed ',inpopt,ch10 265 if (inpopt==1 .or. inpopt ==2 .or. inpopt==3 .or. inpopt==4) okline=1 266 end do 267 268 !In the case of a line defined by its two extreme points 269 if (inpopt==1) then 270 write(std_out,*) ' Type the first point coordinates (Bohrs):' 271 write(std_out,*) ' -> X-dir Y-dir Z-dir:' 272 read(std_in,*) x1 273 write(std_out,'(a,3es16.6,a)') ' You typed ',x1,ch10 274 call reduce(r1,x1,rprimd) 275 276 write(std_out,*) ' Type the second point coordinates (Bohrs):' 277 write(std_out,*) ' -> X-dir Y-dir Z-dir:' 278 read(std_in,*) x2 279 write(std_out,'(a,3es16.6,a)') ' You typed ',x2,ch10 280 call reduce(r2,x2,rprimd) 281 end if 282 283 if (inpopt==2) then 284 write(std_out,*) ' Type the first point coordinates (fractional):' 285 write(std_out,*) ' -> X-dir Y-dir Z-dir:' 286 read(std_in,*) r1 287 write(std_out,'(a,3es16.6,a)') ' You typed ',r1,ch10 288 289 write(std_out,*) ' Type the second point coordinates (fractional):' 290 write(std_out,*) ' -> X-dir Y-dir Z-dir:' 291 read(std_in,*) r2 292 write(std_out,'(a,3es16.6,a)') ' You typed ',r2,ch10 293 end if 294 295 if(inpopt==3 .or. inpopt==4 )then 296 297 write(std_out,*) 'Please enter now the line direction:' 298 write(std_out,*) ' -> X-dir Y-dir Z-dir:' 299 read(std_in,*) x2 300 write(std_out,'(a,3es16.6,a)') 'The line direction is:',x2(1),x2(2),x2(3),ch10 301 302 if (inpopt == 4) then 303 rcart=matmul(x2,rprimd) 304 x2(:)=rcart(:) 305 write(std_out,'(a,3es16.6,a)') 'Expressed in cartesian coordinates: ',x2(1),x2(2),x2(3),ch10 306 end if 307 308 call normalize(x2) 309 310 write(std_out,*) 'Enter now the central point of line:' 311 write(std_out,*) 'Type 1) for cartesian coordinates' 312 write(std_out,*) ' or 2) for crystallographic coordinates' 313 read(std_in,*) inpopt2 314 if (inpopt2==1 .or. inpopt2==2) then 315 write(std_out,*) 'Type the point coordinates:' 316 write(std_out,*) ' -> X-Coord Y-Coord Z-Coord:' 317 read(std_in,*) cent 318 write(std_out,'(a,3es16.6,a)') 'Central point coordinates:', cent(1),cent(2),cent(3),ch10 319 if (inpopt2==2) then 320 rcart=matmul(cent,rprimd) 321 cent(:)=rcart(:) 322 write(std_out,'(a,3es16.6,a)') 'Expressed in cartesian coordinates:',cent(1),cent(2),cent(3),ch10 323 end if 324 write(std_out,*) 'Enter line length (in cartesian coordinates, in Bohr):' 325 read(std_in,*) length 326 327 ! Compute the extremal points in cartesian coordinates 328 x1(:)=cent(:)-length*x2(:)*half 329 x2(:)=cent(:)+length*x2(:)*half 330 331 ! Transfer to crystallographic coordinates 332 call reduce(r1,x1,rprimd) 333 call reduce(r2,x2,rprimd) 334 335 end if 336 337 end if ! inpopt 338 339 write(std_out,*) 340 write(std_out,'(a,3es16.6)' ) ' Crystallographic coordinates of the first point :',r1 341 write(std_out,'(a,3es16.6)' ) ' Crystallographic coordinates of the second point :',r2 342 write(std_out,*) 343 344 write(std_out,*) ' Enter line resolution: (integer, number of points on the line)' 345 read(std_in,*) nresol 346 write(std_out,*) ' You typed',nresol,ch10 347 348 !At this moment the code knows everything about the geometric input, the data and 349 !the line direction. It will further calculate the values along this line using 350 !an interpolation 351 352 write(std_out,*) ch10,' Enter the name of an output file:' 353 if (read_string(filnam, unit=std_in) /= 0) then 354 MSG_ERROR("Fatal error!") 355 end if 356 write(std_out,*) ' The name of your file is : ',trim(filnam),ch10 357 358 if (open_file(filnam,msg,newunit=unt,status='unknown') /= 0) then 359 MSG_ERROR(msg) 360 end if 361 362 dx=(r2(1)-r1(1))/nresol 363 dy=(r2(2)-r1(2))/nresol 364 dz=(r2(3)-r1(3))/nresol 365 366 !DEBUG 367 !write(std_out,*)' nspden=',nspden 368 !ENDDEBUG 369 370 if(nspden==1)then 371 write(std_out,*)' Index of point value ' 372 else if (nspden==2)then 373 write(std_out,*)' Index of point non-spin-polarized spin up spin down difference ' 374 else if (nspden==4)then 375 write(std_out,*)' Index of point non-spin-polarized x y z ' 376 end if 377 378 do k2=0,nresol 379 380 rr(1)=r1(1)+k2*dx 381 rr(2)=r1(2)+k2*dy 382 rr(3)=r1(3)+k2*dz 383 384 rr(1)=mod(mod(rr(1),1._dp)+1._dp,1._dp) 385 rr(2)=mod(mod(rr(2),1._dp)+1._dp,1._dp) 386 rr(3)=mod(mod(rr(3),1._dp)+1._dp,1._dp) 387 388 denvaltt = interpol3d(rr,nr1,nr2,nr3,gridtt) 389 if(nspden==1)then 390 write(unt, '(i13,es22.12)' ) k2,denvaltt 391 write(std_out,'(i13,es22.12)' ) k2,denvaltt 392 393 else if(nspden==2 .or. nspden==4)then 394 denvalux = interpol3d(rr,nr1,nr2,nr3,gridux) 395 denvaldy = interpol3d(rr,nr1,nr2,nr3,griddy) 396 denvalmz = interpol3d(rr,nr1,nr2,nr3,gridmz) 397 write(unt, '(i13,4(es22.12))' ) k2,denvaltt,denvalux,denvaldy,denvalmz 398 write(std_out,'(i13,4es22.12)' ) k2,denvaltt,denvalux,denvaldy,denvalmz 399 end if 400 end do 401 402 close(unt) 403 404 end subroutine cut3d_lineint
m_cut3d/cut3d_planeint [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_planeint
FUNCTION
Computes the values within a plane
INPUTS
gridtt(nr1,nr2,nr3)=Total density gridux(nr1,nr2,nr3)=spin-Up density, or magnetization density in X direction griddy(nr1,nr2,nr3)=spin-Down density, or magnetization density in Y direction gridmz(nr1,nr2,nr3)=spin-polarization density or magnetization density in Z direction natom=integer number of atoms nr1=grid size along x nr2=grid size along y nr3=grid size along z nspden=number of spin-density components rprimd(3,3)=orientation of the unit cell in 3D tau(3,nat)=atomic positions in 3D cartesian space (from XMOL format)
OUTPUT
only writing
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
505 subroutine cut3d_planeint(gridtt,gridux,griddy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,tau) 506 507 508 !This section has been created automatically by the script Abilint (TD). 509 !Do not modify the following lines by hand. 510 #undef ABI_FUNC 511 #define ABI_FUNC 'cut3d_planeint' 512 !End of the abilint section 513 514 implicit none 515 516 !Arguments ------------------------------------ 517 !scalars 518 integer,intent(in) :: natom,nr1,nr2,nr3,nspden 519 !arrays 520 real(dp),intent(in) :: griddy(nr1,nr2,nr3) 521 real(dp),intent(in) :: gridmz(nr1,nr2,nr3),gridtt(nr1,nr2,nr3) 522 real(dp),intent(in) :: gridux(nr1,nr2,nr3),rprimd(3,3),tau(3,natom) 523 524 !Local variables ------------------------- 525 !scalars 526 integer :: iat,idir,ii,inpopt,itypat,k2,k3,mu,nresoll,nresolw,okhkl,okinp 527 integer :: okparam,oksure,unt 528 real(dp) :: denvaldy,denvalmz,denvaltt,denvalux,length 529 real(dp) :: width,xcoord,ycoord 530 character(len=fnlen) :: filnam 531 character(len=500) :: msg 532 !arrays 533 integer :: hkl(3) 534 real(dp) :: cent(3),mminv(3,3),r1(3),r2(3),r3(3),rcart(3),rr(3),x1(3) 535 real(dp) :: x2(3),x3(3),xcart(3) 536 537 ! ********************************************************************* 538 539 !Several lines to compute the transformation matrix from crystallographic to cartesian 540 541 call matr3inv(rprimd,mminv) 542 543 !Start of the real input of the plane orientation 544 545 okinp=0 546 do while (okinp==0) 547 write(std_out,*) 548 write(std_out,*) ' Type 1) for a plane passing through 3 atoms' 549 write(std_out,*) ' or 2) for a plane passing through 3 cartesian points' 550 write(std_out,*) ' or 3) for a plane passing through 3 crystallographic points' 551 write(std_out,*) ' or 4) for a plane parallel to a crystallographic plane' 552 write(std_out,*) ' or 5) for a plane orthogonal to a cartesian direction' 553 write(std_out,*) ' or 6) for a plane orthogonal to a crystallographic direction' 554 write(std_out,*) ' or 0) to stop' 555 read(std_in,*) itypat 556 select case (itypat) 557 558 case (0) 559 stop 560 561 ! A plane passing through 3 atoms 562 case (1) 563 write(std_out,*) ' The X axis will be through atms: 1,2 ' 564 write(std_out,*) ' Define each atom by its species and its number:' 565 write(std_out,*) ' -> atom 1 (iat):' 566 read(std_in,*) iat 567 x1(1)=tau(1,iat) 568 x1(2)=tau(2,iat) 569 x1(3)=tau(3,iat) 570 write(std_out,'(a,3f10.6)') ' position: ',x1 571 write(std_out,*) 572 write(std_out,*) ' -> atom 2 (iat):' 573 read(std_in,*) iat 574 x2(1)=tau(1,iat) 575 x2(2)=tau(2,iat) 576 x2(3)=tau(3,iat) 577 write(std_out,'(a,3f10.6)') ' position: ',x2 578 write(std_out,*) 579 write(std_out,*) ' -> atom 3 (iat):' 580 read(std_in,*) iat 581 x3(1)=tau(1,iat) 582 x3(2)=tau(2,iat) 583 x3(3)=tau(3,iat) 584 write(std_out,'(a,3f10.6)') ' position: ',x3 585 write(std_out,*) 586 587 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 588 do idir=1,3 589 x2(idir)=x2(idir)-x1(idir) 590 x3(idir)=x3(idir)-x1(idir) 591 end do 592 call normalize(x2) 593 call vdot(x3,x2,x1) 594 call normalize(x1) 595 call vdot(x2,x1,x3) 596 call normalize(x3) 597 okinp=1 598 599 ! A plane passing through 3 cartesian points 600 case (2) 601 write(std_out,*) ' The X axis will be through points: 1,2 ' 602 write(std_out,*) ' Define each :point coordinates' 603 write(std_out,*) ' -> point 1: X-coord Y-coord Z-coord:' 604 read(std_in,*) xcart 605 x1(:)=xcart(:) 606 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x1 607 write(std_out,*) 608 write(std_out,*) ' -> point 2: X-coord Y-coord Z-coord:' 609 read(std_in,*) xcart 610 x2(:)=xcart(:) 611 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x2 612 write(std_out,*) 613 write(std_out,*) ' -> point 3: X-coord Y-coord Z-coord:' 614 read(std_in,*) xcart 615 x3(:)=xcart(:) 616 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x3 617 write(std_out,*) 618 619 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 620 do idir=1,3 621 x2(idir)=x2(idir)-x1(idir) 622 x3(idir)=x3(idir)-x1(idir) 623 end do 624 call normalize(x2) 625 call vdot(x3,x2,x1) 626 call normalize(x1) 627 call vdot(x2,x1,x3) 628 call normalize(x3) 629 okinp=1 630 631 ! A plane passing through 3 crystallographic points 632 case (3) 633 write(std_out,*) ' The X axis will be through points: 1,2 ' 634 write(std_out,*) ' Define each :point coordinates' 635 write(std_out,*) ' -> point 1: X-coord Y-coord Z-coord:' 636 read(std_in,*) r1 637 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r1 638 write(std_out,*) 639 write(std_out,*) ' -> point 2: X-coord Y-coord Z-coord:' 640 read(std_in,*) r2 641 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r2 642 write(std_out,*) 643 write(std_out,*) ' -> point 3: X-coord Y-coord Z-coord:' 644 read(std_in,*) r3 645 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r3 646 write(std_out,*) 647 648 ! Transforms the points coordinates into cartesian 649 do mu=1,3 650 x1(mu)=rprimd(mu,1)*r1(1)+rprimd(mu,2)*r1(2)+rprimd(mu,3)*r1(3) 651 x2(mu)=rprimd(mu,1)*r2(1)+rprimd(mu,2)*r2(2)+rprimd(mu,3)*r2(3) 652 x3(mu)=rprimd(mu,1)*r3(1)+rprimd(mu,2)*r3(2)+rprimd(mu,3)*r3(3) 653 end do 654 655 write(std_out,*) ' Cartesian positions:' 656 write(std_out,*) x1 657 write(std_out,*) x2 658 write(std_out,*) x3 659 660 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 661 do idir=1,3 662 x2(idir)=x2(idir)-x1(idir) 663 x3(idir)=x3(idir)-x1(idir) 664 end do 665 call normalize(x2) 666 call vdot(x3,x2,x1) 667 call normalize(x1) 668 call vdot(x2,x1,x3) 669 call normalize(x3) 670 okinp=1 671 672 ! A plane parallel to a crystallographic plane 673 case (4) 674 okhkl=0 675 do while (okhkl==0) 676 write(std_out,*) ' Enter plane coordinates:' 677 write(std_out,*) ' -> H K L ' 678 read(std_in,*) hkl 679 if (.not. (hkl(1)==0 .and. hkl(2)==0 .and. hkl(3)==0)) okhkl=1 680 end do 681 write(std_out,*) ' Miller indices are:',hkl 682 683 do ii=1,3 684 x1(ii)=mminv(ii,1)*hkl(1) + mminv(ii,2)*hkl(2) + mminv(ii,3)*hkl(3) 685 end do 686 write(std_out,*) ' Orthogonal vector to the plane',x1 687 688 call normalize(x1) 689 if((x1(1).ne.0).or.(x1(2).ne.0)) then 690 x2(1)=-x1(2) 691 x2(2)=x1(1) 692 x2(3)=0 693 call normalize(x2) 694 else 695 x2(1)=1 696 x2(2)=0 697 x2(3)=0 698 end if 699 call vdot(x2,x1,x3) 700 call normalize(x3) 701 okinp=1 702 703 ! A plane orthogonal to a cartesian direction 704 case (5) 705 write(std_out,*) ' Enter the cartesian coordinates of the vector orthogonal to plane:' 706 write(std_out,*) ' -> X-dir Y-dir Z-dir (Angstroms or Bohrs):' 707 read(std_in,*) x1 708 call normalize(x1) 709 if((x1(1).ne.0).or.(x1(2).ne.0)) then 710 x2(1)=-x1(2) 711 x2(2)=x1(1) 712 x2(3)=0 713 call normalize(x2) 714 else 715 x2(1)=1 716 x2(2)=0 717 x2(3)=0 718 end if 719 call vdot(x2,x1,x3) 720 call normalize(x3) 721 okinp=1 722 723 ! A plane orthogonal to a crystallographic direction 724 case (6) 725 write(std_out,*) ' Enter crystallographic vector orthogonal to plane:' 726 write(std_out,*) ' -> X-dir Y-dir Z-dir (Fractional coordinates):' 727 read(std_in,*) r1 728 okinp=1 729 do mu=1,3 730 x1(mu)=rprimd(mu,1)*r1(1)+rprimd(mu,2)*r1(2)+rprimd(mu,3)*r1(3) 731 end do 732 call normalize(x1) 733 if((x1(1).ne.0).or.(x1(2).ne.0)) then 734 x2(1)=-x1(2) 735 x2(2)=x1(1) 736 x2(3)=0 737 call normalize(x2) 738 else 739 x2(1)=1 740 x2(2)=0 741 x2(3)=0 742 end if 743 call vdot(x2,x1,x3) 744 call normalize(x3) 745 okinp=1 746 747 case default 748 okinp=0 749 write(std_out,*) 'Input option do not correspond to the available options' 750 write(std_out,*) 'Please try again' 751 end select 752 753 end do 754 755 !At this moment the family of planes was defined 756 !The code knows also some of the geometric input 757 !It will proceed to the anchorage of the plane onto a point and then 758 !to the effective calculation 759 760 write(std_out,*) ' Vectors: (orthogonal & normalized) ' 761 write(std_out,'(11x,a,3f10.6)') ' X-dir in the plot ',x2 762 write(std_out,'(11x,a,3f10.6)') ' Y-dir in the plot ',x3 763 write(std_out,'(11x,a,3f10.6)') ' Z-dir (orth. to the plot) ',x1 764 765 write(std_out,*) 766 write(std_out,*) ' Enter central point of plane (Bohrs):' 767 write(std_out,*) ' Type 1) for Cartesian coordinates.' 768 write(std_out,*) ' or 2) for Crystallographic coordinates.' 769 read(std_in,*) inpopt 770 write(std_out,*) ' -> X-Coord Y-Coord Z-Coord:' 771 read(std_in,*) cent 772 773 if (inpopt==2) then 774 775 do mu=1,3 776 rcart(mu)=rprimd(mu,1)*cent(1)+rprimd(mu,2)*cent(2)+rprimd(mu,3)*cent(3) 777 end do 778 779 cent(:)=rcart(:) 780 write(std_out,'(a,3f16.6)' ) ' Expressed in cartesian coordinates: ',cent(1),cent(2),cent(3) 781 782 end if 783 784 okparam=0 785 do while(okparam==0) 786 write(std_out,*) 787 write(std_out,*) ' Enter plane width:' 788 read(std_in,*) width 789 write(std_out,*) ' Enter plane length:' 790 read(std_in,*) length 791 write(std_out,*) 792 write(std_out,*) ' Enter plane resolution in width:' 793 read(std_in,*) nresolw 794 write(std_out,*) ' Enter plane resolution in length:' 795 read(std_in,*) nresoll 796 write(std_out,*) ch10,' Enter the name of an output file:' 797 if (read_string(filnam, unit=std_in) /= 0) then 798 MSG_ERROR("Fatal error!") 799 end if 800 write(std_out,*) ' The name of your file is : ',trim(filnam) 801 write(std_out,*) 802 write(std_out,*) ' You asked for a plane of ',length,' x ',width 803 write(std_out,*) ' With a resolution of ',nresoll,' x ',nresolw 804 write(std_out,*) ' The result will be redirected to the file: ',trim(filnam) 805 write(std_out,*) ' These parameters may still be changed.' 806 write(std_out,*) ' Are you sure you want to keep them? (1=default=yes,2=no) ' 807 read(std_in,*) oksure 808 if (oksure/=2) okparam=1 809 end do 810 811 if (open_file(filnam,msg,newunit=unt,status='unknown') /= 0) then 812 MSG_ERROR(msg) 813 end if 814 815 do k2=-nresoll/2,nresoll/2 816 do k3=-nresolw/2,nresolw/2 817 rcart(1)=cent(1) + k2*x2(1)*length/nresoll + k3*x3(1)*width/nresolw 818 rcart(2)=cent(2) + k2*x2(2)*length/nresoll + k3*x3(2)*width/nresolw 819 rcart(3)=cent(3) + k2*x2(3)*length/nresoll + k3*x3(3)*width/nresolw 820 xcoord=k2*length/nresoll 821 ycoord=k3*width/nresolw 822 call reduce(rr,rcart,rprimd) 823 rr(1)=mod(mod(rr(1),1._dp)+1._dp,1._dp) 824 rr(2)=mod(mod(rr(2),1._dp)+1._dp,1._dp) 825 rr(3)=mod(mod(rr(3),1._dp)+1._dp,1._dp) 826 denvaltt = interpol3d(rr,nr1,nr2,nr3,gridtt) 827 if(nspden==2 .or. nspden==4)then 828 denvalux = interpol3d(rr,nr1,nr2,nr3,gridux) 829 denvaldy = interpol3d(rr,nr1,nr2,nr3,griddy) 830 denvalmz = interpol3d(rr,nr1,nr2,nr3,gridmz) 831 end if 832 if(nspden==1)then 833 write(unt, '(3e16.8)' ) xcoord,ycoord,denvaltt 834 else 835 write(unt, '(3e16.8)' ) xcoord,ycoord,denvaltt,denvalux,denvaldy,denvalmz 836 end if 837 end do 838 end do 839 840 close(unt) 841 842 end subroutine cut3d_planeint
m_cut3d/cut3d_pointint [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_pointint
FUNCTION
Computes the values at any point rr (this point is input from keyboard)
INPUTS
gridtt(nr1,nr2,nr3)=Total density gridux(nr1,nr2,nr3)=spin-Up density, or magnetization density in X direction griddy(nr1,nr2,nr3)=spin-Down density, or magnetization density in Y direction gridmz(nr1,nr2,nr3)=spin-polarization density or magnetization density in Z direction nr1=grid size along x nr2=grid size along y nr3=grid size along z nspden=number of spin-density components rprimd(3,3)=orientation of the unit cell in 3D
OUTPUT
only writing
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
878 subroutine cut3d_pointint(gridt,gridu,gridd,gridm,nr1,nr2,nr3,nspden,rprimd) 879 880 881 !This section has been created automatically by the script Abilint (TD). 882 !Do not modify the following lines by hand. 883 #undef ABI_FUNC 884 #define ABI_FUNC 'cut3d_pointint' 885 !End of the abilint section 886 887 implicit none 888 889 !Arguments-------------------------------------------------------------- 890 !scalars 891 integer,intent(in) :: nr1,nr2,nr3,nspden 892 !arrays 893 real(dp),intent(in) :: gridd(nr1,nr2,nr3),gridm(nr1,nr2,nr3) 894 real(dp),intent(in) :: gridt(nr1,nr2,nr3),gridu(nr1,nr2,nr3),rprimd(3,3) 895 896 !Local variables-------------------------------------------------------- 897 !scalars 898 integer :: inpopt,mu,okinp 899 real(dp) :: denvaldy,denvalmz,denvaltt,denvalux 900 !arrays 901 real(dp) :: rcart(3),rr(3) 902 903 ! ************************************************************************* 904 905 okinp=0 906 do while (okinp==0) 907 write(std_out,*) ' Select the coordinate system:' 908 write(std_out,*) ' Type 1) for cartesian coordinates' 909 write(std_out,*) ' or 2) for crystallographic coordinates' 910 read(std_in,*) inpopt 911 if (inpopt==1 .or. inpopt==2) okinp=1 912 end do 913 914 if (inpopt==1) then 915 916 write(std_out,*) ' Input point Cartesian Coord: X Y Z' 917 read(std_in,*) rcart(1),rcart(2),rcart(3) 918 call reduce(rr,rcart,rprimd) 919 write(std_out,'(a,3es16.6)' ) ' Crystallographic coordinates: ',rr(1:3) 920 921 else 922 923 write(std_out,*) ' Input point Crystallographic Coord: X Y Z' 924 read(std_in,*) rr(1),rr(2),rr(3) 925 926 do mu=1,3 927 rcart(mu)=rprimd(mu,1)*rr(1)+rprimd(mu,2)*rr(2)+rprimd(mu,3)*rr(3) 928 end do 929 930 write(std_out,*) ' Cartesian coordinates : ' 931 write(std_out,'(3es16.6)' ) rcart(1),rcart(2),rcart(3) 932 933 end if 934 935 !At this moment the code knows everything needed about the geometric input 936 !It will further proceed to calculate the interpolation at the demanded point 937 938 rr(1)=mod(mod(rr(1),1._dp)+1._dp,1._dp) 939 rr(2)=mod(mod(rr(2),1._dp)+1._dp,1._dp) 940 rr(3)=mod(mod(rr(3),1._dp)+1._dp,1._dp) 941 942 write(std_out,'(a,es16.6)' ) ' X coordinate, r1 is:',rr(1) 943 write(std_out,'(a,es16.6)' ) ' Y coordinate, r2 is:',rr(2) 944 write(std_out,'(a,es16.6)' ) ' Z coordinate, r3 is:',rr(3) 945 946 !devalt = total density value 947 !devalu = spin-up density value 948 !devald = spin-down density value 949 !devalm = magnetization density value 950 denvaltt = interpol3d(rr,nr1,nr2,nr3,gridt) 951 if(nspden==2 .or. nspden==4)then 952 denvalux = interpol3d(rr,nr1,nr2,nr3,gridu) 953 denvaldy = interpol3d(rr,nr1,nr2,nr3,gridd) 954 denvalmz = interpol3d(rr,nr1,nr2,nr3,gridm) 955 end if 956 write(std_out,*) 957 write(std_out,*)'---------------------------------------------' 958 write(std_out,'(a,es16.6)') ' Non-spin-polarized value= ',denvaltt 959 if(nspden==2)then 960 write(std_out,'(a,es16.6)')' Spin-up value = ',denvalux 961 write(std_out,'(a,es16.6)')' Spin-down value = ',denvaldy 962 write(std_out,'(a,es16.6)')' Spin difference value = ',denvalmz 963 else if(nspden==4)then 964 write(std_out,'(a,es16.6)')' x component = ',denvalux 965 write(std_out,'(a,es16.6)')' y component = ',denvaldy 966 write(std_out,'(a,es16.6)')' z component = ',denvalmz 967 end if 968 write(std_out,*)'---------------------------------------------' 969 970 end subroutine cut3d_pointint
m_cut3d/cut3d_rrho [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_rrho
FUNCTION
Reads in the charge in mkdens3D format The file was opened in the calling program, unit number 19. The header was already read in the case of the unformatted file
INPUTS
path=File name varname=Name of the netcdf variable to be read. iomode=flag specifying the IO library. nr1=grid_full size along x nr2=grid_full size along y nr3=grid_full size along z nspden=number of spin polartized densities (1 for non-spin polarized, 2 for spin-polarized)
OUTPUT
grid_full(nr1,nr2,nr3)=grid_full matrix
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
1064 subroutine cut3d_rrho(path,varname,iomode,grid_full,nr1,nr2,nr3,nspden) 1065 1066 1067 !This section has been created automatically by the script Abilint (TD). 1068 !Do not modify the following lines by hand. 1069 #undef ABI_FUNC 1070 #define ABI_FUNC 'cut3d_rrho' 1071 !End of the abilint section 1072 1073 implicit none 1074 1075 !Arguments------------------------------------------------------------- 1076 !scalars 1077 integer,intent(in) :: iomode,nr1,nr2,nr3,nspden 1078 character(len=*),intent(in) :: path,varname 1079 !arrays 1080 real(dp),intent(out),target :: grid_full(nr1,nr2,nr3,nspden) 1081 1082 !Local variables-------------------------------------------------------- 1083 !scalars 1084 integer :: ispden,unt,fform 1085 #ifdef HAVE_NETCDF 1086 integer :: varid 1087 #endif 1088 character(len=500) :: msg 1089 type(hdr_type) :: hdr 1090 1091 ! ************************************************************************* 1092 1093 select case (iomode) 1094 case (IO_MODE_FORTRAN) 1095 !Unformatted, on one record 1096 if (open_file(path, msg, newunit=unt, form='unformatted', status='old', action="read") /= 0) then 1097 MSG_ERROR(msg) 1098 end if 1099 call hdr_fort_read(hdr, unt, fform) 1100 ABI_CHECK(fform /= 0, sjoin("Error while reading:", path)) 1101 call hdr_free(hdr) 1102 1103 do ispden=1,nspden 1104 read(unit=unt) grid_full(1:nr1,1:nr2,1:nr3,ispden) 1105 end do 1106 1107 close(unt) 1108 1109 case (IO_MODE_ETSF) 1110 ! ETSF case 1111 #ifdef HAVE_NETCDF 1112 NCF_CHECK(nctk_open_read(unt, path, xmpi_comm_self)) 1113 NCF_CHECK(nf90_inq_varid(unt, varname, varid)) 1114 ! [cplex, n1, n2, n3, nspden] 1115 NCF_CHECK(nf90_get_var(unt, varid, grid_full, start=[1,1,1,1,1], count=[1, nr1,nr2,nr3,nspden])) 1116 NCF_CHECK(nf90_close(unt)) 1117 #else 1118 MSG_ERROR('netcdf support is not compiled. Reconfigure with --enable-netcdf.') 1119 #endif 1120 1121 case default 1122 MSG_BUG(sjoin("invalid iomode:", itoa(iomode))) 1123 end select 1124 1125 end subroutine cut3d_rrho
m_cut3d/cut3d_volumeint [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_volumeint
FUNCTION
Computes the values within a volume
INPUTS
gridtt(nr1,nr2,nr3)=Total density gridux(nr1,nr2,nr3)=spin-Up density, or magnetization density in X direction griddy(nr1,nr2,nr3)=spin-Down density, or magnetization density in Y direction gridmz(nr1,nr2,nr3)=spin-polarization density or magnetization density in Z direction natom=integer number of atoms nr1=grid size along x nr2=grid size along y nr3=grid size along z nspden=number of spin-density components rprimd(3,3)=orientation of the unit cell in 3D tau(3,natom)=list of atoms in cartesian coordinates
OUTPUT
only writing
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
1216 subroutine cut3d_volumeint(gridtt,gridux,griddy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,tau) 1217 1218 1219 !This section has been created automatically by the script Abilint (TD). 1220 !Do not modify the following lines by hand. 1221 #undef ABI_FUNC 1222 #define ABI_FUNC 'cut3d_volumeint' 1223 !End of the abilint section 1224 1225 implicit none 1226 1227 !Arguments ------------------------------------ 1228 !scalars 1229 integer,intent(in) :: natom,nr1,nr2,nr3,nspden 1230 !arrays 1231 real(dp),intent(in) :: griddy(nr1,nr2,nr3) 1232 real(dp),intent(in) :: gridmz(nr1,nr2,nr3),gridtt(nr1,nr2,nr3) 1233 real(dp),intent(in) :: gridux(nr1,nr2,nr3),rprimd(3,3),tau(3,natom) 1234 1235 !Local variables ------------------------- 1236 !scalars 1237 integer :: fileformattype,iat,idir,ii,inpopt,itypat,k1,k2,k3,mu,nresolh 1238 integer :: nresoll,nresolw,okhkl,okparam,planetype 1239 integer :: referenceposition,unt 1240 real(dp) :: denvaldy,denvalmz,denvaltt,denvalux,height 1241 real(dp) :: length,width 1242 real(dp) :: xm,xp,ym,yp,zm,zp 1243 character(len=fnlen) :: filnam 1244 character(len=500) :: msg 1245 !arrays 1246 integer :: hkl(3) 1247 real(dp) :: cent(3),centpl(3),mminv(3,3),r1(3),r2(3),r3(3),rcart(3) 1248 real(dp) :: rr(3),x1(3),x2(3),x3(3),xcart(3) 1249 real(dp),allocatable :: rhomacudy(:,:),rhomacumz(:,:),rhomacutt(:,:) 1250 real(dp),allocatable :: rhomacuux(:,:) 1251 1252 ! ********************************************************************* 1253 1254 call matr3inv(rprimd,mminv) 1255 !Start of the real input of the volume orientation 1256 1257 write(std_out,*) 1258 write(std_out,*) ' The volume is an orthogonal prism, that is defined by: ' 1259 write(std_out,*) ' the basal plane and' 1260 write(std_out,*) ' the height perpendicular to the basal plane' 1261 write(std_out,*) 1262 write(std_out,*) ' First you will define the basal plane ' 1263 write(std_out,*) ' second you will define the height' 1264 write(std_out,*) ' and third you will define the basal plane position ' 1265 write(std_out,*) ' along the height vector' 1266 1267 do 1268 write(std_out,*) 1269 write(std_out,*) ' Type 1) for a plane passing through 3 atoms' 1270 write(std_out,*) ' or 2) for a plane passing through 3 cartesian points' 1271 write(std_out,*) ' or 3) for a plane passing through 3 crystallographic points' 1272 write(std_out,*) ' or 4) for a plane parallel to a crystallographic plane' 1273 write(std_out,*) ' or 5) for a plane orthogonal to a cartesian direction' 1274 write(std_out,*) ' or 6) for a plane orthogonal to a crystallographic direction' 1275 write(std_out,*) ' or 0) to stop' 1276 read(std_in,*) itypat 1277 1278 select case (itypat) 1279 1280 case (0) 1281 stop 1282 1283 ! A plane passing through 3 atoms 1284 case (1) 1285 write(std_out,*) ' The X axis will be through atoms: 1,2 ' 1286 write(std_out,*) ' Define each atom by its species and its number:' 1287 write(std_out,*) ' -> atom 1 (iat):' 1288 read(std_in,*) iat 1289 x1(1)=tau(1,iat) 1290 x1(2)=tau(2,iat) 1291 x1(3)=tau(3,iat) 1292 write(std_out,'(a,3f10.6)') ' position: ',x1 1293 write(std_out,*) 1294 write(std_out,*) ' -> atom 2 (iat):' 1295 read(std_in,*) iat 1296 x2(1)=tau(1,iat) 1297 x2(2)=tau(2,iat) 1298 x2(3)=tau(3,iat) 1299 write(std_out,'(a,3f10.6)') ' position: ',x2 1300 write(std_out,*) 1301 write(std_out,*) ' -> atom 3 (iat):' 1302 read(std_in,*) iat 1303 x3(1)=tau(1,iat) 1304 x3(2)=tau(2,iat) 1305 x3(3)=tau(3,iat) 1306 write(std_out,'(a,3f10.6)') ' position: ',x3 1307 write(std_out,*) 1308 1309 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 1310 do idir=1,3 1311 x2(idir)=x2(idir)-x1(idir) 1312 x3(idir)=x3(idir)-x1(idir) 1313 end do 1314 call normalize(x2) 1315 call vdot(x3,x2,x1) 1316 call normalize(x1) 1317 call vdot(x2,x1,x3) 1318 call normalize(x3) 1319 exit 1320 1321 ! A plane passing through 3 cartesian points 1322 case (2) 1323 write(std_out,*) ' The X axis will be through points: 1,2 ' 1324 write(std_out,*) ' Define each :point coordinates' 1325 write(std_out,*) ' -> point 1: X-coord Y-coord Z-coord:' 1326 read(std_in,*) xcart 1327 x1(:)=xcart(:) 1328 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x1 1329 write(std_out,*) 1330 write(std_out,*) ' -> point 2: X-coord Y-coord Z-coord:' 1331 read(std_in,*) xcart 1332 x2(:)=xcart(:) 1333 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x2 1334 write(std_out,*) 1335 write(std_out,*) ' -> point 3: X-coord Y-coord Z-coord:' 1336 read(std_in,*) xcart 1337 x3(:)=xcart(:) 1338 write(std_out,'(a,3f10.6)') ' crystallographic position: ',x3 1339 write(std_out,*) 1340 1341 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 1342 do idir=1,3 1343 x2(idir)=x2(idir)-x1(idir) 1344 x3(idir)=x3(idir)-x1(idir) 1345 end do 1346 call normalize(x2) 1347 call vdot(x3,x2,x1) 1348 call normalize(x1) 1349 call vdot(x2,x1,x3) 1350 call normalize(x3) 1351 exit 1352 1353 ! A plane passing through 3 crystallographic points 1354 case (3) 1355 write(std_out,*) ' The X axis will be through points: 1,2 ' 1356 write(std_out,*) ' Define each :point coordinates' 1357 write(std_out,*) ' -> point 1: X-coord Y-coord Z-coord:' 1358 read(std_in,*) r1 1359 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r1 1360 write(std_out,*) 1361 write(std_out,*) ' -> point 2: X-coord Y-coord Z-coord:' 1362 read(std_in,*) r2 1363 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r2 1364 write(std_out,*) 1365 write(std_out,*) ' -> point 3: X-coord Y-coord Z-coord:' 1366 read(std_in,*) r3 1367 write(std_out,'(a,3f10.6)') ' crystallographic position: ',r3 1368 write(std_out,*) 1369 1370 ! Transforms the points coordinates into cartesian 1371 do mu=1,3 1372 x1(mu)=rprimd(mu,1)*r1(1)+rprimd(mu,2)*r1(2)+rprimd(mu,3)*r1(3) 1373 x2(mu)=rprimd(mu,1)*r2(1)+rprimd(mu,2)*r2(2)+rprimd(mu,3)*r2(3) 1374 x3(mu)=rprimd(mu,1)*r3(1)+rprimd(mu,2)*r3(2)+rprimd(mu,3)*r3(3) 1375 end do 1376 1377 write(std_out,*) ' Cartesian positions:' 1378 write(std_out,*) x1 1379 write(std_out,*) x2 1380 write(std_out,*) x3 1381 1382 ! Compute the 3 orthogonal normalized vectors from x2-x1, x3-x1 1383 do idir=1,3 1384 x2(idir)=x2(idir)-x1(idir) 1385 x3(idir)=x3(idir)-x1(idir) 1386 end do 1387 call normalize(x2) 1388 call vdot(x3,x2,x1) 1389 call normalize(x1) 1390 call vdot(x2,x1,x3) 1391 call normalize(x3) 1392 exit 1393 1394 ! A plane parallel to a crystallographic plane 1395 case (4) 1396 okhkl=0 1397 do while (okhkl==0) 1398 write(std_out,*) ' Enter plane coordinates:' 1399 write(std_out,*) ' ->H K L ' 1400 read(std_in,*) hkl 1401 if (.not. (hkl(1)==0 .and. hkl(2)==0 .and. hkl(3)==0)) okhkl=1 1402 end do 1403 write(std_out,*) ' Miller indices are:',hkl 1404 1405 do ii=1,3 1406 x1(ii)=mminv(ii,1)*hkl(1) + mminv(ii,2)*hkl(2) + mminv(ii,3)*hkl(3) 1407 end do 1408 write(std_out,*) ' Orthogonal vector to the plane',x1 1409 1410 call normalize(x1) 1411 if((x1(1).ne.0).or.(x1(2).ne.0)) then 1412 x2(1)=-x1(2) 1413 x2(2)=x1(1) 1414 x2(3)=0 1415 call normalize(x2) 1416 else 1417 x2(1)=1 1418 x2(2)=0 1419 x2(3)=0 1420 end if 1421 call vdot(x2,x1,x3) 1422 call normalize(x3) 1423 exit 1424 1425 ! A plane orthogonal to a cartesian direction 1426 case (5) 1427 write(std_out,*) ' Enter cartesian vector orthogonal to plane:' 1428 write(std_out,*) ' -> X-dir Y-dir Z-dir (Angstroms or Bohrs):' 1429 read(std_in,*) x1 1430 call normalize(x1) 1431 if((x1(1).ne.0).or.(x1(2).ne.0)) then 1432 x2(1)=-x1(2) 1433 x2(2)=x1(1) 1434 x2(3)=0 1435 call normalize(x2) 1436 else 1437 x2(1)=1 1438 x2(2)=0 1439 x2(3)=0 1440 end if 1441 call vdot(x1,x2,x3) 1442 call normalize(x3) 1443 exit 1444 1445 ! A plane orthogonal to a crystallographic direction 1446 case (6) 1447 write(std_out,*) ' Enter crystallographic vector orthogonal to plane:' 1448 write(std_out,*) ' -> X-dir Y-dir Z-dir (Fractional coordinates):' 1449 read(std_in,*) r1 1450 do mu=1,3 1451 x1(mu)=rprimd(mu,1)*r1(1)+rprimd(mu,2)*r1(2)+rprimd(mu,3)*r1(3) 1452 end do 1453 call normalize(x1) 1454 if(abs(x1(1))<tol10 .or. abs(x1(2)) < tol10) then 1455 x2(1)=-x1(2) 1456 x2(2)= x1(1) 1457 x2(3)= 0 1458 call normalize(x2) 1459 else 1460 x2(1)=1 1461 x2(2)=0 1462 x2(3)=0 1463 end if 1464 call vdot(x1,x2,x3) 1465 call normalize(x3) 1466 exit 1467 1468 case default 1469 write(std_out,*) ' Input option does not correspond to one available option' 1470 write(std_out,*) ' Please try again' 1471 cycle 1472 1473 end select 1474 end do 1475 1476 !At this moment the family of planes was defined 1477 !The code knows also some of the geometric input 1478 !It will proceed to the anchorage of the plane onto a point and then 1479 !to the effective calculation 1480 1481 write(std_out,*) ' Vectors: (orthogonal & normalized) ' 1482 write(std_out,'(11x,a,3f10.6)') ' X-dir in the plot ',x2 1483 write(std_out,'(11x,a,3f10.6)') ' Y-dir in the plot ',x3 1484 write(std_out,'(11x,a,3f10.6)') ' Z-dir (orth. to the plot) ',x1 1485 1486 do 1487 write(std_out,*) 1488 write(std_out,*) ' Enter reference point of plane (Bohr):' 1489 write(std_out,*) ' Type 1) for Cartesian coordinates.' 1490 write(std_out,*) ' or 2) for Crystallographic coordinates.' 1491 read(std_in,*) inpopt 1492 1493 select case (inpopt) 1494 1495 case(1) 1496 write(std_out,*) ' -> X-Coord Y-Coord Z-Coord:' 1497 read(std_in,*) cent 1498 exit 1499 case(2) 1500 write(std_out,*) ' -> X-Coord Y-Coord Z-Coord:' 1501 read(std_in,*) cent 1502 do mu=1,3 1503 rcart(mu)=rprimd(mu,1)*cent(1)+rprimd(mu,2)*cent(2)+rprimd(mu,3)*cent(3) 1504 end do 1505 cent(:)=rcart(:) 1506 write(std_out,'(a,3es16.6)' ) ' Expressed in cartesian coordinates: ',cent(1:3) 1507 exit 1508 case (3) 1509 cycle 1510 1511 end select 1512 end do 1513 1514 !End of basal plane orientation 1515 1516 !Input box dimensions now 1517 1518 write(std_out,*) 1519 write(std_out,*) ' It is now time to input the 3D box dimensions.' 1520 write(std_out,*) ' and the position of the basal plane in the box.' 1521 1522 do 1523 write(std_out,*) 1524 write(std_out,*) ' Enter in-plane width:' 1525 read(std_in,*) width 1526 write(std_out,*) ' Enter in-plane length:' 1527 read(std_in,*) length 1528 write(std_out,*) ' Enter box height:' 1529 read(std_in,*) height 1530 write(std_out,*) 1531 write(std_out,*) ' Enter the position of the basal plane in the box:' 1532 do 1533 write(std_out,*) 1534 write(std_out,*) ' Type 1) for DOWN' 1535 write(std_out,*) ' Type 2) for MIDDLE' 1536 write(std_out,*) ' Type 3) for UP' 1537 read(std_in,*) planetype 1538 1539 select case(planetype) 1540 1541 case (1) 1542 exit 1543 case (2) 1544 exit 1545 case (3) 1546 exit 1547 case default 1548 cycle 1549 1550 end select 1551 end do 1552 1553 write(std_out,*) ' Enter the position of the reference point in the basal plane ' 1554 1555 do 1556 write(std_out,*) 1557 write(std_out,*) ' Type 1) for CENTRAL position ' 1558 write(std_out,*) ' Type 2) for CORNER(0,0) position ' 1559 read(std_in,*) referenceposition 1560 1561 select case(referenceposition) 1562 1563 case (1) 1564 exit 1565 case (2) 1566 exit 1567 case default 1568 cycle 1569 1570 end select 1571 end do 1572 1573 write(std_out,*) 1574 write(std_out,*) ' Enter the box grid values:' 1575 write(std_out,*) ' Enter plane resolution in width:' 1576 read(std_in,*) nresolw 1577 write(std_out,*) ' Enter plane resolution in lenth:' 1578 read(std_in,*) nresoll 1579 write(std_out,*) ' Enter height resolution:' 1580 read(std_in,*) nresolh 1581 write(std_out,*) 1582 write(std_out,*) ch10,' Enter the name of an output file:' 1583 if (read_string(filnam, unit=std_in) /= 0) then 1584 MSG_ERROR("Fatal error!") 1585 end if 1586 write(std_out,*) ' The name of your file is : ',trim(filnam) 1587 1588 do 1589 write(std_out,*) ' Enter the format of the output file:' 1590 write(std_out,*) ' Type 1=> ASCII formatted' 1591 write(std_out,*) ' Type 2=> 3D index + data, ASCII formatted' 1592 write(std_out,*) ' Type 3=> Molekel formatted' 1593 read(std_in,*) fileformattype 1594 if (fileformattype>=1 .and. fileformattype<=3) then 1595 exit 1596 else 1597 cycle 1598 end if 1599 end do 1600 1601 write(std_out,*) ' You asked for a 3d box of:' 1602 write(std_out,*) length,' x ',width,' x ',height 1603 write(std_out,*) ' With a resolution of ;' 1604 write(std_out,*) nresoll,' x ',nresolw,' x ',nresolh 1605 write(std_out,*) ' The result will be redirected to the file: ',trim(filnam) 1606 if (fileformattype==1) then 1607 write(std_out,*) ' ASCII formatted' 1608 else if (fileformattype==2) then 1609 write(std_out,*) ' 3d index + data, ASCII formatted' 1610 else if (fileformattype==3) then 1611 write(std_out,*) ' Molekel formatted' 1612 end if 1613 write(std_out,*) ' These parameters may still be changed.' 1614 write(std_out,*) ' Are you sure you want to keep them? (1=default=yes,2=no) ' 1615 read(std_in,*) okparam 1616 if (okparam==2) then 1617 cycle 1618 else 1619 exit 1620 end if 1621 end do 1622 1623 !Write the header of the Molekel input file 1624 1625 if (fileformattype==1 .or. fileformattype==2) then 1626 if (open_file(filnam,msg,newunit=unt,status='unknown') /= 0) then 1627 MSG_ERROR(msg) 1628 end if 1629 else if (fileformattype==3) then 1630 if (open_file(filnam,msg,newunit=unt,form='unformatted') /= 0) then 1631 MSG_ERROR(msg) 1632 end if 1633 1634 xm=0 1635 xp=length 1636 ym=0 1637 yp=width 1638 zm=0 1639 zp=height 1640 1641 write(std_out,'(a)' )& 1642 & ' Extremas of the cube in which the system is placed (x,y,z), in Angs.:' 1643 write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp 1644 write(std_out,'(a,a,3i5)' ) ch10,& 1645 & ' Number of points per side: ',nresolw+1,nresoll+1,nresolh+1 1646 write(std_out,'(a,a,i10,a,a)' ) ch10,& 1647 & ' Total number of points: ',(nresolw+1)*(nresoll+1)*(nresolh+1),& 1648 & ch10,ch10 1649 1650 write(unt) xm,xp,ym,yp,zm,zp,nresolw+1,nresoll+1,nresolh+1 1651 1652 end if 1653 1654 !Allocate rhomacu in case of molekel output format 1655 ABI_ALLOCATE(rhomacutt,(nresoll+1,nresolw+1)) 1656 ABI_ALLOCATE(rhomacuux,(nresoll+1,nresolw+1)) 1657 ABI_ALLOCATE(rhomacudy,(nresoll+1,nresolw+1)) 1658 ABI_ALLOCATE(rhomacumz,(nresoll+1,nresolw+1)) 1659 1660 do k1=0,nresolh 1661 1662 select case (planetype) 1663 1664 ! Basal plane at the bottom 1665 case (1) 1666 centpl(1)=cent(1)+k1*x1(1)*height/nresolh 1667 centpl(2)=cent(2)+k1*x1(2)*height/nresolh 1668 centpl(3)=cent(3)+k1*x1(3)*height/nresolh 1669 1670 ! Basal plane in the middle 1671 case (2) 1672 centpl(1)=cent(1)+(k1-nresolh/2)*x1(1)*height/nresolh 1673 centpl(2)=cent(2)+(k1-nresolh/2)*x1(2)*height/nresolh 1674 centpl(3)=cent(3)+(k1-nresolh/2)*x1(3)*height/nresolh 1675 1676 ! Basal plane on the top 1677 case (3) 1678 centpl(1)=cent(1)+(k1-nresolh)*x1(1)*height/nresolh 1679 centpl(2)=cent(2)+(k1-nresolh)*x1(2)*height/nresolh 1680 centpl(3)=cent(3)+(k1-nresolh)*x1(3)*height/nresolh 1681 1682 end select 1683 1684 do k3=0,nresolw 1685 do k2=0,nresoll 1686 1687 select case(referenceposition) 1688 1689 ! Reference point in the middle of the basal plane 1690 case (1) 1691 rcart(1)=centpl(1) + (k2-nresoll/2)*x2(1)*length/nresoll + (k3-nresolw/2)*x3(1)*width/nresolw 1692 rcart(2)=centpl(2) + (k2-nresoll/2)*x2(2)*length/nresoll + (k3-nresolw/2)*x3(2)*width/nresolw 1693 rcart(3)=centpl(3) + (k2-nresoll/2)*x2(3)*length/nresoll + (k3-nresolw/2)*x3(3)*width/nresolw 1694 1695 ! Reference point in the corner of the basal plane 1696 case (2) 1697 rcart(1)=centpl(1) + k2*x2(1)*length/nresoll + k3*x3(1)*width/nresolw 1698 rcart(2)=centpl(2) + k2*x2(2)*length/nresoll + k3*x3(2)*width/nresolw 1699 rcart(3)=centpl(3) + k2*x2(3)*length/nresoll + k3*x3(3)*width/nresolw 1700 1701 end select 1702 1703 call reduce(rr,rcart,rprimd) 1704 rr(1)=mod(mod(rr(1),1._dp)+1._dp,1._dp) 1705 rr(2)=mod(mod(rr(2),1._dp)+1._dp,1._dp) 1706 rr(3)=mod(mod(rr(3),1._dp)+1._dp,1._dp) 1707 1708 denvaltt = interpol3d(rr,nr1,nr2,nr3,gridtt) 1709 if(nspden==2 .or. nspden==4)then 1710 denvalux = interpol3d(rr,nr1,nr2,nr3,gridux) 1711 denvaldy = interpol3d(rr,nr1,nr2,nr3,griddy) 1712 denvalmz = interpol3d(rr,nr1,nr2,nr3,gridmz) 1713 end if 1714 1715 if (fileformattype==1) then 1716 if(nspden==1)then 1717 write(unt, '(es22.12)' ) denvaltt 1718 else if(nspden==2 .or. nspden==4)then 1719 write(unt, '(4(es22.12))' ) denvaltt,denvalux,denvaldy,denvalmz 1720 end if 1721 1722 else if (fileformattype==2) then 1723 if(nspden==1)then 1724 write(unt, '(4es22.12)' ) rcart, denvaltt 1725 else if(nspden==2 .or. nspden==4)then 1726 write(unt, '(3(e22.12), 4(es22.12))' ) rcart, denvaltt,denvalux,denvaldy,denvalmz 1727 end if 1728 1729 else if (fileformattype==3) then 1730 rhomacutt(k2+1,k3+1)=denvaltt 1731 if(nspden==2 .or. nspden==4)then 1732 rhomacuux(k2+1,k3+1)=denvalux 1733 rhomacudy(k2+1,k3+1)=denvaldy 1734 rhomacumz(k2+1,k3+1)=denvalmz 1735 end if 1736 end if 1737 1738 end do ! resoll 1739 write(unt, * ) 1740 end do ! resolw 1741 1742 if (fileformattype==3) then 1743 write(unt) rhomacutt(:,:) 1744 if(nspden==2 .or. nspden==4)then 1745 write(unt) rhomacuux(:,:) 1746 write(unt) rhomacudy(:,:) 1747 write(unt) rhomacumz(:,:) 1748 end if 1749 end if 1750 1751 end do 1752 1753 close(unt) 1754 1755 ABI_DEALLOCATE(rhomacutt) 1756 ABI_DEALLOCATE(rhomacuux) 1757 ABI_DEALLOCATE(rhomacudy) 1758 ABI_DEALLOCATE(rhomacumz) 1759 1760 end subroutine cut3d_volumeint
m_cut3d/cut3d_wffile [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
cut3d_wffile
FUNCTION
Part of cut3d that gives the wavefunction for one kpt,one band and one spin polarisation in real space. The output depends on the chosen option.
INPUTS
wfk_fname=Name of the WFK file. Needs an unformatted wave function from abinit. ecut= effective ecut (ecut*dilatmx**2) exchn2n3d= if 1, n2 and n3 are exchanged istwfk= input variable indicating the storage option of each k-point natom = number of atoms in the unit cell nband= size of e_kpt nkpt= number of k-points npwarr= array holding npw for each k point nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension) nspinor= number of spinorial components of the wavefunctions nsppol= number of spin polarization ntypat = number of atom type paral_kgb= parallization option, it is set to 0 in the parent subroutine rprim = orientation of the unit cell axes xcart = cartesian coordinates typat= input variable typat(natom) znucl= znucltypat(ntypat) from alchemy
OUTPUT
Depends on the option chosen. It is the wave function for the k point, band and spin polarisation chosen. It can be written in different ways. The option are describe with the option list. It is possible to output a Data Explorer file.
PARENTS
cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
1810 subroutine cut3d_wffile(wfk_fname,ecut,exchn2n3d,istwfk,kpt,natom,nband,nkpt,npwarr,& 1811 & nr1,nr2,nr3,nspinor,nsppol,ntypat,paral_kgb,rprimd,xcart,typat,znucl) 1812 1813 1814 !This section has been created automatically by the script Abilint (TD). 1815 !Do not modify the following lines by hand. 1816 #undef ABI_FUNC 1817 #define ABI_FUNC 'cut3d_wffile' 1818 !End of the abilint section 1819 1820 implicit none 1821 1822 !Arguments ----------------------------------- 1823 !scalars 1824 integer,intent(in) :: exchn2n3d,natom,nkpt,nr1,nr2,nr3,nspinor,nsppol 1825 integer,intent(in) :: ntypat,paral_kgb 1826 real(dp),intent(in) :: ecut 1827 character(len=*),intent(in) :: wfk_fname 1828 !arrays 1829 integer,intent(in) :: istwfk(nkpt),nband(nkpt),npwarr(nkpt),typat(natom) 1830 real(dp),intent(in) :: kpt(3,nkpt),rprimd(3,3),znucl(ntypat) 1831 real(dp),intent(in) :: xcart(3,natom) 1832 1833 !Local variables------------------------------- 1834 !scalars 1835 integer,parameter :: tim_fourwf0=0,tim_rwwf0=0,ndat1=1,formeig0=0 1836 integer :: cband,cgshift,ckpt,cplex,cspinor,csppol,gridshift1 1837 integer :: gridshift2,gridshift3,ia,iatom,iband,ichoice,ifile,iomode 1838 integer :: ii1,ii2,ii3,ikpt,ilang,ioffkg,iout,iprompt,ipw 1839 integer :: ir1,ir2,ir3,ivect,ixint,mband,mbess,mcg,mgfft 1840 integer :: mkmem,mlang,mpw,n4,n5,n6,nfit,npw_k 1841 integer :: nradintmax,oldcband,oldckpt,oldcspinor,oldcsppol 1842 integer :: prtsphere,select_exit,unout,iunt,rc_ylm 1843 integer :: ikpt_qps,nkpt_qps,nband_qps,iscf_qps 1844 real(dp) :: arg,bessargmax,bessint_delta,kpgmax,ratsph,tmpi,tmpr,ucvol,weight,eig_k_qps 1845 character(len=*), parameter :: INPUTfile='cut.in' 1846 character(len=1) :: outputchar 1847 character(len=10) :: string 1848 character(len=4) :: mode_paral 1849 character(len=500) :: msg 1850 character(len=fnlen) :: output,output1 1851 type(MPI_type) :: mpi_enreg 1852 type(wfk_t) :: Wfk 1853 type(jlspline_t) :: jlspl 1854 !arrays 1855 integer :: atindx(natom),iatsph(natom),ngfft(18),nradint(natom),mlang_type(ntypat) 1856 integer,allocatable :: gbound(:,:),iindex(:),kg(:,:),kg_dum(:,:),kg_k(:,:) 1857 integer,allocatable :: npwarr1(:),npwarrk1(:),npwtot1(:) 1858 real(dp) :: cmax(natom),gmet(3,3),gprimd(3,3) 1859 real(dp) :: phkxred(2,natom),ratsph_arr(natom),rmet(3,3),shift_tau(3) 1860 real(dp) :: tau2(3,natom),xred(3,natom),kpt_qps(3) 1861 real(dp) :: znucl_atom(natom) 1862 integer :: znucl_atom_int(natom) 1863 real(dp),allocatable :: bess_fit(:,:,:) 1864 real(dp),allocatable :: cg_k(:,:),cgcband(:,:),denpot(:,:,:),eig_k(:) 1865 real(dp),allocatable :: fofgout(:,:),fofr(:,:,:,:),k1(:,:) 1866 real(dp),allocatable :: kpgnorm(:),occ_k(:),ph1d(:,:),ph3d(:,:,:),rint(:) 1867 real(dp),allocatable :: sum_1atom_1ll(:,:,:),sum_1atom_1lm(:,:,:) 1868 real(dp),allocatable :: xfit(:),yfit(:),ylm_k(:,:) 1869 real(dp),allocatable :: ylmgr_dum(:,:,:) 1870 character(len=fnlen) :: fileqps 1871 character(len=fnlen),allocatable :: filename(:) 1872 complex(dp),allocatable :: ccoeff(:,:),wfg(:,:),wfg_qps(:) 1873 real(dp) :: spinvec(3) 1874 1875 ! *********************************************************************** 1876 1877 call initmpi_seq(mpi_enreg) 1878 mband=maxval(nband) 1879 ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt,mband,nsppol)) 1880 mpi_enreg%proc_distrb=0 1881 mpi_enreg%me_g0 = 1 1882 oldckpt=0 1883 oldcband=0 1884 oldcsppol=0 1885 oldcspinor=0 1886 1887 iout=-1 1888 call metric(gmet,gprimd,iout,rmet,rprimd,ucvol) 1889 1890 !get xred 1891 call xcart2xred(natom,rprimd,xcart,xred) 1892 1893 !znucl indexed by atoms 1894 znucl_atom = znucl(typat(1:natom)) 1895 znucl_atom_int = INT(znucl(typat(1:natom))) 1896 1897 do iatom=1,natom 1898 iatsph(iatom) = iatom 1899 atindx(iatom) = iatom 1900 end do 1901 1902 !max ang mom + 1 1903 mlang = 5 1904 1905 ABI_ALLOCATE(kg_dum,(3,0)) 1906 1907 ABI_ALLOCATE(ph1d,(2,(2*nr1+1+2*nr2+1+2*nr3+1)*natom)) 1908 call getph(atindx,natom,nr1,nr2,nr3,ph1d,xred) 1909 1910 do 1911 ! Get k-point, band and spin polarisation for the output 1912 if(nkpt/=1)then 1913 write(std_out,*) 1914 write(std_out,'(a,i4,a)') ' For which k-points? (1 to ',nkpt,')' 1915 read(std_in,*)ckpt 1916 ! Check if kpt exist 1917 if(ckpt<1 .or. ckpt>nkpt) then 1918 write(msg,'(a,i0)') 'Invalid k-point ',ckpt 1919 MSG_ERROR(msg) 1920 end if 1921 else 1922 ckpt=nkpt 1923 end if 1924 write(std_out,*) ' => Your k-point is : ',ckpt 1925 write(std_out,*) 1926 1927 if(nband(ckpt)/=1)then 1928 write(std_out,*) 1929 write(std_out,'(a,i5,a)') ' For which band ? (1 to ',nband(ckpt),')' 1930 read(std_in,*)cband 1931 ! Check if band number exist 1932 1933 if(cband<1 .or. cband>nband(ckpt)) then 1934 write(msg,'(a,i0)')'Invalid band number',cband 1935 MSG_ERROR(msg) 1936 end if 1937 else 1938 cband=nband(ckpt) 1939 end if 1940 write(std_out,*) ' => Your band number is : ',(cband) 1941 write(std_out,*) 1942 1943 if(nsppol/=1)then 1944 write(std_out,*) 1945 write(std_out,*) ' For which spin polarisation ?' 1946 read(std_in,*)csppol 1947 ! Check if spin polarisation exist 1948 if(csppol<1 .or. csppol>nsppol) then 1949 write(msg,'(a,i0)')'Invalid spin polarisation ',csppol 1950 MSG_ERROR(msg) 1951 end if 1952 else 1953 csppol=1 1954 end if 1955 1956 write(std_out,*) ' => Your spin polarisation number is : ',(csppol) 1957 write(std_out,*) 1958 1959 if(nspinor/=1) then 1960 write(std_out,*) ' nspinor = ', nspinor 1961 write(std_out,*) 1962 write(std_out,*) ' For which spinor component ?' 1963 read(std_in,*) cspinor 1964 ! Check if spin polarisation exist 1965 if(cspinor<1 .or. cspinor>nspinor) then 1966 write(msg,'(a,i0)')'Invalid spinor index ',cspinor 1967 MSG_ERROR(msg) 1968 end if 1969 write(std_out,*) ' => Your spinor component is : ',(cspinor) 1970 write(std_out,*) 1971 else 1972 cspinor=1 1973 end if 1974 1975 ! Reading of the data if the value of ckpt and csppol are different from oldckpt and oldcsppol 1976 ! formeig=0 gstate calculation 1977 ! formeig=1 for response function calculation 1978 if(csppol/=oldcsppol .or. ckpt/=oldckpt)then 1979 mband=maxval(nband) 1980 mpw=maxval(npwarr) 1981 mcg=mpw*nspinor*mband 1982 if (allocated (cg_k)) then 1983 ABI_DEALLOCATE(cg_k) 1984 ABI_DEALLOCATE(eig_k) 1985 ABI_DEALLOCATE(occ_k) 1986 end if 1987 ABI_ALLOCATE(cg_k,(2,mcg)) 1988 ABI_ALLOCATE(eig_k,((2*mband)**formeig0*mband)) 1989 ABI_ALLOCATE(occ_k,(mband)) 1990 1991 ! FIXME 1992 ! nband depends on (kpt,spin) 1993 iomode = iomode_from_fname(wfk_fname) 1994 call wfk_open_read(Wfk,wfk_fname,formeig0,iomode,get_unit(),xmpi_comm_self) 1995 call wfk_read_band_block(Wfk,[1,nband(ckpt)],ckpt,csppol,xmpio_single,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 1996 call wfk_close(Wfk) 1997 end if 1998 1999 if (csppol/=oldcsppol .or. ckpt/=oldckpt .or. cband/=oldcband .or. cspinor/=oldcspinor ) then 2000 ! The data of ckpt,cnsspol are in cg_k. Now we have to do the Fourier Transform of the datas 2001 2002 ngfft(1)=nr1 2003 ngfft(2)=nr2 2004 ngfft(3)=nr3 2005 ! ngfft(4) and ngfft(5) can not be even (see getng.f) 2006 if (mod(nr1,2)==0)then 2007 ngfft(4)=nr1+1 2008 else 2009 ngfft(4)=nr1 2010 end if 2011 if (mod(nr2,2)==0)then 2012 ngfft(5)=nr2+1 2013 else 2014 ngfft(5)=nr2 2015 end if 2016 ngfft(6)=nr3 2017 ! XG 020829 : 112 does not work yet for all istwfk values 2018 ngfft(7)=111 2019 ngfft(8)=256 2020 ngfft(9)=0 2021 ngfft(10)=1 2022 ngfft(11)=0 2023 ngfft(12)=ngfft(2) 2024 ngfft(13)=ngfft(3) 2025 ngfft(14)=0 2026 2027 ! if iout<0, the output of metric will not be print 2028 mode_paral='PERS' 2029 mkmem=nkpt 2030 mgfft=maxval(ngfft(1:3)) 2031 ABI_ALLOCATE(npwarr1,(nkpt)) 2032 ABI_ALLOCATE(kg,(3,mpw*mkmem)) 2033 ABI_ALLOCATE(npwtot1,(nkpt)) 2034 call init_distribfft_seq(mpi_enreg%distribfft,'c',ngfft(2),ngfft(3),'all') 2035 2036 ! Create positions index for pw 2037 call kpgio(ecut,exchn2n3d,gmet,istwfk,kg,kpt,mkmem,nband,nkpt,& 2038 & mode_paral,mpi_enreg,mpw,npwarr1,npwtot1,nsppol) 2039 2040 ioffkg=0 2041 do ikpt=1,ckpt-1 2042 ioffkg=ioffkg+npwarr1(ikpt) 2043 end do 2044 npw_k=npwarr(ckpt) 2045 ABI_ALLOCATE(gbound,(2*mgfft+8,2)) 2046 ABI_ALLOCATE(kg_k,(3,npw_k)) 2047 kg_k(:,1:npw_k)=kg(:,1+ioffkg:npw_k+ioffkg) 2048 2049 ABI_ALLOCATE(ylm_k,(npw_k,mlang*mlang)) 2050 ABI_ALLOCATE(ylmgr_dum,(npw_k,3,mlang*mlang)) 2051 2052 ! call for only the kpoint we are interested in ! 2053 ABI_ALLOCATE(k1,(3,1)) 2054 k1(:,1)=kpt(:,ckpt) 2055 ABI_ALLOCATE(npwarrk1,(1)) 2056 npwarrk1 = (/npw_k/) 2057 call initylmg(gprimd,kg_k,k1,1,mpi_enreg,mlang,npw_k,nband,1,& 2058 & npwarrk1,nsppol,0,rprimd,ylm_k,ylmgr_dum) 2059 ABI_DEALLOCATE(ylmgr_dum) 2060 ABI_DEALLOCATE(k1) 2061 ABI_DEALLOCATE(npwarrk1) 2062 2063 ! Compute the norms of the k+G vectors 2064 ABI_ALLOCATE(kpgnorm,(npw_k)) 2065 call getkpgnorm(gprimd,kpt(:,ckpt),kg_k,kpgnorm,npw_k) 2066 2067 call sphereboundary(gbound,istwfk(ckpt),kg_k,mgfft,npw_k) 2068 ! Do the Fourier Transform 2069 n4=ngfft(4) 2070 n5=ngfft(5) 2071 n6=ngfft(6) 2072 ! cplex=0 2073 cplex=1 2074 ! Complex can be set to 0 with this option(0) of fourwf 2075 2076 ! Read the QPS file if GW wavefunctions are to be analysed 2077 write(std_out,*) 'Do you want to analyze a GW wavefunction? (1=yes,0=no)' 2078 read(std_in,*) ii1 2079 write(std_out,*) '=> Your choice is :',ii1 2080 write(std_out,*) 2081 2082 if(ii1==1) then 2083 write(std_out,*) 'What is the name of the QPS file?' 2084 if (read_string(fileqps, unit=std_in) /= 0) then 2085 MSG_ERROR("Fatal error!") 2086 end if 2087 ! Checking the existence of data file 2088 if (.not. file_exists(fileqps)) then 2089 MSG_ERROR(sjoin('Missing data file:', fileqps)) 2090 end if 2091 2092 if (open_file(fileqps, msg, newunit=iunt, status='old',form='formatted') /= 0) then 2093 MSG_ERROR(msg) 2094 end if 2095 2096 read(iunt,*) iscf_qps 2097 read(iunt,*) nkpt_qps 2098 read(iunt,*) nband_qps 2099 read(iunt,*) ikpt_qps 2100 2101 ABI_ALLOCATE(ccoeff,(nband_qps,nband_qps)) 2102 do ikpt=1,ckpt ! nkpt_qps 2103 read(iunt,*) kpt_qps(:) 2104 do iband=1,nband_qps 2105 read(iunt,*) eig_k_qps 2106 read(iunt,*) ccoeff(:,iband) 2107 end do 2108 end do 2109 close(iunt) 2110 2111 ABI_ALLOCATE(wfg,(npw_k,nband_qps)) 2112 ABI_ALLOCATE(wfg_qps,(npw_k)) 2113 do iband=1,nband_qps 2114 cgshift=(iband-1)*npw_k*nspinor + (cspinor-1)*npw_k 2115 wfg(:,iband) = dcmplx( cg_k(1,cgshift+1:cgshift+npw_k),cg_k(2,cgshift+1:cgshift+npw_k) ) 2116 end do 2117 2118 wfg_qps = matmul( wfg(:,:) , ccoeff(:,cband) ) 2119 2120 ! write(std_out,*) 'norm',SUM( abs(wfg(:,cband))**2 ) 2121 ! write(std_out,*) 'norm',SUM( abs(wfg_qps(:))**2 ) 2122 ABI_DEALLOCATE(ccoeff) 2123 ABI_DEALLOCATE(wfg) 2124 ABI_ALLOCATE(cgcband,(2,npw_k*nspinor)) 2125 cgcband = zero 2126 cgcband(1,:)= real(wfg_qps(:)) 2127 cgcband(2,:)= aimag(wfg_qps(:)) 2128 ABI_DEALLOCATE(wfg_qps) 2129 2130 else ! not a GW wavefunction 2131 2132 ! get spin vector for present state 2133 cgshift=(cband-1)*npw_k*nspinor 2134 ABI_ALLOCATE(cgcband,(2,npw_k*nspinor)) 2135 cgcband(:,1:npw_k*nspinor)=cg_k(:,cgshift+1:cgshift+nspinor*npw_k) 2136 end if ! test QPS wavefunction from GW 2137 2138 if (nspinor == 2) then 2139 call cg_getspin(cgcband, npw_k, spinvec) 2140 write(std_out,'(a,6E20.10)' ) ' spin vector for this state = ', (spinvec) 2141 end if 2142 2143 ! Fix the phase of cgcband, for portability reasons 2144 ! call fxphas(cgcband,cgcband,0,npw_k,1,npw_k,0) 2145 2146 ABI_ALLOCATE(denpot,(cplex*n4,n5,n6)) 2147 ABI_ALLOCATE(fofgout,(2,npw_k)) 2148 ABI_ALLOCATE(fofr,(2,n4,n5,n6)) 2149 2150 call fourwf(cplex,denpot,cgcband(:,(cspinor-1)*npw_k+1:cspinor*npw_k),fofgout,fofr,gbound,gbound,& 2151 & istwfk(ckpt),kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,& 2152 & npw_k,n4,n5,n6,0,paral_kgb,tim_fourwf0,weight,weight) 2153 2154 ! TODO 2155 ! call fft_ug_dp(npw_k,nfft,nspinor,ndat1,mgfft,ngfft,istwf_k(ckpt),kg_k,gbound,cgcband,fofr) 2156 2157 ! Analyse wavefunction inside atomic sphere 2158 2159 write(std_out,'(a)' ) ' Do you want the atomic analysis for this state : ' 2160 write(std_out,'(a,2i5,a)' ) ' (kpt,band)= (',ckpt,cband,')? ' 2161 write(std_out,'(a)' ) ' If yes, enter the radius of the atomic spheres, in bohr ' 2162 write(std_out,'(a)' ) ' If no, enter 0 ' 2163 read (std_in,*) ratsph 2164 write(std_out,'(a,f16.8,a)' ) ' You entered ratsph=',ratsph,' Bohr ' 2165 2166 if (ratsph >= tol10) then 2167 2168 write(std_out,'(3a)' ) ch10,' Atomic sphere analysis ',ch10 2169 2170 ! Init bessel function integral for recip_ylm: max ang mom + 1 2171 mlang = 5 2172 bessint_delta = 0.1_dp 2173 kpgmax = sqrt(ecut) 2174 bessargmax = ratsph*two_pi*kpgmax 2175 mbess = int (bessargmax / bessint_delta) + 1 2176 bessargmax = bessint_delta*mbess 2177 2178 ! Intervals in radial integration 2179 nradintmax = mbess 2180 nradint(1:natom)=nradintmax 2181 2182 write(std_out,'(a,2es16.6,i6)')' wffile : kpgmax, bessargmax, nradint = ', kpgmax, bessargmax,nradintmax 2183 2184 ! Initialize general Bessel function array on uniform grid xx, from 0 to (2 \pi |k+G|_{max} |r_{max}|) 2185 ABI_ALLOCATE(rint,(nradintmax)) 2186 2187 jlspl = jlspline_new(mbess, bessint_delta, mlang) 2188 2189 ABI_ALLOCATE(bess_fit,(mpw,nradintmax,mlang)) 2190 ABI_ALLOCATE(xfit,(npw_k)) 2191 ABI_ALLOCATE(yfit,(npw_k)) 2192 ABI_ALLOCATE(iindex,(npw_k)) 2193 nfit = npw_k 2194 2195 do ixint=1,nradintmax 2196 rint(ixint) = (ixint-1)*ratsph / (nradintmax-1) 2197 2198 do ipw=1,npw_k 2199 xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint) 2200 iindex(ipw) = ipw 2201 end do 2202 call sort_dp (npw_k,xfit,iindex,tol14) 2203 do ilang=1,mlang 2204 call splint(mbess,jlspl%xx,jlspl%bess_spl(:,ilang),jlspl%bess_spl_der(:,ilang),nfit,xfit,yfit) 2205 ! Re-order results for different G vectors 2206 do ipw=1,npw_k 2207 bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw) 2208 end do 2209 end do ! ipw 2210 end do ! ixint 2211 2212 ! Construct phases ph3d for all G vectors in present sphere make phkred for all atoms 2213 do ia=1,natom 2214 iatom=atindx(ia) 2215 arg=two_pi*( kpt(1,ckpt)*xred(1,ia) + kpt(2,ckpt)*xred(2,ia) + kpt(3,ckpt)*xred(3,ia)) 2216 phkxred(1,iatom)=cos(arg) 2217 phkxred(2,iatom)=sin(arg) 2218 end do 2219 2220 ABI_ALLOCATE(ph3d,(2,npw_k,natom)) 2221 ! Get full phases exp (2 pi i (k+G).x_tau) in ph3d 2222 call ph1d3d(1,natom,kg_k,natom,natom,npw_k,nr1,nr2,nr3,phkxred,ph1d,ph3d) 2223 2224 ABI_ALLOCATE(sum_1atom_1ll,(nspinor**2,mlang,natom)) 2225 ABI_ALLOCATE(sum_1atom_1lm,(nspinor**2,mlang**2,natom)) 2226 prtsphere=1 2227 ratsph_arr(:)=ratsph 2228 2229 rc_ylm = 1 ! Real or Complex spherical harmonics. 2230 mlang_type = 5 2231 2232 call recip_ylm (bess_fit,cgcband,istwfk(ckpt),mpi_enreg,& 2233 & nradint,nradintmax,mlang,mpw,natom,typat,mlang_type,npw_k,nspinor,ph3d,prtsphere,rint,& 2234 & ratsph_arr,rc_ylm,sum_1atom_1ll,sum_1atom_1lm,ucvol,ylm_k,znucl_atom) 2235 2236 call dens_in_sph(cmax,cgcband(:,(cspinor-1)*npw_k+1:cspinor*npw_k),gmet,istwfk(ckpt),& 2237 & kg_k,natom,ngfft,mpi_enreg,npw_k,paral_kgb,ph1d,ratsph_arr,ucvol) 2238 2239 write(std_out,'(a)' )' Charge in the sphere around each atom ' 2240 do iatom=1,natom 2241 write(std_out,'(a,i4,a,f14.8)' ) ' Atom number ',iatom,' : charge =',cmax(iatom) 2242 end do 2243 2244 ABI_DEALLOCATE(sum_1atom_1ll) 2245 ABI_DEALLOCATE(sum_1atom_1lm) 2246 ABI_DEALLOCATE(ph3d) 2247 ABI_DEALLOCATE(iindex) 2248 ABI_DEALLOCATE(yfit) 2249 ABI_DEALLOCATE(xfit) 2250 ABI_DEALLOCATE(bess_fit) 2251 call jlspline_free(jlspl) 2252 ABI_DEALLOCATE(rint) 2253 end if ! ratsph < 0 = end if for atomic sphere analysis 2254 2255 ABI_DEALLOCATE(cgcband) 2256 ABI_DEALLOCATE(fofgout) 2257 ABI_DEALLOCATE(denpot) 2258 ABI_DEALLOCATE(gbound) 2259 ABI_DEALLOCATE(kg_k) 2260 ABI_DEALLOCATE(npwarr1) 2261 ABI_DEALLOCATE(kg) 2262 ABI_DEALLOCATE(npwtot1) 2263 ABI_DEALLOCATE(kpgnorm) 2264 ABI_DEALLOCATE(ylm_k) 2265 call destroy_distribfft(mpi_enreg%distribfft) 2266 end if 2267 2268 write(std_out,*) 2269 write(std_out,*) ' 3D wave function was read. ','Ready for further treatment.' 2270 write(std_out,*) 2271 write(std_out,*) '============================','===============================' 2272 write(std_out,*) 2273 2274 ! ------------------------------------------------------------------------ 2275 2276 ! At this moment all the input is done 2277 ! The code knows the geometry of the system, and the data file. 2278 2279 2280 select_exit = 0 2281 do while (select_exit == 0) 2282 write(std_out,*) ' What is your choice ? Type:' 2283 write(std_out,*) ' 0 => exit to k-point / band / spin-pol loop' 2284 write(std_out,*) ' 1 => 3D formatted real and imaginary data' 2285 write(std_out,*) ' (output the bare 3D data - two column,R,I)' 2286 write(std_out,*) ' 2 => 3D formatted real data' 2287 write(std_out,*) ' (output the bare 3D data - one column)' 2288 write(std_out,*) ' 3 => 3D formatted imaginary data' 2289 write(std_out,*) ' (output the bare 3D data - one column)' 2290 write(std_out,*) ' 4 => 3D indexed real and imaginary data' 2291 write(std_out,*) ' (3D data, preceeded by 3D index)' 2292 write(std_out,*) ' 5 => 3D indexed real data' 2293 write(std_out,*) ' (bare 3D data, preceeded by 3D index)' 2294 write(std_out,*) ' 6 => 3D indexed imaginary data' 2295 write(std_out,*) ' (bare 3D data, preceeded by 3D index)' 2296 write(std_out,*) ' 7 => 3D Data Explorer formatted data ' 2297 write(std_out,*) ' (Real file and Imaginary file)' 2298 write(std_out,*) ' 8 => 3D Data Explorer formatted data ' 2299 write(std_out,*) ' (Only the Real file)' 2300 write(std_out,*) ' 9 => 3D Data Explorer formatted data ' 2301 write(std_out,*) ' (Only the Imaginary file)' 2302 write(std_out,*) ' 10 => 3D Data Explorer formatted data and position files' 2303 write(std_out,*) ' 11 => XCrysden formatted data (norm of wf) and position files' 2304 write(std_out,*) ' 12 => NetCDF data and position file' 2305 write(std_out,*) ' 13 => XCrysden/VENUS wavefunction (real part of data)' 2306 write(std_out,*) ' 14 => Gaussian/cube wavefunction module' 2307 read(std_in,*) ichoice 2308 write(std_out,'(a,a,i2,a)' ) ch10,' Your choice is ',ichoice,char(10) 2309 2310 if (ichoice>0 .and. ichoice<15)then 2311 write(std_out,*) ch10,' Enter the root of an output file:' 2312 if (read_string(output1, unit=std_in) /= 0) then 2313 MSG_ERROR("Fatal error!") 2314 end if 2315 write(std_out,*) ' The root of your file is : ',trim(output1) 2316 output=trim(output1) 2317 call int2char10(ckpt,string) 2318 output=trim(output)//'_k'//trim(string) 2319 call int2char10(cband,string) 2320 output=trim(output)//'_b'//trim(string) 2321 if (nsppol > 1) then 2322 call int2char10(csppol,string) 2323 output=trim(output)//'_sppol'//trim(string) 2324 end if 2325 if (nspinor > 1) then 2326 call int2char10(cspinor,string) 2327 output=trim(output)//'_spinor'//trim(string) 2328 end if 2329 2330 write(std_out,*) ' The corresponding filename is : ',trim(output) 2331 end if 2332 2333 select case (ichoice) 2334 2335 case (1) ! data R,I 2336 write(std_out,*) 2337 write(std_out,*) 'Give 1 file of 3D formatted real and imaginary data' 2338 write(std_out,*) 'The first column is the real data' 2339 write(std_out,*) 'The second column is the imaginary data' 2340 write(std_out,*) 2341 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2342 MSG_ERROR(msg) 2343 end if 2344 call print_fofr_ri("RI",nr1,nr2,nr3,n4,n5,n6,fofr,unit=unout) 2345 close(unout) 2346 exit 2347 2348 case (2) ! data R 2349 write(std_out,*) 2350 write(std_out,*) 'Give 1 file of 3D formatted real data' 2351 write(std_out,*) 'The only column is the real data' 2352 write(std_out,*) 2353 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2354 MSG_ERROR(msg) 2355 end if 2356 call print_fofr_ri("R",nr1,nr2,nr3,n4,n5,n6,fofr,unit=unout) 2357 close(unout) 2358 exit 2359 2360 case (3) ! data I 2361 write(std_out,*) 2362 write(std_out,*) 'Give 1 file of 3D formatted real data' 2363 write(std_out,*) 'The only column is the imaginary data' 2364 write(std_out,*) 2365 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2366 MSG_ERROR(msg) 2367 end if 2368 call print_fofr_ri("I",nr1,nr2,nr3,n4,n5,n6,fofr,unit=unout) 2369 close(unout) 2370 exit 2371 2372 case (4) ! coord(x,y,z) data R,I 2373 write(std_out,*) 2374 write(std_out,*) 'Give 1 file of 3D formatted data' 2375 write(std_out,*) 'The first three columns are the x,y,z positions(Angstrom)' 2376 write(std_out,*) 'The fourth column is the real data' 2377 write(std_out,*) 'The fifth column is the imaginary data' 2378 write(std_out,*) 2379 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2380 MSG_ERROR(msg) 2381 end if 2382 call print_fofr_xyzri("RI",nr1,nr2,nr3,n4,n5,n6,fofr,rprimd,conv_fact=Bohr_Ang,unit=unout) 2383 close(unout) 2384 exit 2385 2386 case (5) ! coord(x,y,z) data R 2387 write(std_out,*) 2388 write(std_out,*) 'Give 1 file of 3D formatted data' 2389 write(std_out,*) 'The first three columns are the x,y,z positions(Angstrom)' 2390 write(std_out,*) 'The fourth column is the real data' 2391 write(std_out,*) 2392 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2393 MSG_ERROR(msg) 2394 end if 2395 call print_fofr_xyzri("R",nr1,nr2,nr3,n4,n5,n6,fofr,rprimd,conv_fact=Bohr_Ang,unit=unout) 2396 close(unout) 2397 exit 2398 2399 case(6) ! coord(x,y,z) data I 2400 write(std_out,*) 2401 write(std_out,*) 'Give 1 file of 3D formatted data' 2402 write(std_out,*) 'The first three columns are the x,y,z positions(Angstrom)' 2403 write(std_out,*) 'The fourth column is the imaginary data' 2404 write(std_out,*) 2405 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2406 MSG_ERROR(msg) 2407 end if 2408 call print_fofr_xyzri("I",nr1,nr2,nr3,n4,n5,n6,fofr,rprimd,conv_fact=Bohr_Ang,unit=unout) 2409 close(unout) 2410 exit 2411 2412 case(7) !OpenDX format, data R and data I 2413 write(std_out,*) 2414 write(std_out,*) 'Give 2 files of 3D formatted data' 2415 write(std_out,*) 'The file is ready to be use with OpenDX' 2416 write(std_out,*) 'The eig_kvalues and occupations numbers are in comments' 2417 write(std_out,*) 2418 ABI_ALLOCATE(filename,(2)) 2419 filename(1)=trim(output)//'Real.dx' 2420 filename(2)=trim(output)//'Imag.dx' 2421 write(std_out,*) ' The name of your files is : ' 2422 write(std_out,*) trim(filename(1)),' for the real part,' 2423 write(std_out,*) trim(filename(2)),' for the imaginary part.' 2424 write(std_out,*) 2425 2426 do ifile=1,2 2427 if (open_file(filename(ifile), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2428 MSG_ERROR(msg) 2429 end if 2430 rewind(unout) 2431 write(unout,*)'# band, eig_kvalues and occupations' 2432 do iband=1,nband(ckpt) 2433 write(unout,'(a,i4,2f20.16)')'#',iband,eig_k(iband),occ_k(iband) 2434 end do 2435 write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',nr1*nr2*nr3,' data follows' 2436 do ir3=1,nr3 2437 do ir2=1,nr2 2438 do ir1=1,nr1 2439 write(unout,'(f20.16)')fofr(ifile,ir1,ir2,ir3) 2440 end do 2441 end do 2442 end do 2443 2444 write(unout,'(a)')'# this is the object defining the grid connections' 2445 write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',nr3,nr2,nr1 2446 write(unout,*) 2447 write(unout,*) 2448 write(unout,'(a)')'# this is the object defining the grid' 2449 write(unout,'(a,3i5)')'object "positions" class gridpositions counts',nr3,nr2,nr1 2450 2451 write(unout,'(a)') 'origin 0 0 0' 2452 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3) 2453 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3) 2454 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3) 2455 2456 write(unout,'(a)')'# this is the collective object, one for each grid ' 2457 write(unout,'(a)')'object "densite" class field ' 2458 write(unout,'(a)')'component "positions" value "positions"' 2459 write(unout,'(a)')'component "connections" value "gridconnections" ' 2460 write(unout,'(a)')'component "data" value "donnees"' 2461 2462 close(unit=unout) 2463 end do 2464 ABI_DEALLOCATE(filename) 2465 exit 2466 2467 case(8) ! OpenDX format, data R and data I 2468 write(std_out,*) 2469 write(std_out,*) 'Give 2 files of 3D formatted data' 2470 write(std_out,*) 'The file is ready to be use with OpenDX' 2471 write(std_out,*) 'The eig_kvalues and occupations numbers are in comments' 2472 write(std_out,*) 2473 ABI_ALLOCATE(filename,(1)) 2474 filename(1)=trim(output)//'Real.dx' 2475 write(std_out,*) ' The name of your file is : ' 2476 write(std_out,*) trim(filename(1)),' for the real part,' 2477 write(std_out,*) 2478 2479 if (open_file(filename(1), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2480 MSG_ERROR(msg) 2481 end if 2482 rewind(unout) 2483 write(unout,*)'# band, eig_kvalues and occupations' 2484 do iband=1,nband(ckpt) 2485 write(unout,'(a,i4,2f20.16)')'#',iband,eig_k(iband),occ_k(iband) 2486 end do 2487 write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',nr1*nr2*nr3,' data follows' 2488 do ir3=1,nr3 2489 do ir2=1,nr2 2490 do ir1=1,nr1 2491 write(unout,'(f20.16)')fofr(1,ir1,ir2,ir3) 2492 end do 2493 end do 2494 end do 2495 2496 write(unout,'(a)')'# this is the object defining the grid connections' 2497 write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',nr3,nr2,nr1 2498 write(unout,*) 2499 write(unout,*) 2500 write(unout,'(a)')'# this is the object defining the grid' 2501 write(unout,'(a,3i5)')'object "positions" class gridpositions counts',nr3,nr2,nr1 2502 2503 write(unout,'(a)') 'origin 0 0 0' 2504 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3) 2505 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3) 2506 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3) 2507 2508 write(unout,'(a)')'# this is the collective object, one for each grid ' 2509 write(unout,'(a)')'object "densite" class field ' 2510 write(unout,'(a)')'component "positions" value "positions"' 2511 write(unout,'(a)')'component "connections" value "gridconnections" ' 2512 write(unout,'(a)')'component "data" value "donnees"' 2513 2514 close(unit=unout) 2515 ABI_DEALLOCATE(filename) 2516 exit 2517 2518 case(9) !OpenDX format, data R and data I 2519 write(std_out,*) 2520 write(std_out,*) 'Give 2 files of 3D formatted data' 2521 write(std_out,*) 'The file is ready to be use with OpenDX' 2522 write(std_out,*) 'The eig_kvalues and occupations numbers are in comments' 2523 write(std_out,*) 2524 ABI_ALLOCATE(filename,(1)) 2525 filename(1)=trim(output)//'Imag.dx' 2526 write(std_out,*) ' The name of your file is : ' 2527 write(std_out,*) trim(filename(1)),' for the imaginary part.' 2528 write(std_out,*) 2529 2530 if (open_file(filename(1), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2531 MSG_ERROR(msg) 2532 end if 2533 rewind(unout) 2534 write(unout,*)'# band, eig_kvalues and occupations' 2535 do iband=1,nband(ckpt) 2536 write(unout,'(a,i4,2f20.16)')'#',iband,eig_k(iband),occ_k(iband) 2537 end do 2538 write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',nr1*nr2*nr3,' data follows' 2539 do ir3=1,nr3 2540 do ir2=1,nr2 2541 do ir1=1,nr1 2542 write(unout,'(f20.16)')fofr(2,ir1,ir2,ir3) 2543 end do 2544 end do 2545 end do 2546 2547 write(unout,'(a)')'# this is the object defining the grid connections' 2548 write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',nr3,nr2,nr1 2549 write(unout,*) 2550 write(unout,*) 2551 write(unout,'(a)')'# this is the object defining the grid' 2552 write(unout,'(a,3i5)')'object "positions" class gridpositions counts',nr3,nr2,nr1 2553 2554 write(unout,'(a)') 'origin 0 0 0' 2555 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3) 2556 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3) 2557 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3) 2558 2559 write(unout,'(a)')'# this is the collective object, one for each grid ' 2560 write(unout,'(a)')'object "densite" class field ' 2561 write(unout,'(a)')'component "positions" value "positions"' 2562 write(unout,'(a)')'component "connections" value "gridconnections" ' 2563 write(unout,'(a)')'component "data" value "donnees"' 2564 2565 close(unit=unout) 2566 ABI_DEALLOCATE(filename) 2567 exit 2568 2569 case(10) !OpenDX format, data R and data I, atoms positions, lattice and cell 2570 write(std_out,*) 2571 write(std_out,*) 'Give 5 files of formatted data' 2572 write(std_out,*) 'The files are ready to be use with Data Explorer' 2573 write(std_out,*) 'The eig_kvalues and occupations numbers are in comments' 2574 write(std_out,*) 'of the two data files' 2575 write(std_out,*) 2576 ABI_ALLOCATE(filename,(2)) 2577 filename(1)=trim(output)//'Real.dx' 2578 filename(2)=trim(output)//'Imag.dx' 2579 write(std_out,*) ' The name of your data files is : ' 2580 write(std_out,*) trim(filename(1)),' for the real part,' 2581 write(std_out,*) trim(filename(2)),' for the imaginary part.' 2582 write(std_out,*) 2583 2584 do ifile=1,2 2585 if (open_file(filename(ifile), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2586 MSG_ERROR(msg) 2587 end if 2588 rewind(unout) 2589 do iband=1,nband(ckpt) 2590 write(unout,'(a,2f20.16)')'#', eig_k(iband),occ_k(iband) 2591 end do 2592 write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',nr1*nr2*nr3,' data follows' 2593 do ir3=1,nr3 2594 do ir2=1,nr2 2595 do ir1=1,nr1 2596 write(unout,'(f20.16)')fofr(ifile,ir1,ir2,ir3) 2597 end do 2598 end do 2599 end do 2600 2601 write(unout,'(a)')'# this is the object defining the grid connections' 2602 write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',nr3,nr2,nr1 2603 write(unout,*) 2604 write(unout,*) 2605 write(unout,'(a)')'# this is the object defining the grid' 2606 write(unout,'(a,3i5)')'object "positions" class gridpositions counts',nr3,nr2,nr1 2607 2608 write(unout,'(a)') 'origin 0 0 0' 2609 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3) 2610 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3) 2611 write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3) 2612 2613 write(unout,'(a)')'# this is the collective object, one for each grid ' 2614 write(unout,'(a)')'object "densite" class field ' 2615 write(unout,'(a)')'component "positions" value "positions"' 2616 write(unout,'(a)')'component "connections" value "gridconnections" ' 2617 write(unout,'(a)')'component "data" value "donnees"' 2618 2619 close(unit=unout) 2620 end do 2621 ABI_DEALLOCATE(filename) 2622 ! 2623 ! write LATTICE_VEC.dx file 2624 ! 2625 ABI_ALLOCATE(filename,(3)) 2626 filename(1)=trim(output1)//'_LATTICE_VEC.dx' 2627 filename(2)=trim(output1)//'_ATOM_POS.dx' 2628 filename(3)=trim(output1)//'_UCELL_FRAME.dx' 2629 write(std_out,*) 2630 write(std_out,*)'Give the lattice file, ', trim(filename(1)) 2631 if (open_file(filename(1), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2632 MSG_ERROR(msg) 2633 end if 2634 2635 write(unout,'("#",/,"#",/,"# LATTICE VECTOR INFO:",/,"#",/,"#")') 2636 write(unout,'(a)') 'object "lattices" class array type float rank 1 shape 3 items 3 data follows' 2637 do ivect=1,3 2638 write(unout,'(3f16.10)') Bohr_Ang*rprimd(1,ivect),Bohr_Ang*rprimd(2,ivect),Bohr_Ang*rprimd(3,ivect) 2639 end do 2640 write(unout,'(a,a)') 'object "lattices_location" class array type float ','rank 1 shape 3 items 3 data follows' 2641 do ivect=1,3 2642 write(unout,'(3f16.10)') 0_dp,0_dp,0_dp 2643 end do 2644 write(unout,'("object 3 class field")') 2645 write(unout,'(a)') 'component "data" value "lattices"' 2646 write(unout,'(a)') 'component "positions" value "lattices_location"' 2647 close(unout) 2648 ! 2649 ! write ATOM_POS.dx file 2650 ! 2651 write(std_out,*)'Give the atoms positions file, ', trim(filename(2)) 2652 2653 if (open_file(filename(2), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2654 MSG_ERROR(msg) 2655 end if 2656 2657 write(unout,'("#",/,"#",/,"# BALL AND STICK INFO:",/,"#",/,"#")') 2658 write(unout,'(a,i5,a)') 'object "atomcoord" array type float rank 1 shape 3 items ',natom,' data follows' 2659 do iatom=1,natom 2660 write(unout,'(3f16.10)') Bohr_Ang*xcart(1:3,iatom) 2661 end do 2662 ! write(unout,'(a,i5,a)') 'object "data" array type string rank 0 shape 2 items ',natom,' data follows' 2663 write(unout,'(a,i5,a)') 'object "colorcode" array type float rank 0 items ',natom,' data follows' 2664 do iatom=1,natom 2665 write(unout,'(f10.4)') znucl(typat(iatom)) 2666 end do 2667 write(unout,'(a)') 'object "molecule" field' 2668 write(unout,'(a)') 'component "positions" value "atomcoord"' 2669 write(unout,'(a)') 'component "data" value "colorcode"' 2670 close(unout) 2671 2672 ! 2673 ! write UCELL_FRAME.dx file 2674 ! 2675 write(std_out,*)'Give the enveloppe of the cell file, ',trim(filename(3)) 2676 if (open_file(filename(3), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2677 MSG_ERROR(msg) 2678 end if 2679 2680 write(unout,'("#",/,"#",/,"# UNIT CELL FRAME INFO:",/,"#",/,"#")') 2681 write(unout,'(a)')'object 3 class array type int rank 1 shape 2 items 12 data follows' 2682 write(unout,'(" 0 1",/," 0 2",/," 0 3",/," 1 4",/," 1 5",/," 3 5")') 2683 write(unout,'(" 3 6",/," 2 6",/," 2 4",/," 7 5",/," 7 6",/," 7 4")') 2684 write(unout,'(a)') 'attribute "element type" string "lines"' 2685 write(unout,'("object 4 class array type float rank 1 shape 3 items 8 data follows")') 2686 write(unout,'(" .00000000 .00000000 .00000000")') 2687 write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,1) 2688 write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,2) 2689 write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,3) 2690 write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,2)) 2691 write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,3)) 2692 write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,2)+rprimd(:,3)) 2693 write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,2)+rprimd(:,3)) 2694 write(unout,'("object 5 array type float rank 0 items 12 data follows")') 2695 do ivect=1,12 2696 write(unout,'("1.0")') 2697 end do 2698 write(unout,'(a)') 'attribute "dep" string "connections"' 2699 write(unout,'("object 6 class field")') 2700 write(unout,'(a)') 'component "data" value 5' 2701 write(unout,'(a)') 'component "positions" value 4' 2702 write(unout,'(a)') 'component "connections" value 3' 2703 close(unout) 2704 ABI_DEALLOCATE(filename) 2705 2706 write(std_out,*) 2707 exit 2708 2709 case(11) 2710 write(std_out,*) 2711 write(std_out,*) 'Give 1 files of formatted data' 2712 write(std_out,*) 'The files are ready to be used with XCrysDen' 2713 write(std_out,*) 2714 gridshift1 = 0 2715 gridshift2 = 0 2716 gridshift3 = 0 2717 write(std_out,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?' 2718 write(std_out,*) 2719 shift_tau(:) = 0.0 2720 if (read_string(outputchar, unit=std_in) /= 0) then 2721 MSG_ERROR("Fatal error!") 2722 end if 2723 if (outputchar == 'y' .or. outputchar == 'Y') then 2724 MSG_ERROR("Shift is buggy, don't use it") 2725 write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :' 2726 write(std_out,*) 2727 read (std_in,*) gridshift1, gridshift2, gridshift3 2728 shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1) 2729 end if 2730 2731 ABI_ALLOCATE(filename,(1)) 2732 filename(1)=trim(output) 2733 write(std_out,*) ' The name of your data files is : ' 2734 write(std_out,*) trim(filename(1)),' for the density (norm of the wfk),' 2735 write(std_out,*) 2736 2737 if (open_file(filename(1), msg, newunit=unout, status='replace',form='formatted') /= 0) then 2738 MSG_ERROR(msg) 2739 end if 2740 rewind(unout) 2741 do iband=1,nband(ckpt) 2742 write(unout,'(a,2f20.16)')'#', eig_k(iband),occ_k(iband) 2743 end do 2744 2745 write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1 2746 write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1) 2747 write(std_out,*) ' znucl = ', znucl, ' typat = ', typat, ' ntypat = ', ntypat 2748 2749 write(unout,'(1X,A)') 'DIM-GROUP' 2750 write(unout,*) '3 1' 2751 write(unout,'(1X,A)') 'PRIMVEC' 2752 do ir1 = 1,3 2753 write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3) 2754 end do 2755 write(unout,'(1X,A)') 'PRIMCOORD' 2756 write(unout,*) natom, ' 1' 2757 ! 2758 ! generate translated coordinates to match density shift 2759 ! 2760 do iatom = 1,natom 2761 tau2 (:,iatom) = xcart(:,iatom) - shift_tau(:) 2762 end do 2763 2764 do iatom = 1,natom 2765 write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), & 2766 & Bohr_Ang*tau2(1,iatom), & 2767 & Bohr_Ang*tau2(2,iatom), & 2768 & Bohr_Ang*tau2(3,iatom) 2769 end do 2770 write(unout,'(1X,A)') 'ATOMS' 2771 do iatom = 1,natom 2772 write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), & 2773 & Bohr_Ang*tau2(1,iatom), & 2774 & Bohr_Ang*tau2(2,iatom), & 2775 & Bohr_Ang*tau2(3,iatom) 2776 end do 2777 2778 ! write(unout,'(1X,A)') 'FRAMES' 2779 write(unout,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D' 2780 write(unout,*) 'datagrids' 2781 write(unout,'(1X,A)') 'DATAGRID_3D_DENSITY' 2782 write(unout,*) nr1+1,nr2+1,nr3+1 2783 write(unout,*) '0.0 0.0 0.0 ' 2784 do ir1 = 1,3 2785 write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3) 2786 end do 2787 2788 do ir3=gridshift3+1,nr3+1 2789 ii3=mod(ir3-1,nr3) + 1 2790 do ir2=gridshift2+1,nr2+1 2791 ii2=mod(ir2-1,nr2) + 1 2792 do ir1=gridshift1+1,nr1+1 2793 ii1=mod(ir1-1,nr1) + 1 2794 tmpr=fofr(1,ii1,ii2,ii3) 2795 tmpi=fofr(2,ii1,ii2,ii3) 2796 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2797 end do 2798 do ir1=1,gridshift1 2799 ii1=mod(ir1-1,nr1) + 1 2800 tmpr=fofr(1,ii1,ii2,ii3) 2801 tmpi=fofr(2,ii1,ii2,ii3) 2802 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2803 end do 2804 end do 2805 do ir2=1,gridshift2 2806 ii2=mod(ir2-1,nr2) + 1 2807 do ir1=gridshift1+1,nr1+1 2808 ii1=mod(ir1-1,nr1) + 1 2809 tmpr=fofr(1,ii1,ii2,ii3) 2810 tmpi=fofr(2,ii1,ii2,ii3) 2811 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2812 end do 2813 do ir1=1,gridshift1 2814 ii1=mod(ir1-1,nr1) + 1 2815 tmpr=fofr(1,ii1,ii2,ii3) 2816 tmpi=fofr(2,ii1,ii2,ii3) 2817 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2818 end do 2819 end do 2820 end do 2821 do ir3=1,gridshift3 2822 ii3=mod(ir3-1,nr3) + 1 2823 do ir2=gridshift2+1,nr2+1 2824 ii2=mod(ir2-1,nr2) + 1 2825 do ir1=gridshift1+1,nr1+1 2826 ii1=mod(ir1-1,nr1) + 1 2827 tmpr=fofr(1,ii1,ii2,ii3) 2828 tmpi=fofr(2,ii1,ii2,ii3) 2829 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2830 end do 2831 do ir1=1,gridshift1 2832 ii1=mod(ir1-1,nr1) + 1 2833 tmpr=fofr(1,ii1,ii2,ii3) 2834 tmpi=fofr(2,ii1,ii2,ii3) 2835 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2836 end do 2837 end do 2838 do ir2=1,gridshift2 2839 ii2=mod(ir2-1,nr2) + 1 2840 do ir1=gridshift1+1,nr1+1 2841 ii1=mod(ir1-1,nr1) + 1 2842 tmpr=fofr(1,ii1,ii2,ii3) 2843 tmpi=fofr(2,ii1,ii2,ii3) 2844 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2845 end do 2846 do ir1=1,gridshift1 2847 ii1=mod(ir1-1,nr1) + 1 2848 tmpr=fofr(1,ii1,ii2,ii3) 2849 tmpi=fofr(2,ii1,ii2,ii3) 2850 write(unout,'(e12.5)') tmpr*tmpr + tmpi*tmpi 2851 end do 2852 end do 2853 end do 2854 2855 2856 write(unout,*) 2857 write(unout,'(1X,A)') 'END_DATAGRID_3D' 2858 write(unout,'(1X,A)') 'END_BLOCK_DATAGRID3D' 2859 close(unout) 2860 2861 ABI_DEALLOCATE(filename) 2862 2863 write(std_out,*) 2864 exit 2865 2866 2867 case(12) 2868 write(std_out,*)"NetCDF output is not available anymore" 2869 exit 2870 2871 ! ************************************************************ 2872 2873 case(13) 2874 write(std_out,*) 2875 write(std_out,*) 'Give 1 files of formatted data' 2876 write(std_out,*) 'The files are ready to be used with XCrysDen' 2877 write(std_out,*) 2878 gridshift1 = 0 2879 gridshift2 = 0 2880 gridshift3 = 0 2881 write(std_out,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?' 2882 write(std_out,*) 2883 shift_tau(:) = 0.0 2884 if (read_string(outputchar, unit=std_in) /= 0) then 2885 MSG_ERROR("Fatal error!") 2886 end if 2887 if (outputchar == 'y' .or. outputchar == 'Y') then 2888 MSG_ERROR("Shift is buggy, don't use it") 2889 write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :' 2890 write(std_out,*) 2891 read (std_in,*) gridshift1, gridshift2, gridshift3 2892 shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1) 2893 end if 2894 2895 ABI_ALLOCATE(filename,(1)) 2896 filename(1)=trim(output) 2897 write(std_out,*) ' The name of your data files is : ' 2898 write(std_out,*) trim(filename(1)),' for the density (norm of the wfk),' 2899 write(std_out,*) 2900 2901 if (open_file(filename(1), msg, newunit=unout, status='unknown',form='formatted') /= 0) then 2902 MSG_ERROR(msg) 2903 end if 2904 rewind(unout) 2905 2906 do iband=1,nband(ckpt) 2907 write(unout,'(a,2f20.16)')'#', eig_k(iband),occ_k(iband) 2908 end do 2909 2910 write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1 2911 write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1) 2912 write(std_out,*) ' znucl = ', znucl, ' typat = ', typat, ' ntypat = ', ntypat 2913 2914 write(unout,'(1X,A)') 'DIM-GROUP' 2915 write(unout,*) '3 1' 2916 write(unout,'(1X,A)') 'PRIMVEC' 2917 do ir1 = 1,3 2918 write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3) 2919 end do 2920 write(unout,'(1X,A)') 'PRIMCOORD' 2921 write(unout,*) natom, ' 1' 2922 ! 2923 ! generate translated coordinates to match density shift 2924 ! 2925 do iatom = 1,natom 2926 tau2 (:,iatom) = xcart(:,iatom) - shift_tau(:) 2927 end do 2928 2929 do iatom = 1,natom 2930 write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), & 2931 & Bohr_Ang*tau2(1,iatom), & 2932 & Bohr_Ang*tau2(2,iatom), & 2933 & Bohr_Ang*tau2(3,iatom) 2934 end do 2935 write(unout,'(1X,A)') 'ATOMS' 2936 do iatom = 1,natom 2937 write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), & 2938 & Bohr_Ang*tau2(1,iatom), & 2939 & Bohr_Ang*tau2(2,iatom), & 2940 & Bohr_Ang*tau2(3,iatom) 2941 end do 2942 2943 ! write(unout,'(1X,A)') 'FRAMES' 2944 write(unout,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D' 2945 write(unout,*) 'datagrids' 2946 write(unout,'(1X,A)') 'DATAGRID_3D_DENSITY' 2947 write(unout,*) nr1+1,nr2+1,nr3+1 2948 write(unout,*) '0.0 0.0 0.0 ' 2949 do ir1 = 1,3 2950 write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3) 2951 end do 2952 2953 do ir3=1,nr3+1 2954 ii3=mod(ir3-1+gridshift3, nr3) + 1 2955 do ir2=1,nr2+1 2956 ii2=mod(ir2-1+gridshift2, nr2) + 1 2957 do ir1=1,nr1+1 2958 ii1=mod(ir1-1+gridshift1, nr1) + 1 2959 write(unout,'(ES17.10)') fofr(1,ii1,ii2,ii3) 2960 end do 2961 end do 2962 end do 2963 write(unout,*) 2964 write(unout,'(1X,A)') 'END_DATAGRID_3D' 2965 write(unout,'(1X,A)') 'END_BLOCK_DATAGRID3D' 2966 close(unout) 2967 2968 ABI_DEALLOCATE(filename) 2969 2970 write(std_out,*) 2971 exit 2972 2973 case(14) ! CUBE file format from GAUSSIAN 2974 2975 write(std_out,*) 2976 write(std_out,*) 'Output a cube file of 3D volumetric data' 2977 write(std_out,*) 2978 2979 if (open_file(output, msg, newunit=unout, status='replace',form='formatted') /= 0) then 2980 MSG_ERROR(msg) 2981 end if 2982 call print_fofr_cube(nr1,nr2,nr3,n4,n5,n6,fofr,rprimd,natom,znucl_atom_int,xcart,unit=unout) 2983 close(unout) 2984 exit 2985 2986 case(0) 2987 write(std_out,*)' Exit inner loop' 2988 select_exit = 1 2989 2990 case default 2991 write(std_out,*) ' This choice is not valid.' 2992 write(std_out,*) 2993 cycle 2994 2995 end select 2996 2997 end do 2998 2999 ckpt=oldckpt 3000 cband=oldcband 3001 csppol=oldcsppol 3002 cspinor=oldcspinor 3003 ! deallocate the datas 3004 ABI_DEALLOCATE(fofr) 3005 3006 write(std_out,*) ' Task ',ichoice,' has been done !' 3007 write(std_out,*) 3008 write(std_out,*) ' Run interpolation again? (1=default=yes,0=no)' 3009 read(std_in,*) iprompt 3010 if(iprompt==0) then 3011 exit 3012 else 3013 cycle 3014 end if 3015 end do 3016 3017 !Deallocate the datas 3018 ABI_DEALLOCATE(cg_k) 3019 ABI_DEALLOCATE(eig_k) 3020 ABI_DEALLOCATE(kg_dum) 3021 ABI_DEALLOCATE(ph1d) 3022 ABI_DEALLOCATE(occ_k) 3023 3024 call destroy_mpi_enreg(mpi_enreg) 3025 3026 end subroutine cut3d_wffile
m_cut3d/normalize [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
normalize
FUNCTION
Normalizes the value of v
INPUTS
v = on entry, vector to be normalized
OUTPUT
v = on exit, vector normalized
SIDE EFFECTS
v=value to be normalized
PARENTS
m_cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
435 subroutine normalize(v) 436 437 438 !This section has been created automatically by the script Abilint (TD). 439 !Do not modify the following lines by hand. 440 #undef ABI_FUNC 441 #define ABI_FUNC 'normalize' 442 !End of the abilint section 443 444 implicit none 445 446 !Arguments------------------------------------------------------------- 447 !arrays 448 real(dp),intent(inout) :: v(3) 449 450 !Local variables-------------------------------------------------------- 451 !scalars 452 integer :: idir 453 real(dp) :: norm 454 455 ! ************************************************************************* 456 457 norm=0.0 458 do idir=1,3 459 norm=norm+v(idir)**2 460 end do 461 norm=sqrt(norm) 462 463 do idir=1,3 464 v(idir)=v(idir)/norm 465 end do 466 467 end subroutine normalize
m_cut3d/reduce [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
reduce
FUNCTION
Transforms coordinates of an input point from cartesian to crystallographic
INPUTS
rcart(3)=position vector in crystallographic coordinates rprimd(3,3)=orientation of the unit cell in 3D
OUTPUT
r(3)=position vector in cartesian coordinates
PARENTS
m_cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
1000 subroutine reduce(r,rcart,rprimd) 1001 1002 1003 !This section has been created automatically by the script Abilint (TD). 1004 !Do not modify the following lines by hand. 1005 #undef ABI_FUNC 1006 #define ABI_FUNC 'reduce' 1007 !End of the abilint section 1008 1009 implicit none 1010 1011 !Arguments------------------------------------------------------------- 1012 !arrays 1013 real(dp),intent(in) :: rcart(3),rprimd(3,3) 1014 real(dp),intent(out) :: r(3) 1015 1016 !Local variables-------------------------------------------------------- 1017 !scalars 1018 !arrays 1019 real(dp) :: mminv(3,3) 1020 1021 ! ************************************************************************* 1022 1023 call matr3inv(rprimd,mminv) 1024 r(1)=rcart(1)*mminv(1,1)+rcart(2)*mminv(2,1)+rcart(3)*mminv(3,1) 1025 r(2)=rcart(1)*mminv(1,2)+rcart(2)*mminv(2,2)+rcart(3)*mminv(3,2) 1026 r(3)=rcart(1)*mminv(1,3)+rcart(2)*mminv(2,3)+rcart(3)*mminv(3,3) 1027 1028 end subroutine reduce
m_cut3d/vdot [ Functions ]
[ Top ] [ m_cut3d ] [ Functions ]
NAME
vdot
FUNCTION
Computes the cross product of two vectors
INPUTS
x1(3)=first vector x2(3)=second vector
OUTPUT
x3(3)=cross product of x1 * x2
PARENTS
m_cut3d
CHILDREN
cg_getspin,dens_in_sph,destroy_distribfft,destroy_mpi_enreg,fourwf getkpgnorm,getph,init_distribfft_seq,initmpi_seq,initylmg,int2char10 jlspline_free,kpgio,metric,ph1d3d,print_fofr_cube,print_fofr_ri print_fofr_xyzri,recip_ylm,sort_dp,sphereboundary,splint,wfk_close wfk_open_read,wfk_read_band_block,xcart2xred
SOURCE
1154 subroutine vdot(x1,x2,x3) 1155 1156 1157 !This section has been created automatically by the script Abilint (TD). 1158 !Do not modify the following lines by hand. 1159 #undef ABI_FUNC 1160 #define ABI_FUNC 'vdot' 1161 !End of the abilint section 1162 1163 implicit none 1164 1165 !Arguments------------------------------------------------------------- 1166 !arrays 1167 real(dp),intent(in) :: x1(3),x2(3) 1168 real(dp),intent(out) :: x3(3) 1169 1170 !Local variables------------------------------- 1171 1172 ! ************************************************************************* 1173 1174 x3(1)=x1(2)*x2(3)-x2(2)*x1(3) 1175 x3(2)=x1(3)*x2(1)-x2(3)*x1(1) 1176 x3(3)=x1(1)*x2(2)-x2(1)*x1(2) 1177 1178 end subroutine vdot