TABLE OF CONTENTS
- m_mkrho/atomden
- m_mkrho/initro
- m_mkrho/m_mkrho
- m_mkrho/mkrho
- m_mkrho/prtrhomxmn
- m_mkrho/read_atomden
m_mkrho/atomden [ Functions ]
[ Top ] [ m_mkrho ] [ Functions ]
NAME
atomden
FUNCTION
Construct atomic proto-bulk density (i.e. the superposed density
from neutral, isolated atoms at the bulk atomic positions).
This is useful if one wants to construct the bonding density:
rho^{bnd} = rho^{bulk}(r)
- \sum_{\alpha}\rho^{atm}_{\alpha}(r-R_{\alpha})
Where rho^{bulk} is the bulk density, rho^{atm} the atomic density
and the index \alpha sums over all atoms. the R_{\alpha} are the
atomic positions in the bulk. This routine calculates the sum over
rho^{atm}_{\alpha}(r-R_{\alpha}) on a grid.
Units are atomic.
INPUTS
calctype : type of calculation
'replace' zero the input/output density array
'add' add to the input/output density array
natom : number of atoms in cell
ntypat : number of different types of atoms in cell
typat(natom) : type of each atom
ngrid : number of gridpoints
r_vec_grid(3,ngrid) : real (non-reduced) coordinates for grid points
rho(ngrid) : input/output density array
a(3),b(3),c(3) : real-space basis vectors
atom_pos(3,natom) : reduced coordinates for atomic positions
natomgr(ntypat) : number of gridpoints for each atomic density grid
natomgrmax : max(natomgr(ntypat))
atomrgrid(natomgrmax,ntypat)
density(natomgrmax,ntypat)
OUTPUT
rho(ngrid) : input/output density array
SIDE EFFECTS
NOTES
There are two ways to compile the proto density in real space for a solid. One alternative is that the density is calculated for an extended grid encompassing the sphere of points around one atom, and the results are folded back into the unit cell. On the other hand one can, around each grid point, identify the number of atoms in a sphere equivalent to the length of the radial grid for each type of atom. The second approach, with some modification, is taken here. The numer of atoms in a supercell cell are listed such that the supercell encompasses the atoms which could contribute to any point in the grid. That list is kept and cycled through, to avoid recalculating it at each point. Note that the density calculated from the atom is the spherical average, since there is no preferred direction without any external field (and it's simpler)
SOURCE
2253 subroutine atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,atom_pos, & 2254 & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype) 2255 2256 !Arguments ------------------------------------ 2257 !scalars 2258 integer,intent(in) :: natom,ntypat,ngrid,natomgrmax,prtvol 2259 character(len=7),intent(in) :: calctype 2260 !arrays 2261 type(MPI_type),intent(in) :: MPI_enreg 2262 integer,intent(in) :: typat(natom),natomgr(ntypat) 2263 real(dp),intent(in) :: r_vec_grid(3,ngrid),a(3),b(3),c(3) 2264 real(dp),intent(in) :: atom_pos(3,natom),atomrgrid(natomgrmax,ntypat) 2265 real(dp),intent(in) :: density(natomgrmax,ntypat) 2266 real(dp),intent(inout) :: rho(ngrid) 2267 2268 !Local variables------------------------------- 2269 !scalars 2270 character(len=500) :: message 2271 integer :: cnt,delta,i,l,m,n,iatom,itypat,igrid,ncells,n_grid_p 2272 integer :: ierr,spaceComm,nprocs,master,rank,remainder 2273 real(dp) :: a_norm,b_norm,c_norm 2274 real(dp) :: r_max,R_sphere_max,dp_dummy,ybcbeg,ybcend 2275 !arrays 2276 integer :: n_equiv_atoms(ntypat),grid_index(ngrid) 2277 integer :: my_start_equiv_atoms(ntypat) 2278 integer :: my_end_equiv_atoms(ntypat) 2279 integer :: l_min(ntypat),m_min(ntypat),n_min(ntypat) 2280 integer :: l_max(ntypat),m_max(ntypat),n_max(ntypat) 2281 real(dp) :: center(3),dp_vec_dummy(3),delta_a(3),delta_b(3),delta_c(3) 2282 real(dp) :: r_atom(3),grid_distances(ngrid) 2283 integer, allocatable :: new_index(:),i_1d_dummy(:) 2284 real(dp),allocatable :: equiv_atom_dist(:,:),equiv_atom_pos(:,:,:),rho_temp(:,:) 2285 real(dp),allocatable :: dp_1d_dummy(:),dp_2d_dummy(:,:),ypp(:) 2286 real(dp),allocatable :: x_fit(:),y_fit(:) 2287 2288 2289 ! ************************************************************************ 2290 2291 !initialise and check parallel execution 2292 spaceComm=MPI_enreg%comm_cell 2293 nprocs=xmpi_comm_size(spaceComm) 2294 rank=MPI_enreg%me_kpt 2295 2296 master=0 2297 2298 !initialise variables and vectors 2299 a_norm = sqrt(dot_product(a,a)) 2300 b_norm = sqrt(dot_product(b,b)) 2301 c_norm = sqrt(dot_product(c,c)) 2302 center = (a+b+c)*half 2303 dp_dummy = dot_product(a,b)/(b_norm*b_norm) 2304 dp_vec_dummy = dp_dummy*b 2305 delta_a = a - dp_vec_dummy 2306 dp_dummy = dot_product(b,a)/(a_norm*a_norm) 2307 dp_vec_dummy = dp_dummy*a 2308 delta_b = b - dp_vec_dummy 2309 dp_dummy = dot_product(c,(a+b))/(dot_product((a+b),(a+b))) 2310 dp_vec_dummy = dp_dummy*(a+b) 2311 delta_c = c - dp_vec_dummy 2312 ABI_MALLOC(rho_temp,(ngrid,ntypat)) 2313 rho_temp = zero 2314 2315 !write(std_out,*) '*** --- In atomden --- ***' 2316 !write(std_out,*) ' a_norm:',a_norm,' b_norm:',b_norm,' c_norm:',c_norm 2317 !write(std_out,*) 'delta_a:',delta_a,'delta_b:',delta_b,'delta_c:',delta_c 2318 !write(std_out,*) ' center:',center 2319 2320 !Find supercell which will contain all possible contributions 2321 !for all atoms, and enumerate positions for all atoms 2322 !TODO list of atoms can be "pruned", i.e identify all atoms 2323 !that can't possibly contribute and remove from list. 2324 !Should be most important for very oblique cells 2325 do itypat=1,ntypat 2326 R_sphere_max = atomrgrid(natomgr(itypat),itypat) 2327 l_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_a,delta_a))) 2328 l_max(itypat) = -l_min(itypat) 2329 m_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_b,delta_b))) 2330 m_max(itypat) = -m_min(itypat) 2331 n_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_c,delta_c))) 2332 n_max(itypat) = -n_min(itypat) 2333 ncells = (l_max(itypat)-l_min(itypat)+1) & 2334 & *(m_max(itypat)-m_min(itypat)+1) & 2335 & *(n_max(itypat)-n_min(itypat)+1) 2336 n_equiv_atoms(itypat) = 0 2337 do iatom=1,natom 2338 if (typat(iatom)==itypat) then 2339 n_equiv_atoms(itypat) = n_equiv_atoms(itypat) + ncells 2340 end if ! if type=itypat 2341 end do ! number of atoms per cell 2342 if ((rank==master).and.(prtvol>9)) then 2343 write(message,'(a)') '*** --- In atomden --- find box ***' 2344 call wrtout(std_out,message,'COLL') 2345 write(message,'(a,I4)') ' itypat:',itypat 2346 call wrtout(std_out,message,'COLL') 2347 write(message,'(2(a,I4))') ' l_min:',l_min(itypat),' l_max:',l_max(itypat) 2348 call wrtout(std_out,message,'COLL') 2349 write(message,'(2(a,I4))') ' m_min:',m_min(itypat),' m_max:',m_max(itypat) 2350 call wrtout(std_out,message,'COLL') 2351 write(message,'(2(a,I4))') ' n_min:',n_min(itypat),' n_max:',n_max(itypat) 2352 call wrtout(std_out,message,'COLL') 2353 write(message,'(2(a,I4))') ' n_equiv_atoms:',n_equiv_atoms(itypat) 2354 call wrtout(std_out,message,'COLL') 2355 end if 2356 end do !atom type 2357 2358 !allocate arrays 2359 n = maxval(n_equiv_atoms) 2360 ABI_MALLOC(equiv_atom_pos,(3,n,ntypat)) 2361 ABI_MALLOC(equiv_atom_dist,(n,ntypat)) 2362 equiv_atom_pos = zero 2363 equiv_atom_dist = zero 2364 2365 !Find positions and distance of atoms from center of cell 2366 do itypat=1,ntypat 2367 i = 1 2368 do l=l_min(itypat),l_max(itypat) 2369 do m=m_min(itypat),m_max(itypat) 2370 do n=n_min(itypat),n_max(itypat) 2371 do iatom=1,natom 2372 if (typat(iatom)==itypat) then 2373 if (i>n_equiv_atoms(itypat)) then 2374 ABI_ERROR('atomden: i>n_equiv_atoms') 2375 end if 2376 equiv_atom_pos(:,i,itypat) = (atom_pos(1,iatom)+dble(l))*a & 2377 & + (atom_pos(2,iatom)+dble(m))*b & 2378 & + (atom_pos(3,iatom)+dble(n))*c 2379 dp_vec_dummy = equiv_atom_pos(:,i,itypat)-center 2380 equiv_atom_dist(i,itypat) = & 2381 & sqrt(dot_product(dp_vec_dummy,dp_vec_dummy)) 2382 i = i + 1 2383 end if 2384 end do 2385 end do !n 2386 end do !m 2387 end do !l 2388 ! write(std_out,*) '*** --- In atomden --- find equiv ***' 2389 ! write(std_out,*) ' itypat:',itypat 2390 ! write(std_out,*) ' equiv_atom_pos:' 2391 ! write(std_out,*) equiv_atom_pos(:,:,itypat) 2392 ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat) 2393 end do !atom type 2394 2395 !Sort the atoms after distance so that the density from the ones 2396 !furthest away can be added first. This is to prevent truncation error. 2397 do itypat=1,ntypat 2398 n = n_equiv_atoms(itypat) 2399 ABI_MALLOC(dp_1d_dummy,(n)) 2400 ABI_MALLOC(new_index,(n)) 2401 ABI_MALLOC(dp_2d_dummy,(3,n)) 2402 dp_1d_dummy = equiv_atom_dist(1:n,itypat) 2403 dp_2d_dummy = equiv_atom_pos(1:3,1:n,itypat) 2404 do i=1,n 2405 new_index(i) = i 2406 end do 2407 call sort_dp(n,dp_1d_dummy,new_index,tol14) 2408 do i=1,n 2409 !write(std_out,*) i,' -> ',new_index(i) 2410 equiv_atom_pos(1:3,n+1-i,itypat) = dp_2d_dummy(1:3,new_index(i)) 2411 equiv_atom_dist(1:n,itypat) = dp_1d_dummy 2412 end do 2413 ABI_FREE(dp_1d_dummy) 2414 ABI_FREE(new_index) 2415 ABI_FREE(dp_2d_dummy) 2416 ! write(std_out,*) '*** --- In atomden --- sorting atoms ***' 2417 ! write(std_out,*) ' itypat:',itypat 2418 ! write(std_out,*) ' equiv_atom_pos:' 2419 ! write(std_out,*) equiv_atom_pos(:,:,itypat) 2420 ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat) 2421 end do ! atom type 2422 2423 !Divide the work in case of parallel execution 2424 if (nprocs==1) then ! Make sure everything runs with one proc 2425 if (prtvol>9) then 2426 write(message,'(a)') ' In atomden - number of processors: 1' 2427 call wrtout(std_out,message,'COLL') 2428 write(message,'(a)') ' Calculation of proto-atomic density done in serial' 2429 call wrtout(std_out,message,'COLL') 2430 end if 2431 do itypat=1,ntypat 2432 if (prtvol>9) then 2433 write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat) 2434 call wrtout(std_out,message,'COLL') 2435 end if 2436 my_start_equiv_atoms(itypat) = 1 2437 my_end_equiv_atoms(itypat) = n_equiv_atoms(itypat) 2438 end do 2439 else 2440 if (rank==master.and.prtvol>9) then 2441 write(message,'(a,I5)') ' In atomden - number of processors:',nprocs 2442 call wrtout(std_out,message,'COLL') 2443 write(message,'(a)') ' Calculation of proto-atomic density done in parallel' 2444 call wrtout(std_out,message,'COLL') 2445 end if 2446 do itypat=1,ntypat 2447 if (rank==master.and.prtvol>9) then 2448 write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat) 2449 call wrtout(std_out,message,'COLL') 2450 end if 2451 ! Divide the atoms among the processors by shuffling indices 2452 delta = int(floor(real(n_equiv_atoms(itypat))/real(nprocs))) 2453 remainder = n_equiv_atoms(itypat)-nprocs*delta 2454 my_start_equiv_atoms(itypat) = 1+rank*delta 2455 my_end_equiv_atoms(itypat) = (rank+1)*delta 2456 ! Divide the remainder points among the processors 2457 ! by shuffling indices 2458 if ((rank+1)>remainder) then 2459 my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + remainder 2460 my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + remainder 2461 else 2462 my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + rank 2463 my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + rank + 1 2464 end if 2465 if (prtvol>9) then 2466 write(message,'(a,I3)') ' For atom type: ',itypat 2467 call wrtout(std_out,message,'PERS') 2468 ! write(message,'(a,I6)') ' I''ll take atoms from: ',my_start_equiv_atoms(itypat) 2469 ! call wrtout(std_out,message,'PERS') 2470 ! write(message,'(a,I6)') ' total for me: ',my_end_equiv_atoms(itypat) 2471 ! call wrtout(std_out,message,'PERS') 2472 write(message,'(a,I6)') ' total for me: ', & 2473 & my_end_equiv_atoms(itypat)+1-my_start_equiv_atoms(itypat) 2474 call wrtout(std_out,message,'PERS') 2475 end if 2476 end do 2477 end if 2478 2479 !Loop over types of atoms and equivalent atoms and 2480 !interpolate density onto grid 2481 do itypat=1,ntypat 2482 ! do iatom=my_start_equiv_atoms(itypat),my_end_equiv_atoms(itypat) 2483 2484 cnt = 0 2485 iatom = rank+1 - nprocs 2486 ! Parallel execution of loop 2487 do 2488 cnt = cnt + 1 2489 iatom = iatom + nprocs 2490 if (iatom>n_equiv_atoms(itypat)) exit ! Exit if index is too large 2491 2492 if (mod(cnt,100)==0.and.prtvol>0) then 2493 write(message,'(2(a,I6))') ' atoms so far',cnt,' of: ',n_equiv_atoms(itypat)/nprocs 2494 call wrtout(std_out,message,'PERS') 2495 end if 2496 2497 r_max = atomrgrid(natomgr(itypat),itypat) 2498 r_atom = equiv_atom_pos(:,iatom,itypat) 2499 2500 ! Set up an array with the gridpoint distances 2501 i = 1 2502 grid_distances = zero 2503 grid_index = 0 2504 do igrid=1,ngrid 2505 dp_vec_dummy(:) = r_vec_grid(:,igrid) - r_atom(:) 2506 dp_dummy = sqrt(dot_product(dp_vec_dummy,dp_vec_dummy)) 2507 if (dp_dummy <= r_max) then 2508 grid_distances(i) = dp_dummy 2509 grid_index(i) = igrid 2510 i = i + 1 2511 else 2512 cycle ! cycle if point is too far away 2513 end if 2514 end do 2515 n_grid_p = i - 1 2516 2517 if (n_grid_p==0) cycle ! Cycle if no point needs 2518 ! to be interpolated 2519 2520 ! Sort points to be interpolated in ascending order 2521 ABI_MALLOC(dp_1d_dummy,(n_grid_p)) 2522 ABI_MALLOC(new_index,(n_grid_p)) 2523 ABI_MALLOC(i_1d_dummy,(n_grid_p)) 2524 dp_1d_dummy = grid_distances(1:n_grid_p) 2525 do i=1,n_grid_p 2526 new_index(i) = i 2527 end do 2528 call sort_dp(n_grid_p,dp_1d_dummy,new_index,tol16) 2529 grid_distances(1:n_grid_p) = dp_1d_dummy 2530 i_1d_dummy = grid_index(1:n_grid_p) 2531 do i=1,n_grid_p 2532 !write(std_out,*) i_1d_dummy(i),' -> ',i_1d_dummy(new_index(i)) 2533 grid_index(i) = i_1d_dummy(new_index(i)) 2534 end do 2535 ABI_FREE(dp_1d_dummy) 2536 ABI_FREE(new_index) 2537 ABI_FREE(i_1d_dummy) 2538 2539 ! Interpolate density onto all grid points 2540 ABI_MALLOC(ypp,(natomgr(itypat))) 2541 ABI_MALLOC(x_fit,(n_grid_p)) 2542 ABI_MALLOC(y_fit,(n_grid_p)) 2543 ypp = zero; y_fit = zero 2544 ybcbeg = zero; ybcend = zero 2545 x_fit = grid_distances(1:n_grid_p) 2546 call spline(atomrgrid(1:natomgr(itypat),itypat), & 2547 & density(1:natomgr(itypat),itypat), & 2548 & natomgr(itypat),ybcbeg,ybcend,ypp) 2549 call splint(natomgr(itypat),atomrgrid(1:natomgr(itypat),itypat), & 2550 & density(1:natomgr(itypat),itypat),ypp,n_grid_p, & 2551 & x_fit,y_fit) 2552 2553 ! Save the interpolated points to grid 2554 do i=1,n_grid_p 2555 rho_temp(grid_index(i),itypat) = rho_temp(grid_index(i),itypat) + y_fit(i) 2556 end do 2557 ABI_FREE(ypp) 2558 ABI_FREE(x_fit) 2559 ABI_FREE(y_fit) 2560 2561 end do ! n equiv atoms 2562 end do ! type of atom 2563 2564 ! Collect all contributions to rho_temp if we are running in parallel 2565 if (nprocs>1) then 2566 call xmpi_barrier(spaceComm) 2567 call xmpi_sum_master(rho_temp,master,spaceComm,ierr) 2568 call xmpi_barrier(spaceComm) 2569 if (prtvol>9) then 2570 write(message,'(a)') ' In atomden - contributions to rho_temp collected' 2571 call wrtout(std_out,message,'PERS') 2572 end if 2573 end if 2574 2575 !Now rho_temp contains the atomic protodensity for each atom. 2576 !Check whether this is to replace or be added to the input/output array 2577 !and sum up contributions 2578 if (trim(calctype)=='replace') rho = zero 2579 do itypat=1,ntypat 2580 rho(:) = rho(:) + rho_temp(:,itypat) 2581 end do 2582 2583 ! deallocations 2584 ABI_SFREE(rho_temp) 2585 ABI_SFREE(equiv_atom_pos) 2586 ABI_SFREE(equiv_atom_dist) 2587 2588 end subroutine atomden
m_mkrho/initro [ Functions ]
[ Top ] [ m_mkrho ] [ Functions ]
NAME
initro
FUNCTION
Initialize the density using either: - a gaussian of adjustable decay length (norm-conserving psp) - PS atomic valence density from psp file (PAW or NC psps with valence change in the pp file)
INPUTS
atindx(natom)=index table for atoms (see gstate.f) densty(ntypat,4)=parameters for initialisation of the density of each atom type gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box (larger sphere). izero=if 1, unbalanced components of rho(g) have to be set to zero mgfft=maximum size of 1D FFTs mpi_enreg=information about mpi parallelization mqgrid=number of grid pts in q array for n^AT(q) spline. natom=number of atoms in cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ntypat=number of types of atoms in cell. nspden=number of spin-density components psps<type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase information for given atom coordinates. qgrid(mqgrid)=q grid for spline atomic valence density n^AT(q) from 0 to qmax. spinat(3,natom)=initial spin of each atom, in unit of hbar/2. ucvol=unit cell volume (Bohr**3). usepaw= 0 for non paw calculation; =1 for paw calculation zion(ntypat)=charge on each type of atom (real number) znucl(ntypat)=atomic number, for each type of atom
OUTPUT
rhog(2,nfft)=initialized total density in reciprocal space rhor(nfft,nspden)=initialized total density in real space. as well as spin-up part if spin-polarized
SOURCE
1004 subroutine initro(atindx,densty,gmet,gsqcut,izero,mgfft,mpi_enreg,mqgrid,natom,nattyp,& 1005 & nfft,ngfft,nspden,ntypat,psps,pawtab,ph1d,qgrid,rhog,rhor,spinat,ucvol,usepaw,zion,znucl) 1006 1007 !Arguments ------------------------------------ 1008 !scalars 1009 integer,intent(in) :: izero,mgfft,mqgrid,natom,nfft,nspden,ntypat 1010 integer,intent(in) :: usepaw 1011 real(dp),intent(in) :: gsqcut,ucvol 1012 type(mpi_type),intent(in) :: mpi_enreg 1013 type(pseudopotential_type),intent(in) :: psps 1014 !arrays 1015 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 1016 real(dp),intent(in) :: densty(ntypat,4),gmet(3,3),ph1d(2,3*(2*mgfft+1)*natom) 1017 real(dp),intent(in) :: qgrid(mqgrid),spinat(3,natom),zion(ntypat) 1018 real(dp),intent(in) :: znucl(ntypat) 1019 real(dp),intent(out) :: rhog(2,nfft),rhor(nfft,nspden) 1020 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 1021 1022 !Local variables------------------------------- 1023 !scalars 1024 integer,parameter :: im=2,re=1 1025 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,ispden 1026 integer :: itypat,jj,jtemp,me_fft,n1,n2,n3,nproc_fft 1027 real(dp),parameter :: tolfix=1.000000001_dp 1028 real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqm1,fact,fact0,gmag 1029 real(dp) :: gsquar,rhoat,sfi,sfr 1030 real(dp) :: xnorm 1031 character(len=500) :: message 1032 !arrays 1033 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:) 1034 real(dp) :: length(ntypat) 1035 real(dp),allocatable :: work(:), spinat_indx(:,:) 1036 logical :: use_gaussian(ntypat) 1037 1038 ! ************************************************************************* 1039 1040 if (nspden==4) then 1041 ABI_COMMENT('initro: might work yet for nspden=4 (not checked)') 1042 !write(std_out,*)' spinat',spinat(1:3,1:natom) 1043 end if 1044 1045 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3); me_fft=ngfft(11); nproc_fft=ngfft(10) 1046 1047 ABI_MALLOC(work,(nfft)) 1048 ABI_MALLOC(spinat_indx,(3,natom)) 1049 1050 ! Get the distrib associated with this fft_grid 1051 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1052 1053 ! Transfer the spinat array to an array in which the atoms have the proper order, type by type. 1054 do ia=1,natom 1055 spinat_indx(:,atindx(ia))=spinat(:,ia) 1056 end do 1057 1058 ! Check whether the values of spinat are acceptable 1059 if (nspden==2)then 1060 ia1=1 1061 do itypat=1,ntypat 1062 !ia1,ia2 sets range of loop over atoms: 1063 ia2=ia1+nattyp(itypat)-1 1064 do ia=ia1,ia2 1065 if( sqrt(spinat_indx(1,ia)**2+spinat_indx(2,ia)**2+spinat_indx(3,ia)**2) & 1066 > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then 1067 write(message, '(a,i0,a,a,3es11.4,a,a,a,es11.4)' )& 1068 ' For type-ordered atom number ',ia,ch10,& 1069 ' input spinat=',spinat_indx(:,ia),' is larger, in magnitude,',ch10,& 1070 ' than zion(ia)=',zion(itypat) 1071 call wrtout([std_out, ab_out], message) 1072 end if 1073 end do 1074 ia1=ia2+1 1075 end do 1076 end if 1077 1078 ! Compute the decay length of each type of atom depending on data available in pseudos. 1079 jtemp=0 1080 do itypat=1,ntypat 1081 use_gaussian(itypat)=.true. 1082 if (usepaw==0) use_gaussian(itypat) = .not. psps%nctab(itypat)%has_tvale 1083 if (usepaw==1) use_gaussian(itypat)=(pawtab(itypat)%has_tvale==0) 1084 if (.not.use_gaussian(itypat)) jtemp=jtemp+1 1085 1086 if (use_gaussian(itypat)) then 1087 length(itypat) = atom_length(densty(itypat,1),zion(itypat),znucl(itypat)) 1088 write(message,'(a,i3,a,f12.4,a,a,a,f12.4,a,i3,a,es12.4,a)' )& 1089 ' initro: for itypat=',itypat,', take decay length=',length(itypat),',',ch10,& 1090 ' initro: indeed, coreel=',znucl(itypat)-zion(itypat),', nval=',int(zion(itypat)),' and densty=',densty(itypat,1),'.' 1091 call wrtout(std_out,message) 1092 else 1093 write(message,"(a,i3,a)")' initro: for itypat=',itypat,", take pseudo charge density from pp file" 1094 call wrtout(std_out,message) 1095 end if 1096 end do 1097 1098 if (jtemp>0) then 1099 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 1100 dqm1=1.0_dp/dq 1101 dq2div6=dq**2/6.0_dp 1102 end if 1103 1104 cutoff=gsqcut*tolfix 1105 xnorm=1.0_dp/ucvol 1106 1107 id1=n1/2+2 1108 id2=n2/2+2 1109 id3=n3/2+2 1110 1111 if(nspden /= 4) then 1112 1113 do ispden=nspden,1,-1 1114 ! This loop overs spins will actually be as follows: 1115 ! ispden=2 for spin up 1116 ! ispden=1 for total spin (also valid for non-spin-polarized calculations) 1117 ! 1118 ! The reverse ispden order is chosen, in order to end up with 1119 ! rhog containing the proper total density. 1120 rhog = zero 1121 1122 ia1=1 1123 do itypat=1,ntypat 1124 if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2 1125 1126 ! ia1,ia2 sets range of loop over atoms: 1127 ia2=ia1+nattyp(itypat)-1 1128 ii=0 1129 jtemp=0 1130 1131 do i3=1,n3 1132 ig3=i3-(i3/id3)*n3-1 1133 do i2=1,n2 1134 ig2=i2-(i2/id2)*n2-1 1135 if (fftn2_distrib(i2)==me_fft) then 1136 do i1=1,n1 1137 1138 ig1=i1-(i1/id1)*n1-1 1139 ii=ii+1 1140 gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+& 1141 dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+& 1142 dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1) 1143 1144 ! Skip G**2 outside cutoff: 1145 if (gsquar<=cutoff) then 1146 1147 ! Assemble structure factor over all atoms of given type, 1148 ! also taking into account the spin-charge on each atom: 1149 1150 sfr=zero;sfi=zero 1151 if (ispden==1) then 1152 do ia=ia1,ia2 1153 sfr=sfr+phre_ini(ig1,ig2,ig3,ia) 1154 sfi=sfi-phimag_ini(ig1,ig2,ig3,ia) 1155 end do 1156 if (use_gaussian(itypat)) then 1157 sfr=sfr*zion(itypat) 1158 sfi=sfi*zion(itypat) 1159 end if 1160 else 1161 fact0=half;if (.not.use_gaussian(itypat)) fact0=half/zion(itypat) 1162 do ia=ia1,ia2 1163 ! Here, take care only of the z component 1164 fact=fact0*(zion(itypat)+spinat_indx(3,ia)) 1165 sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact 1166 sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact 1167 end do 1168 end if 1169 1170 ! Charge density integrating to one 1171 if (use_gaussian(itypat)) then 1172 rhoat=xnorm*exp(-gsquar*alf2pi2) 1173 ! Multiply structure factor times rhoat (atomic density in reciprocal space) 1174 rhog(re,ii)=rhog(re,ii)+sfr*rhoat 1175 rhog(im,ii)=rhog(im,ii)+sfi*rhoat 1176 else 1177 gmag=sqrt(gsquar) 1178 jj=1+int(gmag*dqm1) 1179 diff=gmag-qgrid(jj) 1180 bb = diff*dqm1 1181 aa = one-bb 1182 cc = aa*(aa**2-one)*dq2div6 1183 dd = bb*(bb**2-one)*dq2div6 1184 if (usepaw == 1) then 1185 rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+& 1186 & cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm 1187 else if (usepaw == 0) then 1188 rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+& 1189 cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm 1190 else 1191 ABI_BUG('Initialization of density is non consistent.') 1192 end if 1193 ! Multiply structure factor times rhoat (atomic density in reciprocal space) 1194 rhog(re,ii)=rhog(re,ii)+sfr*rhoat 1195 rhog(im,ii)=rhog(im,ii)+sfi*rhoat 1196 end if 1197 1198 else 1199 jtemp=jtemp+1 1200 end if 1201 1202 end do ! i1 1203 end if 1204 end do ! i2 1205 end do ! i3 1206 ia1=ia2+1 1207 1208 end do ! itypat 1209 1210 ! Set contribution of unbalanced components to zero 1211 if (izero==1) then 1212 call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 1213 end if 1214 !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1) 1215 1216 ! Note, we end with ispden=1, so that rhog contains the total density 1217 call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0) 1218 rhor(:,ispden)=work(:) 1219 end do ! End loop on spins 1220 1221 else if (nspden==4) then 1222 1223 do ispden=nspden,1,-1 1224 ! This loop overs spins will actually be as follows: 1225 ! ispden=2,3,4 for mx,my,mz 1226 ! ispden=1 for total spin (also valid for non-spin-polarized calculations) 1227 ! The reverse ispden order is chosen, in order to end up with 1228 ! rhog containing the proper total density. 1229 1230 rhog(:,:)=zero 1231 1232 ia1=1 1233 do itypat=1,ntypat 1234 1235 if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2 1236 1237 ! ia1,ia2 sets range of loop over atoms: 1238 ia2=ia1+nattyp(itypat)-1 1239 ii=0 1240 jtemp=0 1241 do i3=1,n3 1242 ig3=i3-(i3/id3)*n3-1 1243 do i2=1,n2 1244 ig2=i2-(i2/id2)*n2-1 1245 if (fftn2_distrib(i2)==me_fft) then 1246 do i1=1,n1 1247 1248 ig1=i1-(i1/id1)*n1-1 1249 ii=ii+1 1250 gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+& 1251 dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+& 1252 dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1) 1253 1254 ! Skip G**2 outside cutoff: 1255 if (gsquar<=cutoff) then 1256 1257 ! Assemble structure factor over all atoms of given type, 1258 ! also taking into account the spin-charge on each atom: 1259 sfr=zero;sfi=zero 1260 if(ispden==1)then 1261 do ia=ia1,ia2 1262 sfr=sfr+phre_ini(ig1,ig2,ig3,ia) 1263 sfi=sfi-phimag_ini(ig1,ig2,ig3,ia) 1264 end do 1265 if (use_gaussian(itypat)) then 1266 sfr=sfr*zion(itypat) 1267 sfi=sfi*zion(itypat) 1268 end if 1269 else 1270 fact0=one;if (.not.use_gaussian(itypat)) fact0=one/zion(itypat) 1271 do ia=ia1,ia2 1272 ! Here, take care of the components of m 1273 fact=fact0*spinat_indx(ispden-1,ia) 1274 sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact 1275 sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact 1276 end do 1277 end if 1278 1279 ! Charge density integrating to one 1280 if (use_gaussian(itypat)) then 1281 rhoat=xnorm*exp(-gsquar*alf2pi2) 1282 else 1283 gmag=sqrt(gsquar) 1284 jj=1+int(gmag*dqm1) 1285 diff=gmag-qgrid(jj) 1286 bb = diff*dqm1 1287 aa = one-bb 1288 cc = aa*(aa**2-one)*dq2div6 1289 dd = bb*(bb**2-one)*dq2div6 1290 if (usepaw == 1) then 1291 rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+& 1292 & cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm 1293 else if (usepaw == 0) then 1294 rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+& 1295 cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm 1296 else 1297 ABI_BUG('Initialization of density is non consistent.') 1298 end if 1299 end if 1300 1301 ! Multiply structure factor times rhoat (atomic density in reciprocal space) 1302 rhog(re,ii)=rhog(re,ii)+sfr*rhoat 1303 rhog(im,ii)=rhog(im,ii)+sfi*rhoat 1304 else 1305 jtemp=jtemp+1 1306 end if 1307 1308 end do ! i1 1309 end if 1310 end do ! i2 1311 end do ! i3 1312 ia1=ia2+1 1313 end do ! itypat 1314 1315 ! Set contribution of unbalanced components to zero 1316 if (izero==1) then 1317 call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 1318 end if 1319 !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1) 1320 1321 ! Note, we end with ispden=1, so that rhog contains the total density 1322 call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0) 1323 rhor(:,ispden)=work(:) 1324 end do ! ispden 1325 1326 ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities 1327 ! Add a small real to the magnetization 1328 if (all(abs(spinat(:,:))<tol10)) rhor(:,4)=rhor(:,4)+tol14 1329 1330 end if ! nspden==4 1331 1332 ABI_FREE(spinat_indx) 1333 ABI_FREE(work) 1334 1335 contains 1336 1337 !Real and imaginary parts of phase. 1338 function phr_ini(x1,y1,x2,y2,x3,y3) 1339 1340 real(dp) :: phr_ini 1341 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1342 phr_ini=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 1343 end function phr_ini 1344 1345 function phi_ini(x1,y1,x2,y2,x3,y3) 1346 1347 real(dp) :: phi_ini 1348 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1349 phi_ini=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 1350 end function phi_ini 1351 1352 function ph1_ini(nri,ig1,ia) 1353 1354 real(dp) :: ph1_ini 1355 integer,intent(in) :: nri,ig1,ia 1356 ph1_ini=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 1357 end function ph1_ini 1358 1359 function ph2_ini(nri,ig2,ia) 1360 1361 real(dp) :: ph2_ini 1362 integer,intent(in) :: nri,ig2,ia 1363 ph2_ini=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 1364 end function ph2_ini 1365 1366 function ph3_ini(nri,ig3,ia) 1367 1368 real(dp) :: ph3_ini 1369 integer,intent(in) :: nri,ig3,ia 1370 ph3_ini=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 1371 end function ph3_ini 1372 1373 function phre_ini(ig1,ig2,ig3,ia) 1374 1375 real(dp) :: phre_ini 1376 integer,intent(in) :: ig1,ig2,ig3,ia 1377 phre_ini=phr_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),& 1378 & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia)) 1379 end function phre_ini 1380 1381 function phimag_ini(ig1,ig2,ig3,ia) 1382 1383 real(dp) :: phimag_ini 1384 integer,intent(in) :: ig1,ig2,ig3,ia 1385 phimag_ini=phi_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),& 1386 & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia)) 1387 end function phimag_ini 1388 1389 end subroutine initro
m_mkrho/m_mkrho [ Modules ]
NAME
m_mkrho
FUNCTION
Procedures for computing densities from KS orbitals.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, LSI, AR, MB, MT, SM, VR, FJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mkrho 23 24 use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64 25 26 use defs_basis 27 use defs_wvltypes 28 use m_abicore 29 use m_xmpi 30 use m_errors 31 use m_dtset 32 use m_extfpmd 33 34 use defs_abitypes, only : MPI_type 35 use m_fstrings, only : sjoin, itoa 36 use m_time, only : timab 37 use m_fftcore, only : sphereboundary 38 use m_fft, only : fftpac, zerosym, fourwf, fourdp 39 use m_bandfft_kpt, only : bandfft_kpt_set_ikpt 40 use m_paw_dmft, only : paw_dmft_type 41 use m_spacepar, only : symrhg 42 use defs_datatypes, only : pseudopotential_type 43 use m_atomdata, only : atom_length 44 use m_mpinfo, only : ptabs_fourdp, proc_distrb_cycle 45 use m_pawtab, only : pawtab_type 46 use m_io_tools, only : open_file 47 use m_splines, only : spline, splint 48 use m_sort, only : sort_dp 49 use m_prep_kgb, only : prep_fourwf 50 use m_wvl_rho, only : wvl_mkrho 51 use m_rot_cg, only : rot_cg 52 53 #if defined HAVE_YAKL 54 use gator_mod 55 #endif 56 57 #ifdef HAVE_OPENMP_OFFLOAD 58 use m_ompgpu_fourwf 59 #endif 60 61 implicit none 62 63 private
m_mkrho/mkrho [ Functions ]
[ Top ] [ m_mkrho ] [ Functions ]
NAME
mkrho
FUNCTION
Depending on option argument value: --Compute charge density rho(r) and rho(G) in electrons/bohr**3 from input wavefunctions, band occupations, and k point wts. --Compute kinetic energy density tau(r) and tau(G) in bohr**-5 from input wavefunctions, band occupations, and k point wts. --Compute a given element of the kinetic energy density tensor tau_{alpha,beta}(r) and tau_{alpha,beta}(G) in bohr**-5 from input wavefunctions, band occupations, and k point wts.
INPUTS
cg(2,mcg)=wf in G space dtset <type(dataset_type)>=all input variables for this dataset | istwfk(nkpt)=input option parameter that describes the storage of wfs | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem=Number of k points treated by this node | mpw=maximum allowed value for npw | nband(nkpt*nsppol)=number of bands to be included in summation | at each k point for each spin channel | nfft=(effective) number of FFT grid points (for this processor) | ngfft(18)=contain all needed information about 3D FFT, | see ~abinit/doc/variables/vargs.htm#ngfft | nkpt=number of k points | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetry elements in group (at least 1 for identity) | symafm(nsym)=(anti)ferromagnetic part of symmetry operations | symrel(3,3,nsym)=symmetry matrices in real space (integers) | wtk(nkpt)=k point weights (they sum to 1.0) gprimd(3,3)=dimensional reciprocal space primitive translations irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data kg(3,mpw*mkmem)=reduced planewave coordinates mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization npwarr(nkpt)=number of planewaves and boundary planewaves at each k occ(mband*nkpt*nsppol)= occupation numbers for each band (usually 2.0) at each k point option if 0: compute rhor (electron density) if 1: compute taur (kinetic energy density) (i.e. Trace over the kinetic energy density tensor) if 2: compute taur_{alpha,beta} !!NOT YET IMPLEMENTED (a given element of the kinetic energy density tensor) paw_dmft <type(paw_dmft_type)>= paw+dmft related data phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases rprimd(3,3)=dimensional real space primitive translations tim_mkrho=timing code of the calling routine(can be set to 0 if not attributed) ucvol=unit cell volume (Bohr**3) wvl_den <type(wvl_denspot_type)>=density information for wavelets wvl_wfs <type(wvl_projector_type)>=wavefunctions information for wavelets
OUTPUT
rhog(2,nfft)=total electron density in G space rhor(nfft,nspden)=electron density in r space (if spin polarized, array contains total density in first half and spin-up density in second half) (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2)
SOURCE
138 subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 139 & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wvl_wfs,& 140 & extfpmd,option) !optional 141 142 !Arguments ------------------------------------ 143 !scalars 144 integer,intent(in) :: mcg,tim_mkrho 145 integer,intent(in),optional :: option 146 real(dp),intent(in) :: ucvol 147 type(extfpmd_type),intent(in),pointer,optional :: extfpmd 148 type(MPI_type),intent(inout) :: mpi_enreg 149 type(dataset_type),intent(in) :: dtset 150 type(paw_dmft_type), intent(in) :: paw_dmft 151 type(wvl_wf_type),intent(inout) :: wvl_wfs 152 type(wvl_denspot_type), intent(inout) :: wvl_den 153 !no_abirules 154 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 155 integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2, & 156 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 157 integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 158 real(dp), intent(in) :: gprimd(3,3) 159 real(dp), intent(in) :: cg(2,mcg) 160 real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 161 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 162 real(dp), intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym), & 163 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 164 real(dp), intent(in) :: rprimd(3,3) 165 real(dp), intent(out) :: rhor(dtset%nfft,dtset%nspden),rhog(2,dtset%nfft) 166 167 !Local variables------------------------------- 168 !scalars 169 integer,save :: nskip=0 170 integer :: alpha,use_nondiag_occup_dmft,bdtot_index,beta,blocksize,iband,iband1,ibandc1,ib,iblock,icg,ierr 171 integer :: ifft,ikg,ikpt,ioption,ipw,ipwbd,ipwsp,ishf,ispden,ispinor,ispinor_index 172 integer :: isppol,istwf_k,jspinor_index 173 integer :: me,my_nspinor,n1,n2,n3,n4,n5,n6,nalpha,nband_k,nband_occ,nbandc1,nbdblock,nbeta 174 integer :: ndat,nfftot,npw_k,spaceComm,tim_fourwf 175 integer :: iband_me 176 integer :: mband_mem 177 real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp 178 real(dp) :: weight,weight_i 179 !character(len=500) :: message 180 !arrays 181 integer,allocatable :: gbound(:,:) 182 integer, ABI_CONTIGUOUS pointer :: kg_k(:,:) => null() 183 logical :: locc_test,nspinor1TreatedByThisProc,nspinor2TreatedByThisProc 184 real(dp) :: dummy(2,1) = reshape( (/0.0, 0.0/), shape(dummy)) 185 real(dp) :: tsec(2) 186 real(dp),allocatable :: cwavef_rot(:,:,:,:),occ_diag(:),occ_k(:) 187 real(dp),allocatable :: kg_k_cart_block(:),taur_alphabeta(:,:,:,:),weight_t(:) 188 real(dp), ABI_CONTIGUOUS pointer :: cwavef(:,:,:) => null() 189 real(dp), ABI_CONTIGUOUS pointer :: cwavefb(:,:,:) => null() 190 real(dp), ABI_CONTIGUOUS pointer :: cwavef_x(:,:) => null() 191 real(dp), ABI_CONTIGUOUS pointer :: cwavef_y(:,:) => null() 192 real(dp), ABI_CONTIGUOUS pointer :: cwavefb_x(:,:) => null() ! only use when paral_kgb=0 193 real(dp), ABI_CONTIGUOUS pointer :: cwavefb_y(:,:) => null() ! only use when paral_kgb=0 194 real(dp), ABI_CONTIGUOUS pointer :: rhoaug(:,:,:) => null() 195 real(dp), ABI_CONTIGUOUS pointer :: rhoaug_down(:,:,:) => null() 196 real(dp), ABI_CONTIGUOUS pointer :: rhoaug_up(:,:,:) => null() 197 real(dp), ABI_CONTIGUOUS pointer :: rhoaug_mx(:,:,:) => null() 198 real(dp), ABI_CONTIGUOUS pointer :: rhoaug_my(:,:,:) => null() 199 real(dp), ABI_CONTIGUOUS pointer :: wfraug(:,:,:,:) => null() 200 201 ! ************************************************************************* 202 203 DBG_ENTER("COLL") 204 205 call timab(790+tim_mkrho,1,tsec) 206 call timab(799,1,tsec) 207 208 if(mpi_enreg%paralbd==0) tim_fourwf=3 209 if(mpi_enreg%paralbd==1) tim_fourwf=6 210 211 if(.not.(present(option))) then 212 ioption=0 213 else 214 ioption=option 215 end if 216 217 if(ioption/=0.and.paw_dmft%use_sc_dmft==1) then 218 ABI_ERROR('option argument value of this routines should be 0 if usedmft=1.') 219 end if 220 if(paw_dmft%use_sc_dmft/=0) then 221 nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1 222 else 223 nbandc1=1 224 end if 225 use_nondiag_occup_dmft=0 226 227 !if(dtset%nspinor==2.and.paw_dmft%use_sc_dmft==1) then 228 !write(message, '(a,a,a,a)' )ch10,& 229 !& ' mkrho : ERROR -',ch10,& 230 !& ' nspinor argument value of this routines should be 1 if usedmft=1. ' 231 !call wrtout(std_out,message,'COLL') 232 !end if 233 234 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 235 if (mpi_enreg%paral_spinor==0) then 236 ispinor_index=1;jspinor_index=1 237 nspinor1TreatedByThisProc=.true. 238 nspinor2TreatedByThisProc=(dtset%nspinor==2) 239 else 240 ispinor_index=mpi_enreg%me_spinor+1;jspinor_index=3-ispinor_index 241 nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0) 242 nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1) 243 end if 244 245 !Set local variable which depend on option argument 246 247 !nalpha*nbeta is the number of element of the kinetic energy density tensor 248 !to be computed in the irreducible Brillouin Zone (BZ) to get the result in the full BZ. 249 !In case of electron density calculation, nalpha=nbeta=1 250 select case (ioption) 251 case (0) 252 nalpha = 1 253 nbeta = 1 254 case (1) 255 nalpha = 3 256 nbeta = 1 257 ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,1)) 258 case (2) 259 nalpha = 3 260 nbeta = 3 261 ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,3)) 262 case default 263 ABI_BUG(sjoin('ioption argument value should be 0,1 or 2 while got:', itoa(ioption))) 264 end select 265 266 !Init me 267 me=mpi_enreg%me_kpt 268 !zero the charge density array in real space 269 !$OMP PARALLEL DO COLLAPSE(2) 270 do ispden=1,dtset%nspden 271 do ifft=1,dtset%nfft 272 rhor(ifft,ispden)=zero 273 end do 274 end do 275 276 !WVL - Branching with a separate mkrho procedure in wavelet. 277 if (dtset%usewvl == 1) then 278 select case(ioption) 279 case (0) 280 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den) 281 return 282 case (1) 283 !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den) 284 ABI_ERROR("kinetic energy density (taur) is not yet implemented in wavelet formalism.") 285 case (2) 286 !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den) 287 ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented in wavelet formalism.') 288 end select 289 end if 290 !WVL - Following is done in plane waves. 291 292 !start loop over alpha and beta 293 294 do alpha=1,nalpha 295 do beta=1,nbeta 296 297 ! start loop over spin and k points 298 bdtot_index=0 299 icg=0 300 301 ! n4,n5,n6 are FFT dimensions, modified to avoid cache trashing 302 n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3) 303 n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6) 304 305 ndat = 1 306 if (mpi_enreg%paral_kgb==1) then 307 ndat = mpi_enreg%bandpp 308 else if (dtset%gpu_option/=ABI_GPU_DISABLED) then 309 ndat = dtset%mband 310 end if 311 312 if (dtset%gpu_option == ABI_GPU_KOKKOS) then 313 #if defined HAVE_GPU && defined HAVE_YAKL 314 ABI_MALLOC_MANAGED(cwavef,(/2,dtset%mpw*ndat,my_nspinor/)) 315 #endif 316 else 317 if (dtset%gpu_option/=ABI_GPU_DISABLED) then 318 ABI_MALLOC(cwavef,(2,dtset%mpw*ndat,my_nspinor)) 319 else 320 ABI_MALLOC(cwavef,(2,dtset%mpw,my_nspinor)) 321 end if 322 end if 323 324 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 325 #if defined HAVE_GPU && defined HAVE_YAKL 326 ABI_MALLOC_MANAGED(rhoaug, (/n4,n5,n6/)) 327 ABI_MALLOC_MANAGED(wfraug, (/2,n4,n5,n6*ndat/)) 328 ABI_MALLOC_MANAGED(cwavefb, (/2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor/)) 329 if(dtset%nspden==4) then 330 ABI_MALLOC_MANAGED(rhoaug_up, (/n4,n5,n6/)) 331 ABI_MALLOC_MANAGED(rhoaug_down,(/n4,n5,n6/)) 332 ABI_MALLOC_MANAGED(rhoaug_mx, (/n4,n5,n6/)) 333 ABI_MALLOC_MANAGED(rhoaug_my, (/n4,n5,n6/)) 334 rhoaug_up(:,:,:)=zero 335 rhoaug_down(:,:,:)=zero 336 rhoaug_mx(:,:,:)=zero 337 rhoaug_my(:,:,:)=zero 338 end if 339 #endif 340 else 341 ABI_MALLOC(rhoaug, (n4,n5,n6)) 342 ABI_MALLOC(wfraug, (2,n4,n5,n6*ndat)) 343 ABI_MALLOC(cwavefb, (2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor)) 344 if(dtset%nspden==4) then 345 ABI_MALLOC(rhoaug_up, (n4,n5,n6)) 346 ABI_MALLOC(rhoaug_down,(n4,n5,n6)) 347 ABI_MALLOC(rhoaug_mx, (n4,n5,n6)) 348 ABI_MALLOC(rhoaug_my, (n4,n5,n6)) 349 rhoaug_up(:,:,:)=zero 350 rhoaug_down(:,:,:)=zero 351 rhoaug_mx(:,:,:)=zero 352 rhoaug_my(:,:,:)=zero 353 end if 354 end if 355 356 do isppol=1,dtset%nsppol 357 ikg=0 358 359 rhoaug(:,:,:)=zero 360 do ikpt=1,dtset%nkpt 361 362 nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 363 mband_mem = nband_k 364 if (dtset%paral_kgb==0) mband_mem = nband_k/mpi_enreg%nproc_band 365 npw_k=npwarr(ikpt) 366 istwf_k = dtset%istwfk(ikpt) 367 368 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 369 bdtot_index=bdtot_index+nband_k 370 cycle 371 end if 372 373 ABI_MALLOC(gbound,(2*dtset%mgfft+8,2)) 374 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 375 #if defined HAVE_GPU && defined HAVE_YAKL 376 ABI_MALLOC_MANAGED(kg_k, (/3,npw_k/)) 377 #endif 378 else 379 ABI_MALLOC(kg_k,(3,npw_k)) 380 end if 381 382 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 383 call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k) 384 385 ! Loop over bands to fft and square for rho(r) 386 ! Should be changed to treat bands by batch always 387 388 if (mpi_enreg%paral_kgb /= 1) then ! Not yet parallelized on spinors 389 390 #ifdef HAVE_OPENMP_OFFLOAD 391 ! With OpenMP GPU, uploading kg_k when paral_kgb==0 392 !$OMP TARGET ENTER DATA MAP(to:kg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP) 393 #endif 394 395 if (dtset%gpu_option /= ABI_GPU_DISABLED) then 396 !On GPU, treat all bands at once 397 ABI_MALLOC(weight_t,(nband_k)) 398 nband_occ = 0 399 do iband=1,nband_k 400 ipwsp = (iband-1)*npw_k*my_nspinor + icg 401 weight_t(iband) = occ(iband+bdtot_index) * dtset%wtk(ikpt)/ucvol 402 locc_test = abs(occ(iband+bdtot_index))>tol8 403 if (locc_test) then 404 nband_occ = nband_occ +1 405 ipwbd = (nband_occ-1) * npw_k 406 cwavef(:,ipwbd+1:ipwbd+npw_k,1) = cg(:,ipwsp+1:ipwsp+npw_k) 407 if (my_nspinor==2) cwavef(:,ipwbd+1:ipwbd+npw_k,2) = cg(:,ipwsp+npw_k+1:ipwsp+npw_k+npw_k) 408 if (ioption==1) then ! Multiplication by 2pi i (k+G)_alpha 409 gp2pi2 = gprimd(alpha,2)*two_pi ; gp2pi2 = gprimd(alpha,2)*two_pi ; gp2pi3 = gprimd(alpha,3)*two_pi 410 kpt_cart = gp2pi1*dtset%kptns(1,ikpt) + gp2pi2*dtset%kptns(2,ikpt) + gp2pi3*dtset%kptns(3,ikpt) 411 do ispinor=1,my_nspinor 412 do ipw=1,npw_k 413 kg_k_cart = gp2pi1*kg_k(1,ipw) + gp2pi2*kg_k(2,ipw) + gp2pi3*kg_k(3,ipw) + kpt_cart 414 cwftmp = -cwavef(2,ipwbd+ipw,ispinor)*kg_k_cart 415 cwavef(2,ipwbd+ipw,ispinor) = cwavef(1,ipwbd+ipw,ispinor)*kg_k_cart 416 cwavef(1,ipwbd+ipw,ispinor) = cwftmp 417 end do 418 end do 419 else if (ioption==2) then 420 ABI_ERROR('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.') 421 end if ! end if ioption==1 422 end if ! end if locc_test 423 end do ! end iband=1,nband_k 424 if (nband_occ>0) then 425 call fourwf(1,rhoaug,cwavef(:,1:nband_occ*npw_k,1),dummy,wfraug(:,:,:,1:n6*nband_occ),& 426 & gbound,gbound,istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,dtset%ngfft,& 427 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 428 & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),& 429 & gpu_option=dtset%gpu_option) 430 end if 431 ABI_FREE(weight_t) 432 433 else ! CPU version 434 435 iband_me = 0 436 do iband=1,nband_k 437 if(mpi_enreg%paralbd==1)then 438 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle 439 end if 440 iband_me = iband_me + 1 441 do ibandc1=1,nbandc1 ! in case of DMFT 442 ! Check if DMFT and only treat occupied states (check on occ.) 443 if(paw_dmft%use_sc_dmft == 1) then 444 iband1 = paw_dmft%include_bands(ibandc1) 445 if(paw_dmft%band_in(iband)) then 446 if(.not. paw_dmft%band_in(iband1)) stop 447 use_nondiag_occup_dmft = 1 448 locc_test = abs(paw_dmft%occnd(1,iband,iband1,ikpt,isppol)) +& 449 & abs(paw_dmft%occnd(2,iband,iband1,ikpt,isppol))>tol8 450 else 451 use_nondiag_occup_dmft = 0 452 locc_test = abs(occ(iband+bdtot_index))>tol8 453 if(ibandc1 /=1 .and. .not. paw_dmft%band_in(iband)) cycle 454 end if 455 else 456 use_nondiag_occup_dmft = 0 457 locc_test = abs(occ(iband+bdtot_index))>tol8 458 end if 459 if (locc_test) then 460 ! Obtain Fourier transform in fft box and accumulate the density or 461 ! the kinetic energy density 462 ! Not yet parallelized on nspinor if paral_kgb/=1 463 ipwsp=(iband_me-1)*npw_k*my_nspinor +icg 464 cwavef(:,1:npw_k,1) = cg(:,1+ipwsp :ipwsp+npw_k) 465 if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,1+ipwsp+npw_k:ipwsp+2*npw_k) 466 if(ioption==1)then 467 ! Multiplication by 2pi i (k+G)_alpha 468 gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi 469 kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt) 470 do ispinor=1,my_nspinor 471 do ipw=1,npw_k 472 kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart 473 cwftmp=-cwavef(2,ipw,ispinor)*kg_k_cart 474 cwavef(2,ipw,ispinor)=cwavef(1,ipw,ispinor)*kg_k_cart 475 cwavef(1,ipw,ispinor)=cwftmp 476 end do 477 end do 478 else if(ioption==2)then 479 ABI_ERROR('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.') 480 end if 481 ! Non diag occupation in DMFT. 482 ! TODO : this will break in full distrib of band memory 483 if(use_nondiag_occup_dmft==1) then 484 ipwsp=(iband1-1)*npw_k*my_nspinor +icg 485 cwavefb(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k) 486 if (my_nspinor==2) cwavefb(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k) 487 weight =paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol 488 weight_i=paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol 489 else 490 weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol 491 weight_i=weight 492 end if 493 494 ! The same section of code is also found in vtowfk.F90 : should be rationalized ! 495 496 call fourwf(1,rhoaug,cwavef(:,:,1),dummy,wfraug,gbound,gbound,& 497 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 498 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 499 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,1),& 500 & gpu_option=dtset%gpu_option) 501 if(dtset%nspinor==2)then 502 if(dtset%nspden==1) then 503 ! We need only the total density : accumulation continues on top of rhoaug 504 call fourwf(1,rhoaug,cwavef(:,:,2),dummy,wfraug,gbound,gbound,& 505 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 506 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 507 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),& 508 & gpu_option=dtset%gpu_option) 509 else if(dtset%nspden==4) then 510 ! Build the four components of rho. We use only norm quantities and, so fourwf. 511 ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$ 512 ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$ 513 ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$ 514 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 515 #if defined HAVE_GPU && defined HAVE_YAKL 516 ABI_MALLOC_MANAGED(cwavef_x, (/2,npw_k/)) 517 ABI_MALLOC_MANAGED(cwavef_y, (/2,npw_k/)) 518 ABI_MALLOC_MANAGED(cwavefb_x,(/2,npw_k*paw_dmft%use_sc_dmft/)) 519 ABI_MALLOC_MANAGED(cwavefb_y,(/2,npw_k*paw_dmft%use_sc_dmft/)) 520 #endif 521 else 522 ABI_MALLOC(cwavef_x,(2,npw_k)) 523 ABI_MALLOC(cwavef_y,(2,npw_k)) 524 ABI_MALLOC(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft)) 525 ABI_MALLOC(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft)) 526 end if 527 ! $(\Psi^{1}+\Psi^{2})$ 528 cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2) 529 ! $(\Psi^{1}-i \Psi^{2})$ 530 cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2) 531 cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2) 532 if(use_nondiag_occup_dmft==1) then 533 cwavefb_x(:,:)=cwavefb(:,1:npw_k,1)+cwavefb(:,1:npw_k,2) 534 cwavefb_y(1,:)=cwavefb(1,1:npw_k,1)+cwavefb(2,1:npw_k,2) 535 cwavefb_y(2,:)=cwavefb(2,1:npw_k,1)-cwavefb(1,1:npw_k,2) 536 end if 537 rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed 538 call fourwf(1,rhoaug_down,cwavef(:,:,2),dummy,wfraug,gbound,gbound,& 539 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 540 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 541 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),& 542 & gpu_option=dtset%gpu_option) 543 544 call fourwf(1,rhoaug_mx,cwavef_x,dummy,wfraug,gbound,gbound,& 545 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 546 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 547 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,& 548 & gpu_option=dtset%gpu_option) 549 550 call fourwf(1,rhoaug_my,cwavef_y,dummy,wfraug,gbound,gbound,& 551 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 552 & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,& 553 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,& 554 & gpu_option=dtset%gpu_option) 555 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 556 #if defined HAVE_GPU && defined HAVE_YAKL 557 ABI_FREE_MANAGED(cwavef_x) 558 ABI_FREE_MANAGED(cwavef_y) 559 ABI_FREE_MANAGED(cwavefb_x) 560 ABI_FREE_MANAGED(cwavefb_y) 561 #endif 562 else 563 ABI_FREE(cwavef_x) 564 ABI_FREE(cwavef_y) 565 ABI_FREE(cwavefb_x) 566 ABI_FREE(cwavefb_y) 567 end if 568 end if ! dtset%nspden/=4 569 end if 570 else 571 ! Accumulate the number of one-way 3D ffts skipped 572 nskip=nskip+1 573 end if ! abs(occ(iband+bdtot_index))>tol8 574 end do ! iband1=1,(nband_k-1)*paw_dmft%use_sc_dmft+1 575 end do ! iband=1,nband_k 576 577 end if ! gpu_option 578 579 #ifdef HAVE_OPENMP_OFFLOAD 580 !$OMP TARGET EXIT DATA MAP(delete:kg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP) 581 #endif 582 else !paral_kgb==1 583 584 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 585 586 call bandfft_kpt_set_ikpt(ikpt,mpi_enreg) 587 nbdblock=nband_k/(mpi_enreg%nproc_band * mpi_enreg%bandpp) 588 blocksize=nband_k/nbdblock 589 590 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 591 #if defined HAVE_GPU && defined HAVE_YAKL 592 if(associated(cwavef)) then 593 ABI_FREE_MANAGED(cwavef) 594 end if 595 ABI_MALLOC_MANAGED(cwavef,(/2,npw_k*blocksize,dtset%nspinor/)) 596 #endif 597 else 598 if(associated(cwavef)) then 599 ABI_FREE(cwavef) 600 end if 601 ABI_MALLOC(cwavef,(2,npw_k*blocksize,dtset%nspinor)) 602 end if 603 if(ioption==1) then 604 ABI_MALLOC(kg_k_cart_block,(npw_k)) 605 end if 606 ABI_MALLOC(occ_k,(nband_k)) 607 occ_k(:)=occ(bdtot_index+1:bdtot_index+nband_k) 608 609 ! ---------- DMFT 610 if(allocated(cwavef_rot)) then 611 ABI_FREE(cwavef_rot) 612 ABI_FREE(occ_diag) 613 ! ABI_FREE(occ_nd) 614 end if 615 if(paw_dmft%use_sc_dmft==1) then 616 ! Allocation of DMFT temporaries arrays 617 ABI_MALLOC(cwavef_rot,(2,npw_k,blocksize,dtset%nspinor)) 618 ABI_MALLOC(occ_diag,(blocksize)) 619 ! ABI_MALLOC(occ_nd,(2, blocksize, blocksize, dtset%nspinor)) 620 end if 621 ! ---------- END DMFT 622 623 do iblock=1,nbdblock 624 if (dtset%nspinor==1) then 625 cwavef(:,1:npw_k*blocksize,1)=cg(:,1+(iblock-1)*npw_k*blocksize+icg:iblock*npw_k*blocksize+icg) 626 else 627 if (mpi_enreg%paral_spinor==0) then 628 ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg 629 do ib=1,blocksize 630 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,1)=cg(:,1+(2*ib-2)*npw_k+ishf:(2*ib-1)*npw_k+ishf) 631 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,2)=cg(:,1+(2*ib-1)*npw_k+ishf:ib*2*npw_k+ishf) 632 end do 633 else 634 ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg 635 do ib=1,blocksize 636 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,ispinor_index)=& 637 & cg(:,1+(ib-1)*npw_k+ishf:ib*npw_k+ishf) 638 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,jspinor_index)=zero 639 end do 640 call xmpi_sum(cwavef,mpi_enreg%comm_spinor,ierr) 641 end if 642 end if 643 644 if(ioption==1)then 645 ! Multiplication by 2pi i (k+G)_alpha 646 gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi 647 kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt) 648 kg_k_cart_block(1:npw_k)=gp2pi1*kg_k(1,1:npw_k)+gp2pi2*kg_k(2,1:npw_k)+gp2pi3*kg_k(3,1:npw_k)+kpt_cart 649 do ib=1,blocksize 650 do ipw=1,npw_k 651 cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw) 652 cwavef(2,ipw,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw) 653 cwavef(1,ipw,1)=cwftmp 654 if (my_nspinor==2) then 655 cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw) 656 cwavef(2,ipw,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw) 657 cwavef(1,ipw,2)=cwftmp 658 end if 659 end do 660 end do 661 else if(ioption==2)then 662 ABI_ERROR("kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.") 663 end if 664 665 ! ---------- DMFT 666 if(paw_dmft%use_sc_dmft==1) then 667 ! initialisation of DMFT arrays 668 cwavef_rot(:,:,:,:) = zero 669 occ_diag(:) = zero 670 ! occ_nd(:,:,:,:) = paw_dmft%occnd(:,:,:,ikpt,:) 671 672 do ib=1,blocksize 673 cwavef_rot(:, :, ib, :) = cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :) 674 end do 675 676 call rot_cg(paw_dmft%occnd(:,:,:,ikpt,isppol), cwavef_rot, npw_k, nband_k, blocksize,& 677 & dtset%nspinor, paw_dmft%include_bands(1), paw_dmft%mbandc, occ_diag) 678 do ib=1,blocksize 679 cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :) = cwavef_rot(:, :, ib, :) 680 end do 681 682 occ_k(:) = occ_diag(:) 683 end if 684 ! ---------- END DMFT 685 686 call timab(538,1,tsec) 687 if (nspinor1TreatedByThisProc) then 688 call prep_fourwf(rhoaug,blocksize,cwavef(:,:,1),wfraug,iblock,istwf_k,dtset%mgfft,mpi_enreg,& 689 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 690 & dtset%wtk(ikpt),gpu_option=dtset%gpu_option) 691 end if 692 call timab(538,2,tsec) 693 if(dtset%nspinor==2)then 694 if (dtset%nspden==1) then 695 if (nspinor2TreatedByThisProc) then 696 call prep_fourwf(rhoaug,blocksize,cwavef(:,:,2),wfraug,& 697 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 698 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 699 & dtset%wtk(ikpt),gpu_option=dtset%gpu_option) 700 end if 701 else if(dtset%nspden==4 ) then 702 703 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 704 #if defined HAVE_GPU && defined HAVE_YAKL 705 ABI_MALLOC_MANAGED(cwavef_x,(/2,npw_k*blocksize/)) 706 ABI_MALLOC_MANAGED(cwavef_y,(/2,npw_k*blocksize/)) 707 #endif 708 else 709 ABI_MALLOC(cwavef_x,(2,npw_k*blocksize)) 710 ABI_MALLOC(cwavef_y,(2,npw_k*blocksize)) 711 end if 712 713 cwavef_x(:,:)=cwavef(:,:,1)+cwavef(:,:,2) 714 cwavef_y(1,:)=cwavef(1,:,1)+cwavef(2,:,2) 715 cwavef_y(2,:)=cwavef(2,:,1)-cwavef(1,:,2) 716 717 call timab(538,1,tsec) 718 if (nspinor1TreatedByThisProc) then 719 call prep_fourwf(rhoaug_down,blocksize,cwavef(:,:,2),wfraug,& 720 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 721 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 722 & dtset%wtk(ikpt),gpu_option=dtset%gpu_option) 723 end if 724 if (nspinor2TreatedByThisProc) then 725 call prep_fourwf(rhoaug_mx,blocksize,cwavef_x,wfraug,& 726 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 727 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 728 & dtset%wtk(ikpt),gpu_option=dtset%gpu_option) 729 call prep_fourwf(rhoaug_my,blocksize,cwavef_y,wfraug,& 730 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 731 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 732 & dtset%wtk(ikpt),gpu_option=dtset%gpu_option) 733 end if 734 call timab(538,2,tsec) 735 736 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 737 #if defined HAVE_GPU && defined HAVE_YAKL 738 ABI_FREE_MANAGED(cwavef_x) 739 ABI_FREE_MANAGED(cwavef_y) 740 #endif 741 else 742 ABI_FREE(cwavef_x) 743 ABI_FREE(cwavef_y) 744 end if 745 746 end if 747 end if 748 end do !iblock 749 if(ioption==1) then 750 ABI_FREE(kg_k_cart_block) 751 end if 752 753 if (associated(cwavef)) then 754 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 755 #if defined HAVE_GPU && defined HAVE_YAKL 756 ABI_FREE_MANAGED(cwavef) 757 #endif 758 else 759 ABI_FREE(cwavef) 760 end if 761 end if 762 763 ABI_FREE(occ_k) 764 end if 765 766 ABI_FREE(gbound) 767 768 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 769 #if defined HAVE_GPU && defined HAVE_YAKL 770 ABI_FREE_MANAGED(kg_k) 771 #endif 772 else 773 ABI_FREE(kg_k) 774 end if 775 776 bdtot_index=bdtot_index+nband_k 777 778 if (dtset%mkmem/=0) then 779 icg=icg+npw_k*my_nspinor*mband_mem !iband_me 780 ikg=ikg+npw_k 781 end if 782 783 end do ! ikpt 784 785 if(mpi_enreg%paral_kgb == 1) then 786 call bandfft_kpt_set_ikpt(-1,mpi_enreg) 787 if (dtset%nspden==4) then 788 ! Sum the contribution of the band and of the FFT 789 call xmpi_sum(rhoaug ,mpi_enreg%comm_bandspinorfft, ierr) 790 call xmpi_sum(rhoaug_down,mpi_enreg%comm_bandspinorfft, ierr) 791 call xmpi_sum(rhoaug_mx ,mpi_enreg%comm_bandspinorfft, ierr) 792 call xmpi_sum(rhoaug_my ,mpi_enreg%comm_bandspinorfft, ierr) 793 rhoaug_up(:,:,:) = rhoaug(:,:,:) 794 else 795 call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) 796 end if 797 end if 798 799 ! Transfer density on augmented fft grid to normal fft grid in real space 800 ! Take also into account the spin, to place it correctly in rhor. 801 if(dtset%nspden==1 .or. dtset%nspden==2) then 802 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1) 803 else if(dtset%nspden==4) then 804 ispden=1 805 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1) 806 ispden=2 807 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1) 808 ispden=3 809 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1) 810 ispden=4 811 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1) 812 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 813 #if defined HAVE_GPU && defined HAVE_YAKL 814 ABI_FREE_MANAGED(rhoaug_up) 815 ABI_FREE_MANAGED(rhoaug_down) 816 ABI_FREE_MANAGED(rhoaug_mx) 817 ABI_FREE_MANAGED(rhoaug_my) 818 #endif 819 else 820 ABI_FREE(rhoaug_up) 821 ABI_FREE(rhoaug_down) 822 ABI_FREE(rhoaug_mx) 823 ABI_FREE(rhoaug_my) 824 end if 825 end if 826 827 end do ! isppol 828 829 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 830 #if defined HAVE_GPU && defined HAVE_YAKL 831 if(associated(cwavef)) then 832 ABI_FREE_MANAGED(cwavef) 833 end if 834 if(associated(cwavefb)) then 835 ABI_FREE_MANAGED(cwavefb) 836 end if 837 ABI_FREE_MANAGED(rhoaug) 838 ABI_FREE_MANAGED(wfraug) 839 #endif 840 else 841 if(associated(cwavef)) then 842 ABI_FREE(cwavef) 843 end if 844 if(associated(cwavefb)) then 845 ABI_FREE(cwavefb) 846 endif 847 ABI_FREE(rhoaug) 848 ABI_FREE(wfraug) 849 end if 850 851 if(allocated(cwavef_rot)) then 852 ABI_FREE(cwavef_rot) 853 ABI_FREE(occ_diag) 854 ! ABI_FREE(occ_nd) 855 end if 856 857 ! Recreate full rhor on all proc. 858 call timab(48,1,tsec) 859 call timab(71,1,tsec) 860 spaceComm=mpi_enreg%comm_cell 861 if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt 862 if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt 863 call xmpi_sum(rhor,spaceComm,ierr) 864 call timab(71,2,tsec) 865 call timab(48,2,tsec) 866 867 call timab(799,2,tsec) 868 call timab(549,1,tsec) 869 870 if(ioption==1 .or. ioption==2) then 871 !$OMP PARALLEL DO COLLAPSE(2) 872 do ispden=1,dtset%nspden 873 do ifft=1,dtset%nfft 874 taur_alphabeta(ifft,ispden,alpha,beta) = rhor(ifft,ispden) 875 end do 876 end do 877 end if 878 879 end do ! beta=1,nbeta 880 end do ! alpha=1,nalpha 881 882 !Compute the trace over the kinetic energy density tensor. i.e. Sum of the 3 diagonal elements. 883 if(ioption==1)then 884 ! zero rhor array in real space 885 do ispden=1,dtset%nspden 886 !$OMP PARALLEL DO 887 do ifft=1,dtset%nfft 888 rhor(ifft,ispden)=zero 889 end do 890 end do 891 do alpha = 1, nalpha 892 !$OMP PARALLEL DO COLLAPSE(2) 893 do ispden=1,dtset%nspden 894 do ifft=1,dtset%nfft 895 rhor(ifft,ispden) = rhor(ifft,ispden) + taur_alphabeta(ifft,ispden,alpha,1) 896 end do 897 end do 898 end do 899 end if 900 901 nfftot=dtset%ngfft(1) * dtset%ngfft(2) * dtset%ngfft(3) 902 903 !Add extfpmd free electrons contribution to density 904 if(present(extfpmd)) then 905 if(associated(extfpmd)) then 906 if(extfpmd%version==1.or.extfpmd%version==2.or.extfpmd%version==3.or.extfpmd%version==4) then 907 rhor(:,:)=rhor(:,:)+extfpmd%nelect/ucvol/dtset%nspden 908 else if(extfpmd%version==10) then 909 do ispden=1,dtset%nspden 910 do ifft=1,dtset%nfft 911 rhor(ifft,ispden)=rhor(ifft,ispden)+extfpmd%nelectarr(ifft,ispden)/ucvol/dtset%nspden 912 end do 913 end do 914 end if 915 rhog(1,1)=rhog(1,1)+extfpmd%nelect/ucvol/dtset%nspden 916 end if 917 end if 918 919 select case (ioption) 920 case(0, 1) 921 call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 922 phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons) 923 if(ioption==1)then 924 !$OMP PARALLEL DO 925 do ifft=1,dtset%nfft 926 do ispden=1,dtset%nspden 927 rhor(ifft,ispden)=1.0d0/2.0d0*rhor(ifft,ispden) 928 end do 929 rhog(:,ifft)=1.0d0/2.0d0*rhog(:,ifft) 930 end do 931 end if 932 case(2) 933 ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.') 934 !call symtaug(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 935 !dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel) 936 end select 937 938 call timab(549,2,tsec) 939 940 !We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2 941 !we also have the spin-up density, symmetrized, in rhor(:,2). 942 !In case of non collinear magnetism, we have rho,mx,my,mz. No symmetry is applied 943 944 call timab(799,1,tsec) 945 946 if(ioption==1 .or. ioption==2) then 947 ABI_FREE(taur_alphabeta) 948 end if 949 950 !Find and print minimum and maximum total electron density 951 !(or total kinetic energy density, or total element of kinetic energy density tensor) and locations 952 call wrtout(std_out,' mkrho: echo density (plane-wave part only)','COLL') 953 call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol) 954 955 call timab(799,2,tsec) 956 call timab(790+tim_mkrho,2,tsec) 957 958 DBG_EXIT("COLL") 959 960 end subroutine mkrho
m_mkrho/prtrhomxmn [ Functions ]
[ Top ] [ m_mkrho ] [ Functions ]
NAME
prtrhomxmn
FUNCTION
If option==1, compute the maximum and minimum of the density (and spin-polarization if nspden==2), and print it. If option==2, also compute and print the second maximum or minimum
INPUTS
iout=unit for output file mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components option, see above optrhor=option for rhor (If optrhor==0, rhor is expected to be the electron density) (If optrhor==1, rhor is expected to be the kinetic energy density (taur)) (If optrhor==2, rhor is expected to be the gradient of the electron density (grhor)) (If optrhor==3, rhor is expected to be the laplacian of the electron density (lrhor)) (If optrhor==4, rhor is expected to be the ELF (elfr)) rhor(nfft,nspden)=electron density (electrons/bohr^3)
NOTES
The tolerance tol12 aims at giving a machine-independent ordering. (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
SOURCE
1420 subroutine prtrhomxmn(iout,mpi_enreg,nfft,ngfft,nspden,option,rhor,optrhor,ucvol) 1421 1422 !Arguments ------------------------------------ 1423 !scalars 1424 integer,intent(in) :: iout,nfft,nspden,option 1425 type(MPI_type),intent(in) :: mpi_enreg 1426 integer,intent(in),optional :: optrhor 1427 real(dp),intent(in),optional :: ucvol 1428 !arrays 1429 integer,intent(in) :: ngfft(18) 1430 real(dp),intent(in) :: rhor(nfft,nspden) 1431 1432 !Local variables------------------------------- 1433 !scalars 1434 integer :: i1,i2,i3,ierr,ifft,ii,iisign,iitems,index1,ioptrhor 1435 integer :: index2,indsign,iproc,istart,me,n1,n2,n3,nitems 1436 integer :: nfft_,nfftot,nproc,spaceComm 1437 real(dp) :: temp,value1,value2 1438 character(len=500) :: message,txt1_in_mssg,txt2_in_mssg,txt3_in_mssg 1439 logical :: reduce=.false. 1440 !arrays 1441 integer,allocatable :: iindex(:,:,:),index_fft(:,:,:,:) 1442 real(dp) :: rhomn1(4),rhomn2(4),rhomx1(4),rhomx2(4),ri_rhomn1(3,4) 1443 real(dp) :: ri_rhomn2(3,4),ri_rhomx1(3,4),ri_rhomx2(3,4),ri_zetmn1(3,2) 1444 real(dp) :: ri_zetmn2(3,2),ri_zetmx1(3,2),ri_zetmx2(3,2),zetmn1(2) 1445 real(dp) :: zetmn2(2),zetmx1(2),zetmx2(2) 1446 real(dp),allocatable :: array(:),coord(:,:,:,:),value(:,:,:),integrated(:) 1447 real(dp),allocatable :: value_fft(:,:,:) 1448 1449 ! ************************************************************************* 1450 1451 if(.not.(present(optrhor))) then 1452 ioptrhor=0 1453 else 1454 ioptrhor=optrhor 1455 end if 1456 1457 if(option/=1 .and. option/=2)then 1458 ABI_BUG(sjoin(' Option must be 1 or 2, while it is:', itoa(option))) 1459 end if 1460 1461 if (mpi_enreg%nproc_wvl>1) then 1462 ! nfft is always the potential size (in GGA, the density has buffers). 1463 nfft_ = ngfft(1) * ngfft(2) * mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 2) 1464 n1 = ngfft(1) 1465 n2 = ngfft(2) 1466 n3 = sum(mpi_enreg%nscatterarr(:, 2)) 1467 istart = mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 4) 1468 else 1469 nfft_ = nfft 1470 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1471 istart = 0 1472 end if 1473 1474 !-------------------------------------------------------------------------- 1475 !One has to determine the maximum and minimum (etc...) values 1476 !over all space, and then output it, as well as to identify 1477 !the point at which it occurs ... 1478 !This will require a bit of data exchange, and correct indirect indexing ... 1479 1480 !For the local processor, find different items : 1481 !maximum and minimum total electron density and locations 1482 !and also spin-polarisation and magnetization 1483 !also keep the second maximal or minimal value 1484 if(nspden==1)nitems=1 ! Simply the total density 1485 if(nspden==2)nitems=5 ! Total density, spin up, spin down, magnetization, zeta 1486 if(nspden==4)nitems=6 ! Total density, x, y, z, magnetization, zeta 1487 1488 ABI_MALLOC(value,(2,2,nitems)) 1489 ABI_MALLOC(iindex,(2,2,nitems)) 1490 ABI_MALLOC(array,(nfft)) 1491 ABI_MALLOC(integrated,(nitems)) 1492 1493 do iitems=1,nitems 1494 1495 ! Copy the correct values into the array 1496 ! First set of items : the density, for each spin component 1497 if(iitems<=nspden)then 1498 array(:)=rhor(:,iitems) 1499 end if 1500 ! Case nspden==2, some computation to be done 1501 if(nspden==2)then 1502 if(iitems==3)then ! Spin down 1503 array(:)=rhor(:,1)-rhor(:,2) 1504 else if(iitems==4)then ! Magnetization 1505 array(:)=2*rhor(:,2)-rhor(:,1) 1506 else if(iitems==5)then ! zeta = relative magnetization 1507 ! Avoid 0/0: the limit of (x - y) / (x+ y) depends on the direction. 1508 array(:)=zero 1509 where (abs(rhor(:,1)) > tol12) array(:)=(2*rhor(:,2)-rhor(:,1))/rhor(:,1) 1510 end if 1511 ! Case nspden==4, some other computation to be done 1512 else if(nspden==4)then 1513 if(iitems==5)then ! Magnetization 1514 array(:)=sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2) 1515 else if(iitems==6)then ! zeta = relative magnetization 1516 array(:)=(sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2))/rhor(:,1) 1517 end if 1518 end if 1519 1520 ! Zero all the absolute values that are lower than tol8, for portability reasons. 1521 do ifft = 1, nfft_ 1522 if(abs(array(ifft))<tol8)array(ifft)=zero 1523 end do 1524 1525 ! DEBUG 1526 ! write(std_out,*) ' iitems,array(1:2)=',iitems,array(1:2) 1527 ! ENDDEBUG 1528 1529 do indsign=1,2 ! Find alternatively the maximum and the minimum 1530 iisign=3-2*indsign 1531 1532 if (nfft_ > 1) then 1533 ! Initialize the two first values 1534 value1=array(istart + 1) ; value2=array(istart + 2) 1535 index1=1 ; index2=2 1536 1537 ! Ordering, if needed 1538 if( iisign*(value2+tol12) > iisign*(value1)) then 1539 temp=value2 ; value2=value1 ; value1=temp 1540 index1=2 ; index2=1 1541 end if 1542 1543 ! Integration, if relevant 1544 if(present(ucvol).and. indsign==1)then 1545 integrated(iitems) = array(istart + 1)+array(istart + 2) 1546 end if 1547 else 1548 value1 = zero; value2 = zero 1549 index1 = 0; index2 = 0 1550 end if 1551 1552 ! DEBUG 1553 ! write(std_out,*) ' value1,value2,index1,index2=',value1,value2,index1,index2 1554 ! ENDDEBUG 1555 1556 ! Loop over all points 1557 do ifft = 3, nfft_ 1558 1559 temp=array(istart + ifft) 1560 if(present(ucvol).and. indsign==1)integrated(iitems) = integrated(iitems)+temp 1561 ! Compares it to the second value 1562 if( iisign*(temp+tol12) > iisign*value2 ) then 1563 ! Compare it to the first value 1564 if( iisign*(temp+tol12) > iisign*value1 ) then 1565 value2=value1 ; index2=index1 1566 value1=temp ; index1=ifft 1567 else 1568 value2=temp ; index2=ifft 1569 end if 1570 end if 1571 1572 end do ! ifft 1573 1574 value(1,indsign,iitems)=value1 1575 value(2,indsign,iitems)=value2 1576 iindex(1,indsign,iitems)=index1 1577 iindex(2,indsign,iitems)=index2 1578 1579 ! DEBUG 1580 ! write(std_out,*) ' it,v1,i1=',iitems, value1,index1 1581 ! write(std_out,*) ' it,v2,i2=',iitems, value2,index2 1582 ! ENDDEBUG 1583 1584 end do ! indsign 1585 1586 if(present(ucvol))then 1587 nfftot=ngfft(1) * ngfft(2) * ngfft(3) 1588 integrated(iitems)=integrated(iitems)*ucvol/nfftot 1589 end if 1590 1591 ! Integrate the array 1592 ! integrated(iitems)=zero 1593 ! do ifft=1,nfft_ 1594 ! integrated(iitems) = integrated(iitems) + array(istart + ifft) 1595 ! enddo 1596 ! if(present(ucvol))integrated(iitems) = integrated(iitems)*ucvol/nfft_ 1597 ! write(std_err,*)present(ucvol) 1598 ! if(present(ucvol))then 1599 ! write(std_err,*)ucvol 1600 ! endif 1601 1602 end do ! iitems 1603 1604 ABI_FREE(array) 1605 1606 !------------------------------------------------------------------- 1607 !Enter section for FFT parallel case 1608 !if(mpi_enreg%paral_kgb>1) spaceComm=mpi_enreg%comm_fft; reduce=.true. 1609 spaceComm=mpi_enreg%comm_fft; reduce=.false. 1610 if(mpi_enreg%nproc_fft>1) then 1611 spaceComm=mpi_enreg%comm_fft; reduce=.true. 1612 else if(mpi_enreg%nproc_wvl>1) then 1613 spaceComm=mpi_enreg%comm_wvl; reduce=.true. 1614 end if 1615 nproc=xmpi_comm_size(spaceComm) 1616 me=xmpi_comm_rank(spaceComm) 1617 1618 if (reduce) then 1619 1620 ! Communicate all data to all processors with only two global communications 1621 ABI_MALLOC(value_fft,(5,nitems,nproc)) 1622 ABI_MALLOC(index_fft,(2,2,nitems,nproc)) 1623 value_fft(:,:,:)=zero 1624 index_fft(:,:,:,:)=0 1625 value_fft(1,:,me + 1)=value(1,1,:) 1626 value_fft(2,:,me + 1)=value(2,1,:) 1627 value_fft(3,:,me + 1)=value(1,2,:) 1628 value_fft(4,:,me + 1)=value(2,2,:) 1629 if(present(ucvol))value_fft(5,:,me + 1)=integrated(:) 1630 index_fft(:,:,:,me + 1)=iindex(:,:,:) 1631 call xmpi_sum(value_fft,spaceComm,ierr) 1632 call xmpi_sum(index_fft,spaceComm,ierr) 1633 1634 ! Determine the global optimum and second optimum for each item 1635 ! Also, the integrated quantities, if relevant. 1636 do iitems=1,nitems 1637 1638 if(present(ucvol))integrated(iitems)=sum(value_fft(5,iitems,1:nproc)) 1639 1640 do indsign=1,2 ! Find alternatively the maximum and the minimum 1641 iisign=3-2*indsign 1642 1643 ! Initialisation 1644 value1=value_fft(2*indsign-1,iitems,1) 1645 value2=value_fft(2*indsign ,iitems,1) 1646 index1=index_fft(1,indsign,iitems,1) 1647 index2=index_fft(2,indsign,iitems,1) 1648 1649 ! Loop 1650 do iproc=1, nproc, 1 1651 do ii=1,2 1652 if(iproc>1 .or. ii==2)then 1653 1654 temp=value_fft(ii+2*(indsign-1),iitems,iproc) 1655 ! Compares it to the second value 1656 if( iisign*(temp+tol12) > iisign*value2 ) then 1657 ! Compare it to the first value 1658 if( iisign*(temp+tol12) > iisign*value1 ) then 1659 value2=value1 ; index2=index1 1660 value1=temp ; index1=index_fft(ii,indsign,iitems,iproc) 1661 else 1662 value2=temp ; index2=index_fft(ii,indsign,iitems,iproc) 1663 end if 1664 end if 1665 1666 end if ! if(iproc>1 .or. ii==2) 1667 end do ! ii 1668 end do ! iproc 1669 1670 value(1,indsign,iitems)=value1 1671 value(2,indsign,iitems)=value2 1672 iindex(1,indsign,iitems)=index1 1673 iindex(2,indsign,iitems)=index2 1674 1675 end do ! iisign 1676 end do ! iitems 1677 1678 ABI_FREE(value_fft) 1679 ABI_FREE(index_fft) 1680 1681 end if !if(reduce) 1682 1683 !------------------------------------------------------------------- 1684 1685 !Determines the reduced coordinates of the min and max for each item 1686 ABI_MALLOC(coord,(3,2,2,nitems)) 1687 do iitems=1,nitems 1688 do indsign=1,2 1689 do ii=1,2 1690 index1=iindex(ii,indsign,iitems) 1691 i3=(index1-1)/n1/n2 1692 i2=(index1-1-i3*n1*n2)/n1 1693 i1=index1-1-i3*n1*n2-i2*n1 1694 coord(1,ii,indsign,iitems)=dble(i1)/dble(n1)+tol12 1695 coord(2,ii,indsign,iitems)=dble(i2)/dble(n2)+tol12 1696 coord(3,ii,indsign,iitems)=dble(i3)/dble(n3)+tol12 1697 ! DEBUG 1698 ! write(std_out,*)' ii,indsign,iitems,coord(1:3)=',ii,indsign,iitems,coord(:,ii,indsign,iitems) 1699 ! write(std_out,*)' value ', value(ii, indsign, iitems) 1700 ! ENDDEBUG 1701 end do 1702 end do 1703 end do 1704 1705 !------------------------------------------------------------------------- 1706 !Output 1707 if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then 1708 if(.true.)then 1709 do iitems=1,nitems 1710 1711 if(ioptrhor==4 .and. iitems>2)exit 1712 1713 select case (ioptrhor) 1714 case(0) 1715 1716 if(iitems==1) write(message,'(a)')' Total charge density [el/Bohr^3]' 1717 if(nspden==2)then 1718 if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^3]' 1719 if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^3]' 1720 if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^3]' 1721 if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1722 else if(nspden==4)then 1723 if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^3]' 1724 if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^3]' 1725 if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^3]' 1726 if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^3]' 1727 if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1728 end if 1729 1730 case(1) 1731 1732 if(iitems==1) write(message,'(a)')' Total kinetic energy density [Ha/Bohr^3]' 1733 if(nspden==2)then 1734 if(iitems==2) write(message,'(a)')' Spin up density [Ha/Bohr^3]' 1735 if(iitems==3) write(message,'(a)')' Spin down density [Ha/Bohr^3]' 1736 if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [Ha/Bohr^3]' 1737 if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1738 else if(nspden==4)then 1739 if(iitems==2) write(message,'(a)')' x component of magnetization [Ha/Bohr^3]' 1740 if(iitems==3) write(message,'(a)')' y component of magnetization [Ha/Bohr^3]' 1741 if(iitems==4) write(message,'(a)')' z component of magnetization [Ha/Bohr^3]' 1742 if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [Ha/Bohr^3]' 1743 if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1744 end if 1745 1746 case(2) 1747 1748 if(iitems==1) write(message,'(a)')' Gradient of the electronic density [el/Bohr^4]' 1749 if(nspden==2)then 1750 if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^4]' 1751 if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^4]' 1752 if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]' 1753 if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1754 else if(nspden==4)then 1755 if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]' 1756 if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]' 1757 if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]' 1758 if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^4]' 1759 if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1760 end if 1761 1762 case(3) 1763 1764 if(iitems==1) write(message,'(a)')' Laplacian of the electronic density [el/Bohr^5]' 1765 if(nspden==2)then 1766 if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^5]' 1767 if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^5]' 1768 if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^5]' 1769 if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1770 else if(nspden==4)then 1771 if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^5]' 1772 if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^5]' 1773 if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^5]' 1774 if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^5]' 1775 if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1776 end if 1777 1778 case(4) 1779 1780 if(iitems==1) write(message,'(a)')' Electron Localization Function (ELF) [min:0;max:1]' 1781 if(nspden==2)then 1782 if(iitems==2) write(message,'(a)')' Spin up ELF [min:0;max:1]' 1783 ! if(iitems==3) write(message,'(a)')' Spin down ELF [min:0;max:1]' 1784 ! if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]' 1785 ! if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1786 else if(nspden==4)then 1787 ! if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]' 1788 ! if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]' 1789 ! if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]' 1790 ! if(iitems==5) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]' 1791 ! if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) ' 1792 end if 1793 end select 1794 1795 call wrtout(iout,message,'COLL') 1796 1797 write(message,'(a,es13.4,a,3f10.4)') ') Maximum= ',& 1798 & value(1,1,iitems),' at reduced coord.',coord(:,1,1,iitems) 1799 call wrtout(iout,message,'COLL') 1800 if(option==2)then 1801 write(message,'(a,es13.4,a,3f10.4)') ')Next maximum= ',& 1802 & value(2,1,iitems),' at reduced coord.',coord(:,2,1,iitems) 1803 call wrtout(iout,message,'COLL') 1804 end if 1805 write(message,'(a,es13.4,a,3f10.4)') ') Minimum= ',& 1806 & value(1,2,iitems),' at reduced coord.',coord(:,1,2,iitems) 1807 call wrtout(iout,message,'COLL') 1808 if(option==2)then 1809 write(message,'(a,es13.4,a,3f10.4)') ')Next minimum= ',& 1810 & value(2,2,iitems),' at reduced coord.',coord(:,2,2,iitems) 1811 call wrtout(iout,message,'COLL') 1812 end if 1813 if(present(ucvol))then 1814 if(.not.(nspden==2.and.iitems==5) .and. .not.(nspden==4.and.iitems==6))then 1815 if(abs(integrated(iitems))<tol10)integrated(iitems)=zero 1816 write(message,'(a,es13.4)')' Integrated= ',integrated(iitems) 1817 call wrtout(iout,message,'COLL') 1818 end if 1819 end if 1820 1821 end do ! iitems 1822 end if 1823 1824 if(.false.)then 1825 1826 select case(optrhor) 1827 case(0) 1828 write(txt1_in_mssg, '(a)')" Min el dens=" 1829 write(txt2_in_mssg, '(a)')" el/bohr^3 at reduced coord." 1830 write(txt3_in_mssg, '(a)')" Max el dens=" 1831 case(1) 1832 write(txt1_in_mssg, '(a)')" Min kin energy dens=" 1833 write(txt2_in_mssg, '(a)')" bohr^(-5) at reduced coord." 1834 write(txt3_in_mssg, '(a)')" Max kin energy dens=" 1835 end select 1836 1837 write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,& 1838 & trim(txt1_in_mssg),value(1,2,1),& 1839 & trim(txt2_in_mssg),coord(:,1,2,1) 1840 call wrtout(iout,message,'COLL') 1841 if(option==2)then 1842 write(message, '(a,1p,e12.4,a,0p,3f8.4)' ) & 1843 & ', next min=',value(2,2,1),& 1844 & trim(txt2_in_mssg),coord(:,2,2,1) 1845 call wrtout(iout,message,'COLL') 1846 end if 1847 write(message, '(a,1p,e12.4,a,0p,3f8.4)' )& 1848 & trim(txt3_in_mssg),value(1,1,1),& 1849 & trim(txt2_in_mssg),coord(:,1,1,1) 1850 call wrtout(iout,message,'COLL') 1851 if(option==2)then 1852 write(message, '(a,1p,e12.4,a,0p,3f8.4)' )& 1853 & ', next max=',value(2,1,1),& 1854 & trim(txt2_in_mssg),coord(:,2,1,1) 1855 call wrtout(iout,message,'COLL') 1856 end if 1857 1858 if(nspden>=2)then 1859 write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,& 1860 & ',Min spin pol zeta=',value(1,2,4+nspden/2),& 1861 & ' at reduced coord.',coord(:,1,2,4+nspden/2) 1862 call wrtout(iout,message,'COLL') 1863 if(option==2)then 1864 write(message, '(a,1p,e12.4,a,0p,3f8.4)' )& 1865 & ', next min=',value(2,2,4+nspden/2),& 1866 & ' at reduced coord.',coord(:,2,2,4+nspden/2) 1867 call wrtout(iout,message,'COLL') 1868 end if 1869 write(message, '(a,1p,e12.4,a,0p,3f8.4)' )& 1870 & ',Max spin pol zeta=',value(1,1,4+nspden/2),& 1871 & ' at reduced coord.',coord(:,1,1,4+nspden/2) 1872 call wrtout(iout,message,'COLL') 1873 if(option==2)then 1874 write(message, '(a,1p,e12.4,a,0p,3f8.4)' )& 1875 & ', next max=',value(2,1,4+nspden/2),& 1876 & ' at reduced coord.',coord(:,2,1,4+nspden/2) 1877 call wrtout(iout,message,'COLL') 1878 end if 1879 end if ! nspden 1880 1881 end if ! second section always true 1882 1883 if(nspden==2 .and. .false.)then 1884 write(message,'(a)')& 1885 & ' Position in reduced coord. ( x y z )' 1886 call wrtout(iout,message,'COLL') 1887 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Total el-den) : [el/Bohr^3]',& 1888 & rhomn1(1),' at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1) 1889 call wrtout(iout,message,'COLL') 1890 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin-up den) : [el/Bohr^3]',& 1891 & rhomn1(2),' at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2) 1892 call wrtout(iout,message,'COLL') 1893 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin-down den) : [el/Bohr^3]',& 1894 & zetmn1(1),' at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1) 1895 call wrtout(iout,message,'COLL') 1896 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin pol zeta) : [m/|m|] ',& 1897 & zetmn1(2),' at',ri_zetmn1(1,2),ri_zetmn1(2,2),ri_zetmn1(3,2) 1898 call wrtout(iout,message,'COLL') 1899 if(option==2)then 1900 write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Total el-den) : [el/Bohr^3]',& 1901 & rhomn2(1),' at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1) 1902 call wrtout(iout,message,'COLL') 1903 write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-up den) : [el/Bohr^3]',& 1904 & rhomn2(2),' at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2) 1905 call wrtout(iout,message,'COLL') 1906 write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-down den) : [el/Bohr^3]',& 1907 & zetmn2(1),' at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1) 1908 call wrtout(iout,message,'COLL') 1909 write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin pol zeta) : [m/|m|] ',& 1910 & zetmn2(2),' at',ri_zetmn2(1,2),ri_zetmn2(2,2),ri_zetmn2(3,2) 1911 call wrtout(iout,message,'COLL') 1912 end if 1913 write(message,*)' ' 1914 call wrtout(iout,message,'COLL') 1915 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Total el-den) : [el/Bohr^3]',& 1916 & rhomx1(1),' at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1) 1917 call wrtout(iout,message,'COLL') 1918 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin-up den) : [el/Bohr^3]',& 1919 & rhomx1(2),' at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2) 1920 call wrtout(iout,message,'COLL') 1921 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin-down den) : [el/Bohr^3]',& 1922 & zetmx1(1),' at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1) 1923 call wrtout(iout,message,'COLL') 1924 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin pol zeta) : [m/|m|] ',& 1925 & zetmx1(2),' at',ri_zetmx1(1,2),ri_zetmx1(2,2),ri_zetmx1(3,2) 1926 call wrtout(iout,message,'COLL') 1927 if(option==2)then 1928 write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Total el-den) : [el/Bohr^3]',& 1929 & rhomx2(1),' at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1) 1930 call wrtout(iout,message,'COLL') 1931 write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-up den) : [el/Bohr^3]',& 1932 & rhomx2(2),' at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2) 1933 call wrtout(iout,message,'COLL') 1934 write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-down den) : [el/Bohr^3]',& 1935 & zetmx2(1),' at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1) 1936 call wrtout(iout,message,'COLL') 1937 write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin pol zeta) : [m/|m|] ',& 1938 & zetmx2(2),' at',ri_zetmx2(1,2),ri_zetmx2(2,2),ri_zetmx2(3,2) 1939 call wrtout(iout,message,'COLL') 1940 end if 1941 end if 1942 1943 if(nspden==4 .and. .false.)then 1944 write(message,'(a)')& 1945 & ' Position in reduced coord. ( x y z )' 1946 call wrtout(iout,message,'COLL') 1947 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Total el-den) : [el/Bohr^3]',& 1948 & rhomn1(1),' at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1) 1949 call wrtout(iout,message,'COLL') 1950 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-x) : [m/|m|] ',& 1951 & rhomn1(2),' at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2) 1952 call wrtout(iout,message,'COLL') 1953 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-y) : [m/|m|] ',& 1954 & rhomn1(3),' at',ri_rhomn1(1,3),ri_rhomn1(2,3),ri_rhomn1(3,3) 1955 call wrtout(iout,message,'COLL') 1956 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-z) : [m/|m|] ',& 1957 & rhomn1(4),' at',ri_rhomn1(1,4),ri_rhomn1(2,4),ri_rhomn1(3,4) 1958 call wrtout(iout,message,'COLL') 1959 write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin pol zeta) : [m/|m|] ',& 1960 & zetmn1(1),' at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1) 1961 call wrtout(iout,message,'COLL') 1962 if(option==2)then 1963 write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Total el-den) : [el/Bohr^3]',& 1964 & rhomn2(1),' at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1) 1965 call wrtout(iout,message,'COLL') 1966 write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-x) : [m/|m|] ',& 1967 & rhomn2(2),' at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2) 1968 call wrtout(iout,message,'COLL') 1969 write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-y) : [m/|m|] ',& 1970 & rhomn2(3),' at',ri_rhomn2(1,3),ri_rhomn2(2,3),ri_rhomn2(3,3) 1971 call wrtout(iout,message,'COLL') 1972 write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-z) : [m/|m|] ',& 1973 & rhomn2(4),' at',ri_rhomn2(1,4),ri_rhomn2(2,4),ri_rhomn2(3,4) 1974 call wrtout(iout,message,'COLL') 1975 write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Spin pol zeta) : [m/|m|] ',& 1976 & zetmn2(1),' at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1) 1977 call wrtout(iout,message,'COLL') 1978 end if 1979 write(message,*)' ' 1980 call wrtout(iout,message,'COLL') 1981 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Total el-den) : [el/Bohr^3]',& 1982 & rhomx1(1),' at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1) 1983 call wrtout(iout,message,'COLL') 1984 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-x) : [m/|m|] ',& 1985 & rhomx1(2),' at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2) 1986 call wrtout(iout,message,'COLL') 1987 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-y) : [m/|m|] ',& 1988 & rhomx1(3),' at',ri_rhomx1(1,3),ri_rhomx1(2,3),ri_rhomx1(3,3) 1989 call wrtout(iout,message,'COLL') 1990 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-z) : [m/|m|] ',& 1991 & rhomx1(4),' at',ri_rhomx1(1,4),ri_rhomx1(2,4),ri_rhomx1(3,4) 1992 call wrtout(iout,message,'COLL') 1993 write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin pol zeta) : [m/|m|] ',& 1994 & zetmx1(1),' at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1) 1995 call wrtout(iout,message,'COLL') 1996 if(option==2)then 1997 write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Total el-den) : [el/Bohr^3]',& 1998 & rhomx2(1),' at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1) 1999 call wrtout(iout,message,'COLL') 2000 write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-x) : [m/|m|] ',& 2001 & rhomx2(2),' at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2) 2002 call wrtout(iout,message,'COLL') 2003 write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-y) : [m/|m|] ',& 2004 & rhomx2(3),' at',ri_rhomx2(1,3),ri_rhomx2(2,3),ri_rhomx2(3,3) 2005 call wrtout(iout,message,'COLL') 2006 write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-z) : [m/|m|] ',& 2007 & rhomx2(4),' at',ri_rhomx2(1,4),ri_rhomx2(2,4),ri_rhomx2(3,4) 2008 call wrtout(iout,message,'COLL') 2009 write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Spin pol zeta) : [m/|m|] ',& 2010 & zetmx2(1),' at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1) 2011 call wrtout(iout,message,'COLL') 2012 end if 2013 end if 2014 end if 2015 2016 ABI_FREE(coord) 2017 ABI_FREE(value) 2018 ABI_FREE(iindex) 2019 ABI_FREE(integrated) 2020 2021 end subroutine prtrhomxmn
m_mkrho/read_atomden [ Functions ]
[ Top ] [ m_mkrho ] [ Functions ]
NAME
read_atomden
FUNCTION
INPUTS
natom : number of atoms in cell nfft=(effective) number of FFT grid points (for this processor) - fine grid ngfft(18)=contain all needed information about 3D FFT, nspden : number of spin densities ntypat : number of types of atoms in the cell typat(natom) : list of atom types
OUTPUT
rhor_atm(nfft,nspden) : full electron density on the (fine) grid
SOURCE
2043 subroutine read_atomden(MPI_enreg,natom,nfft,ngfft,nspden,ntypat, & 2044 & rhor_atm,typat,rprimd,xred,prtvol,file_prefix) 2045 2046 !Arguments ------------------------------------ 2047 !scalars 2048 integer,intent(in) :: natom,nfft,nspden,ntypat,prtvol 2049 !arrays 2050 type(MPI_type),intent(in) :: MPI_enreg 2051 integer,intent(in) :: ngfft(18),typat(natom) 2052 real(dp), intent(in) :: rprimd(3,3),xred(3,natom) 2053 real(dp),intent(inout) :: rhor_atm(nfft,nspden) 2054 character(len=7), intent(in) :: file_prefix 2055 2056 !Local variables------------------------------- 2057 !scalars 2058 character(len=500) :: message 2059 character(len=120) :: filename 2060 character(len=7) :: calctype='replace' 2061 integer :: igrid,i,i1,i2,i3,io_err,itypat,unt 2062 integer :: natomgrmax,nlines,ngrid,n1,n2,n3 2063 real(dp) :: difx,dify,difz,ucvol!,norm 2064 !arrays 2065 integer :: natomgr(ntypat) 2066 real(dp) :: a(3),b(3),c(3) 2067 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2068 real(dp),allocatable :: atomrgrid(:,:),r_vec_grid(:,:),density(:,:) 2069 real(dp),allocatable :: rho(:) 2070 2071 ! ************************************************************************ 2072 2073 !Initialise various variables 2074 ngrid = nfft 2075 a(:) = rprimd(:,1) 2076 b(:) = rprimd(:,2) 2077 c(:) = rprimd(:,3) 2078 ABI_MALLOC(rho,(ngrid)) 2079 if (nspden/=1) then 2080 ABI_ERROR('read_atomden: Only nspden=1 allowed.') 2081 end if 2082 rho = rhor_atm(:,1) 2083 gmet=zero;gprimd=zero;rmet=zero;ucvol=zero 2084 2085 2086 !Calculate the r vector (reduced coord.) of the fine gridpoints 2087 ABI_MALLOC(r_vec_grid,(3,ngrid)) 2088 igrid = 0 2089 n1 = ngfft(1) 2090 n2 = ngfft(2) 2091 n3 = ngfft(3) 2092 do i3=0,n3-1 2093 difz=dble(i3)/dble(n3) 2094 do i2=0,n2-1 2095 dify=dble(i2)/dble(n2) 2096 do i1=0,n1-1 2097 difx=dble(i1)/dble(n1) 2098 igrid = igrid + 1 2099 r_vec_grid(1,igrid)=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3) 2100 r_vec_grid(2,igrid)=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3) 2101 r_vec_grid(3,igrid)=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3) 2102 end do 2103 end do 2104 end do 2105 if (igrid/=ngrid) then 2106 ABI_ERROR('read_atomden: igrid not equal to ngrid') 2107 end if 2108 2109 !Read in atomic density data for each atom type 2110 !first check how many datapoints are in each file 2111 do itypat=1,ntypat 2112 filename='';io_err=0; 2113 if (itypat>0) write(filename,'(a,a,i1,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat' 2114 if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat' 2115 if (open_file(filename, message, newunit=unt, status='old',action='read') /= 0) then 2116 write(std_out,*) 'ERROR in read_atomden: Could not open file: ',filename 2117 write(std_out,*) ' Current implementation requires this file to be present' 2118 write(std_out,*) ' for each type of atom.' 2119 write(std_out,*)trim(message) 2120 ABI_ERROR("Cannot continue") 2121 end if 2122 ! Check number of lines in file 2123 nlines = 1;io_err=0; 2124 do 2125 read(unt,*,iostat=io_err) 2126 if (io_err<0) exit 2127 nlines = nlines + 1 2128 end do 2129 close(unt) 2130 natomgr(itypat) = nlines - 2 2131 end do ! Atom type 2132 !Allocate arrays and read in data 2133 natomgrmax = maxval(natomgr) 2134 ABI_MALLOC(atomrgrid,(natomgrmax,ntypat)) 2135 ABI_MALLOC(density,(natomgrmax,ntypat)) 2136 atomrgrid = zero ; density = zero 2137 do itypat=1,ntypat 2138 filename='';io_err=0; 2139 if (itypat>0) write(filename,'(a,a,i1,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat' 2140 if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat' 2141 if (open_file(filename,message,newunit=unt,status='old',action='read') /= 0) then 2142 ABI_ERROR(message) 2143 end if 2144 read(unt,*) ! Skip comment line 2145 do i=1,natomgr(itypat) 2146 read(unt,*) atomrgrid(i,itypat),density(i,itypat) 2147 end do 2148 close(unt) 2149 if (atomrgrid(1,itypat)/=zero) then 2150 write(std_out,*) 'ERROR in read_atomden, in file: ',filename 2151 write(std_out,*) ' First gridpoint has to be the origin.' 2152 ABI_ERROR("Cannot continue") 2153 end if 2154 end do ! Atom type 2155 2156 !write(std_out,*) '*** --- In read_atomden before call--- ***' 2157 !write(std_out,*) ' calctype:',calctype,' natom:',natom 2158 !write(std_out,*) ' ntypat:',ntypat,' typat:',typat 2159 !write(std_out,*) ' ngrid:',ngrid 2160 !write(std_out,*) ' a:',a 2161 !write(std_out,*) ' b:',b 2162 !write(std_out,*) ' c:',c 2163 !write(std_out,*) ' xred:',xred 2164 !write(std_out,*) ' natomgr:',natomgr 2165 !write(std_out,*) 'natomgrmax:',natomgrmax 2166 !write(std_out,*) ' atomrgrid:',atomrgrid 2167 !write(std_out,*) ' density:',density 2168 !write(std_out,*) 'r_vec_grid:' 2169 !write(std_out,*) r_vec_grid 2170 2171 !Call atomden 2172 call atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,xred, & 2173 & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype) 2174 2175 !if (prtvol>9) then ! calculate norm 2176 !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2177 !norm = SUM(rho(:))*ucvol/dble(n1*n2*n3) 2178 !write(message,'(a,F8.4)') ' In read_atomden - NORM OF DENSITY: ',norm 2179 !call wrtout(std_out,message,'COLL') 2180 !end if 2181 2182 rhor_atm(:,1) = rho 2183 2184 ABI_SFREE(atomrgrid) 2185 ABI_SFREE(density) 2186 ABI_SFREE(r_vec_grid) 2187 ABI_SFREE(rho) 2188 2189 end subroutine read_atomden