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