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
|
