TABLE OF CONTENTS
- ABINIT/m_dens
- ABINIT/printmagvtk
- m_dens/add_atomic_fcts
- m_dens/calcdenmagsph
- m_dens/constrained_dft_free
- m_dens/constrained_dft_ini
- m_dens/constrained_dft_t
- m_dens/constrained_residual
- m_dens/dens_hirsh
- m_dens/mag_penalty
- m_dens/mag_penalty_e
- m_dens/prtdenmagsph
- m_dens/radsmear
ABINIT/m_dens [ Modules ]
NAME
m_dens
FUNCTION
Module containing the definition of the constrained_dft_t data type and methods used to handle it, and also includes the computation of integrated atomic charge and magnetization, as well as Hirshfeld charges.
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_dens 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_splines 30 31 use defs_abitypes, only : MPI_type 32 use m_fft, only : fourdp 33 use m_time, only : timab 34 use m_numeric_tools, only : wrap2_zero_one 35 use m_io_tools, only : open_file 36 use m_geometry, only : dist2, xcart2xred, metric 37 use m_mpinfo, only : ptabs_fourdp 38 39 implicit none 40 41 private 42 43 public :: dens_hirsh ! Compute the Hirshfeld charges 44 public :: add_atomic_fcts ! Add atomic functions to real space function 45 public :: constrained_dft_ini ! Initialize the constrained_dft datastructure 46 public :: constrained_dft_free ! Free the constrained_dft datastructure 47 public :: constrained_residual ! Recompute the potential residual, to account for constraints 48 public :: mag_penalty ! Compute the potential corresponding to constrained magnetic moments (using add_atomic_fcts) with the penalty function. 49 public :: mag_penalty_e ! Compute the energy corresponding to constrained magnetic moments. 50 public :: calcdenmagsph ! Compute integral of total density and magnetization inside spheres around atoms. 51 public :: prtdenmagsph ! Print integral of total density and magnetization inside spheres around atoms.
ABINIT/printmagvtk [ Functions ]
NAME
printmagvtk
FUNCTION
Auxiliary routine for printing out magnetization density in VTK format. Output file name is DEN.vtk
INPUTS
mpi_enreg = information about adopted parallelization strategy
nspden = number of components of density matrix (possible values re 1,2, or 4)
nspden: 1 -> rho
nspden: 2 -> rho_up,rho_dwn
nspden: 4 -> rho,mx,my,mz
nfft = number of fft points per FFT processor
ngfft = full information about FFT mesh
rhor = density array stored in the memory of current FFT processor
rprimd = array of lattice vectors
OUTPUT
SIDE EFFECTS
NOTES
At the moment this routine is mainly used for development and debugging
of gs and dfpt calculations with non-collinear spins. If needed, can be used
to print final density in vtk format.
IMPORTANT: implementation is thoroughly checked only for npspinor = 1,
for other case might need to change the part gathering
the FFT mesh info
SOURCE
2317 subroutine printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,fname) 2318 2319 !Arguments ------------------------------------ 2320 !scalars 2321 type(MPI_type),intent(in) :: mpi_enreg 2322 integer,intent(in) :: nfft,nspden,cplex 2323 !arrays 2324 integer,intent(in) :: ngfft(18) 2325 real(dp),intent(in) :: rhor(cplex*nfft,nspden),rprimd(3,3) 2326 character(len=*),intent(in) :: fname 2327 2328 !Local variables------------------------------- 2329 !scalars 2330 integer :: denvtk,denxyz,denxyz_im,nfields 2331 integer :: nx,ny,nz,nfft_tot 2332 integer :: ii,jj,kk,ind,ispden 2333 integer :: mpi_comm,mpi_head,mpi_rank,ierr 2334 real(dp) :: rx,ry,rz 2335 integer :: nproc_fft,ir 2336 character(len=500) :: msg 2337 character(len=10) :: outformat 2338 character(len=50) :: fname_vtk 2339 character(len=50) :: fname_xyz 2340 character(len=50) :: fname_xyz_re 2341 character(len=50) :: fname_xyz_im 2342 !arrays 2343 real(dp),allocatable :: rhorfull(:,:) 2344 2345 2346 ! ************************************************************************* 2347 2348 DBG_ENTER("COLL") 2349 2350 ! if (option/=1 .and. option/=2 ) then 2351 ! write(msg,'(3a,i0)')& 2352 !& 'The argument option should be 1 or 2,',ch10,& 2353 !& 'however, option=',option 2354 ! ABI_BUG(msg) 2355 ! end if 2356 ! 2357 ! if (sizein<1) then 2358 ! write(msg,'(3a,i0)')& 2359 !& 'The argument sizein should be a positive number,',ch10,& 2360 !& 'however, sizein=',sizein 2361 ! ABI_ERROR(msg) 2362 ! end if 2363 2364 DBG_EXIT("COLL") 2365 2366 fname_vtk=adjustl(adjustr(fname)//".vtk") 2367 fname_xyz=adjustl(adjustr(fname)//".xyz") 2368 fname_xyz_re=adjustl(adjustr(fname)//"_re.xyz") 2369 fname_xyz_im=adjustl(adjustr(fname)//"_im.xyz") 2370 !write(std_out,*) ' Writing out .vtk file: ',fname_vtk 2371 !write(std_out,*) ' Writing out .xyz file: ',fname_xyz 2372 2373 !if 1 or two component density then write out either 1 or 2 scalar density fields 2374 !if 4, then write one scalar field (density) and one vector field (magnetization density) 2375 if(nspden/=4)then 2376 nfields=nspden 2377 else 2378 nfields=2 2379 end if 2380 2381 nfields=nfields*cplex 2382 2383 ! FFT mesh specifications: full grid 2384 nx=ngfft(1) ! number of points along 1st lattice vector 2385 ny=ngfft(2) ! number of points along 2nd lattice vector 2386 nz=ngfft(3) ! number of points along 3rd lattice vector 2387 nfft_tot=nx*ny*nz ! total number of fft mesh points (can be different from nfft in case of distributed memory of nproc_fft processors) 2388 2389 2390 ! Gather information about memory distribution 2391 mpi_head=0 2392 mpi_comm = mpi_enreg%comm_fft 2393 mpi_rank = xmpi_comm_rank(mpi_comm) 2394 nproc_fft=ngfft(10) 2395 2396 ! Create array to host full FFT mesh 2397 if(mpi_rank==mpi_head)then 2398 ABI_MALLOC(rhorfull,(cplex*nfft_tot,nspden)) 2399 end if 2400 2401 ! Fill in the full mesh 2402 if(nproc_fft==1)then 2403 rhorfull=rhor 2404 else 2405 do ir=1,nspden 2406 call xmpi_gather(rhor(:,ir),cplex*nfft,rhorfull(:,ir),cplex*nfft,mpi_head,mpi_comm,ierr) 2407 end do 2408 end if 2409 2410 if(mpi_rank==mpi_head)then 2411 2412 ! Open the output vtk file 2413 if (open_file(fname_vtk,msg,newunit=denvtk,status='replace',form='formatted') /=0) then 2414 ABI_WARNING(msg) 2415 RETURN 2416 end if 2417 2418 if(cplex==1) then 2419 if (open_file(fname_xyz,msg,newunit=denxyz,status='replace',form='formatted') /=0) then 2420 ABI_WARNING(msg) 2421 RETURN 2422 end if 2423 else if (cplex==2) then 2424 if (open_file(fname_xyz_re,msg,newunit=denxyz,status='replace',form='formatted') /=0) then 2425 ABI_WARNING(msg) 2426 RETURN 2427 end if 2428 if (open_file(fname_xyz_im,msg,newunit=denxyz_im,status='replace',form='formatted') /=0) then 2429 ABI_WARNING(msg) 2430 RETURN 2431 end if 2432 end if 2433 2434 ! Write the header of the output vtk file 2435 write(denvtk,"(a)") '# vtk DataFile Version 2.0' 2436 write(denvtk,"(a)") 'Electron density components' 2437 write(denvtk,"(a)") 'ASCII' 2438 write(denvtk,"(a)") 'DATASET STRUCTURED_GRID' 2439 write(denvtk,"(a,3i6)") 'DIMENSIONS ', nx,ny,nz 2440 write(denvtk,"(a,i18,a)") 'POINTS ',nfft_tot,' double' 2441 2442 if (nspden==1) then 2443 outformat="(4e16.8)" 2444 else if (nspden==2) then 2445 outformat="(5e16.8)" 2446 else 2447 outformat="(7e16.8)" 2448 end if 2449 2450 ! Write out information about grid points 2451 do kk=0,nz-1 2452 do jj=0,ny-1 2453 do ii=0,nx-1 2454 2455 rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 2456 ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 2457 rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 2458 write(denvtk,'(3f16.8)') rx,ry,rz !coordinates of the grid point 2459 ind=1+ii+nx*(jj+ny*kk) 2460 if (cplex==1) then 2461 write(denxyz,outformat) rx,ry,rz,(rhorfull(ind,ispden),ispden=1,nspden) 2462 else 2463 write(denxyz,outformat) rx,ry,rz,(rhorfull(2*ind-1,ispden),ispden=1,nspden) 2464 write(denxyz_im,outformat) rx,ry,rz,(rhorfull(2*ind ,ispden),ispden=1,nspden) 2465 end if 2466 end do 2467 end do 2468 end do 2469 2470 if(cplex==1) then 2471 close(denxyz) 2472 else 2473 close(denxyz) 2474 close(denxyz_im) 2475 end if 2476 2477 ! Write out information about field defined on the FFT mesh 2478 write(denvtk,"(a,i18)") 'POINT_DATA ',nfft_tot 2479 write(denvtk,"(a,i6)") 'FIELD Densities ',nfields 2480 2481 2482 ! Write out different fields depending on the number of density matrix components 2483 if(nspden==1)then 2484 2485 !single component, so just write out the density 2486 if(cplex==1) then 2487 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 2488 do kk=0,nz-1 2489 do jj=0,ny-1 2490 do ii=0,nx-1 2491 ind=1+ii+nx*(jj+ny*kk) 2492 write(denvtk,'(f16.8)') rhorfull(ind,1) 2493 end do 2494 end do 2495 end do 2496 else 2497 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 2498 do kk=0,nz-1 2499 do jj=0,ny-1 2500 do ii=0,nx-1 2501 ind=2*(1+ii+nx*(jj+ny*kk))-1 2502 write(denvtk,'(f16.8)') rhorfull(ind,1) 2503 end do 2504 end do 2505 end do 2506 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 2507 do kk=0,nz-1 2508 do jj=0,ny-1 2509 do ii=0,nx-1 2510 ind=2*(1+ii+nx*(jj+ny*kk)) 2511 write(denvtk,'(f16.8)') rhorfull(ind,1) 2512 end do 2513 end do 2514 end do 2515 end if 2516 2517 else if(nspden==2)then 2518 2519 !two component, write the density for spin_up and spin_down channels 2520 if(cplex==1) then 2521 2522 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 2523 do kk=0,nz-1 2524 do jj=0,ny-1 2525 do ii=0,nx-1 2526 ind=1+ii+nx*(jj+ny*kk) 2527 write(denvtk,'(f16.8)') rhorfull(ind,1) 2528 end do 2529 end do 2530 end do 2531 write(denvtk,"(a,i18,a)") 'mag 1 ',nfft_tot,' double' 2532 do kk=0,nz-1 2533 do jj=0,ny-1 2534 do ii=0,nx-1 2535 ind=1+ii+nx*(jj+ny*kk) 2536 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 2537 end do 2538 end do 2539 end do 2540 2541 else 2542 2543 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 2544 do kk=0,nz-1 2545 do jj=0,ny-1 2546 do ii=0,nx-1 2547 ind=2*(1+ii+nx*(jj+ny*kk))-1 2548 write(denvtk,'(f16.8)') rhorfull(ind,1) 2549 end do 2550 end do 2551 end do 2552 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 2553 do kk=0,nz-1 2554 do jj=0,ny-1 2555 do ii=0,nx-1 2556 ind=2*(1+ii+nx*(jj+ny*kk)) 2557 write(denvtk,'(f16.8)') rhorfull(ind,1) 2558 end do 2559 end do 2560 end do 2561 write(denvtk,"(a,i18,a)") 'Re_mag 1 ',nfft_tot,' double' 2562 do kk=0,nz-1 2563 do jj=0,ny-1 2564 do ii=0,nx-1 2565 ind=2*(1+ii+nx*(jj+ny*kk))-1 2566 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 2567 end do 2568 end do 2569 end do 2570 write(denvtk,"(a,i18,a)") 'Im_mag 1 ',nfft_tot,' double' 2571 do kk=0,nz-1 2572 do jj=0,ny-1 2573 do ii=0,nx-1 2574 ind=2*(1+ii+nx*(jj+ny*kk)) 2575 write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1) 2576 end do 2577 end do 2578 end do 2579 2580 end if 2581 2582 else !here is the last option: nspden==4 2583 2584 if(cplex==1) then 2585 2586 !four component, write the density (scalar field) and magnetization density (vector field) 2587 write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double' 2588 do kk=0,nz-1 2589 do jj=0,ny-1 2590 do ii=0,nx-1 2591 ind=1+ii+nx*(jj+ny*kk) 2592 write(denvtk,'(f16.8)') rhorfull(ind,1) 2593 end do 2594 end do 2595 end do 2596 write(denvtk,"(a,i18,a)") 'mag 3 ',nfft_tot,' double' 2597 do kk=0,nz-1 2598 do jj=0,ny-1 2599 do ii=0,nx-1 2600 ind=1+ii+nx*(jj+ny*kk) 2601 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 2602 end do 2603 end do 2604 end do 2605 2606 else 2607 2608 write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double' 2609 do kk=0,nz-1 2610 do jj=0,ny-1 2611 do ii=0,nx-1 2612 ind=2*(1+ii+nx*(jj+ny*kk))-1 2613 write(denvtk,'(f16.8)') rhorfull(ind,1) 2614 end do 2615 end do 2616 end do 2617 write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double' 2618 do kk=0,nz-1 2619 do jj=0,ny-1 2620 do ii=0,nx-1 2621 ind=2*(1+ii+nx*(jj+ny*kk)) 2622 write(denvtk,'(f16.8)') rhorfull(ind,1) 2623 end do 2624 end do 2625 end do 2626 write(denvtk,"(a,i18,a)") 'Re_mag 3 ',nfft_tot,' double' 2627 do kk=0,nz-1 2628 do jj=0,ny-1 2629 do ii=0,nx-1 2630 ind=2*(1+ii+nx*(jj+ny*kk))-1 2631 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 2632 end do 2633 end do 2634 end do 2635 write(denvtk,"(a,i18,a)") 'Im_mag 3 ',nfft_tot,' double' 2636 do kk=0,nz-1 2637 do jj=0,ny-1 2638 do ii=0,nx-1 2639 ind=2*(1+ii+nx*(jj+ny*kk)) 2640 write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4) 2641 end do 2642 end do 2643 end do 2644 2645 end if 2646 2647 end if ! nspden options condition 2648 2649 close (denvtk) 2650 2651 !clean up the gathered FFT mesh 2652 ABI_FREE(rhorfull) 2653 2654 end if 2655 2656 end subroutine printmagvtk
m_dens/add_atomic_fcts [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
add_atomic_fcts
FUNCTION
This routine is called to assemble the atomic spherical functions, and, if option/=0, to add it to some input function (usually an input potential residual). The contributions from each atomic sphere are governed by parameters coeff_constr_dft, input to the present routine.
INPUTS
natom=number of atoms nspden = number of spin densities (1 2 or 4) option= if 0, the sum of the atomic spherical functions is returned in nv_constr_dft_r; if non-zero, they are added to nv_constr_dft_r rprimd=lattice vectors (dimensionful) mpi_enreg=mpi structure with communicator info nfft=number of points in standard fft grid ngfft=FFT grid dimensions ntypat=number of types of atoms ratsph(ntypat)=radii for muffin tin spheres of each atom typat(natom)=types of atoms xred(3,natom)=reduced atomic positions
SIDE EFFECTS
nv_constr_dft_r=the constrained potential or density in real space
SOURCE
454 subroutine add_atomic_fcts(natom,nspden,rprimd,mpi_enreg,nfft,ngfft,ntypat,option,ratsph, & 455 ratsm, typat,coeffs_constr_dft,nv_constr_dft_r,xred) 456 457 !Arguments ------------------------------------ 458 !scalars 459 integer,intent(in) :: natom,nfft,nspden,ntypat,option 460 type(MPI_type),intent(in) :: mpi_enreg 461 !arrays 462 integer,intent(in) :: typat(natom) 463 integer,intent(in) :: ngfft(18) 464 real(dp),intent(in) :: coeffs_constr_dft(nspden,natom) 465 real(dp),intent(inout) :: nv_constr_dft_r(nfft,nspden) 466 real(dp),intent(in) :: ratsph(ntypat) 467 real(dp),intent(in) :: ratsm ! Ben change 468 real(dp),intent(in) :: rprimd(3,3) 469 real(dp),intent(in) :: xred(3,natom) 470 471 !Local variables------------------------------- 472 !scalars 473 integer,parameter :: ishift=5 474 integer :: iatom, ierr 475 integer :: n1a, n1b, n3a, n3b, n2a, n2b 476 integer :: n1, n2, n3 477 integer :: ifft_local 478 integer :: i1,i2,i3,ix,iy,iz,izloc 479 real(dp) :: dfsm,dify,difz,fsm,r2atsph,rr1,rr2,rr3,ratsm2,rx23,ry23,rz23 !Ben change: remove ratsm 480 real(dp) :: r2,r2_11,r2_123,r2_23 481 real(dp) :: ucvol 482 real(dp),parameter :: delta=0.99_dp 483 !arrays 484 real(dp), allocatable :: difx(:) 485 real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3) 486 real(dp) :: tsec(2) 487 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 488 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 489 490 ! *********************************************************************************************** 491 492 !We need the metric because it is needed to compute the "box" around each atom 493 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 494 495 n1 = ngfft(1) 496 n2 = ngfft(2) 497 n3 = ngfft(3) 498 499 !Ben change: comment out ratsm line 500 !ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later 501 502 if(option==0)then 503 nv_constr_dft_r = zero 504 endif 505 506 !Get the distrib associated with this fft_grid 507 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 508 509 !Loop over atoms 510 !------------------------------------------- 511 do iatom=1,natom 512 513 if(sum(coeffs_constr_dft(1:nspden,iatom)**2)<tol12)then 514 cycle 515 endif 516 517 ! Define a "box" around the atom 518 r2atsph=1.0000001_dp*ratsph(typat(iatom))**2 519 rr1=sqrt(r2atsph*gmet(1,1)) 520 rr2=sqrt(r2atsph*gmet(2,2)) 521 rr3=sqrt(r2atsph*gmet(3,3)) 522 523 n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1 524 n1b=int((xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1 525 n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2 526 n2b=int((xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2 527 n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3 528 n3b=int((xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3 529 530 ratsm2 = -(ratsm**2 - 2*ratsph(typat(iatom))*ratsm) 531 532 ABI_MALLOC(difx,(n1a:n1b)) 533 do i1=n1a,n1b 534 difx(i1)=dble(i1)/dble(n1)-xred(1,iatom) 535 enddo ! i1 536 537 do i3=n3a,n3b 538 iz=mod(i3+ishift*n3,n3) 539 if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then 540 izloc = ffti3_local(iz+1) - 1 541 difz=dble(i3)/dble(n3)-xred(3,iatom) 542 do i2=n2a,n2b 543 iy=mod(i2+ishift*n2,n2) 544 dify=dble(i2)/dble(n2)-xred(2,iatom) 545 rx23=dify*rprimd(1,2)+difz*rprimd(1,3) 546 ry23=dify*rprimd(2,2)+difz*rprimd(2,3) 547 rz23=dify*rprimd(3,2)+difz*rprimd(3,3) 548 r2_23=rx23**2+ry23**2+rz23**2 549 r2_11=rprimd(1,1)**2+rprimd(2,1)**2+rprimd(3,1)**2 550 r2_123=2*(rprimd(1,1)*rx23+rprimd(2,1)*ry23+rprimd(3,1)*rz23) 551 do i1=n1a,n1b 552 r2=(difx(i1)*r2_11+r2_123)*difx(i1)+r2_23 553 if (r2 > r2atsph) cycle 554 call radsmear(dfsm,fsm,r2,r2atsph,ratsm2) 555 ix=mod(i1+ishift*n1,n1) 556 ! Identify the fft indexes of the rectangular grid around the atom 557 ifft_local=1+ix+n1*(iy+n2*izloc) 558 nv_constr_dft_r(ifft_local,1:nspden)=nv_constr_dft_r(ifft_local,1:nspden) + fsm*coeffs_constr_dft(1:nspden,iatom) 559 560 end do ! i1 561 end do ! i2 562 end if ! if this is my fft slice 563 end do ! i3 564 ABI_FREE(difx) 565 566 ! end loop over atoms 567 end do 568 569 !MPI parallelization 570 !TODO: test if xmpi_sum does the correct stuff for a slice of nv_constr_dft_r 571 if(mpi_enreg%nproc_fft>1)then 572 call timab(48,1,tsec) 573 call xmpi_sum(nv_constr_dft_r,mpi_enreg%comm_fft,ierr) 574 call timab(48,2,tsec) 575 end if 576 577 ! write (201,*) '# potential 1' 578 ! write (201,*) nv_constr_dft_r(:,1) 579 580 ! write (202,*) '# potential 2' 581 ! write (202,*) nv_constr_dft_r(:,2) 582 583 ! if (nspden > 2) then 584 ! write (203,*) '# potential 3' 585 ! write (203,*) nv_constr_dft_r(:,3) 586 587 ! write (204,*) '# potential 4' 588 ! write (204,*) nv_constr_dft_r(:,4) 589 ! end if 590 591 end subroutine add_atomic_fcts
m_dens/calcdenmagsph [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
calcdenmagsph
FUNCTION
Compute and print integral of total density inside spheres around atoms, or optionally of integral of potential residual. Also can compute the contributions to forces and stresses due to density-magnetization type constraints.
INPUTS
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 option = if not larger than 10, then a density is input , if larger than 10 then a potential residual is input. ratsm=smearing width for ratsph 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 xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
dentot(nspden)=integrated density (magnetization...) over full u.c. vol. Optional argument gr_intgden(3,nspden,natom)=grad wrt atomic positions, of integrated density (magnetization...) for each atom in a sphere. Optional arg intgden(nspden, natom)=integrated density (magnetization...) for each atom in a sphere of radius ratsph. Optional arg Note that when intgden is present, the definition of the spherical integration function changes, as it is smoothed. intgf2(natom,natom)=overlaps of the spherical integration functions for each atom in a sphere of radius ratsph. Optional arg rhomag(2,nspden)=integrated complex density (magnetization...) over full u.c. vol. Optional argument In collinear case component 1 is total density and 2 is _magnetization_ up-down In non collinear case component 1 is total density, and 2:4 are the magnetization vector strs_intgden(6,nspden,natom)=stress contribution due to constrained integrated density (magnetization...), due to each atom. Optional arg Rest is printing
SOURCE
1517 subroutine calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,typat,xred,& 1518 & option,cplex,dentot,gr_intgden,intgden,intgf2,rhomag,strs_intgden) 1519 1520 !Arguments --------------------------------------------- 1521 !scalars 1522 integer,intent(in) :: natom,nfft,nspden,ntypat 1523 real(dp),intent(in) :: ratsm 1524 type(MPI_type),intent(in) :: mpi_enreg 1525 integer ,intent(in) :: option 1526 integer, intent(in) :: cplex 1527 !arrays 1528 integer,intent(in) :: ngfft(18),typat(natom) 1529 real(dp),intent(in) :: ratsph(ntypat),rhor(cplex*nfft,nspden),rprimd(3,3) 1530 real(dp),intent(in) :: xred(3,natom) 1531 real(dp),intent(out),optional :: dentot(nspden) 1532 real(dp),intent(out),optional :: gr_intgden(3,nspden,natom) 1533 real(dp),intent(out),optional :: intgden(nspden,natom) 1534 real(dp),intent(out),optional :: intgf2(natom,natom) 1535 real(dp),intent(out),optional :: rhomag(2,nspden) 1536 real(dp),intent(out),optional :: strs_intgden(6,nspden,natom) 1537 !Local variables ------------------------------ 1538 !scalars 1539 integer,parameter :: ndir=3,ishift=5 1540 integer :: i1,i2,i3,iatom,ierr,ifft_local,ii,isp,ispden,ix,iy,iz,izloc,jatom,n1,n1a,n1b,n2,ifft 1541 integer :: neighbor_overlap,n2a,n2b,n3,n3a,n3b,nfftot 1542 integer :: jfft 1543 real(dp),parameter :: delta=0.99_dp 1544 real(dp) :: difx,dify,difz,r2,r2atsph,rr1,rr2,rr3,rx,ry,rz 1545 real(dp) :: dfsm,fact,fsm,ratsm2,ucvol 1546 logical :: grid_found 1547 !arrays 1548 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1549 integer :: overlap_ij(natom,natom) 1550 real(dp) :: gmet(3,3),gprimd(3,3),gr_intg(3,4) 1551 real(dp) :: intg(4),rhomag_(2,nspden) 1552 real(dp) :: strs(3,3),strs_cartred(3,3),strs_intg(6,4),tsec(2) 1553 real(dp) :: dist_ij(natom,natom),intgden_(nspden,natom) 1554 real(dp) :: my_xred(3, natom), rmet(3,3),xshift(3, natom) 1555 real(dp), allocatable :: fsm_atom(:,:) 1556 1557 !real(dp) :: rprimd_mod(3,3),strain 1558 1559 ! ************************************************************************* 1560 1561 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 1562 nfftot=n1*n2*n3 1563 intgden_=zero 1564 if(present(intgden))then 1565 intgden=zero 1566 endif 1567 if(present(gr_intgden))then 1568 gr_intgden=zero 1569 endif 1570 if(present(strs_intgden))then 1571 strs_intgden=zero 1572 endif 1573 1574 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1575 1576 ! This routine is not able to handle xred positions that are "far" from the 1577 ! first unit cell so wrap xred into [0, 1[ interval here. 1578 call wrap2_zero_one(xred, my_xred, xshift) 1579 1580 !If intgf2 present, check that the spheres do not overlap. If they overlap, one needs to treat explicitly the neighbors. 1581 if(present(intgf2))then 1582 neighbor_overlap=0 1583 dist_ij(:,:)=zero 1584 intgf2(:,:)=zero 1585 overlap_ij(:,:)=0 1586 dist_ij(1,1)=dist2(xred(:,1),xred(:,1),rprimd,-1) 1587 do iatom=1,natom 1588 dist_ij(iatom,iatom)=dist_ij(1,1) 1589 overlap_ij(iatom,iatom)=1 1590 do jatom=iatom,natom 1591 if(iatom/=jatom)dist_ij(iatom,jatom)=dist2(xred(:,iatom),xred(:,jatom),rprimd,1) 1592 if(dist_ij(iatom,jatom)-ratsph(typat(iatom))-ratsph(typat(jatom))<tol10)then 1593 overlap_ij(iatom,jatom)=1 1594 neighbor_overlap=1 1595 endif 1596 enddo 1597 enddo 1598 if(neighbor_overlap==1)then 1599 ABI_MALLOC(fsm_atom,(nfft,natom)) 1600 fsm_atom(:,:)=zero 1601 endif 1602 endif 1603 1604 !Get the distrib associated with this fft_grid 1605 grid_found=.false. 1606 1607 if(n2 == mpi_enreg%distribfft%n2_coarse ) then 1608 if(n3== size(mpi_enreg%distribfft%tab_fftdp3_distrib)) then 1609 fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3_distrib 1610 ffti3_local => mpi_enreg%distribfft%tab_fftdp3_local 1611 grid_found=.true. 1612 end if 1613 end if 1614 1615 if(n2 == mpi_enreg%distribfft%n2_fine ) then 1616 if(n3 == size(mpi_enreg%distribfft%tab_fftdp3dg_distrib)) then 1617 fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3dg_distrib 1618 ffti3_local => mpi_enreg%distribfft%tab_fftdp3dg_local 1619 grid_found = .true. 1620 end if 1621 end if 1622 1623 if(.not.(grid_found)) then 1624 ABI_BUG("Unable to find an allocated distrib for this fft grid") 1625 end if 1626 1627 !Loop over atoms 1628 !------------------------------------------- 1629 do iatom=1,natom 1630 1631 ! Define a "box" around the atom 1632 r2atsph=1.0000001_dp*ratsph(typat(iatom))**2 1633 rr1=sqrt(r2atsph*gmet(1,1)) 1634 rr2=sqrt(r2atsph*gmet(2,2)) 1635 rr3=sqrt(r2atsph*gmet(3,3)) 1636 1637 n1a=int((my_xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1 1638 n1b=int((my_xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1 1639 n2a=int((my_xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2 1640 n2b=int((my_xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2 1641 n3a=int((my_xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3 1642 n3b=int((my_xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3 1643 1644 !This is the "width" of the zone of smearing, in term of the square of radius 1645 ratsm2 = (2*ratsph(typat(iatom))-ratsm)*ratsm 1646 1647 intg(:)=zero 1648 gr_intg(:,:)=zero 1649 strs_intg(:,:)=zero 1650 1651 do i3=n3a,n3b 1652 iz=mod(i3+ishift*n3,n3) 1653 1654 if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then 1655 1656 izloc = ffti3_local(iz+1) - 1 1657 difz=dble(i3)/dble(n3)-my_xred(3,iatom) 1658 do i2=n2a,n2b 1659 iy=mod(i2+ishift*n2,n2) 1660 dify=dble(i2)/dble(n2)-my_xred(2,iatom) 1661 do i1=n1a,n1b 1662 ix=mod(i1+ishift*n1,n1) 1663 1664 difx=dble(i1)/dble(n1)-my_xred(1,iatom) 1665 !DEBUG 1666 ! if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then 1667 ! difx=dble(i1)/dble(n1)-(my_xred(1,iatom)+0.00005) 1668 ! endif 1669 !ENDDEBUG 1670 rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3) 1671 1672 !DEBUG 1673 ! if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then 1674 ! strain=-0.001 1675 ! rx=difx*rprimd(1,1)*(one+strain)+dify*rprimd(1,2)+difz*rprimd(1,3) 1676 ! endif 1677 !ENDDEBUG 1678 1679 ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3) 1680 rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3) 1681 r2=rx**2+ry**2+rz**2 1682 1683 1684 ! Identify the fft indexes of the rectangular grid around the atom 1685 if(r2 > r2atsph) then 1686 cycle 1687 end if 1688 1689 call radsmear(dfsm,fsm,r2,r2atsph,ratsm2) 1690 1691 ifft_local=1+ix+n1*(iy+n2*izloc) 1692 1693 if(present(intgf2))then 1694 ! intgden_(1,iatom)= integral of the square of the spherical integrating function 1695 if(neighbor_overlap==0)intgf2(iatom,iatom)=intgf2(iatom,iatom)+fsm*fsm 1696 if(neighbor_overlap==1)fsm_atom(ifft_local,iatom)=fsm_atom(ifft_local,iatom)+fsm 1697 endif 1698 ! Integral of density or potential residual 1699 intg(1:nspden)=intg(1:nspden)+fsm*rhor(ifft_local,1:nspden) 1700 if((present(gr_intgden).or.present(strs_intgden)).and. option<10 .and. ratsm2>tol12)then 1701 do ispden=1,nspden 1702 fact=dfsm*rhor(ifft_local,ispden) 1703 if(present(gr_intgden))then 1704 gr_intg(1,ispden)=gr_intg(1,ispden)+difx*fact 1705 gr_intg(2,ispden)=gr_intg(2,ispden)+dify*fact 1706 gr_intg(3,ispden)=gr_intg(3,ispden)+difz*fact 1707 endif 1708 if(present(strs_intgden))then 1709 strs_intg(1,ispden)=strs_intg(1,ispden)+difx*difx*fact 1710 strs_intg(2,ispden)=strs_intg(2,ispden)+dify*dify*fact 1711 strs_intg(3,ispden)=strs_intg(3,ispden)+difz*difz*fact 1712 strs_intg(4,ispden)=strs_intg(4,ispden)+dify*difz*fact 1713 strs_intg(5,ispden)=strs_intg(5,ispden)+difx*difz*fact 1714 strs_intg(6,ispden)=strs_intg(6,ispden)+difx*dify*fact 1715 endif 1716 enddo 1717 endif 1718 end do 1719 end do 1720 end if 1721 end do 1722 1723 if(present(intgf2) .and. neighbor_overlap==0)then 1724 intgf2(iatom,iatom)=intgf2(iatom,iatom)*ucvol/dble(nfftot) 1725 endif 1726 1727 intg(:)=intg(:)*ucvol/dble(nfftot) 1728 1729 if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then 1730 ! Convert to gradient in reduced coordinates 1731 gr_intg=matmul(rmet,gr_intg) 1732 gr_intg(:,:)=-gr_intg(:,:)*two*ucvol/dble(nfftot) 1733 !DEBUG 1734 ! write(6,*)' calcdenmagsph : intg(1)=',intg(1) 1735 ! write(6,*)' calcdenmagsph : iatom,gr_intg(1,1)=',iatom,gr_intg(1,1) 1736 ! call flush(6) 1737 ! stop 1738 !ENDDEBUG 1739 endif 1740 1741 if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then 1742 ! Convert to stress in cartesian coordinates, for each spin constraint 1743 do isp=1,4 1744 ! First change the representation 1745 strs(1,1)=strs_intg(1,isp) ; strs(2,2)=strs_intg(2,isp) ; strs(3,3)=strs_intg(3,isp) 1746 strs(2,3)=strs_intg(4,isp) ; strs(3,2)=strs_intg(4,isp) 1747 strs(1,3)=strs_intg(5,isp) ; strs(3,1)=strs_intg(5,isp) 1748 strs(1,2)=strs_intg(6,isp) ; strs(2,1)=strs_intg(6,isp) 1749 1750 ! Then perform representation change, following Eq.(25) in Hamann2005 1751 do ii=1,3 1752 strs_cartred(:,ii)=matmul(rprimd,strs(:,ii)) 1753 enddo 1754 do ii=1,3 1755 strs(ii,:)=matmul(rprimd,strs_cartred(ii,:)) 1756 enddo 1757 1758 !DEBUG 1759 ! rprimd_mod=rprimd 1760 ! rprimd_mod(1,1)=rprimd(1,1)*(one+strain) 1761 ! do ii=1,3 1762 ! strs_cartred(:,ii)=matmul(rprimd_mod,strs(:,ii)) 1763 ! enddo 1764 ! do ii=1,3 1765 ! strs(ii,:)=matmul(rprimd_mod,strs_cartred(ii,:)) 1766 ! enddo 1767 !ENDDEBUG 1768 1769 strs_intg(1,isp)=two*strs(1,1) ; strs_intg(2,isp)=two*strs(2,2) ; strs_intg(3,isp)=two*strs(3,3) 1770 strs_intg(4,isp)=strs(2,3)+strs(3,2) 1771 strs_intg(5,isp)=strs(1,3)+strs(3,1) 1772 strs_intg(6,isp)=strs(1,2)+strs(2,1) 1773 enddo 1774 strs_intg(:,:)=strs_intg(:,:)/dble(nfftot) 1775 !DEBUG 1776 ! strs_intg(:,:)=zero 1777 ! write(6,*)' calcdenmagsph : iatom,-strs_intg(1,1)*ucvol=',iatom,strs_intg(1,1)*ucvol 1778 ! call flush(6) 1779 ! stop 1780 !ENDDEBUG 1781 endif 1782 1783 if(nspden==2 .and. option/=11)then 1784 ! Specific treatment of collinear density, due to the storage mode. 1785 ! intgden_(1,iatom)= integral of up density 1786 ! intgden_(2,iatom)= integral of dn density 1787 intgden_(1,iatom)=intg(2) 1788 intgden_(2,iatom)=intg(1)-intg(2) 1789 if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then 1790 gr_intgden(:,1,iatom)=gr_intg(:,2) 1791 gr_intgden(:,2,iatom)=gr_intg(:,1)-gr_intg(:,2) 1792 endif 1793 if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then 1794 strs_intgden(:,1,iatom)=strs_intg(:,2) 1795 strs_intgden(:,2,iatom)=strs_intg(:,1)-strs_intg(:,2) 1796 endif 1797 else 1798 intgden_(1:nspden,iatom)=intg(1:nspden) 1799 if(present(gr_intgden).and. option<10 .and. ratsm2>tol12)then 1800 gr_intgden(:,1:nspden,iatom)=gr_intg(:,1:nspden) 1801 endif 1802 if(present(strs_intgden).and. option<10 .and. ratsm2>tol12)then 1803 strs_intgden(:,1:nspden,iatom)=strs_intg(:,1:nspden) 1804 endif 1805 endif 1806 1807 end do ! iatom 1808 !------------------------------------------- 1809 ! 1810 ! In case intgf2 must be computed, while the atoms overlap, a double loop over atoms is needed 1811 if(present(intgf2) .and. neighbor_overlap==1)then 1812 do iatom=1,natom 1813 do jatom=iatom,natom 1814 if(overlap_ij(iatom,jatom)/=0)then 1815 do i3=1,n3 1816 iz=mod(i3,n3) 1817 if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then 1818 izloc = ffti3_local(iz+1) - 1 1819 do i2=1,n2 1820 iy=mod(i2,n2) 1821 do i1=1,n1 1822 ix=mod(i1,n1) 1823 ifft_local=1+ix+n1*(iy+n2*izloc) 1824 intgf2(iatom,jatom)=intgf2(iatom,jatom)+fsm_atom(ifft_local,iatom)*fsm_atom(ifft_local,jatom) 1825 enddo 1826 enddo 1827 endif 1828 enddo ! i3 1829 intgf2(iatom,jatom)=intgf2(iatom,jatom)*ucvol/dble(nfftot) 1830 endif 1831 enddo 1832 if(iatom/=1)then 1833 do jatom=1,iatom-1 1834 intgf2(iatom,jatom)=intgf2(jatom,iatom) 1835 enddo 1836 endif 1837 enddo 1838 ABI_FREE(fsm_atom) 1839 endif 1840 1841 !MPI parallelization 1842 if(present(intgf2)) then 1843 if(mpi_enreg%nproc_fft>1)then 1844 call timab(48,1,tsec) 1845 call xmpi_sum(intgf2,mpi_enreg%comm_fft,ierr) 1846 call timab(48,2,tsec) 1847 end if 1848 end if 1849 1850 !------------------------------------------- 1851 1852 !MPI parallelization 1853 if(present(intgden) .or. option/=0) then 1854 if(mpi_enreg%nproc_fft>1)then 1855 call timab(48,1,tsec) 1856 call xmpi_sum(intgden_,mpi_enreg%comm_fft,ierr) 1857 call timab(48,2,tsec) 1858 end if 1859 if(present(intgden))intgden = intgden_ 1860 end if 1861 1862 if(present(gr_intgden) .and. option<10 .and. ratsm2>tol12) then 1863 if(mpi_enreg%nproc_fft>1)then 1864 call timab(48,1,tsec) 1865 call xmpi_sum(gr_intgden,mpi_enreg%comm_fft,ierr) 1866 call timab(48,2,tsec) 1867 end if 1868 end if 1869 1870 if(present(strs_intgden) .and. option<10 .and. ratsm2>tol12) then 1871 if(mpi_enreg%nproc_fft>1)then 1872 call timab(48,1,tsec) 1873 call xmpi_sum(strs_intgden,mpi_enreg%comm_fft,ierr) 1874 call timab(48,2,tsec) 1875 end if 1876 end if 1877 1878 !EB - Compute magnetization of the whole cell 1879 if(present(dentot) .or. present(rhomag))then 1880 rhomag_(:,:)=zero 1881 if(nspden==2) then 1882 do ifft=1,nfft 1883 jfft=1+cplex*(ifft-1) 1884 rhomag_(1:cplex,1)=rhomag_(1:cplex,1)+rhor(jfft:jfft+cplex-1,1) ! real & imag part of density 1885 rhomag_(1:cplex,2)=rhomag_(1:cplex,2)+2*rhor(jfft:jfft+cplex-1,2)-rhor(jfft:jfft+cplex-1,1) ! real & imag part of magnetization 1886 end do 1887 else if(nspden==4) then 1888 do ifft=1,nfft 1889 jfft=1+cplex*(ifft-1) 1890 rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)+rhor(jfft:jfft+cplex-1,1:nspden) 1891 end do 1892 end if 1893 1894 rhomag_(1:cplex,1:nspden)=rhomag_(1:cplex,1:nspden)*ucvol/dble(nfftot) 1895 1896 !MPI parallelization 1897 if(mpi_enreg%nproc_fft>1)then 1898 call timab(48,1,tsec) 1899 call xmpi_sum(rhomag_,mpi_enreg%comm_fft,ierr) 1900 call timab(48,2,tsec) 1901 end if 1902 1903 if(present(dentot))then 1904 dentot(:)=rhomag_(1,:) 1905 endif 1906 1907 if(present(rhomag))then 1908 rhomag(:,:)=rhomag_(:,:) 1909 endif 1910 endif 1911 1912 !DEBUG BUT KEEP 1913 if(.false.)then 1914 call printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,'DEN.vtk') 1915 endif 1916 !ENDDEBUG 1917 1918 end subroutine calcdenmagsph
m_dens/constrained_dft_free [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
constrained_dft_free
FUNCTION
Free the constrained_dft datastructure.
INPUTS
constrained_dft=datastructure that contain the needed information to enforce the density and magnetization constraints
OUTPUT
SOURCE
726 subroutine constrained_dft_free(constrained_dft) 727 728 !Arguments ------------------------------------ 729 !scalars 730 type(constrained_dft_t),intent(inout):: constrained_dft 731 732 !Local variables------------------------------- 733 734 ! *********************************************************************************************** 735 736 ABI_SFREE(constrained_dft%chrgat) 737 ABI_SFREE(constrained_dft%constraint_kind) 738 ABI_SFREE(constrained_dft%intgf2) 739 ABI_SFREE(constrained_dft%ratsph) 740 ABI_SFREE(constrained_dft%spinat) 741 ABI_SFREE(constrained_dft%typat) 742 ABI_SFREE(constrained_dft%ziontypat) 743 744 end subroutine constrained_dft_free
m_dens/constrained_dft_ini [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
constrained_dft_ini
FUNCTION
Initialize the constrained_dft datastructure. Mostly copying already available (dtset) information, but also computing intgf2
INPUTS
chrgat(natom) = target charge for each atom. Not always used, it depends on the value of constraint_kind
constraint_kind(ntypat)=for each type of atom, 0=no constraint,
1=fix only the magnetization direction, following spinat direction,
2=fix the magnetization vector to be the spinat one,
3=fix the magnetization amplitude to be the spinat one, but does not fix its direction
other future values will constrain the local atomic charge and possibly mix constraints if needed.
magconon=type of penalty function (so, not constrained DFT).
magcon_lambda=strength of the atomic spherical constraint
mpi_enreg=mpi structure with communicator info
natom=number of atoms
nfft=number of points in standard fft grid
ngfft=FFT grid dimensions
nspden = number of spin densities (1 2 or 4)
ntypat=number of types of atoms
ratsm=smearing width for ratsph
ratsph(ntypat)=radii for muffin tin spheres of each atom
rprimd=lattice vectors (dimensioned)
spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind
typat(natom)=types of atoms
xred(3,natom)=reduced atomic positions
ziontypat(ntypat)=ionic charge, per type of atom
OUTPUT
constrained_dft=datastructure that contain the needed information to enforce the density and magnetization constraints
Most of the data are simply copied from dtset, but also constrained_dft%intgf2(natom,natom) is computed from the available data.
SOURCE
632 subroutine constrained_dft_ini(chrgat,constrained_dft,constraint_kind,& 633 & magconon,magcon_lambda,mpi_enreg,natom,nfftf,ngfftf,nspden,ntypat,& 634 & ratsm,ratsph,rprimd,spinat,typat,xred,ziontypat) 635 636 !Arguments ------------------------------------ 637 !scalars 638 integer,intent(in) :: magconon,natom,nfftf,nspden,ntypat 639 real(dp),intent(in) :: magcon_lambda,ratsm 640 type(MPI_type),intent(in) :: mpi_enreg 641 type(constrained_dft_t),intent(out):: constrained_dft 642 !arrays 643 integer,intent(in) :: constraint_kind(ntypat) 644 integer,intent(in) :: ngfftf(18) 645 integer,intent(in) :: typat(natom) 646 real(dp),intent(in) :: chrgat(natom) 647 real(dp),intent(in) :: ratsph(ntypat) 648 real(dp),intent(in) :: rprimd(3,3) 649 real(dp),intent(in) :: spinat(3,natom) 650 real(dp),intent(in) :: xred(3,natom) 651 real(dp),intent(in) :: ziontypat(ntypat) 652 653 !Local variables------------------------------- 654 !scalars 655 integer :: cplex1=1 656 real(dp) :: ucvol 657 !arrays 658 real(dp), allocatable :: intgf2(:,:) ! natom,natom 659 real(dp), allocatable :: rhor_dum(:,:) ! nfftf,nspden 660 real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3) 661 662 ! *********************************************************************************************** 663 664 ABI_MALLOC(intgf2,(natom,natom)) 665 666 !We need the metric because it is needed in calcdenmagsph.F90 667 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 668 669 if(any(constraint_kind(:)/=0))then 670 !We need to precompute intgf2 671 ABI_MALLOC(rhor_dum,(nfftf,nspden)) 672 rhor_dum(:,:)=zero 673 call calcdenmagsph(mpi_enreg,natom,nfftf,ngfftf,nspden,ntypat,& 674 & ratsm,ratsph,rhor_dum,rprimd,typat,xred,0,cplex1,intgf2=intgf2) 675 ABI_FREE(rhor_dum) 676 else 677 intgf2=zero 678 endif 679 680 constrained_dft%magconon =magconon 681 constrained_dft%magcon_lambda =magcon_lambda 682 constrained_dft%natom =natom 683 constrained_dft%nfftf =nfftf 684 constrained_dft%ngfftf =ngfftf 685 constrained_dft%nspden =nspden 686 constrained_dft%ntypat =ntypat 687 constrained_dft%ratsm =ratsm 688 constrained_dft%rprimd =rprimd 689 690 ABI_MALLOC(constrained_dft%chrgat,(natom)) 691 ABI_MALLOC(constrained_dft%constraint_kind,(ntypat)) 692 ABI_MALLOC(constrained_dft%intgf2,(natom,natom)) 693 ABI_MALLOC(constrained_dft%ratsph,(ntypat)) 694 ABI_MALLOC(constrained_dft%spinat,(3,natom)) 695 ABI_MALLOC(constrained_dft%typat,(natom)) 696 ABI_MALLOC(constrained_dft%ziontypat,(ntypat)) 697 698 constrained_dft%chrgat =chrgat 699 constrained_dft%constraint_kind =constraint_kind 700 constrained_dft%intgf2 =intgf2 701 constrained_dft%ratsph =ratsph 702 constrained_dft%spinat =spinat 703 constrained_dft%typat =typat 704 constrained_dft%ziontypat =ziontypat 705 706 ABI_FREE(intgf2) 707 708 end subroutine constrained_dft_ini
m_dens/constrained_dft_t [ Types ]
NAME
constrained_dft_t
FUNCTION
Structure gathering the relevant information for constrained DFT calculations
SOURCE
65 type,public :: constrained_dft_t 66 67 !scalars 68 integer :: natom ! Number of atoms 69 integer :: nfftf ! Number of FFT grid points (for this processor) for the "fine" grid 70 integer :: nspden ! Number of spin-density components 71 integer :: ntypat ! Number of type of atoms 72 73 real(dp) :: magcon_lambda ! Strength of the atomic spherical constraint 74 real(dp) :: ratsm ! Smearing width for ratsph 75 76 integer :: magconon ! Turn on the penalty function constraint instead of the more powerful constrainedDFT algorithm 77 78 !arrays 79 80 integer :: ngfftf(18) ! Number of FFT grid points (for this processor) for the "fine" grid 81 82 integer,allocatable :: typat(:) 83 ! typat(natom) 84 ! Type of each natom 85 86 integer,allocatable :: constraint_kind(:) 87 ! constraint_kind(ntypat) 88 ! Constraint kind to be applied to each type of atom. See corresponding input variable 89 90 real(dp) :: rprimd(3,3) 91 ! Direct lattice vectors, Bohr units. 92 93 real(dp),allocatable :: chrgat(:) 94 ! chrgat(natom) 95 ! Target charge for each atom. Not always used, it depends on the value of constraint_kind 96 97 real(dp),allocatable :: intgf2(:,:) 98 ! intgf2(natom,natom) 99 ! Overlap of the spherical integrating functions, for each atom. 100 ! Initialized using some xred values, will not change during the SCF cycles, except for exotic algorithms, not in production, 101 102 real(dp),allocatable :: ratsph(:) 103 ! ratsph(ntypat) 104 ! Radius of the atomic sphere for each type of atom 105 106 real(dp),allocatable :: spinat(:,:) 107 ! spinat(3,natom) 108 ! Target magnetization for each atom. Possibly only the direction or the magnitude, depending on constraint_kind 109 110 real(dp),allocatable :: ziontypat(:) 111 ! ziontypat(ntypat) 112 ! Ionic charge, per type of atom 113 114 end type constrained_dft_t
m_dens/constrained_residual [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
constrained_residual
FUNCTION
Recompute the residual to take into account the constraints, within constrained DFT. The kind of constraint is given by constraint_kind, and the target values are given by spinat, for the local atomic magnetization, and chrgat minus the ionic charge, for the local atomic charge.
INPUTS
c_dft <type(constrained_dft_t)>=datastructure for the information related to constrained DFT ! chrgat(natom) = target charge for each atom. Not always used, it depends on the value of constraint_kind ! constraint_kind(ntypat)=for each type of atom, 0=no constraint, ! 1=fix only the magnetization direction, following spinat direction, ! 2=fix the magnetization vector to be the spinat one, ! 3=fix the magnetization amplitude to be the spinat one, but does not fix its direction ! other future values will constrain the local atomic charge and possibly mix constraints if needed. ! intgf2(natom,natom)=(precomputed) overlap of the spherical integration functions for each atom in a sphere of radius ratsph. ! magcon_lambda=strength of the atomic spherical constraint ! natom=number of atoms ! nfftf=number of points in fine fft grid ! ngfftf=FFT grid dimensions ! nspden = number of spin densities (1 2 or 4) ! ntypat=number of types of atoms ! ratsm=smearing width for ratsph ! ratsph(ntypat)=radii for muffin tin spheres of each atom ! rprimd=lattice vectors (dimensioned) ! spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind ! typat(natom)=types of atoms ! ziontypat(ntypat)= ionic charge, per type of atom mpi_enreg=mpi structure with communicator info rhor(nfft,nspden)=array for electron density in el./bohr**3. At output it will be constrained. xred(3,natom)=reduced atomic positions
OUTPUT
e_constrained_dft=correction to the total energy, to make it variational grcondft(3,natom)=d(E_constrained_DFT)/d(xred) (hartree) intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints. strscondft(6)=stress due to constraints = -d(E_constrained_DFT)/d(strain) / ucvol (hartree/Bohr^3)
SIDE EFFECTS
vresid(nfft,nspden)==array for potential residual in real space At output it will be modified: projected onto the space orthogonal to the atomic spherical functions (if there is a related constrained, and augmented by such atomic spherical functions multiplied by the difference between the actual integrated charge or magnetization and the target ones.
SOURCE
797 subroutine constrained_residual(c_dft,e_constrained_dft,grcondft,intgres,mpi_enreg,rhor,strscondft,vresid,xred) 798 799 !Arguments ------------------------------------ 800 !scalars 801 real(dp),intent(out) :: e_constrained_dft 802 type(constrained_dft_t),intent(in) :: c_dft 803 type(MPI_type),intent(in) :: mpi_enreg 804 !arrays 805 real(dp),intent(out) :: grcondft(:,:) ! 3,natom 806 real(dp),intent(out) :: intgres(:,:) ! nspden,natom 807 real(dp),intent(in) :: rhor(c_dft%nfftf,c_dft%nspden) 808 real(dp),intent(out) :: strscondft(6) 809 real(dp),intent(inout) :: vresid(c_dft%nfftf,c_dft%nspden) 810 real(dp),intent(in) :: xred(3,c_dft%natom) 811 812 !Local variables------------------------------- 813 !scalars 814 integer :: conkind,iatom,ii,jatom,info,natom,nfftf,nspden,ntypat,option 815 integer :: cplex1=1 816 real(dp) :: intgd,intgden_norm,intgden_proj,intgres_proj,norm,scprod 817 !arrays 818 integer :: ipiv(c_dft%natom) 819 real(dp) :: corr_denmag(4),gr_intgd(3),strs_intgd(6) 820 real(dp), allocatable :: coeffs_constr_dft(:,:) ! nspden,natom 821 real(dp), allocatable :: gr_intgden(:,:,:) ! 3,nspden,natom 822 real(dp), allocatable :: intgden(:,:) ! nspden,natom 823 real(dp), allocatable :: intgden_delta(:,:) ! nspden,natom 824 real(dp), allocatable :: intgres_tmp(:,:) ! nspden,natom 825 real(dp), allocatable :: intgr(:,:) ! natom,nspden 826 real(dp), allocatable :: strs_intgden(:,:,:) ! 6,nspden,natom 827 real(dp) :: intgf2(c_dft%natom,c_dft%natom),rhomag(2,c_dft%nspden),work(2*c_dft%natom) 828 real(dp) :: intgden_normed(3) 829 real(dp) :: spinat_normed(3) 830 831 ! *********************************************************************************************** 832 833 !DEBUG 834 !write(std_out,*) ' constrained_residual : enter ' 835 !ENDDEBUG 836 837 natom=c_dft%natom 838 nfftf=c_dft%nfftf 839 nspden=c_dft%nspden 840 ntypat=c_dft%ntypat 841 842 !We need the integrated magnetic moments 843 ABI_MALLOC(intgden,(nspden,natom)) 844 ABI_MALLOC(gr_intgden,(3,nspden,natom)) 845 ABI_MALLOC(strs_intgden,(6,nspden,natom)) 846 call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,c_dft%ratsm,c_dft%ratsph,rhor,c_dft%rprimd,c_dft%typat,& 847 & xred,1,cplex1,intgden=intgden,gr_intgden=gr_intgden,rhomag=rhomag,strs_intgden=strs_intgden) 848 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,std_out,1,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat) 849 850 !DEBUG 851 !write(std_out,*) ' intgden(1:nspden,1:natom)=',intgden(1:nspden,1:natom) 852 !ENDDEBUG 853 854 !We need the integrated residuals 855 ABI_MALLOC(intgres_tmp,(nspden,natom)) 856 intgres_tmp(:,:)=zero 857 call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,& 858 & c_dft%ratsm,c_dft%ratsph,vresid,c_dft%rprimd,c_dft%typat,xred,11,cplex1,intgden=intgres_tmp,rhomag=rhomag) 859 860 !DEBUG 861 !write(std_out,*) ' intgres_tmp(1:nspden,1:natom)=',intgres_tmp(1:nspden,1:natom) 862 !ENDDEBUG 863 864 !Make the proper combination of intgres_tmp, to single out the scalar potential residual and the magnetic field potential residuals for x,y,z. 865 intgres(:,:)=zero 866 do iatom=1,natom 867 if(nspden==1)then 868 intgres(1,iatom)=intgres_tmp(1,iatom) 869 else if(nspden==2)then 870 intgres(1,iatom)=half*(intgres_tmp(1,iatom)+intgres_tmp(2,iatom)) 871 intgres(2,iatom)=half*(intgres_tmp(1,iatom)-intgres_tmp(2,iatom)) 872 else if(nspden==4)then 873 !Change the potential residual to the density+magnetization convention 874 intgres(1,iatom)=half*(intgres_tmp(1,iatom)+intgres_tmp(2,iatom)) 875 intgres(2,iatom)= intgres_tmp(3,iatom) 876 intgres(3,iatom)=-intgres_tmp(4,iatom) 877 intgres(4,iatom)=half*(intgres_tmp(1,iatom)-intgres_tmp(2,iatom)) 878 endif 879 conkind=c_dft%constraint_kind(c_dft%typat(iatom)) 880 if(conkind <10)intgres(1,iatom)=zero 881 if( mod(conkind,10)==0 .and. nspden>1)intgres(2:nspden,iatom)=zero 882 enddo 883 !Print the potential residuals 884 call prtdenmagsph(cplex1,intgres,natom,nspden,ntypat,std_out,11,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat) 885 ABI_FREE(intgres_tmp) 886 887 !Also exchanges the spin and atom indices to prepare the solution of the linear system of equation 888 ABI_MALLOC(intgr,(natom,nspden)) 889 do iatom=1,natom 890 intgr(iatom,1:nspden)=intgres(1:nspden,iatom) 891 enddo 892 893 !In case there is an overlap between spheres, must solve the linear system of equations 894 !and take into account non-diagonal elements only for the set of atoms for which there is a constraint. 895 !This can be different for the 896 !charge residual or for the spin residual (but for the spin constraints, the set of atoms is inclusive of all constraints) 897 do ii=1,2 ! Charge, then spin 898 intgf2(:,:)=zero 899 do iatom=1,natom 900 intgf2(iatom,iatom)=c_dft%intgf2(iatom,iatom) 901 enddo 902 do iatom=1,natom-1 903 do jatom=iatom,natom 904 if(c_dft%intgf2(iatom,jatom)>tol8)then 905 !In the charge case, must have both atoms with constraints bigger than 10 906 if(ii==1)then 907 if(c_dft%constraint_kind(c_dft%typat(iatom))>=10 .and. & 908 & c_dft%constraint_kind(c_dft%typat(jatom))>=10 )then 909 intgf2(iatom,jatom)=c_dft%intgf2(iatom,jatom) 910 intgf2(jatom,iatom)=c_dft%intgf2(iatom,jatom) 911 endif 912 endif 913 !In the spin case, must have both atoms with constraints not ending with 0 914 if(ii==2)then 915 if(mod(c_dft%constraint_kind(c_dft%typat(iatom)),10)/=0 .and. & 916 & mod(c_dft%constraint_kind(c_dft%typat(jatom)),10)/=0 )then 917 intgf2(iatom,jatom)=c_dft%intgf2(iatom,jatom) 918 intgf2(jatom,iatom)=c_dft%intgf2(iatom,jatom) 919 endif 920 endif 921 endif 922 enddo 923 enddo 924 925 !Solve the linear system of equation, for the different spins 926 call dsytrf('U',natom,intgf2,natom,ipiv,work,2*natom,info) 927 ! call dsytri('U',natom,intgf2,natom,ipiv,work,info) 928 if(ii==1)then 929 call dsytrs('U',natom,1,intgf2,natom,ipiv,intgr,natom,info) 930 else if(ii==2 .and. nspden>1)then 931 call dsytrs('U',natom,nspden-1,intgf2,natom,ipiv,intgr(1:natom,2:nspden),natom,info) 932 endif 933 934 !Store the new residuals 935 do iatom=1,natom 936 if(ii==1)intgres(1,iatom)=intgr(iatom,1) 937 if(ii==2)intgres(2:nspden,iatom)=intgr(iatom,2:nspden) 938 enddo 939 940 enddo 941 942 !DEBUG 943 !write(std_out,*) ' after multiplication by ftt-1 , so, torque :' 944 !write(std_out,*) ' intgres(1:nspden,1:natom)=',intgres(1:nspden,1:natom) 945 !ENDDEBUG 946 947 !Compute the delta of the integrated dens with respect to the target 948 !Compute the energy correction, to make the energy functional variational 949 !Also projects the residual in case constraint_kind 2 950 e_constrained_dft=zero 951 grcondft=zero 952 strscondft=zero 953 ABI_MALLOC(intgden_delta,(nspden,natom)) 954 intgden_delta(:,:)=zero 955 do iatom=1,natom 956 957 ! The integrated density must be in the total density+magnetization representation 958 if(nspden==2)then 959 intgd =intgden(1,iatom)+intgden(2,iatom) 960 intgden(2,iatom)=intgden(1,iatom)-intgden(2,iatom) 961 intgden(1,iatom)=intgd 962 do ii=1,3 963 gr_intgd(ii) =gr_intgden(ii,1,iatom)+gr_intgden(ii,2,iatom) 964 gr_intgden(ii,2,iatom)=gr_intgden(ii,1,iatom)-gr_intgden(ii,2,iatom) 965 gr_intgden(ii,1,iatom)=gr_intgd(ii) 966 enddo 967 do ii=1,6 968 strs_intgd(ii) =strs_intgden(ii,1,iatom)+strs_intgden(ii,2,iatom) 969 strs_intgden(ii,2,iatom)=strs_intgden(ii,1,iatom)-strs_intgden(ii,2,iatom) 970 strs_intgden(ii,1,iatom)=strs_intgd(ii) 971 enddo 972 endif 973 974 !Comparison with the target value, and computation of the correction in terms of density and magnetization coefficients. 975 conkind=c_dft%constraint_kind(c_dft%typat(iatom)) 976 977 if(conkind >=10)then 978 !The electronic constraint is such that the ziontypat charge minus (the electronic charge is negative) the atomic electronic density 979 !intgden gives the target charge chrgat. 980 intgden_delta(1,iatom)=intgden(1,iatom)+c_dft%chrgat(iatom)-c_dft%ziontypat(c_dft%typat(iatom)) 981 ! Uses the usual electronic charge definition, instead of the total nucleus-electronic charge 982 ! intgden_delta(1,iatom)=intgden(1,iatom)-c_dft%chrgat(iatom) 983 endif 984 985 if( mod(conkind,10)==1 .and. nspden>1)then 986 !Fix the different components of the magnetization vector 987 if(nspden==2)intgden_delta(2,iatom)=intgden(2,iatom)-c_dft%spinat(3,iatom) 988 if(nspden==4)intgden_delta(2:4,iatom)=intgden(2:4,iatom)-c_dft%spinat(1:3,iatom) 989 else if( ( mod(conkind,10)>=2 .and. mod(conkind,10)<=4) .and. nspden>1)then 990 norm = sqrt(sum(c_dft%spinat(:,iatom)**2)) 991 if (norm > tol10) then 992 if( mod(conkind,10)==2 )then 993 !Fix the axis of the magnetization vector 994 if(nspden==4)then 995 spinat_normed(:) = c_dft%spinat(:,iatom) / norm 996 !Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector 997 !This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector 998 intgden_proj=spinat_normed(1)*intgden(2,iatom)+ & 999 & spinat_normed(2)*intgden(3,iatom)+ & 1000 & spinat_normed(3)*intgden(4,iatom) 1001 intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)-spinat_normed(1:3)*intgden_proj 1002 !Also projects the residual, so that the usual optimization is done is the largest possible space 1003 intgres_proj=spinat_normed(1)*intgres(2,iatom)+ & 1004 & spinat_normed(2)*intgres(3,iatom)+ & 1005 & spinat_normed(3)*intgres(4,iatom) 1006 intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-spinat_normed(1:3)*intgres_proj 1007 else if(nspden==2)then 1008 !The direction must be correct, collinear, so no change. 1009 intgden_delta(2,iatom)=zero 1010 endif 1011 else if( mod(conkind,10)==3 )then 1012 !Fix the amplitude of the magnetization vector 1013 intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2)) 1014 intgden_delta(2:nspden,iatom)=(one-norm/intgden_norm)*intgden(2:nspden,iatom) 1015 else if( mod(conkind,10)==4 )then 1016 !Fix the direction of the magnetization vector 1017 if(nspden==4)then 1018 spinat_normed(:) = c_dft%spinat(:,iatom) / norm 1019 intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2)) 1020 intgden_normed(:) = intgden(2:nspden,iatom)/intgden_norm 1021 !Calculate the difference vector between the actual magnetization vectoro, times the scalar product between the 1022 !directions of present magnetization and target one, and the target magnetization direction renormalized by the magnetization length. 1023 !See notes 12 October 2021 1024 !DEBUG First possibility 1025 ! intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)-spinat_normed(1:3)*intgden_norm 1026 !ENDDEBUG 1027 !DEBUG Second possibility 1028 scprod=sum(spinat_normed(:)*intgden_normed(:)) 1029 intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom)*scprod-spinat_normed(1:3)*intgden_norm 1030 !ENDDEBUG 1031 !Also projects the residual. There might be some misalignement of the intgden_delta with the correct space of allowed variations 1032 !(not exactly perpendicular to spinat_normed or intgden_normed, or even frankly not at all perpendicular), 1033 !but when close to fulfilling the constraint, this difference becomes negligible. 1034 !The present choice, couple with the above definition of intgden_delta aligns both. 1035 !DEBUG First possibility 1036 intgres_proj=intgden_normed(1)*intgres(2,iatom)+ & 1037 & intgden_normed(2)*intgres(3,iatom)+ & 1038 & intgden_normed(3)*intgres(4,iatom) 1039 intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-intgden_normed(1:3)*intgres_proj 1040 !ENDDEBUG 1041 !DEBUG Second possibility 1042 ! intgres_proj=spinat_normed(1)*intgres(2,iatom)+ & 1043 !& spinat_normed(2)*intgres(3,iatom)+ & 1044 !& spinat_normed(3)*intgres(4,iatom) 1045 ! intgres(2:nspden,iatom)=intgres(2:nspden,iatom)-spinat_normed(1:3)*intgres_proj 1046 !ENDDEBUG 1047 else if(nspden==2)then 1048 !This case is ill-defined ... What follows is rather arbitrary. 1049 !When the actual magnetization and the target magnetization have same sign, the intgden_delta is zero. 1050 ! Otherwise, intgres is set to the difference between the present intgden and the target 1051 intgden_delta(2,iatom)=zero 1052 if(c_dft%spinat(2,iatom)*intgden(2,iatom)<-tol10)intgden_delta(2,iatom)=c_dft%spinat(2,iatom)-intgden(2,iatom) 1053 endif 1054 endif 1055 else 1056 !In this case (norm of constraint vanishes), we set the atomic magnetization to zero. 1057 intgden_delta(2:nspden,iatom)=intgden(2:nspden,iatom) 1058 endif 1059 end if 1060 ! Lagrange energy contribution. Note that intgres is the derivative with respect to the chrgat and spinat constraints, 1061 ! because chrgat and spinat have directly been used in the definition of intgden_delta. 1062 ! WOOPS : chrgat comes with a POSITIVE sign in intgden_delta ?!?! 1063 e_constrained_dft=e_constrained_dft-sum(intgden_delta(:,iatom)*intgres(:,iatom)) 1064 do ii=1,3 1065 grcondft(ii,iatom)=grcondft(ii,iatom)-sum(gr_intgden(ii,:,iatom)*intgres(:,iatom)) 1066 enddo 1067 ! For the stress, this is the place where the summation over atoms is performed. 1068 do ii=1,6 1069 strscondft(ii)=strscondft(ii)-sum(strs_intgden(ii,:,iatom)*intgres(:,iatom)) 1070 enddo 1071 1072 !DEBUG 1073 ! write(6,*)' calcdenmagsph/constrained_residual, line 1058 : iatom=',iatom 1074 ! write(6,*)' e_constrained_dft,intgden_delta(:,iatom),intgres(:,iatom)=',e_constrained_dft,intgden_delta(:,iatom),intgres(:,iatom) 1075 ! write(6,*)' grcondft(1,iatom),gr_intgden(1,:,iatom),intgres(:,iatom)=',grcondft(1,iatom),gr_intgden(1,:,iatom),intgres(:,iatom) 1076 ! write(6,*)' strscondft(1),strs_intgden(1,:,iatom),intgres(:,iatom)=',strscondft(1),strs_intgden(1,:,iatom),intgres(:,iatom) 1077 ! call flush(6) 1078 ! stop 1079 !ENDDEBUG 1080 1081 enddo 1082 1083 ABI_FREE(gr_intgden) 1084 ABI_FREE(strs_intgden) 1085 1086 ABI_MALLOC(coeffs_constr_dft,(nspden,natom)) 1087 coeffs_constr_dft=zero 1088 1089 !With the delta of the integrated density and the atomic residual, compute the atomic correction to be applied to the potential 1090 !See the Eqs.(31), (33) and (34) of notes. 1091 do iatom=1,natom 1092 1093 !Computation of the correction in terms of density and magnetization coefficients. 1094 conkind=c_dft%constraint_kind(c_dft%typat(iatom)) 1095 corr_denmag(:)=zero 1096 1097 if(conkind >=10)then 1098 corr_denmag(1)=intgden_delta(1,iatom)*c_dft%magcon_lambda - intgres(1,iatom) 1099 endif 1100 1101 if( mod(conkind,10)==1 .and. nspden>1)then 1102 1103 !Fix the different components of the magnetization vector 1104 if(nspden==2)corr_denmag(2)=intgden_delta(2,iatom)*c_dft%magcon_lambda - intgres(2,iatom) 1105 if(nspden==4)corr_denmag(2:4)=intgden_delta(2:4,iatom)*c_dft%magcon_lambda - intgres(2:4,iatom) 1106 1107 else if( ( mod(conkind,10)>=2 .and. mod(conkind,10)<=4) .and. nspden>1)then 1108 1109 norm = sqrt(sum(c_dft%spinat(:,iatom)**2)) 1110 if (norm > tol10) then 1111 1112 if( mod(conkind,10)==2 .or. mod(conkind,10)==4)then 1113 !Fix the axis/direction of the magnetization vector 1114 if(nspden==4 .or. mod(conkind,10)==4)then 1115 corr_denmag(2:nspden)=intgden_delta(2:nspden,iatom)*c_dft%magcon_lambda -intgres(2:nspden,iatom) 1116 else if(nspden==2 .and. mod(conkind,10)==2)then 1117 !The direction must be correct, collinear, so no change. 1118 corr_denmag(2)=zero 1119 endif 1120 1121 else if( mod(conkind,10)==3 )then 1122 1123 !Fix the amplitude of the magnetization vector 1124 !This is a special case, one does not work (at present) with the intgden_delta, while one should ... 1125 intgden(2:nspden,iatom)=intgden(2:nspden,iatom) - intgres(2:nspden,iatom)/c_dft%magcon_lambda 1126 intgden_norm = sqrt(sum(intgden(2:nspden,iatom)**2)) 1127 corr_denmag(2:nspden)=(one-norm/intgden_norm)*intgden(2:nspden,iatom) 1128 corr_denmag(2:nspden)=corr_denmag(2:nspden)*c_dft%magcon_lambda 1129 1130 endif 1131 1132 else 1133 !In this case, we set the atomic magnetization to zero. 1134 corr_denmag(2:nspden)=intgden_delta(2:nspden,iatom)*c_dft%magcon_lambda - intgres(2,iatom) 1135 endif 1136 1137 end if 1138 1139 !Convert from density/magnetization constraint residual to actual coefficient that will multiply the spherical function for the potential 1140 if(nspden==1)then 1141 !From charge to potential 1142 coeffs_constr_dft(1,iatom)=corr_denmag(1) 1143 else if(nspden==2 .or. nspden==4)then 1144 !From charge and magnetization to potential 1145 coeffs_constr_dft(1,iatom)=corr_denmag(1)+corr_denmag(nspden) 1146 coeffs_constr_dft(2,iatom)=corr_denmag(1)-corr_denmag(nspden) 1147 if(nspden==4)then 1148 coeffs_constr_dft(3,iatom)= corr_denmag(2) 1149 coeffs_constr_dft(4,iatom)=-corr_denmag(3) 1150 endif 1151 endif 1152 1153 enddo 1154 1155 !Now compute the new residual, by adding the spherical functions 1156 option=1 1157 call add_atomic_fcts(natom,nspden,c_dft%rprimd,mpi_enreg,nfftf,c_dft%ngfftf,ntypat,option,& 1158 & c_dft%ratsph,c_dft%ratsm,c_dft%typat,coeffs_constr_dft,vresid,xred) ! Ben change: add ratsm 1159 1160 ABI_FREE(coeffs_constr_dft) 1161 ABI_FREE(intgden) 1162 ABI_FREE(intgden_delta) 1163 ABI_FREE(intgr) 1164 1165 end subroutine constrained_residual
m_dens/dens_hirsh [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
dens_hirsh
FUNCTION
Compute the Hirshfeld charges
INPUTS
mpoint=Maximum number of points in radial meshes. radii(mpoint, ntypat)=Radial meshes for each type aeden(mpoint, nytpat)=All-electron densities. npoint(ntypat)=The number of the last point with significant density is stored in npoint(itypat) minimal_den=Tolerance on the minum value of the density grid_den(nrx,nry,nrz)= density on the grid natom = number of atoms in the unit cell nrx,nry,nrz= number of points in the grid for the three directions ntypat=number of types of atoms in unit cell. rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(natom)=type of each atom xcart(3,natom) = different positions of the atoms in the unit cell zion=(ntypat)gives the ionic charge for each type of atom prtcharge=1 to write the Hirshfeld charge decomposition
OUTPUT
hcharge(natom), hden(natom), hweight(natom)= Hirshfeld charges, densities, weights.
SOURCE
150 subroutine dens_hirsh(mpoint,radii,aeden,npoint,minimal_den,grid_den, & 151 natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,prtcharge,hcharge,hden,hweight) 152 153 !Arguments ------------------------------------ 154 !scalars 155 integer,intent(in) :: natom,nrx,nry,nrz,ntypat,prtcharge,mpoint 156 real(dp),intent(in) :: minimal_den 157 !arrays 158 integer,intent(in) :: typat(natom),npoint(ntypat) 159 real(dp),intent(in) :: grid_den(nrx,nry,nrz),rprimd(3,3),zion(ntypat) 160 real(dp),intent(in) :: xcart(3,natom) 161 real(dp),intent(in) :: radii(mpoint,ntypat),aeden(mpoint,ntypat) 162 real(dp),intent(out) :: hcharge(natom),hden(natom),hweight(natom) 163 164 !Local variables ------------------------- 165 !scalars 166 integer :: i1,i2,i3,iatom,icell,igrid,ii,inmax,inmin,istep,itypat 167 integer :: k1,k2,k3,mcells,nfftot,ngoodpoints,npt 168 real(dp) :: aa,bb,coeff1,coeff2,coeff3,den,factor,h_inv,hh,maxrad 169 real(dp) :: rr,rr2,total_charge,total_weight,total_zion,ucvol 170 real(dp) :: yp1,ypn 171 !arrays 172 integer :: highest(3),lowest(3) 173 integer,allocatable :: ncells(:) 174 real(dp) :: coordat(3),coord23_1,coord23_2,coord23_3,diff1,diff2,diff3,gmet(3,3),gprimd(3,3),rmet(3,3) 175 real(dp) :: vperp(3),width(3) 176 real(dp),allocatable :: coord1(:,:),local_den(:,:,:,:) 177 real(dp),allocatable :: step(:,:),sum_den(:,:,:) 178 real(dp),allocatable :: xcartcells(:,:,:),xred(:,:),yder2(:) 179 180 ! ********************************************************************* 181 182 !1. Read the 1D all-electron atomic files 183 !Store the radii in radii(:,itypat), and the all-electron 184 !densities in aeden(:,itypat). The number of the last 185 !point with significant density is stored in npoint(itypat) 186 187 !2. Compute the list of atoms that are sufficiently close to the cell 188 189 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 190 nfftot=nrx*nry*nrz 191 192 ABI_MALLOC(xred,(3,natom)) 193 call xcart2xred(natom,rprimd,xcart,xred) 194 195 !Compute the widths of the cell 196 !First width : perpendicular vector length 197 vperp(:)=rprimd(:,1)-rprimd(:,2)*rmet(1,2)/rmet(2,2) -rprimd(:,3)*rmet(1,3)/rmet(3,3) 198 width(1)=sqrt(dot_product(vperp,vperp)) 199 !Second width 200 vperp(:)=rprimd(:,2)-rprimd(:,1)*rmet(2,1)/rmet(1,1) -rprimd(:,3)*rmet(2,3)/rmet(3,3) 201 width(2)=sqrt(dot_product(vperp,vperp)) 202 !Third width 203 vperp(:)=rprimd(:,3)-rprimd(:,1)*rmet(3,1)/rmet(1,1) -rprimd(:,2)*rmet(3,2)/rmet(2,2) 204 width(3)=sqrt(dot_product(vperp,vperp)) 205 206 !Compute the number of cells that will make up the supercell 207 ABI_MALLOC(ncells,(natom)) 208 mcells=0 209 do iatom=1,natom 210 itypat=typat(iatom) 211 maxrad=radii(npoint(itypat),itypat) 212 ! Compute the lower and higher indices of the supercell 213 ! for this atom 214 do ii=1,3 215 lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii)) 216 highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1) 217 ! Next coding, still incorrect 218 ! lowest(ii)=floor(xred(ii,iatom)-maxrad/width(ii))-1 219 ! highest(ii)=ceiling(xred(ii,iatom)+maxrad/width(ii))+1 220 ! Old coding incorrect 221 ! lowest(ii)=ceiling(-xred(ii,iatom)-maxrad/width(ii)) 222 ! highest(ii)=floor(-xred(ii,iatom)+maxrad/width(ii)+1) 223 end do 224 ncells(iatom)=(highest(1)-lowest(1)+1)* & 225 & (highest(2)-lowest(2)+1)* & 226 & (highest(3)-lowest(3)+1) 227 ! DEBUG 228 ! write(std_out,*)' maxrad=',maxrad 229 ! write(std_out,*)' lowest(:)=',lowest(:) 230 ! write(std_out,*)' highest(:)=',highest(:) 231 ! write(std_out,*)' ncells(iatom)=',ncells(iatom) 232 ! ENDDEBUG 233 end do 234 mcells=maxval(ncells(:)) 235 236 !Compute, for each atom, the set of image atoms in the whole supercell 237 ABI_MALLOC(xcartcells,(3,mcells,natom)) 238 do iatom=1,natom 239 itypat=typat(iatom) 240 maxrad=radii(npoint(itypat),itypat) 241 ! Compute the lower and higher indices of the supercell 242 ! for this atom 243 244 do ii=1,3 245 lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii)) 246 highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1) 247 end do 248 icell=0 249 do i1=lowest(1),highest(1) 250 do i2=lowest(2),highest(2) 251 do i3=lowest(3),highest(3) 252 icell=icell+1 253 xcartcells(:,icell,iatom)=xcart(:,iatom)+i1*rprimd(:,1)+i2*rprimd(:,2)+i3*rprimd(:,3) 254 end do 255 end do 256 end do 257 end do 258 259 !Compute, for each atom, the all-electron pro-atom density 260 !at each point in the primitive cell 261 ABI_MALLOC(local_den,(nrx,nry,nrz,natom)) 262 ABI_MALLOC(step,(2,mpoint)) 263 ABI_MALLOC(yder2,(mpoint)) 264 ABI_MALLOC(coord1,(3,nrx)) 265 coeff1=one/nrx 266 coeff2=one/nry 267 coeff3=one/nrz 268 269 do iatom=1,natom 270 itypat=typat(iatom) 271 npt=npoint(itypat) 272 maxrad=radii(npt,itypat) 273 ! write(std_out,*) 274 ! write(std_out,'(a,i4)' )' hirsh : accumulating density for atom ',iatom 275 ! write(std_out,*)' ncells(iatom)=',ncells(iatom) 276 do istep=1,npt-1 277 step(1,istep)=radii(istep+1,itypat) - radii(istep,itypat) 278 step(2,istep)=one/step(1,istep) 279 end do 280 ! Approximate first derivative for small radii 281 yp1=(aeden(2,itypat)-aeden(1,itypat))/(radii(2,itypat)-radii(1,itypat)) 282 ypn=zero 283 call spline(radii(1:npt,itypat),aeden(1:npt,itypat),npt,yp1,ypn,yder2) 284 285 local_den(:,:,:,iatom)=zero 286 287 ! Big loop on the cells 288 do icell=1,ncells(iatom) 289 ! write(std_out,*)' icell=',icell 290 coordat(:)=xcartcells(:,icell,iatom) 291 292 ! Big loop on the grid points 293 do k1 = 1,nrx 294 coord1(:,k1)=rprimd(:,1)*(k1-1)*coeff1 295 end do 296 do k3 = 1, nrz 297 do k2 = 1, nry 298 coord23_1=rprimd(1,2)*(k2-1)*coeff2+rprimd(1,3)*(k3-1)*coeff3-coordat(1) 299 coord23_2=rprimd(2,2)*(k2-1)*coeff2+rprimd(2,3)*(k3-1)*coeff3-coordat(2) 300 coord23_3=rprimd(3,2)*(k2-1)*coeff2+rprimd(3,3)*(k3-1)*coeff3-coordat(3) 301 do k1 = 1, nrx 302 diff1=coord1(1,k1)+coord23_1 303 diff2=coord1(2,k1)+coord23_2 304 diff3=coord1(3,k1)+coord23_3 305 rr2=diff1**2+diff2**2+diff3**2 306 if(rr2<maxrad**2)then 307 308 rr=sqrt(rr2) 309 ! Find the index of the radius by bissection 310 if (rr < radii(1,itypat)) then 311 ! Linear extrapolation 312 den=aeden(1,itypat)+(rr-radii(1,itypat))/(radii(2,itypat)-radii(1,itypat))& 313 & *(aeden(2,itypat)-aeden(1,itypat)) 314 else 315 ! Use the spline interpolation 316 ! Find the index of the radius by bissection 317 inmin=1 318 inmax=npt 319 igrid=1 320 do 321 if(inmax-inmin==1)exit 322 igrid=(inmin+inmax)/2 323 if(rr>=radii(igrid,itypat))then 324 inmin=igrid 325 else 326 inmax=igrid 327 end if 328 end do 329 igrid=inmin 330 ! write(std_out,*)' igrid',igrid 331 332 hh=step(1,igrid) 333 h_inv=step(2,igrid) 334 aa= (radii(igrid+1,itypat)-rr)*h_inv 335 bb= (rr-radii(igrid,itypat))*h_inv 336 den = aa*aeden(igrid,itypat) + bb*aeden(igrid+1,itypat) & 337 & +( (aa*aa*aa-aa)*yder2(igrid) & 338 & +(bb*bb*bb-bb)*yder2(igrid+1) ) *hh*hh*sixth 339 end if ! Select small radius or spline 340 341 local_den(k1,k2,k3,iatom)=local_den(k1,k2,k3,iatom)+den 342 end if ! dist2<maxrad 343 344 end do ! k1 345 end do ! k2 346 end do ! k3 347 348 end do ! icell 349 end do ! iatom 350 351 !Compute, the total all-electron density at each point in the primitive cell 352 ABI_MALLOC(sum_den,(nrx,nry,nrz)) 353 sum_den(:,:,:)=zero 354 do iatom=1,natom 355 sum_den(:,:,:)=sum_den(:,:,:)+local_den(:,:,:,iatom) 356 end do 357 358 !Accumulate the integral of the density, to get Hirshfeld charges 359 !There is a minus sign because the electron has a negative charge 360 ngoodpoints = 0 361 hcharge(:)=zero 362 hweight(:)=zero 363 do k3=1,nrz 364 do k2=1,nry 365 do k1=1,nrx 366 ! Use minimal_den in order to avoid divide by zero 367 if (abs(sum_den(k1,k2,k3)) > minimal_den) then 368 ngoodpoints = ngoodpoints+1 369 factor=grid_den(k1,k2,k3)/(sum_den(k1,k2,k3)+minimal_den) 370 do iatom=1,natom 371 hden(iatom)=hden(iatom)+local_den(k1,k2,k3,iatom) 372 hcharge(iatom)=hcharge(iatom)-local_den(k1,k2,k3,iatom)*factor 373 hweight(iatom)=hweight(iatom)+local_den(k1,k2,k3,iatom)/(sum_den(k1,k2,k3)+minimal_den) 374 end do 375 end if 376 end do 377 end do 378 end do 379 380 !DEBUG 381 !do iatom=1,natom 382 !write(std_out,'(i9,3es17.6)' )iatom,hden(iatom),hcharge(iatom),hweight(iatom) 383 !end do 384 !ENDDEBUG 385 386 hcharge(:)=hcharge(:)*ucvol/dble(nfftot) 387 hweight(:)=hweight(:)/dble(nfftot) 388 389 !Check on the total charge 390 total_zion=sum(zion(typat(1:natom))) 391 total_charge=sum(hcharge(1:natom)) 392 total_weight=sum(hweight(1:natom)) 393 394 !DEBUG 395 !write(std_out,*)' ngoodpoints = ', ngoodpoints, ' out of ', nfftot 396 !write(std_out,*)' total_weight=',total_weight 397 !write(std_out,*)' total_weight=',total_weight 398 !ENDDEBUG 399 400 !Output 401 if (prtcharge == 1) then 402 write(std_out,*) 403 write(std_out,*)' Hirshfeld analysis' 404 write(std_out,*)' Atom Zion Electron Charge Net charge ' 405 write(std_out,*) 406 do iatom=1,natom 407 write(std_out,'(i9,3es17.6)' )& 408 & iatom,zion(typat(iatom)),hcharge(iatom),hcharge(iatom)+zion(typat(iatom)) 409 end do 410 write(std_out,*) 411 write(std_out,'(a,3es17.6)')' Total',total_zion,total_charge,total_charge+total_zion 412 write(std_out,*) 413 end if 414 415 ABI_FREE(coord1) 416 ABI_FREE(local_den) 417 ABI_FREE(ncells) 418 ABI_FREE(step) 419 ABI_FREE(sum_den) 420 ABI_FREE(xcartcells) 421 ABI_FREE(xred) 422 ABI_FREE(yder2) 423 424 end subroutine dens_hirsh
m_dens/mag_penalty [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
mag_penalty
FUNCTION
This routine is called to compute the potential corresponding to constrained magnetic moments using the penalty function algorithm.
INPUTS
c_dft <type(constrained_dft_t)>=datastructure for the information related to constrained DFT ! magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size ! magcon_lambda=strength of the atomic spherical constraint ! natom=number of atoms ! nfftf=number of points in fine fft grid ! ngfftf=FFT grid dimensions ! nspden = number of spin densities (1 2 or 4) ! ntypat=number of types of atoms ! ratsm=smearing width for ratsph ! ratsph(ntypat)=radii for muffin tin spheres of each atom ! rprimd=lattice vectors (dimensioned) ! spinat(3,natom)=magnetic moments vectors, possible targets according to the value of constraint_kind ! typat(natom)=types of atoms mpi_enreg=mpi structure with communicator info rhor=density in real space xred=reduced atomic positions
OUTPUT
nv_constr_dft_r=the constrained potential
NOTES
based on html notes for the VASP implementation at http://cms.mpi.univie.ac.at/vasp/vasp/Constraining_direction_magnetic_moments.html
SOURCE
1202 subroutine mag_penalty(c_dft,mpi_enreg,rhor,nv_constr_dft_r,xred) 1203 1204 !Arguments ------------------------------------ 1205 !scalars 1206 type(constrained_dft_t),intent(in) :: c_dft 1207 real(dp),intent(out) :: nv_constr_dft_r(c_dft%nfftf,c_dft%nspden) 1208 type(MPI_type),intent(in) :: mpi_enreg 1209 !arrays 1210 real(dp),intent(in) :: rhor(c_dft%nfftf,c_dft%nspden) 1211 real(dp),intent(in) :: xred(3,c_dft%natom) 1212 1213 !Local variables------------------------------- 1214 !scalars 1215 integer :: iatom,magconon,natom,nfftf,nspden,ntypat,option 1216 integer :: cplex1=1 1217 real(dp):: cmm_x,cmm_y,cmm_z 1218 real(dp) :: intgden_proj,norm 1219 !arrays 1220 real(dp), allocatable :: coeffs_constr_dft(:,:) ! nspden,natom 1221 real(dp), allocatable :: intgden(:,:) ! nspden,natom 1222 real(dp) :: rhomag(2,c_dft%nspden),spinat_normed(3) 1223 1224 ! *********************************************************************************************** 1225 1226 magconon=c_dft%magconon 1227 natom=c_dft%natom 1228 nfftf=c_dft%nfftf 1229 nspden=c_dft%nspden 1230 ntypat=c_dft%ntypat 1231 1232 ABI_MALLOC(coeffs_constr_dft,(nspden,natom)) 1233 ABI_MALLOC(intgden,(nspden,natom)) 1234 1235 !We need the integrated magnetic moments and the smoothing function 1236 call calcdenmagsph(mpi_enreg,natom,nfftf,c_dft%ngfftf,nspden,ntypat,& 1237 & c_dft%ratsm,c_dft%ratsph,rhor,c_dft%rprimd,c_dft%typat,xred,1,cplex1,intgden=intgden,rhomag=rhomag) 1238 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,std_out,1,c_dft%ratsm,c_dft%ratsph,rhomag,c_dft%typat) 1239 1240 !Loop over atoms 1241 !------------------------------------------- 1242 do iatom=1,natom 1243 1244 norm = sqrt(sum(c_dft%spinat(:,iatom)**2)) 1245 spinat_normed(:) = zero 1246 if (norm > tol10) then 1247 spinat_normed(:) = c_dft%spinat(:,iatom) / norm 1248 else if (magconon == 1) then 1249 ! if spinat = 0 and we are imposing the direction only, skip this atom 1250 cycle 1251 end if 1252 1253 ! Calculate the x- and y-components of the square bracket term 1254 cmm_x = zero 1255 cmm_y = zero 1256 cmm_z = zero 1257 intgden_proj = zero 1258 if (nspden == 4) then 1259 if (magconon==1) then 1260 ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector 1261 ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector 1262 intgden_proj=spinat_normed(1)*intgden(2,iatom)+ & 1263 & spinat_normed(2)*intgden(3,iatom)+ & 1264 & spinat_normed(3)*intgden(4,iatom) 1265 1266 cmm_x=intgden(2,iatom) 1267 cmm_x=cmm_x-spinat_normed(1)*intgden_proj 1268 1269 cmm_y=intgden(3,iatom) 1270 cmm_y=cmm_y-spinat_normed(2)*intgden_proj 1271 1272 else if (magconon==2 .and. nspden == 4) then 1273 cmm_x=intgden(2,iatom)-c_dft%spinat(1,iatom) 1274 cmm_y=intgden(3,iatom)-c_dft%spinat(2,iatom) 1275 end if 1276 1277 ! Calculate the constraining potential for x- and y- components of the mag. mom. vector 1278 ! Eric Bousquet has derived the relationship between spin components and potential spin matrix elements: 1279 ! 1 = up up = +z 1280 ! 2 = down down = -z 1281 ! 3 = up down = +x 1282 ! 4 = down up = -y 1283 coeffs_constr_dft(3,iatom)= 2*c_dft%magcon_lambda*cmm_x 1284 coeffs_constr_dft(4,iatom)=-2*c_dft%magcon_lambda*cmm_y 1285 end if ! nspden 4 1286 1287 ! Calculate the z-component of the square bracket term 1288 if (magconon==1) then 1289 if (nspden == 4) then 1290 ! This apparently enforces the axis of magnetization, not its direction. 1291 ! m_z - spinat_z * <m | spinat> 1292 cmm_z = intgden(4,iatom) - spinat_normed(3)*intgden_proj 1293 else if (nspden == 2) then 1294 ! This apparently enforces the direction of the magnetization. So, the behaviour differs in nspden=4 and 2 cases . 1295 ! this will be just a sign +/- : are we in the same direction as spinat_z? 1296 ! need something more continuous??? To make sure the gradient pushes the state towards FM/AFM? 1297 cmm_z = -sign(one, (intgden(1,iatom)-intgden(2,iatom))*spinat_normed(3)) 1298 end if 1299 else if (magconon==2) then 1300 if (nspden == 4) then 1301 cmm_z=intgden(4,iatom)-c_dft%spinat(3,iatom) 1302 else if (nspden == 2) then 1303 ! this is up spins - down spins - requested moment ~ 0 1304 ! EB: note that intgden comes from calcdenmagsph, which, in nspden=2 case, returns 1305 ! intgden(1)=rho_up=n+m 1306 ! intgden(2)=rho_dn=n-m 1307 ! Then, is the following line be 1308 ! cmm_z=half*(intgden(1,iatom)-intgden(2,iatom)) - spinat(3,iatom) 1309 ! ?? 1310 cmm_z=intgden(1,iatom)-intgden(2,iatom) - c_dft%spinat(3,iatom) 1311 end if 1312 endif 1313 1314 ! Calculate the constraining potential for z-component of the mag. mom. vector 1315 coeffs_constr_dft(1,iatom)= 2*c_dft%magcon_lambda*cmm_z 1316 coeffs_constr_dft(2,iatom)=-2*c_dft%magcon_lambda*cmm_z 1317 1318 enddo ! iatom 1319 1320 !Now compute the potential in real space 1321 option=0 1322 call add_atomic_fcts(natom,nspden,c_dft%rprimd,mpi_enreg,nfftf,c_dft%ngfftf,ntypat,option,c_dft%ratsph, & 1323 c_dft%ratsm,c_dft%typat,coeffs_constr_dft,nv_constr_dft_r,xred) ! Ben change: add ratsm 1324 1325 ABI_FREE(coeffs_constr_dft) 1326 ABI_FREE(intgden) 1327 1328 end subroutine mag_penalty
m_dens/mag_penalty_e [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
mag_penalty_e
FUNCTION
Compute the energy corresponding to constrained magnetic moments.
INPUTS
magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size spinat=fixed magnetic moments vectors magcon_lambda=the size of the penalty terms
OUTPUT
Epen=penalty contribution to the total energy corresponding to the constrained potential Econstr=??? Eexp=???
SOURCE
1351 subroutine mag_penalty_e(magconon,magcon_lambda,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,spinat,typat,xred) 1352 1353 !Arguments ------------------------------------ 1354 !scalars 1355 integer,intent(in) :: natom,magconon,nspden,nfft,ntypat 1356 real(dp),intent(in) :: magcon_lambda,ratsm 1357 !arrays 1358 integer, intent(in) :: ngfft(18),typat(natom) 1359 real(dp),intent(in) :: spinat(3,natom), rprimd(3,3) 1360 real(dp),intent(in) :: ratsph(ntypat),rhor(nfft,nspden),xred(3,natom) 1361 type(MPI_type),intent(in) :: mpi_enreg 1362 1363 !Local variables------------------------------- 1364 !scalars 1365 integer :: iatom,ii 1366 integer :: cplex1=1 ! dummy argument for calcdenmagsph 1367 real(dp) :: intgden_proj, Epen,Econstr,lVp, norm 1368 !arrays 1369 real(dp) :: intmm(3), mag_1atom(3) 1370 real(dp), allocatable :: intgden(:,:) 1371 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),ucvol 1372 real(dp) :: rhomag(2,nspden),spinat_normed(3) 1373 character(len=500) :: msg 1374 1375 ! ********************************************************************* 1376 1377 !We need the metric because it is needed in calcdenmagsph.F90 1378 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1379 1380 ABI_MALLOC(intgden, (nspden,natom)) 1381 1382 !We need the integrated magnetic moments 1383 cplex1=1 1384 call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsm,ratsph,rhor,rprimd,typat,xred,& 1385 & 1,cplex1,intgden=intgden,rhomag=rhomag) 1386 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,std_out,1,ratsm,ratsph,rhomag,typat) 1387 1388 Epen=0 1389 Econstr=0 1390 lVp=0 1391 1392 !Loop over atoms 1393 !------------------------------------------- 1394 do iatom=1,natom 1395 1396 norm = sqrt(sum(spinat(:,iatom)**2)) 1397 spinat_normed(:) = zero 1398 if (norm > tol10) then 1399 spinat_normed(:) = spinat(:,iatom) / norm 1400 else if (magconon == 1) then 1401 ! if spinat = 0 and we are imposing the direction only, skip this atom 1402 cycle 1403 end if 1404 ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector 1405 ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector 1406 1407 ! for the collinear spin case, set up a fictitious 3D vector along z 1408 if (nspden == 4) then 1409 mag_1atom(1:3) = intgden(2:4,iatom) 1410 else if (nspden == 2) then 1411 mag_1atom = zero 1412 mag_1atom(3) = intgden(1,iatom)-intgden(2,iatom) 1413 end if 1414 1415 intgden_proj = zero 1416 intmm = zero 1417 ! Calculate the square bracket term 1418 if (magconon==1) then 1419 intgden_proj=spinat_normed(1)*mag_1atom(1)+ & 1420 & spinat_normed(2)*mag_1atom(2)+ & 1421 & spinat_normed(3)*mag_1atom(3) 1422 1423 do ii=1,3 1424 intmm(ii)=mag_1atom(ii)-spinat_normed(ii)*intgden_proj 1425 end do 1426 1427 ! Calculate the energy Epen corresponding to the constraining potential 1428 ! Econstr and lVp do not have a clear meaning (yet) 1429 Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3)) 1430 Econstr=Econstr-magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 1431 lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 1432 1433 else if (magconon==2) then 1434 do ii=1,3 1435 intmm(ii)=mag_1atom(ii)-spinat(ii,iatom) 1436 end do 1437 1438 ! Calculate the energy Epen corresponding to the constraining potential 1439 ! Epen = -Econstr - lVp 1440 ! Econstr = -M**2 + spinat**2 1441 ! lVp = +2 M \cdot spinat 1442 Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3)) 1443 Econstr=Econstr-magcon_lambda*(mag_1atom(1)*mag_1atom(1)+& 1444 & mag_1atom(2)*mag_1atom(2)+& 1445 & mag_1atom(3)*mag_1atom(3)) & 1446 & +magcon_lambda*(spinat(1,iatom)*spinat(1,iatom)+& 1447 & spinat(2,iatom)*spinat(2,iatom)+& 1448 & spinat(3,iatom)*spinat(3,iatom)) 1449 lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 1450 end if 1451 1452 write(msg, *) 'atom constraining magnetic field' 1453 call wrtout(std_out,msg,'COLL') 1454 write(msg, '(I3,A2,E12.5,A2,E12.5,A2,E12.5)') & 1455 iatom,' ',magcon_lambda*intmm(1),' ',magcon_lambda*intmm(2),' ',magcon_lambda*intmm(3) 1456 call wrtout(std_out,msg,'COLL') 1457 1458 ! End loop over atoms 1459 ! ------------------------------------------- 1460 end do 1461 1462 !Printing 1463 write(msg, '(A17,E10.3)' ) ' magcon_lambda = ',magcon_lambda 1464 call wrtout(std_out,msg,'COLL') 1465 write(msg, '(A17,E12.5)' ) ' Lagrange penalty = ',Epen 1466 call wrtout(std_out,msg,'COLL') 1467 write(msg, '(A17,E12.5)' ) ' E_constraint = ',Econstr 1468 call wrtout(std_out,msg,'COLL') 1469 write(msg, '(A17,E12.5)' ) ' lVp = ',lVp 1470 call wrtout(std_out,msg,'COLL') 1471 1472 ABI_FREE(intgden) 1473 1474 end subroutine mag_penalty_e
m_dens/prtdenmagsph [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
prtdenmagsph
FUNCTION
Print integral of total density inside spheres around atoms, and optionally integral of potential residual (also gradient of the constraint energy wrt constraint).
INPUTS
intgden(nspden, natom)=integrated rhor or potential residual, for each atom in a sphere of radius ratsph. if option <10, intgden is a density (+magnetization) Representation differs according to nspden : if nspden=1, total density if nspden=2, spin up, then spin down if nspden=4, total density, then mag_x, mag_y, mag_z if 20>option>=10, intgden is a potential residual (+spin magnetic field residual). Representation is: first, mean potential; then (if nspden>=2) B_z for nspden=2, B_x, B_y and B_z for nspden=4. if option>=20, intgden is a a gradient wrt target (=torque), also potential residual (+spin magnetic field residual) multiplied by f-1. natom=number of atoms in cell. nspden=number of spin-density components ntypat=number of atom types nunit=number of the unit for printing option = if not larger than 10, then a density is input , if between 10 and 19 then a potential residual is input, if beyond, a torque is input.. When 1, 11, 21, the default printing is on (to unit nunit), if -1, 2, 3, 4, special printing options, if 0 no printing. ratsm=smearing width for ratsph ratsph(ntypat)=radius of spheres around atoms rhomag(2,nspden)=integral of charge or magnetization over the whole cell (also taking into account a possible imaginary part for DFPT). typat(natom)=type of each atom ziontypat(ntypat)= --optional-- ionic charge of each atomic type
OUTPUT
Printing
SOURCE
1956 subroutine prtdenmagsph(cplex,intgden,natom,nspden,ntypat,nunit,option,ratsm,ratsph,rhomag,typat,ziontypat) 1957 1958 !Arguments --------------------------------------------- 1959 !scalars 1960 integer,intent(in) :: natom,nspden,ntypat,nunit 1961 real(dp),intent(in) :: ratsm 1962 integer ,intent(in) :: option 1963 integer, intent(in) :: cplex 1964 !arrays 1965 integer,intent(in) :: typat(natom) 1966 real(dp),intent(in) :: intgden(nspden,natom) 1967 real(dp),intent(in) :: ratsph(ntypat),rhomag(2,nspden) 1968 real(dp),intent(in),optional :: ziontypat(ntypat) 1969 !Local variables ------------------------------ 1970 !scalars 1971 integer :: iatom,ix 1972 real(dp) :: mag_coll , mag_x, mag_y, mag_z ! EB 1973 real(dp) :: mag_coll_im, mag_x_im, mag_y_im, mag_z_im ! SPr 1974 real(dp) :: rho_tot, rho_tot_im 1975 real(dp) :: sum_mag, sum_mag_x,sum_mag_y,sum_mag_z,sum_rho_up,sum_rho_dn,sum_rho_tot ! EB 1976 character(len=500) :: msg,msg1 1977 1978 ! ************************************************************************* 1979 1980 !DEBUG 1981 !write(ab_out,*)' prtdenmagsph : enter, rhomag(1,2)=',rhomag(1,2) 1982 !ENDDEBUG 1983 1984 if(nspden==2)then 1985 rho_tot=rhomag(1,1) ; mag_coll=rhomag(1,2) 1986 if(cplex==2)then 1987 rho_tot_im=rhomag(2,1) ; mag_coll_im=rhomag(2,2) 1988 endif 1989 else if(nspden==4)then 1990 rho_tot=rhomag(1,1) ; mag_x=rhomag(1,2) ; mag_y=rhomag(1,3) ; mag_z=rhomag(1,4) 1991 if(cplex==2)then 1992 rho_tot_im=rhomag(2,1) ; mag_x_im=rhomag(2,2) ; mag_y_im=rhomag(2,3) ; mag_z_im=rhomag(2,4) 1993 endif 1994 endif 1995 1996 if(option/=0)then 1997 1998 !Printing 1999 sum_mag=zero 2000 sum_mag_x=zero 2001 sum_mag_y=zero 2002 sum_mag_z=zero 2003 sum_rho_up=zero 2004 sum_rho_dn=zero 2005 sum_rho_tot=zero 2006 2007 if(option==1 .or. option==11 .or. option==21) then 2008 2009 if(nspden==1) then 2010 if(option== 1)msg1=' Integrated electronic density in atomic spheres:' 2011 if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres:' 2012 if(option==21)msg1=ch10//' Gradient with respect to target (=torque) :' 2013 write(msg, '(3a)' ) trim(msg1),ch10,' ------------------------------------------------' 2014 call wrtout(nunit,msg,'COLL') 2015 if(ratsm>tol8)then 2016 write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,'.' 2017 call wrtout(nunit,msg,'COLL') 2018 endif 2019 if(option== 1)then 2020 msg=' Atom Sphere_radius Integrated_density' 2021 if(present(ziontypat)) write(msg,'(a,a)')trim(msg),' Atomic charge' 2022 endif 2023 if(option==11)msg=' Atom Sphere_radius Integrated_potresid' 2024 if(option==21)msg=' Atom Sphere_radius Torque' 2025 call wrtout(nunit,msg,'COLL') 2026 do iatom=1,natom 2027 write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),intgden(1,iatom) 2028 if(option==21)then 2029 ! There is a change of sign to get the gradient wrt chrgat. 2030 write(msg, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),-intgden(1,iatom) 2031 endif 2032 !If option=1, print atomic charge 2033 if(option==1 .and. present(ziontypat))then 2034 write(msg, '(a,f20.8)' ) trim(msg),ziontypat(typat(iatom))-intgden(1,iatom) 2035 endif 2036 call wrtout(nunit,msg,'COLL') 2037 end do 2038 endif 2039 2040 if(nspden==2 .or. nspden==4) then 2041 2042 if(option== 1)msg1=' Integrated electronic and magnetization densities in atomic spheres:' 2043 if(option==11)msg1=ch10//' Integrated potential residual in atomic spheres (scalar + magnetic field):' 2044 if(option==21)msg1=ch10//' Gradient with respect to target (=torque) (scalar + magnetic field):' 2045 write(msg, '(3a)' ) trim(msg1),ch10,' ---------------------------------------------------------------------' 2046 call wrtout(nunit,msg,'COLL') 2047 2048 if(option== 1 .and. nspden==2) msg1='. Diff(up-dn)=approximate z local magnetic moment.' 2049 if(option== 1 .and. nspden==4) msg1='. mag(i)=approximate local magnetic moment.' 2050 if(option==11 .or. option==21) msg1='.' 2051 write(msg, '(a,f8.4,a)' ) ' Radius=ratsph(iatom), smearing ratsm=',ratsm,trim(msg1) 2052 call wrtout(nunit,msg,'COLL') 2053 2054 if(option==1)then 2055 if(nspden==2) msg=' Atom Radius up_density dn_density Total(up+dn) Diff(up-dn)' 2056 if(nspden==4) msg=' Atom Radius Total density mag(x) mag(y) mag(z) ' 2057 if(present(ziontypat))msg=trim(msg)//' Atomic charge' 2058 else if(option==11)then 2059 if(nspden==2) msg=' Atom Radius Potential B(z) up pot down pot' 2060 if(nspden==4) msg=' Atom Radius Potential B(x) B(y) B(z) ' 2061 else if(option==21)then 2062 if(nspden==2) msg=' Atom Radius grchrg T(z) up torque down torque' 2063 if(nspden==4) msg=' Atom Radius Torque T(x) T(y) T(z) ' 2064 endif 2065 call wrtout(nunit,msg,'COLL') 2066 2067 endif 2068 2069 if(nspden==2)then 2070 do iatom=1,natom 2071 if(option/=21)then 2072 write(msg,'(i5,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,ratsph(typat(iatom)),intgden(1,iatom),intgden(2,iatom) 2073 else 2074 write(msg,'(i5,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,ratsph(typat(iatom)),-intgden(1,iatom),intgden(2,iatom) 2075 endif 2076 write(msg,'(a,a,f12.6,a,f12.6)')trim(msg),' ',(intgden(1,iatom)+intgden(2,iatom)),' ',(intgden(1,iatom)-intgden(2,iatom)) 2077 if(option==1 .and. present(ziontypat))& 2078 & write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-(intgden(1,iatom)+intgden(2,iatom)) 2079 call wrtout(nunit,msg,'COLL') 2080 ! Compute the sum of the magnetization 2081 sum_mag=sum_mag+intgden(1,iatom)-intgden(2,iatom) 2082 sum_rho_up=sum_rho_up+intgden(1,iatom) 2083 sum_rho_dn=sum_rho_dn+intgden(2,iatom) 2084 sum_rho_tot=sum_rho_tot+intgden(1,iatom)+intgden(2,iatom) 2085 end do 2086 write(msg, '(a)') ' ---------------------------------------------------------------------' 2087 call wrtout(nunit,msg,'COLL') 2088 write(msg, '(a,2f13.6,a,f12.6,a,f12.6)') ' Sum: ', sum_rho_up,sum_rho_dn,' ',sum_rho_tot,' ',sum_mag 2089 call wrtout(nunit,msg,'COLL') 2090 2091 if(option==1)then 2092 write(msg, '(a,f14.6)') ' Total magnetization (from the atomic spheres): ', sum_mag 2093 call wrtout(nunit,msg,'COLL') 2094 write(msg, '(a,f14.6)') ' Total magnetization (exact up - dn): ', mag_coll 2095 call wrtout(nunit,msg,'COLL') 2096 endif 2097 2098 elseif(nspden==4) then 2099 2100 do iatom=1,natom 2101 if(option/=21)then 2102 write(msg, '(i5,f10.5,f16.6,a,3f12.6)' ) iatom,ratsph(typat(iatom)),intgden(1,iatom),' ',(intgden(ix,iatom),ix=2,4) 2103 else 2104 write(msg, '(i5,f10.5,f16.6,a,3f12.6)' ) iatom,ratsph(typat(iatom)),-intgden(1,iatom),' ',(intgden(ix,iatom),ix=2,4) 2105 endif 2106 if(option==1 .and. present(ziontypat))& 2107 & write(msg, '(a,f14.6)') trim(msg),ziontypat(typat(iatom))-intgden(1,iatom) 2108 call wrtout(nunit,msg,'COLL') 2109 ! Compute the sum of the magnetization in x, y and z directions 2110 sum_mag_x=sum_mag_x+intgden(2,iatom) 2111 sum_mag_y=sum_mag_y+intgden(3,iatom) 2112 sum_mag_z=sum_mag_z+intgden(4,iatom) 2113 end do 2114 write(msg, '(a)') ' ---------------------------------------------------------------------' 2115 call wrtout(nunit,msg,'COLL') 2116 2117 if(option==1)then 2118 write(msg, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (spheres) ', sum_mag_x,sum_mag_y,sum_mag_z 2119 call wrtout(nunit,msg,'COLL') 2120 write(msg, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (exact) ', mag_x,mag_y,mag_z 2121 call wrtout(nunit,msg,'COLL') 2122 endif 2123 2124 end if 2125 2126 elseif(option==-1) then 2127 2128 write(msg, '(2a)') ch10,' ------------------------------------------------------------------------' 2129 call wrtout(nunit,msg,'COLL') 2130 2131 if(nspden==1) then 2132 write(msg, '(4a)' ) & 2133 & ' Fermi level charge density n_f:',ch10,& 2134 & ' ------------------------------------------------------------------------',ch10 2135 else 2136 write(msg, '(4a)' ) & 2137 & ' Fermi level charge density n_f and magnetization m_f:',ch10,& 2138 & ' ------------------------------------------------------------------------',ch10 2139 end if 2140 call wrtout(nunit,msg,'COLL') 2141 2142 if(cplex==1) then 2143 write(msg, '(a,f13.8)') ' n_f = ',rho_tot 2144 else 2145 write(msg, '(a,f13.8,a,f13.8)') ' Re[n_f]= ', rho_tot," Im[n_f]= ",rho_tot_im 2146 end if 2147 call wrtout(nunit,msg,'COLL') 2148 if(nspden==2) then 2149 if(cplex==1) then 2150 write(msg, '(a,f13.8)') ' m_f = ', mag_coll 2151 else 2152 write(msg, '(a,f13.8,a,f13.8)') ' Re[m_f]= ', mag_coll," Im[m_f]= ",mag_coll_im 2153 end if 2154 call wrtout(nunit,msg,'COLL') 2155 elseif (nspden==4) then 2156 write(msg, '(a,f13.8)') ' mx_f = ',mag_x 2157 call wrtout(nunit,msg,'COLL') 2158 write(msg, '(a,f13.8)') ' my_f = ',mag_y 2159 call wrtout(nunit,msg,'COLL') 2160 write(msg, '(a,f13.8)') ' mz_f = ',mag_z 2161 call wrtout(nunit,msg,'COLL') 2162 end if 2163 2164 write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10 2165 call wrtout(nunit,msg,'COLL') 2166 2167 2168 else if (option==2 .or. option==3 .or. option==4) then 2169 ! Used in the DFPT case, option=idir+1 2170 2171 if(abs(rho_tot)<1.0d-10) then 2172 rho_tot=0 2173 end if 2174 2175 write(msg, '(2a)') ch10,' ------------------------------------------------------------------------' 2176 call wrtout(nunit,msg,'COLL') 2177 2178 if(nspden==1) then 2179 write(msg, '(4a)' ) & 2180 & ' Integral of the first order density n^(1):',ch10,& 2181 & ' ------------------------------------------------------------------------',ch10 2182 else 2183 write(msg, '(4a)' ) & 2184 & ' Integrals of the first order density n^(1) and magnetization m^(1):',ch10,& 2185 & ' ------------------------------------------------------------------------',ch10 2186 end if 2187 call wrtout(nunit,msg,'COLL') 2188 2189 if(cplex==1) then 2190 write(msg, '(a,e16.8)') ' n^(1) = ', rho_tot 2191 else 2192 write(msg, '(a,e16.8,a,e16.8)') ' Re[n^(1)] = ', rho_tot," Im[n^(1)] = ",rho_tot_im 2193 end if 2194 call wrtout(nunit,msg,'COLL') 2195 2196 if(nspden==2) then 2197 2198 if(cplex==1) then 2199 write(msg, '(a,e16.8)') ' m^(1) = ', mag_coll 2200 else 2201 write(msg, '(a,e16.8,a,e16.8)') ' Re[m^(1)] = ', mag_coll," Im[m^(1)] = ",mag_coll_im 2202 end if 2203 call wrtout(nunit,msg,'COLL') 2204 2205 elseif (nspden==4) then 2206 if(cplex==1) then 2207 write(msg, '(a,e16.8)') ' mx^(1) = ', mag_x 2208 call wrtout(nunit,msg,'COLL') 2209 write(msg, '(a,e16.8)') ' my^(1) = ', mag_y 2210 call wrtout(nunit,msg,'COLL') 2211 write(msg, '(a,e16.8)') ' mz^(1) = ', mag_z 2212 call wrtout(nunit,msg,'COLL') 2213 else 2214 write(msg, '(a,e16.8,a,e16.8)') ' Re[mx^(1)]= ', mag_x, " Im[mx^(1)]= ", mag_x_im 2215 call wrtout(nunit,msg,'COLL') 2216 write(msg, '(a,e16.8,a,e16.8)') ' Re[my^(1)]= ', mag_y, " Im[my^(1)]= ", mag_y_im 2217 call wrtout(nunit,msg,'COLL') 2218 write(msg, '(a,e16.8,a,e16.8)') ' Re[mz^(1)]= ', mag_z, " Im[mz^(1)]= ", mag_z_im 2219 call wrtout(nunit,msg,'COLL') 2220 end if 2221 end if 2222 2223 write(msg, '(3a)') ch10,' ------------------------------------------------------------------------',ch10 2224 call wrtout(nunit,msg,'COLL') 2225 2226 end if 2227 2228 end if ! option/=0 2229 2230 end subroutine prtdenmagsph
m_dens/radsmear [ Functions ]
[ Top ] [ m_dens ] [ Functions ]
NAME
radsmear
FUNCTION
As a function of the argument xarg (a positive number), return a function fsm that is zero beyond some cut-off value xcut, one for xarg smaller than xcut-xsmear, and interpolates smoothly between one and zero in the region from xcut-xsmear to xcut. Also returns the derivative of this function, called dfsm. The function fsm is twice differentiable at xcut (the first derivative is continuous, not the second), and three times differentiable at xcut-xsmear (the second derivative is continuous, not the third).
INPUTS
xarg=argument of the function (should be positive) xcut=largest value for which the function is non-zero xsmear=defined the smearing region, between xcut-xsmear and xcut
OUTPUT
fsm=value of the function dfsm=derivative of the function with respect to xarg (zero, except in the smearing region).
SOURCE
2257 subroutine radsmear(dfsm,fsm,xarg,xcut,xsmear) 2258 2259 !Arguments ------------------------------------ 2260 !scalars 2261 real(dp), intent(out) :: dfsm,fsm 2262 real(dp), intent(in) :: xarg, xcut, xsmear 2263 2264 !Local variables ------------------------------ 2265 !scalars 2266 real(dp) :: xsmearinv,xx 2267 2268 !****************************************************************** 2269 2270 fsm = zero 2271 dfsm=zero 2272 if (xarg < xcut - xsmear - tol12) then 2273 fsm = one 2274 else if (xarg < xcut - tol12) then 2275 xsmearinv=one/xsmear 2276 xx = (xcut - xarg) * xsmearinv 2277 fsm = xx**2*(3+xx*(1+xx*(-6+3*xx))) 2278 dfsm = -(xx*(6+xx*(3+xx*(-24+15*xx))))*xsmearinv 2279 end if 2280 2281 end subroutine radsmear