TABLE OF CONTENTS
- ABINIT/m_bader
- defs_aimprom/aim_shutdown
- m_bader/addout
- m_bader/adini
- m_bader/aim_dataset_type
- m_bader/aim_follow
- m_bader/bcp_type
- m_bader/bschg1
- m_bader/bschg2
- m_bader/consist
- m_bader/cpdrv
- m_bader/critic
- m_bader/critics
- m_bader/defad
- m_bader/drvaim
- m_bader/graph
- m_bader/initaim
- m_bader/inpar
- m_bader/inspln
- m_bader/integrho
- m_bader/integvol
- m_bader/mprod
- m_bader/onestep
- m_bader/ordr
- m_bader/plint
- m_bader/rsurf
- m_bader/surf
- m_bader/vec_prod
- m_bader/vgh_rho
- m_bader/vnorm
ABINIT/m_bader [ Modules ]
NAME
m_bader
FUNCTION
Procedures used by AIM code.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (PCasek,FF,XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_bader 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_xmpi 28 use m_sort 29 use m_hdr 30 use m_splines 31 #ifdef HAVE_NETCDF 32 use netcdf 33 #endif 34 35 use m_time, only : timein 36 use m_geometry, only : metric 37 use m_parser, only : inread 38 use m_numeric_tools, only : coeffs_gausslegint 39 use m_hide_lapack, only : jacobi, lubksb, ludcmp 40 41 implicit none 42 43 !private 44 public
defs_aimprom/aim_shutdown [ Functions ]
NAME
aim_shutdown
FUNCTION
Free memory allocated in the module. Close units. Mainly used to pass the abirules
SOURCE
221 subroutine aim_shutdown() 222 223 !Local variables------------------------------- 224 integer :: ii 225 logical :: is_open 226 integer :: all_units(12) 227 228 ! ********************************************************************* 229 230 !if (allocated(typat)) then 231 ! ABI_FREE(typat) 232 !end if 233 !if (allocated(corlim)) then 234 ! ABI_FREE(corlim) 235 !end if 236 !if (allocated(xred)) then 237 ! ABI_FREE(xred) 238 !end if 239 !if (allocated(xatm)) then 240 ! ABI_FREE(rminl) 241 !end if 242 243 all_units(:) = [unt0,unto,unt,untc,unts,untd,untl,untg,unta,untad,untp,untout] 244 do ii=1,size(all_units) 245 inquire(unit=all_units(ii), opened=is_open) 246 if (is_open) close(all_units(ii)) 247 end do 248 249 end subroutine aim_shutdown
m_bader/addout [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
addout
FUNCTION
Output density and laplacian (see input variables denout and lapout)
INPUTS
aim_dtset=the structured entity containing all input variables also, uses the variables saved in the module "defs_aimprom"
OUTPUT
(print) WARNING This file does not follow the ABINIT coding rules (yet) : the use of a module to transfer data should be avoided
SOURCE
760 subroutine addout(aim_dtset) 761 762 !Arguments ------------------------------------ 763 !scalars 764 type(aim_dataset_type),intent(in) :: aim_dtset 765 766 !Local variables ------------------------------ 767 !scalars 768 integer :: cod,dims,iat,ii,ipos,jj,nn,tgrd 769 real(dp) :: alfa,rho,rr,xx,yy 770 !arrays 771 real(dp) :: grho(3),hrho(3,3),orig(3),vv(3) 772 real(dp),allocatable :: dfld(:),lfld(:),nr(:),stp(:),uu(:,:) 773 774 !************************************************************************ 775 orig(:)=aim_dtset%vpts(:,1) 776 if (aim_dtset%denout > 0) then 777 dims=aim_dtset%denout 778 elseif (aim_dtset%lapout > 0) then 779 dims=aim_dtset%lapout 780 end if 781 782 select case (aim_dtset%dltyp) 783 case (1) 784 cod=1 785 case (2) 786 cod=2 787 case default 788 cod=0 789 end select 790 791 ABI_MALLOC(uu,(3,dims)) 792 ABI_MALLOC(nr,(dims)) 793 ABI_MALLOC(stp,(dims)) 794 795 write(std_out,*) 'grid:', aim_dtset%ngrid(1:dims) 796 write(std_out,*) 'kod :', cod 797 tgrd=1 798 do ii=1,dims 799 tgrd=tgrd*aim_dtset%ngrid(ii) 800 uu(:,ii)=aim_dtset%vpts(:,ii+1)-aim_dtset%vpts(:,1) 801 nr(ii)=vnorm(uu(:,ii),0) 802 stp(ii)=nr(ii)/(aim_dtset%ngrid(ii)-1) 803 uu(:,ii)=uu(:,ii)/nr(ii) 804 end do 805 write(std_out,*) 'tgrd :', tgrd 806 do ii=1,dims 807 write(std_out,*) 'uu :', uu(1:3,ii) 808 end do 809 810 if (aim_dtset%denout > 0) then 811 ABI_MALLOC(dfld,(tgrd+1)) 812 dfld(:)=0._dp 813 end if 814 if (aim_dtset%lapout > 0) then 815 ABI_MALLOC(lfld,(tgrd+1)) 816 end if 817 818 select case (dims) 819 case (1) 820 nn=0 821 do ii=0,aim_dtset%ngrid(1)-1 822 nn=nn+1 823 vv(:)=orig(:)+ii*stp(1)*uu(:,1) 824 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,cod) 825 if (aim_dtset%denout > 0) dfld(nn)=rho 826 if (aim_dtset%lapout > 0) lfld(nn)=hrho(1,1)+hrho(2,2)+hrho(3,3) 827 end do 828 if (aim_dtset%denout==1) then 829 do ii=0,aim_dtset%ngrid(1)-1 830 xx=ii*stp(1) 831 write(untd,'(2E16.8)') xx, dfld(ii+1) 832 end do 833 end if 834 if (aim_dtset%lapout==1) then 835 do ii=0,aim_dtset%ngrid(1)-1 836 xx=ii*stp(1) 837 write(untl,'(2E16.8)') xx, lfld(ii+1) 838 end do 839 end if 840 case (2) 841 nn=0 842 alfa=dot_product(uu(:,1),uu(:,2)) 843 alfa=acos(alfa) 844 do ii=0,aim_dtset%ngrid(2)-1 845 do jj=0,aim_dtset%ngrid(1)-1 846 nn=nn+1 847 vv(:)=orig(:)+jj*uu(:,2)*stp(2)+ii*stp(1)*uu(:,1) 848 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,cod) 849 if (aim_dtset%denout > 0) dfld(nn)=rho 850 if (aim_dtset%lapout > 0) lfld(nn)=hrho(1,1)+hrho(2,2)+hrho(3,3) 851 end do 852 end do 853 write(std_out,*) 'generace hotova', nn 854 nn=0 855 if (aim_dtset%denout==2) then 856 do ii=0,aim_dtset%ngrid(2)-1 857 do jj=0,aim_dtset%ngrid(1)-1 858 nn=nn+1 859 xx=jj*stp(1)+cos(alfa)*ii*stp(2) 860 yy=sin(alfa)*ii*stp(2) 861 write(untd,'(3E16.8)') xx, yy, dfld(nn) 862 end do 863 write(untd,*) ' ' 864 end do 865 end if 866 nn=0 867 if (aim_dtset%lapout==2) then 868 write(std_out,*) 'lezes sem?' 869 do ii=0,aim_dtset%ngrid(2)-1 870 do jj=0,aim_dtset%ngrid(1)-1 871 nn=nn+1 872 xx=jj*stp(1)+cos(alfa)*ii*stp(2) 873 yy=sin(alfa)*ii*stp(2) 874 write(untl,'(3E16.8)') xx, yy, lfld(nn) 875 end do 876 write(untl,*) ' ' 877 end do 878 end if 879 end select 880 ABI_FREE(uu) 881 ABI_FREE(stp) 882 ABI_FREE(nr) 883 if(aim_dtset%denout>0) then 884 ABI_FREE(dfld) 885 end if 886 if(aim_dtset%lapout>0) then 887 ABI_FREE(lfld) 888 end if 889 890 end subroutine addout
m_bader/adini [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
adini
FUNCTION
Analysis of the input string "inpstr" (the content of input file) and setting of the corresponding input variables
INPUTS
inpstr=character string containing the input data, to be treated lenstr=actual length of the string contained in inpstr
OUTPUT
aim_dtset=the structured entity containing all input variables
SOURCE
269 subroutine adini(aim_dtset,inpstr,lenstr) 270 271 !Arguments ------------------------------------ 272 !scalars 273 integer,intent(in) :: lenstr 274 character(len=*),intent(in) :: inpstr 275 !no_abirules 276 type(aim_dataset_type), intent(inout) :: aim_dtset !vz_i 277 278 !Local variables ------------------------------ 279 !scalars 280 integer :: errcod,ii,inxh,ipos,jj,lenc,ll,outi,tstngr=0,tstvpt=0 !vz_z 281 real(dp) :: outr 282 logical :: nbtst,try 283 character(len=20) :: cmot 284 285 ! ********************************************************************* 286 287 if (iachar(inpstr(1:1)) < 32) then 288 ipos=2 289 else 290 ipos=1 291 end if 292 293 write(std_out,*) 'ECHO of the INPUT' 294 write(std_out,*) '************************' 295 write(untout,*) 'ECHO of the INPUT' 296 write(untout,*) '************************' 297 298 mread: do ii=1,lenstr 299 try=.false. 300 nbtst=.true. 301 inxh=index(inpstr(ipos:lenstr),' ') 302 if ((ipos >= lenstr)) exit 303 if ((inxh==2).or.(inxh==1)) then 304 ipos=ipos+inxh 305 cycle 306 end if 307 lenc=inxh-1 308 cmot(1:lenc)=inpstr(ipos:ipos+inxh-2) 309 ipos=ipos+inxh 310 ! write(std_out,*) cmot(1:lenc), lenc 311 312 select case (cmot(1:lenc)) 313 314 ! DRIVER SPECIFICATIONS 315 316 case ('SURF') 317 inxh=index(inpstr(ipos:lenstr),' ') 318 if ((inxh /= 2).and.(inpstr(ipos:ipos)/='-')) then 319 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 320 ABI_ERROR("Aborting now") 321 end if 322 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 323 aim_dtset%isurf=outi 324 write(std_out,*) cmot(1:lenc),' ', aim_dtset%isurf 325 write(untout,*) cmot(1:lenc),' ', aim_dtset%isurf 326 ipos=ipos+inxh 327 328 case ('CRIT') 329 inxh=index(inpstr(ipos:lenstr),' ') 330 if ((inxh /= 2).and.(inpstr(ipos:ipos)/='-')) then 331 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 332 ABI_ERROR("Aborting now") 333 end if 334 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 335 aim_dtset%crit=outi 336 write(std_out,*) cmot(1:lenc),' ', aim_dtset%crit 337 write(untout,*) cmot(1:lenc),' ', aim_dtset%crit 338 ipos=ipos+inxh 339 340 case ('RSURF') 341 inxh=index(inpstr(ipos:lenstr),' ') 342 if (inxh /= 2) then 343 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 344 ABI_ERROR("Aborting now") 345 end if 346 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 347 aim_dtset%irsur=outi 348 write(std_out,*) cmot(1:lenc),' ', aim_dtset%irsur 349 write(untout,*) cmot(1:lenc),' ', aim_dtset%irsur 350 ipos=ipos+inxh 351 352 case ('FOLLOW') 353 inxh=index(inpstr(ipos:lenstr),' ') 354 if (inxh /= 2) then 355 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 356 ABI_ERROR("Aborting now") 357 end if 358 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 359 aim_dtset%foll=outi 360 write(std_out,*) cmot(1:lenc),' ', aim_dtset%foll 361 write(untout,*) cmot(1:lenc),' ', aim_dtset%foll 362 ipos=ipos+inxh 363 364 case ('IRHO') 365 inxh=index(inpstr(ipos:lenstr),' ') 366 if (inxh /= 2) then 367 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 368 ABI_ERROR("Aborting now") 369 end if 370 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 371 aim_dtset%irho=outi 372 write(std_out,*) cmot(1:lenc),' ', aim_dtset%irho 373 write(untout,*) cmot(1:lenc),' ', aim_dtset%irho 374 ipos=ipos+inxh 375 376 case ('PLDEN') 377 inxh=index(inpstr(ipos:lenstr),' ') 378 if (inxh /= 2) then 379 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 380 ABI_ERROR("Aborting now") 381 end if 382 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 383 aim_dtset%plden=outi 384 write(std_out,*) cmot(1:lenc),' ', aim_dtset%plden 385 write(untout,*) cmot(1:lenc),' ', aim_dtset%plden 386 ipos=ipos+inxh 387 388 389 case ('IVOL') 390 inxh=index(inpstr(ipos:lenstr),' ') 391 if (inxh /= 2) then 392 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 393 ABI_ERROR("Aborting now") 394 end if 395 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 396 aim_dtset%ivol=outi 397 write(std_out,*) cmot(1:lenc),' ', aim_dtset%ivol 398 write(untout,*) cmot(1:lenc),' ', aim_dtset%ivol 399 ipos=ipos+inxh 400 401 case ('DENOUT') 402 inxh=index(inpstr(ipos:lenstr),' ') 403 if (inxh /= 2) then 404 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 405 ABI_ERROR("Aborting now") 406 end if 407 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 408 aim_dtset%denout=outi 409 if ((aim_dtset%denout < -1).or.(aim_dtset%denout>3)) then 410 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 411 ABI_ERROR("Aborting now") 412 end if 413 write(std_out,*) cmot(1:lenc),' ', aim_dtset%denout 414 write(untout,*) cmot(1:lenc),' ', aim_dtset%denout 415 ipos=ipos+inxh 416 417 case ('LAPOUT') 418 inxh=index(inpstr(ipos:lenstr),' ') 419 if (inxh /= 2) then 420 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 421 ABI_ERROR("Aborting now") 422 end if 423 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 424 aim_dtset%lapout=outi 425 if ((aim_dtset%lapout < -1).or.(aim_dtset%lapout>3)) then 426 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 427 ABI_ERROR("Aborting now") 428 end if 429 write(std_out,*) cmot(1:lenc),' ', aim_dtset%lapout 430 write(untout,*) cmot(1:lenc),' ', aim_dtset%lapout 431 ipos=ipos+inxh 432 433 case ('DLTYP') 434 inxh=index(inpstr(ipos:lenstr),' ') 435 if (inxh /= 2) then 436 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 437 ABI_ERROR("Aborting now") 438 end if 439 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 440 aim_dtset%dltyp=outi 441 write(std_out,*) cmot(1:lenc),' ', aim_dtset%dltyp 442 write(untout,*) cmot(1:lenc),' ', aim_dtset%dltyp 443 ipos=ipos+inxh 444 445 case ('GPSURF') 446 inxh=index(inpstr(ipos:lenstr),' ') 447 if (inxh /= 2) then 448 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 449 ABI_ERROR("Aborting now") 450 end if 451 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 452 aim_dtset%gpsurf=outi 453 write(std_out,*) cmot(1:lenc),' ', aim_dtset%gpsurf 454 write(untout,*) cmot(1:lenc),' ', aim_dtset%gpsurf 455 ipos=ipos+inxh 456 457 458 ! END OF THE DRIVER SPECIFICATIONS 459 460 case ('ATOM') 461 inxh=index(inpstr(ipos:lenstr),' ') 462 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 463 aim_dtset%batom=outi 464 write(std_out,*) cmot(1:lenc),' ', aim_dtset%batom 465 write(untout,*) cmot(1:lenc),' ', aim_dtset%batom 466 ipos=ipos+inxh 467 468 case ('NSA') 469 inxh=index(inpstr(ipos:lenstr),' ') 470 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 471 aim_dtset%nsa=outi 472 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsa 473 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsa 474 ipos=ipos+inxh 475 476 case ('NSB') 477 inxh=index(inpstr(ipos:lenstr),' ') 478 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 479 aim_dtset%nsb=outi 480 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsb 481 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsb 482 ipos=ipos+inxh 483 484 case ('NSC') 485 inxh=index(inpstr(ipos:lenstr),' ') 486 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 487 aim_dtset%nsc=outi 488 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsc 489 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsc 490 ipos=ipos+inxh 491 492 case ('INPT') 493 inxh=index(inpstr(ipos:lenstr),' ') 494 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 495 aim_dtset%npt=outi 496 write(std_out,*) cmot(1:lenc),' ', aim_dtset%npt 497 write(untout,*) cmot(1:lenc),' ', aim_dtset%npt 498 ipos=ipos+inxh 499 500 case ('NTHETA') 501 inxh=index(inpstr(ipos:lenstr),' ') 502 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 503 aim_dtset%nth=outi 504 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nth 505 write(untout,*) cmot(1:lenc),' ', aim_dtset%nth 506 ipos=ipos+inxh 507 508 case ('NPHI') 509 inxh=index(inpstr(ipos:lenstr),' ') 510 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 511 aim_dtset%nph=outi 512 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nph 513 write(untout,*) cmot(1:lenc),' ', aim_dtset%nph 514 ipos=ipos+inxh 515 516 case ('THETAMIN') 517 inxh=index(inpstr(ipos:lenstr),' ') 518 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 519 aim_dtset%themin=outr 520 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themin 521 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themin 522 ipos=ipos+inxh 523 524 case ('THETAMAX') 525 inxh=index(inpstr(ipos:lenstr),' ') 526 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 527 aim_dtset%themax=outr 528 write(std_out, '(1x,a,a,es17.10)' ) cmot(1:lenc),' ', aim_dtset%themax 529 write(untout,'(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themax 530 ipos=ipos+inxh 531 532 case ('PHIMIN') 533 inxh=index(inpstr(ipos:lenstr),' ') 534 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 535 aim_dtset%phimin=outr 536 write(std_out, '(1x,a,a,es17.10)' ) cmot(1:lenc),' ', aim_dtset%phimin 537 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimin 538 ipos=ipos+inxh 539 540 case ('PHIMAX') 541 inxh=index(inpstr(ipos:lenstr),' ') 542 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 543 aim_dtset%phimax=outr 544 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimax 545 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimax 546 ipos=ipos+inxh 547 548 case ('ATRAD') 549 inxh=index(inpstr(ipos:lenstr),' ') 550 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 551 aim_dtset%atrad=outr 552 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%atrad 553 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%atrad 554 ipos=ipos+inxh 555 556 case ('RADSTP') 557 inxh=index(inpstr(ipos:lenstr),' ') 558 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 559 aim_dtset%dr0=outr 560 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dr0 561 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dr0 562 ipos=ipos+inxh 563 564 case ('FOLSTP') 565 inxh=index(inpstr(ipos:lenstr),' ') 566 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 567 aim_dtset%folstp=outr 568 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%folstp 569 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%folstp 570 ipos=ipos+inxh 571 572 case ('RATMIN') 573 inxh=index(inpstr(ipos:lenstr),' ') 574 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 575 aim_dtset%rmin=outr 576 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%rmin 577 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%rmin 578 ipos=ipos+inxh 579 580 case ('COFF1') 581 inxh=index(inpstr(ipos:lenstr),' ') 582 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 583 aim_dtset%coff1=outr 584 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff1 585 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff1 586 ipos=ipos+inxh 587 588 case ('COFF2') 589 inxh=index(inpstr(ipos:lenstr),' ') 590 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 591 aim_dtset%coff2=outr 592 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff2 593 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff2 594 ipos=ipos+inxh 595 596 case ('DPCLIM') 597 inxh=index(inpstr(ipos:lenstr),' ') 598 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 599 aim_dtset%dpclim=outr 600 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dpclim 601 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dpclim 602 ipos=ipos+inxh 603 604 case ('LGRAD') 605 inxh=index(inpstr(ipos:lenstr),' ') 606 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 607 aim_dtset%lgrad=outr 608 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad 609 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad 610 ipos=ipos+inxh 611 612 case ('LGRAD2') 613 inxh=index(inpstr(ipos:lenstr),' ') 614 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 615 aim_dtset%lgrad2=outr 616 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad2 617 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad2 618 ipos=ipos+inxh 619 620 case ('LSTEP') 621 inxh=index(inpstr(ipos:lenstr),' ') 622 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 623 aim_dtset%lstep=outr 624 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep 625 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep 626 ipos=ipos+inxh 627 628 case ('LSTEP2') 629 inxh=index(inpstr(ipos:lenstr),' ') 630 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 631 aim_dtset%lstep2=outr 632 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep2 633 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep2 634 ipos=ipos+inxh 635 636 case ('RSURDIR') 637 inxh=index(inpstr(ipos:lenstr),' ') 638 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 639 aim_dtset%th0=outr 640 ipos=ipos+inxh 641 inxh=index(inpstr(ipos:lenstr),' ') 642 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 643 aim_dtset%phi0=outr 644 ipos=ipos+inxh 645 write(std_out, '(1x,a,a,2es17.10)') cmot(1:lenc),' ', aim_dtset%th0, aim_dtset%phi0 646 write(untout, '(1x,a,a,2es17.10)') cmot(1:lenc),' ', aim_dtset%th0, aim_dtset%phi0 647 648 case ('FOLDEP') 649 do jj=1,3 650 inxh=index(inpstr(ipos:lenstr),' ') 651 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 652 aim_dtset%foldep(jj)=outr 653 ipos=ipos+inxh 654 end do 655 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%foldep 656 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%foldep 657 658 case ('SCAL') 659 do jj=1,3 660 inxh=index(inpstr(ipos:lenstr),' ') 661 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 662 aim_dtset%scal(jj)=outr 663 ipos=ipos+inxh 664 end do 665 write(std_out,*) cmot(1:lenc),' ', aim_dtset%scal 666 write(untout,*) cmot(1:lenc),' ', aim_dtset%scal 667 668 case ('NGRID') 669 try=.true. 670 do jj=1,3 671 inxh=index(inpstr(ipos:lenstr),' ') 672 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 673 aim_dtset%ngrid(jj)=outi 674 if (.not.nbtst) then 675 tstngr=jj-1 676 cycle mread 677 end if 678 if (inxh==0) then 679 tstvpt=jj 680 exit mread 681 end if 682 ipos=ipos+inxh 683 if (ipos==lenstr-1) then 684 tstngr=jj 685 exit mread 686 end if 687 end do 688 ! Why no echo ?? XG 030218 689 tstngr=3 690 691 case ('VPTS') 692 do jj=1,4 693 do ll=1,3 694 try=.true. 695 inxh=index(inpstr(ipos:lenstr),' ') 696 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 697 aim_dtset%vpts(ll,jj)=outr 698 if (.not.nbtst) then 699 tstvpt=jj-1 700 cycle mread 701 end if 702 ipos=ipos+inxh 703 if (ipos>=lenstr) then 704 tstvpt=jj 705 exit mread 706 end if 707 end do 708 end do 709 ! Why no echo ?? XG 030218 710 tstvpt=4 711 712 case ('MAXATD') 713 inxh=index(inpstr(ipos:lenstr),' ') 714 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 715 aim_dtset%maxatd=outr 716 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxatd 717 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxatd 718 ipos=ipos+inxh 719 720 case ('MAXCPD') 721 inxh=index(inpstr(ipos:lenstr),' ') 722 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 723 aim_dtset%maxcpd=outr 724 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxcpd 725 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxcpd 726 ipos=ipos+inxh 727 728 case default 729 write(std_out,*) 'ERROR Bad key word ! ',cmot(1:lenc) 730 end select 731 end do mread 732 733 write(std_out,*) '************************' 734 735 call consist(aim_dtset,tstngr,tstvpt) 736 737 end subroutine adini
m_bader/aim_dataset_type [ Types ]
NAME
aim_dataset_type
FUNCTION
The aim_dataset_type structured datatype gathers all the input variables for the aim code
SOURCE
57 type aim_dataset_type 58 59 ! Since all these input variables are described in the aim_help.html 60 ! file, they are not described in length here ... 61 62 ! Integer 63 integer :: crit 64 integer :: denout 65 integer :: dltyp 66 integer :: gpsurf 67 integer :: irho 68 integer :: ivol 69 integer :: lapout 70 integer :: nsa 71 integer :: nsb 72 integer :: nsc 73 74 integer :: batom ! Warning : corresponds to the input variable atom 75 integer :: foll ! Warning : corresponds to the input variable follow 76 integer :: isurf ! Warning : corresponds to the input variable surf 77 integer :: irsur ! Warning : corresponds to the input variable rsurf 78 integer :: nph ! Warning : corresponds to the input variable nphi 79 integer :: npt ! Warning : corresponds to the input variable inpt 80 integer :: nth ! Warning : corresponds to the input variable ntheta 81 integer :: plden ! Warning : not documented in help file ?! 82 83 integer :: ngrid(3) 84 85 ! Real 86 real(dp) :: atrad 87 real(dp) :: coff1 88 real(dp) :: coff2 89 real(dp) :: dpclim 90 real(dp) :: folstp 91 real(dp) :: lgrad 92 real(dp) :: lgrad2 93 real(dp) :: lstep 94 real(dp) :: lstep2 95 real(dp) :: maxatd 96 real(dp) :: maxcpd 97 real(dp) :: phimax 98 real(dp) :: phimin 99 100 real(dp) :: dr0 ! Warning : correspond to the input variable radstp 101 real(dp) :: phi0 ! Warning : correspond to the input variable rsurdir(2) 102 real(dp) :: rmin ! Warning : correspond to the input variable ratmin 103 real(dp) :: th0 ! Warning : correspond to the input variable rsurdir(1) 104 real(dp) :: themax ! Warning : correspond to the input variable thetamax 105 real(dp) :: themin ! Warning : correspond to the input variable thetamin 106 107 real(dp) :: foldep(3) 108 real(dp) :: scal(3) 109 real(dp) :: vpts(3,4) 110 111 end type aim_dataset_type
m_bader/aim_follow [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
aim_follow
FUNCTION
This routine follows the gradient line starting from the point vv. It stop when it arrives to the atom (nearer than rminl(iat)) or - if srch=true - also if it arrives under the already known part of Bader surface
INPUTS
aim_dtset= the structured entity containing all input variables iatinit,iposinit= indexes of initial atom npmax= maximum number of division in each step
OUTPUT
iat,ipos= index of final atom nstep= returns the number of step needed
SIDE EFFECTS
srch= (true/false) check if the line is outside or
inside the atomic surface.
vv(3)= initial point in orthogonal coordinates
SOURCE
919 subroutine aim_follow(aim_dtset,vv,npmax,srch,iatinit,iposinit,iat,ipos,nstep) 920 921 !Arguments ------------------------------------ 922 !scalars 923 integer,intent(in) :: iatinit,iposinit,npmax 924 integer,intent(out) :: iat,ipos,nstep 925 logical,intent(inout) :: srch 926 type(aim_dataset_type),intent(in) :: aim_dtset 927 !arrays 928 real(dp),intent(inout) :: vv(3) 929 930 !Local variables ------------------------------ 931 !scalars 932 integer :: i1,i2,i3,ii,iph,ires,ith,jj,kk,nit,np,nph,nsi,nth 933 real(dp) :: deltar,dg,dist,dph,dth,fac2,facf,h0old,hh,hold,rho,rr,rsmed 934 real(dp) :: t1,t2,t3,vcth,vph,vth,wall,xy,xyz 935 logical :: fin,ldebold,srchold,stemp,stemp2 936 character(len=50) :: formpc 937 character(len=500) :: msg 938 !arrays 939 real(dp) :: ev(3),grho(3),hrho(3,3),pom(3),vold(3),vt(3),vt1(3) 940 real(dp) :: zz(3,3) 941 942 !************************************************************************ 943 formpc='(":CP",2I5,3F12.8,3E12.4,I4,2E12.4)' 944 945 946 fin=.false. 947 948 srchold=srch 949 ldebold=ldeb 950 h0old=h0 951 952 nth=aim_dtset%nth 953 nph=aim_dtset%nph 954 955 if (slc==0) then 956 rminl(:)=aim_dtset%rmin 957 end if 958 959 if (deb) then 960 ldeb=.true. 961 end if 962 963 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 964 965 !Initial tests 966 967 if (iat/=0) then 968 if (rr<rminl(iat)) then 969 fin=.true. 970 write(std_out,*) 'rr < rmin iat=',iat,' ipos=',ipos 971 elseif (rho<aim_rhomin) then 972 fin=.true. 973 write(std_out,*) 'CHARGE LT rhomin ',rho,' < ',aim_rhomin 974 if (rho<zero) then 975 ABI_ERROR('RHO < 0 !!!') 976 end if 977 end if 978 end if 979 980 facf=aim_fac0 981 hh=aim_hmax 982 983 call timein(t1,wall) 984 nstep=0 985 nsi=0 986 987 !the principal cycle 988 989 madw : do while(.not.fin) 990 hold=hh 991 992 dg=vnorm(grho,0) 993 if (ldeb.or.deb) write(std_out,*) 'dg= ',dg 994 995 ! the time test 996 997 call timein(t3,wall) 998 t2=t3-t1 999 if (t2>300.0) then 1000 write(std_out,*) 'TIME EXCEEDED 5 min IN FOLLOW' 1001 write(std_out,*) 'h0 =',h0,' h =',hh,' h0old =',h0old,' dg =',dg 1002 write(std_out,*) 'facf =',facf 1003 msg = 'TIME EXCEEDED 5 min IN FOLLOW' 1004 ABI_ERROR(msg) 1005 end if 1006 1007 if (dg<aim_dgmin) then 1008 write(std_out,*) 'gradient < dgmin ',dg,' < ',aim_dgmin 1009 fin=.true. 1010 iat=0 1011 ipos=0 1012 ! testing for the CP 1013 if (npc>0) then 1014 call critic(aim_dtset,vv,ev,zz,aim_dmaxcrit,ires,0) 1015 if (ires==0) then 1016 do jj=1,npc 1017 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 1018 dist=vnorm(pom,0) 1019 if (dist<aim_tiny) cycle madw 1020 end do 1021 write(std_out,*) 'C.P. found !!' 1022 npc=npc+1 1023 do jj=1,3 1024 pc(jj,npc)=vv(jj) 1025 evpc(jj,npc)=ev(jj) 1026 do kk=1,3 1027 zpc(kk,jj,npc)=zz(kk,jj) 1028 end do 1029 end do 1030 i1=ev(1)/abs(ev(1)) 1031 i2=ev(2)/abs(ev(2)) 1032 i3=ev(3)/abs(ev(3)) 1033 icpc(npc)=i1+i2+i3 1034 if (icpc(npc)==-3) then ! pseudoatom handling 1035 npcm3=npcm3+1 1036 write(std_out,*) 'Pseudo-atom found !!' 1037 end if 1038 1039 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 1040 write(22,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 1041 & ev(1)+ev(2)+ev(3),rho 1042 write(std_out,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 1043 & ev(1)+ev(2)+ev(3),rho 1044 else 1045 write(std_out,*) 'C.P. not found !!' 1046 end if 1047 end if 1048 1049 cycle madw 1050 end if 1051 1052 hh=h0/dg 1053 if (ldeb.or.deb) write(std_out,*) 'h= ',hh,' h0= ',h0,' dg= ',dg 1054 if (hh>aim_hmax) hh=aim_hmax 1055 ! step modifications 1056 1057 hh=hh*facf 1058 if (hh>(hold*aim_hmult)) then 1059 hh=hold*aim_hmult 1060 end if 1061 1062 do ii=1,3 1063 vold(ii)=vv(ii) 1064 end do 1065 1066 nit=0 1067 hold=hh 1068 1069 ! one step following the gradient line 1070 ! 1071 call onestep(vv,rho,grho,hh,np,npmax,deltar) 1072 do while (((np>npmax).or.(deltar>aim_stmax)).and.(deltar>aim_dmin)) 1073 nit=nit+1 1074 if (nit>5) then 1075 if (deltar>aim_stmax) then 1076 write(std_out,*) 'nit > 5 and deltar > stmax nit=',nit 1077 else 1078 write(std_out,*) 'nit > 5 and np > npmax nit=',nit 1079 end if 1080 end if 1081 do ii=1,3 1082 vv(ii)=vold(ii) 1083 end do 1084 hh=hh*0.3 1085 call onestep(vv,rho,grho,hh,np,npmax,deltar) 1086 end do 1087 1088 1089 nstep=nstep+1 1090 if (ldeb.or.deb) write(std_out,*) 'h= ',hh 1091 1092 fac2=hh/hold 1093 if (fac2>=1._dp) then 1094 facf=facf*1.2 1095 else 1096 if (fac2>=aim_facmin) then 1097 facf=fac2 1098 else 1099 facf=aim_facmin 1100 end if 1101 end if 1102 1103 if (deb.or.ldeb) then 1104 write(std_out,*) ':POS ',vv 1105 write(std_out,*) ':RBPOS ',vt1 1106 write(std_out,*) ':GRAD ',grho 1107 end if 1108 1109 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 1110 dg=vnorm(grho,0) 1111 pom(:)=vv(:)-xatm(:,iatinit)-atp(:,iposinit) 1112 1113 if (iat /= 0) then 1114 fin=.true. 1115 write(std_out,*) 'r < rmin iat=',iat,' ipos=',ipos 1116 cycle madw 1117 end if 1118 1119 if (rho<aim_rhomin) then 1120 fin=.true. 1121 write(std_out,*) 'charge < rhomin ',rho,' < ',aim_rhomin 1122 if (rho<zero) then 1123 ABI_ERROR('RHO < 0 !!!') 1124 end if 1125 iat=0 1126 ipos=0 1127 cycle madw 1128 end if 1129 1130 if (npcm3>0) then 1131 do jj=1,npc 1132 if (icpc(jj)==(-3)) then 1133 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 1134 dist=vnorm(pom,0) 1135 if (dist<(aim_dtset%rmin**2*0.1)) then 1136 iat=0 1137 ipos=0 1138 fin=.true. 1139 write(std_out,*) 'We are inside a pseudo-atom' 1140 cycle madw 1141 end if 1142 end if 1143 end do 1144 end if 1145 1146 nsi=nsi+1 1147 1148 ! surface checking 1149 1150 if (srch.and.(nsi>=nsimax)) then 1151 nsi=0 1152 ith=0 1153 iph=0 1154 do ii=1,3 1155 vt(ii)=vv(ii)-xatm(ii,iatinit) 1156 end do 1157 xy=vt(1)*vt(1)+vt(2)*vt(2) 1158 xyz=xy+vt(3)*vt(3) 1159 xyz=sqrt(xyz) 1160 if (xy<aim_snull) then 1161 vcth=1._dp 1162 if (vt(3)<0._dp) vcth=-vcth 1163 vph=0._dp 1164 else 1165 vcth=vt(3)/xyz 1166 vph=atan2(vt(2),vt(1)) 1167 end if 1168 vth=acos(vcth) 1169 if (vth<th(1)) then 1170 ith=0 1171 else 1172 if (vth>th(nth)) then 1173 ith=nth 1174 else 1175 do ii=2,nth 1176 if (vth<th(ii)) then 1177 ith=ii-1 1178 exit 1179 end if 1180 end do 1181 end if 1182 end if 1183 1184 if (vph<ph(1)) then 1185 iph=0 1186 else 1187 if (vph>ph(nph)) then 1188 iph=nph 1189 else 1190 do ii=2,nph 1191 if (vph<ph(ii)) then 1192 iph=ii-1 1193 exit 1194 end if 1195 end do 1196 end if 1197 end if 1198 1199 stemp=(iph>0).and.(iph<nph) 1200 stemp=stemp.and.(ith>0).and.(ith<nth) 1201 1202 if (stemp) then 1203 stemp2=rs(ith,iph)>0._dp 1204 stemp2=stemp2.and.(rs(ith+1,iph)>0._dp) 1205 stemp2=stemp2.and.(rs(ith+1,iph+1)>0._dp) 1206 stemp2=stemp2.and.(rs(ith,iph+1)>0._dp) 1207 if (stemp2) then 1208 dth=th(ith+1)-th(ith) 1209 dph=ph(iph+1)-ph(iph) 1210 rsmed=rs(ith,iph)*(th(ith+1)-vth)/dth*(ph(iph+1)-vph)/dph 1211 rsmed=rsmed+rs(ith+1,iph)*(vth-th(ith))/dth*(ph(iph+1)-vph)/dph 1212 rsmed=rsmed+rs(ith+1,iph+1)*(vth-th(ith))/dth*(vph-ph(iph))/dph 1213 rsmed=rsmed+rs(ith,iph+1)*(th(ith+1)-vth)/dth*(vph-ph(iph))/dph 1214 if (rsmed>xyz) then 1215 write(std_out,*) 'We are inside the surface' 1216 iat=iatinit 1217 ipos=iposinit 1218 else 1219 write(std_out,*) 'We are outside the surface' 1220 iat=0 1221 ipos=0 1222 end if 1223 fin=.true. 1224 cycle madw 1225 end if 1226 end if 1227 end if 1228 1229 end do madw 1230 1231 1232 srch=srchold 1233 ldeb=ldebold 1234 h0=h0old 1235 1236 1237 end subroutine aim_follow
m_bader/bcp_type [ Types ]
NAME
bcp_type
FUNCTION
a "bonding critical point" for aim
SOURCE
189 type, private :: bcp_type 190 191 ! Integer 192 integer :: iat ! number of the bonding atom inside a primitive cell 193 integer :: ipos ! number of the primitive cell of the bonding atom 194 195 ! Real 196 real(dp) :: chg ! charge at the critical point 197 real(dp) :: diff(3) ! three distances : AT-CP,BAT-CP,AT-BAT 198 real(dp) :: ev(3) ! eigenvalues of the Hessian 199 real(dp) :: pom(3) ! position of the bonding atom 200 real(dp) :: rr(3) ! position of the bcp 201 real(dp) :: vec(3,3) ! eigenvectors of the Hessian 202 real(dp) :: vv(3) ! position of the bcp relative to the central atom 203 204 end type bcp_type
m_bader/bschg1 [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
bschg1
FUNCTION
bschg1: Vector transformation of coordinates
SOURCE
5813 subroutine bschg1(vv,dir) 5814 5815 !Arguments ------------------------------------ 5816 !scalars 5817 integer,intent(in) :: dir 5818 !arrays 5819 real(dp),intent(inout) :: vv(3) 5820 5821 !Local variables ------------------------------ 5822 !scalars 5823 integer :: ii 5824 !arrays 5825 real(dp) :: vt(3) 5826 5827 ! ********************************************************************* 5828 5829 if (dir==1) then 5830 do ii=1,3 5831 vt(ii)=rprimd(ii,1)*vv(1)+rprimd(ii,2)*vv(2)+rprimd(ii,3)*vv(3) 5832 end do 5833 elseif (dir==-1) then 5834 do ii=1,3 5835 vt(ii)=ivrprim(ii,1)*vv(1)+ivrprim(ii,2)*vv(2)+ivrprim(ii,3)*vv(3) 5836 end do 5837 elseif (dir==2) then 5838 do ii=1,3 5839 vt(ii)=trivrp(ii,1)*vv(1)+trivrp(ii,2)*vv(2)+trivrp(ii,3)*vv(3) 5840 end do 5841 else 5842 ABI_ERROR('Transformation of coordinates') 5843 end if 5844 vv(:)=vt(:) 5845 5846 end subroutine bschg1
m_bader/bschg2 [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
bschg2
FUNCTION
bschg2: Matrix transformation of coordinates
SOURCE
5858 subroutine bschg2(aa,dir) 5859 5860 !Arguments ------------------------------------ 5861 !scalars 5862 integer,intent(in) :: dir 5863 !arrays 5864 real(dp),intent(inout) :: aa(3,3) 5865 5866 !Local variables ------------------------------ 5867 !arrays 5868 real(dp) :: bb(3,3) 5869 5870 ! ********************************************************************* 5871 5872 if (dir==1) then 5873 call mprod(aa,ivrprim,bb) 5874 call mprod(rprimd,bb,aa) 5875 elseif (dir==2) then 5876 call mprod(aa,ivrprim,bb) 5877 call mprod(trivrp,bb,aa) 5878 elseif (dir==-1) then 5879 call mprod(aa,rprimd,bb) 5880 call mprod(ivrprim,bb,aa) 5881 else 5882 ABI_ERROR("transformation of coordinates") 5883 end if 5884 end subroutine bschg2
m_bader/consist [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
consist
FUNCTION
Checking of the consistency between the values of input variables
INPUTS
aim_dtset= the structured entity containing all input variables tstngr= information about the test on the ngrid input variable tstvpt= information about the test on the vpts input variable
OUTPUT
(only checking : print error message and stop if there is a problem) WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
1260 subroutine consist(aim_dtset,tstngr,tstvpt) 1261 1262 !Arguments ------------------------------------ 1263 !scalars 1264 integer,intent(in) :: tstngr,tstvpt 1265 type(aim_dataset_type),intent(in) :: aim_dtset 1266 1267 !Local variables ------------------------------ 1268 1269 ! ********************************************************************* 1270 1271 !write(std_out,*) tstngr, tstvpt 1272 1273 if (((aim_dtset%denout/=0).or.(aim_dtset%lapout/=0)).and.((tstngr < 1).or.(tstvpt < 2))) then 1274 ABI_ERROR('in input1 - I cannot do the output !') 1275 end if 1276 if ((aim_dtset%denout > 0).and.(aim_dtset%lapout>0)) then 1277 if (aim_dtset%denout/=aim_dtset%lapout) then 1278 write(std_out,*) 'ERROR in input - when both denout and lapout are positive non-zero,' 1279 write(std_out,*) 'they must be equal.' 1280 ABI_ERROR("Aborting now") 1281 end if 1282 if ((tstvpt < aim_dtset%denout+1).or.(tstngr < aim_dtset%denout)) then 1283 write(std_out,*) 'ERROR in input2 - I cannot do the output !' 1284 ABI_ERROR("Aborting now") 1285 end if 1286 elseif (aim_dtset%denout > 0) then 1287 if ((tstvpt < aim_dtset%denout+1).or.(tstngr < aim_dtset%denout)) then 1288 write(std_out,*) 'ERROR in input - I cannot do the output !' 1289 ABI_ERROR("Aborting now") 1290 end if 1291 elseif (aim_dtset%lapout > 0) then 1292 if ((tstvpt < aim_dtset%lapout+1).or.(tstngr < aim_dtset%lapout)) then 1293 write(std_out,*) 'ERROR in input - I cannot do the output !' 1294 ABI_ERROR("Aborting now") 1295 end if 1296 end if 1297 1298 if ((aim_dtset%isurf==1).and.(aim_dtset%crit==0)) then 1299 write(std_out,*) 'ERROR in input - must have crit/=0 for isurf==1' 1300 ABI_ERROR("Aborting now") 1301 end if 1302 1303 if (((aim_dtset%ivol/=0).or.(aim_dtset%irho/=0)).and.(aim_dtset%isurf==0)) then 1304 ABI_ERROR('in input - I cannot integrate without surface !') 1305 end if 1306 1307 end subroutine consist
m_bader/cpdrv [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
cpdrv
FUNCTION
Critical points (CPs) searching driver First Bond CPs are searched for each pair atom-its neighbor (distance cutoff=maxatdst) then Ring CPs for each pair of BCPs and finally Cage CPs for each pair of RCPs.
INPUTS
aim_dtset= the structured entity containing all input variables
OUTPUT
SIDE EFFECTS
this routine treat information contained in the aim_prom module WARNING This file does not follow the ABINIT coding rules (yet)
TODO
Should combine parts of code that are similar ...
SOURCE
1337 subroutine cpdrv(aim_dtset) 1338 1339 !Arguments ------------------------------------ 1340 !scalars 1341 type(aim_dataset_type),intent(in) :: aim_dtset 1342 1343 !Local variables ------------------------------ 1344 !scalars 1345 integer :: iat,iatinit,ii,inxat,inxcell,ipair,ipos,iposinit,ires,jj,kk,nb,nb_now 1346 integer :: nn,nstep,nvs,me,nproc,ierr 1347 real(dp) :: candidate,chg,diff1,diff2,diff3,dist,prj,rtdiff,ss,tt0,wall 1348 logical :: srch=.false. 1349 !arrays 1350 integer :: ibat(nnpos*natom),inatm(nnpos*natom),incell(nnpos*natom) 1351 integer :: ipibat(nnpos*natom) 1352 integer,allocatable :: indexcp(:),nr(:) 1353 real(dp) :: bmin(natom),dif(3),dists(nnpos*natom),ev(3),evec(3,3),grho(3) 1354 real(dp) :: hrho(3,3),pom(3),rr(3),uu(3),vv(3),xorig(3) 1355 real(dp),allocatable :: buffer(:,:),sortguide(:) 1356 !no_abirules 1357 !Warning : bcp_type should be transformed to cp_type 1358 type(bcp_type),allocatable :: bcp(:),ccp(:),cp_tmp(:),rcp(:) 1359 1360 !************************************************************************ 1361 1362 me=xmpi_comm_rank(xmpi_world) 1363 nproc=xmpi_comm_size(xmpi_world) 1364 1365 !Consider the critical points starting from atom #batom 1366 inxat=aim_dtset%batom 1367 slc=-1 1368 rminl(:)=aim_dtset%rmin 1369 bmin(:)=0._dp 1370 ttcp=0._dp 1371 1372 write(std_out,*) 1373 write(std_out,*) "CRITICAL POINTS ANALYSIS" 1374 write(std_out,*) "========================" 1375 write(std_out,*) 1376 1377 write(untout,*) 1378 write(untout,*) "CRITICAL POINTS ANALYSIS" 1379 write(untout,*) "========================" 1380 write(untout,*) 1381 1382 1383 xorig(:)=xatm(:,inxat) 1384 1385 call timein(tt0,wall) 1386 1387 !Searching the neighbouring atoms 1388 1389 if (aim_dtset%crit > 0) then 1390 nvs=0 1391 do ii=1,nnpos 1392 do jj=1,natom 1393 dist=0._dp 1394 dif(:)=xatm(:,inxat)-xatm(:,jj)-atp(:,ii) 1395 dif(:)=dif(:)/aim_dtset%scal(:) 1396 dist=vnorm(dif,0) 1397 if (dist < tol6 ) then 1398 inxcell=ii 1399 elseif (dist < maxatdst) then 1400 nvs=nvs+1 1401 dists(nvs)=dist 1402 inatm(nvs)=jj 1403 incell(nvs)=ii 1404 end if 1405 end do 1406 end do 1407 1408 write(std_out,*) "ATOM:" 1409 write(std_out,*) 'inxat :', inxat, 'inxcell :', inxcell 1410 write(std_out, '(3es16.6)' ) (xorig(ii),ii=1,3) 1411 write(std_out,*) 1412 1413 write(untout,*) "ATOM:" 1414 write(untout,*) 'inxat :', inxat, 'inxcell :', inxcell 1415 write(untout, '(3es16.6)') (xorig(ii),ii=1,3) 1416 write(untout,*) 1417 1418 ABI_MALLOC(nr,(nvs)) 1419 do ii=1,nvs 1420 nr(ii)=ii 1421 end do 1422 1423 ! Ordering of the nearest neighbouring atoms 1424 call sort_dp(nvs,dists,nr,tol14) 1425 1426 nb=0 1427 write(std_out,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 1428 write(untout,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 1429 do ii=1,nvs 1430 nn=nr(ii) 1431 if (dists(ii) < maxatdst) then 1432 nb=nb+1 1433 ibat(nb)=inatm(nn) 1434 ipibat(nb)=incell(nn) 1435 write(std_out,*) ':NEIG ',inatm(nn),incell(nn),dists(ii) 1436 write(untout,'(" ",2I6,F16.8)')inatm(nn),incell(nn),dists(ii) 1437 else 1438 exit 1439 end if 1440 end do 1441 1442 ! SEARCHING BCP 1443 ABI_MALLOC(bcp,(nb)) 1444 nbcp=0 1445 iatinit=inxat 1446 iposinit=inxcell 1447 bcp(:)%iat=0 1448 bcp(:)%ipos=0 1449 1450 write(std_out,*) 1451 write(std_out,*) "BONDING CRITICAL POINTS (BCP)" 1452 write(std_out,*) "=============================" 1453 write(std_out,*) 1454 1455 write(untout,*) 1456 write(untout,*) "BONDING CRITICAL POINTS (BCP)" 1457 write(untout,*) "=============================" 1458 write(untout,*) 1459 1460 srbcp: do ii=1,nb 1461 1462 ! Start the search for BCP from the midistance between the atom 1463 ! and his neighbor. 1464 vv(:)=(xatm(:,inxat)+xatm(:,ibat(ii))+atp(:,ipibat(ii)))/2._dp 1465 1466 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,-1) 1467 1468 if (ires==0) then 1469 ! Testing if CP is already known 1470 if (nbcp > 0) then 1471 do jj=1,nbcp 1472 pom(:)=vv(:)-bcp(jj)%rr(:)-xorig(:) 1473 dist=vnorm(pom,0) 1474 if (dist < aim_dtset%dpclim) then 1475 write(std_out,*) 'BCP already known !' 1476 cycle srbcp 1477 end if 1478 end do 1479 end if 1480 rr(:)=vv(:)-xorig(:) 1481 ss=vnorm(rr,0) 1482 if (ss > maxcpdst) then 1483 write(std_out, '(a,es16.6,a,es16.6)' ) 'BCP distance from atom,',ss,', exceed maxcpdst =',maxcpdst 1484 cycle srbcp 1485 end if 1486 nn=0 1487 do jj=1,3 1488 nn=nn+ev(jj)/abs(ev(jj)) 1489 end do 1490 write(std_out, '(a,3es16.6,i4)') ' vv(1:3), nn',(vv(jj), jj=1,3), nn 1491 write(std_out, '(a,3es16.6)') 'ev: ', (ev(jj), jj=1,3) 1492 if (nn /= -1) then 1493 write(std_out,*) ' The trial critical point is not a BCP !' 1494 cycle srbcp 1495 end if 1496 write(std_out, '(a,3es16.6)' ) 'evec(:,1): ',(evec(jj,1), jj=1,3) 1497 pom(:)=evec(:,1) 1498 dist=vnorm(pom,0) 1499 prj=dot_product(evec(:,1),rr) 1500 write(std_out,*) 'prj:', prj, vnorm(evec(:,1),0) 1501 dist=vnorm(evec(:,1),0) 1502 uu(:)=vv(:)-sign(aim_epstep,prj)*evec(:,1)/dist 1503 1504 ! Testing whether this BCP "is bonded" to the considered atom 1505 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1506 ! write(std_out,*) 'do', iat, ipos 1507 ! if ((iat==0).or.(ipos==0)) cycle 1508 ! write(std_out,*) 'APOS: ',(xatm(jj,iat)+atp(jj,ipos), jj=1,3) 1509 if ((iat/=inxat).or.(inxcell/=ipos)) then 1510 write(std_out,*) ' The trial BCP is not bonded to the Bader atom' 1511 cycle srbcp 1512 end if 1513 1514 ! A new BCP has been found ! 1515 nbcp=nbcp+1 1516 1517 ! Searching for the second bonded atom 1518 ss=vnorm(rr,0) 1519 diff1=ss 1520 diff3=dists(ii) 1521 uu(:)=vv(:)+sign(aim_epstep,prj)*evec(:,1)/dist 1522 if ((abs(bmin(iat))<1.0d-12).or.( ss<bmin(iat))) then 1523 bmin(iat)=ss 1524 end if 1525 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1526 if ((iat==0).or.(ipos==0)) then 1527 write(std_out,*) ' The trial BCP is not bonded to a bonding atom !' 1528 ! cycle srbcp 1529 end if 1530 pom(:)=vv(:)-xatm(:,iat)-atp(:,ipos) 1531 ss=vnorm(pom,0) 1532 diff2=ss 1533 pom(:)=xorig(:)-xatm(:,iat)-atp(:,ipos) 1534 diff3=vnorm(pom,0) 1535 rtdiff=diff1/diff3 1536 if ((abs(bmin(iat))<1.0d-12).or.(ss<bmin(iat))) then 1537 bmin(iat)=ss 1538 end if 1539 pom(:)=xatm(:,iat)+atp(:,ipos) 1540 1541 ! Store more results, for coherent, and portable output 1542 bcp(nbcp)%iat=iat 1543 bcp(nbcp)%ipos=ipos 1544 bcp(nbcp)%chg=chg 1545 bcp(nbcp)%diff(1)=diff1 1546 bcp(nbcp)%diff(2)=diff2 1547 bcp(nbcp)%diff(3)=diff3 1548 bcp(nbcp)%ev(:)=ev(:) 1549 bcp(nbcp)%pom(:)=pom(:) 1550 bcp(nbcp)%rr(:)=rr(:) 1551 bcp(nbcp)%vec(:,:)=evec(:,:) 1552 bcp(nbcp)%vv(:)=vv(:) 1553 ! Warning : iat, ipos might be modified by this call 1554 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1555 bcp(nbcp)%chg=chg 1556 1557 end if ! ires==0 1558 end do srbcp 1559 1560 if(nbcp>0)then 1561 1562 ! Order the BCP. CPs should appear by increasing values of x,y,z , the latter 1563 ! varying the fastest 1564 ABI_MALLOC(sortguide,(nbcp)) 1565 ABI_MALLOC(indexcp,(nbcp)) 1566 ABI_MALLOC(cp_tmp,(nbcp)) 1567 do ii=3,1,-1 1568 ! DEBUG 1569 ! write(std_out,*)' cpdrv : sort on index ii=',ii 1570 ! ENDDEBUG 1571 1572 do jj=1,nbcp 1573 ! DEBUG 1574 ! write(std_out,*)bcp(jj)%vv(:) 1575 ! ENDDEBUG 1576 sortguide(jj)=bcp(jj)%vv(ii) 1577 indexcp(jj)=jj 1578 end do 1579 1580 ! Try to be platform-independent. Might need a larger tolerance. 1581 call sort_dp(nbcp,sortguide,indexcp,tol3) 1582 do jj=1,nbcp 1583 cp_tmp(jj)=bcp(indexcp(jj)) 1584 end do 1585 do jj=1,nbcp 1586 bcp(jj)=cp_tmp(jj) 1587 end do 1588 end do 1589 ! DEBUG 1590 ! write(std_out,*)' cpdrv : after the sort ' 1591 ! do jj=1,nbcp 1592 ! write(std_out,*)bcp(jj)%vv(:) 1593 ! end do 1594 ! ENDDEBUG 1595 1596 1597 ! Output the info about the BCP 1598 do jj=1,nbcp 1599 write(untout,'(" Bonded atom (BAT) (indxatm,indxcell,position): ",/,2I6,3F16.8)')& 1600 & bcp(jj)%iat,bcp(jj)%ipos,bcp(jj)%pom(:) 1601 write(untout,'("%Bonding CP: ",3F16.8)') bcp(jj)%vv(:) 1602 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') bcp(jj)%ev(:) 1603 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1604 & ' Eigenvec. of Hessian:',char(10),& 1605 & '-',bcp(jj)%vec(1,:),char(10),& 1606 & '-',bcp(jj)%vec(2,:),char(10),& 1607 & '-',bcp(jj)%vec(3,:),char(10) 1608 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1609 & bcp(jj)%chg, bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 1610 write(untout,'("%Relative position of BCP (AT-CP,BAT-CP,AT-BAT,relative(AT): ",/,4F16.8)') & 1611 & bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 1612 write(untout,*) "********************************************************************" 1613 write(std_out,'(/," BCP: ",3F10.6,3E12.4,E12.4,/)') & 1614 & bcp(jj)%rr(:),bcp(jj)%ev(:),bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 1615 write(std_out,'(":DISPC ",4F12.6)') bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 1616 end do 1617 1618 ABI_FREE(cp_tmp) 1619 ABI_FREE(indexcp) 1620 ABI_FREE(sortguide) 1621 1622 end if ! nbcp>0 1623 1624 if (abs(bmin(inxat))>1.0d-12) then 1625 rminl(inxat)=aim_dtset%coff1*bmin(inxat) 1626 r0=bmin(inxat) 1627 else 1628 r0=0._dp 1629 end if 1630 1631 ! !AD-HOC PARAMETER 1632 1633 do ii=1,natom 1634 if ((abs(bmin(ii))>1.0d-12).and.(ii /= inxat)) rminl(ii)=aim_dtset%coff2*bmin(ii) 1635 end do 1636 1637 ! END WARNING 1638 1639 write(std_out,*) ' number of BCP:', nbcp 1640 write(untout,'(" Number of BCP found: ",I4)') nbcp 1641 nn=nbcp*(nbcp-1)*(nbcp-2)/6 1642 if (bit_size(ii) <= nbcp+1) then 1643 ABI_ERROR("b-test!") 1644 end if 1645 1646 ! SEARCHING RCP 1647 1648 write(std_out,*) 1649 write(std_out,*) "RING CRITICAL POINTS (RCP)" 1650 write(std_out,*) "=============================" 1651 write(std_out,*) 1652 1653 write(untout,*) 1654 write(untout,*) "RING CRITICAL POINTS (RCP)" 1655 write(untout,*) "=============================" 1656 write(untout,*) 1657 1658 nrcp=0 1659 if(aim_dtset%crit==1)nb_now=nbcp 1660 if(aim_dtset%crit==2)nb_now=nb 1661 ! DEBUG 1662 ! nb_now=nbcp 1663 ! ENDDEBUG 1664 nn=nb_now*(nb_now-1)/2 1665 ABI_MALLOC(rcp,(nn)) 1666 1667 ! Loop on pairs of BCP or atoms 1668 ipair=0 1669 ABI_MALLOC(buffer,(16,nn)) 1670 buffer=zero 1671 1672 ! DEBUG 1673 ! write(std_out,*)ch10,ch10,' drvcpr : enter loop to search for RCPs,nb_now,nn=',nb_now,nn 1674 ! ENDDEBUG 1675 1676 do ii=1,nb_now-1 1677 srcp1: do jj=ii+1,nb_now 1678 ipair=ipair+1 1679 if(mod(ipair,nproc)==me)then 1680 if (aim_dtset%crit==1) then 1681 vv(:)=xorig(:)+(bcp(ii)%rr(:)+bcp(jj)%rr(:))/2._dp 1682 else if (aim_dtset%crit==2) then 1683 vv(:)=xorig(:)*half+(xatm(:,ibat(ii))+atp(:,ipibat(ii))+xatm(:,ibat(jj))+atp(:,ipibat(jj)))*quarter 1684 end if 1685 1686 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,1) 1687 1688 if(ires==1)then 1689 cycle srcp1 1690 end if 1691 1692 ! Check that it is within the maximum allowed distance for a CP 1693 rr(:)=vv(:)-xorig(:) 1694 ss=vnorm(rr,0) 1695 if (ss > maxcpdst) then 1696 write(std_out,*) 'RCP distance from atom exceed maxcpdst !' 1697 cycle srcp1 1698 end if 1699 ! Check that it is a RCP 1700 nn=0 1701 do kk=1,3 1702 nn=nn+ev(kk)/abs(ev(kk)) 1703 end do 1704 if (nn /= 1) then 1705 write(std_out,*) ' the critical point that is found is not a RCP ' 1706 cycle srcp1 1707 end if 1708 ! Might be the same RCP than one already found on the same processor 1709 if (nrcp > 0) then 1710 do kk=1,nrcp 1711 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 1712 dist=vnorm(pom,0) 1713 if (dist < aim_dtset%dpclim) then 1714 write(std_out,*) ':RCP already known' 1715 cycle srcp1 1716 end if 1717 end do 1718 end if 1719 ! If crit==2, check that it is on the Bader surface 1720 if (aim_dtset%crit==2) then 1721 uu(:)=vv(:)-aim_epstep*rr(:)/ss 1722 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1723 if ((iat/=inxat).or.(inxcell/=ipos))then 1724 write(std_out,*) ' RCP is not on the Bader surface (outside of it)' 1725 cycle srcp1 1726 end if 1727 end if 1728 nrcp=nrcp+1 1729 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 1730 1731 buffer(1:3,ipair)=vv 1732 buffer(4:6,ipair)=ev 1733 buffer(7:9,ipair)=evec(:,1) 1734 buffer(10:12,ipair)=evec(:,2) 1735 buffer(13:15,ipair)=evec(:,3) 1736 buffer(16,ipair)=one 1737 1738 ! DEBUG 1739 ! write(std_out,*)ch10,ch10,' drvcpr : ipair,candidate=',ipair,candidate 1740 ! ENDDEBUG 1741 end if 1742 end do srcp1 1743 end do 1744 call xmpi_sum(buffer,xmpi_world,ierr) 1745 1746 nrcp=0 1747 ipair=0 1748 do ii=1,nb_now-1 1749 srcp: do jj=ii+1,nb_now 1750 ipair=ipair+1 1751 candidate=buffer(16,ipair) 1752 1753 ! One CP has been found, must make tests to see whether it is a new RCP 1754 if (nint(candidate)==1) then 1755 1756 vv=buffer(1:3,ipair) 1757 ev=buffer(4:6,ipair) 1758 evec(:,1)=buffer(7:9,ipair) 1759 evec(:,2)=buffer(10:12,ipair) 1760 evec(:,3)=buffer(13:15,ipair) 1761 1762 ! Check that it is not the same as a previous one 1763 if (nrcp > 0) then 1764 do kk=1,nrcp 1765 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 1766 dist=vnorm(pom,0) 1767 if (dist < aim_dtset%dpclim) then 1768 write(std_out,*) ':RCP already known' 1769 cycle srcp 1770 end if 1771 end do 1772 end if 1773 1774 ! A new RCP has been found ! 1775 nrcp=nrcp+1 1776 1777 ! DEBUG 1778 ! write(std_out,*)' drvcpr : A new RCP has been found, for kk=',kk 1779 ! ENDDEBUG 1780 1781 1782 rcp(nrcp)%iat=iat 1783 rcp(nrcp)%ipos=ipos 1784 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 1785 rcp(nrcp)%vec(:,:)=evec(:,:) 1786 rcp(nrcp)%ev(:)=ev(:) 1787 rcp(nrcp)%vv(:)=vv(:) 1788 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1789 rcp(nrcp)%chg=chg 1790 1791 end if ! ires==0 1792 end do srcp ! jj=ii+2,nb_now 1793 end do ! ii=1,nb_now-1 1794 1795 ABI_FREE(buffer) 1796 1797 if(nrcp>0)then 1798 1799 ! Order the RCP. CPs should appear by increasing values of x,y,z , the latter 1800 ! varying the fastest 1801 ABI_MALLOC(sortguide,(nrcp)) 1802 ABI_MALLOC(indexcp,(nrcp)) 1803 ABI_MALLOC(cp_tmp,(nrcp)) 1804 do ii=3,1,-1 1805 ! DEBUG 1806 ! write(std_out,*)' cpdrv : sort on index ii=',ii 1807 ! ENDDEBUG 1808 do jj=1,nrcp 1809 1810 ! DEBUG 1811 ! write(std_out,*)rcp(jj)%vv(:) 1812 ! ENDDEBUG 1813 1814 ! Try to be platform-independent. Might need a larger tolerance. 1815 sortguide(jj)=rcp(jj)%vv(ii) 1816 indexcp(jj)=jj 1817 end do 1818 call sort_dp(nrcp,sortguide,indexcp,tol3) 1819 do jj=1,nrcp 1820 cp_tmp(jj)=rcp(indexcp(jj)) 1821 end do 1822 do jj=1,nrcp 1823 rcp(jj)=cp_tmp(jj) 1824 end do 1825 end do 1826 1827 ! DEBUG 1828 ! write(std_out,*)' cpdrv : after the sort ' 1829 ! do jj=1,nrcp 1830 ! write(std_out,*)rcp(jj)%vv(:) 1831 ! end do 1832 ! ENDDEBUG 1833 1834 1835 ! Write the Ring Critical Point information 1836 do jj=1,nrcp 1837 write(untout,'(";Ring CP: ",3F16.8)') rcp(jj)%vv(:) 1838 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') rcp(jj)%ev(:) 1839 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1840 & ' Eigenvec. of Hessian:',char(10),& 1841 & '-',rcp(jj)%vec(1,:),char(10),& 1842 & '-',rcp(jj)%vec(2,:),char(10),& 1843 & '-',rcp(jj)%vec(3,:),char(10) 1844 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1845 & rcp(jj)%chg, rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 1846 write(untout,*) "********************************************************************" 1847 write(std_out,'(/," RCP: ",3F10.6,3E12.4,E12.4,/)') & 1848 & rcp(jj)%rr(:),rcp(jj)%ev(:),rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 1849 end do 1850 1851 ABI_FREE(cp_tmp) 1852 ABI_FREE(indexcp) 1853 ABI_FREE(sortguide) 1854 1855 end if ! nrcp>0 1856 1857 write(untout,'(" Number of RCP found: ",I4)') nrcp 1858 write(std_out,*) ' Number of RCP:', nrcp 1859 1860 ! SEARCHING CCP 1861 1862 write(std_out,*) 1863 write(std_out,*) "CAGE CRITICAL POINTS (CCP)" 1864 write(std_out,*) "=============================" 1865 write(std_out,*) 1866 1867 write(untout,*) 1868 write(untout,*) "CAGE CRITICAL POINTS (CCP)" 1869 write(untout,*) "=============================" 1870 write(untout,*) 1871 1872 1873 nn=nrcp*(nrcp-1)/2 1874 ABI_MALLOC(ccp,(nn)) 1875 1876 nccp=0 1877 do ii=1,nrcp-1 1878 srccp: do jj=ii+1,nrcp 1879 vv(:)=xorig(:)+(rcp(ii)%rr(:)+rcp(jj)%rr(:))/2._dp 1880 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,3) 1881 if (ires==0) then 1882 rr(:)=vv(:)-xorig(:) 1883 ss=vnorm(rr,0) 1884 if (ss > maxcpdst) then 1885 write(std_out,*) 'CCP distance from atom exceed maxcpdst !' 1886 cycle srccp 1887 end if 1888 nn=0 1889 do kk=1,3 1890 nn=nn+ev(kk)/abs(ev(kk)) 1891 end do 1892 if (nn /= 3) then 1893 write(std_out,*) ' the critical point that is found is not a CCP ' 1894 cycle srccp 1895 end if 1896 1897 if (nccp > 0) then 1898 do kk=1,nccp 1899 pom(:)=vv(:)-ccp(kk)%rr(:)-xorig(:) 1900 dist=vnorm(pom,0) 1901 if (dist < aim_dtset%dpclim) then 1902 write(std_out,*) ':CCP already known' 1903 cycle srccp 1904 end if 1905 end do 1906 end if 1907 if (aim_dtset%crit==2) then 1908 uu(:)=vv(:)-aim_epstep*rr(:)/ss 1909 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1910 if ((iat/=inxat).or.(inxcell/=ipos)) then 1911 write(std_out,*) ' This CCP is not on the Bader surface (outside of it)' 1912 cycle srccp 1913 end if 1914 end if 1915 1916 nccp=nccp+1 1917 1918 ccp(nccp)%iat=iat 1919 ccp(nccp)%ipos=ipos 1920 ccp(nccp)%rr(:)=vv(:)-xorig(:) 1921 ccp(nccp)%vec(:,:)=evec(:,:) 1922 ccp(nccp)%ev(:)=ev(:) 1923 ccp(nccp)%vv(:)=vv(:) 1924 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1925 ccp(nccp)%chg=chg 1926 1927 end if 1928 end do srccp 1929 end do 1930 1931 if(nccp>0)then 1932 1933 ! Order the CCP. CPs should appear by increasing values of x,y,z , the latter 1934 ! varying the fastest 1935 ABI_MALLOC(sortguide,(nccp)) 1936 ABI_MALLOC(indexcp,(nccp)) 1937 ABI_MALLOC(cp_tmp,(nccp)) 1938 do ii=3,1,-1 1939 do jj=1,nccp 1940 ! Try to be platform-independent. Might need a larger tolerance. 1941 sortguide(jj)=ccp(jj)%vv(ii) 1942 indexcp(jj)=jj 1943 end do 1944 call sort_dp(nccp,sortguide,indexcp,tol3) 1945 do jj=1,nccp 1946 cp_tmp(jj)=ccp(indexcp(jj)) 1947 end do 1948 do jj=1,nccp 1949 ccp(jj)=cp_tmp(jj) 1950 end do 1951 end do 1952 1953 ! Write the Cage Critical Point information 1954 do jj=1,nccp 1955 write(untout,'("%Cage CP: ",3F16.8)') ccp(jj)%vv(:) 1956 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') ccp(jj)%ev(:) 1957 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1958 & ' Eigenvec. of Hessian:',char(10),& 1959 & '-',ccp(jj)%vec(1,:),char(10),& 1960 & '-',ccp(jj)%vec(2,:),char(10),& 1961 & '-',ccp(jj)%vec(3,:),char(10) 1962 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1963 & ccp(jj)%chg, ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 1964 write(untout,*) "********************************************************************" 1965 write(std_out,'(/," CCP: ",3F10.6,3E12.4,E12.4,/)') & 1966 & ccp(jj)%rr(:),ccp(jj)%ev(:),ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 1967 end do 1968 1969 ABI_FREE(sortguide) 1970 ABI_FREE(indexcp) 1971 ABI_FREE(cp_tmp) 1972 1973 end if ! nccp>0 1974 1975 write(untout,'(" Number of CCP found: ",I4)') nccp 1976 write(std_out,*) 'Number of CCP:', nccp 1977 write(std_out,*) 1978 write(untout,*) 1979 write(std_out, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 1980 write(untout, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 1981 1982 write(std_out,*) 1983 write(std_out,*) "===============================" 1984 write(std_out,*) "END OF CRITICAL POINTS ANALYSIS" 1985 write(std_out,*) 1986 1987 write(untout,*) 1988 write(untout,*) "===============================" 1989 write(untout,*) "END OF CRITICAL POINTS ANALYSIS" 1990 write(untout,*) 1991 1992 1993 ! Output of the CPs 1994 1995 write(untc,'(I4, " :BCP''s, coordinates, laplacian eigs, type of bonding at., sum of lap.eigs., density")') nbcp 1996 do ii=1,nbcp 1997 write(untc,'(3F10.6,3E12.4,I4,2E12.4)') & 1998 & bcp(ii)%rr(:),bcp(ii)%ev(:),bcp(ii)%iat,bcp(ii)%ev(1)+bcp(ii)%ev(2)+bcp(ii)%ev(3),bcp(ii)%chg 1999 end do 2000 2001 write(untc,'(I4, " :RCP''s, coordinates, laplacian eigenvalues, sum of these, density")') nrcp 2002 do ii=1,nrcp 2003 vv(:)=rcp(ii)%rr(:)+xorig(:) 2004 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 2005 write(untc,'(3F10.6,3E12.4,2E12.4)') & 2006 & rcp(ii)%rr(:),rcp(ii)%ev(:),rcp(ii)%ev(1)+rcp(ii)%ev(2)+rcp(ii)%ev(3),rcp(ii)%chg 2007 end do 2008 2009 write(untc,'(I4, " :CCP''s coordinates, laplacian eigenvalues, sum of these, density")') nccp 2010 do ii=1,nccp 2011 vv(:)=ccp(ii)%rr(:)+xorig(:) 2012 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 2013 write(untc,'(3F10.6,3E12.4,2E12.4)') & 2014 & ccp(ii)%rr(:),ccp(ii)%ev(:),ccp(ii)%ev(1)+ccp(ii)%ev(2)+ccp(ii)%ev(3),ccp(ii)%chg 2015 end do 2016 2017 end if ! End the condition on aim_dtset%crit > 0 2018 2019 !Reading of the CPs from the file 2020 2021 if (aim_dtset%crit==-1) then 2022 read(untc,*) nbcp 2023 ABI_MALLOC(bcp,(nbcp)) 2024 do ii=1,nbcp 2025 read(untc,*) bcp(ii)%rr(:) 2026 end do 2027 read(untc,*) nrcp 2028 ABI_MALLOC(rcp,(nrcp)) 2029 do ii=1,nrcp 2030 read(untc,*) rcp(ii)%rr(:) 2031 end do 2032 read(untc,*) nccp 2033 ABI_MALLOC(ccp,(nccp)) 2034 do ii=1,nccp 2035 read(untc,*) ccp(ii)%rr(:) 2036 end do 2037 end if 2038 2039 do ii=1,nbcp 2040 pc(:,ii)=bcp(ii)%rr(:) 2041 icpc(ii)=-1 2042 end do 2043 do ii=1,nrcp 2044 pc(:,nbcp+ii)=rcp(ii)%rr(:) 2045 icpc(nbcp+ii)=1 2046 end do 2047 do ii=1,nccp 2048 pc(:,nbcp+nrcp+ii)=ccp(ii)%rr(:) 2049 icpc(nbcp+nrcp+ii)=3 2050 end do 2051 npc=nbcp+nrcp+nccp 2052 2053 !Checking 2054 2055 if (allocated(bcp)) then 2056 do ii=1,nbcp 2057 do jj=1,npc 2058 iat=bcp(ii)%iat 2059 ipos=bcp(ii)%ipos 2060 if ((iat/=0).and.(ipos/=0)) then 2061 pom(:)=pc(:,jj)+xorig(:)-xatm(:,iat)-atp(:,ipos) 2062 ss=aim_dtset%coff2*vnorm(pom,0) 2063 if (rminl(iat) >= ss) rminl(iat)=ss 2064 end if 2065 end do 2066 end do 2067 ABI_FREE(bcp) 2068 end if 2069 do ii=1,natom 2070 write(std_out,*) 'atom: ', ii, rminl(ii) 2071 end do 2072 2073 if(allocated(rcp)) then 2074 ABI_FREE(rcp) 2075 end if 2076 if(allocated(ccp)) then 2077 ABI_FREE(ccp) 2078 end if 2079 2080 !END CP ANALYSIS 2081 2082 call timein(ttcp,wall) 2083 ttcp=ttcp-tt0 2084 2085 end subroutine cpdrv
m_bader/critic [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
critic
FUNCTION
Search for a critical point starting from point vv
INPUTS
aim_dtset= the structured entity containing all input variables
dmax= maximal step
sort= 0(default) general CP searching (Newton-Raphson)
-1,1,3 searching of specific type CP (Popelier)
OUTPUT
ev= eigenvalues (ordered) of the Hessian in the final point
zz= eigenvectors of the Hessian in the final point
ires= if ires==0 => CP found
if ires==1 => CP not found within the maximum steps
SIDE EFFECTS
vv(3)= starting point and final point
SOURCE
2112 subroutine critic(aim_dtset,vv,ev,zz,dmax,ires,sort) 2113 2114 !Arguments ------------------------------------ 2115 !scalars 2116 integer,intent(in) :: sort 2117 integer,intent(out) :: ires 2118 real(dp),intent(in) :: dmax 2119 !arrays 2120 real(dp),intent(inout) :: vv(3) 2121 real(dp),intent(out) :: ev(3),zz(3,3) 2122 !no_abirules 2123 type(aim_dataset_type), intent(in) :: aim_dtset 2124 2125 !Local variables ------------------------------ 2126 !scalars 2127 integer :: iat,id,ii,info,ipos,istep,jii,jj,nrot 2128 real(dp),parameter :: evol=1.d-3 2129 real(dp) :: chg,dg,dltcmax,dv,dvold,rr,ss 2130 logical :: oscl,outof 2131 !arrays 2132 integer :: ipiv(3) 2133 real(dp) :: dc(3),ff(3),grho(3),hrho(3,3),lp(3),vold(3),vt(3),yy(3,3) 2134 real(dp),allocatable :: lamb(:),pom(:,:),pom2(:,:) 2135 2136 !************************************************************************ 2137 2138 !DEBUG 2139 !write(std_out,*)' critic : enter ' 2140 !ENDDEBUG 2141 oscl=.false. 2142 if (sort==3) then 2143 ABI_MALLOC(pom,(4,4)) 2144 ABI_MALLOC(pom2,(4,4)) 2145 ABI_MALLOC(lamb,(4)) 2146 elseif (sort/=0) then 2147 ABI_MALLOC(pom,(3,3)) 2148 ABI_MALLOC(pom2,(3,3)) 2149 ABI_MALLOC(lamb,(3)) 2150 end if 2151 2152 2153 deb=.false. 2154 istep=0 2155 ires=0 2156 2157 !DEBUG 2158 !write(std_out,'(":POSIN ",3F16.8)') vv 2159 !do jj=1,3 2160 !vt(jj)=rprimd(1,jj)*vv(1)+rprimd(2,jj)*vv(2)+rprimd(3,jj)*vv(3) 2161 !end do 2162 !write(std_out,'(":RBPOSIN ",3F16.8)') vt 2163 !ENDDEBUG 2164 2165 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2166 2167 !write(std_out,'(":GRAD ",3F16.8)') grho 2168 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 2169 2170 dg=1.0_dp 2171 dv=1.0_dp 2172 dg = vnorm(grho,0) 2173 2174 if (chg < aim_rhomin) then 2175 ires=1 2176 ! DEBUG 2177 ! write(std_out,*)' critic : exit, ires=1' 2178 ! ENDDEBUG 2179 return 2180 end if 2181 2182 !main cycle => limits (adhoc): 2183 !aim_dtset%lstep - minimal step 2184 !aim_dtset%lgrad - minimal norm of gradient 2185 !aim_maxstep - max number of steps 2186 2187 do while ((dv>aim_dtset%lstep).and.(dg>aim_dtset%lgrad).and.(istep<aim_maxstep)) 2188 istep=istep+1 2189 vold(:)=vv(:) 2190 dvold=dv 2191 ev(:)=0._dp 2192 yy(:,:)=0._dp 2193 call jacobi(hrho,3,3,ev,yy,nrot) ! eigenval of Hessian 2194 call ordr(ev,yy,3,-1) ! ordering 2195 2196 ! modification of the Newton-Raphson step to searching 2197 ! specific type of CP (Popelier algorithm) 2198 2199 ff(:)=0._dp 2200 lp(:)=0._dp 2201 dc(:)=0._dp 2202 outof=.false. 2203 do ii=1,3 2204 do jj=1,3 2205 ff(ii)=ff(ii)+yy(jj,ii)*grho(jj) 2206 end do 2207 end do 2208 id=sign(1._dp,ev(1))+sign(1._dp,ev(2))+sign(1._dp,ev(3)) 2209 if (id /= sort) then 2210 outof=.true. 2211 select case (sort) 2212 case (-1) 2213 lp(3)=0.5_dp*(ev(3)-sqrt(ev(3)*ev(3)+4._dp*ff(3)*ff(3))) 2214 pom(:,:)=0._dp 2215 pom2(:,:)=0._dp 2216 lamb(:)=0._dp 2217 do ii=1,2 2218 pom(ii,ii)=ev(ii) 2219 pom(ii,3)=ff(ii) 2220 pom(3,ii)=ff(ii) 2221 end do 2222 call jacobi(pom,3,3,lamb,pom2,nrot) 2223 call ordr(lamb,pom2,3,1) 2224 do ii=1,3 2225 lp(1)=lamb(ii) 2226 if (abs(pom2(3,ii))>1.0d-24) exit 2227 end do 2228 lp(2)=lp(1) 2229 2230 ! write(std_out,*) (ev(ii),ii=1,3) 2231 ! write(std_out,*) (lamb(ii),ii=1,3) 2232 ! write(std_out,*) ':ID ',id,lp(1),lp(3) 2233 2234 case (1) 2235 lp(1)=0.5_dp*(ev(1)+sqrt(ev(1)*ev(1)+4._dp*ff(1)*ff(1))) 2236 pom(:,:)=0._dp 2237 pom2(:,:)=0._dp 2238 lamb(:)=0._dp 2239 do ii=2,3 2240 pom(ii-1,ii-1)=ev(ii) 2241 pom(ii-1,3)=ff(ii) 2242 pom(3,ii-1)=ff(ii) 2243 end do 2244 call jacobi(pom,3,3,lamb,pom2,nrot) 2245 call ordr(lamb,pom2,3,1) 2246 do ii=3,1,-1 2247 lp(2)=lamb(ii) 2248 if (abs(pom2(3,ii))>1.0d-24) exit 2249 end do 2250 lp(3)=lp(2) 2251 2252 case (3) 2253 pom(:,:)=0._dp 2254 pom2(:,:)=0._dp 2255 lamb(:)=0._dp 2256 do ii=1,3 2257 pom(ii,ii)=ev(ii) 2258 pom(ii,4)=ff(ii) 2259 pom(4,ii)=ff(ii) 2260 end do 2261 call jacobi(pom,4,4,lamb,pom2,nrot) 2262 call ordr(lamb,pom2,4,1) 2263 do ii=4,1,-1 2264 lp(1)=lamb(ii) 2265 if (abs(pom2(4,ii))>1.0d-24) exit 2266 end do 2267 lp(2)=lp(1); lp(3)=lp(1) 2268 case default 2269 lp(:)=0._dp 2270 end select 2271 end if 2272 2273 do ii=1,3 2274 if (abs(ev(ii)-lp(ii))<1.0d-24) then 2275 outof=.false. 2276 exit 2277 end if 2278 end do 2279 do ii=1,3 ! SEARCHING STEP 2280 do jj=1,3 2281 if (outof) then 2282 dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/(ev(jj)-lp(jj)) 2283 elseif (abs(ev(jj))>1.0d-24) then 2284 dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/ev(jj) 2285 else 2286 ABI_ERROR("zero eigval of Hessian") 2287 end if 2288 end do 2289 end do 2290 2291 dltcmax = vnorm(dc,0) 2292 if (dltcmax>dmax) then ! STEP RESTRICTION 2293 do ii=1,3 2294 dc(ii)=dc(ii)*dmax/dltcmax 2295 end do 2296 end if ! primitive handling of oscillations 2297 ss=vnorm(dc,0) ! usually not needed 2298 ss=abs(ss-dv)/ss 2299 if ((ss < evol).and.(oscl)) then 2300 dc(:)=dc(:)/2._dp 2301 end if 2302 2303 2304 do ii=1,3 2305 vv(ii) = vv(ii) - dc(ii) 2306 end do 2307 2308 ! DEBUG 2309 ! write(std_out,'(":POSIN ",3F16.8)') vv 2310 ! ENDDEBUG 2311 2312 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2313 dg = vnorm(grho,0) 2314 2315 if (deb) then ! DEBUGG OUTPUT 2316 write(std_out,'("AFTER STEP ===================================")') 2317 write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((yy(ii,jii),jii=1,3),ii=1,3) 2318 write(std_out,'(":DC ",3F16.8)') dc 2319 write(std_out,*) 'STEP ',istep 2320 write(std_out,'(":POS ",3F16.8)') vv 2321 write(std_out,'(":GRAD ",3F16.8)') grho 2322 write(std_out,*) ':DGRAD,CHG ',dg,chg 2323 write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jii),jii=1,3),ii=1,3) 2324 end if 2325 vt(:)=vv(:)-vold(:) 2326 dv=vnorm(vt,0) 2327 ss=abs(dvold-dv)/dv 2328 if (ss < evol) oscl=.true. 2329 end do 2330 2331 !end of main cycle 2332 2333 !the final output 2334 2335 write(std_out,*) 'iste:',istep, dv, dg 2336 if (istep>=aim_maxstep)then 2337 write(std_out,*) ' istep=MAXSTEP ! Examine lstep2 and lgrad2 .' 2338 if ( (dv>aim_dtset%lstep2) .and. (dg>aim_dtset%lgrad2 )) then 2339 write(std_out,'(":POSOUT ",3F16.8)') vv 2340 ires=1 2341 end if 2342 end if 2343 2344 vt(:)=vv(:) 2345 2346 !write(std_out,'(":POSOUT ",3F16.8)') vv 2347 !write(std_out,'(":RBPOSOUT ",3F16.8)') vt 2348 2349 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2350 2351 !write(std_out,'(":GRAD ",3F16.8)') grho 2352 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)')& 2353 !& ((hrho(ii,jii),jii=1,3),ii=1,3) 2354 2355 2356 !FINAL INVERSION OF HESSIAN 2357 2358 call ludcmp(hrho,3,3,ipiv,id,info) 2359 if (info /= 0) then 2360 write(std_out,*) 'Error inverting hrho:' 2361 do ii=1,3 2362 write(std_out,*) (hrho(ii,jii),jii=1,3) 2363 end do 2364 ires=1 2365 ! DEBUG 2366 ! write(std_out,*)' critic : exit, ires=1' 2367 ! ENDDEBUG 2368 return 2369 ! stop 'ERROR INVERTING HESSIAN' 2370 end if 2371 do ii=1,3 2372 yy(ii,1:3)=0. 2373 yy(ii,ii)=1. 2374 end do 2375 do jii=1,3 2376 call lubksb(hrho,3,3,ipiv,yy(1,jii)) 2377 end do 2378 2379 2380 !write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((y(ii,jii),jii=1,3),ii=1,3) 2381 2382 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2383 2384 !write(std_out,'("LAPLAC:",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 2385 2386 call jacobi(hrho,3,3,ev,yy,nrot) 2387 call ordr(ev,yy,3,1) 2388 zz(:,:)=yy(:,:) 2389 2390 !do ii=1,3 2391 !do jii=1,3 2392 !zz(ii,jii)=yy(jii,ii) 2393 !end do 2394 !end do 2395 2396 !write(std_out,'(":AUTOVAL ",3F16.8)') (ev(ii),ii=1,3) 2397 !write(std_out,'(":AUTOVEC ",/,3F16.8,/,3F16.8,/,3F16.8)') ((zz(ii,jii),ii=1,3),jii=1,3) 2398 2399 if (sort/=0) then 2400 ABI_FREE(pom) 2401 ABI_FREE(pom2) 2402 ABI_FREE(lamb) 2403 end if 2404 2405 !DEBUG 2406 !write(std_out,*)' critic : exit, ires= ',ires 2407 !ENDDEBUG 2408 end subroutine critic
m_bader/critics [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
critics
FUNCTION
Search for critical points starting between
atom inxat and its neighbors.
INPUTS
aim_dtset= the structured entity containing all input variables
dstmax=maximum distance to search for neighbors
stwo, sthree, sfour: logical switches (TRUE/FALSE) indicating
to search CP starting in the middle point
of two, three or four atoms. One of these
atoms is inxat.
OUTPUT
(see side effects)
SIDE EFFECTS
This routines acts primarily on the data contained in the aim_prom module WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
2497 subroutine critics(aim_dtset,inxat,stwo,sthree,sfour,dstmax) 2498 2499 !Arguments ------------------------------------ 2500 !scalars 2501 integer,intent(in) :: inxat 2502 real(dp),intent(in) :: dstmax 2503 logical,intent(in) :: sfour,sthree,stwo 2504 !no_abirules 2505 type(aim_dataset_type), intent(in) :: aim_dtset 2506 2507 !Local variables ------------------------------ 2508 !scalars 2509 integer :: i1,i2,i3,iat,ii,ipos,ires,jii,jj,kjj,kk,ll,n1,n2,n3,nb,nc 2510 integer :: nshell 2511 real(dp) :: chg,dif1,dif2,diff,dist,olddist,rr 2512 ! real(dp) :: ss,uu 2513 logical :: found,inter 2514 !arrays 2515 integer :: ibat(nnpos*natom),inat(nnpos*natom),ipibat(nnpos*natom) 2516 integer :: nnat(nnpos*natom),nr(nnpos*natom) 2517 real(dp) :: dif(3),dists(nnpos*natom),ev(3),grho(3),hrho(3,3) 2518 real(dp) :: pom(3),v1(3),v2(3),v3(3),v4(3),vi(3),vt(3),zz(3,3) 2519 2520 !************************************************************************ 2521 vi(:)=xatm(:,inxat) 2522 2523 nc=0 2524 do jii=1,nnpos 2525 do kjj=1,natom 2526 dist=0._dp 2527 dif(:)=xatm(:,inxat)-xatm(:,kjj)-atp(:,jii) 2528 2529 ! do ii=1,3 2530 ! dif(ii)=xatm(ii,inxat)-xatm(ii,kjj)-atp(ii,jii) 2531 ! end do 2532 dist=vnorm(dif,0) 2533 if (.not.((dist>dstmax).or.(dist<0.001))) then 2534 nc=nc+1 2535 dists(nc)=dist 2536 nnat(nc)=kjj 2537 inat(nc)=jii 2538 end if 2539 end do 2540 end do 2541 do n1=1,nc 2542 nr(n1)=n1 2543 end do 2544 call sort_dp(nc,dists,nr,tol14) 2545 nb=0 2546 olddist=0._dp 2547 nshell=0 2548 !write(std_out,*) ':ORIAT ', (xatm(ii,inxat),ii=1,3) 2549 do n1=1,nc 2550 n2=nr(n1) 2551 n3=nnat(n2) 2552 if (dists(n1)<(2*dists(1))) then 2553 if ((dists(n1)-olddist)>aim_dlimit) then 2554 nshell=nshell+1 2555 olddist=dists(n1) 2556 if (nshell==5) exit 2557 end if 2558 nb=nb+1 2559 ibat(nb)=n3 2560 ipibat(nb)=inat(n2) 2561 write(std_out,*) ':NEIG ',inxat,n3,inat(n2),dists(n1) 2562 ! write(std_out,*) ':POSAT',(xatm(ii,ibat(nb))+atp(ii,ipibat(nb)),ii=1,3) 2563 else 2564 exit 2565 end if 2566 end do 2567 2568 npc=0 2569 npcm3=0 2570 2571 ! 2572 !.....SEARCH BETWEEN EACH PAIR OF ATOMS 2573 ! 2574 2575 if (stwo) then 2576 do jii=1,nb 2577 do ii=1,3 2578 v1(ii)=xatm(ii,inxat) 2579 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2580 vt(ii)=(v1(ii)+v2(ii))/2._dp 2581 end do 2582 inter=.true. 2583 diff=0._dp 2584 pom(:)=vt(:) 2585 pom(:)=pom(:)-vi(:) 2586 diff=vnorm(pom,0) 2587 if (diff > maxcpdst) inter=.false. 2588 if (inter) then 2589 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2590 if (ires==0) then 2591 found=.false. 2592 if (npc > 0) then 2593 do jj=1,npc 2594 pom(:)=vt(:)-pc(:,jj) 2595 dist=vnorm(pom,0) 2596 if (dist < aim_dtset%dpclim) found=.true. 2597 end do 2598 end if 2599 if (.not.found) then 2600 pom(:)=vt(:) 2601 call bschg1(pom,-1) 2602 pcrb(:,npc+1)=pom(:) 2603 pom(:)=pom(:)-vi(:) 2604 diff=vnorm(pom,0) 2605 if (abs(diff) > maxcpdst) found=.true. 2606 end if 2607 if (.not.found) then 2608 npc=npc+1 2609 do jj=1,3 2610 pc(jj,npc)=vt(jj) 2611 evpc(jj,npc)=ev(jj) 2612 do kk=1,3 2613 zpc(kk,jj,npc)=zz(kk,jj) 2614 end do 2615 end do 2616 i1=ev(1)/abs(ev(1)) 2617 i2=ev(2)/abs(ev(2)) 2618 i3=ev(3)/abs(ev(3)) 2619 icpc(npc)=i1+i2+i3 2620 if (icpc(npc)==-3) then 2621 npcm3=npcm3+1 2622 end if 2623 write(std_out,*) 'New critical point found' 2624 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2625 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2626 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2627 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2628 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2629 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2630 write(22,'(":PC2",3F10.6,3E12.4,I4,2E12.4)') & 2631 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2632 write(std_out,'(":PC2",3F10.6,3E12.4,I4,2E12.4)') & 2633 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2634 pom(:)=vt(:)-v1(:) 2635 dif1=vnorm(pom,0) 2636 pom(:)=vt(:)-v2(:) 2637 dif2=vnorm(pom,0) 2638 write(std_out,'(":DISPC ",2F12.8)') dif1,dif2 2639 end if 2640 end if 2641 end if 2642 end do 2643 end if 2644 ! 2645 !.....SEARCH BETWEEN EACH THREE ATOMS 2646 ! 2647 if(sthree) then 2648 do jii=1,nb 2649 do kjj=jii+1,nb 2650 do ii=1,3 2651 v1(ii)=xatm(ii,inxat) 2652 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2653 v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj)) 2654 vt(ii)=(v1(ii)+v2(ii)+v3(ii))/3._dp 2655 end do 2656 inter=.true. 2657 pom(:)=vt(:) 2658 pom(:)=pom(:)-vi(:) 2659 dist=vnorm(pom,0) 2660 if (abs(diff)>maxcpdst) then 2661 inter=.false. 2662 exit 2663 end if 2664 if (inter) then 2665 do jj=1,npc 2666 pom(:)=pc(:,jj)-vt(:) 2667 diff=vnorm(pom,0) 2668 if (diff<aim_dpc0) then 2669 inter=.false. 2670 exit 2671 end if 2672 end do 2673 end if 2674 if (inter) then 2675 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2676 if (ires==0) then 2677 found=.false. 2678 if (npc>0) then 2679 do jj=1,npc 2680 pom(:)=vt(:)-pc(:,jj) 2681 dist=vnorm(pom,0) 2682 if (dist<aim_dtset%dpclim) then 2683 found=.true. 2684 exit 2685 end if 2686 end do 2687 end if 2688 if (.not.found) then 2689 pom(:)=vt(:) 2690 call bschg1(pom,-1) 2691 pcrb(:,npc+1)=pom(:) 2692 pom(:)=pom(:)-vi(:) 2693 diff=vnorm(pom,0) 2694 if (abs(diff)>maxcpdst) found=.true. 2695 end if 2696 if (.not.found) then 2697 npc=npc+1 2698 do jj=1,3 2699 pc(jj,npc)=vt(jj) 2700 evpc(jj,npc)=ev(jj) 2701 do kk=1,3 2702 zpc(kk,jj,npc)=zz(kk,jj) 2703 end do 2704 end do 2705 i1=ev(1)/abs(ev(1)) 2706 i2=ev(2)/abs(ev(2)) 2707 i3=ev(3)/abs(ev(3)) 2708 icpc(npc)=i1+i2+i3 2709 if (icpc(npc)==-3) then 2710 npcm3=npcm3+1 2711 end if 2712 write(std_out,*) 'New critical point found' 2713 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2714 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2715 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2716 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2717 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2718 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2719 write(22,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') & 2720 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2721 write(std_out,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') & 2722 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2723 end if 2724 end if 2725 end if 2726 end do 2727 end do 2728 end if 2729 2730 ! 2731 !.....SEARCH BETWEEN EACH FOUR ATOMS 2732 ! 2733 if (sfour) then 2734 do jii=1,nb 2735 do kjj=jii+1,nb 2736 do ll=jii+1,nb 2737 do ii=1,3 2738 v1(ii)=xatm(ii,inxat) 2739 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2740 v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj)) 2741 v4(ii)=xatm(ii,ibat(ll))+atp(ii,ipibat(ll)) 2742 vt(ii)=(v1(ii)+v2(ii)+v3(ii)+v4(ii))/4._dp 2743 end do 2744 inter=.true. 2745 pom(:)=vt(:) 2746 pom(:)=pom(:)-vi(:) 2747 diff=vnorm(pom,0) 2748 if (abs(diff)>maxcpdst) then 2749 inter=.false. 2750 exit 2751 end if 2752 if (inter) then 2753 do jj=1,npc 2754 pom(:)=pc(:,jj)-vt(:) 2755 diff=vnorm(pom,0) 2756 if (diff < aim_dpc0) then 2757 inter=.false. 2758 exit 2759 end if 2760 end do 2761 end if 2762 if (inter) then 2763 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2764 if (ires==0) then 2765 found=.false. 2766 if (npc>0) then 2767 do jj=1,npc 2768 pom(:)=vt(:)-pc(:,jj) 2769 dist=vnorm(pom,0) 2770 if (dist < aim_dtset%dpclim) found=.true. 2771 end do 2772 end if 2773 if (.not.found) then 2774 pom(:)=vt(:) 2775 pcrb(:,npc+1)=pom(:) 2776 pom(:)=pom(:)-vi(:) 2777 diff=vnorm(pom,0) 2778 if (abs(diff)>maxcpdst) found=.true. 2779 end if 2780 if (.not.found) then 2781 npc=npc+1 2782 do jj=1,3 2783 pc(jj,npc)=vt(jj) 2784 evpc(jj,npc)=ev(jj) 2785 do kk=1,3 2786 zpc(kk,jj,npc)=zz(kk,jj) 2787 end do 2788 end do 2789 i1=ev(1)/abs(ev(1)) 2790 i2=ev(2)/abs(ev(2)) 2791 i3=ev(3)/abs(ev(3)) 2792 icpc(npc)=i1+i2+i3 2793 if (icpc(npc)==-3) then 2794 npcm3=npcm3+1 2795 end if 2796 write(std_out,*) 'New critical point found' 2797 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2798 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2799 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2800 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2801 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2802 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2803 write(22,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') & 2804 & (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2805 write(std_out,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') & 2806 & (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2807 end if 2808 end if 2809 end if 2810 end do 2811 end do 2812 end do 2813 end if 2814 2815 write(std_out,*) npc 2816 end subroutine critics
m_bader/defad [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
defad
FUNCTION
Initialisation of aim input variables to their default values.
INPUTS
(no input : initialisation by default values)
OUTPUT
aim_dtset = the structured entity containing all input variables
SOURCE
2834 subroutine defad(aim_dtset) 2835 2836 !Arguments ------------------------------------ 2837 !scalars 2838 type(aim_dataset_type),intent(out) :: aim_dtset 2839 2840 !Local variables ------------------------------ 2841 2842 ! ********************************************************************* 2843 2844 aim_dtset%isurf=0 2845 aim_dtset%crit=0 2846 aim_dtset%irsur=0 2847 aim_dtset%foll=0 2848 aim_dtset%irho=0 2849 aim_dtset%ivol=0 2850 aim_dtset%denout=0 2851 aim_dtset%lapout=0 2852 aim_dtset%gpsurf=0 2853 aim_dtset%plden=0 2854 aim_dtset%dltyp=0 2855 2856 aim_dtset%batom=1 2857 aim_dtset%nsa=3 2858 aim_dtset%nsb=3 2859 aim_dtset%nsc=3 2860 aim_dtset%npt=100 2861 aim_dtset%nth=32 2862 aim_dtset%nph=48 2863 2864 aim_dtset%themax=pi 2865 aim_dtset%themin=zero 2866 aim_dtset%phimin=zero 2867 aim_dtset%phimax=two_pi 2868 aim_dtset%phi0=zero 2869 aim_dtset%th0=zero 2870 aim_dtset%folstp=5.d-2 2871 aim_dtset%dr0=5.d-2 2872 aim_dtset%atrad=one 2873 aim_dtset%rmin=one 2874 2875 aim_dtset%foldep(:)=zero 2876 aim_dtset%vpts(:,:)=zero 2877 aim_dtset%ngrid(:)=30 2878 aim_dtset%scal(:)=one 2879 aim_dtset%maxatd=1.d1 2880 aim_dtset%maxcpd=5.d1 2881 2882 aim_dtset%dpclim=1.d-2 2883 aim_dtset%lstep=1.d-10 2884 aim_dtset%lstep2=1.d-5 2885 aim_dtset%lgrad=1.d-12 2886 aim_dtset%lgrad2=1.d-5 2887 aim_dtset%coff1=0.98_dp 2888 aim_dtset%coff2=0.95_dp 2889 2890 end subroutine defad
m_bader/drvaim [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
drvaim
FUNCTION
Main driver for the Bader analysis it looks the values of the input variables and calls corresponding procedures
INPUTS
aim_dtset = the structured entity containing all input variables tcpui=initial CPU time twalli=initial wall clock time
OUTPUT
(see side effects)
SIDE EFFECTS
this routine acts primarily on the data contained in the aimprom module WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
2918 subroutine drvaim(aim_dtset,tcpui,twalli) 2919 2920 !Arguments ------------------------------------ 2921 !scalars 2922 real(dp) :: tcpui,twalli 2923 type(aim_dataset_type),intent(in) :: aim_dtset 2924 2925 !Local variables ------------------------------ 2926 !scalars 2927 integer :: iat,iatinit,inxat,ipos,iposinit 2928 integer :: me,npmax,nproc,nstep 2929 real(dp) :: dstlim,rr,ss,t1,t2,tf,wall 2930 real(dp) :: tcpu,twall,znucl_batom 2931 logical :: debold,sfour,srch,sthree,stwo 2932 !arrays 2933 real(dp) :: tsec(2) 2934 real(dp) :: grho(3),xstart(3) 2935 2936 ! ********************************************************************* 2937 2938 me=xmpi_comm_rank(xmpi_world) 2939 nproc=xmpi_comm_size(xmpi_world) 2940 2941 !These input variables might be modified during what follows, 2942 !so, they are copied outside of aim_dtset. 2943 inxat=aim_dtset%batom 2944 r0=aim_dtset%atrad 2945 h0=aim_dtset%folstp 2946 maxatdst=aim_dtset%maxatd 2947 maxcpdst=aim_dtset%maxcpd 2948 2949 dstlim=maxcpdst 2950 2951 !Flags from the old version 2952 !to be remove later 2953 deb=.false. 2954 stwo=.true. 2955 sthree=.true. 2956 sfour=.false. 2957 srch=.false. 2958 2959 npmax=aim_npmaxin 2960 2961 !Main initialisation procedure - 2962 !- it reads ABINIT density file and files 2963 !with core densities and initialises the fields for 2964 !spline interpolation 2965 2966 call initaim(aim_dtset,znucl_batom) 2967 2968 2969 !CP SEARCHING 2970 2971 if (aim_dtset%crit /= 0) then 2972 2973 call timein(tcpu,twall) 2974 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 2975 write(std_out, '(5a,f13.1,a,f13.1)' ) & 2976 & '-',ch10,'- Before searching the CP ',ch10,& 2977 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 2978 2979 if (aim_dtset%crit==3) then 2980 ! old version of the driver for searching CPs (original code) 2981 call critics(aim_dtset,inxat,stwo,sthree,sfour,dstlim) 2982 else 2983 ! driver for searching CPs with Popellier algorithm 2984 call cpdrv(aim_dtset) 2985 end if 2986 2987 call timein(tcpu,twall) 2988 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 2989 write(std_out, '(5a,f13.1,a,f13.1)' ) & 2990 & '-',ch10,'- After searching the CP ',ch10,& 2991 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 2992 2993 end if 2994 2995 ! 2996 !BADER SURFACE CALCULATION 2997 ! 2998 2999 if (aim_dtset%isurf==1) then 3000 ! driver for determination of the Bader surface 3001 3002 call timein(tcpu,twall) 3003 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3004 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3005 & '-',ch10,'- Before determinating the Bader surface ',ch10,& 3006 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3007 3008 call surf(aim_dtset) 3009 3010 call timein(tcpu,twall) 3011 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3012 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3013 & '-',ch10,'- After determinating the Bader surface ',ch10,& 3014 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3015 3016 end if 3017 3018 ! 3019 !CHARGE INTEGRATIOM 3020 ! 3021 3022 if (aim_dtset%irho==1) then 3023 call integrho(aim_dtset,znucl_batom) 3024 end if 3025 3026 ! 3027 !VOLUME INTEGRATION OF THE BADER ATOM 3028 ! 3029 3030 if (aim_dtset%ivol==1) then 3031 call integvol() 3032 end if 3033 3034 ! 3035 !ONE RADIUS OF THE BADER SURFACE 3036 ! 3037 3038 if (aim_dtset%irsur==1) then 3039 if (aim_dtset%isurf/=0) srch=.true. 3040 iat=aim_dtset%batom 3041 ss=r0 3042 call timein(t1,wall) 3043 call rsurf(aim_dtset,rr,grho,aim_dtset%th0,aim_dtset%phi0,ss,iat,npmax,srch) 3044 call timein(t2,wall) 3045 t2=t2-t1 3046 write(unts,'(2F12.8,F15.10)') aim_dtset%th0,aim_dtset%phi0,rr 3047 write(std_out,'(":RSUR ",2F12.8,2F15.10)') aim_dtset%th0,aim_dtset%phi0,rr,t2 3048 end if 3049 3050 ! 3051 !FOLLOW THE GRADIENT PATH FROM ONE POINT 3052 ! 3053 3054 if (aim_dtset%foll==1) then 3055 iatinit=aim_dtset%batom 3056 iposinit=batcell 3057 if (aim_dtset%isurf/=0) srch=.true. 3058 debold=deb 3059 xstart(:)=aim_dtset%foldep(:) 3060 call timein(t1,wall) 3061 call aim_follow(aim_dtset,xstart,npmax,srch,iatinit,iposinit,iat,ipos,nstep) 3062 call timein(t2,wall) 3063 tf=t2-t1 3064 write(std_out,'(":TIME in aim_follow:", F12.4)') tf 3065 end if 3066 3067 if (aim_dtset%plden == 1) then 3068 ! profile of the density integrated in plane xy 3069 ! belong the z-axes - not finished - cut3d better ! 3070 call plint() 3071 end if 3072 3073 if ((aim_dtset%denout > 0).or.(aim_dtset%lapout > 0)) then 3074 ! additional outputs of density and laplacian fields 3075 ! in the plane or line 3076 call addout(aim_dtset) 3077 end if 3078 3079 if (aim_dtset%gpsurf == 1) then 3080 ! script for gnuplot - simple demonstration of the 3081 ! computed surface 3082 call graph(unts,untg) 3083 end if 3084 3085 !Deallocation of global variables allocated in initaim 3086 !and declared in defs_aimfields. 3087 ABI_FREE(dig1) 3088 ABI_FREE(dig2) 3089 ABI_FREE(dig3) 3090 ABI_FREE(llg1) 3091 ABI_FREE(llg2) 3092 ABI_FREE(llg3) 3093 ABI_FREE(cdig1) 3094 ABI_FREE(cdig2) 3095 ABI_FREE(cdig3) 3096 ABI_FREE(ddx) 3097 ABI_FREE(ddy) 3098 ABI_FREE(ddz) 3099 ABI_FREE(rrad) 3100 ABI_FREE(crho) 3101 ABI_FREE(sp2) 3102 ABI_FREE(sp3) 3103 ABI_FREE(sp4) 3104 ABI_FREE(corlim) 3105 ABI_FREE(dvl) 3106 ABI_FREE(ndat) 3107 ABI_FREE(rminl) 3108 !Deallocation of global variables allocated in initaim 3109 !and declared in defs_aimprom. 3110 ABI_FREE(typat) 3111 ABI_FREE(xred) 3112 ABI_FREE(xatm) 3113 3114 end subroutine drvaim
m_bader/graph [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
graph
FUNCTION
Writing of the gnuplot script to show the computed part of Bader surface with lines
INPUTS
untg = unit number of the file on which the info is written unts = unit number of the file from which the Bader surface is read
OUTPUT
(written in the untg file)
SOURCE
3134 subroutine graph(unts,untg) 3135 3136 !Arguments ------------------------------------ 3137 !scalars 3138 integer,intent(in) :: untg,unts 3139 3140 !Local variables ------------------------------ 3141 !scalars 3142 integer :: ii,indx,jj,nphi,nth 3143 real(dp),parameter :: snull=1.d-6 3144 real(dp) :: phimax,phimin,ss,thmax,thmin 3145 !arrays 3146 real(dp) :: xorig(3) 3147 real(dp),allocatable :: phi(:),rr(:,:),th(:) 3148 3149 ! ********************************************************************* 3150 3151 rewind(unts) 3152 read(unts,*) indx, xorig(1:3) 3153 read(unts,*) nth, thmin, thmax 3154 read(unts,*) nphi, phimin, phimax 3155 ABI_MALLOC(th,(nth)) 3156 ABI_MALLOC(phi,(nphi)) 3157 ABI_MALLOC(rr,(nth,nphi)) 3158 do ii=1,nth 3159 do jj=1,nphi 3160 read(unts,*) th(ii),phi(jj),rr(ii,jj),ss 3161 end do 3162 end do 3163 3164 !end of reading 3165 3166 write(untg,*) 'reset' 3167 write(untg,*) 'set st d l' 3168 write(untg,*) 'set ticslevel 0' 3169 write(untg,*) 'set title ''Bader surface'' ' 3170 write(untg,*) 'splot ''-'' using ($3*sin($1)*cos($2)):($3*sin($1)*sin($2)):($3*cos($1)) notitle' 3171 do ii=1,nth 3172 do jj=1,nphi 3173 write(untg,'(2F12.8,E16.8)') th(ii),phi(jj),rr(ii,jj) 3174 end do 3175 if ((ii==nth).and.(jj==nphi)) then 3176 cycle 3177 else 3178 write(untg,*) 3179 end if 3180 end do 3181 3182 end subroutine graph
m_bader/initaim [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
initaim
FUNCTION
Initialization for the 3D interpolation for the AIM code:
- this procedure reads the charge density of the electrons of valence on
the equidistant 3D grid (*_DEN output file of ABINIT) and the core charge
density of electrons from *.fc files (fhi package)
- the Cholesky decomposition of the general matrix for
the computation of the 1D spline coeficients in each direction is done.
Warning - the procedure is modified to use periodic boundary conditions
already during the decomposition
- the second derivations of valence density in three directions are computed
and stored in the real space grid of the density for interpolation.
- the core density is stored separately in the radial grid together with th
second radial derivation
INPUTS
aim_dtset= the structured entity containing all input variables
OUTPUT
znucl_batom= the nuclear charge of the Bader atom (see side effects)
SIDE EFFECTS
thie routine works on the data contained in the aim_fields and aim_prom modules WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
3218 subroutine initaim(aim_dtset,znucl_batom) 3219 3220 !Arguments ------------------------------------ 3221 !scalars 3222 type(aim_dataset_type),intent(in) :: aim_dtset 3223 3224 !Local variables ------------------------------ 3225 !scalars 3226 integer,parameter :: master=0 3227 integer :: fform0,id,ierr,ii,info,jj,kk,kod,mm,ndtmax,nn,nsa,nsb,nsc,nsym,me,nproc,npsp 3228 integer :: unth,comm 3229 #ifdef HAVE_NETCDF 3230 integer :: den_id 3231 #endif 3232 real(dp) :: ss,ucvol,znucl_batom 3233 real(dp) :: zz 3234 type(hdr_type) :: hdr 3235 !arrays 3236 integer :: ipiv(3) 3237 integer,allocatable :: symrel(:,:,:) 3238 real(dp) :: aa(3),bb(3),gmet(3,3),gprimd(3,3),rmet(3,3),yy(3,3) 3239 real(dp),allocatable :: tnons(:,:),znucl(:),zionpsp(:) 3240 real(dp),pointer :: ptc(:),ptd(:),ptf(:),ptp(:),ptsd(:) 3241 3242 ! ********************************************************************* 3243 3244 !DEBUG 3245 !write(std_out,*) ' initaim : enter ' 3246 !ENDDEBUG 3247 3248 comm = xmpi_world 3249 me=xmpi_comm_rank(comm) 3250 nproc=xmpi_comm_size(comm) 3251 3252 slc=0 ! code for follow 3253 3254 !The use of the "hdr" routines is much better for the future 3255 !maintenance of the code. Indeed, the content of the header 3256 !will continue to change from time to time, and the associated 3257 !changes should be done in one place only. 3258 3259 !Read ABINIT header ---------------------------------------------------------- 3260 if(me==master)then 3261 if (aim_iomode == IO_MODE_ETSF) then 3262 call hdr_ncread(hdr, untad, fform0) 3263 else 3264 call hdr_fort_read(hdr, untad, fform0) 3265 end if 3266 end if 3267 ABI_CHECK(fform0 /= 0, "hdr_read returned fform == 0") 3268 call hdr%bcast(master, me, comm) 3269 3270 !Echo part of the header 3271 call hdr%echo(fform0, 4, unit=std_out) 3272 call hdr%echo(fform0, 4, unit=untout) 3273 3274 natom=hdr%natom 3275 ngfft(1:3)=hdr%ngfft(:) 3276 nsym=hdr%nsym 3277 npsp=hdr%npsp 3278 ntypat=hdr%ntypat 3279 rprimd(:,:)=hdr%rprimd(:,:) 3280 3281 ABI_MALLOC(zionpsp,(npsp)) 3282 ABI_MALLOC(znucl,(ntypat)) 3283 ABI_MALLOC(typat,(natom)) 3284 ABI_MALLOC(xred,(3,natom)) 3285 ABI_MALLOC(symrel,(3,3,nsym)) 3286 ABI_MALLOC(tnons,(3,nsym)) 3287 ABI_MALLOC(xatm,(3,natom)) 3288 3289 symrel(:,:,:)=hdr%symrel(:,:,:) 3290 typat(:)=hdr%typat(:) 3291 tnons(:,:)=hdr%tnons(:,:) 3292 znucl(:)=hdr%znucltypat(:) 3293 zionpsp(:)=hdr%zionpsp(:) 3294 xred(:,:)=hdr%xred(:,:) 3295 3296 call hdr%free() 3297 3298 !------------------------------------------------------------------------------- 3299 3300 ABI_MALLOC(dvl,(ngfft(1),ngfft(2),ngfft(3))) 3301 3302 if(me==master)then 3303 if (aim_iomode == IO_MODE_ETSF) then 3304 #ifdef HAVE_NETCDF 3305 ! netcdf array has shape [cplex, n1, n2, n3, nspden]), here we read only the total density. 3306 NCF_CHECK(nf90_inq_varid(untad, "density", den_id)) 3307 NCF_CHECK(nf90_get_var(untad, den_id, dvl, start=[1,1,1,1], count=[1, ngfft(1), ngfft(2), ngfft(3), 1])) 3308 #endif 3309 else 3310 read(untad,iostat=nn) dvl(1:ngfft(1),1:ngfft(2),1:ngfft(3)) 3311 ABI_CHECK(nn==0,"error of reading !") 3312 end if 3313 end if 3314 call xmpi_bcast(dvl, master, comm, ierr) 3315 3316 write(std_out,*)ch10,' initaim : the valence density has been read' ,ch10 3317 3318 !INITIALISATION OF SOME IMPORTANT FIELDS 3319 3320 !Only interpolation is computed (inside vgh_rho) in reduced 3321 !coordinates. In all other routines the cart. coordinates (CC) are used. 3322 3323 !transformation of the atom positions to CC 3324 do ii=1,natom 3325 xatm(:,ii)=xred(:,ii) 3326 call bschg1(xatm(:,ii),1) 3327 end do 3328 3329 !Generation of the neighbouring cells + transf to CC 3330 nn=0 3331 nsa=aim_dtset%nsa ; nsb=aim_dtset%nsb ; nsc=aim_dtset%nsc 3332 do ii=-nsa,nsa 3333 do jj=-nsb,nsb 3334 do kk=-nsc,nsc 3335 nn=nn+1 3336 atp(1,nn)=ii*1._dp 3337 atp(2,nn)=jj*1._dp 3338 atp(3,nn)=kk*1._dp 3339 call bschg1(atp(:,nn),1) 3340 end do 3341 end do 3342 end do 3343 nnpos=nn 3344 3345 !DEBUG 3346 !write(std_out,*)' initaim : nnpos=',nnpos 3347 !ENDDEBUG 3348 3349 batcell=nsa*(2*nsb+1)*(2*nsc+1)+(2*nsc+1)*nsb+nsc+1 3350 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 3351 maxatdst=min(maxatdst, nsa*sqrt(rmet(1,1)), nsb*sqrt(rmet(2,2)), nsc*sqrt(rmet(3,3)) ) 3352 if (maxcpdst > maxatdst) maxcpdst=0.75*maxatdst 3353 3354 3355 !RPRIM ITS INVERSE AND TRANSPOSE 3356 3357 do ii=1,3 3358 do jj=1,3 3359 yy(ii,jj)=rprimd(ii,jj) 3360 end do 3361 end do 3362 call ludcmp(yy,3,3,ipiv,id,info) 3363 ABI_CHECK(info==0,'Error inverting rprimd') 3364 3365 do ii=1,3 3366 do jj=1,3 3367 ivrprim(ii,jj)=0._dp 3368 end do 3369 ivrprim(ii,ii)=1._dp 3370 end do 3371 do ii=1,3 3372 call lubksb(yy,3,3,ipiv,ivrprim(:,ii)) 3373 end do 3374 do ii=1,3 3375 do jj=1,3 3376 trivrp(ii,jj)=ivrprim(jj,ii) 3377 end do 3378 end do 3379 3380 write(std_out,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 3381 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 3382 write(untout,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 3383 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 3384 3385 write(std_out,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 3386 write(std_out,*) "=======================================" 3387 do ii=1,natom 3388 jj=typat(ii) 3389 write(std_out,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 3390 end do 3391 write(untout,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 3392 write(untout,*) "=======================================" 3393 do ii=1,natom 3394 jj=typat(ii) 3395 write(untout,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 3396 end do 3397 3398 !STEPS IN REAL SPACE GRID (REDUCED) 3399 do ii=1,3 3400 dix(ii)=1._dp/ngfft(ii) 3401 end do 3402 3403 !READING OF THE CORE DENSITY 3404 write(std_out,*)ch10,' initaim : will read the core densities' ,ch10 3405 3406 ABI_MALLOC(ndat,(ntypat)) 3407 ABI_MALLOC(rminl,(natom)) 3408 ndtmax=0 3409 if(me==master)then 3410 do ii=1,ntypat 3411 unth=unt+ii 3412 ! DEBUG 3413 ! write(std_out,*)' read from unit ',unth 3414 ! call flush(std_out) 3415 ! stop 3416 ! ENDDEBUG 3417 read(unth,*) ndat(ii),ss 3418 if (ndat(ii)>ndtmax) ndtmax=ndat(ii) 3419 end do 3420 end if 3421 call xmpi_bcast(ndat,master,comm,ierr) 3422 call xmpi_bcast(ndtmax,master,comm,ierr) 3423 call xmpi_bcast(ss,master,comm,ierr) 3424 3425 !FIELDS FOR STORING CORE DENSITY 3426 3427 ABI_MALLOC(rrad,(ndtmax,ntypat)) 3428 ABI_MALLOC(crho,(ndtmax,ntypat)) 3429 ABI_MALLOC(sp2,(ndtmax,ntypat)) 3430 ABI_MALLOC(sp3,(ndtmax,ntypat)) 3431 ABI_MALLOC(sp4,(ndtmax,ntypat)) 3432 ABI_MALLOC(corlim,(ntypat)) 3433 3434 sp2(:,:)=zero 3435 sp3(:,:)=zero 3436 sp4(:,:)=zero 3437 3438 !Reading of the core densities 3439 corlim(:)=0 3440 kod=0 3441 if(me==master)then 3442 do ii=1,ntypat 3443 unth=unt+ii 3444 do jj=1,ndat(ii) 3445 read(unth,*) rrad(jj,ii),crho(jj,ii),sp2(jj,ii),sp3(jj,ii) 3446 ! this is the integral of the core charge read in 3447 crho(jj,ii) = crho(jj,ii)/4._dp/pi 3448 if ((crho(jj,ii) < aim_rhocormin) .and. (corlim(ii)==0)) corlim(ii)=jj 3449 sp2(jj,ii)=sp2(jj,ii)/4._dp/pi 3450 sp3(jj,ii)=sp3(jj,ii)/4._dp/pi ! ATENTION!!! in sp3 is just second derivation 3451 end do 3452 do jj=1,ndat(ii)-1 3453 sp4(jj,ii)=(sp3(jj+1,ii)-sp3(jj,ii))/(6._dp*(rrad(jj+1,ii)-rrad(jj,ii))) 3454 end do 3455 ! 3456 zz = crho(1,ii) * rrad(1,ii)**2 * (rrad(2,ii)-rrad(1,ii)) 3457 do jj=2,ndat(ii)-1 3458 zz = zz + crho(jj,ii) * rrad(jj,ii)**2 * (rrad(jj+1,ii)-rrad(jj-1,ii)) 3459 end do 3460 zz = zz * half * 4._dp * pi 3461 if (corlim(ii)==0) corlim(ii)=ndat(ii) 3462 3463 ! add check on zion wrt FHI .fc file 3464 ! compare zion to zionpsp(typat(aim_dtset%batom)) 3465 if (abs(znucl(ii) - zz - zionpsp(ii)) > 1.e-1_dp) then 3466 write (std_out,*) 'error: your core charge ', zz, ' does not correspond to the correct number' 3467 write (std_out,*) ' of valence electrons', zionpsp(ii), ' and the nuclear charge ', znucl(ii) 3468 write (std_out,*) ' You have probably used a pseudopotential which has more valence electrons than the' 3469 write (std_out,*) ' original FHI ones. ACTION: make a .fc file with the correct core charge' 3470 stop 3471 end if 3472 3473 end do 3474 end if 3475 call xmpi_bcast(rrad,master,comm,ierr) 3476 call xmpi_bcast(crho,master,comm,ierr) 3477 call xmpi_bcast(sp2,master,comm,ierr) 3478 call xmpi_bcast(sp3,master,comm,ierr) 3479 call xmpi_bcast(sp4,master,comm,ierr) 3480 call xmpi_bcast(corlim,master,comm,ierr) 3481 3482 write(std_out,*)ch10,' initaim : the core densities have been read' ,ch10 3483 3484 3485 !CORRECTION OF THE CORE DENSITY NORMALISATION 3486 crho(:,:)=1.0003*crho(:,:) 3487 sp2(:,:)=1.0003*sp2(:,:) 3488 sp3(:,:)=1.0003*sp3(:,:) 3489 sp4(:,:)=1.0003*sp4(:,:) 3490 3491 !FIELDS FOR INTERPOLATIONS OF THE VALENCE DENSITY 3492 3493 ABI_MALLOC(dig1,(ngfft(1))) 3494 ABI_MALLOC(dig2,(ngfft(2))) 3495 ABI_MALLOC(dig3,(ngfft(3))) 3496 ABI_MALLOC(llg1,(ngfft(1))) 3497 ABI_MALLOC(llg2,(ngfft(2))) 3498 ABI_MALLOC(llg3,(ngfft(3))) 3499 ABI_MALLOC(cdig1,(ngfft(1)-1)) 3500 ABI_MALLOC(cdig2,(ngfft(2)-1)) 3501 ABI_MALLOC(cdig3,(ngfft(3)-1)) 3502 ABI_MALLOC(ddx,(ngfft(1),ngfft(2),ngfft(3))) 3503 ABI_MALLOC(ddy,(ngfft(1),ngfft(2),ngfft(3))) 3504 ABI_MALLOC(ddz,(ngfft(1),ngfft(2),ngfft(3))) 3505 3506 !DECOMPOSITION OF THE MATRIX FOR THE DETERMINATION OF COEFFICIENTS 3507 !FOR CUBIC SPLINE INTERPOLATION (using the periodic boundary conditions) 3508 3509 !MAIN DIAGONAL (aa) AND SECONDARY DIAGONAL (bb) MATRIX ELEMENTS 3510 3511 nmax=ngfft(1) 3512 do ii=2,3 3513 if (ngfft(ii) > nmax) nmax=ngfft(ii) 3514 end do 3515 nullify(ptf,ptsd) 3516 nullify(ptd,ptc,ptp) 3517 aa(:)=2.0*dix(:)**2/3.0 3518 bb(:)=dix(:)**2/6.0 3519 3520 do ii=1,3 3521 if(ii==1) then 3522 ptd=>dig1;ptc=>cdig1;ptp=>llg1 3523 elseif (ii==2) then 3524 ptd=>dig2;ptc=>cdig2;ptp=>llg2 3525 else 3526 ptd=>dig3;ptc=>cdig3;ptp=>llg3 3527 end if 3528 ptd(1)=sqrt(aa(ii)) 3529 ptc(1)=bb(ii)/ptd(1) 3530 ptp(1)=ptc(1) 3531 do jj=2,ngfft(ii)-1 3532 ptd(jj)=aa(ii)-ptc(jj-1)**2 3533 if(ptd(jj)<zero) then 3534 ABI_ERROR('Matrix is not positive definite !') 3535 end if 3536 ptd(jj)=sqrt(ptd(jj)) 3537 if (jj==ngfft(ii)-1) then 3538 ptc(jj)=(bb(ii)-ptp(jj-1)*ptc(jj-1))/ptd(jj) 3539 ptp(jj)=ptc(jj) 3540 exit 3541 end if 3542 ptc(jj)=bb(ii)/ptd(jj) 3543 ptp(jj)=-ptp(jj-1)*ptc(jj-1)/ptd(jj) 3544 end do 3545 ss=0._dp 3546 do jj=1,ngfft(ii)-1 3547 ss=ss+ptp(jj)**2 3548 end do 3549 ss=aa(ii)-ss 3550 if(ss<zero) then 3551 ABI_ERROR('Matrix is not positive definite !') 3552 end if 3553 ptd(ngfft(ii))=sqrt(ss) 3554 ptp(ngfft(ii))=ptd(ngfft(ii)) 3555 3556 3557 ! INITIALISATION OF THE SECOND DERIVATIVE FIELDS 3558 3559 nn=ii+1 3560 if (nn>3) nn=nn-3 3561 mm=ii+2 3562 if (mm>3) mm=mm-3 3563 do jj=1,ngfft(nn) 3564 do kk=1,ngfft(mm) 3565 ! The calcul of the second derivations on the grid 3566 call inspln(ii,jj,kk) 3567 end do 3568 end do 3569 nullify(ptd,ptc,ptp) 3570 end do 3571 nullify(ptd,ptc,ptp) 3572 3573 znucl_batom=znucl(typat(aim_dtset%batom)) 3574 3575 ABI_FREE(znucl) 3576 ABI_FREE(zionpsp) 3577 ABI_FREE(symrel) 3578 ABI_FREE(tnons) 3579 3580 !the pointers are obsolete - to remove later 3581 3582 end subroutine initaim
m_bader/inpar [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
inpar
FUNCTION
Parser for the aim utility (shorter than the one of ABINIT)
INPUTS
This routine uses data from the defs_aimprom module
OUTPUT
instr=string of character containing the input data lenstr=actual length of the character string WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
3604 subroutine inpar(instr,lenstr) 3605 3606 !Arguments ------------------------------------ 3607 !scalars 3608 integer,intent(out) :: lenstr 3609 character(len=*),intent(out) :: instr 3610 3611 !Local variables ------------------------------ 3612 character(len=1),parameter :: space=' ' 3613 character(len=26),parameter :: uplett='ABCDEFGHIJKLMNOPQRSTUVWXYZ', lolett='abcdefghijklmnopqrstuvwxyz' 3614 !scalars 3615 integer,parameter :: nline=100 3616 integer :: ii,inxh,inxl,ios,jj,kk,ll 3617 character(len=fnlen) :: line 3618 3619 ! ********************************************************************* 3620 3621 lenstr=0 3622 3623 do ii=1,26 3624 inxh=index(lolett,uplett(ii:ii)) 3625 if (inxh > 0) then 3626 write(std_out,*) 'ERROR The ', uplett(ii:ii) ,' is considered come lowcase !' 3627 ABI_ERROR("Aborting now") 3628 end if 3629 end do 3630 rewind(unt0) 3631 do ii=1,nline 3632 read(unt0,'(A)',iostat=ios) line(1:fnlen) 3633 if (ios/=0) exit 3634 inxh=index(line,'#') 3635 if (inxh == 1) then 3636 cycle 3637 elseif (inxh > 0) then 3638 inxl=inxh-1 3639 line(inxh:inxh)=space 3640 else 3641 inxl=len_trim(line) 3642 if (inxl==0) cycle 3643 end if 3644 inxh=index(line(1:inxl),char(9)) 3645 if (inxh/=0) line(inxh:inxh)=space 3646 do ll=1,inxl 3647 if (iachar(line(ll:ll)) < 32) line(ll:ll)=space 3648 end do 3649 inxh=index(line(1:inxl),'- ') 3650 if (inxh/=0) then 3651 write(std_out,*) 'ERROR sign minus with white space in input file' 3652 ABI_ERROR("Aborting now") 3653 end if 3654 line(1:inxl)=adjustl(line(1:inxl)) 3655 inxl=len_trim(line(1:inxl))+1 3656 jj=2;kk=0 3657 line(1:inxl)=adjustl(line(1:inxl)) 3658 kk=len_trim(line(1:inxl))+1 3659 do ll=1,inxl 3660 inxh=index(line(jj:kk),space) 3661 if ((inxh==0).or.((jj+inxh-1)==kk)) exit 3662 line(inxh+jj:kk)=adjustl(line(inxh+jj:kk)) 3663 kk=len_trim(line(1:inxl)) 3664 if (kk == inxl) then 3665 exit 3666 end if 3667 jj=jj+inxh 3668 end do 3669 inxl=len_trim(line(1:inxl))+1 3670 do ll=1,inxl-1 3671 inxh=index(lolett,line(ll:ll)) 3672 if (inxh/=0) line(ll:ll)=uplett(inxh:inxh) 3673 end do 3674 if ((lenstr+inxl) > strlen ) then 3675 write(std_out,*) 'ERROR Too large input !' 3676 ABI_ERROR("Aborting now") 3677 else 3678 instr(lenstr+1:lenstr+inxl)=line(1:inxl) 3679 lenstr=lenstr+inxl 3680 end if 3681 end do 3682 end subroutine inpar
m_bader/inspln [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
inspln
FUNCTION
This procedure gives the values of the spline coefficients (second derivatives) in the 1D grid with periodic boundary conditions at rsid - the values of the unknown functions specified in the vector valf of direction idir
INPUTS
idir= direction following which the derivatives are evaluated
snn, tnn=remaining bi-dimensional coordinates of the line along which
the derivative is to be computed
OUTPUT
(see side effects)
SIDE EFFECTS
This routine works on the data contained in the aimfields module WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
3711 subroutine inspln(idir,snn,tnn) 3712 3713 !Arguments ------------------------------------ 3714 !scalars 3715 integer,intent(in) :: idir,snn,tnn 3716 3717 !Local variables------------------------------- 3718 !scalars 3719 integer :: dim,ii 3720 real(dp) :: ss 3721 !arrays 3722 real(dp) :: rsid(ngfft(idir)),valf(ngfft(idir)) 3723 real(dp),pointer :: ptc(:),ptd(:),ptp(:) 3724 3725 ! ************************************************************************* 3726 3727 !POINTER INITIALIZATION 3728 3729 if (idir==1) then 3730 valf(:)=dvl(:,snn,tnn) 3731 elseif (idir==2) then 3732 valf(:)=dvl(tnn,:,snn) 3733 else 3734 valf(:)=dvl(snn,tnn,:) 3735 end if 3736 3737 nullify(ptd,ptc,ptp) 3738 if(idir==1) then 3739 ptd=>dig1;ptc=>cdig1;ptp=>llg1 3740 elseif (idir==2) then 3741 ptd=>dig2;ptc=>cdig2;ptp=>llg2 3742 else 3743 ptd=>dig3;ptc=>cdig3;ptp=>llg3 3744 end if 3745 3746 dim=ngfft(idir) 3747 3748 !FIRST CYCLE OF RECURRENCE 3749 3750 rsid(1)=valf(2)+valf(dim)-2.*valf(1) 3751 rsid(1)=rsid(1)/ptd(1) 3752 do ii=2,dim-1 3753 rsid(ii)=valf(ii+1)+valf(ii-1)-2.*valf(ii) 3754 rsid(ii)=(rsid(ii)-ptc(ii-1)*rsid(ii-1))/ptd(ii) 3755 end do 3756 ss=0._dp 3757 do ii=1,dim-1 3758 ss=ss+rsid(ii)*ptp(ii) 3759 end do 3760 rsid(dim)=valf(1)+valf(dim-1)-2.*valf(dim) 3761 rsid(dim)=(rsid(dim)-ss)/ptd(dim) 3762 3763 !SECOND CYCLE WITH TRANSPOSED MATRIX 3764 3765 rsid(dim)=rsid(dim)/ptd(dim) 3766 rsid(dim-1)=(rsid(dim-1)-ptc(dim-1)*rsid(dim))/ptd(dim-1) 3767 do ii=dim-2,1,-1 3768 rsid(ii)=(rsid(ii)-ptc(ii)*rsid(ii+1)-ptp(ii)*rsid(dim))/ptd(ii) 3769 end do 3770 3771 if (idir==1) then 3772 ddx(:,snn,tnn)=rsid(:) 3773 elseif (idir==2) then 3774 ddy(tnn,:,snn)=rsid(:) 3775 else 3776 ddz(snn,tnn,:)=rsid(:) 3777 end if 3778 3779 end subroutine inspln
m_bader/integrho [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
integrho
FUNCTION
This routine integrates the electron density inside the atomic surface already calculated - it reads the file *.surf The radial integration is always performed with splines and the two angular integrations with Gauss quadrature
INPUTS
aim_dtset = the structured entity containing all input variables znucl_batom=the nuclear charge of the Bader atom
OUTPUT
(see side effects)
SIDE EFFECTS
This routine works primarily on the data contained in the aimfields and aimprom modules WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
3807 subroutine integrho(aim_dtset,znucl_batom) 3808 3809 !Arguments ------------------------------------ 3810 !scalars 3811 type(aim_dataset_type),intent(in) :: aim_dtset 3812 3813 !Local variables ------------------------------ 3814 !scalars 3815 integer :: batom,chs,iat,ii,inx,inxf,ipos,jj,kk,ll,nn,nph,nth 3816 real(dp) :: chg,chgint,cintr,ct1,ct2,lder,nsphe,phimax,phimin,rder 3817 real(dp) :: rsmax,rsmin,ss,stp,themax,themin,uu 3818 real(dp) :: znucl_batom,zz 3819 logical :: gaus,weit 3820 !arrays 3821 real(dp) :: grho(3),hrho(3,3),shift(3),unvec(3),vv(3) 3822 real(dp),allocatable :: ncrho(:),nsp2(:),nsp3(:),nsp4(:),rdint(:,:),rr(:) 3823 real(dp),allocatable :: vdd(:),vrho(:),wgrs(:,:),work(:) 3824 3825 ! ********************************************************************* 3826 3827 gaus=.true. 3828 weit=.true. 3829 3830 write(std_out,*) 'npt = ',aim_dtset%npt 3831 3832 rewind(unts) 3833 read(unts,*) batom,shift ! Warning : batom is read, instead of coming from aim_dtset 3834 read(unts,*) nth,themin,themax ! Warning : these numbers are read, instead of coming from aim_dtset 3835 read(unts,*) nph,phimin,phimax ! Warning : these numbers are read, instead of coming from aim_dtset 3836 3837 write(std_out,*) 'NTH NPH ',nth,nph 3838 3839 ABI_MALLOC(wgrs,(nth,nph)) 3840 ABI_MALLOC(rdint,(nth,nph)) 3841 3842 do ii=1,nth 3843 do jj=1,nph 3844 if (weit) then 3845 read(unts,*) th(ii),ph(jj),rs(ii,jj),wgrs(ii,jj) 3846 else 3847 read(unts,*) th(ii),ph(jj),rs(ii,jj) 3848 end if 3849 end do 3850 end do 3851 read(unts,*) rsmin,rsmax 3852 3853 3854 if (gaus) then 3855 ct1=cos(themin) 3856 ct2=cos(themax) 3857 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 3858 call coeffs_gausslegint(phimin,phimax,ph,wph,nph) 3859 end if 3860 3861 do ii=1,nth 3862 do jj=1,nph 3863 if (.not.weit) then 3864 if (gaus) then 3865 wgrs(ii,jj)=wcth(ii)*wph(jj) 3866 else 3867 wgrs(ii,jj)=1._dp 3868 end if 3869 end if 3870 end do 3871 end do 3872 3873 3874 do ii=1,nth 3875 do jj=1,nph 3876 if (rs(ii,jj) < rsmin) rsmin=rs(ii,jj) 3877 end do 3878 end do 3879 3880 3881 !INTEGRATION OF THE CORE DENSITY 3882 3883 nn=typat(batom) 3884 kk=ndat(nn) 3885 3886 3887 !spherical integration of the core density in the sphere 3888 !of the minimal Bader radius 3889 3890 !COEF. FOR SPHERICAL INTEGRATION 3891 3892 ABI_MALLOC(nsp2,(kk)) 3893 ABI_MALLOC(nsp3,(kk)) 3894 ABI_MALLOC(nsp4,(kk)) 3895 ABI_MALLOC(ncrho,(kk)) 3896 3897 do ii=1,kk 3898 ncrho(ii)=crho(ii,nn)*4._dp*pi*rrad(ii,nn)*rrad(ii,nn) 3899 nsp3(ii)=4._dp*pi*(2._dp*crho(ii,nn)+2._dp*rrad(ii,nn)*sp2(ii,nn)+& 3900 & rrad(ii,nn)*rrad(ii,nn)*sp3(ii,nn)) 3901 end do 3902 3903 if (rsmin < rrad(ndat(nn),nn)) then ! search index 3904 inx=0 3905 if (rsmin < rrad(1,nn)) then 3906 ABI_ERROR('absurd') 3907 elseif (rsmin > rrad(ndat(nn),nn)) then 3908 inx=ndat(nn) 3909 else 3910 do while (rsmin >= rrad(inx+1,nn)) 3911 inx=inx+1 3912 end do 3913 end if 3914 else 3915 inx=ndat(nn) 3916 end if 3917 3918 cintr=4._dp/3._dp*pi*rrad(1,nn)**3*crho(1,nn) 3919 3920 !spline integration 3921 3922 do ii=1,inx-1 3923 uu=rrad(ii+1,nn)-rrad(ii,nn) 3924 cintr=cintr+(ncrho(ii)+ncrho(ii+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(ii)+nsp3(ii+1)) 3925 end do 3926 if (inx/=ndat(nn)) then 3927 uu=rsmin-rrad(inx,nn) 3928 zz=rrad(inx+1,nn)-rsmin 3929 ss=rrad(inx+1,nn)-rrad(inx,nn) 3930 cintr=cintr+ncrho(inx)/2._dp*(ss-zz*zz/ss)+ncrho(inx+1)/2._dp*uu*uu/ss+& 3931 nsp3(inx)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 3932 nsp3(inx+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 3933 end if 3934 3935 3936 !INTEGRATION OF THE REST OF THE CORE DENSITY 3937 !(for gauss quadrature) 3938 !For the Gauss quadrature it is added 3939 !to the radial integrated valence density 3940 3941 rdint(:,:)=0._dp 3942 nsphe=0._dp 3943 do ii=1,nth 3944 do jj=1,nph 3945 if (inx==ndat(nn)) cycle 3946 inxf=inx 3947 if (rs(ii,jj) < rsmin) then 3948 write(std_out,*) rs(ii,jj),rsmin 3949 ABI_ERROR('in surface') 3950 elseif (rs(ii,jj) > rrad(ndat(nn),nn)) then 3951 inxf=ndat(nn) 3952 else 3953 do while (rs(ii,jj) >= rrad(inxf+1,nn)) 3954 inxf=inxf+1 3955 end do 3956 end if 3957 3958 if (inxf==inx) then 3959 uu=rrad(inx+1,nn)-rs(ii,jj) 3960 zz=rrad(inx+1,nn)-rsmin 3961 ss=rrad(inx+1,nn)-rrad(inx,nn) 3962 3963 rdint(ii,jj)=(ncrho(inx)/2._dp/ss-nsp3(inx)/1.2d1*ss)*(zz*zz-uu*uu)+& 3964 nsp3(inx)/2.4d1/ss*(zz**4-uu**4) 3965 uu=rs(ii,jj)-rrad(inx,nn) 3966 zz=rsmin-rrad(inx,nn) 3967 rdint(ii,jj)=rdint(ii,jj)+(uu*uu-zz*zz)*(ncrho(inx+1)/2._dp/ss-nsp3(inx+1)/1.2d1*ss)+& 3968 nsp3(inx+1)/2.4d1/ss*(uu**4-zz**4) 3969 else 3970 uu=rrad(inx+1,nn)-rsmin 3971 zz=rsmin-rrad(inx,nn) 3972 3973 rdint(ii,jj)=ncrho(inx)/2._dp/ss*uu*uu+ncrho(inx+1)/2._dp*(ss-zz*zz/ss)+& 3974 nsp3(inx)/1.2d1*(uu**4/2._dp/ss-uu*uu*ss)+nsp3(inx+1)/1.2d1*(zz*zz*ss-ss**3/2._dp-zz**4/2._dp/ss) 3975 if (inxf > inx+1) then 3976 do kk=inx+1,inxf-1 3977 uu=rrad(kk+1,nn)-rrad(kk,nn) 3978 rdint(ii,jj)=rdint(ii,jj)+(ncrho(kk)+ncrho(kk+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(kk)+nsp3(kk+1)) 3979 end do 3980 end if 3981 3982 if (inxf/=ndat(nn)) then 3983 uu=rs(ii,jj)-rrad(inxf,nn) 3984 zz=rrad(inxf+1,nn)-rs(ii,jj) 3985 ss=rrad(inxf+1,nn)-rrad(inxf,nn) 3986 rdint(ii,jj)=rdint(ii,jj)+ncrho(inxf)/2._dp*(ss-zz*zz/ss)+ncrho(inxf+1)/2._dp*uu*uu/ss+& 3987 nsp3(inxf)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 3988 nsp3(inxf+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 3989 end if 3990 end if 3991 rdint(ii,jj)=rdint(ii,jj)/4._dp/pi 3992 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 3993 end do 3994 end do 3995 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 3996 3997 write(untout,*) 3998 write(untout,*) "CHARGE INTEGRATION" 3999 write(untout,*) "==================" 4000 write(untout,'(" Core density contribution: ",/,/," ",F16.8)') cintr+nsphe 4001 4002 write(std_out,*) ':INTECOR ', cintr+nsphe 4003 4004 ABI_FREE(ncrho) 4005 ABI_FREE(nsp2) 4006 ABI_FREE(nsp3) 4007 ABI_FREE(nsp4) 4008 4009 !INTEGRATION OF THE VALENCE DENSITY 4010 4011 ABI_MALLOC(rr,(aim_dtset%npt+1)) 4012 ABI_MALLOC(vrho,(aim_dtset%npt+1)) 4013 ABI_MALLOC(vdd,(aim_dtset%npt+1)) 4014 4015 !in the case of the only irho appelation 4016 4017 nn=0 4018 do ii=-3,3 4019 do jj=-3,3 4020 do kk=-3,3 4021 nn=nn+1 4022 atp(1,nn)=ii*1._dp 4023 atp(2,nn)=jj*1._dp 4024 atp(3,nn)=kk*1._dp 4025 call bschg1(atp(:,nn),1) 4026 if ((ii==0).and.(jj==0).and.(kk==0)) ipos=nn 4027 end do 4028 end do 4029 end do 4030 nnpos=nn 4031 iat=batom 4032 4033 !XG020629 There is a problem with this routine 4034 !(or vgh_rho), when one uses the PGI compiler : 4035 !The following line is needed, otherwise, iat and ipos 4036 !are set to 0 inside vgh_now. Why ???? 4037 write(std_out,*)' integrho : iat,ipos=',iat,ipos 4038 ! 4039 4040 nsphe=0._dp 4041 ABI_MALLOC(work,(aim_dtset%npt+1)) 4042 do ii=1,nth 4043 do jj=1,nph 4044 4045 stp=rs(ii,jj)/aim_dtset%npt 4046 unvec(1)=sin(th(ii))*cos(ph(jj)) 4047 unvec(2)=sin(th(ii))*sin(ph(jj)) 4048 unvec(3)=cos(th(ii)) 4049 do kk=0,aim_dtset%npt 4050 rr(kk+1)=kk*stp 4051 vv(:)=xatm(:,batom)+kk*stp*unvec(:) 4052 chs=-2 4053 call vgh_rho(vv,chg,grho,hrho,uu,iat,ipos,chs) 4054 vrho(kk+1)=chg*rr(kk+1)*rr(kk+1) 4055 if (kk==aim_dtset%npt) then 4056 rder=0._dp 4057 do ll=1,3 4058 rder=rder+grho(ll)*unvec(ll) 4059 end do 4060 rder=rder*rr(kk+1)*rr(kk+1)+2._dp*rr(kk+1)*chg 4061 end if 4062 end do 4063 lder=0._dp 4064 kk=aim_dtset%npt+1 4065 call spline(rr,vrho,kk,lder,rder,vdd) 4066 4067 ! INTEGRATION 4068 4069 do kk=1,aim_dtset%npt 4070 rdint(ii,jj)=rdint(ii,jj)+stp/2._dp*(vrho(kk)+vrho(kk+1))& 4071 & -stp*stp*stp/24._dp*(vdd(kk)+vdd(kk+1)) 4072 end do 4073 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 4074 end do 4075 end do 4076 ABI_FREE(work) 4077 4078 if (gaus.or.weit) then 4079 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 4080 else 4081 nsphe=nsphe/(nth*nph)*2.0*two_pi 4082 end if 4083 chgint=cintr+nsphe 4084 4085 write(untout,'(/," Different density contributions: Core (only spherical part) and the rest ",/,/," ",2F16.8)') & 4086 & cintr, nsphe 4087 write(untout,'(/,a,i4,a,f14.8)') ' For atom number ',batom,', the number of electrons in the Bader volume is ',chgint 4088 write(untout,'(a,f15.7,a,f17.8)') ' The nuclear charge is',znucl_batom,', so that the Bader charge is ',znucl_batom-chgint 4089 write(untout,*) 4090 write(std_out,*) ':INTEPAR ', cintr, nsphe 4091 write(std_out,*) ':RHOTOT ',batom,chgint 4092 4093 end subroutine integrho
m_bader/integvol [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
integvol
FUNCTION
This routine integrates the volume of the Bader atom
INPUTS
(see side effects)
OUTPUT
(see side effects)
SIDE EFFECTS
This routine works on the data contained in the aimfields and aimprom modules WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
4117 subroutine integvol() 4118 4119 !Arguments ------------------------------------ 4120 4121 !Local variables ------------------------------ 4122 !scalars 4123 integer :: batom,ii,jj,nph,nth 4124 real(dp) :: chgint,ct1,ct2,nsphe,phimax,phimin 4125 real(dp) :: rsmax,rsmin,themax,themin 4126 logical :: gaus,weit 4127 !arrays 4128 real(dp) :: shift(3) 4129 real(dp),allocatable :: rdint(:,:) 4130 real(dp),allocatable :: wgrs(:,:) 4131 4132 ! ********************************************************************* 4133 4134 tpi=two_pi 4135 gaus=.true. 4136 weit=.true. 4137 4138 4139 rewind(unts) 4140 read(unts,*) batom,shift 4141 read(unts,*) nth,themin,themax 4142 read(unts,*) nph,phimin,phimax 4143 4144 write(std_out,*) 'NTH NPH ',nth,nph 4145 4146 ABI_MALLOC(wgrs,(nth,nph)) 4147 ABI_MALLOC(rdint,(nth,nph)) 4148 4149 do ii=1,nth 4150 do jj=1,nph 4151 if (weit) then 4152 read(unts,*) th(ii),ph(jj),rs(ii,jj),wgrs(ii,jj) 4153 else 4154 read(unts,*) th(ii),ph(jj),rs(ii,jj) 4155 end if 4156 end do 4157 end do 4158 read(unts,*) rsmin,rsmax 4159 4160 4161 if (gaus) then 4162 ct1=cos(themin) 4163 ct2=cos(themax) 4164 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 4165 call coeffs_gausslegint(phimin,phimax,ph,wph,nph) 4166 end if 4167 4168 do ii=1,nth 4169 do jj=1,nph 4170 if (.not.weit) then 4171 if (gaus) then 4172 wgrs(ii,jj)=wcth(ii)*wph(jj) 4173 else 4174 wgrs(ii,jj)=1._dp 4175 end if 4176 end if 4177 end do 4178 end do 4179 4180 nsphe=0._dp 4181 do ii=1,nth 4182 do jj=1,nph 4183 nsphe=nsphe+rs(ii,jj)**3/3._dp*wgrs(ii,jj) 4184 end do 4185 end do 4186 if (gaus.or.weit) then 4187 nsphe=nsphe*(pi/(themin-themax))*(tpi/(phimax-phimin)) 4188 else 4189 nsphe=nsphe/(nth*nph)*2.0*tpi 4190 end if 4191 chgint=nsphe 4192 4193 write(std_out,*) ':VOLTOT ',batom,chgint 4194 write(untout,'("Volume of the Bader atom: ", I6, F16.8)') batom,chgint 4195 4196 end subroutine integvol
m_bader/mprod [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
mprod
FUNCTION
Matrix multiplication cc=aa*bb
SOURCE
5779 subroutine mprod(aa,bb,cc) 5780 5781 !Arguments ------------------------------------ 5782 !arrays 5783 real(dp),intent(in) :: aa(3,3),bb(3,3) 5784 real(dp),intent(out) :: cc(3,3) 5785 5786 !Local variables------------------------------- 5787 !scalars 5788 integer :: ii,jj,kk 5789 5790 ! ************************************************************************* 5791 5792 do ii=1,3 5793 do jj=1,3 5794 cc(ii,jj)=0._dp 5795 do kk=1,3 5796 cc(ii,jj)=cc(ii,jj)+aa(ii,kk)*bb(kk,jj) 5797 end do 5798 end do 5799 end do 5800 5801 end subroutine mprod
m_bader/onestep [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
onestep
FUNCTION
Advance one step following the gradient from vv(3). It returns a new point in vv(3) and the value and gradient of the electron density at this point in chg and grho(3)
INPUTS
npmax= maximum number of divisions hh= determines the initial value of the step (to be multiplied by grho)
OUTPUT
chg= value of electron density deltar= the length of the step thaty was needed grho(3)= gradient of electron density np= returns the number of divisions that were needed
SIDE EFFECTS
vv(3)=starting and updated point WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
4226 subroutine onestep(vv,chg,grho,hh,np,npmax,deltar) 4227 4228 !Arguments ------------------------------------ 4229 !scalars 4230 integer,intent(in) :: npmax 4231 integer,intent(out) :: np 4232 real(dp),intent(in) :: hh 4233 real(dp),intent(out) :: chg,deltar 4234 !arrays 4235 real(dp),intent(inout) :: vv(3) 4236 real(dp),intent(out) :: grho(3) 4237 4238 !Local variables ------------------------------ 4239 !scalars 4240 integer :: iat,ii,ipos,jj 4241 real(dp) :: dt,rr 4242 !arrays 4243 real(dp) :: hrho(3,3),pom(3),vinter(3,200),vk(3),vkold(3) 4244 4245 !************************************************************************ 4246 dt=hh 4247 np=1 4248 deltar=1._dp 4249 vk(1:3)=vv(1:3) 4250 4251 4252 do while((np<3).or.((np<=npmax).and.(deltar>aim_deltarmin))) 4253 np=np*2 4254 dt=dt*0.5_dp 4255 vkold(1:3)=vk(1:3) 4256 call vgh_rho(vk,chg,grho,hrho,rr,iat,ipos,0) 4257 vinter(1:3,1)=vv(1:3)+dt*grho(1:3) 4258 do jj=2,np 4259 call vgh_rho(vinter(1,jj-1),chg,grho,hrho,rr,iat,ipos,0) 4260 if(jj.eq.2) then 4261 vinter(1:3,2)=vv(1:3)+2.0*dt*grho(1:3) 4262 else 4263 vinter(1:3,jj)=vinter(1:3,jj-2)+2.0*dt*grho(1:3) 4264 end if 4265 end do 4266 4267 call vgh_rho(vinter(1,np),chg,grho,hrho,rr,iat,ipos,0) 4268 vinter(1:3,np+1)=vinter(1:3,np-1)+dt*grho(1:3) 4269 4270 deltar=0._dp 4271 do ii=1,3 4272 vk(ii)=(vinter(ii,np)+vinter(ii,np+1))*0.5_dp 4273 deltar=deltar+(vkold(ii)-vk(ii))*(vkold(ii)-vk(ii)) 4274 end do 4275 end do 4276 4277 pom(:)=vk(:)-vv(:) 4278 deltar=vnorm(pom,0) 4279 vv(1:3)=vk(1:3) 4280 4281 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 4282 if(deb) write(std_out,*) ':VKf ',np,vk 4283 4284 end subroutine onestep
m_bader/ordr [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
ordr
FUNCTION
INPUTS
(to be filled)
OUTPUT
(to be filled)
SOURCE
2426 subroutine ordr(aa,dd,nn,cff) 2427 2428 !Arguments ---------------------------- 2429 !scalars 2430 integer,intent(in) :: cff,nn 2431 !arrays 2432 real(dp),intent(inout) :: aa(nn),dd(nn,nn) 2433 2434 !Local variables ---------------------- 2435 !scalars 2436 integer :: ii,jj,kk 2437 real(dp) :: uu 2438 2439 ! ********************************************************************* 2440 2441 do ii=1,nn-1 2442 kk=ii 2443 uu=aa(ii) 2444 do jj=ii+1,nn 2445 if (cff==1) then 2446 if (aa(jj) >= uu+tol12) then 2447 kk=jj 2448 uu=aa(jj) 2449 end if 2450 else 2451 if (aa(jj) <= uu-tol12) then 2452 kk=jj 2453 uu=aa(jj) 2454 end if 2455 end if 2456 end do 2457 if (kk /= ii) then 2458 aa(kk)=aa(ii) 2459 aa(ii)=uu 2460 do jj=1,nn 2461 uu=dd(jj,ii) 2462 dd(jj,ii)=dd(jj,kk) 2463 dd(jj,kk)=uu 2464 end do 2465 end if 2466 end do 2467 end subroutine ordr
m_bader/plint [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
plint
FUNCTION
This simple routine gives the profile of the density integrated in xy plane belong the z-axes (it works only for orthogonal coordinates at present - it is better to use cut3d) integration in plane - with equilateral triangles (not really finished and not tested!)
INPUTS
(this routine works on the data in the aimprom module)
OUTPUT
(this routine works on the data in the aimprom module) WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
4309 subroutine plint() 4310 4311 !Arguments ------------------------------------ 4312 4313 !Local variables ------------------------------ 4314 !scalars 4315 integer,parameter :: nd=150,ng=300 4316 integer :: cod,iat,ii,ipos,jj,kk,nn 4317 real(dp) :: dd,ee,ff,gg,hh,igr,rho,ss 4318 logical :: prep 4319 !arrays 4320 real(dp) :: grho(3),hrho(3,3),vv(3),xl(nd+1),xs(nd) 4321 real(dp),allocatable :: uu(:) 4322 4323 ! ********************************************************************* 4324 4325 ff=rprimd(1,1)/nd 4326 ss=2._dp/sqrt(3._dp)*rprimd(2,2)/rprimd(1,1)*nd 4327 nn=int(ss) 4328 gg=sqrt(3._dp)/2.*ff 4329 hh=rprimd(2,2)-nn/nd*sqrt(3._dp)/2.*rprimd(1,1) 4330 ee=hh/sqrt(3._dp) 4331 hh=hh/2. 4332 ss=sqrt(3._dp)*ff*ff/4. 4333 dd=ee*ff/2. 4334 4335 do ii=1,nd 4336 xl(ii)=ii*ff 4337 xs(ii)=ff/2.+ii*ff 4338 end do 4339 xl(nd+1)=rprimd(1,1) 4340 4341 ABI_MALLOC(uu,(nn+3)) 4342 4343 uu(1)=0._dp 4344 uu(nn+3)=rprimd(2,2) 4345 do ii=2,nn+2 4346 uu(ii)=hh+(ii-1)*gg 4347 end do 4348 igr=0._dp 4349 prep=.true. 4350 do kk=1,ng 4351 igr=0._dp 4352 vv(3)=(kk-1)*rprimd(3,3)/ng 4353 do ii=1,nn+3 4354 vv(2)=uu(ii) 4355 do jj=1,nd 4356 if (prep) then 4357 vv(1)=xl(jj) 4358 prep=.false. 4359 else 4360 vv(1)=xs(jj) 4361 prep=.true. 4362 end if 4363 call vgh_rho(vv,rho,grho,hrho,dd,iat,ipos,cod) 4364 if ((ii==1).or.(ii==nn+3)) then 4365 igr=igr+dd*rho 4366 elseif ((ii==2).or.(ii==nn+2)) then 4367 igr=igr+(dd+ss)*rho 4368 else 4369 igr=igr+ss*2*rho 4370 end if 4371 end do 4372 end do 4373 write(untp,'(2E16.8)') vv(3), igr 4374 end do 4375 ABI_FREE(uu) 4376 4377 end subroutine plint
m_bader/rsurf [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
rsurf
FUNCTION
Basic routine for determination of the radius of Bader surface for spherical rayon theta,phi the bassin is tested by following the gradient line If srch==true (in general for calls from surf) the routine aim_follow is called to stop when it arrives under already known part of surface Simple bissection method is used to obtain the radius WARNING This file does not follow the ABINIT coding rules (yet)
INPUTS
aim_dtset= the structured entity containing all input variables rr0= starting radius theta,phi = the spherical direction iatinit= the atom index srch= see above npmax= maximum number of divisions in one step for follow
OUTPUT
rr= radius grho(3)= gradient on the surface
SOURCE
4409 subroutine rsurf(aim_dtset,rr,grho,theta,phi,rr0,iatinit,npmax,srch) 4410 4411 !Arguments ------------------------------------ 4412 !scalars 4413 integer,intent(in) :: iatinit,npmax 4414 real(dp),intent(in) :: phi,rr0,theta 4415 real(dp),intent(out) :: rr 4416 logical,intent(in) :: srch 4417 !arrays 4418 real(dp),intent(out) :: grho(3) 4419 !no_abirules 4420 type(aim_dataset_type),intent(in) :: aim_dtset 4421 4422 !Local variables ------------------------------ 4423 !scalars 4424 integer :: iat,ii,ipos,iposinit,jj,nstep 4425 real(dp),parameter :: mfkt=1.d1 4426 real(dp) :: aa,dmax,dr,drr,rho,rr1,rr2,t1,t2,wall 4427 logical :: cross,deb_tmp,in,in1,in2,low,srch_tmp 4428 !arrays 4429 real(dp) :: hrho(3,3),unvec(3),vv(3) 4430 4431 ! ********************************************************************* 4432 4433 srch_tmp=srch 4434 deb_tmp=deb 4435 4436 !unity vecteur in the direction (theta,phi) 4437 4438 unvec(1)=sin(theta)*cos(phi) 4439 unvec(2)=sin(theta)*sin(phi) 4440 unvec(3)=cos(theta) 4441 4442 4443 rr=rr0 4444 rr1=rr 4445 rr2=rr 4446 drr=1._dp 4447 if (abs(rr0-r0)<1.0d-12) then 4448 dr=aim_dtset%dr0*mfkt 4449 else 4450 dr=aim_dtset%dr0 4451 end if 4452 4453 vv(1)=xatm(1,aim_dtset%batom) 4454 vv(2)=xatm(2,aim_dtset%batom) 4455 vv(3)=xatm(3,aim_dtset%batom) 4456 4457 4458 iposinit=batcell 4459 write(std_out,'("ATOM iat=",i4," ipos=",i4)') aim_dtset%batom,batcell 4460 jj=0 4461 4462 cross=.false. 4463 4464 in=.true. 4465 low=.false. 4466 4467 dmax=h0 4468 4469 in1=.true. 4470 in2=in1 4471 4472 do while((drr>aim_drmin).or.(jj<2)) 4473 call timein(t1,wall) 4474 jj=jj+1 4475 do ii=1,3 4476 vv(ii)=xatm(ii,aim_dtset%batom)+rr*unvec(ii) 4477 end do 4478 4479 ! VACUUM CONDITION 4480 4481 call vgh_rho(vv,rho,grho,hrho,aa,iat,ipos,0) 4482 if (rho < aim_rhomin) exit 4483 4484 ldeb=.false. 4485 4486 call aim_follow(aim_dtset,vv,npmax,srch_tmp,iatinit,iposinit,iat,ipos,nstep) 4487 4488 call timein(t2,wall) 4489 t2=t2-t1 4490 4491 write(std_out,'(a,i4,a,f12.8,a,i4,a,i4,a,f10.5,a,i4)') & 4492 & ' :STEP ',jj,' r=',rr,' iat=',iat,' ipos=',ipos,' time(sec)=',t2,' nstep=',nstep 4493 4494 if ((iat.eq.iatinit).and.(ipos.eq.iposinit)) then 4495 in=.true. 4496 else 4497 in=.false. 4498 end if 4499 4500 ! 4501 ! NEW RADIUS 4502 ! 4503 4504 if ((jj.eq.1).or.((in1.eqv.in).and.(.not.cross))) then 4505 if (in) then 4506 rr2=rr1 4507 rr1=rr 4508 rr=rr+dr 4509 else 4510 rr2=rr1 4511 rr1=rr 4512 rr=rr-dr 4513 end if 4514 if ((jj>2).and.(dr<(0.6))) then 4515 ! modification of the step 4516 dr=dr*aim_fac 4517 if (deb_tmp) write(std_out,*) ':DR ',dr 4518 end if 4519 else 4520 if (.not.cross) then 4521 cross=.true. 4522 rr2=rr1 4523 else 4524 if (in2) then 4525 if (in) then 4526 rr2=rr1 4527 else 4528 in1=in2 4529 end if 4530 else 4531 if (in) then 4532 in1=in2 4533 else 4534 rr2=rr1 4535 end if 4536 end if 4537 end if 4538 rr1=rr 4539 rr=(rr2+rr1)/2.0 4540 end if 4541 4542 in2=in1 4543 in1=in 4544 drr=abs(rr2-rr1)/rr 4545 if (deb_tmp) write(std_out,*) ':DRR ',jj,rr2,rr1,drr 4546 end do 4547 4548 end subroutine rsurf
m_bader/surf [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
surf
FUNCTION
Determination of the Bader surface. Use rsurf to determine radius for one direction simple bisection method is used the bassin is tested following the gradient (follow) = = the most time consuming follow stops if the gradient line is near the atom or if it is under already known part of surface - this is why the surface is not computed row by row.
INPUTS
aim_dtset= the structured entity containing all input variables
OUTPUT
(see side effects)
SIDE EFFECTS
This routine works primarily on the data contained in the defs_aimprom module WARNING This file does not follow the ABINIT coding rules (yet)
SOURCE
4579 subroutine surf(aim_dtset) 4580 4581 !Arguments ------------------------------------ 4582 !scalars 4583 type(aim_dataset_type) :: aim_dtset 4584 4585 !Local variables ------------------------------ 4586 !scalars 4587 integer :: ierr,ii,ijj,ijj_exist,incr,init,iph,iph2,ith,ith2,jj,jj_exist,kk,level,me,mm,nn,nph,npmax,nproc,nth,comm 4588 real(dp) :: ct1,ct2,phi,rr,rsmax,rsmin,rthe,rthe0,t1,t2,theta,tt0,vcth,vph,vth 4589 real(dp) :: wall,xy,xyz 4590 logical :: srch,stemp 4591 !arrays 4592 real(dp) :: grho(3),vr(3),vv(3) 4593 real(dp),allocatable :: rs_computed(:,:) 4594 4595 !************************************************************************ 4596 4597 comm = xmpi_world 4598 me=xmpi_comm_rank(comm) 4599 nproc=xmpi_comm_size(comm) 4600 4601 ttsrf=zero 4602 4603 rewind(unts) 4604 4605 nth=aim_dtset%nth 4606 nph=aim_dtset%nph 4607 4608 !Coefficients for spherical Gauss quadrature 4609 4610 ct1=cos(aim_dtset%themin) 4611 ct2=cos(aim_dtset%themax) 4612 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 4613 call coeffs_gausslegint(aim_dtset%phimin,aim_dtset%phimax,ph,wph,nph) 4614 4615 !DEBUG 4616 !write(std_out,*)' surf : wcth=',wcth(1:nth) 4617 !write(std_out,*)' surf : wph=',wph(1:nth) 4618 !ENDDEBUG 4619 4620 do ijj=1,nth 4621 th(ijj)=acos(cth(ijj)) 4622 if (aim_dtset%isurf/=-1) then 4623 do jj=1,nph 4624 rs(ijj,jj)=zero 4625 end do 4626 end if 4627 end do 4628 4629 npmax=aim_npmaxin 4630 rsmax=0.0 4631 rsmin=100.0 4632 rthe0=r0 4633 srch=.false. 4634 4635 do ijj=1,3 4636 vv(ijj)=xatm(ijj,aim_dtset%batom) 4637 end do 4638 4639 4640 write(std_out,*) 4641 write(std_out,*) "BADER SURFACE DETERMINATION" 4642 write(std_out,*) "===========================" 4643 write(std_out,*) 4644 4645 write(untout,*) 4646 write(untout,*) "BADER SURFACE DETERMINATION" 4647 write(untout,*) "===========================" 4648 write(untout,*) 4649 4650 write(std_out,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 4651 write(std_out,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4652 write(std_out,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4653 4654 write(untout,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 4655 write(untout,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4656 write(untout,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4657 4658 write(unts,'(i3,3F15.10)') aim_dtset%batom,vv 4659 write(unts,'(i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4660 write(unts,'(i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4661 4662 !write(std_out,*) 'npmax in surf= ',npmax 4663 4664 ith=0 4665 iph=0 4666 tt0=0._dp 4667 call timein(tt0,wall) 4668 4669 write(untout,*) 4670 write(untout,*) "DEVELOPMENT OF THE RADII DETERMINATIONS" 4671 write(untout,*) "========================================" 4672 write(untout,*) 4673 write(untout,*) "Determination near the CPs:" 4674 4675 !Determination of the CP neighbouring radii 4676 4677 if (aim_dtset%isurf/=-1) then 4678 4679 ! Precomputation of the value of the radii (for parallelisation) 4680 ! To make the output independent of the number of processors, but still 4681 ! cut down the CPU time, use a multigrid technique 4682 srch=.true. 4683 ABI_MALLOC(rs_computed,(nth,nph)) 4684 rs(:,:)=zero 4685 rs_computed(:,:)=zero 4686 kk=0 ; init=0 4687 do level=3,0,-1 4688 incr=2**level 4689 if(incr<nth .and. incr<nph)then 4690 rs_computed(:,:)=rs(1:nth,1:nph) 4691 rs(1:nth,1:nph)=zero 4692 do ijj=1,nth,incr 4693 do jj=1,nph,incr 4694 if(rs_computed(ijj,jj)<1.0d-12) then 4695 kk=kk+1 4696 if(mod(kk,nproc)==me)then 4697 ! Find an approximate starting radius, from the already computed ones 4698 if(init==0)then 4699 rthe=r0 4700 else 4701 ijj_exist=ijj ; if(mod(ijj-1,2*incr)>=incr)ijj_exist=ijj-incr 4702 jj_exist=jj ; if(mod(jj-1,2*incr)>=incr)jj_exist=jj-incr 4703 rthe=rs_computed(ijj_exist,jj_exist) 4704 if(rthe<1.0d-12)then 4705 write(std_out,*)' surf : there is a bug ! rthe=',rthe 4706 ABI_ERROR("Aborting now") 4707 end if 4708 end if 4709 call timein(t1,wall) ; t2=zero 4710 call rsurf(aim_dtset,rr,grho,th(ijj),ph(jj),rthe,aim_dtset%batom,npmax,srch) 4711 rs(ijj,jj)=rr 4712 if (deb) then 4713 call timein(t2,wall) ; t2=t2-t1 4714 write(std_out,*) ':CALCULATED NP',ijj,jj,th(ijj),ph(jj),rthe,npmax,rs(ijj,jj),t2 4715 end if 4716 end if 4717 end if 4718 end do ! jj 4719 end do ! ijj 4720 call xmpi_sum(rs,comm,ierr) 4721 ! Combine the set of already computed radii and the set of the newly computed, to obtain all computed. 4722 rs(1:nth,1:nph)=rs(1:nth,1:nph)+rs_computed(:,:) 4723 init=1 4724 end if 4725 end do 4726 ABI_FREE(rs_computed) 4727 4728 srch=.true. 4729 4730 do ijj=1,nbcp 4731 ! if ((icpc(ijj) == -1)) then 4732 rthe0=vnorm(pc(:,ijj),0) 4733 do jj=1,3 4734 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 4735 end do 4736 xy=vr(1)*vr(1)+vr(2)*vr(2) 4737 xyz=xy+vr(3)*vr(3) 4738 xyz=sqrt(xyz) 4739 4740 if (xy < aim_xymin) then 4741 vcth=1._dp 4742 if (vr(3) < 0._dp) vcth=-vcth 4743 vph=0._dp 4744 else 4745 vcth=vr(3)/xyz 4746 vph=atan2(vr(2),vr(1)) 4747 end if 4748 4749 vth=acos(vcth) 4750 write(untout,'(/," BCP: (index,theta,phi)",I4,2E16.8)') ijj,vth,vph 4751 4752 if (vth < th(1)) then 4753 ith=0 4754 else 4755 if (vth > th(nth)) then 4756 ith=nth 4757 else 4758 do ii=2,nth 4759 if (vth < th(ii)) then 4760 ith=ii-1 4761 exit 4762 end if 4763 end do 4764 end if 4765 end if 4766 4767 if (vph < ph(1)) then 4768 iph=0 4769 else 4770 if (vph > ph(nph)) then 4771 iph=nph 4772 else 4773 do ii=2,nph 4774 if (vph < ph(ii)) then 4775 iph=ii-1 4776 exit 4777 end if 4778 end do 4779 end if 4780 end if 4781 4782 write(untout,*) "ATOMIC RADII (ith,iph,theta,phi,radius)" 4783 do jj=-1,2 4784 do kk=-1,2 4785 ith2=ith+jj 4786 iph2=iph+kk 4787 stemp=(iph2 > 0).and.(iph2 < nph+1) 4788 stemp=(stemp.and.((ith2 > 0).and.(ith2 < nth+1))) 4789 if (stemp) then 4790 theta=th(ith2) 4791 phi=ph(iph2) 4792 if (abs(rs(ith2,iph2))<1.0d-12) then 4793 rthe=rthe0 4794 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 4795 call timein(t1,wall) 4796 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 4797 call timein(t2,wall) 4798 t2=t2-t1 4799 rs(ith2,iph2)=rr 4800 end if 4801 rr=rs(ith2,iph2) 4802 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 4803 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 4804 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 4805 rthe0=rr 4806 end if 4807 4808 end do ! kk 4809 end do ! jj 4810 4811 ! end if 4812 4813 end do ! ijj (loop on BCP) 4814 4815 ! DEBUG 4816 ! write(std_out,*)' surf : near BCP ' 4817 ! do ijj=1,nth 4818 ! do jj=1,nph 4819 ! write(std_out,*)ijj,jj,rs(ijj,jj) 4820 ! end do 4821 ! end do 4822 ! ENDDEBUG 4823 4824 4825 srch=.true. 4826 do ijj=nbcp+1,nbcp+nrcp ! Loop on RCP 4827 ! if ((icpc(ijj) == 1)) then 4828 rthe0=max(rminl(aim_dtset%batom),r0) 4829 do jj=1,3 4830 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 4831 end do 4832 xy=vr(1)*vr(1)+vr(2)*vr(2) 4833 xyz=xy+vr(3)*vr(3) 4834 xyz=sqrt(xyz) 4835 4836 if (xy < aim_xymin) then 4837 vcth=1._dp 4838 if (vr(3) < 0._dp) vcth=-vcth 4839 vph=0._dp 4840 else 4841 vcth=vr(3)/xyz 4842 vph=atan2(vr(2),vr(1)) 4843 end if 4844 vth=acos(vcth) 4845 write(untout,'(/,";RCP: (index,theta,phi)",I4,2E16.8)') ijj-nbcp,vth,vph 4846 4847 if (vth < th(1)) then 4848 ith=0 4849 else 4850 if (vth > th(nth)) then 4851 ith=nth 4852 else 4853 do ii=2,nth 4854 if (vth < th(ii)) then 4855 ith=ii-1 4856 exit 4857 end if 4858 end do 4859 end if 4860 end if 4861 4862 if (vph < ph(1)) then 4863 iph=0 4864 else 4865 if (vph > ph(nph)) then 4866 iph=nph 4867 else 4868 do ii=2,nph 4869 if (vph < ph(ii)) then 4870 iph=ii-1 4871 exit 4872 end if 4873 end do 4874 end if 4875 end if 4876 4877 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 4878 do jj=-1,2 4879 do kk=-1,2 4880 ith2=ith+jj 4881 iph2=iph+kk 4882 stemp=(iph2 > 0).and.(iph2 < nph+1) 4883 stemp=stemp.and.(ith2 > 0).and.(ith2 < nth+1) 4884 4885 if (stemp) then 4886 theta=th(ith2) 4887 phi=ph(iph2) 4888 if ((abs(rs(ith2,iph2))<1.0d-12)) then 4889 rthe=rthe0 4890 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 4891 call timein(t1,wall) 4892 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 4893 call timein(t2,wall) 4894 t2=t2-t1 4895 rs(ith2,iph2)=rr 4896 end if 4897 rr=rs(ith2,iph2) 4898 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 4899 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 4900 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 4901 rthe0=rr 4902 end if 4903 4904 end do ! kk 4905 end do ! jj 4906 ! end if 4907 4908 end do ! ijj (Loop on RCP) 4909 4910 ! DEBUG 4911 ! write(std_out,*)' surf : near RCP ' 4912 ! do ijj=1,nth 4913 ! do jj=1,nph 4914 ! write(std_out,*)ijj,jj,rs(ijj,jj) 4915 ! end do 4916 ! end do 4917 ! ENDDEBUG 4918 4919 ! Boundary angles 4920 rthe0=r0 4921 srch=.true. 4922 write(untout,*) 4923 write(untout,*) "The boundary angles:" 4924 write(untout,*) "====================" 4925 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 4926 4927 ! Must have sufficient angular sampling 4928 if ((nth > 8).and.(nph > 8)) then 4929 rthe=r0 4930 do ijj=1,2 4931 theta=th(ijj) 4932 if (ijj==2) rthe=rs(1,1) 4933 do jj=1,nph 4934 phi=ph(jj) 4935 call timein(t1,wall) 4936 if (abs(rs(ijj,jj))<1.0d-12) then 4937 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 4938 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 4939 rs(ijj,jj)=rr 4940 end if 4941 rr=rs(ijj,jj) 4942 call timein(t2,wall) 4943 t2=t2-t1 4944 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 4945 write(untout,'(a,2i3,3E16.8)') '- ',ijj,jj,theta,phi,rr 4946 rthe=rs(ijj,jj) 4947 end do ! jj 4948 end do ! ijj 4949 4950 write(untout,*) 4951 4952 rthe=rs(2,1) 4953 do jj=1,2 4954 phi=ph(jj) 4955 if (jj==2) rthe=rs(2,2) 4956 do ijj=3,nth 4957 theta=th(ijj) 4958 t2=0.0 4959 call timein(t1,wall) 4960 if (abs(rs(ijj,jj))<1.0d-12) then 4961 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 4962 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 4963 rs(ijj,jj)=rr 4964 end if 4965 rr=rs(ijj,jj) 4966 call timein(t2,wall) 4967 t2=t2-t1 4968 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 4969 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 4970 rthe=rs(ijj,jj) 4971 end do ! ijj 4972 end do ! jj 4973 4974 write(untout,*) 4975 4976 rthe=rs(nth-1,2) 4977 do ijj=nth-1,nth 4978 theta=th(ijj) 4979 if (ijj==nth) rthe=rs(nth,2) 4980 do jj=3,nph 4981 phi=ph(jj) 4982 call timein(t1,wall) 4983 if (abs(rs(ijj,jj))<1.0d-12) then 4984 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 4985 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 4986 rs(ijj,jj)=rr 4987 end if 4988 rr=rs(ijj,jj) 4989 call timein(t2,wall) 4990 t2=t2-t1 4991 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 4992 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 4993 rthe=rs(ijj,jj) 4994 end do ! jj 4995 end do ! ijj 4996 4997 rthe=rs(2,nph-1) 4998 do jj=nph-1,nph 4999 phi=ph(jj) 5000 if (jj==nph) rthe=rs(2,nph) 5001 do ijj=3,nth-2 5002 theta=th(ijj) 5003 t2=0.0 5004 call timein(t1,wall) 5005 if (abs(rs(ijj,jj))<1.0d-12) then 5006 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5007 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5008 rs(ijj,jj)=rr 5009 end if 5010 rr=rs(ijj,jj) 5011 call timein(t2,wall) 5012 t2=t2-t1 5013 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5014 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 5015 rthe=rs(ijj,jj) 5016 end do ! ijj 5017 end do ! jj 5018 write(untout,*) 5019 5020 ! Complementary bands for boundary angles 5021 nn=int(real(nth)/1.4d1) 5022 if (nn > 1) then 5023 do ii=1,nn-1 5024 mm=int(nth/nn)*ii 5025 do kk=0,1 5026 mm=mm+kk 5027 theta=th(mm) 5028 rthe=rs(mm,2) 5029 do jj=3,nph-2 5030 phi=ph(jj) 5031 call timein(t1,wall) 5032 if (abs(rs(mm,jj))<1.0d-12) then 5033 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5034 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5035 rs(mm,jj)=rr 5036 end if 5037 rr=rs(mm,jj) 5038 call timein(t2,wall) 5039 t2=t2-t1 5040 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(mm)*wph(jj),t2 5041 write(untout,'(2i3,3E16.8)') mm,jj,theta,phi,rr 5042 rthe=rs(mm,jj) 5043 end do ! jj 5044 end do ! kk 5045 end do ! ii 5046 end if ! nn>1 5047 5048 write(untout,*) 5049 5050 nn=nint(real(nph)/1.2d1) 5051 if (nn > 1) then 5052 do ii=1,nn-1 5053 mm=int(nph/nn)*ii 5054 do kk=0,1 5055 mm=mm+kk 5056 phi=ph(mm) 5057 rthe=rs(2,mm) 5058 5059 do jj=3,nth-2 5060 theta=th(jj) 5061 call timein(t1,wall) 5062 if (abs(rs(jj,mm))<1.0d-12) then 5063 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5064 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5065 rs(jj,mm)=rr 5066 end if 5067 rr=rs(mm,jj) 5068 call timein(t2,wall) 5069 t2=t2-t1 5070 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(jj)*wph(mm),t2 5071 write(untout,'(2i3,3E16.8)') jj,mm,theta,phi,rr 5072 rthe=rs(jj,mm) 5073 end do ! jj 5074 5075 end do ! kk 5076 end do ! ii 5077 end if ! nn>1 5078 5079 end if ! sufficient sampling to determine boundary angles 5080 5081 write(untout,*) 5082 5083 ! DEBUG 5084 ! write(std_out,*)' surf : after boundary angles ' 5085 ! do ijj=1,nth 5086 ! do jj=1,nph 5087 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5088 ! end do 5089 ! end do 5090 ! ENDDEBUG 5091 5092 ! Output the complete Bader surface 5093 5094 write(untout,*) "The complete Bader surface:" 5095 write(untout,*) "===========================" 5096 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 5097 rthe0=r0 5098 srch=.true. 5099 5100 ! Write all the values 5101 5102 do ijj=1,nth 5103 theta=th(ijj) 5104 do jj=1,nph 5105 phi=ph(jj) 5106 rr=rs(ijj,jj) 5107 write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 5108 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5109 write(untout,'(a,2i3,3E16.8)') ' ',ijj,jj,theta,phi,rr 5110 if (rr < rsmin) rsmin=rr 5111 if (rr> rsmax) rsmax=rr 5112 end do ! jj 5113 end do ! ijj 5114 write(unts,'(2F15.10)') rsmin,rsmax 5115 write(untout,'(/," The minimal and maximal radii:",/,/," ",2F15.10)') rsmin,rsmax 5116 5117 ! DEBUG 5118 ! write(std_out,*)' surf : final output ' 5119 ! do ijj=1,nth 5120 ! do jj=1,nph 5121 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5122 ! end do 5123 ! end do 5124 ! ENDDEBUG 5125 5126 end if ! determination of the critical surface 5127 5128 call timein(ttsrf,wall) 5129 ttsrf=ttsrf-tt0 5130 5131 end subroutine surf
m_bader/vec_prod [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
vec_prod
FUNCTION
Vector product
INPUTS
vv,uu = vectors to compute vector product
OUTPUT
(return the value of the vector product)
SOURCE
5752 function vec_prod(uu,vv) 5753 5754 !Arguments ------------------------------------ 5755 !arrays 5756 real(dp) :: vec_prod(3) 5757 real(dp),intent(in) :: uu(3),vv(3) 5758 5759 !Local variables------------------------------- 5760 5761 ! ************************************************************************* 5762 5763 vec_prod(1)=uu(2)*vv(3)-vv(2)*uu(3) 5764 vec_prod(2)=uu(3)*vv(1)-vv(3)*uu(1) 5765 vec_prod(3)=uu(1)*vv(2)-vv(1)*uu(2) 5766 5767 end function vec_prod
m_bader/vgh_rho [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
vgh_rho
FUNCTION
The general procedure to obtain the value, the gradient and the hessian of the density of electrons in the point vv (in cart.coord). WARNING This file does not follow the ABINIT coding rules (yet)
INPUTS
vv(3)=position
chs 1 only valence density
2 only core density
0 total density
-2 iat, ipos are nulify and ignored
-1 iat,ipos = index of atom if vv is inside
the "core sphere (rminl)", 0 otherwise
OUTPUT
rho,grho(3),hrho(3,3) - density, gradient of density, hessian of density
(cart. coord)
iat, ipos - index of the nearest atom (except chs < 0 see above)
rdmin - the distance to the nearest atom
SIDE EFFECTS
This routine also works on the data contained in the defs_aimprom and defs_aimfields modules
SOURCE
5165 subroutine vgh_rho(vv,rho,grho,hrho,rdmin,iat,ipos,chs) 5166 5167 !Arguments ------------------------------------ 5168 !scalars 5169 integer,intent(in) :: chs 5170 integer,intent(inout) :: iat,ipos 5171 real(dp),intent(out) :: rdmin,rho 5172 !arrays 5173 real(dp),intent(in) :: vv(3) 5174 real(dp),intent(out) :: grho(3),hrho(3,3) 5175 5176 !Local variables ------------------------------ 5177 !scalars 5178 integer :: ii,inmax,inmin,inx,jj,kk,ll,nn,oii,omm,onn 5179 integer :: selct 5180 ! real(dp),save :: cumul_cpu=0.0_dp,cumul_cpu_old=0.0_dp 5181 real(dp),save :: tcpui,tcpuo,twalli 5182 real(dp),save :: twallo 5183 real(dp) :: aa,bb,cc,cgrad1_rr_inv,coeff,dd,rr,rr2,rr_inv 5184 real(dp) :: rrad2_nn,rrad_nn,ss,uu,uu_inv,val,vt1,vt2,vt3,vw1,vw2 5185 ! real(dp) :: ss_inv 5186 real(dp) :: vw3 5187 !arrays 5188 integer :: indx(3),inii(4,3) 5189 real(dp) :: cgrad(3),ches(3,3),cof(2,3),ddstar(6),ddu(2),grd(4) 5190 real(dp) :: hh(4,2),hrh(2),lder(4),pom2sq(2,3),pomsq(2,3) 5191 real(dp) :: rhstar(6),sqder(6,4),sqvlr(6,4),trsf(3,3),xx(3) 5192 real(dp),pointer :: ptddx(:,:,:),ptddy(:,:,:),ptddz(:,:,:),ptrho(:,:,:) 5193 5194 !************************************************************************ 5195 tcpui=0.0_dp 5196 tcpuo=0.0_dp 5197 twalli=0.0_dp 5198 twallo=0.0_dp 5199 5200 nullify(ptddx,ptddy,ptddz,ptrho) 5201 5202 selct=chs 5203 5204 if (selct/=2) then 5205 5206 ! call timein(tcpui,twalli) 5207 5208 ! TRANSFORMATION TO THE REDUCED COORD. 5209 5210 xx(:)=vv(:) 5211 call bschg1(xx,-1) 5212 5213 ! call timein(tcpuo,twallo) 5214 ! cumul_cpu=cumul_cpu+(tcpuo-tcpui) 5215 5216 ! REDUCTION TO THE PRIMITIVE CELL 5217 5218 do ii=1,3 5219 if (xx(ii) >= one-tol12 ) then 5220 xx(ii)=xx(ii)-aint(xx(ii)) 5221 elseif (xx(ii) < -tol12 ) then 5222 xx(ii)=xx(ii)-floor(xx(ii)) 5223 end if 5224 end do 5225 5226 5227 ! DETERMINATION OF THE INDEX IN THE GRID 5228 5229 do ii=1,3 5230 indx(ii)=aint(xx(ii)*ngfft(ii)) 5231 bb=(xx(ii)-indx(ii)*dix(ii))*ngfft(ii) 5232 if (indx(ii)==ngfft(ii)) then 5233 indx(ii)=1 5234 xx(ii)=0._dp 5235 else 5236 indx(ii)=indx(ii)+1 5237 end if 5238 5239 ! Explicit handling to avoid numeric problems 5240 5241 if (bb > 1._dp+tol12 ) then 5242 cof(1,ii)=0._dp 5243 cof(2,ii)=1._dp 5244 elseif (bb < -tol12 ) then 5245 cof(1,ii)=1._dp 5246 cof(2,ii)=0._dp 5247 else 5248 cof(1,ii)=1._dp-bb 5249 cof(2,ii)=bb 5250 end if 5251 end do 5252 5253 ! 3D INTERPOLATION OF THE VALENCE DENSITY 5254 5255 ! determination of the values of density and of its second derivative 5256 ! at the "star" = constructed at vv with primitive directions 5257 ! To interpolation the values at the faces of the grid cell are needed 5258 5259 rhstar(:)=0._dp 5260 sqder(:,:)=0._dp 5261 sqvlr(:,:)=0._dp 5262 ddstar(:)=0._dp 5263 pomsq(:,:)=0._dp 5264 pom2sq(:,:)=0._dp 5265 5266 oii=1; onn=1; omm=1 5267 if (indx(1)==ngfft(1)) oii=1-ngfft(1) 5268 if (indx(2)==ngfft(2)) onn=1-ngfft(2) 5269 if (indx(3)==ngfft(3)) omm=1-ngfft(3) 5270 5271 ! the values in the corners of the grid cell 5272 5273 ptddx=>ddx(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5274 ptddy=>ddy(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5275 ptddz=>ddz(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5276 ptrho=>dvl(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5277 5278 ! the coefficients for spline interpolation of density and its derivation 5279 do ii=1,3 5280 do jj=1,2 5281 pomsq(jj,ii)=(cof(jj,ii)*cof(jj,ii)*cof(jj,ii)-cof(jj,ii))/6._dp*dix(ii)*dix(ii) 5282 pom2sq(jj,ii)=(3._dp*cof(jj,ii)*cof(jj,ii)-1._dp)/6._dp*dix(ii) 5283 if (jj==1) pom2sq(jj,ii)=-pom2sq(jj,ii) 5284 end do 5285 end do 5286 5287 5288 do ii=1,2 5289 do jj=1,2 5290 do kk=1,2 5291 ddstar(ii)=ddstar(ii)+cof(jj,2)*cof(kk,3)*ptddx(ii,jj,kk) 5292 ddstar(ii+2)=ddstar(ii+2)+cof(jj,3)*cof(kk,1)*ptddy(kk,ii,jj) 5293 ddstar(ii+4)=ddstar(ii+4)+cof(jj,1)*cof(kk,2)*ptddz(jj,kk,ii) 5294 sqder(ii,jj)=sqder(ii,jj)+cof(kk,2)*ptddz(ii,kk,jj) 5295 sqder(ii,jj+2)=sqder(ii,jj+2)+cof(kk,3)*ptddy(ii,jj,kk) 5296 sqder(ii+2,jj)=sqder(ii+2,jj)+cof(kk,3)*ptddx(jj,ii,kk) 5297 sqder(ii+2,jj+2)=sqder(ii+2,jj+2)+cof(kk,1)*ptddz(kk,ii,jj) 5298 sqder(ii+4,jj)=sqder(ii+4,jj)+cof(kk,1)*ptddy(kk,jj,ii) 5299 sqder(ii+4,jj+2)=sqder(ii+4,jj+2)+cof(kk,2)*ptddx(jj,kk,ii) 5300 sqvlr(ii,jj)=sqvlr(ii,jj)+cof(kk,2)*ptrho(ii,kk,jj)+pomsq(kk,2)*ptddy(ii,kk,jj) 5301 sqvlr(ii,jj+2)=sqvlr(ii,jj+2)+cof(kk,3)*ptrho(ii,jj,kk)+pomsq(kk,3)*ptddz(ii,jj,kk) 5302 sqvlr(ii+2,jj+2)=sqvlr(ii+2,jj+2)+cof(kk,1)*ptrho(kk,ii,jj)+pomsq(kk,1)*ptddx(kk,ii,jj) 5303 end do 5304 end do 5305 end do 5306 5307 do ii=1,2 5308 do jj=1,2 5309 sqvlr(ii+2,jj)=sqvlr(jj,ii+2) 5310 sqvlr(ii+4,jj)=sqvlr(jj+2,ii+2) 5311 sqvlr(ii+4,jj+2)=sqvlr(jj,ii) 5312 end do 5313 end do 5314 5315 do ii=1,2 5316 do jj=1,2 5317 rhstar(ii)=rhstar(ii)+cof(jj,3)*sqvlr(ii,jj)+pomsq(jj,3)*sqder(ii,jj)+& 5318 & cof(jj,2)*sqvlr(ii,jj+2)+pomsq(jj,2)*sqder(ii,jj+2) 5319 rhstar(ii+2)=rhstar(ii+2)+cof(jj,1)*sqvlr(ii+2,jj)+pomsq(jj,1)*sqder(ii+2,jj)+& 5320 & cof(jj,3)*sqvlr(ii+2,jj+2)+pomsq(jj,3)*sqder(ii+2,jj+2) 5321 rhstar(ii+4)=rhstar(ii+4)+cof(jj,2)*sqvlr(ii+4,jj)+pomsq(jj,2)*sqder(ii+4,jj)+& 5322 & cof(jj,1)*sqvlr(ii+4,jj+2)+pomsq(jj,1)*sqder(ii+4,jj+2) 5323 end do 5324 end do 5325 rhstar(:)=rhstar(:)/2._dp 5326 5327 rho=0._dp 5328 grho(:)=0._dp 5329 hrho(:,:)=0._dp 5330 kk=1; nn=1 5331 do ii=1,5,2 5332 do jj=1,2 5333 nn=-nn 5334 rho=rho+cof(jj,kk)*rhstar(ii+jj-1)+pomsq(jj,kk)*ddstar(ii+jj-1) 5335 grho(kk)=grho(kk)+pom2sq(jj,kk)*ddstar(ii+jj-1) 5336 hrho(kk,kk)=hrho(kk,kk)+cof(jj,kk)*ddstar(ii+jj-1) 5337 grho(kk)=grho(kk)+nn*rhstar(ii+jj-1)/dix(kk) 5338 end do 5339 kk=kk+1 5340 end do 5341 rho=rho/3._dp 5342 5343 ! Off-diagonal elements of the hessian 5344 5345 ! for the speed reasons the polynomial interpolation 5346 ! for second derivation fields is used in this case 5347 ! but the last step is always done by spline interpolation. 5348 5349 5350 do ii=1,3 5351 do jj=-1,2 5352 inii(jj+2,ii)=indx(ii)+jj 5353 if (inii(jj+2,ii) < 1) inii(jj+2,ii)=inii(jj+2,ii)+ngfft(ii) 5354 if (inii(jj+2,ii) > ngfft(ii)) inii(jj+2,ii)=inii(jj+2,ii)-ngfft(ii) 5355 end do 5356 end do 5357 5358 ! Not very nice 5359 5360 do ii=1,3 5361 select case (ii) 5362 case (1) 5363 do jj=1,4 5364 ddu(1)=cof(1,2)*ddz(inii(jj,1),inii(2,2),inii(2,3))+cof(2,2)*ddz(inii(jj,1),inii(3,2),inii(2,3)) 5365 ddu(2)=cof(1,2)*ddz(inii(jj,1),inii(2,2),inii(3,3))+cof(2,2)*ddz(inii(jj,1),inii(3,2),inii(3,3)) 5366 hrh(1)=cof(1,2)*dvl(inii(jj,1),inii(2,2),inii(2,3))+cof(2,2)*dvl(inii(jj,1),inii(3,2),inii(2,3))+& 5367 & pomsq(1,2)*ddy(inii(jj,1),inii(2,2),inii(2,3))+pomsq(2,2)*ddy(inii(jj,1),inii(3,2),inii(2,3)) 5368 hrh(2)=cof(1,2)*dvl(inii(jj,1),inii(2,2),inii(3,3))+cof(2,2)*dvl(inii(jj,1),inii(3,2),inii(3,3))+& 5369 & pomsq(1,2)*ddy(inii(jj,1),inii(2,2),inii(3,3))+pomsq(2,2)*ddy(inii(jj,1),inii(3,2),inii(3,3)) 5370 hh(jj,2)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 5371 5372 ddu(1)=cof(1,3)*ddy(inii(jj,1),inii(2,2),inii(2,3))+cof(2,3)*ddy(inii(jj,1),inii(2,2),inii(3,3)) 5373 ddu(2)=cof(1,3)*ddy(inii(jj,1),inii(3,2),inii(2,3))+cof(2,3)*ddy(inii(jj,1),inii(3,2),inii(3,3)) 5374 hrh(1)=cof(1,3)*dvl(inii(jj,1),inii(2,2),inii(2,3))+cof(2,3)*dvl(inii(jj,1),inii(2,2),inii(3,3))+& 5375 & pomsq(1,3)*ddz(inii(jj,1),inii(2,2),inii(2,3))+pomsq(2,3)*ddz(inii(jj,1),inii(2,2),inii(3,3)) 5376 hrh(2)=cof(1,3)*dvl(inii(jj,1),inii(3,2),inii(2,3))+cof(2,3)*dvl(inii(jj,1),inii(3,2),inii(3,3))+& 5377 & pomsq(1,3)*ddz(inii(jj,1),inii(3,2),inii(2,3))+pomsq(2,3)*ddz(inii(jj,1),inii(3,2),inii(3,3)) 5378 hh(jj,1)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 5379 end do 5380 case (2) 5381 do jj=1,4 5382 ddu(1)=cof(1,3)*ddx(inii(2,1),inii(jj,2),inii(2,3))+cof(2,3)*ddx(inii(2,1),inii(jj,2),inii(3,3)) 5383 ddu(2)=cof(1,3)*ddx(inii(3,1),inii(jj,2),inii(2,3))+cof(2,3)*ddx(inii(3,1),inii(jj,2),inii(3,3)) 5384 hrh(1)=cof(1,3)*dvl(inii(2,1),inii(jj,2),inii(2,3))+cof(2,3)*dvl(inii(2,1),inii(jj,2),inii(3,3))+& 5385 & pomsq(1,3)*ddz(inii(2,1),inii(jj,2),inii(2,3))+pomsq(2,3)*ddz(inii(2,1),inii(jj,2),inii(3,3)) 5386 hrh(2)=cof(1,3)*dvl(inii(3,1),inii(jj,2),inii(2,3))+cof(2,3)*dvl(inii(3,1),inii(jj,2),inii(3,3))+& 5387 & pomsq(1,3)*ddz(inii(3,1),inii(jj,2),inii(2,3))+pomsq(2,3)*ddz(inii(3,1),inii(jj,2),inii(3,3)) 5388 hh(jj,2)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 5389 5390 ddu(1)=cof(1,1)*ddz(inii(2,1),inii(jj,2),inii(2,3))+cof(2,1)*ddz(inii(3,1),inii(jj,2),inii(2,3)) 5391 ddu(2)=cof(1,1)*ddz(inii(2,1),inii(jj,2),inii(3,3))+cof(2,1)*ddz(inii(3,1),inii(jj,2),inii(3,3)) 5392 hrh(1)=cof(1,1)*dvl(inii(2,1),inii(jj,2),inii(2,3))+cof(2,1)*dvl(inii(3,1),inii(jj,2),inii(2,3))+& 5393 & pomsq(1,1)*ddx(inii(2,1),inii(jj,2),inii(2,3))+pomsq(2,1)*ddx(inii(3,1),inii(jj,2),inii(2,3)) 5394 hrh(2)=cof(1,1)*dvl(inii(2,1),inii(jj,2),inii(3,3))+cof(2,1)*dvl(inii(3,1),inii(jj,2),inii(3,3))+& 5395 & pomsq(1,1)*ddx(inii(2,1),inii(jj,2),inii(3,3))+pomsq(2,1)*ddx(inii(3,1),inii(jj,2),inii(3,3)) 5396 hh(jj,1)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 5397 end do 5398 case (3) 5399 do jj=1,4 5400 ddu(1)=cof(1,1)*ddy(inii(2,1),inii(2,2),inii(jj,3))+cof(2,1)*ddy(inii(3,1),inii(2,2),inii(jj,3)) 5401 ddu(2)=cof(1,1)*ddy(inii(2,1),inii(3,2),inii(jj,3))+cof(2,1)*ddy(inii(3,1),inii(3,2),inii(jj,3)) 5402 hrh(1)=cof(1,1)*dvl(inii(2,1),inii(2,2),inii(jj,3))+cof(2,1)*dvl(inii(3,1),inii(2,2),inii(jj,3))+& 5403 & pomsq(1,1)*ddx(inii(2,1),inii(2,2),inii(jj,3))+pomsq(2,1)*ddx(inii(3,1),inii(2,2),inii(jj,3)) 5404 hrh(2)=cof(1,1)*dvl(inii(2,1),inii(3,2),inii(jj,3))+cof(2,1)*dvl(inii(3,1),inii(3,2),inii(jj,3))+& 5405 & pomsq(1,1)*ddx(inii(2,1),inii(3,2),inii(jj,3))+pomsq(2,1)*ddx(inii(3,1),inii(3,2),inii(jj,3)) 5406 hh(jj,2)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 5407 5408 ddu(1)=cof(1,2)*ddx(inii(2,1),inii(2,2),inii(jj,3))+cof(2,2)*ddx(inii(2,1),inii(3,2),inii(jj,3)) 5409 ddu(2)=cof(1,2)*ddx(inii(3,1),inii(2,2),inii(jj,3))+cof(2,2)*ddx(inii(3,1),inii(3,2),inii(jj,3)) 5410 hrh(1)=cof(1,2)*dvl(inii(2,1),inii(2,2),inii(jj,3))+cof(2,2)*dvl(inii(2,1),inii(3,2),inii(jj,3))+& 5411 & pomsq(1,2)*ddy(inii(2,1),inii(2,2),inii(jj,3))+pomsq(2,2)*ddy(inii(2,1),inii(3,2),inii(jj,3)) 5412 hrh(2)=cof(1,2)*dvl(inii(3,1),inii(2,2),inii(jj,3))+cof(2,2)*dvl(inii(3,1),inii(3,2),inii(jj,3))+& 5413 & pomsq(1,2)*ddy(inii(3,1),inii(2,2),inii(jj,3))+pomsq(2,2)*ddy(inii(3,1),inii(3,2),inii(jj,3)) 5414 hh(jj,1)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 5415 end do 5416 end select 5417 do jj=-2,1 5418 grd(jj+3)=(indx(ii)+jj)*dix(ii) 5419 end do 5420 5421 ! write(std_out,'("hh: ",/,4F16.8,/,4F16.8)') ((hh(kk,jj),kk=1,4),jj=1,2) 5422 ! write(std_out,'("grad: ",3F16.8)') (grho(kk),kk=1,3) 5423 ! write(std_out,'("dix: ",3F16.8)') (dix(kk),kk=1,3) 5424 ! write(std_out,'("grd: ",4F16.8)') (grd(kk),kk=1,4) 5425 ! write(std_out,'("inii: ",4I4)') (inii(kk,ii),kk=1,4) 5426 5427 do jj=1,2 5428 5429 ! polynomial interpolation 5430 5431 do kk=1,3 5432 do ll=4,kk+1,-1 5433 hh(ll,jj)=(hh(ll,jj)-hh(ll-1,jj))/(grd(ll)-grd(ll-1)) 5434 end do 5435 end do 5436 lder(4)=hh(4,jj) 5437 do kk=3,1,-1 5438 lder(kk)=hh(kk,jj)+(xx(ii)-grd(kk))*lder(kk+1) 5439 end do 5440 do kk=1,2 5441 do ll=3,kk+1,-1 5442 lder(ll)=lder(ll)+(xx(ii)-grd(ll-kk))*lder(ll+1) 5443 end do 5444 end do 5445 nn=ii+jj 5446 if (nn > 3) nn=nn-3 5447 hrho(ii,nn)=hrho(ii,nn)+lder(2) 5448 hrho(nn,ii)=hrho(nn,ii)+lder(2) 5449 end do 5450 end do 5451 5452 ! averaging of the mixed derivations obtained in different order 5453 5454 do ii=1,3 5455 do jj=1,3 5456 if (ii /= jj) hrho(ii,jj)=hrho(ii,jj)/2._dp 5457 end do 5458 end do 5459 5460 5461 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 5462 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 5463 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 5464 ! stop 5465 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 5466 ! write(std_out,'(":GRAD pred tr ",3F16.8)') grho 5467 ! write(std_out,'(":HESSIAN pred tr",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5468 5469 5470 ! Transformation back to Cart. coordonnes 5471 5472 call bschg1(grho,2) 5473 call bschg2(hrho,2) 5474 5475 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 5476 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 5477 ! stop 5478 5479 nullify(ptddx,ptddy,ptddz,ptrho) 5480 5481 if (selct==1) return 5482 5483 end if 5484 5485 !write(51,'(":GRADv ",3F16.8)') grho 5486 !write(52,'(":LAPv ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 5487 !write(52,'(":HESNv ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5488 5489 !INTERPOLATION OF THE CORE DENSITY 5490 5491 if (selct/=1) then 5492 5493 if (selct==2) then 5494 grho(:)=0._dp 5495 hrho(:,:)=0._dp 5496 rho=0._dp 5497 end if 5498 5499 ! SEARCH OF THE NEIGHBOUR ATOMS 5500 5501 if (selct /= -2) then 5502 iat=0 5503 ipos=0 5504 end if 5505 rdmin=20._dp 5506 5507 do jj=1,natom 5508 nn=typat(jj) 5509 rrad_nn=rrad(corlim(nn),nn) 5510 rrad2_nn=rrad_nn*rrad_nn 5511 vw1=vv(1)-xatm(1,jj) 5512 vw2=vv(2)-xatm(2,jj) 5513 vw3=vv(3)-xatm(3,jj) 5514 5515 do kk=1,nnpos 5516 5517 vt1=vw1-atp(1,kk) 5518 vt2=vw2-atp(2,kk) 5519 vt3=vw3-atp(3,kk) 5520 rr2=vt1*vt1+vt2*vt2+vt3*vt3 5521 5522 ! rr=vnorm(vt,0) 5523 5524 ! Only contribution > rhocormin (adhoc.f90) are considered 5525 5526 if (rr2 < rrad2_nn .and.(.not.((selct==-2).and.(iat==jj).and.(ipos==kk)))) then 5527 ! if (rr /= 0.0_dp) then ! XG020629 : never test a real number against zero (not portable) 5528 if (rr2 > 1.0d-28) then ! SEARCH INDEX 5529 5530 rr=sqrt(rr2) 5531 rr_inv=1.0_dp/rr 5532 5533 if (rr < rrad(1,nn)) then 5534 inx=-1 5535 elseif (rr > rrad(ndat(nn),nn)) then 5536 inx=ndat(nn) 5537 else 5538 ! Find the index of the radius by bissection 5539 inmin=1 5540 inmax=ndat(nn) 5541 inx=1 5542 do 5543 if(inmax-inmin==1)exit 5544 inx=(inmin+inmax)/2 5545 if(rr>=rrad(inx,nn))then 5546 inmin=inx 5547 else 5548 inmax=inx 5549 end if 5550 end do 5551 inx=inmin 5552 5553 ! XG020629 : old coding, slower 5554 ! inx=0 5555 ! do while (rr >= rrad(inx+1,nn)) 5556 ! inx=inx+1 5557 ! end do 5558 5559 end if 5560 5561 ! Transformation matrix radial -> cart. coord 5562 ss=sqrt(vt1*vt1+vt2*vt2) 5563 ! if (ss /=0._dp) then ! XG020629 : never test a real number against zero (not portable) 5564 if (ss*ss > 1.0d-28) then ! ss non-zero 5565 ! XG020629 : very strange : only trsf(:,1) is needed in what follows ? ! 5566 ! ss_inv=1.0_dp/ss 5567 trsf(1,1)=vt1*rr_inv 5568 ! trsf(1,2)=-vt2*ss_inv 5569 ! trsf(1,3)=vt3*vt1*rr_inv*ss_inv 5570 trsf(2,1)=vt2*rr_inv 5571 ! trsf(2,2)=vt1*ss_inv 5572 ! trsf(2,3)=vt3*vt2*rr_inv*ss_inv 5573 trsf(3,1)=vt3*rr_inv 5574 ! trsf(3,2)=0._dp 5575 ! trsf(3,3)=-ss*rr_inv 5576 ! XG020629 Not needed 5577 ! do ii=1,3 5578 ! do ll=1,3 5579 ! ches(ii,ll)=0._dp 5580 ! end do 5581 ! cgrad(ii)=0._dp 5582 ! end do 5583 else ! ss zero 5584 do ii=1,3 5585 do ll=1,3 5586 trsf(ii,ll)=0._dp 5587 end do 5588 trsf(ii,4-ii)=1._dp 5589 end do 5590 end if ! ss zero or non-zero 5591 5592 if (inx == -1) then ! LEFT EXTRAPOLATION y=a*x^2+b (a<0) 5593 val=sp2(1,nn)*0.5_dp*rr*rr/rrad(1,nn)+crho(1,nn)-sp2(1,nn)*rrad(1,nn)*0.5_dp 5594 cgrad(1)=sp2(1,nn)*rr/rrad(1,nn) 5595 ches(1,1)=sp2(1,nn)/rrad(1,nn) 5596 elseif (inx == ndat(nn) ) then ! RIGHT EXTRAPOLATION y=a*exp(b*x) 5597 val=rrad(ndat(nn),nn)*exp(sp2(ndat(nn),nn)*(rr-rrad(ndat(nn),nn))/crho(ndat(nn),nn)) 5598 cgrad(1)=val*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 5599 ches(1,1)=cgrad(1)*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 5600 else ! INTERPOLATION 5601 uu=rrad(inx+1,nn)-rrad(inx,nn) 5602 uu_inv=1.0_dp/uu 5603 aa=(rrad(inx+1,nn)-rr)*uu_inv 5604 bb=(rr-rrad(inx,nn))*uu_inv 5605 cc=(aa*aa*aa-aa)*uu*uu*0.16666666666666666_dp 5606 dd=(bb*bb*bb-bb)*uu*uu*0.16666666666666666_dp 5607 val=aa*crho(inx,nn)+bb*crho(inx+1,nn)+cc*sp3(inx,nn)+dd*sp3(inx+1,nn) 5608 cgrad(1)=(crho(inx+1,nn)-crho(inx,nn))*uu_inv& 5609 & -(3._dp*aa*aa-1._dp)*uu*0.16666666666666666_dp*sp3(inx,nn)+& 5610 & (3._dp*bb*bb-1._dp)*uu*0.16666666666666666_dp*sp3(inx+1,nn) 5611 ches(1,1)=aa*sp3(inx,nn)+bb*sp3(inx+1,nn) 5612 5613 end if ! TRANSFORMATION TO CARTEZ. COORD. 5614 5615 cgrad1_rr_inv=cgrad(1)*rr_inv 5616 coeff=(ches(1,1)-cgrad1_rr_inv)*rr_inv*rr_inv 5617 cgrad(3)=trsf(3,1)*cgrad(1) 5618 cgrad(2)=trsf(2,1)*cgrad(1) 5619 cgrad(1)=trsf(1,1)*cgrad(1) 5620 ches(1,1)=coeff*vt1*vt1+cgrad1_rr_inv 5621 ches(2,2)=coeff*vt2*vt2+cgrad1_rr_inv 5622 ches(3,3)=coeff*vt3*vt3+cgrad1_rr_inv 5623 ches(1,2)=coeff*vt1*vt2 ; ches(2,1)=coeff*vt1*vt2 5624 ches(1,3)=coeff*vt1*vt3 ; ches(3,1)=coeff*vt1*vt3 5625 ches(2,3)=coeff*vt2*vt3 ; ches(3,2)=coeff*vt2*vt3 5626 5627 else ! case rr==0 5628 5629 val=crho(1,nn)-sp2(1,nn)*rrad(1,nn)/2._dp 5630 do ii=1,3 5631 do ll=1,3 5632 ches(ii,ll)=0._dp 5633 end do 5634 cgrad(ii)=0._dp 5635 ches(ii,ii)=sp2(1,nn)/rrad(1,nn) 5636 end do 5637 5638 end if ! rr>0 or rr==0 5639 5640 do ii=1,3 5641 do ll=1,3 5642 hrho(ii,ll)=hrho(ii,ll)+ches(ii,ll) 5643 end do 5644 grho(ii)=grho(ii)+cgrad(ii) 5645 end do 5646 rho=rho+val 5647 5648 end if ! rr2< rrad_nn*rrad_nn 5649 5650 if (selct==-1) then 5651 if (rr2 < rminl(jj)*rminl(jj) ) then 5652 iat=jj 5653 ipos=kk 5654 rdmin=sqrt(rr2) 5655 end if 5656 elseif (selct==-2) then 5657 cycle 5658 else 5659 if (rr2 < rdmin*rdmin) then 5660 iat=jj 5661 ipos=kk 5662 rdmin=sqrt(rr2) 5663 end if 5664 end if 5665 5666 end do 5667 end do 5668 5669 end if 5670 5671 !write(51,'(":GRADt ",3F16.8)') grho 5672 !write(52,'(":LAPt ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 5673 !write(52,'(":HESNt ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5674 5675 !if(abs(cumul_cpu-cumul_cpu_old)>0.499)then 5676 !write(std_out,'(a,f7.1)' )' vgh_rho : cumul_cpu=',cumul_cpu 5677 !cumul_cpu_old=cumul_cpu 5678 !end if 5679 5680 end subroutine vgh_rho
m_bader/vnorm [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
vnorm
FUNCTION
Default declarations, and interfaces for the aim.f utility.
INPUTS
vector norm ->dir==1: vector in reduced coordinates
dir==0: vector in cartes. coordinates
OUTPUT
(see side effects)
SIDE EFFECTS
vv = on entry, vector to normalized
= on exit, normalized vector
SOURCE
5703 function vnorm(vv,dir) 5704 5705 !Arguments ------------------------------------ 5706 !scalars 5707 integer,intent(in) :: dir 5708 real(dp) :: vnorm 5709 !arrays 5710 real(dp),intent(in) :: vv(3) 5711 5712 !Local variables------------------------------- 5713 !scalars 5714 integer :: ii 5715 !arrays 5716 real(dp) :: vt(3) 5717 5718 ! ************************************************************************* 5719 5720 vnorm=zero 5721 if (dir==1) then 5722 do ii=1,3 5723 vt(ii)=rprimd(ii,1)*vv(1)+rprimd(ii,2)*vv(2)+rprimd(ii,3)*vv(3) 5724 vnorm=vnorm+vt(ii)*vt(ii) 5725 end do 5726 elseif (dir==0) then 5727 do ii=1,3 5728 vnorm=vnorm+vv(ii)*vv(ii) 5729 end do 5730 else 5731 ABI_ERROR('vnorm calcul') 5732 end if 5733 vnorm=sqrt(vnorm) 5734 end function vnorm