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