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