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