TABLE OF CONTENTS
ABINIT/calcdensph [ Functions ]
NAME
calcdensph
FUNCTION
Compute and print integral of total density inside spheres around atoms.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT,ILuk,MVer,EB,SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
gmet(3,3)=reciprocal space metric tensor in bohr**-2 mpi_enreg=information about MPI parallelization natom=number of atoms 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 nspden=number of spin-density components ntypat=number of atom types nunit=number of the unit for writing ratsph(ntypat)=radius of spheres around atoms rhor(nfft,nspden)=array for electron density in electrons/bohr**3. (total in first half and spin-up in second half if nspden=2) (total in first comp. and magnetization in comp. 2 to 4 if nspden=4) rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(natom)=type of each atom ucvol=unit cell volume in bohr**3 xred(3,natom)=reduced dimensionless atomic coordinates prtopt = if 1, the default printing is on, otherwise special printing options
OUTPUT
intgden(nspden, natom)=intgrated density (magnetization...) for each atom in a sphere of radius ratsph. Optional arg dentot(nspden)=integrated density (magnetization...) over full u.c. vol, optional argument Rest is printing
PARENTS
dfpt_scfcv,mag_constr,mag_constr_e,outscfcv
CHILDREN
timab,wrtout,xmpi_sum
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,nunit,ratsph,rhor,rprimd,typat,ucvol,xred,& 55 & prtopt,cplex,intgden,dentot) 56 57 use defs_basis 58 use defs_abitypes 59 use m_profiling_abi 60 use m_errors 61 62 use m_xmpi, only : xmpi_sum 63 64 !This section has been created automatically by the script Abilint (TD). 65 !Do not modify the following lines by hand. 66 #undef ABI_FUNC 67 #define ABI_FUNC 'calcdensph' 68 use interfaces_14_hidewrite 69 use interfaces_18_timing 70 use interfaces_32_util 71 !End of the abilint section 72 73 implicit none 74 75 !Arguments --------------------------------------------- 76 !scalars 77 integer,intent(in) :: natom,nfft,nspden,ntypat,nunit 78 real(dp),intent(in) :: ucvol 79 type(MPI_type),intent(in) :: mpi_enreg 80 integer ,intent(in) :: prtopt 81 integer, intent(in) :: cplex 82 !arrays 83 integer,intent(in) :: ngfft(18),typat(natom) 84 real(dp),intent(in) :: gmet(3,3),ratsph(ntypat),rhor(cplex*nfft,nspden),rprimd(3,3) 85 real(dp),intent(in) :: xred(3,natom) 86 !integer,intent(out),optional :: atgridpts(nfft) 87 real(dp),intent(out),optional :: intgden(nspden,natom) 88 real(dp),intent(out),optional :: dentot(nspden) 89 !Local variables ------------------------------ 90 !scalars 91 integer,parameter :: ishift=5 92 integer :: i1,i2,i3,iatom,ierr,ifft_local,ix,iy,iz,izloc,n1,n1a,n1b,n2,ifft 93 integer :: n2a,n2b,n3,n3a,n3b,nd3,nfftot 94 integer :: ii 95 !integer :: is,npts(natom) 96 integer :: cmplex_den,jfft 97 real(dp),parameter :: delta=0.99_dp 98 real(dp) :: difx,dify,difz,r2,r2atsph,rr1,rr2,rr3,rx,ry,rz 99 real(dp) :: fsm, ratsm, ratsm2 100 real(dp) :: mag_coll , mag_x, mag_y, mag_z ! EB 101 real(dp) :: mag_coll_im, mag_x_im, mag_y_im, mag_z_im ! SPr 102 real(dp) :: sum_mag, sum_mag_x,sum_mag_y,sum_mag_z,sum_rho_up,sum_rho_dn,sum_rho_tot ! EB 103 real(dp) :: rho_tot, rho_tot_im 104 ! real(dp) :: rho_up,rho_dn,rho_tot !EB 105 logical :: grid_found 106 character(len=500) :: message 107 !arrays 108 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 109 real(dp) :: intgden_(nspden,natom),tsec(2) 110 111 ! ************************************************************************* 112 113 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nd3=n3/mpi_enreg%nproc_fft 114 nfftot=n1*n2*n3 115 intgden_=zero 116 117 ratsm = zero 118 if(present(intgden)) then 119 120 ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later 121 122 end if 123 124 !Get the distrib associated with this fft_grid 125 grid_found=.false. 126 127 if(n2 == mpi_enreg%distribfft%n2_coarse ) then 128 129 if(n3== size(mpi_enreg%distribfft%tab_fftdp3_distrib)) then 130 131 fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3_distrib 132 ffti3_local => mpi_enreg%distribfft%tab_fftdp3_local 133 grid_found=.true. 134 135 end if 136 137 end if 138 139 if(n2 == mpi_enreg%distribfft%n2_fine ) then 140 141 if(n3 == size(mpi_enreg%distribfft%tab_fftdp3dg_distrib)) then 142 143 fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3dg_distrib 144 ffti3_local => mpi_enreg%distribfft%tab_fftdp3dg_local 145 grid_found = .true. 146 147 end if 148 149 end if 150 if(.not.(grid_found)) then 151 152 MSG_BUG("Unable to find an allocated distrib for this fft grid") 153 154 end if 155 156 !Loop over atoms 157 !------------------------------------------- 158 ii=0 159 160 do iatom=1,natom 161 162 !if (present(atgridpts)) then 163 ! npts(iatom)=0 !SPr: initialize the number of grid points within atomic sphere around atom i 164 ! ii=ii+1 !SPr: initialize running index for constructing an array of grid point indexes 165 !end if ! within atomic spheres 166 167 ! Define a "box" around the atom 168 r2atsph=1.0000001_dp*ratsph(typat(iatom))**2 169 rr1=sqrt(r2atsph*gmet(1,1)) 170 rr2=sqrt(r2atsph*gmet(2,2)) 171 rr3=sqrt(r2atsph*gmet(3,3)) 172 173 n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1 174 n1b=int((xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1 175 n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2 176 n2b=int((xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2 177 n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3 178 n3b=int((xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3 179 180 ratsm2 = (2*ratsph(typat(iatom))-ratsm)*ratsm 181 182 do i3=n3a,n3b 183 iz=mod(i3+ishift*n3,n3) 184 185 if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then 186 187 izloc = ffti3_local(iz+1) - 1 188 difz=dble(i3)/dble(n3)-xred(3,iatom) 189 do i2=n2a,n2b 190 iy=mod(i2+ishift*n2,n2) 191 dify=dble(i2)/dble(n2)-xred(2,iatom) 192 do i1=n1a,n1b 193 ix=mod(i1+ishift*n1,n1) 194 195 difx=dble(i1)/dble(n1)-xred(1,iatom) 196 rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3) 197 ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3) 198 rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3) 199 r2=rx**2+ry**2+rz**2 200 201 202 ! Identify the fft indexes of the rectangular grid around the atom 203 if(r2 > r2atsph) then 204 cycle 205 end if 206 207 fsm = radsmear(r2, r2atsph, ratsm2) 208 209 ifft_local=1+ix+n1*(iy+n2*izloc) 210 211 !if(present(atgridpts)) then 212 ! ii=ii+1 213 ! atgridpts(ii)=ifft_local !SPr: save grid point index (dbg: to check whether ifft_local is a valid "global" index ) 214 !end if 215 216 if(nspden==1) then 217 ! intgden_(1,iatom)= integral of total density 218 intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,1) 219 elseif(nspden==2) then 220 ! intgden_(1,iatom)= integral of up density 221 ! intgden_(2,iatom)= integral of dn density 222 intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,2) 223 intgden_(2,iatom)=intgden_(2,iatom)+fsm*rhor(ifft_local,1)-rhor(ifft_local,2) 224 else 225 ! intgden_(1,iatom)= integral of total density 226 ! intgden_(2,iatom)= integral of magnetization, x-component 227 ! intgden_(3,iatom)= integral of magnetization, y-component 228 ! intgden_(4,iatom)= integral of magnetization, z-component 229 intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,1) 230 intgden_(2,iatom)=intgden_(2,iatom)+fsm*rhor(ifft_local,2) 231 intgden_(3,iatom)=intgden_(3,iatom)+fsm*rhor(ifft_local,3) 232 intgden_(4,iatom)=intgden_(4,iatom)+fsm*rhor(ifft_local,4) 233 end if 234 235 end do 236 end do 237 end if 238 end do 239 240 intgden_(:,iatom)=intgden_(:,iatom)*ucvol/dble(nfftot) 241 242 !if (present(atgridpts)) then 243 ! npts(iatom)=ii-1 244 ! do is=1,iatom-1,1 245 ! npts(iatom)=npts(iatom)-npts(is)-1 246 ! end do 247 ! atgridpts(ii-npts(iatom))=npts(iatom) !SPr: save number of grid points around atom i 248 !end if 249 250 251 ! End loop over atoms 252 ! ------------------------------------------- 253 254 end do 255 256 ! EB - Compute magnetization of the whole cell 257 ! SPr - in case of complex density array set cmplex_den to one 258 cmplex_den=0 259 260 if(cplex==2) then 261 262 cmplex_den=1 263 264 end if 265 266 if(nspden==2) then 267 mag_coll=zero 268 mag_coll_im=zero 269 rho_tot=zero 270 rho_tot_im=zero 271 do ifft=1,nfft 272 jfft=(cmplex_den+1)*ifft 273 ! rho_up=rho_up+rhor(ifft,2) 274 ! rho_dn=rho_dn+(rhor(ifft,2)-rhor(ifft,1)) 275 rho_tot=rho_tot+rhor(jfft-cmplex_den,1) ! real parts of density and magnetization 276 mag_coll=mag_coll+2*rhor(jfft-cmplex_den,2)-rhor(jfft-cmplex_den,1) 277 278 rho_tot_im=rho_tot_im+rhor(jfft,1) ! imaginary parts of density and magnetization 279 mag_coll_im=mag_coll_im+2*rhor(jfft,2)-rhor(jfft,1) ! in case of real array both are equal to corresponding real parts 280 end do 281 mag_coll=mag_coll*ucvol/dble(nfftot) 282 rho_tot =rho_tot *ucvol/dble(nfftot) 283 284 mag_coll_im=mag_coll_im*ucvol/dble(nfftot) 285 rho_tot_im =rho_tot_im *ucvol/dble(nfftot) 286 ! rho_up=rho_up*ucvol/dble(nfftot) 287 ! rho_dn=rho_dn*ucvol/dble(nfftot) 288 ! rho_tot=rho_tot*ucvol/dble(nfftot) 289 else if(nspden==4) then 290 rho_tot=0 291 rho_tot_im=0 292 mag_x=0 293 mag_y=0 294 mag_z=0 295 mag_x_im=0 296 mag_y_im=0 297 mag_z_im=0 298 do ifft=1,nfft 299 jfft=(cmplex_den+1)*ifft 300 rho_tot=rho_tot+rhor(jfft-cmplex_den,1) 301 mag_x=mag_x+rhor(jfft-cmplex_den,2) 302 mag_y=mag_y+rhor(jfft-cmplex_den,3) 303 mag_z=mag_z+rhor(jfft-cmplex_den,4) 304 rho_tot_im=rho_tot_im+rhor(jfft,1) 305 mag_x_im=mag_x_im+rhor(jfft,2) 306 mag_y_im=mag_y_im+rhor(jfft,3) 307 mag_z_im=mag_z_im+rhor(jfft,4) 308 end do 309 rho_tot=rho_tot*ucvol/dble(nfftot) 310 mag_x=mag_x*ucvol/dble(nfftot) 311 mag_y=mag_y*ucvol/dble(nfftot) 312 mag_z=mag_z*ucvol/dble(nfftot) 313 rho_tot_im=rho_tot_im*ucvol/dble(nfftot) 314 mag_x_im=mag_x_im*ucvol/dble(nfftot) 315 mag_y_im=mag_y_im*ucvol/dble(nfftot) 316 mag_z_im=mag_z_im*ucvol/dble(nfftot) 317 end if 318 319 !MPI parallelization 320 if(mpi_enreg%nproc_fft>1)then 321 call timab(48,1,tsec) 322 call xmpi_sum(intgden_,mpi_enreg%comm_fft,ierr) 323 call xmpi_sum(rho_tot,mpi_enreg%comm_fft,ierr) ! EB 324 call xmpi_sum(mag_coll,mpi_enreg%comm_fft,ierr) ! EB 325 326 call xmpi_sum(rho_tot_im,mpi_enreg%comm_fft,ierr) 327 call xmpi_sum(mag_coll_im,mpi_enreg%comm_fft,ierr) 328 ! call xmpi_sum(rho_up,mpi_enreg%comm_fft,ierr) ! EB 329 ! call xmpi_sum(rho_dn,mpi_enreg%comm_fft,ierr) ! EB 330 call xmpi_sum(mag_x,mpi_enreg%comm_fft,ierr) ! EB 331 call xmpi_sum(mag_y,mpi_enreg%comm_fft,ierr) ! EB 332 call xmpi_sum(mag_z,mpi_enreg%comm_fft,ierr) ! EB 333 call xmpi_sum(mag_x_im,mpi_enreg%comm_fft,ierr) ! EB 334 call xmpi_sum(mag_y_im,mpi_enreg%comm_fft,ierr) ! EB 335 call xmpi_sum(mag_z_im,mpi_enreg%comm_fft,ierr) ! EB 336 call timab(48,2,tsec) 337 end if 338 339 !Printing 340 sum_mag=zero 341 sum_mag_x=zero 342 sum_mag_y=zero 343 sum_mag_z=zero 344 sum_rho_up=zero 345 sum_rho_dn=zero 346 sum_rho_tot=zero 347 348 if(prtopt==1) then 349 350 if(nspden==1) then 351 write(message, '(4a)' ) & 352 & ' Integrated electronic density in atomic spheres:',ch10,& 353 & ' ------------------------------------------------' 354 call wrtout(nunit,message,'COLL') 355 write(message, '(a)' ) ' Atom Sphere_radius Integrated_density' 356 call wrtout(nunit,message,'COLL') 357 do iatom=1,natom 358 write(message, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom) 359 call wrtout(nunit,message,'COLL') 360 end do 361 elseif(nspden==2) then 362 write(message, '(4a)' ) & 363 & ' Integrated electronic and magnetization densities in atomic spheres:',ch10,& 364 & ' ---------------------------------------------------------------------' 365 call wrtout(nunit,message,'COLL') 366 write(message, '(3a)' ) ' Note: Diff(up-dn) is a rough ',& 367 & 'approximation of local magnetic moment' 368 call wrtout(nunit,message,'COLL') 369 write(message, '(a)' ) ' Atom Radius up_density dn_density Total(up+dn) Diff(up-dn)' 370 call wrtout(nunit,message,'COLL') 371 do iatom=1,natom 372 write(message, '(i5,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom),intgden_(2,iatom),& 373 & ' ',(intgden_(1,iatom)+intgden_(2,iatom)),' ',(intgden_(1,iatom)-intgden_(2,iatom)) 374 call wrtout(nunit,message,'COLL') 375 ! Compute the sum of the magnetization 376 sum_mag=sum_mag+intgden_(1,iatom)-intgden_(2,iatom) 377 sum_rho_up=sum_rho_up+intgden_(1,iatom) 378 sum_rho_dn=sum_rho_dn+intgden_(2,iatom) 379 sum_rho_tot=sum_rho_tot+intgden_(1,iatom)+intgden_(2,iatom) 380 end do 381 write(message, '(a)') ' ---------------------------------------------------------------------' 382 call wrtout(nunit,message,'COLL') 383 write(message, '(a,2f13.6,a,f12.6,a,f12.6)') ' Sum: ', sum_rho_up,sum_rho_dn,' ',sum_rho_tot,' ',sum_mag 384 call wrtout(nunit,message,'COLL') 385 write(message, '(a,f14.6)') ' Total magnetization (from the atomic spheres): ', sum_mag 386 call wrtout(nunit,message,'COLL') 387 write(message, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll 388 call wrtout(nunit,message,'COLL') 389 ! EB for testing purpose print rho_up, rho_dn and rho_tot 390 ! write(message, '(a,3f14.4,2i8)') ' rho_up, rho_dn, rho_tot, nfftot, nfft: ', rho_up,rho_dn,rho_tot,nfft,nfft 391 ! call wrtout(nunit,message,'COLL') 392 393 elseif(nspden==4) then 394 395 write(message, '(4a)' ) & 396 & ' Integrated electronic and magnetization densities in atomic spheres:',ch10,& 397 & ' ---------------------------------------------------------------------' 398 call wrtout(nunit,message,'COLL') 399 write(message, '(3a)' ) ' Note: this is a rough approximation of local magnetic moments' 400 call wrtout(nunit,message,'COLL') 401 write(message, '(a)' ) ' Atom Radius Total density mag(x) mag(y) mag(z) ' 402 call wrtout(nunit,message,'COLL') 403 do iatom=1,natom 404 write(message, '(i5,f10.5,f16.6,a,3f12.6)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom),' ',(intgden_(ix,iatom),ix=2,4) 405 call wrtout(nunit,message,'COLL') 406 ! Compute the sum of the magnetization in x, y and z directions 407 sum_mag_x=sum_mag_x+intgden_(2,iatom) 408 sum_mag_y=sum_mag_y+intgden_(3,iatom) 409 sum_mag_z=sum_mag_z+intgden_(4,iatom) 410 end do 411 write(message, '(a)') ' ---------------------------------------------------------------------' 412 call wrtout(nunit,message,'COLL') 413 ! write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization : ', sum_mag_x,sum_mag_y,sum_mag_z 414 write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (spheres) ', sum_mag_x,sum_mag_y,sum_mag_z 415 call wrtout(nunit,message,'COLL') 416 write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (exact) ', mag_x,mag_y,mag_z 417 call wrtout(nunit,message,'COLL') 418 ! SPr for dfpt debug 419 ! write(message, '(a,f12.6)') ' Total density (exact) ', rho_tot 420 ! call wrtout(nunit,message,'COLL') 421 end if 422 423 elseif(prtopt==-1) then 424 425 write(message, '(2a)') ch10,' ------------------------------------------------------------------------' 426 call wrtout(nunit,message,'COLL') 427 428 if(nspden==1) then 429 write(message, '(4a)' ) & 430 & ' Fermi level charge density n_f:',ch10,& 431 & ' ------------------------------------------------------------------------',ch10 432 else 433 write(message, '(4a)' ) & 434 & ' Fermi level charge density n_f and magnetization m_f:',ch10,& 435 & ' ------------------------------------------------------------------------',ch10 436 end if 437 call wrtout(nunit,message,'COLL') 438 439 if(cmplex_den==0) then 440 write(message, '(a,f13.8)') ' n_f = ',rho_tot 441 else 442 write(message, '(a,f13.8,a,f13.8)') ' Re[n_f]= ', rho_tot," Im[n_f]= ",rho_tot_im 443 end if 444 call wrtout(nunit,message,'COLL') 445 if(nspden==2) then 446 if(cmplex_den==0) then 447 write(message, '(a,f13.8)') ' m_f = ', mag_coll 448 else 449 write(message, '(a,f13.8,a,f13.8)') ' Re[m_f]= ', mag_coll," Im[m_f]= ",mag_coll_im 450 end if 451 call wrtout(nunit,message,'COLL') 452 elseif (nspden==4) then 453 write(message, '(a,f13.8)') ' mx_f = ',mag_x 454 call wrtout(nunit,message,'COLL') 455 write(message, '(a,f13.8)') ' my_f = ',mag_y 456 call wrtout(nunit,message,'COLL') 457 write(message, '(a,f13.8)') ' mz_f = ',mag_z 458 call wrtout(nunit,message,'COLL') 459 end if 460 461 write(message, '(3a)') ch10,' ------------------------------------------------------------------------',ch10 462 call wrtout(nunit,message,'COLL') 463 464 465 else 466 ! prtopt different from 1 and -1 (either 2,3 or 4) 467 468 if(abs(rho_tot)<1.0d-10) then 469 rho_tot=0 470 end if 471 472 write(message, '(2a)') ch10,' ------------------------------------------------------------------------' 473 call wrtout(nunit,message,'COLL') 474 475 if(nspden==1) then 476 write(message, '(4a)' ) & 477 & ' Integral of the first order density n^(1):',ch10,& 478 & ' ------------------------------------------------------------------------',ch10 479 else 480 write(message, '(4a)' ) & 481 & ' Integrals of the first order density n^(1) and magnetization m^(1):',ch10,& 482 & ' ------------------------------------------------------------------------',ch10 483 end if 484 call wrtout(nunit,message,'COLL') 485 486 if(cmplex_den==0) then 487 write(message, '(a,e16.8)') ' n^(1) = ', rho_tot 488 else 489 write(message, '(a,e16.8,a,e16.8)') ' Re[n^(1)] = ', rho_tot," Im[n^(1)] = ",rho_tot_im 490 end if 491 call wrtout(nunit,message,'COLL') 492 493 if(nspden==2) then 494 495 if(cmplex_den==0) then 496 write(message, '(a,e16.8)') ' m^(1) = ', mag_coll 497 else 498 write(message, '(a,e16.8,a,e16.8)') ' Re[m^(1)] = ', mag_coll," Im[m^(1)] = ",mag_coll_im 499 end if 500 call wrtout(nunit,message,'COLL') 501 502 elseif (nspden==4) then 503 if(cmplex_den==0) then 504 write(message, '(a,e16.8)') ' mx^(1) = ', mag_x 505 call wrtout(nunit,message,'COLL') 506 write(message, '(a,e16.8)') ' my^(1) = ', mag_y 507 call wrtout(nunit,message,'COLL') 508 write(message, '(a,e16.8)') ' mz^(1) = ', mag_z 509 call wrtout(nunit,message,'COLL') 510 else 511 write(message, '(a,e16.8,a,e16.8)') ' Re[mx^(1)]= ', mag_x, " Im[mx^(1)]= ", mag_x_im 512 call wrtout(nunit,message,'COLL') 513 write(message, '(a,e16.8,a,e16.8)') ' Re[my^(1)]= ', mag_y, " Im[my^(1)]= ", mag_y_im 514 call wrtout(nunit,message,'COLL') 515 write(message, '(a,e16.8,a,e16.8)') ' Re[mz^(1)]= ', mag_z, " Im[mz^(1)]= ", mag_z_im 516 call wrtout(nunit,message,'COLL') 517 end if 518 end if 519 520 write(message, '(3a)') ch10,' ------------------------------------------------------------------------',ch10 521 call wrtout(nunit,message,'COLL') 522 523 end if 524 525 if(present(intgden)) then 526 intgden = intgden_ 527 end if 528 529 if(present(dentot)) then 530 if(nspden==2) then 531 dentot(1)=rho_tot 532 dentot(2)=mag_coll 533 elseif(nspden==4) then 534 dentot(1)=rho_tot 535 dentot(2)=mag_x 536 dentot(3)=mag_y 537 dentot(4)=mag_z 538 end if 539 end if 540 541 end subroutine calcdensph