Branch data 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 : : !!***
|
