TABLE OF CONTENTS
- ABINIT/m_bader
- defs_aimprom/aim_shutdown
- m_bader/addout
- m_bader/adini
- 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-2018 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 .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_bader 28 29 use defs_basis 30 use defs_abitypes 31 use defs_datatypes 32 use m_errors 33 use m_abicore 34 use m_xmpi 35 use m_sort 36 use m_hdr 37 use m_splines 38 #ifdef HAVE_NETCDF 39 use netcdf 40 #endif 41 42 use m_time, only : timein 43 use m_geometry, only : metric 44 use m_parser, only : inread 45 use m_numeric_tools, only : coeffs_gausslegint 46 use m_hide_lapack, only : jacobi, lubksb, ludcmp 47 48 implicit none 49 50 !private 51 public
defs_aimprom/aim_shutdown [ Functions ]
NAME
aim_shutdown
FUNCTION
Free memory allocated in the module. Close units. Mainly used to pass the abirules
PARENTS
aim
CHILDREN
SOURCE
166 subroutine aim_shutdown() 167 168 169 !This section has been created automatically by the script Abilint (TD). 170 !Do not modify the following lines by hand. 171 #undef ABI_FUNC 172 #define ABI_FUNC 'aim_shutdown' 173 !End of the abilint section 174 175 implicit none 176 177 !This section has been created automatically by the script Abilint (TD). 178 !Do not modify the following lines by hand. 179 #undef ABI_FUNC 180 #define ABI_FUNC 'aim_shutdown' 181 !End of the abilint section 182 183 !Local variables------------------------------- 184 integer :: ii 185 logical :: is_open 186 integer :: all_units(12) 187 188 ! ********************************************************************* 189 190 !if (allocated(typat)) then 191 ! ABI_FREE(typat) 192 !end if 193 !if (allocated(corlim)) then 194 ! ABI_FREE(corlim) 195 !end if 196 !if (allocated(xred)) then 197 ! ABI_FREE(xred) 198 !end if 199 !if (allocated(xatm)) then 200 ! ABI_FREE(rminl) 201 !end if 202 203 all_units(:) = [unt0,unto,unt,untc,unts,untd,untl,untg,unta,untad,untp,untout] 204 do ii=1,size(all_units) 205 inquire(unit=all_units(ii), opened=is_open) 206 if (is_open) close(all_units(ii)) 207 end do 208 209 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
PARENTS
drvaim
CHILDREN
vgh_rho
SOURCE
741 subroutine addout(aim_dtset) 742 743 744 !This section has been created automatically by the script Abilint (TD). 745 !Do not modify the following lines by hand. 746 #undef ABI_FUNC 747 #define ABI_FUNC 'addout' 748 !End of the abilint section 749 750 implicit none 751 752 !Arguments ------------------------------------ 753 !scalars 754 type(aim_dataset_type),intent(in) :: aim_dtset 755 756 !Local variables ------------------------------ 757 !scalars 758 integer :: cod,dims,iat,ii,ipos,jj,nn,tgrd 759 real(dp) :: alfa,rho,rr,xx,yy 760 !arrays 761 real(dp) :: grho(3),hrho(3,3),orig(3),vv(3) 762 real(dp),allocatable :: dfld(:),lfld(:),nr(:),stp(:),uu(:,:) 763 764 !************************************************************************ 765 orig(:)=aim_dtset%vpts(:,1) 766 if (aim_dtset%denout > 0) then 767 dims=aim_dtset%denout 768 elseif (aim_dtset%lapout > 0) then 769 dims=aim_dtset%lapout 770 end if 771 772 select case (aim_dtset%dltyp) 773 case (1) 774 cod=1 775 case (2) 776 cod=2 777 case default 778 cod=0 779 end select 780 781 ABI_ALLOCATE(uu,(3,dims)) 782 ABI_ALLOCATE(nr,(dims)) 783 ABI_ALLOCATE(stp,(dims)) 784 785 write(std_out,*) 'grid:', aim_dtset%ngrid(1:dims) 786 write(std_out,*) 'kod :', cod 787 tgrd=1 788 do ii=1,dims 789 tgrd=tgrd*aim_dtset%ngrid(ii) 790 uu(:,ii)=aim_dtset%vpts(:,ii+1)-aim_dtset%vpts(:,1) 791 nr(ii)=vnorm(uu(:,ii),0) 792 stp(ii)=nr(ii)/(aim_dtset%ngrid(ii)-1) 793 uu(:,ii)=uu(:,ii)/nr(ii) 794 end do 795 write(std_out,*) 'tgrd :', tgrd 796 do ii=1,dims 797 write(std_out,*) 'uu :', uu(1:3,ii) 798 end do 799 800 if (aim_dtset%denout > 0) then 801 ABI_ALLOCATE(dfld,(tgrd+1)) 802 dfld(:)=0._dp 803 end if 804 if (aim_dtset%lapout > 0) then 805 ABI_ALLOCATE(lfld,(tgrd+1)) 806 end if 807 808 select case (dims) 809 case (1) 810 nn=0 811 do ii=0,aim_dtset%ngrid(1)-1 812 nn=nn+1 813 vv(:)=orig(:)+ii*stp(1)*uu(:,1) 814 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,cod) 815 if (aim_dtset%denout > 0) dfld(nn)=rho 816 if (aim_dtset%lapout > 0) lfld(nn)=hrho(1,1)+hrho(2,2)+hrho(3,3) 817 end do 818 if (aim_dtset%denout==1) then 819 do ii=0,aim_dtset%ngrid(1)-1 820 xx=ii*stp(1) 821 write(untd,'(2E16.8)') xx, dfld(ii+1) 822 end do 823 end if 824 if (aim_dtset%lapout==1) then 825 do ii=0,aim_dtset%ngrid(1)-1 826 xx=ii*stp(1) 827 write(untl,'(2E16.8)') xx, lfld(ii+1) 828 end do 829 end if 830 case (2) 831 nn=0 832 alfa=dot_product(uu(:,1),uu(:,2)) 833 alfa=acos(alfa) 834 do ii=0,aim_dtset%ngrid(2)-1 835 do jj=0,aim_dtset%ngrid(1)-1 836 nn=nn+1 837 vv(:)=orig(:)+jj*uu(:,2)*stp(2)+ii*stp(1)*uu(:,1) 838 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,cod) 839 if (aim_dtset%denout > 0) dfld(nn)=rho 840 if (aim_dtset%lapout > 0) lfld(nn)=hrho(1,1)+hrho(2,2)+hrho(3,3) 841 end do 842 end do 843 write(std_out,*) 'generace hotova', nn 844 nn=0 845 if (aim_dtset%denout==2) then 846 do ii=0,aim_dtset%ngrid(2)-1 847 do jj=0,aim_dtset%ngrid(1)-1 848 nn=nn+1 849 xx=jj*stp(1)+cos(alfa)*ii*stp(2) 850 yy=sin(alfa)*ii*stp(2) 851 write(untd,'(3E16.8)') xx, yy, dfld(nn) 852 end do 853 write(untd,*) ' ' 854 end do 855 end if 856 nn=0 857 if (aim_dtset%lapout==2) then 858 write(std_out,*) 'lezes sem?' 859 do ii=0,aim_dtset%ngrid(2)-1 860 do jj=0,aim_dtset%ngrid(1)-1 861 nn=nn+1 862 xx=jj*stp(1)+cos(alfa)*ii*stp(2) 863 yy=sin(alfa)*ii*stp(2) 864 write(untl,'(3E16.8)') xx, yy, lfld(nn) 865 end do 866 write(untl,*) ' ' 867 end do 868 end if 869 end select 870 ABI_DEALLOCATE(uu) 871 ABI_DEALLOCATE(stp) 872 ABI_DEALLOCATE(nr) 873 if(aim_dtset%denout>0) then 874 ABI_DEALLOCATE(dfld) 875 end if 876 if(aim_dtset%lapout>0) then 877 ABI_DEALLOCATE(lfld) 878 end if 879 880 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
PARENTS
aim
CHILDREN
consist,inread
SOURCE
235 subroutine adini(aim_dtset,inpstr,lenstr) 236 237 238 !This section has been created automatically by the script Abilint (TD). 239 !Do not modify the following lines by hand. 240 #undef ABI_FUNC 241 #define ABI_FUNC 'adini' 242 !End of the abilint section 243 244 implicit none 245 246 !Arguments ------------------------------------ 247 !scalars 248 integer,intent(in) :: lenstr 249 character(len=*),intent(in) :: inpstr 250 !no_abirules 251 type(aim_dataset_type), intent(inout) :: aim_dtset !vz_i 252 253 !Local variables ------------------------------ 254 !scalars 255 integer :: errcod,ii,inxh,ipos,jj,lenc,ll,outi,tstngr=0,tstvpt=0 !vz_z 256 real(dp) :: outr 257 logical :: nbtst,try 258 character(len=20) :: cmot 259 260 ! ********************************************************************* 261 262 if (iachar(inpstr(1:1)) < 32) then 263 ipos=2 264 else 265 ipos=1 266 end if 267 268 write(std_out,*) 'ECHO of the INPUT' 269 write(std_out,*) '************************' 270 write(untout,*) 'ECHO of the INPUT' 271 write(untout,*) '************************' 272 273 mread: do ii=1,lenstr 274 try=.false. 275 nbtst=.true. 276 inxh=index(inpstr(ipos:lenstr),' ') 277 if ((ipos >= lenstr)) exit 278 if ((inxh==2).or.(inxh==1)) then 279 ipos=ipos+inxh 280 cycle 281 end if 282 lenc=inxh-1 283 cmot(1:lenc)=inpstr(ipos:ipos+inxh-2) 284 ipos=ipos+inxh 285 ! write(std_out,*) cmot(1:lenc), lenc 286 287 select case (cmot(1:lenc)) 288 289 ! DRIVER SPECIFICATIONS 290 291 case ('SURF') 292 inxh=index(inpstr(ipos:lenstr),' ') 293 if ((inxh /= 2).and.(inpstr(ipos:ipos)/='-')) then 294 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 295 MSG_ERROR("Aborting now") 296 end if 297 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 298 aim_dtset%isurf=outi 299 write(std_out,*) cmot(1:lenc),' ', aim_dtset%isurf 300 write(untout,*) cmot(1:lenc),' ', aim_dtset%isurf 301 ipos=ipos+inxh 302 303 case ('CRIT') 304 inxh=index(inpstr(ipos:lenstr),' ') 305 if ((inxh /= 2).and.(inpstr(ipos:ipos)/='-')) then 306 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 307 MSG_ERROR("Aborting now") 308 end if 309 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 310 aim_dtset%crit=outi 311 write(std_out,*) cmot(1:lenc),' ', aim_dtset%crit 312 write(untout,*) cmot(1:lenc),' ', aim_dtset%crit 313 ipos=ipos+inxh 314 315 case ('RSURF') 316 inxh=index(inpstr(ipos:lenstr),' ') 317 if (inxh /= 2) then 318 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 319 MSG_ERROR("Aborting now") 320 end if 321 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 322 aim_dtset%irsur=outi 323 write(std_out,*) cmot(1:lenc),' ', aim_dtset%irsur 324 write(untout,*) cmot(1:lenc),' ', aim_dtset%irsur 325 ipos=ipos+inxh 326 327 case ('FOLLOW') 328 inxh=index(inpstr(ipos:lenstr),' ') 329 if (inxh /= 2) then 330 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 331 MSG_ERROR("Aborting now") 332 end if 333 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 334 aim_dtset%foll=outi 335 write(std_out,*) cmot(1:lenc),' ', aim_dtset%foll 336 write(untout,*) cmot(1:lenc),' ', aim_dtset%foll 337 ipos=ipos+inxh 338 339 case ('IRHO') 340 inxh=index(inpstr(ipos:lenstr),' ') 341 if (inxh /= 2) then 342 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 343 MSG_ERROR("Aborting now") 344 end if 345 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 346 aim_dtset%irho=outi 347 write(std_out,*) cmot(1:lenc),' ', aim_dtset%irho 348 write(untout,*) cmot(1:lenc),' ', aim_dtset%irho 349 ipos=ipos+inxh 350 351 case ('PLDEN') 352 inxh=index(inpstr(ipos:lenstr),' ') 353 if (inxh /= 2) then 354 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 355 MSG_ERROR("Aborting now") 356 end if 357 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 358 aim_dtset%plden=outi 359 write(std_out,*) cmot(1:lenc),' ', aim_dtset%plden 360 write(untout,*) cmot(1:lenc),' ', aim_dtset%plden 361 ipos=ipos+inxh 362 363 364 case ('IVOL') 365 inxh=index(inpstr(ipos:lenstr),' ') 366 if (inxh /= 2) then 367 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 368 MSG_ERROR("Aborting now") 369 end if 370 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 371 aim_dtset%ivol=outi 372 write(std_out,*) cmot(1:lenc),' ', aim_dtset%ivol 373 write(untout,*) cmot(1:lenc),' ', aim_dtset%ivol 374 ipos=ipos+inxh 375 376 case ('DENOUT') 377 inxh=index(inpstr(ipos:lenstr),' ') 378 if (inxh /= 2) then 379 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 380 MSG_ERROR("Aborting now") 381 end if 382 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 383 aim_dtset%denout=outi 384 if ((aim_dtset%denout < -1).or.(aim_dtset%denout>3)) then 385 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 386 MSG_ERROR("Aborting now") 387 end if 388 write(std_out,*) cmot(1:lenc),' ', aim_dtset%denout 389 write(untout,*) cmot(1:lenc),' ', aim_dtset%denout 390 ipos=ipos+inxh 391 392 case ('LAPOUT') 393 inxh=index(inpstr(ipos:lenstr),' ') 394 if (inxh /= 2) then 395 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 396 MSG_ERROR("Aborting now") 397 end if 398 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 399 aim_dtset%lapout=outi 400 if ((aim_dtset%lapout < -1).or.(aim_dtset%lapout>3)) then 401 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 402 MSG_ERROR("Aborting now") 403 end if 404 write(std_out,*) cmot(1:lenc),' ', aim_dtset%lapout 405 write(untout,*) cmot(1:lenc),' ', aim_dtset%lapout 406 ipos=ipos+inxh 407 408 case ('DLTYP') 409 inxh=index(inpstr(ipos:lenstr),' ') 410 if (inxh /= 2) then 411 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 412 MSG_ERROR("Aborting now") 413 end if 414 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 415 aim_dtset%dltyp=outi 416 write(std_out,*) cmot(1:lenc),' ', aim_dtset%dltyp 417 write(untout,*) cmot(1:lenc),' ', aim_dtset%dltyp 418 ipos=ipos+inxh 419 420 case ('GPSURF') 421 inxh=index(inpstr(ipos:lenstr),' ') 422 if (inxh /= 2) then 423 write(std_out,*) 'ERROR in specif. of ', cmot(1:lenc) 424 MSG_ERROR("Aborting now") 425 end if 426 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 427 aim_dtset%gpsurf=outi 428 write(std_out,*) cmot(1:lenc),' ', aim_dtset%gpsurf 429 write(untout,*) cmot(1:lenc),' ', aim_dtset%gpsurf 430 ipos=ipos+inxh 431 432 433 ! END OF THE DRIVER SPECIFICATIONS 434 435 case ('ATOM') 436 inxh=index(inpstr(ipos:lenstr),' ') 437 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 438 aim_dtset%batom=outi 439 write(std_out,*) cmot(1:lenc),' ', aim_dtset%batom 440 write(untout,*) cmot(1:lenc),' ', aim_dtset%batom 441 ipos=ipos+inxh 442 443 case ('NSA') 444 inxh=index(inpstr(ipos:lenstr),' ') 445 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 446 aim_dtset%nsa=outi 447 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsa 448 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsa 449 ipos=ipos+inxh 450 451 case ('NSB') 452 inxh=index(inpstr(ipos:lenstr),' ') 453 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 454 aim_dtset%nsb=outi 455 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsb 456 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsb 457 ipos=ipos+inxh 458 459 case ('NSC') 460 inxh=index(inpstr(ipos:lenstr),' ') 461 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 462 aim_dtset%nsc=outi 463 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nsc 464 write(untout,*) cmot(1:lenc),' ', aim_dtset%nsc 465 ipos=ipos+inxh 466 467 case ('INPT') 468 inxh=index(inpstr(ipos:lenstr),' ') 469 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 470 aim_dtset%npt=outi 471 write(std_out,*) cmot(1:lenc),' ', aim_dtset%npt 472 write(untout,*) cmot(1:lenc),' ', aim_dtset%npt 473 ipos=ipos+inxh 474 475 case ('NTHETA') 476 inxh=index(inpstr(ipos:lenstr),' ') 477 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 478 aim_dtset%nth=outi 479 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nth 480 write(untout,*) cmot(1:lenc),' ', aim_dtset%nth 481 ipos=ipos+inxh 482 483 case ('NPHI') 484 inxh=index(inpstr(ipos:lenstr),' ') 485 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 486 aim_dtset%nph=outi 487 write(std_out,*) cmot(1:lenc),' ', aim_dtset%nph 488 write(untout,*) cmot(1:lenc),' ', aim_dtset%nph 489 ipos=ipos+inxh 490 491 case ('THETAMIN') 492 inxh=index(inpstr(ipos:lenstr),' ') 493 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 494 aim_dtset%themin=outr 495 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themin 496 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themin 497 ipos=ipos+inxh 498 499 case ('THETAMAX') 500 inxh=index(inpstr(ipos:lenstr),' ') 501 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 502 aim_dtset%themax=outr 503 write(std_out, '(1x,a,a,es17.10)' ) cmot(1:lenc),' ', aim_dtset%themax 504 write(untout,'(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%themax 505 ipos=ipos+inxh 506 507 case ('PHIMIN') 508 inxh=index(inpstr(ipos:lenstr),' ') 509 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 510 aim_dtset%phimin=outr 511 write(std_out, '(1x,a,a,es17.10)' ) cmot(1:lenc),' ', aim_dtset%phimin 512 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimin 513 ipos=ipos+inxh 514 515 case ('PHIMAX') 516 inxh=index(inpstr(ipos:lenstr),' ') 517 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 518 aim_dtset%phimax=outr 519 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimax 520 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%phimax 521 ipos=ipos+inxh 522 523 case ('ATRAD') 524 inxh=index(inpstr(ipos:lenstr),' ') 525 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 526 aim_dtset%atrad=outr 527 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%atrad 528 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%atrad 529 ipos=ipos+inxh 530 531 case ('RADSTP') 532 inxh=index(inpstr(ipos:lenstr),' ') 533 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 534 aim_dtset%dr0=outr 535 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dr0 536 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dr0 537 ipos=ipos+inxh 538 539 case ('FOLSTP') 540 inxh=index(inpstr(ipos:lenstr),' ') 541 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 542 aim_dtset%folstp=outr 543 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%folstp 544 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%folstp 545 ipos=ipos+inxh 546 547 case ('RATMIN') 548 inxh=index(inpstr(ipos:lenstr),' ') 549 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 550 aim_dtset%rmin=outr 551 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%rmin 552 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%rmin 553 ipos=ipos+inxh 554 555 case ('COFF1') 556 inxh=index(inpstr(ipos:lenstr),' ') 557 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 558 aim_dtset%coff1=outr 559 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff1 560 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff1 561 ipos=ipos+inxh 562 563 case ('COFF2') 564 inxh=index(inpstr(ipos:lenstr),' ') 565 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 566 aim_dtset%coff2=outr 567 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff2 568 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%coff2 569 ipos=ipos+inxh 570 571 case ('DPCLIM') 572 inxh=index(inpstr(ipos:lenstr),' ') 573 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 574 aim_dtset%dpclim=outr 575 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dpclim 576 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%dpclim 577 ipos=ipos+inxh 578 579 case ('LGRAD') 580 inxh=index(inpstr(ipos:lenstr),' ') 581 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 582 aim_dtset%lgrad=outr 583 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad 584 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad 585 ipos=ipos+inxh 586 587 case ('LGRAD2') 588 inxh=index(inpstr(ipos:lenstr),' ') 589 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 590 aim_dtset%lgrad2=outr 591 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad2 592 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lgrad2 593 ipos=ipos+inxh 594 595 case ('LSTEP') 596 inxh=index(inpstr(ipos:lenstr),' ') 597 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 598 aim_dtset%lstep=outr 599 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep 600 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep 601 ipos=ipos+inxh 602 603 case ('LSTEP2') 604 inxh=index(inpstr(ipos:lenstr),' ') 605 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 606 aim_dtset%lstep2=outr 607 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep2 608 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%lstep2 609 ipos=ipos+inxh 610 611 case ('RSURDIR') 612 inxh=index(inpstr(ipos:lenstr),' ') 613 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 614 aim_dtset%th0=outr 615 ipos=ipos+inxh 616 inxh=index(inpstr(ipos:lenstr),' ') 617 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 618 aim_dtset%phi0=outr 619 ipos=ipos+inxh 620 write(std_out, '(1x,a,a,2es17.10)') cmot(1:lenc),' ', aim_dtset%th0, aim_dtset%phi0 621 write(untout, '(1x,a,a,2es17.10)') cmot(1:lenc),' ', aim_dtset%th0, aim_dtset%phi0 622 623 case ('FOLDEP') 624 do jj=1,3 625 inxh=index(inpstr(ipos:lenstr),' ') 626 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 627 aim_dtset%foldep(jj)=outr 628 ipos=ipos+inxh 629 end do 630 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%foldep 631 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%foldep 632 633 case ('SCAL') 634 do jj=1,3 635 inxh=index(inpstr(ipos:lenstr),' ') 636 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 637 aim_dtset%scal(jj)=outr 638 ipos=ipos+inxh 639 end do 640 write(std_out,*) cmot(1:lenc),' ', aim_dtset%scal 641 write(untout,*) cmot(1:lenc),' ', aim_dtset%scal 642 643 case ('NGRID') 644 try=.true. 645 do jj=1,3 646 inxh=index(inpstr(ipos:lenstr),' ') 647 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"INT",outi,outr,errcod) 648 aim_dtset%ngrid(jj)=outi 649 if (.not.nbtst) then 650 tstngr=jj-1 651 cycle mread 652 end if 653 if (inxh==0) then 654 tstvpt=jj 655 exit mread 656 end if 657 ipos=ipos+inxh 658 if (ipos==lenstr-1) then 659 tstngr=jj 660 exit mread 661 end if 662 end do 663 ! Why no echo ?? XG 030218 664 tstngr=3 665 666 case ('VPTS') 667 do jj=1,4 668 do ll=1,3 669 try=.true. 670 inxh=index(inpstr(ipos:lenstr),' ') 671 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 672 aim_dtset%vpts(ll,jj)=outr 673 if (.not.nbtst) then 674 tstvpt=jj-1 675 cycle mread 676 end if 677 ipos=ipos+inxh 678 if (ipos>=lenstr) then 679 tstvpt=jj 680 exit mread 681 end if 682 end do 683 end do 684 ! Why no echo ?? XG 030218 685 tstvpt=4 686 687 case ('MAXATD') 688 inxh=index(inpstr(ipos:lenstr),' ') 689 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 690 aim_dtset%maxatd=outr 691 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxatd 692 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxatd 693 ipos=ipos+inxh 694 695 case ('MAXCPD') 696 inxh=index(inpstr(ipos:lenstr),' ') 697 call inread(inpstr(ipos:ipos+inxh-2),inxh-1,"DPR",outi,outr,errcod) 698 aim_dtset%maxcpd=outr 699 write(std_out, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxcpd 700 write(untout, '(1x,a,a,es17.10)') cmot(1:lenc),' ', aim_dtset%maxcpd 701 ipos=ipos+inxh 702 703 case default 704 write(std_out,*) 'ERROR Bad key word ! ',cmot(1:lenc) 705 end select 706 end do mread 707 708 write(std_out,*) '************************' 709 710 call consist(aim_dtset,tstngr,tstvpt) 711 712 end subroutine adini
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
PARENTS
cpdrv,drvaim,rsurf
CHILDREN
critic,onestep,timein,vgh_rho
SOURCE
915 subroutine aim_follow(aim_dtset,vv,npmax,srch,iatinit,iposinit,iat,ipos,nstep) 916 917 918 !This section has been created automatically by the script Abilint (TD). 919 !Do not modify the following lines by hand. 920 #undef ABI_FUNC 921 #define ABI_FUNC 'aim_follow' 922 !End of the abilint section 923 924 implicit none 925 926 !Arguments ------------------------------------ 927 !scalars 928 integer,intent(in) :: iatinit,iposinit,npmax 929 integer,intent(out) :: iat,ipos,nstep 930 logical,intent(inout) :: srch 931 type(aim_dataset_type),intent(in) :: aim_dtset 932 !arrays 933 real(dp),intent(inout) :: vv(3) 934 935 !Local variables ------------------------------ 936 !scalars 937 integer :: i1,i2,i3,ii,iph,ires,ith,jj,kk,nit,np,nph,nsi,nth 938 real(dp) :: deltar,dg,dist,dph,dth,fac2,facf,h0old,hh,hold,rho,rr,rsmed 939 real(dp) :: t1,t2,t3,vcth,vph,vth,wall,xy,xyz 940 logical :: fin,ldebold,srchold,stemp,stemp2 941 character(len=50) :: formpc 942 character(len=500) :: msg 943 !arrays 944 real(dp) :: ev(3),grho(3),hrho(3,3),pom(3),vold(3),vt(3),vt1(3) 945 real(dp) :: zz(3,3) 946 947 !************************************************************************ 948 formpc='(":CP",2I5,3F12.8,3E12.4,I4,2E12.4)' 949 950 951 fin=.false. 952 953 srchold=srch 954 ldebold=ldeb 955 h0old=h0 956 957 nth=aim_dtset%nth 958 nph=aim_dtset%nph 959 960 if (slc==0) then 961 rminl(:)=aim_dtset%rmin 962 end if 963 964 if (deb) then 965 ldeb=.true. 966 end if 967 968 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 969 970 !Initial tests 971 972 if (iat/=0) then 973 if (rr<rminl(iat)) then 974 fin=.true. 975 write(std_out,*) 'rr < rmin iat=',iat,' ipos=',ipos 976 elseif (rho<aim_rhomin) then 977 fin=.true. 978 write(std_out,*) 'CHARGE LT rhomin ',rho,' < ',aim_rhomin 979 if (rho<zero) then 980 MSG_ERROR('RHO < 0 !!!') 981 end if 982 end if 983 end if 984 985 facf=aim_fac0 986 hh=aim_hmax 987 988 call timein(t1,wall) 989 nstep=0 990 nsi=0 991 992 !the principal cycle 993 994 madw : do while(.not.fin) 995 hold=hh 996 997 dg=vnorm(grho,0) 998 if (ldeb.or.deb) write(std_out,*) 'dg= ',dg 999 1000 ! the time test 1001 1002 call timein(t3,wall) 1003 t2=t3-t1 1004 if (t2>300.0) then 1005 write(std_out,*) 'TIME EXCEEDED 5 min IN FOLLOW' 1006 write(std_out,*) 'h0 =',h0,' h =',hh,' h0old =',h0old,' dg =',dg 1007 write(std_out,*) 'facf =',facf 1008 msg = 'TIME EXCEEDED 5 min IN FOLLOW' 1009 MSG_ERROR(msg) 1010 end if 1011 1012 if (dg<aim_dgmin) then 1013 write(std_out,*) 'gradient < dgmin ',dg,' < ',aim_dgmin 1014 fin=.true. 1015 iat=0 1016 ipos=0 1017 ! testing for the CP 1018 if (npc>0) then 1019 call critic(aim_dtset,vv,ev,zz,aim_dmaxcrit,ires,0) 1020 if (ires==0) then 1021 do jj=1,npc 1022 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 1023 dist=vnorm(pom,0) 1024 if (dist<aim_tiny) cycle madw 1025 end do 1026 write(std_out,*) 'C.P. found !!' 1027 npc=npc+1 1028 do jj=1,3 1029 pc(jj,npc)=vv(jj) 1030 evpc(jj,npc)=ev(jj) 1031 do kk=1,3 1032 zpc(kk,jj,npc)=zz(kk,jj) 1033 end do 1034 end do 1035 i1=ev(1)/abs(ev(1)) 1036 i2=ev(2)/abs(ev(2)) 1037 i3=ev(3)/abs(ev(3)) 1038 icpc(npc)=i1+i2+i3 1039 if (icpc(npc)==-3) then ! pseudoatom handling 1040 npcm3=npcm3+1 1041 write(std_out,*) 'Pseudo-atom found !!' 1042 end if 1043 1044 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 1045 write(22,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 1046 & ev(1)+ev(2)+ev(3),rho 1047 write(std_out,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 1048 & ev(1)+ev(2)+ev(3),rho 1049 else 1050 write(std_out,*) 'C.P. not found !!' 1051 end if 1052 end if 1053 1054 cycle madw 1055 end if 1056 1057 hh=h0/dg 1058 if (ldeb.or.deb) write(std_out,*) 'h= ',hh,' h0= ',h0,' dg= ',dg 1059 if (hh>aim_hmax) hh=aim_hmax 1060 ! step modifications 1061 1062 hh=hh*facf 1063 if (hh>(hold*aim_hmult)) then 1064 hh=hold*aim_hmult 1065 end if 1066 1067 do ii=1,3 1068 vold(ii)=vv(ii) 1069 end do 1070 1071 nit=0 1072 hold=hh 1073 1074 ! one step following the gradient line 1075 ! 1076 call onestep(vv,rho,grho,hh,np,npmax,deltar) 1077 do while (((np>npmax).or.(deltar>aim_stmax)).and.(deltar>aim_dmin)) 1078 nit=nit+1 1079 if (nit>5) then 1080 if (deltar>aim_stmax) then 1081 write(std_out,*) 'nit > 5 and deltar > stmax nit=',nit 1082 else 1083 write(std_out,*) 'nit > 5 and np > npmax nit=',nit 1084 end if 1085 end if 1086 do ii=1,3 1087 vv(ii)=vold(ii) 1088 end do 1089 hh=hh*0.3 1090 call onestep(vv,rho,grho,hh,np,npmax,deltar) 1091 end do 1092 1093 1094 nstep=nstep+1 1095 if (ldeb.or.deb) write(std_out,*) 'h= ',hh 1096 1097 fac2=hh/hold 1098 if (fac2>=1._dp) then 1099 facf=facf*1.2 1100 else 1101 if (fac2>=aim_facmin) then 1102 facf=fac2 1103 else 1104 facf=aim_facmin 1105 end if 1106 end if 1107 1108 if (deb.or.ldeb) then 1109 write(std_out,*) ':POS ',vv 1110 write(std_out,*) ':RBPOS ',vt1 1111 write(std_out,*) ':GRAD ',grho 1112 end if 1113 1114 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 1115 dg=vnorm(grho,0) 1116 pom(:)=vv(:)-xatm(:,iatinit)-atp(:,iposinit) 1117 1118 if (iat /= 0) then 1119 fin=.true. 1120 write(std_out,*) 'r < rmin iat=',iat,' ipos=',ipos 1121 cycle madw 1122 end if 1123 1124 if (rho<aim_rhomin) then 1125 fin=.true. 1126 write(std_out,*) 'charge < rhomin ',rho,' < ',aim_rhomin 1127 if (rho<zero) then 1128 MSG_ERROR('RHO < 0 !!!') 1129 end if 1130 iat=0 1131 ipos=0 1132 cycle madw 1133 end if 1134 1135 if (npcm3>0) then 1136 do jj=1,npc 1137 if (icpc(jj)==(-3)) then 1138 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 1139 dist=vnorm(pom,0) 1140 if (dist<(aim_dtset%rmin**2*0.1)) then 1141 iat=0 1142 ipos=0 1143 fin=.true. 1144 write(std_out,*) 'We are inside a pseudo-atom' 1145 cycle madw 1146 end if 1147 end if 1148 end do 1149 end if 1150 1151 nsi=nsi+1 1152 1153 ! surface checking 1154 1155 if (srch.and.(nsi>=nsimax)) then 1156 nsi=0 1157 ith=0 1158 iph=0 1159 do ii=1,3 1160 vt(ii)=vv(ii)-xatm(ii,iatinit) 1161 end do 1162 xy=vt(1)*vt(1)+vt(2)*vt(2) 1163 xyz=xy+vt(3)*vt(3) 1164 xyz=sqrt(xyz) 1165 if (xy<aim_snull) then 1166 vcth=1._dp 1167 if (vt(3)<0._dp) vcth=-vcth 1168 vph=0._dp 1169 else 1170 vcth=vt(3)/xyz 1171 vph=atan2(vt(2),vt(1)) 1172 end if 1173 vth=acos(vcth) 1174 if (vth<th(1)) then 1175 ith=0 1176 else 1177 if (vth>th(nth)) then 1178 ith=nth 1179 else 1180 do ii=2,nth 1181 if (vth<th(ii)) then 1182 ith=ii-1 1183 exit 1184 end if 1185 end do 1186 end if 1187 end if 1188 1189 if (vph<ph(1)) then 1190 iph=0 1191 else 1192 if (vph>ph(nph)) then 1193 iph=nph 1194 else 1195 do ii=2,nph 1196 if (vph<ph(ii)) then 1197 iph=ii-1 1198 exit 1199 end if 1200 end do 1201 end if 1202 end if 1203 1204 stemp=(iph>0).and.(iph<nph) 1205 stemp=stemp.and.(ith>0).and.(ith<nth) 1206 1207 if (stemp) then 1208 stemp2=rs(ith,iph)>0._dp 1209 stemp2=stemp2.and.(rs(ith+1,iph)>0._dp) 1210 stemp2=stemp2.and.(rs(ith+1,iph+1)>0._dp) 1211 stemp2=stemp2.and.(rs(ith,iph+1)>0._dp) 1212 if (stemp2) then 1213 dth=th(ith+1)-th(ith) 1214 dph=ph(iph+1)-ph(iph) 1215 rsmed=rs(ith,iph)*(th(ith+1)-vth)/dth*(ph(iph+1)-vph)/dph 1216 rsmed=rsmed+rs(ith+1,iph)*(vth-th(ith))/dth*(ph(iph+1)-vph)/dph 1217 rsmed=rsmed+rs(ith+1,iph+1)*(vth-th(ith))/dth*(vph-ph(iph))/dph 1218 rsmed=rsmed+rs(ith,iph+1)*(th(ith+1)-vth)/dth*(vph-ph(iph))/dph 1219 if (rsmed>xyz) then 1220 write(std_out,*) 'We are inside the surface' 1221 iat=iatinit 1222 ipos=iposinit 1223 else 1224 write(std_out,*) 'We are outside the surface' 1225 iat=0 1226 ipos=0 1227 end if 1228 fin=.true. 1229 cycle madw 1230 end if 1231 end if 1232 end if 1233 1234 end do madw 1235 1236 1237 srch=srchold 1238 ldeb=ldebold 1239 h0=h0old 1240 1241 1242 end subroutine aim_follow
m_bader/bcp_type [ Types ]
NAME
bcp_type
FUNCTION
a "bonding critical point" for aim
SOURCE
129 type, private :: bcp_type 130 131 ! Integer 132 integer :: iat ! number of the bonding atom inside a primitive cell 133 integer :: ipos ! number of the primitive cell of the bonding atom 134 135 ! Real 136 real(dp) :: chg ! charge at the critical point 137 real(dp) :: diff(3) ! three distances : AT-CP,BAT-CP,AT-BAT 138 real(dp) :: ev(3) ! eigenvalues of the Hessian 139 real(dp) :: pom(3) ! position of the bonding atom 140 real(dp) :: rr(3) ! position of the bcp 141 real(dp) :: vec(3,3) ! eigenvectors of the Hessian 142 real(dp) :: vv(3) ! position of the bcp relative to the central atom 143 144 end type bcp_type
m_bader/bschg1 [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
bschg1
FUNCTION
bschg1: Vector transformation of coordinates
PARENTS
critics,initaim,integrho,vgh_rho
CHILDREN
mprod
SOURCE
6123 subroutine bschg1(vv,dir) 6124 6125 6126 !This section has been created automatically by the script Abilint (TD). 6127 !Do not modify the following lines by hand. 6128 #undef ABI_FUNC 6129 #define ABI_FUNC 'bschg1' 6130 !End of the abilint section 6131 6132 implicit none 6133 6134 !Arguments ------------------------------------ 6135 !scalars 6136 integer,intent(in) :: dir 6137 !arrays 6138 real(dp),intent(inout) :: vv(3) 6139 6140 !Local variables ------------------------------ 6141 !scalars 6142 integer :: ii 6143 !arrays 6144 real(dp) :: vt(3) 6145 6146 ! ********************************************************************* 6147 6148 if (dir==1) then 6149 do ii=1,3 6150 vt(ii)=rprimd(ii,1)*vv(1)+rprimd(ii,2)*vv(2)+rprimd(ii,3)*vv(3) 6151 end do 6152 elseif (dir==-1) then 6153 do ii=1,3 6154 vt(ii)=ivrprim(ii,1)*vv(1)+ivrprim(ii,2)*vv(2)+ivrprim(ii,3)*vv(3) 6155 end do 6156 elseif (dir==2) then 6157 do ii=1,3 6158 vt(ii)=trivrp(ii,1)*vv(1)+trivrp(ii,2)*vv(2)+trivrp(ii,3)*vv(3) 6159 end do 6160 else 6161 MSG_ERROR('Transformation of coordinates') 6162 end if 6163 vv(:)=vt(:) 6164 6165 end subroutine bschg1
m_bader/bschg2 [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
bschg2
FUNCTION
bschg2: Matrix transformation of coordinates
PARENTS
vgh_rho
CHILDREN
mprod
SOURCE
6183 subroutine bschg2(aa,dir) 6184 6185 6186 !This section has been created automatically by the script Abilint (TD). 6187 !Do not modify the following lines by hand. 6188 #undef ABI_FUNC 6189 #define ABI_FUNC 'bschg2' 6190 !End of the abilint section 6191 6192 implicit none 6193 6194 !Arguments ------------------------------------ 6195 !scalars 6196 integer,intent(in) :: dir 6197 !arrays 6198 real(dp),intent(inout) :: aa(3,3) 6199 6200 !Local variables ------------------------------ 6201 !arrays 6202 real(dp) :: bb(3,3) 6203 6204 ! ********************************************************************* 6205 6206 if (dir==1) then 6207 call mprod(aa,ivrprim,bb) 6208 call mprod(rprimd,bb,aa) 6209 elseif (dir==2) then 6210 call mprod(aa,ivrprim,bb) 6211 call mprod(trivrp,bb,aa) 6212 elseif (dir==-1) then 6213 call mprod(aa,rprimd,bb) 6214 call mprod(ivrprim,bb,aa) 6215 else 6216 MSG_ERROR("transformation of coordinates") 6217 end if 6218 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)
PARENTS
adini
CHILDREN
SOURCE
1270 subroutine consist(aim_dtset,tstngr,tstvpt) 1271 1272 1273 !This section has been created automatically by the script Abilint (TD). 1274 !Do not modify the following lines by hand. 1275 #undef ABI_FUNC 1276 #define ABI_FUNC 'consist' 1277 !End of the abilint section 1278 1279 implicit none 1280 1281 !Arguments ------------------------------------ 1282 !scalars 1283 integer,intent(in) :: tstngr,tstvpt 1284 type(aim_dataset_type),intent(in) :: aim_dtset 1285 1286 !Local variables ------------------------------ 1287 1288 ! ********************************************************************* 1289 1290 !write(std_out,*) tstngr, tstvpt 1291 1292 if (((aim_dtset%denout/=0).or.(aim_dtset%lapout/=0)).and.((tstngr < 1).or.(tstvpt < 2))) then 1293 MSG_ERROR('in input1 - I cannot do the output !') 1294 end if 1295 if ((aim_dtset%denout > 0).and.(aim_dtset%lapout>0)) then 1296 if (aim_dtset%denout/=aim_dtset%lapout) then 1297 write(std_out,*) 'ERROR in input - when both denout and lapout are positive non-zero,' 1298 write(std_out,*) 'they must be equal.' 1299 MSG_ERROR("Aborting now") 1300 end if 1301 if ((tstvpt < aim_dtset%denout+1).or.(tstngr < aim_dtset%denout)) then 1302 write(std_out,*) 'ERROR in input2 - I cannot do the output !' 1303 MSG_ERROR("Aborting now") 1304 end if 1305 elseif (aim_dtset%denout > 0) then 1306 if ((tstvpt < aim_dtset%denout+1).or.(tstngr < aim_dtset%denout)) then 1307 write(std_out,*) 'ERROR in input - I cannot do the output !' 1308 MSG_ERROR("Aborting now") 1309 end if 1310 elseif (aim_dtset%lapout > 0) then 1311 if ((tstvpt < aim_dtset%lapout+1).or.(tstngr < aim_dtset%lapout)) then 1312 write(std_out,*) 'ERROR in input - I cannot do the output !' 1313 MSG_ERROR("Aborting now") 1314 end if 1315 end if 1316 1317 if ((aim_dtset%isurf==1).and.(aim_dtset%crit==0)) then 1318 write(std_out,*) 'ERROR in input - must have crit/=0 for isurf==1' 1319 MSG_ERROR("Aborting now") 1320 end if 1321 1322 if (((aim_dtset%ivol/=0).or.(aim_dtset%irho/=0)).and.(aim_dtset%isurf==0)) then 1323 MSG_ERROR('in input - I cannot integrate without surface !') 1324 end if 1325 1326 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 ...
PARENTS
drvaim
CHILDREN
aim_follow,critic,sort_dp,timein,vgh_rho,xmpi_sum
SOURCE
1362 subroutine cpdrv(aim_dtset) 1363 1364 1365 !This section has been created automatically by the script Abilint (TD). 1366 !Do not modify the following lines by hand. 1367 #undef ABI_FUNC 1368 #define ABI_FUNC 'cpdrv' 1369 !End of the abilint section 1370 1371 implicit none 1372 1373 !Arguments ------------------------------------ 1374 !scalars 1375 type(aim_dataset_type),intent(in) :: aim_dtset 1376 1377 !Local variables ------------------------------ 1378 !scalars 1379 integer :: iat,iatinit,ii,inxat,inxcell,ipair,ipos,iposinit,ires,jj,kk,nb,nb_now 1380 integer :: nn,nstep,nvs,me,nproc,ierr 1381 real(dp) :: candidate,chg,diff1,diff2,diff3,dist,prj,rtdiff,ss,tt0,wall 1382 logical :: srch=.false. 1383 !arrays 1384 integer :: ibat(nnpos*natom),inatm(nnpos*natom),incell(nnpos*natom) 1385 integer :: ipibat(nnpos*natom) 1386 integer,allocatable :: indexcp(:),nr(:) 1387 real(dp) :: bmin(natom),dif(3),dists(nnpos*natom),ev(3),evec(3,3),grho(3) 1388 real(dp) :: hrho(3,3),pom(3),rr(3),uu(3),vv(3),xorig(3) 1389 real(dp),allocatable :: buffer(:,:),sortguide(:) 1390 !no_abirules 1391 !Warning : bcp_type should be transformed to cp_type 1392 type(bcp_type),allocatable :: bcp(:),ccp(:),cp_tmp(:),rcp(:) 1393 1394 !************************************************************************ 1395 1396 me=xmpi_comm_rank(xmpi_world) 1397 nproc=xmpi_comm_size(xmpi_world) 1398 1399 !Consider the critical points starting from atom #batom 1400 inxat=aim_dtset%batom 1401 slc=-1 1402 rminl(:)=aim_dtset%rmin 1403 bmin(:)=0._dp 1404 ttcp=0._dp 1405 1406 write(std_out,*) 1407 write(std_out,*) "CRITICAL POINTS ANALYSIS" 1408 write(std_out,*) "========================" 1409 write(std_out,*) 1410 1411 write(untout,*) 1412 write(untout,*) "CRITICAL POINTS ANALYSIS" 1413 write(untout,*) "========================" 1414 write(untout,*) 1415 1416 1417 xorig(:)=xatm(:,inxat) 1418 1419 call timein(tt0,wall) 1420 1421 !Searching the neighbouring atoms 1422 1423 if (aim_dtset%crit > 0) then 1424 nvs=0 1425 do ii=1,nnpos 1426 do jj=1,natom 1427 dist=0._dp 1428 dif(:)=xatm(:,inxat)-xatm(:,jj)-atp(:,ii) 1429 dif(:)=dif(:)/aim_dtset%scal(:) 1430 dist=vnorm(dif,0) 1431 if (dist < tol6 ) then 1432 inxcell=ii 1433 elseif (dist < maxatdst) then 1434 nvs=nvs+1 1435 dists(nvs)=dist 1436 inatm(nvs)=jj 1437 incell(nvs)=ii 1438 end if 1439 end do 1440 end do 1441 1442 write(std_out,*) "ATOM:" 1443 write(std_out,*) 'inxat :', inxat, 'inxcell :', inxcell 1444 write(std_out, '(3es16.6)' ) (xorig(ii),ii=1,3) 1445 write(std_out,*) 1446 1447 write(untout,*) "ATOM:" 1448 write(untout,*) 'inxat :', inxat, 'inxcell :', inxcell 1449 write(untout, '(3es16.6)') (xorig(ii),ii=1,3) 1450 write(untout,*) 1451 1452 ABI_ALLOCATE(nr,(nvs)) 1453 do ii=1,nvs 1454 nr(ii)=ii 1455 end do 1456 1457 ! Ordering of the nearest neighbouring atoms 1458 call sort_dp(nvs,dists,nr,tol14) 1459 1460 nb=0 1461 write(std_out,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 1462 write(untout,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 1463 do ii=1,nvs 1464 nn=nr(ii) 1465 if (dists(ii) < maxatdst) then 1466 nb=nb+1 1467 ibat(nb)=inatm(nn) 1468 ipibat(nb)=incell(nn) 1469 write(std_out,*) ':NEIG ',inatm(nn),incell(nn),dists(ii) 1470 write(untout,'(" ",2I6,F16.8)')inatm(nn),incell(nn),dists(ii) 1471 else 1472 exit 1473 end if 1474 end do 1475 1476 ! SEARCHING BCP 1477 ABI_DATATYPE_ALLOCATE(bcp,(nb)) 1478 nbcp=0 1479 iatinit=inxat 1480 iposinit=inxcell 1481 bcp(:)%iat=0 1482 bcp(:)%ipos=0 1483 1484 write(std_out,*) 1485 write(std_out,*) "BONDING CRITICAL POINTS (BCP)" 1486 write(std_out,*) "=============================" 1487 write(std_out,*) 1488 1489 write(untout,*) 1490 write(untout,*) "BONDING CRITICAL POINTS (BCP)" 1491 write(untout,*) "=============================" 1492 write(untout,*) 1493 1494 srbcp: do ii=1,nb 1495 1496 ! Start the search for BCP from the midistance between the atom 1497 ! and his neighbor. 1498 vv(:)=(xatm(:,inxat)+xatm(:,ibat(ii))+atp(:,ipibat(ii)))/2._dp 1499 1500 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,-1) 1501 1502 if (ires==0) then 1503 ! Testing if CP is already known 1504 if (nbcp > 0) then 1505 do jj=1,nbcp 1506 pom(:)=vv(:)-bcp(jj)%rr(:)-xorig(:) 1507 dist=vnorm(pom,0) 1508 if (dist < aim_dtset%dpclim) then 1509 write(std_out,*) 'BCP already known !' 1510 cycle srbcp 1511 end if 1512 end do 1513 end if 1514 rr(:)=vv(:)-xorig(:) 1515 ss=vnorm(rr,0) 1516 if (ss > maxcpdst) then 1517 write(std_out, '(a,es16.6,a,es16.6)' ) 'BCP distance from atom,',ss,', exceed maxcpdst =',maxcpdst 1518 cycle srbcp 1519 end if 1520 nn=0 1521 do jj=1,3 1522 nn=nn+ev(jj)/abs(ev(jj)) 1523 end do 1524 write(std_out, '(a,3es16.6,i4)') ' vv(1:3), nn',(vv(jj), jj=1,3), nn 1525 write(std_out, '(a,3es16.6)') 'ev: ', (ev(jj), jj=1,3) 1526 if (nn /= -1) then 1527 write(std_out,*) ' The trial critical point is not a BCP !' 1528 cycle srbcp 1529 end if 1530 write(std_out, '(a,3es16.6)' ) 'evec(:,1): ',(evec(jj,1), jj=1,3) 1531 pom(:)=evec(:,1) 1532 dist=vnorm(pom,0) 1533 prj=dot_product(evec(:,1),rr) 1534 write(std_out,*) 'prj:', prj, vnorm(evec(:,1),0) 1535 dist=vnorm(evec(:,1),0) 1536 uu(:)=vv(:)-sign(aim_epstep,prj)*evec(:,1)/dist 1537 1538 ! Testing whether this BCP "is bonded" to the considered atom 1539 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1540 ! write(std_out,*) 'do', iat, ipos 1541 ! if ((iat==0).or.(ipos==0)) cycle 1542 ! write(std_out,*) 'APOS: ',(xatm(jj,iat)+atp(jj,ipos), jj=1,3) 1543 if ((iat/=inxat).or.(inxcell/=ipos)) then 1544 write(std_out,*) ' The trial BCP is not bonded to the Bader atom' 1545 cycle srbcp 1546 end if 1547 1548 ! A new BCP has been found ! 1549 nbcp=nbcp+1 1550 1551 ! Searching for the second bonded atom 1552 ss=vnorm(rr,0) 1553 diff1=ss 1554 diff3=dists(ii) 1555 uu(:)=vv(:)+sign(aim_epstep,prj)*evec(:,1)/dist 1556 if ((abs(bmin(iat))<1.0d-12).or.( ss<bmin(iat))) then 1557 bmin(iat)=ss 1558 end if 1559 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1560 if ((iat==0).or.(ipos==0)) then 1561 write(std_out,*) ' The trial BCP is not bonded to a bonding atom !' 1562 ! cycle srbcp 1563 end if 1564 pom(:)=vv(:)-xatm(:,iat)-atp(:,ipos) 1565 ss=vnorm(pom,0) 1566 diff2=ss 1567 pom(:)=xorig(:)-xatm(:,iat)-atp(:,ipos) 1568 diff3=vnorm(pom,0) 1569 rtdiff=diff1/diff3 1570 if ((abs(bmin(iat))<1.0d-12).or.(ss<bmin(iat))) then 1571 bmin(iat)=ss 1572 end if 1573 pom(:)=xatm(:,iat)+atp(:,ipos) 1574 1575 ! Store more results, for coherent, and portable output 1576 bcp(nbcp)%iat=iat 1577 bcp(nbcp)%ipos=ipos 1578 bcp(nbcp)%chg=chg 1579 bcp(nbcp)%diff(1)=diff1 1580 bcp(nbcp)%diff(2)=diff2 1581 bcp(nbcp)%diff(3)=diff3 1582 bcp(nbcp)%ev(:)=ev(:) 1583 bcp(nbcp)%pom(:)=pom(:) 1584 bcp(nbcp)%rr(:)=rr(:) 1585 bcp(nbcp)%vec(:,:)=evec(:,:) 1586 bcp(nbcp)%vv(:)=vv(:) 1587 ! Warning : iat, ipos might be modified by this call 1588 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1589 bcp(nbcp)%chg=chg 1590 1591 end if ! ires==0 1592 end do srbcp 1593 1594 if(nbcp>0)then 1595 1596 ! Order the BCP. CPs should appear by increasing values of x,y,z , the latter 1597 ! varying the fastest 1598 ABI_ALLOCATE(sortguide,(nbcp)) 1599 ABI_ALLOCATE(indexcp,(nbcp)) 1600 ABI_DATATYPE_ALLOCATE(cp_tmp,(nbcp)) 1601 do ii=3,1,-1 1602 ! DEBUG 1603 ! write(std_out,*)' cpdrv : sort on index ii=',ii 1604 ! ENDDEBUG 1605 1606 do jj=1,nbcp 1607 ! DEBUG 1608 ! write(std_out,*)bcp(jj)%vv(:) 1609 ! ENDDEBUG 1610 sortguide(jj)=bcp(jj)%vv(ii) 1611 indexcp(jj)=jj 1612 end do 1613 1614 ! Try to be platform-independent. Might need a larger tolerance. 1615 call sort_dp(nbcp,sortguide,indexcp,tol3) 1616 do jj=1,nbcp 1617 cp_tmp(jj)=bcp(indexcp(jj)) 1618 end do 1619 do jj=1,nbcp 1620 bcp(jj)=cp_tmp(jj) 1621 end do 1622 end do 1623 ! DEBUG 1624 ! write(std_out,*)' cpdrv : after the sort ' 1625 ! do jj=1,nbcp 1626 ! write(std_out,*)bcp(jj)%vv(:) 1627 ! end do 1628 ! ENDDEBUG 1629 1630 1631 ! Output the info about the BCP 1632 do jj=1,nbcp 1633 write(untout,'(" Bonded atom (BAT) (indxatm,indxcell,position): ",/,2I6,3F16.8)')& 1634 & bcp(jj)%iat,bcp(jj)%ipos,bcp(jj)%pom(:) 1635 write(untout,'("%Bonding CP: ",3F16.8)') bcp(jj)%vv(:) 1636 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') bcp(jj)%ev(:) 1637 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1638 & ' Eigenvec. of Hessian:',char(10),& 1639 & '-',bcp(jj)%vec(1,:),char(10),& 1640 & '-',bcp(jj)%vec(2,:),char(10),& 1641 & '-',bcp(jj)%vec(3,:),char(10) 1642 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1643 & bcp(jj)%chg, bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 1644 write(untout,'("%Relative position of BCP (AT-CP,BAT-CP,AT-BAT,relative(AT): ",/,4F16.8)') & 1645 & bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 1646 write(untout,*) "********************************************************************" 1647 write(std_out,'(/," BCP: ",3F10.6,3E12.4,E12.4,/)') & 1648 & bcp(jj)%rr(:),bcp(jj)%ev(:),bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 1649 write(std_out,'(":DISPC ",4F12.6)') bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 1650 end do 1651 1652 ABI_DATATYPE_DEALLOCATE(cp_tmp) 1653 ABI_DEALLOCATE(indexcp) 1654 ABI_DEALLOCATE(sortguide) 1655 1656 end if ! nbcp>0 1657 1658 if (abs(bmin(inxat))>1.0d-12) then 1659 rminl(inxat)=aim_dtset%coff1*bmin(inxat) 1660 r0=bmin(inxat) 1661 else 1662 r0=0._dp 1663 end if 1664 1665 ! !AD-HOC PARAMETER 1666 1667 do ii=1,natom 1668 if ((abs(bmin(ii))>1.0d-12).and.(ii /= inxat)) rminl(ii)=aim_dtset%coff2*bmin(ii) 1669 end do 1670 1671 ! END WARNING 1672 1673 write(std_out,*) ' number of BCP:', nbcp 1674 write(untout,'(" Number of BCP found: ",I4)') nbcp 1675 nn=nbcp*(nbcp-1)*(nbcp-2)/6 1676 if (bit_size(ii) <= nbcp+1) then 1677 MSG_ERROR("b-test!") 1678 end if 1679 1680 ! SEARCHING RCP 1681 1682 write(std_out,*) 1683 write(std_out,*) "RING CRITICAL POINTS (RCP)" 1684 write(std_out,*) "=============================" 1685 write(std_out,*) 1686 1687 write(untout,*) 1688 write(untout,*) "RING CRITICAL POINTS (RCP)" 1689 write(untout,*) "=============================" 1690 write(untout,*) 1691 1692 nrcp=0 1693 if(aim_dtset%crit==1)nb_now=nbcp 1694 if(aim_dtset%crit==2)nb_now=nb 1695 ! DEBUG 1696 ! nb_now=nbcp 1697 ! ENDDEBUG 1698 nn=nb_now*(nb_now-1)/2 1699 ABI_DATATYPE_ALLOCATE(rcp,(nn)) 1700 1701 ! Loop on pairs of BCP or atoms 1702 ipair=0 1703 ABI_ALLOCATE(buffer,(16,nn)) 1704 buffer=zero 1705 1706 ! DEBUG 1707 ! write(std_out,*)ch10,ch10,' drvcpr : enter loop to search for RCPs,nb_now,nn=',nb_now,nn 1708 ! ENDDEBUG 1709 1710 do ii=1,nb_now-1 1711 srcp1: do jj=ii+1,nb_now 1712 ipair=ipair+1 1713 if(mod(ipair,nproc)==me)then 1714 if (aim_dtset%crit==1) then 1715 vv(:)=xorig(:)+(bcp(ii)%rr(:)+bcp(jj)%rr(:))/2._dp 1716 else if (aim_dtset%crit==2) then 1717 vv(:)=xorig(:)*half+(xatm(:,ibat(ii))+atp(:,ipibat(ii))+xatm(:,ibat(jj))+atp(:,ipibat(jj)))*quarter 1718 end if 1719 1720 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,1) 1721 1722 if(ires==1)then 1723 cycle srcp1 1724 end if 1725 1726 ! Check that it is within the maximum allowed distance for a CP 1727 rr(:)=vv(:)-xorig(:) 1728 ss=vnorm(rr,0) 1729 if (ss > maxcpdst) then 1730 write(std_out,*) 'RCP distance from atom exceed maxcpdst !' 1731 cycle srcp1 1732 end if 1733 ! Check that it is a RCP 1734 nn=0 1735 do kk=1,3 1736 nn=nn+ev(kk)/abs(ev(kk)) 1737 end do 1738 if (nn /= 1) then 1739 write(std_out,*) ' the critical point that is found is not a RCP ' 1740 cycle srcp1 1741 end if 1742 ! Might be the same RCP than one already found on the same processor 1743 if (nrcp > 0) then 1744 do kk=1,nrcp 1745 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 1746 dist=vnorm(pom,0) 1747 if (dist < aim_dtset%dpclim) then 1748 write(std_out,*) ':RCP already known' 1749 cycle srcp1 1750 end if 1751 end do 1752 end if 1753 ! If crit==2, check that it is on the Bader surface 1754 if (aim_dtset%crit==2) then 1755 uu(:)=vv(:)-aim_epstep*rr(:)/ss 1756 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1757 if ((iat/=inxat).or.(inxcell/=ipos))then 1758 write(std_out,*) ' RCP is not on the Bader surface (outside of it)' 1759 cycle srcp1 1760 end if 1761 end if 1762 nrcp=nrcp+1 1763 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 1764 1765 buffer(1:3,ipair)=vv 1766 buffer(4:6,ipair)=ev 1767 buffer(7:9,ipair)=evec(:,1) 1768 buffer(10:12,ipair)=evec(:,2) 1769 buffer(13:15,ipair)=evec(:,3) 1770 buffer(16,ipair)=one 1771 1772 ! DEBUG 1773 ! write(std_out,*)ch10,ch10,' drvcpr : ipair,candidate=',ipair,candidate 1774 ! ENDDEBUG 1775 end if 1776 end do srcp1 1777 end do 1778 call xmpi_sum(buffer,xmpi_world,ierr) 1779 1780 nrcp=0 1781 ipair=0 1782 do ii=1,nb_now-1 1783 srcp: do jj=ii+1,nb_now 1784 ipair=ipair+1 1785 candidate=buffer(16,ipair) 1786 1787 ! One CP has been found, must make tests to see whether it is a new RCP 1788 if (nint(candidate)==1) then 1789 1790 vv=buffer(1:3,ipair) 1791 ev=buffer(4:6,ipair) 1792 evec(:,1)=buffer(7:9,ipair) 1793 evec(:,2)=buffer(10:12,ipair) 1794 evec(:,3)=buffer(13:15,ipair) 1795 1796 ! Check that it is not the same as a previous one 1797 if (nrcp > 0) then 1798 do kk=1,nrcp 1799 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 1800 dist=vnorm(pom,0) 1801 if (dist < aim_dtset%dpclim) then 1802 write(std_out,*) ':RCP already known' 1803 cycle srcp 1804 end if 1805 end do 1806 end if 1807 1808 ! A new RCP has been found ! 1809 nrcp=nrcp+1 1810 1811 ! DEBUG 1812 ! write(std_out,*)' drvcpr : A new RCP has been found, for kk=',kk 1813 ! ENDDEBUG 1814 1815 1816 rcp(nrcp)%iat=iat 1817 rcp(nrcp)%ipos=ipos 1818 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 1819 rcp(nrcp)%vec(:,:)=evec(:,:) 1820 rcp(nrcp)%ev(:)=ev(:) 1821 rcp(nrcp)%vv(:)=vv(:) 1822 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1823 rcp(nrcp)%chg=chg 1824 1825 end if ! ires==0 1826 end do srcp ! jj=ii+2,nb_now 1827 end do ! ii=1,nb_now-1 1828 1829 ABI_DEALLOCATE(buffer) 1830 1831 if(nrcp>0)then 1832 1833 ! Order the RCP. CPs should appear by increasing values of x,y,z , the latter 1834 ! varying the fastest 1835 ABI_ALLOCATE(sortguide,(nrcp)) 1836 ABI_ALLOCATE(indexcp,(nrcp)) 1837 ABI_DATATYPE_ALLOCATE(cp_tmp,(nrcp)) 1838 do ii=3,1,-1 1839 ! DEBUG 1840 ! write(std_out,*)' cpdrv : sort on index ii=',ii 1841 ! ENDDEBUG 1842 do jj=1,nrcp 1843 1844 ! DEBUG 1845 ! write(std_out,*)rcp(jj)%vv(:) 1846 ! ENDDEBUG 1847 1848 ! Try to be platform-independent. Might need a larger tolerance. 1849 sortguide(jj)=rcp(jj)%vv(ii) 1850 indexcp(jj)=jj 1851 end do 1852 call sort_dp(nrcp,sortguide,indexcp,tol3) 1853 do jj=1,nrcp 1854 cp_tmp(jj)=rcp(indexcp(jj)) 1855 end do 1856 do jj=1,nrcp 1857 rcp(jj)=cp_tmp(jj) 1858 end do 1859 end do 1860 1861 ! DEBUG 1862 ! write(std_out,*)' cpdrv : after the sort ' 1863 ! do jj=1,nrcp 1864 ! write(std_out,*)rcp(jj)%vv(:) 1865 ! end do 1866 ! ENDDEBUG 1867 1868 1869 ! Write the Ring Critical Point information 1870 do jj=1,nrcp 1871 write(untout,'(";Ring CP: ",3F16.8)') rcp(jj)%vv(:) 1872 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') rcp(jj)%ev(:) 1873 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1874 & ' Eigenvec. of Hessian:',char(10),& 1875 & '-',rcp(jj)%vec(1,:),char(10),& 1876 & '-',rcp(jj)%vec(2,:),char(10),& 1877 & '-',rcp(jj)%vec(3,:),char(10) 1878 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1879 & rcp(jj)%chg, rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 1880 write(untout,*) "********************************************************************" 1881 write(std_out,'(/," RCP: ",3F10.6,3E12.4,E12.4,/)') & 1882 & rcp(jj)%rr(:),rcp(jj)%ev(:),rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 1883 end do 1884 1885 ABI_DATATYPE_DEALLOCATE(cp_tmp) 1886 ABI_DEALLOCATE(indexcp) 1887 ABI_DEALLOCATE(sortguide) 1888 1889 end if ! nrcp>0 1890 1891 write(untout,'(" Number of RCP found: ",I4)') nrcp 1892 write(std_out,*) ' Number of RCP:', nrcp 1893 1894 ! SEARCHING CCP 1895 1896 write(std_out,*) 1897 write(std_out,*) "CAGE CRITICAL POINTS (CCP)" 1898 write(std_out,*) "=============================" 1899 write(std_out,*) 1900 1901 write(untout,*) 1902 write(untout,*) "CAGE CRITICAL POINTS (CCP)" 1903 write(untout,*) "=============================" 1904 write(untout,*) 1905 1906 1907 nn=nrcp*(nrcp-1)/2 1908 ABI_DATATYPE_ALLOCATE(ccp,(nn)) 1909 1910 nccp=0 1911 do ii=1,nrcp-1 1912 srccp: do jj=ii+1,nrcp 1913 vv(:)=xorig(:)+(rcp(ii)%rr(:)+rcp(jj)%rr(:))/2._dp 1914 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,3) 1915 if (ires==0) then 1916 rr(:)=vv(:)-xorig(:) 1917 ss=vnorm(rr,0) 1918 if (ss > maxcpdst) then 1919 write(std_out,*) 'CCP distance from atom exceed maxcpdst !' 1920 cycle srccp 1921 end if 1922 nn=0 1923 do kk=1,3 1924 nn=nn+ev(kk)/abs(ev(kk)) 1925 end do 1926 if (nn /= 3) then 1927 write(std_out,*) ' the critical point that is found is not a CCP ' 1928 cycle srccp 1929 end if 1930 1931 if (nccp > 0) then 1932 do kk=1,nccp 1933 pom(:)=vv(:)-ccp(kk)%rr(:)-xorig(:) 1934 dist=vnorm(pom,0) 1935 if (dist < aim_dtset%dpclim) then 1936 write(std_out,*) ':CCP already known' 1937 cycle srccp 1938 end if 1939 end do 1940 end if 1941 if (aim_dtset%crit==2) then 1942 uu(:)=vv(:)-aim_epstep*rr(:)/ss 1943 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 1944 if ((iat/=inxat).or.(inxcell/=ipos)) then 1945 write(std_out,*) ' This CCP is not on the Bader surface (outside of it)' 1946 cycle srccp 1947 end if 1948 end if 1949 1950 nccp=nccp+1 1951 1952 ccp(nccp)%iat=iat 1953 ccp(nccp)%ipos=ipos 1954 ccp(nccp)%rr(:)=vv(:)-xorig(:) 1955 ccp(nccp)%vec(:,:)=evec(:,:) 1956 ccp(nccp)%ev(:)=ev(:) 1957 ccp(nccp)%vv(:)=vv(:) 1958 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 1959 ccp(nccp)%chg=chg 1960 1961 end if 1962 end do srccp 1963 end do 1964 1965 if(nccp>0)then 1966 1967 ! Order the CCP. CPs should appear by increasing values of x,y,z , the latter 1968 ! varying the fastest 1969 ABI_ALLOCATE(sortguide,(nccp)) 1970 ABI_ALLOCATE(indexcp,(nccp)) 1971 ABI_DATATYPE_ALLOCATE(cp_tmp,(nccp)) 1972 do ii=3,1,-1 1973 do jj=1,nccp 1974 ! Try to be platform-independent. Might need a larger tolerance. 1975 sortguide(jj)=ccp(jj)%vv(ii) 1976 indexcp(jj)=jj 1977 end do 1978 call sort_dp(nccp,sortguide,indexcp,tol3) 1979 do jj=1,nccp 1980 cp_tmp(jj)=ccp(indexcp(jj)) 1981 end do 1982 do jj=1,nccp 1983 ccp(jj)=cp_tmp(jj) 1984 end do 1985 end do 1986 1987 ! Write the Cage Critical Point information 1988 do jj=1,nccp 1989 write(untout,'("%Cage CP: ",3F16.8)') ccp(jj)%vv(:) 1990 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') ccp(jj)%ev(:) 1991 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 1992 & ' Eigenvec. of Hessian:',char(10),& 1993 & '-',ccp(jj)%vec(1,:),char(10),& 1994 & '-',ccp(jj)%vec(2,:),char(10),& 1995 & '-',ccp(jj)%vec(3,:),char(10) 1996 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 1997 & ccp(jj)%chg, ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 1998 write(untout,*) "********************************************************************" 1999 write(std_out,'(/," CCP: ",3F10.6,3E12.4,E12.4,/)') & 2000 & ccp(jj)%rr(:),ccp(jj)%ev(:),ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 2001 end do 2002 2003 ABI_DEALLOCATE(sortguide) 2004 ABI_DEALLOCATE(indexcp) 2005 ABI_DATATYPE_DEALLOCATE(cp_tmp) 2006 2007 end if ! nccp>0 2008 2009 write(untout,'(" Number of CCP found: ",I4)') nccp 2010 write(std_out,*) 'Number of CCP:', nccp 2011 write(std_out,*) 2012 write(untout,*) 2013 write(std_out, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 2014 write(untout, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 2015 2016 write(std_out,*) 2017 write(std_out,*) "===============================" 2018 write(std_out,*) "END OF CRITICAL POINTS ANALYSIS" 2019 write(std_out,*) 2020 2021 write(untout,*) 2022 write(untout,*) "===============================" 2023 write(untout,*) "END OF CRITICAL POINTS ANALYSIS" 2024 write(untout,*) 2025 2026 2027 ! Output of the CPs 2028 2029 write(untc,'(I4, " :BCP''s, coordinates, laplacian eigs, type of bonding at., sum of lap.eigs., density")') nbcp 2030 do ii=1,nbcp 2031 write(untc,'(3F10.6,3E12.4,I4,2E12.4)') & 2032 & bcp(ii)%rr(:),bcp(ii)%ev(:),bcp(ii)%iat,bcp(ii)%ev(1)+bcp(ii)%ev(2)+bcp(ii)%ev(3),bcp(ii)%chg 2033 end do 2034 2035 write(untc,'(I4, " :RCP''s, coordinates, laplacian eigenvalues, sum of these, density")') nrcp 2036 do ii=1,nrcp 2037 vv(:)=rcp(ii)%rr(:)+xorig(:) 2038 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 2039 write(untc,'(3F10.6,3E12.4,2E12.4)') & 2040 & rcp(ii)%rr(:),rcp(ii)%ev(:),rcp(ii)%ev(1)+rcp(ii)%ev(2)+rcp(ii)%ev(3),rcp(ii)%chg 2041 end do 2042 2043 write(untc,'(I4, " :CCP''s coordinates, laplacian eigenvalues, sum of these, density")') nccp 2044 do ii=1,nccp 2045 vv(:)=ccp(ii)%rr(:)+xorig(:) 2046 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 2047 write(untc,'(3F10.6,3E12.4,2E12.4)') & 2048 & ccp(ii)%rr(:),ccp(ii)%ev(:),ccp(ii)%ev(1)+ccp(ii)%ev(2)+ccp(ii)%ev(3),ccp(ii)%chg 2049 end do 2050 2051 end if ! End the condition on aim_dtset%crit > 0 2052 2053 !Reading of the CPs from the file 2054 2055 if (aim_dtset%crit==-1) then 2056 read(untc,*) nbcp 2057 ABI_DATATYPE_ALLOCATE(bcp,(nbcp)) 2058 do ii=1,nbcp 2059 read(untc,*) bcp(ii)%rr(:) 2060 end do 2061 read(untc,*) nrcp 2062 ABI_DATATYPE_ALLOCATE(rcp,(nrcp)) 2063 do ii=1,nrcp 2064 read(untc,*) rcp(ii)%rr(:) 2065 end do 2066 read(untc,*) nccp 2067 ABI_DATATYPE_ALLOCATE(ccp,(nccp)) 2068 do ii=1,nccp 2069 read(untc,*) ccp(ii)%rr(:) 2070 end do 2071 end if 2072 2073 do ii=1,nbcp 2074 pc(:,ii)=bcp(ii)%rr(:) 2075 icpc(ii)=-1 2076 end do 2077 do ii=1,nrcp 2078 pc(:,nbcp+ii)=rcp(ii)%rr(:) 2079 icpc(nbcp+ii)=1 2080 end do 2081 do ii=1,nccp 2082 pc(:,nbcp+nrcp+ii)=ccp(ii)%rr(:) 2083 icpc(nbcp+nrcp+ii)=3 2084 end do 2085 npc=nbcp+nrcp+nccp 2086 2087 !Checking 2088 2089 if (allocated(bcp)) then 2090 do ii=1,nbcp 2091 do jj=1,npc 2092 iat=bcp(ii)%iat 2093 ipos=bcp(ii)%ipos 2094 if ((iat/=0).and.(ipos/=0)) then 2095 pom(:)=pc(:,jj)+xorig(:)-xatm(:,iat)-atp(:,ipos) 2096 ss=aim_dtset%coff2*vnorm(pom,0) 2097 if (rminl(iat) >= ss) rminl(iat)=ss 2098 end if 2099 end do 2100 end do 2101 ABI_DATATYPE_DEALLOCATE(bcp) 2102 end if 2103 do ii=1,natom 2104 write(std_out,*) 'atom: ', ii, rminl(ii) 2105 end do 2106 2107 if(allocated(rcp)) then 2108 ABI_DATATYPE_DEALLOCATE(rcp) 2109 end if 2110 if(allocated(ccp)) then 2111 ABI_DATATYPE_DEALLOCATE(ccp) 2112 end if 2113 2114 !END CP ANALYSIS 2115 2116 call timein(ttcp,wall) 2117 ttcp=ttcp-tt0 2118 2119 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
PARENTS
aim_follow,cpdrv,critics
CHILDREN
SOURCE
2151 subroutine critic(aim_dtset,vv,ev,zz,dmax,ires,sort) 2152 2153 2154 !This section has been created automatically by the script Abilint (TD). 2155 !Do not modify the following lines by hand. 2156 #undef ABI_FUNC 2157 #define ABI_FUNC 'critic' 2158 !End of the abilint section 2159 2160 implicit none 2161 2162 !Arguments ------------------------------------ 2163 !scalars 2164 integer,intent(in) :: sort 2165 integer,intent(out) :: ires 2166 real(dp),intent(in) :: dmax 2167 !arrays 2168 real(dp),intent(inout) :: vv(3) 2169 real(dp),intent(out) :: ev(3),zz(3,3) 2170 !no_abirules 2171 type(aim_dataset_type), intent(in) :: aim_dtset 2172 2173 !Local variables ------------------------------ 2174 !scalars 2175 integer :: iat,id,ii,info,ipos,istep,jii,jj,nrot 2176 real(dp),parameter :: evol=1.d-3 2177 real(dp) :: chg,dg,dltcmax,dv,dvold,rr,ss 2178 logical :: oscl,outof 2179 !arrays 2180 integer :: ipiv(3) 2181 real(dp) :: dc(3),ff(3),grho(3),hrho(3,3),lp(3),vold(3),vt(3),yy(3,3) 2182 real(dp),allocatable :: lamb(:),pom(:,:),pom2(:,:) 2183 2184 !************************************************************************ 2185 2186 !DEBUG 2187 !write(std_out,*)' critic : enter ' 2188 !ENDDEBUG 2189 oscl=.false. 2190 if (sort==3) then 2191 ABI_ALLOCATE(pom,(4,4)) 2192 ABI_ALLOCATE(pom2,(4,4)) 2193 ABI_ALLOCATE(lamb,(4)) 2194 elseif (sort/=0) then 2195 ABI_ALLOCATE(pom,(3,3)) 2196 ABI_ALLOCATE(pom2,(3,3)) 2197 ABI_ALLOCATE(lamb,(3)) 2198 end if 2199 2200 2201 deb=.false. 2202 istep=0 2203 ires=0 2204 2205 !DEBUG 2206 !write(std_out,'(":POSIN ",3F16.8)') vv 2207 !do jj=1,3 2208 !vt(jj)=rprimd(1,jj)*vv(1)+rprimd(2,jj)*vv(2)+rprimd(3,jj)*vv(3) 2209 !end do 2210 !write(std_out,'(":RBPOSIN ",3F16.8)') vt 2211 !ENDDEBUG 2212 2213 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2214 2215 !write(std_out,'(":GRAD ",3F16.8)') grho 2216 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 2217 2218 dg=1.0_dp 2219 dv=1.0_dp 2220 dg = vnorm(grho,0) 2221 2222 if (chg < aim_rhomin) then 2223 ires=1 2224 ! DEBUG 2225 ! write(std_out,*)' critic : exit, ires=1' 2226 ! ENDDEBUG 2227 return 2228 end if 2229 2230 !main cycle => limits (adhoc): 2231 !aim_dtset%lstep - minimal step 2232 !aim_dtset%lgrad - minimal norm of gradient 2233 !aim_maxstep - max number of steps 2234 2235 do while ((dv>aim_dtset%lstep).and.(dg>aim_dtset%lgrad).and.(istep<aim_maxstep)) 2236 istep=istep+1 2237 vold(:)=vv(:) 2238 dvold=dv 2239 ev(:)=0._dp 2240 yy(:,:)=0._dp 2241 call jacobi(hrho,3,3,ev,yy,nrot) ! eigenval of Hessian 2242 call ordr(ev,yy,3,-1) ! ordering 2243 2244 ! modification of the Newton-Raphson step to searching 2245 ! specific type of CP (Popelier algorithm) 2246 2247 ff(:)=0._dp 2248 lp(:)=0._dp 2249 dc(:)=0._dp 2250 outof=.false. 2251 do ii=1,3 2252 do jj=1,3 2253 ff(ii)=ff(ii)+yy(jj,ii)*grho(jj) 2254 end do 2255 end do 2256 id=sign(1._dp,ev(1))+sign(1._dp,ev(2))+sign(1._dp,ev(3)) 2257 if (id /= sort) then 2258 outof=.true. 2259 select case (sort) 2260 case (-1) 2261 lp(3)=0.5_dp*(ev(3)-sqrt(ev(3)*ev(3)+4._dp*ff(3)*ff(3))) 2262 pom(:,:)=0._dp 2263 pom2(:,:)=0._dp 2264 lamb(:)=0._dp 2265 do ii=1,2 2266 pom(ii,ii)=ev(ii) 2267 pom(ii,3)=ff(ii) 2268 pom(3,ii)=ff(ii) 2269 end do 2270 call jacobi(pom,3,3,lamb,pom2,nrot) 2271 call ordr(lamb,pom2,3,1) 2272 do ii=1,3 2273 lp(1)=lamb(ii) 2274 if (abs(pom2(3,ii))>1.0d-24) exit 2275 end do 2276 lp(2)=lp(1) 2277 2278 ! write(std_out,*) (ev(ii),ii=1,3) 2279 ! write(std_out,*) (lamb(ii),ii=1,3) 2280 ! write(std_out,*) ':ID ',id,lp(1),lp(3) 2281 2282 case (1) 2283 lp(1)=0.5_dp*(ev(1)+sqrt(ev(1)*ev(1)+4._dp*ff(1)*ff(1))) 2284 pom(:,:)=0._dp 2285 pom2(:,:)=0._dp 2286 lamb(:)=0._dp 2287 do ii=2,3 2288 pom(ii-1,ii-1)=ev(ii) 2289 pom(ii-1,3)=ff(ii) 2290 pom(3,ii-1)=ff(ii) 2291 end do 2292 call jacobi(pom,3,3,lamb,pom2,nrot) 2293 call ordr(lamb,pom2,3,1) 2294 do ii=3,1,-1 2295 lp(2)=lamb(ii) 2296 if (abs(pom2(3,ii))>1.0d-24) exit 2297 end do 2298 lp(3)=lp(2) 2299 2300 case (3) 2301 pom(:,:)=0._dp 2302 pom2(:,:)=0._dp 2303 lamb(:)=0._dp 2304 do ii=1,3 2305 pom(ii,ii)=ev(ii) 2306 pom(ii,4)=ff(ii) 2307 pom(4,ii)=ff(ii) 2308 end do 2309 call jacobi(pom,4,4,lamb,pom2,nrot) 2310 call ordr(lamb,pom2,4,1) 2311 do ii=4,1,-1 2312 lp(1)=lamb(ii) 2313 if (abs(pom2(4,ii))>1.0d-24) exit 2314 end do 2315 lp(2)=lp(1); lp(3)=lp(1) 2316 case default 2317 lp(:)=0._dp 2318 end select 2319 end if 2320 2321 do ii=1,3 2322 if (abs(ev(ii)-lp(ii))<1.0d-24) then 2323 outof=.false. 2324 exit 2325 end if 2326 end do 2327 do ii=1,3 ! SEARCHING STEP 2328 do jj=1,3 2329 if (outof) then 2330 dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/(ev(jj)-lp(jj)) 2331 elseif (abs(ev(jj))>1.0d-24) then 2332 dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/ev(jj) 2333 else 2334 MSG_ERROR("zero eigval of Hessian") 2335 end if 2336 end do 2337 end do 2338 2339 dltcmax = vnorm(dc,0) 2340 if (dltcmax>dmax) then ! STEP RESTRICTION 2341 do ii=1,3 2342 dc(ii)=dc(ii)*dmax/dltcmax 2343 end do 2344 end if ! primitive handling of oscillations 2345 ss=vnorm(dc,0) ! usually not needed 2346 ss=abs(ss-dv)/ss 2347 if ((ss < evol).and.(oscl)) then 2348 dc(:)=dc(:)/2._dp 2349 end if 2350 2351 2352 do ii=1,3 2353 vv(ii) = vv(ii) - dc(ii) 2354 end do 2355 2356 ! DEBUG 2357 ! write(std_out,'(":POSIN ",3F16.8)') vv 2358 ! ENDDEBUG 2359 2360 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2361 dg = vnorm(grho,0) 2362 2363 if (deb) then ! DEBUGG OUTPUT 2364 write(std_out,'("AFTER STEP ===================================")') 2365 write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((yy(ii,jii),jii=1,3),ii=1,3) 2366 write(std_out,'(":DC ",3F16.8)') dc 2367 write(std_out,*) 'STEP ',istep 2368 write(std_out,'(":POS ",3F16.8)') vv 2369 write(std_out,'(":GRAD ",3F16.8)') grho 2370 write(std_out,*) ':DGRAD,CHG ',dg,chg 2371 write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jii),jii=1,3),ii=1,3) 2372 end if 2373 vt(:)=vv(:)-vold(:) 2374 dv=vnorm(vt,0) 2375 ss=abs(dvold-dv)/dv 2376 if (ss < evol) oscl=.true. 2377 end do 2378 2379 !end of main cycle 2380 2381 !the final output 2382 2383 write(std_out,*) 'iste:',istep, dv, dg 2384 if (istep>=aim_maxstep)then 2385 write(std_out,*) ' istep=MAXSTEP ! Examine lstep2 and lgrad2 .' 2386 if ( (dv>aim_dtset%lstep2) .and. (dg>aim_dtset%lgrad2 )) then 2387 write(std_out,'(":POSOUT ",3F16.8)') vv 2388 ires=1 2389 end if 2390 end if 2391 2392 vt(:)=vv(:) 2393 2394 !write(std_out,'(":POSOUT ",3F16.8)') vv 2395 !write(std_out,'(":RBPOSOUT ",3F16.8)') vt 2396 2397 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2398 2399 !write(std_out,'(":GRAD ",3F16.8)') grho 2400 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)')& 2401 !& ((hrho(ii,jii),jii=1,3),ii=1,3) 2402 2403 2404 !FINAL INVERSION OF HESSIAN 2405 2406 call ludcmp(hrho,3,3,ipiv,id,info) 2407 if (info /= 0) then 2408 write(std_out,*) 'Error inverting hrho:' 2409 do ii=1,3 2410 write(std_out,*) (hrho(ii,jii),jii=1,3) 2411 end do 2412 ires=1 2413 ! DEBUG 2414 ! write(std_out,*)' critic : exit, ires=1' 2415 ! ENDDEBUG 2416 return 2417 ! stop 'ERROR INVERTING HESSIAN' 2418 end if 2419 do ii=1,3 2420 yy(ii,1:3)=0. 2421 yy(ii,ii)=1. 2422 end do 2423 do jii=1,3 2424 call lubksb(hrho,3,3,ipiv,yy(1,jii)) 2425 end do 2426 2427 2428 !write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((y(ii,jii),jii=1,3),ii=1,3) 2429 2430 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 2431 2432 !write(std_out,'("LAPLAC:",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 2433 2434 call jacobi(hrho,3,3,ev,yy,nrot) 2435 call ordr(ev,yy,3,1) 2436 zz(:,:)=yy(:,:) 2437 2438 !do ii=1,3 2439 !do jii=1,3 2440 !zz(ii,jii)=yy(jii,ii) 2441 !end do 2442 !end do 2443 2444 !write(std_out,'(":AUTOVAL ",3F16.8)') (ev(ii),ii=1,3) 2445 !write(std_out,'(":AUTOVEC ",/,3F16.8,/,3F16.8,/,3F16.8)') ((zz(ii,jii),ii=1,3),jii=1,3) 2446 2447 if (sort/=0) then 2448 ABI_DEALLOCATE(pom) 2449 ABI_DEALLOCATE(pom2) 2450 ABI_DEALLOCATE(lamb) 2451 end if 2452 2453 !DEBUG 2454 !write(std_out,*)' critic : exit, ires= ',ires 2455 !ENDDEBUG 2456 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)
PARENTS
drvaim
CHILDREN
bschg1,critic,sort_dp,vgh_rho
SOURCE
2565 subroutine critics(aim_dtset,inxat,stwo,sthree,sfour,dstmax) 2566 2567 2568 !This section has been created automatically by the script Abilint (TD). 2569 !Do not modify the following lines by hand. 2570 #undef ABI_FUNC 2571 #define ABI_FUNC 'critics' 2572 !End of the abilint section 2573 2574 implicit none 2575 2576 !Arguments ------------------------------------ 2577 !scalars 2578 integer,intent(in) :: inxat 2579 real(dp),intent(in) :: dstmax 2580 logical,intent(in) :: sfour,sthree,stwo 2581 !no_abirules 2582 type(aim_dataset_type), intent(in) :: aim_dtset 2583 2584 !Local variables ------------------------------ 2585 !scalars 2586 integer :: i1,i2,i3,iat,ii,ipos,ires,jii,jj,kjj,kk,ll,n1,n2,n3,nb,nc 2587 integer :: nshell 2588 real(dp) :: chg,dif1,dif2,diff,dist,olddist,rr 2589 ! real(dp) :: ss,uu 2590 logical :: found,inter 2591 !arrays 2592 integer :: ibat(nnpos*natom),inat(nnpos*natom),ipibat(nnpos*natom) 2593 integer :: nnat(nnpos*natom),nr(nnpos*natom) 2594 real(dp) :: dif(3),dists(nnpos*natom),ev(3),grho(3),hrho(3,3) 2595 real(dp) :: pom(3),v1(3),v2(3),v3(3),v4(3),vi(3),vt(3),zz(3,3) 2596 2597 !************************************************************************ 2598 vi(:)=xatm(:,inxat) 2599 2600 nc=0 2601 do jii=1,nnpos 2602 do kjj=1,natom 2603 dist=0._dp 2604 dif(:)=xatm(:,inxat)-xatm(:,kjj)-atp(:,jii) 2605 2606 ! do ii=1,3 2607 ! dif(ii)=xatm(ii,inxat)-xatm(ii,kjj)-atp(ii,jii) 2608 ! end do 2609 dist=vnorm(dif,0) 2610 if (.not.((dist>dstmax).or.(dist<0.001))) then 2611 nc=nc+1 2612 dists(nc)=dist 2613 nnat(nc)=kjj 2614 inat(nc)=jii 2615 end if 2616 end do 2617 end do 2618 do n1=1,nc 2619 nr(n1)=n1 2620 end do 2621 call sort_dp(nc,dists,nr,tol14) 2622 nb=0 2623 olddist=0._dp 2624 nshell=0 2625 !write(std_out,*) ':ORIAT ', (xatm(ii,inxat),ii=1,3) 2626 do n1=1,nc 2627 n2=nr(n1) 2628 n3=nnat(n2) 2629 if (dists(n1)<(2*dists(1))) then 2630 if ((dists(n1)-olddist)>aim_dlimit) then 2631 nshell=nshell+1 2632 olddist=dists(n1) 2633 if (nshell==5) exit 2634 end if 2635 nb=nb+1 2636 ibat(nb)=n3 2637 ipibat(nb)=inat(n2) 2638 write(std_out,*) ':NEIG ',inxat,n3,inat(n2),dists(n1) 2639 ! write(std_out,*) ':POSAT',(xatm(ii,ibat(nb))+atp(ii,ipibat(nb)),ii=1,3) 2640 else 2641 exit 2642 end if 2643 end do 2644 2645 npc=0 2646 npcm3=0 2647 2648 ! 2649 !.....SEARCH BETWEEN EACH PAIR OF ATOMS 2650 ! 2651 2652 if (stwo) then 2653 do jii=1,nb 2654 do ii=1,3 2655 v1(ii)=xatm(ii,inxat) 2656 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2657 vt(ii)=(v1(ii)+v2(ii))/2._dp 2658 end do 2659 inter=.true. 2660 diff=0._dp 2661 pom(:)=vt(:) 2662 pom(:)=pom(:)-vi(:) 2663 diff=vnorm(pom,0) 2664 if (diff > maxcpdst) inter=.false. 2665 if (inter) then 2666 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2667 if (ires==0) then 2668 found=.false. 2669 if (npc > 0) then 2670 do jj=1,npc 2671 pom(:)=vt(:)-pc(:,jj) 2672 dist=vnorm(pom,0) 2673 if (dist < aim_dtset%dpclim) found=.true. 2674 end do 2675 end if 2676 if (.not.found) then 2677 pom(:)=vt(:) 2678 call bschg1(pom,-1) 2679 pcrb(:,npc+1)=pom(:) 2680 pom(:)=pom(:)-vi(:) 2681 diff=vnorm(pom,0) 2682 if (abs(diff) > maxcpdst) found=.true. 2683 end if 2684 if (.not.found) then 2685 npc=npc+1 2686 do jj=1,3 2687 pc(jj,npc)=vt(jj) 2688 evpc(jj,npc)=ev(jj) 2689 do kk=1,3 2690 zpc(kk,jj,npc)=zz(kk,jj) 2691 end do 2692 end do 2693 i1=ev(1)/abs(ev(1)) 2694 i2=ev(2)/abs(ev(2)) 2695 i3=ev(3)/abs(ev(3)) 2696 icpc(npc)=i1+i2+i3 2697 if (icpc(npc)==-3) then 2698 npcm3=npcm3+1 2699 end if 2700 write(std_out,*) 'New critical point found' 2701 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2702 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2703 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2704 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2705 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2706 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2707 write(22,'(":PC2",3F10.6,3E12.4,I4,2E12.4)') & 2708 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2709 write(std_out,'(":PC2",3F10.6,3E12.4,I4,2E12.4)') & 2710 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2711 pom(:)=vt(:)-v1(:) 2712 dif1=vnorm(pom,0) 2713 pom(:)=vt(:)-v2(:) 2714 dif2=vnorm(pom,0) 2715 write(std_out,'(":DISPC ",2F12.8)') dif1,dif2 2716 end if 2717 end if 2718 end if 2719 end do 2720 end if 2721 ! 2722 !.....SEARCH BETWEEN EACH THREE ATOMS 2723 ! 2724 if(sthree) then 2725 do jii=1,nb 2726 do kjj=jii+1,nb 2727 do ii=1,3 2728 v1(ii)=xatm(ii,inxat) 2729 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2730 v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj)) 2731 vt(ii)=(v1(ii)+v2(ii)+v3(ii))/3._dp 2732 end do 2733 inter=.true. 2734 pom(:)=vt(:) 2735 pom(:)=pom(:)-vi(:) 2736 dist=vnorm(pom,0) 2737 if (abs(diff)>maxcpdst) then 2738 inter=.false. 2739 exit 2740 end if 2741 if (inter) then 2742 do jj=1,npc 2743 pom(:)=pc(:,jj)-vt(:) 2744 diff=vnorm(pom,0) 2745 if (diff<aim_dpc0) then 2746 inter=.false. 2747 exit 2748 end if 2749 end do 2750 end if 2751 if (inter) then 2752 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2753 if (ires==0) then 2754 found=.false. 2755 if (npc>0) then 2756 do jj=1,npc 2757 pom(:)=vt(:)-pc(:,jj) 2758 dist=vnorm(pom,0) 2759 if (dist<aim_dtset%dpclim) then 2760 found=.true. 2761 exit 2762 end if 2763 end do 2764 end if 2765 if (.not.found) then 2766 pom(:)=vt(:) 2767 call bschg1(pom,-1) 2768 pcrb(:,npc+1)=pom(:) 2769 pom(:)=pom(:)-vi(:) 2770 diff=vnorm(pom,0) 2771 if (abs(diff)>maxcpdst) found=.true. 2772 end if 2773 if (.not.found) then 2774 npc=npc+1 2775 do jj=1,3 2776 pc(jj,npc)=vt(jj) 2777 evpc(jj,npc)=ev(jj) 2778 do kk=1,3 2779 zpc(kk,jj,npc)=zz(kk,jj) 2780 end do 2781 end do 2782 i1=ev(1)/abs(ev(1)) 2783 i2=ev(2)/abs(ev(2)) 2784 i3=ev(3)/abs(ev(3)) 2785 icpc(npc)=i1+i2+i3 2786 if (icpc(npc)==-3) then 2787 npcm3=npcm3+1 2788 end if 2789 write(std_out,*) 'New critical point found' 2790 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2791 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2792 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2793 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2794 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2795 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2796 write(22,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') & 2797 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2798 write(std_out,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') & 2799 & (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2800 end if 2801 end if 2802 end if 2803 end do 2804 end do 2805 end if 2806 2807 ! 2808 !.....SEARCH BETWEEN EACH FOUR ATOMS 2809 ! 2810 if (sfour) then 2811 do jii=1,nb 2812 do kjj=jii+1,nb 2813 do ll=jii+1,nb 2814 do ii=1,3 2815 v1(ii)=xatm(ii,inxat) 2816 v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii)) 2817 v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj)) 2818 v4(ii)=xatm(ii,ibat(ll))+atp(ii,ipibat(ll)) 2819 vt(ii)=(v1(ii)+v2(ii)+v3(ii)+v4(ii))/4._dp 2820 end do 2821 inter=.true. 2822 pom(:)=vt(:) 2823 pom(:)=pom(:)-vi(:) 2824 diff=vnorm(pom,0) 2825 if (abs(diff)>maxcpdst) then 2826 inter=.false. 2827 exit 2828 end if 2829 if (inter) then 2830 do jj=1,npc 2831 pom(:)=pc(:,jj)-vt(:) 2832 diff=vnorm(pom,0) 2833 if (diff < aim_dpc0) then 2834 inter=.false. 2835 exit 2836 end if 2837 end do 2838 end if 2839 if (inter) then 2840 call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0) 2841 if (ires==0) then 2842 found=.false. 2843 if (npc>0) then 2844 do jj=1,npc 2845 pom(:)=vt(:)-pc(:,jj) 2846 dist=vnorm(pom,0) 2847 if (dist < aim_dtset%dpclim) found=.true. 2848 end do 2849 end if 2850 if (.not.found) then 2851 pom(:)=vt(:) 2852 pcrb(:,npc+1)=pom(:) 2853 pom(:)=pom(:)-vi(:) 2854 diff=vnorm(pom,0) 2855 if (abs(diff)>maxcpdst) found=.true. 2856 end if 2857 if (.not.found) then 2858 npc=npc+1 2859 do jj=1,3 2860 pc(jj,npc)=vt(jj) 2861 evpc(jj,npc)=ev(jj) 2862 do kk=1,3 2863 zpc(kk,jj,npc)=zz(kk,jj) 2864 end do 2865 end do 2866 i1=ev(1)/abs(ev(1)) 2867 i2=ev(2)/abs(ev(2)) 2868 i3=ev(3)/abs(ev(3)) 2869 icpc(npc)=i1+i2+i3 2870 if (icpc(npc)==-3) then 2871 npcm3=npcm3+1 2872 end if 2873 write(std_out,*) 'New critical point found' 2874 write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3) 2875 write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3) 2876 write(std_out,'("AUTOVAL: ",3F16.8)') ev 2877 write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 2878 & ((zpc(ii,jj,npc),ii=1,3),jj=1,3) 2879 call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0) 2880 write(22,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') & 2881 & (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2882 write(std_out,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') & 2883 & (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg 2884 end if 2885 end if 2886 end if 2887 end do 2888 end do 2889 end do 2890 end if 2891 2892 write(std_out,*) npc 2893 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
PARENTS
aim
CHILDREN
SOURCE
2916 subroutine defad(aim_dtset) 2917 2918 2919 !This section has been created automatically by the script Abilint (TD). 2920 !Do not modify the following lines by hand. 2921 #undef ABI_FUNC 2922 #define ABI_FUNC 'defad' 2923 !End of the abilint section 2924 2925 implicit none 2926 2927 !Arguments ------------------------------------ 2928 !scalars 2929 type(aim_dataset_type),intent(out) :: aim_dtset 2930 2931 !Local variables ------------------------------ 2932 2933 ! ********************************************************************* 2934 2935 aim_dtset%isurf=0 2936 aim_dtset%crit=0 2937 aim_dtset%irsur=0 2938 aim_dtset%foll=0 2939 aim_dtset%irho=0 2940 aim_dtset%ivol=0 2941 aim_dtset%denout=0 2942 aim_dtset%lapout=0 2943 aim_dtset%gpsurf=0 2944 aim_dtset%plden=0 2945 aim_dtset%dltyp=0 2946 2947 aim_dtset%batom=1 2948 aim_dtset%nsa=3 2949 aim_dtset%nsb=3 2950 aim_dtset%nsc=3 2951 aim_dtset%npt=100 2952 aim_dtset%nth=32 2953 aim_dtset%nph=48 2954 2955 aim_dtset%themax=pi 2956 aim_dtset%themin=zero 2957 aim_dtset%phimin=zero 2958 aim_dtset%phimax=two_pi 2959 aim_dtset%phi0=zero 2960 aim_dtset%th0=zero 2961 aim_dtset%folstp=5.d-2 2962 aim_dtset%dr0=5.d-2 2963 aim_dtset%atrad=one 2964 aim_dtset%rmin=one 2965 2966 aim_dtset%foldep(:)=zero 2967 aim_dtset%vpts(:,:)=zero 2968 aim_dtset%ngrid(:)=30 2969 aim_dtset%scal(:)=one 2970 aim_dtset%maxatd=1.d1 2971 aim_dtset%maxcpd=5.d1 2972 2973 aim_dtset%dpclim=1.d-2 2974 aim_dtset%lstep=1.d-10 2975 aim_dtset%lstep2=1.d-5 2976 aim_dtset%lgrad=1.d-12 2977 aim_dtset%lgrad2=1.d-5 2978 aim_dtset%coff1=0.98_dp 2979 aim_dtset%coff2=0.95_dp 2980 2981 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)
PARENTS
aim
CHILDREN
addout,aim_follow,cpdrv,critics,graph,initaim,integrho,integvol,plint rsurf,surf,timein
SOURCE
3016 subroutine drvaim(aim_dtset,tcpui,twalli) 3017 3018 3019 !This section has been created automatically by the script Abilint (TD). 3020 !Do not modify the following lines by hand. 3021 #undef ABI_FUNC 3022 #define ABI_FUNC 'drvaim' 3023 !End of the abilint section 3024 3025 implicit none 3026 3027 !Arguments ------------------------------------ 3028 !scalars 3029 real(dp) :: tcpui,twalli 3030 type(aim_dataset_type),intent(in) :: aim_dtset 3031 3032 !Local variables ------------------------------ 3033 !scalars 3034 integer :: iat,iatinit,inxat,ipos,iposinit 3035 integer :: me,npmax,nproc,nstep 3036 real(dp) :: dstlim,rr,ss,t1,t2,tf,wall 3037 real(dp) :: tcpu,twall,znucl_batom 3038 logical :: debold,sfour,srch,sthree,stwo 3039 !arrays 3040 real(dp) :: tsec(2) 3041 real(dp) :: grho(3),xstart(3) 3042 3043 ! ********************************************************************* 3044 3045 me=xmpi_comm_rank(xmpi_world) 3046 nproc=xmpi_comm_size(xmpi_world) 3047 3048 !These input variables might be modified during what follows, 3049 !so, they are copied outside of aim_dtset. 3050 inxat=aim_dtset%batom 3051 r0=aim_dtset%atrad 3052 h0=aim_dtset%folstp 3053 maxatdst=aim_dtset%maxatd 3054 maxcpdst=aim_dtset%maxcpd 3055 3056 dstlim=maxcpdst 3057 3058 !Flags from the old version 3059 !to be remove later 3060 deb=.false. 3061 stwo=.true. 3062 sthree=.true. 3063 sfour=.false. 3064 srch=.false. 3065 3066 npmax=aim_npmaxin 3067 3068 !Main initialisation procedure - 3069 !- it reads ABINIT density file and files 3070 !with core densities and initialises the fields for 3071 !spline interpolation 3072 3073 call initaim(aim_dtset,znucl_batom) 3074 3075 3076 !CP SEARCHING 3077 3078 if (aim_dtset%crit /= 0) then 3079 3080 call timein(tcpu,twall) 3081 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3082 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3083 & '-',ch10,'- Before searching the CP ',ch10,& 3084 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3085 3086 if (aim_dtset%crit==3) then 3087 ! old version of the driver for searching CPs (original code) 3088 call critics(aim_dtset,inxat,stwo,sthree,sfour,dstlim) 3089 else 3090 ! driver for searching CPs with Popellier algorithm 3091 call cpdrv(aim_dtset) 3092 end if 3093 3094 call timein(tcpu,twall) 3095 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3096 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3097 & '-',ch10,'- After searching the CP ',ch10,& 3098 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3099 3100 end if 3101 3102 ! 3103 !BADER SURFACE CALCULATION 3104 ! 3105 3106 if (aim_dtset%isurf==1) then 3107 ! driver for determination of the Bader surface 3108 3109 call timein(tcpu,twall) 3110 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3111 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3112 & '-',ch10,'- Before determinating the Bader surface ',ch10,& 3113 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3114 3115 call surf(aim_dtset) 3116 3117 call timein(tcpu,twall) 3118 tsec(1)=tcpu-tcpui ; tsec(2)=twall-twalli 3119 write(std_out, '(5a,f13.1,a,f13.1)' ) & 3120 & '-',ch10,'- After determinating the Bader surface ',ch10,& 3121 & '- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 3122 3123 end if 3124 3125 ! 3126 !CHARGE INTEGRATIOM 3127 ! 3128 3129 if (aim_dtset%irho==1) then 3130 call integrho(aim_dtset,znucl_batom) 3131 end if 3132 3133 ! 3134 !VOLUME INTEGRATION OF THE BADER ATOM 3135 ! 3136 3137 if (aim_dtset%ivol==1) then 3138 call integvol() 3139 end if 3140 3141 ! 3142 !ONE RADIUS OF THE BADER SURFACE 3143 ! 3144 3145 if (aim_dtset%irsur==1) then 3146 if (aim_dtset%isurf/=0) srch=.true. 3147 iat=aim_dtset%batom 3148 ss=r0 3149 call timein(t1,wall) 3150 call rsurf(aim_dtset,rr,grho,aim_dtset%th0,aim_dtset%phi0,ss,iat,npmax,srch) 3151 call timein(t2,wall) 3152 t2=t2-t1 3153 write(unts,'(2F12.8,F15.10)') aim_dtset%th0,aim_dtset%phi0,rr 3154 write(std_out,'(":RSUR ",2F12.8,2F15.10)') aim_dtset%th0,aim_dtset%phi0,rr,t2 3155 end if 3156 3157 ! 3158 !FOLLOW THE GRADIENT PATH FROM ONE POINT 3159 ! 3160 3161 if (aim_dtset%foll==1) then 3162 iatinit=aim_dtset%batom 3163 iposinit=batcell 3164 if (aim_dtset%isurf/=0) srch=.true. 3165 debold=deb 3166 xstart(:)=aim_dtset%foldep(:) 3167 call timein(t1,wall) 3168 call aim_follow(aim_dtset,xstart,npmax,srch,iatinit,iposinit,iat,ipos,nstep) 3169 call timein(t2,wall) 3170 tf=t2-t1 3171 write(std_out,'(":TIME in aim_follow:", F12.4)') tf 3172 end if 3173 3174 if (aim_dtset%plden == 1) then 3175 ! profile of the density integrated in plane xy 3176 ! belong the z-axes - not finished - cut3d better ! 3177 call plint() 3178 end if 3179 3180 if ((aim_dtset%denout > 0).or.(aim_dtset%lapout > 0)) then 3181 ! additional outputs of density and laplacian fields 3182 ! in the plane or line 3183 call addout(aim_dtset) 3184 end if 3185 3186 if (aim_dtset%gpsurf == 1) then 3187 ! script for gnuplot - simple demonstration of the 3188 ! computed surface 3189 call graph(unts,untg) 3190 end if 3191 3192 !Deallocation of global variables allocated in initaim 3193 !and declared in defs_aimfields. 3194 ABI_DEALLOCATE(dig1) 3195 ABI_DEALLOCATE(dig2) 3196 ABI_DEALLOCATE(dig3) 3197 ABI_DEALLOCATE(llg1) 3198 ABI_DEALLOCATE(llg2) 3199 ABI_DEALLOCATE(llg3) 3200 ABI_DEALLOCATE(cdig1) 3201 ABI_DEALLOCATE(cdig2) 3202 ABI_DEALLOCATE(cdig3) 3203 ABI_DEALLOCATE(ddx) 3204 ABI_DEALLOCATE(ddy) 3205 ABI_DEALLOCATE(ddz) 3206 ABI_DEALLOCATE(rrad) 3207 ABI_DEALLOCATE(crho) 3208 ABI_DEALLOCATE(sp2) 3209 ABI_DEALLOCATE(sp3) 3210 ABI_DEALLOCATE(sp4) 3211 ABI_DEALLOCATE(corlim) 3212 ABI_DEALLOCATE(dvl) 3213 ABI_DEALLOCATE(ndat) 3214 ABI_DEALLOCATE(rminl) 3215 !Deallocation of global variables allocated in initaim 3216 !and declared in defs_aimprom. 3217 ABI_DEALLOCATE(typat) 3218 ABI_DEALLOCATE(xred) 3219 ABI_DEALLOCATE(xatm) 3220 3221 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)
PARENTS
drvaim
CHILDREN
SOURCE
3246 subroutine graph(unts,untg) 3247 3248 3249 !This section has been created automatically by the script Abilint (TD). 3250 !Do not modify the following lines by hand. 3251 #undef ABI_FUNC 3252 #define ABI_FUNC 'graph' 3253 !End of the abilint section 3254 3255 implicit none 3256 3257 !Arguments ------------------------------------ 3258 !scalars 3259 integer,intent(in) :: untg,unts 3260 3261 !Local variables ------------------------------ 3262 !scalars 3263 integer :: ii,indx,jj,nphi,nth 3264 real(dp),parameter :: snull=1.d-6 3265 real(dp) :: phimax,phimin,ss,thmax,thmin 3266 !arrays 3267 real(dp) :: xorig(3) 3268 real(dp),allocatable :: phi(:),rr(:,:),th(:) 3269 3270 ! ********************************************************************* 3271 3272 rewind(unts) 3273 read(unts,*) indx, xorig(1:3) 3274 read(unts,*) nth, thmin, thmax 3275 read(unts,*) nphi, phimin, phimax 3276 ABI_ALLOCATE(th,(nth)) 3277 ABI_ALLOCATE(phi,(nphi)) 3278 ABI_ALLOCATE(rr,(nth,nphi)) 3279 do ii=1,nth 3280 do jj=1,nphi 3281 read(unts,*) th(ii),phi(jj),rr(ii,jj),ss 3282 end do 3283 end do 3284 3285 !end of reading 3286 3287 write(untg,*) 'reset' 3288 write(untg,*) 'set st d l' 3289 write(untg,*) 'set ticslevel 0' 3290 write(untg,*) 'set title ''Bader surface'' ' 3291 write(untg,*) 'splot ''-'' using ($3*sin($1)*cos($2)):($3*sin($1)*sin($2)):($3*cos($1)) notitle' 3292 do ii=1,nth 3293 do jj=1,nphi 3294 write(untg,'(2F12.8,E16.8)') th(ii),phi(jj),rr(ii,jj) 3295 end do 3296 if ((ii==nth).and.(jj==nphi)) then 3297 cycle 3298 else 3299 write(untg,*) 3300 end if 3301 end do 3302 3303 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)
PARENTS
drvaim
CHILDREN
bschg1,hdr_bcast,hdr_echo,hdr_fort_read,hdr_free,hdr_ncread,inspln lubksb,ludcmp,metric,xmpi_bcast
SOURCE
3346 subroutine initaim(aim_dtset,znucl_batom) 3347 3348 3349 !This section has been created automatically by the script Abilint (TD). 3350 !Do not modify the following lines by hand. 3351 #undef ABI_FUNC 3352 #define ABI_FUNC 'initaim' 3353 !End of the abilint section 3354 3355 implicit none 3356 3357 !Arguments ------------------------------------ 3358 !scalars 3359 type(aim_dataset_type),intent(in) :: aim_dtset 3360 3361 !Local variables ------------------------------ 3362 !scalars 3363 integer,parameter :: master=0 3364 integer :: fform0,id,ierr,ii,info,jj,kk,kod,mm,ndtmax,nn,nsa,nsb,nsc,nsym,me,nproc,npsp 3365 integer :: unth,comm 3366 #ifdef HAVE_NETCDF 3367 integer :: den_id 3368 #endif 3369 real(dp) :: ss,ucvol,znucl_batom 3370 real(dp) :: zz 3371 type(hdr_type) :: hdr 3372 !arrays 3373 integer :: ipiv(3) 3374 integer,allocatable :: symrel(:,:,:) 3375 real(dp) :: aa(3),bb(3),gmet(3,3),gprimd(3,3),rmet(3,3),yy(3,3) 3376 real(dp),allocatable :: tnons(:,:),znucl(:),zionpsp(:) 3377 real(dp),pointer :: ptc(:),ptd(:),ptf(:),ptp(:),ptsd(:) 3378 3379 ! ********************************************************************* 3380 3381 !DEBUG 3382 !write(std_out,*) ' initaim : enter ' 3383 !ENDDEBUG 3384 3385 comm = xmpi_world 3386 me=xmpi_comm_rank(comm) 3387 nproc=xmpi_comm_size(comm) 3388 3389 slc=0 ! code for follow 3390 3391 !The use of the "hdr" routines is much better for the future 3392 !maintenance of the code. Indeed, the content of the header 3393 !will continue to change from time to time, and the associated 3394 !changes should be done in one place only. 3395 3396 !Read ABINIT header ---------------------------------------------------------- 3397 if(me==master)then 3398 if (aim_iomode == IO_MODE_ETSF) then 3399 call hdr_ncread(hdr, untad, fform0) 3400 else 3401 call hdr_fort_read(hdr, untad, fform0) 3402 end if 3403 end if 3404 ABI_CHECK(fform0 /= 0, "hdr_read returned fform == 0") 3405 call hdr_bcast(hdr,master,me,comm) 3406 3407 !Echo part of the header 3408 call hdr_echo(hdr, fform0, 4, unit=std_out) 3409 call hdr_echo(hdr, fform0, 4, unit=untout) 3410 3411 natom=hdr%natom 3412 ngfft(1:3)=hdr%ngfft(:) 3413 nsym=hdr%nsym 3414 npsp=hdr%npsp 3415 ntypat=hdr%ntypat 3416 rprimd(:,:)=hdr%rprimd(:,:) 3417 3418 ABI_ALLOCATE(zionpsp,(npsp)) 3419 ABI_ALLOCATE(znucl,(ntypat)) 3420 ABI_ALLOCATE(typat,(natom)) 3421 ABI_ALLOCATE(xred,(3,natom)) 3422 ABI_ALLOCATE(symrel,(3,3,nsym)) 3423 ABI_ALLOCATE(tnons,(3,nsym)) 3424 ABI_ALLOCATE(xatm,(3,natom)) 3425 3426 symrel(:,:,:)=hdr%symrel(:,:,:) 3427 typat(:)=hdr%typat(:) 3428 tnons(:,:)=hdr%tnons(:,:) 3429 znucl(:)=hdr%znucltypat(:) 3430 zionpsp(:)=hdr%zionpsp(:) 3431 xred(:,:)=hdr%xred(:,:) 3432 3433 !This is to deallocate records of hdr 3434 call hdr_free(hdr) 3435 3436 !------------------------------------------------------------------------------- 3437 3438 ABI_ALLOCATE(dvl,(ngfft(1),ngfft(2),ngfft(3))) 3439 3440 if(me==master)then 3441 if (aim_iomode == IO_MODE_ETSF) then 3442 #ifdef HAVE_NETCDF 3443 ! netcdf array has shape [cplex, n1, n2, n3, nspden]), here we read only the total density. 3444 NCF_CHECK(nf90_inq_varid(untad, "density", den_id)) 3445 NCF_CHECK(nf90_get_var(untad, den_id, dvl, start=[1,1,1,1], count=[1, ngfft(1), ngfft(2), ngfft(3), 1])) 3446 #endif 3447 else 3448 read(untad,iostat=nn) dvl(1:ngfft(1),1:ngfft(2),1:ngfft(3)) 3449 ABI_CHECK(nn==0,"error of reading !") 3450 end if 3451 end if 3452 call xmpi_bcast(dvl, master, comm, ierr) 3453 3454 write(std_out,*)ch10,' initaim : the valence density has been read' ,ch10 3455 3456 !INITIALISATION OF SOME IMPORTANT FIELDS 3457 3458 !Only interpolation is computed (inside vgh_rho) in reduced 3459 !coordinates. In all other routines the cart. coordinates (CC) are used. 3460 3461 !transformation of the atom positions to CC 3462 do ii=1,natom 3463 xatm(:,ii)=xred(:,ii) 3464 call bschg1(xatm(:,ii),1) 3465 end do 3466 3467 !Generation of the neighbouring cells + transf to CC 3468 nn=0 3469 nsa=aim_dtset%nsa ; nsb=aim_dtset%nsb ; nsc=aim_dtset%nsc 3470 do ii=-nsa,nsa 3471 do jj=-nsb,nsb 3472 do kk=-nsc,nsc 3473 nn=nn+1 3474 atp(1,nn)=ii*1._dp 3475 atp(2,nn)=jj*1._dp 3476 atp(3,nn)=kk*1._dp 3477 call bschg1(atp(:,nn),1) 3478 end do 3479 end do 3480 end do 3481 nnpos=nn 3482 3483 !DEBUG 3484 !write(std_out,*)' initaim : nnpos=',nnpos 3485 !ENDDEBUG 3486 3487 batcell=nsa*(2*nsb+1)*(2*nsc+1)+(2*nsc+1)*nsb+nsc+1 3488 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 3489 maxatdst=min(maxatdst, nsa*sqrt(rmet(1,1)), nsb*sqrt(rmet(2,2)), nsc*sqrt(rmet(3,3)) ) 3490 if (maxcpdst > maxatdst) maxcpdst=0.75*maxatdst 3491 3492 3493 !RPRIM ITS INVERSE AND TRANSPOSE 3494 3495 do ii=1,3 3496 do jj=1,3 3497 yy(ii,jj)=rprimd(ii,jj) 3498 end do 3499 end do 3500 call ludcmp(yy,3,3,ipiv,id,info) 3501 ABI_CHECK(info==0,'Error inverting rprimd') 3502 3503 do ii=1,3 3504 do jj=1,3 3505 ivrprim(ii,jj)=0._dp 3506 end do 3507 ivrprim(ii,ii)=1._dp 3508 end do 3509 do ii=1,3 3510 call lubksb(yy,3,3,ipiv,ivrprim(:,ii)) 3511 end do 3512 do ii=1,3 3513 do jj=1,3 3514 trivrp(ii,jj)=ivrprim(jj,ii) 3515 end do 3516 end do 3517 3518 write(std_out,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 3519 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 3520 write(untout,'(" INVERSE OF RPRIMD: ",/,3F16.8,/,3F16.8,/,3F16.8,/)') & 3521 & ((ivrprim(ii,jj), jj=1,3), ii=1,3) 3522 3523 write(std_out,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 3524 write(std_out,*) "=======================================" 3525 do ii=1,natom 3526 jj=typat(ii) 3527 write(std_out,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 3528 end do 3529 write(untout,*) "ATOMS (index,at.number,Zionic,position(xcart.))" 3530 write(untout,*) "=======================================" 3531 do ii=1,natom 3532 jj=typat(ii) 3533 write(untout,'(I4,2F10.6,3F16.8)') ii, znucl(jj), zionpsp(jj), (xatm(kk,ii),kk=1,3) 3534 end do 3535 3536 !STEPS IN REAL SPACE GRID (REDUCED) 3537 do ii=1,3 3538 dix(ii)=1._dp/ngfft(ii) 3539 end do 3540 3541 !READING OF THE CORE DENSITY 3542 write(std_out,*)ch10,' initaim : will read the core densities' ,ch10 3543 3544 ABI_ALLOCATE(ndat,(ntypat)) 3545 ABI_ALLOCATE(rminl,(natom)) 3546 ndtmax=0 3547 if(me==master)then 3548 do ii=1,ntypat 3549 unth=unt+ii 3550 ! DEBUG 3551 ! write(std_out,*)' read from unit ',unth 3552 ! call flush(std_out) 3553 ! stop 3554 ! ENDDEBUG 3555 read(unth,*) ndat(ii),ss 3556 if (ndat(ii)>ndtmax) ndtmax=ndat(ii) 3557 end do 3558 end if 3559 call xmpi_bcast(ndat,master,comm,ierr) 3560 call xmpi_bcast(ndtmax,master,comm,ierr) 3561 call xmpi_bcast(ss,master,comm,ierr) 3562 3563 !FIELDS FOR STORING CORE DENSITY 3564 3565 ABI_ALLOCATE(rrad,(ndtmax,ntypat)) 3566 ABI_ALLOCATE(crho,(ndtmax,ntypat)) 3567 ABI_ALLOCATE(sp2,(ndtmax,ntypat)) 3568 ABI_ALLOCATE(sp3,(ndtmax,ntypat)) 3569 ABI_ALLOCATE(sp4,(ndtmax,ntypat)) 3570 ABI_ALLOCATE(corlim,(ntypat)) 3571 3572 sp2(:,:)=zero 3573 sp3(:,:)=zero 3574 sp4(:,:)=zero 3575 3576 !Reading of the core densities 3577 corlim(:)=0 3578 kod=0 3579 if(me==master)then 3580 do ii=1,ntypat 3581 unth=unt+ii 3582 do jj=1,ndat(ii) 3583 read(unth,*) rrad(jj,ii),crho(jj,ii),sp2(jj,ii),sp3(jj,ii) 3584 ! this is the integral of the core charge read in 3585 crho(jj,ii) = crho(jj,ii)/4._dp/pi 3586 if ((crho(jj,ii) < aim_rhocormin) .and. (corlim(ii)==0)) corlim(ii)=jj 3587 sp2(jj,ii)=sp2(jj,ii)/4._dp/pi 3588 sp3(jj,ii)=sp3(jj,ii)/4._dp/pi ! ATENTION!!! in sp3 is just second derivation 3589 end do 3590 do jj=1,ndat(ii)-1 3591 sp4(jj,ii)=(sp3(jj+1,ii)-sp3(jj,ii))/(6._dp*(rrad(jj+1,ii)-rrad(jj,ii))) 3592 end do 3593 ! 3594 zz = crho(1,ii) * rrad(1,ii)**2 * (rrad(2,ii)-rrad(1,ii)) 3595 do jj=2,ndat(ii)-1 3596 zz = zz + crho(jj,ii) * rrad(jj,ii)**2 * (rrad(jj+1,ii)-rrad(jj-1,ii)) 3597 end do 3598 zz = zz * half * 4._dp * pi 3599 if (corlim(ii)==0) corlim(ii)=ndat(ii) 3600 3601 ! add check on zion wrt FHI .fc file 3602 ! compare zion to zionpsp(typat(aim_dtset%batom)) 3603 if (abs(znucl(ii) - zz - zionpsp(ii)) > 1.e-1_dp) then 3604 write (std_out,*) 'error: your core charge ', zz, ' does not correspond to the correct number' 3605 write (std_out,*) ' of valence electrons', zionpsp(ii), ' and the nuclear charge ', znucl(ii) 3606 write (std_out,*) ' You have probably used a pseudopotential which has more valence electrons than the' 3607 write (std_out,*) ' original FHI ones. ACTION: make a .fc file with the correct core charge' 3608 stop 3609 end if 3610 3611 end do 3612 end if 3613 call xmpi_bcast(rrad,master,comm,ierr) 3614 call xmpi_bcast(crho,master,comm,ierr) 3615 call xmpi_bcast(sp2,master,comm,ierr) 3616 call xmpi_bcast(sp3,master,comm,ierr) 3617 call xmpi_bcast(sp4,master,comm,ierr) 3618 call xmpi_bcast(corlim,master,comm,ierr) 3619 3620 write(std_out,*)ch10,' initaim : the core densities have been read' ,ch10 3621 3622 3623 !CORRECTION OF THE CORE DENSITY NORMALISATION 3624 crho(:,:)=1.0003*crho(:,:) 3625 sp2(:,:)=1.0003*sp2(:,:) 3626 sp3(:,:)=1.0003*sp3(:,:) 3627 sp4(:,:)=1.0003*sp4(:,:) 3628 3629 !FIELDS FOR INTERPOLATIONS OF THE VALENCE DENSITY 3630 3631 ABI_ALLOCATE(dig1,(ngfft(1))) 3632 ABI_ALLOCATE(dig2,(ngfft(2))) 3633 ABI_ALLOCATE(dig3,(ngfft(3))) 3634 ABI_ALLOCATE(llg1,(ngfft(1))) 3635 ABI_ALLOCATE(llg2,(ngfft(2))) 3636 ABI_ALLOCATE(llg3,(ngfft(3))) 3637 ABI_ALLOCATE(cdig1,(ngfft(1)-1)) 3638 ABI_ALLOCATE(cdig2,(ngfft(2)-1)) 3639 ABI_ALLOCATE(cdig3,(ngfft(3)-1)) 3640 ABI_ALLOCATE(ddx,(ngfft(1),ngfft(2),ngfft(3))) 3641 ABI_ALLOCATE(ddy,(ngfft(1),ngfft(2),ngfft(3))) 3642 ABI_ALLOCATE(ddz,(ngfft(1),ngfft(2),ngfft(3))) 3643 3644 !DECOMPOSITION OF THE MATRIX FOR THE DETERMINATION OF COEFFICIENTS 3645 !FOR CUBIC SPLINE INTERPOLATION (using the periodic boundary conditions) 3646 3647 !MAIN DIAGONAL (aa) AND SECONDARY DIAGONAL (bb) MATRIX ELEMENTS 3648 3649 nmax=ngfft(1) 3650 do ii=2,3 3651 if (ngfft(ii) > nmax) nmax=ngfft(ii) 3652 end do 3653 nullify(ptf,ptsd) 3654 nullify(ptd,ptc,ptp) 3655 aa(:)=2.0*dix(:)**2/3.0 3656 bb(:)=dix(:)**2/6.0 3657 3658 do ii=1,3 3659 if(ii==1) then 3660 ptd=>dig1;ptc=>cdig1;ptp=>llg1 3661 elseif (ii==2) then 3662 ptd=>dig2;ptc=>cdig2;ptp=>llg2 3663 else 3664 ptd=>dig3;ptc=>cdig3;ptp=>llg3 3665 end if 3666 ptd(1)=sqrt(aa(ii)) 3667 ptc(1)=bb(ii)/ptd(1) 3668 ptp(1)=ptc(1) 3669 do jj=2,ngfft(ii)-1 3670 ptd(jj)=aa(ii)-ptc(jj-1)**2 3671 if(ptd(jj)<zero) then 3672 MSG_ERROR('Matrix is not positive definite !') 3673 end if 3674 ptd(jj)=sqrt(ptd(jj)) 3675 if (jj==ngfft(ii)-1) then 3676 ptc(jj)=(bb(ii)-ptp(jj-1)*ptc(jj-1))/ptd(jj) 3677 ptp(jj)=ptc(jj) 3678 exit 3679 end if 3680 ptc(jj)=bb(ii)/ptd(jj) 3681 ptp(jj)=-ptp(jj-1)*ptc(jj-1)/ptd(jj) 3682 end do 3683 ss=0._dp 3684 do jj=1,ngfft(ii)-1 3685 ss=ss+ptp(jj)**2 3686 end do 3687 ss=aa(ii)-ss 3688 if(ss<zero) then 3689 MSG_ERROR('Matrix is not positive definite !') 3690 end if 3691 ptd(ngfft(ii))=sqrt(ss) 3692 ptp(ngfft(ii))=ptd(ngfft(ii)) 3693 3694 3695 ! INITIALISATION OF THE SECOND DERIVATIVE FIELDS 3696 3697 nn=ii+1 3698 if (nn>3) nn=nn-3 3699 mm=ii+2 3700 if (mm>3) mm=mm-3 3701 do jj=1,ngfft(nn) 3702 do kk=1,ngfft(mm) 3703 ! The calcul of the second derivations on the grid 3704 call inspln(ii,jj,kk) 3705 end do 3706 end do 3707 nullify(ptd,ptc,ptp) 3708 end do 3709 nullify(ptd,ptc,ptp) 3710 3711 znucl_batom=znucl(typat(aim_dtset%batom)) 3712 3713 ABI_DEALLOCATE(znucl) 3714 ABI_DEALLOCATE(zionpsp) 3715 ABI_DEALLOCATE(symrel) 3716 ABI_DEALLOCATE(tnons) 3717 3718 !the pointers are obsolete - to remove later 3719 3720 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)
PARENTS
aim
CHILDREN
SOURCE
3747 subroutine inpar(instr,lenstr) 3748 3749 3750 !This section has been created automatically by the script Abilint (TD). 3751 !Do not modify the following lines by hand. 3752 #undef ABI_FUNC 3753 #define ABI_FUNC 'inpar' 3754 !End of the abilint section 3755 3756 implicit none 3757 3758 !Arguments ------------------------------------ 3759 !scalars 3760 integer,intent(out) :: lenstr 3761 character(len=*),intent(out) :: instr 3762 3763 !Local variables ------------------------------ 3764 character(len=1),parameter :: space=' ' 3765 character(len=26),parameter :: uplett='ABCDEFGHIJKLMNOPQRSTUVWXYZ', lolett='abcdefghijklmnopqrstuvwxyz' 3766 !scalars 3767 integer,parameter :: nline=100 3768 integer :: ii,inxh,inxl,ios,jj,kk,ll 3769 character(len=fnlen) :: line 3770 3771 ! ********************************************************************* 3772 3773 lenstr=0 3774 3775 do ii=1,26 3776 inxh=index(lolett,uplett(ii:ii)) 3777 if (inxh > 0) then 3778 write(std_out,*) 'ERROR The ', uplett(ii:ii) ,' is considered come lowcase !' 3779 MSG_ERROR("Aborting now") 3780 end if 3781 end do 3782 rewind(unt0) 3783 do ii=1,nline 3784 read(unt0,'(A)',iostat=ios) line(1:fnlen) 3785 if (ios/=0) exit 3786 inxh=index(line,'#') 3787 if (inxh == 1) then 3788 cycle 3789 elseif (inxh > 0) then 3790 inxl=inxh-1 3791 line(inxh:inxh)=space 3792 else 3793 inxl=len_trim(line) 3794 if (inxl==0) cycle 3795 end if 3796 inxh=index(line(1:inxl),char(9)) 3797 if (inxh/=0) line(inxh:inxh)=space 3798 do ll=1,inxl 3799 if (iachar(line(ll:ll)) < 32) line(ll:ll)=space 3800 end do 3801 inxh=index(line(1:inxl),'- ') 3802 if (inxh/=0) then 3803 write(std_out,*) 'ERROR sign minus with white space in input file' 3804 MSG_ERROR("Aborting now") 3805 end if 3806 line(1:inxl)=adjustl(line(1:inxl)) 3807 inxl=len_trim(line(1:inxl))+1 3808 jj=2;kk=0 3809 line(1:inxl)=adjustl(line(1:inxl)) 3810 kk=len_trim(line(1:inxl))+1 3811 do ll=1,inxl 3812 inxh=index(line(jj:kk),space) 3813 if ((inxh==0).or.((jj+inxh-1)==kk)) exit 3814 line(inxh+jj:kk)=adjustl(line(inxh+jj:kk)) 3815 kk=len_trim(line(1:inxl)) 3816 if (kk == inxl) then 3817 exit 3818 end if 3819 jj=jj+inxh 3820 end do 3821 inxl=len_trim(line(1:inxl))+1 3822 do ll=1,inxl-1 3823 inxh=index(lolett,line(ll:ll)) 3824 if (inxh/=0) line(ll:ll)=uplett(inxh:inxh) 3825 end do 3826 if ((lenstr+inxl) > strlen ) then 3827 write(std_out,*) 'ERROR Too large input !' 3828 MSG_ERROR("Aborting now") 3829 else 3830 instr(lenstr+1:lenstr+inxl)=line(1:inxl) 3831 lenstr=lenstr+inxl 3832 end if 3833 end do 3834 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)
PARENTS
initaim
CHILDREN
SOURCE
3868 subroutine inspln(idir,snn,tnn) 3869 3870 3871 !This section has been created automatically by the script Abilint (TD). 3872 !Do not modify the following lines by hand. 3873 #undef ABI_FUNC 3874 #define ABI_FUNC 'inspln' 3875 !End of the abilint section 3876 3877 implicit none 3878 3879 !Arguments ------------------------------------ 3880 !scalars 3881 integer,intent(in) :: idir,snn,tnn 3882 3883 !Local variables------------------------------- 3884 !scalars 3885 integer :: dim,ii 3886 real(dp) :: ss 3887 !arrays 3888 real(dp) :: rsid(ngfft(idir)),valf(ngfft(idir)) 3889 real(dp),pointer :: ptc(:),ptd(:),ptp(:) 3890 3891 ! ************************************************************************* 3892 3893 !POINTER INITIALIZATION 3894 3895 if (idir==1) then 3896 valf(:)=dvl(:,snn,tnn) 3897 elseif (idir==2) then 3898 valf(:)=dvl(tnn,:,snn) 3899 else 3900 valf(:)=dvl(snn,tnn,:) 3901 end if 3902 3903 nullify(ptd,ptc,ptp) 3904 if(idir==1) then 3905 ptd=>dig1;ptc=>cdig1;ptp=>llg1 3906 elseif (idir==2) then 3907 ptd=>dig2;ptc=>cdig2;ptp=>llg2 3908 else 3909 ptd=>dig3;ptc=>cdig3;ptp=>llg3 3910 end if 3911 3912 dim=ngfft(idir) 3913 3914 !FIRST CYCLE OF RECURRENCE 3915 3916 rsid(1)=valf(2)+valf(dim)-2.*valf(1) 3917 rsid(1)=rsid(1)/ptd(1) 3918 do ii=2,dim-1 3919 rsid(ii)=valf(ii+1)+valf(ii-1)-2.*valf(ii) 3920 rsid(ii)=(rsid(ii)-ptc(ii-1)*rsid(ii-1))/ptd(ii) 3921 end do 3922 ss=0._dp 3923 do ii=1,dim-1 3924 ss=ss+rsid(ii)*ptp(ii) 3925 end do 3926 rsid(dim)=valf(1)+valf(dim-1)-2.*valf(dim) 3927 rsid(dim)=(rsid(dim)-ss)/ptd(dim) 3928 3929 !SECOND CYCLE WITH TRANSPOSED MATRIX 3930 3931 rsid(dim)=rsid(dim)/ptd(dim) 3932 rsid(dim-1)=(rsid(dim-1)-ptc(dim-1)*rsid(dim))/ptd(dim-1) 3933 do ii=dim-2,1,-1 3934 rsid(ii)=(rsid(ii)-ptc(ii)*rsid(ii+1)-ptp(ii)*rsid(dim))/ptd(ii) 3935 end do 3936 3937 if (idir==1) then 3938 ddx(:,snn,tnn)=rsid(:) 3939 elseif (idir==2) then 3940 ddy(tnn,:,snn)=rsid(:) 3941 else 3942 ddz(snn,tnn,:)=rsid(:) 3943 end if 3944 3945 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)
PARENTS
drvaim
CHILDREN
bschg1,coeffs_gausslegint,spline,vgh_rho
SOURCE
3979 subroutine integrho(aim_dtset,znucl_batom) 3980 3981 3982 !This section has been created automatically by the script Abilint (TD). 3983 !Do not modify the following lines by hand. 3984 #undef ABI_FUNC 3985 #define ABI_FUNC 'integrho' 3986 !End of the abilint section 3987 3988 implicit none 3989 3990 !Arguments ------------------------------------ 3991 !scalars 3992 type(aim_dataset_type),intent(in) :: aim_dtset 3993 3994 !Local variables ------------------------------ 3995 !scalars 3996 integer :: batom,chs,iat,ii,inx,inxf,ipos,jj,kk,ll,nn,nph,nth 3997 real(dp) :: chg,chgint,cintr,ct1,ct2,lder,nsphe,phimax,phimin,rder 3998 real(dp) :: rsmax,rsmin,ss,stp,themax,themin,uu 3999 real(dp) :: znucl_batom,zz 4000 logical :: gaus,weit 4001 !arrays 4002 real(dp) :: grho(3),hrho(3,3),shift(3),unvec(3),vv(3) 4003 real(dp),allocatable :: ncrho(:),nsp2(:),nsp3(:),nsp4(:),rdint(:,:),rr(:) 4004 real(dp),allocatable :: vdd(:),vrho(:),wgrs(:,:),work(:) 4005 4006 ! ********************************************************************* 4007 4008 gaus=.true. 4009 weit=.true. 4010 4011 write(std_out,*) 'npt = ',aim_dtset%npt 4012 4013 rewind(unts) 4014 read(unts,*) batom,shift ! Warning : batom is read, instead of coming from aim_dtset 4015 read(unts,*) nth,themin,themax ! Warning : these numbers are read, instead of coming from aim_dtset 4016 read(unts,*) nph,phimin,phimax ! Warning : these numbers are read, instead of coming from aim_dtset 4017 4018 write(std_out,*) 'NTH NPH ',nth,nph 4019 4020 ABI_ALLOCATE(wgrs,(nth,nph)) 4021 ABI_ALLOCATE(rdint,(nth,nph)) 4022 4023 do ii=1,nth 4024 do jj=1,nph 4025 if (weit) then 4026 read(unts,*) th(ii),ph(jj),rs(ii,jj),wgrs(ii,jj) 4027 else 4028 read(unts,*) th(ii),ph(jj),rs(ii,jj) 4029 end if 4030 end do 4031 end do 4032 read(unts,*) rsmin,rsmax 4033 4034 4035 if (gaus) then 4036 ct1=cos(themin) 4037 ct2=cos(themax) 4038 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 4039 call coeffs_gausslegint(phimin,phimax,ph,wph,nph) 4040 end if 4041 4042 do ii=1,nth 4043 do jj=1,nph 4044 if (.not.weit) then 4045 if (gaus) then 4046 wgrs(ii,jj)=wcth(ii)*wph(jj) 4047 else 4048 wgrs(ii,jj)=1._dp 4049 end if 4050 end if 4051 end do 4052 end do 4053 4054 4055 do ii=1,nth 4056 do jj=1,nph 4057 if (rs(ii,jj) < rsmin) rsmin=rs(ii,jj) 4058 end do 4059 end do 4060 4061 4062 !INTEGRATION OF THE CORE DENSITY 4063 4064 nn=typat(batom) 4065 kk=ndat(nn) 4066 4067 4068 !spherical integration of the core density in the sphere 4069 !of the minimal Bader radius 4070 4071 !COEF. FOR SPHERICAL INTEGRATION 4072 4073 ABI_ALLOCATE(nsp2,(kk)) 4074 ABI_ALLOCATE(nsp3,(kk)) 4075 ABI_ALLOCATE(nsp4,(kk)) 4076 ABI_ALLOCATE(ncrho,(kk)) 4077 4078 do ii=1,kk 4079 ncrho(ii)=crho(ii,nn)*4._dp*pi*rrad(ii,nn)*rrad(ii,nn) 4080 nsp3(ii)=4._dp*pi*(2._dp*crho(ii,nn)+2._dp*rrad(ii,nn)*sp2(ii,nn)+& 4081 & rrad(ii,nn)*rrad(ii,nn)*sp3(ii,nn)) 4082 end do 4083 4084 if (rsmin < rrad(ndat(nn),nn)) then ! search index 4085 inx=0 4086 if (rsmin < rrad(1,nn)) then 4087 MSG_ERROR('absurd') 4088 elseif (rsmin > rrad(ndat(nn),nn)) then 4089 inx=ndat(nn) 4090 else 4091 do while (rsmin >= rrad(inx+1,nn)) 4092 inx=inx+1 4093 end do 4094 end if 4095 else 4096 inx=ndat(nn) 4097 end if 4098 4099 cintr=4._dp/3._dp*pi*rrad(1,nn)**3*crho(1,nn) 4100 4101 !spline integration 4102 4103 do ii=1,inx-1 4104 uu=rrad(ii+1,nn)-rrad(ii,nn) 4105 cintr=cintr+(ncrho(ii)+ncrho(ii+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(ii)+nsp3(ii+1)) 4106 end do 4107 if (inx/=ndat(nn)) then 4108 uu=rsmin-rrad(inx,nn) 4109 zz=rrad(inx+1,nn)-rsmin 4110 ss=rrad(inx+1,nn)-rrad(inx,nn) 4111 cintr=cintr+ncrho(inx)/2._dp*(ss-zz*zz/ss)+ncrho(inx+1)/2._dp*uu*uu/ss+& 4112 nsp3(inx)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 4113 nsp3(inx+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 4114 end if 4115 4116 4117 !INTEGRATION OF THE REST OF THE CORE DENSITY 4118 !(for gauss quadrature) 4119 !For the Gauss quadrature it is added 4120 !to the radial integrated valence density 4121 4122 rdint(:,:)=0._dp 4123 nsphe=0._dp 4124 do ii=1,nth 4125 do jj=1,nph 4126 if (inx==ndat(nn)) cycle 4127 inxf=inx 4128 if (rs(ii,jj) < rsmin) then 4129 write(std_out,*) rs(ii,jj),rsmin 4130 MSG_ERROR('in surface') 4131 elseif (rs(ii,jj) > rrad(ndat(nn),nn)) then 4132 inxf=ndat(nn) 4133 else 4134 do while (rs(ii,jj) >= rrad(inxf+1,nn)) 4135 inxf=inxf+1 4136 end do 4137 end if 4138 4139 if (inxf==inx) then 4140 uu=rrad(inx+1,nn)-rs(ii,jj) 4141 zz=rrad(inx+1,nn)-rsmin 4142 ss=rrad(inx+1,nn)-rrad(inx,nn) 4143 4144 rdint(ii,jj)=(ncrho(inx)/2._dp/ss-nsp3(inx)/1.2d1*ss)*(zz*zz-uu*uu)+& 4145 nsp3(inx)/2.4d1/ss*(zz**4-uu**4) 4146 uu=rs(ii,jj)-rrad(inx,nn) 4147 zz=rsmin-rrad(inx,nn) 4148 rdint(ii,jj)=rdint(ii,jj)+(uu*uu-zz*zz)*(ncrho(inx+1)/2._dp/ss-nsp3(inx+1)/1.2d1*ss)+& 4149 nsp3(inx+1)/2.4d1/ss*(uu**4-zz**4) 4150 else 4151 uu=rrad(inx+1,nn)-rsmin 4152 zz=rsmin-rrad(inx,nn) 4153 4154 rdint(ii,jj)=ncrho(inx)/2._dp/ss*uu*uu+ncrho(inx+1)/2._dp*(ss-zz*zz/ss)+& 4155 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) 4156 if (inxf > inx+1) then 4157 do kk=inx+1,inxf-1 4158 uu=rrad(kk+1,nn)-rrad(kk,nn) 4159 rdint(ii,jj)=rdint(ii,jj)+(ncrho(kk)+ncrho(kk+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(kk)+nsp3(kk+1)) 4160 end do 4161 end if 4162 4163 if (inxf/=ndat(nn)) then 4164 uu=rs(ii,jj)-rrad(inxf,nn) 4165 zz=rrad(inxf+1,nn)-rs(ii,jj) 4166 ss=rrad(inxf+1,nn)-rrad(inxf,nn) 4167 rdint(ii,jj)=rdint(ii,jj)+ncrho(inxf)/2._dp*(ss-zz*zz/ss)+ncrho(inxf+1)/2._dp*uu*uu/ss+& 4168 nsp3(inxf)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 4169 nsp3(inxf+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 4170 end if 4171 end if 4172 rdint(ii,jj)=rdint(ii,jj)/4._dp/pi 4173 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 4174 end do 4175 end do 4176 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 4177 4178 write(untout,*) 4179 write(untout,*) "CHARGE INTEGRATION" 4180 write(untout,*) "==================" 4181 write(untout,'(" Core density contribution: ",/,/," ",F16.8)') cintr+nsphe 4182 4183 write(std_out,*) ':INTECOR ', cintr+nsphe 4184 4185 ABI_DEALLOCATE(ncrho) 4186 ABI_DEALLOCATE(nsp2) 4187 ABI_DEALLOCATE(nsp3) 4188 ABI_DEALLOCATE(nsp4) 4189 4190 !INTEGRATION OF THE VALENCE DENSITY 4191 4192 ABI_ALLOCATE(rr,(aim_dtset%npt+1)) 4193 ABI_ALLOCATE(vrho,(aim_dtset%npt+1)) 4194 ABI_ALLOCATE(vdd,(aim_dtset%npt+1)) 4195 4196 !in the case of the only irho appelation 4197 4198 nn=0 4199 do ii=-3,3 4200 do jj=-3,3 4201 do kk=-3,3 4202 nn=nn+1 4203 atp(1,nn)=ii*1._dp 4204 atp(2,nn)=jj*1._dp 4205 atp(3,nn)=kk*1._dp 4206 call bschg1(atp(:,nn),1) 4207 if ((ii==0).and.(jj==0).and.(kk==0)) ipos=nn 4208 end do 4209 end do 4210 end do 4211 nnpos=nn 4212 iat=batom 4213 4214 !XG020629 There is a problem with this routine 4215 !(or vgh_rho), when one uses the PGI compiler : 4216 !The following line is needed, otherwise, iat and ipos 4217 !are set to 0 inside vgh_now. Why ???? 4218 write(std_out,*)' integrho : iat,ipos=',iat,ipos 4219 ! 4220 4221 nsphe=0._dp 4222 ABI_ALLOCATE(work,(aim_dtset%npt+1)) 4223 do ii=1,nth 4224 do jj=1,nph 4225 4226 stp=rs(ii,jj)/aim_dtset%npt 4227 unvec(1)=sin(th(ii))*cos(ph(jj)) 4228 unvec(2)=sin(th(ii))*sin(ph(jj)) 4229 unvec(3)=cos(th(ii)) 4230 do kk=0,aim_dtset%npt 4231 rr(kk+1)=kk*stp 4232 vv(:)=xatm(:,batom)+kk*stp*unvec(:) 4233 chs=-2 4234 call vgh_rho(vv,chg,grho,hrho,uu,iat,ipos,chs) 4235 vrho(kk+1)=chg*rr(kk+1)*rr(kk+1) 4236 if (kk==aim_dtset%npt) then 4237 rder=0._dp 4238 do ll=1,3 4239 rder=rder+grho(ll)*unvec(ll) 4240 end do 4241 rder=rder*rr(kk+1)*rr(kk+1)+2._dp*rr(kk+1)*chg 4242 end if 4243 end do 4244 lder=0._dp 4245 kk=aim_dtset%npt+1 4246 call spline(rr,vrho,kk,lder,rder,vdd) 4247 4248 ! INTEGRATION 4249 4250 do kk=1,aim_dtset%npt 4251 rdint(ii,jj)=rdint(ii,jj)+stp/2._dp*(vrho(kk)+vrho(kk+1))& 4252 & -stp*stp*stp/24._dp*(vdd(kk)+vdd(kk+1)) 4253 end do 4254 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 4255 end do 4256 end do 4257 ABI_DEALLOCATE(work) 4258 4259 if (gaus.or.weit) then 4260 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 4261 else 4262 nsphe=nsphe/(nth*nph)*2.0*two_pi 4263 end if 4264 chgint=cintr+nsphe 4265 4266 write(untout,'(/," Different density contributions: Core (only spherical part) and the rest ",/,/," ",2F16.8)') & 4267 & cintr, nsphe 4268 write(untout,'(/,a,i4,a,f14.8)') ' For atom number ',batom,', the number of electrons in the Bader volume is ',chgint 4269 write(untout,'(a,f15.7,a,f17.8)') ' The nuclear charge is',znucl_batom,', so that the Bader charge is ',znucl_batom-chgint 4270 write(untout,*) 4271 write(std_out,*) ':INTEPAR ', cintr, nsphe 4272 write(std_out,*) ':RHOTOT ',batom,chgint 4273 4274 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)
PARENTS
drvaim
CHILDREN
coeffs_gausslegint
SOURCE
4304 subroutine integvol() 4305 4306 4307 !This section has been created automatically by the script Abilint (TD). 4308 !Do not modify the following lines by hand. 4309 #undef ABI_FUNC 4310 #define ABI_FUNC 'integvol' 4311 !End of the abilint section 4312 4313 implicit none 4314 4315 !Arguments ------------------------------------ 4316 4317 !Local variables ------------------------------ 4318 !scalars 4319 integer :: batom,ii,jj,nph,nth 4320 real(dp) :: chgint,ct1,ct2,nsphe,phimax,phimin 4321 real(dp) :: rsmax,rsmin,themax,themin 4322 logical :: gaus,weit 4323 !arrays 4324 real(dp) :: shift(3) 4325 real(dp),allocatable :: rdint(:,:) 4326 real(dp),allocatable :: wgrs(:,:) 4327 4328 ! ********************************************************************* 4329 4330 tpi=two_pi 4331 gaus=.true. 4332 weit=.true. 4333 4334 4335 rewind(unts) 4336 read(unts,*) batom,shift 4337 read(unts,*) nth,themin,themax 4338 read(unts,*) nph,phimin,phimax 4339 4340 write(std_out,*) 'NTH NPH ',nth,nph 4341 4342 ABI_ALLOCATE(wgrs,(nth,nph)) 4343 ABI_ALLOCATE(rdint,(nth,nph)) 4344 4345 do ii=1,nth 4346 do jj=1,nph 4347 if (weit) then 4348 read(unts,*) th(ii),ph(jj),rs(ii,jj),wgrs(ii,jj) 4349 else 4350 read(unts,*) th(ii),ph(jj),rs(ii,jj) 4351 end if 4352 end do 4353 end do 4354 read(unts,*) rsmin,rsmax 4355 4356 4357 if (gaus) then 4358 ct1=cos(themin) 4359 ct2=cos(themax) 4360 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 4361 call coeffs_gausslegint(phimin,phimax,ph,wph,nph) 4362 end if 4363 4364 do ii=1,nth 4365 do jj=1,nph 4366 if (.not.weit) then 4367 if (gaus) then 4368 wgrs(ii,jj)=wcth(ii)*wph(jj) 4369 else 4370 wgrs(ii,jj)=1._dp 4371 end if 4372 end if 4373 end do 4374 end do 4375 4376 nsphe=0._dp 4377 do ii=1,nth 4378 do jj=1,nph 4379 nsphe=nsphe+rs(ii,jj)**3/3._dp*wgrs(ii,jj) 4380 end do 4381 end do 4382 if (gaus.or.weit) then 4383 nsphe=nsphe*(pi/(themin-themax))*(tpi/(phimax-phimin)) 4384 else 4385 nsphe=nsphe/(nth*nph)*2.0*tpi 4386 end if 4387 chgint=nsphe 4388 4389 write(std_out,*) ':VOLTOT ',batom,chgint 4390 write(untout,'("Volume of the Bader atom: ", I6, F16.8)') batom,chgint 4391 4392 end subroutine integvol
m_bader/mprod [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
mprod
FUNCTION
Matrix multiplication cc=aa*bb
PARENTS
vnorm
CHILDREN
mprod
SOURCE
6074 subroutine mprod(aa,bb,cc) 6075 6076 6077 !This section has been created automatically by the script Abilint (TD). 6078 !Do not modify the following lines by hand. 6079 #undef ABI_FUNC 6080 #define ABI_FUNC 'mprod' 6081 !End of the abilint section 6082 6083 implicit none 6084 6085 !Arguments ------------------------------------ 6086 !arrays 6087 real(dp),intent(in) :: aa(3,3),bb(3,3) 6088 real(dp),intent(out) :: cc(3,3) 6089 6090 !Local variables------------------------------- 6091 !scalars 6092 integer :: ii,jj,kk 6093 6094 ! ************************************************************************* 6095 6096 do ii=1,3 6097 do jj=1,3 6098 cc(ii,jj)=0._dp 6099 do kk=1,3 6100 cc(ii,jj)=cc(ii,jj)+aa(ii,kk)*bb(kk,jj) 6101 end do 6102 end do 6103 end do 6104 6105 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)
PARENTS
aim_follow
CHILDREN
vgh_rho
SOURCE
4428 subroutine onestep(vv,chg,grho,hh,np,npmax,deltar) 4429 4430 4431 !This section has been created automatically by the script Abilint (TD). 4432 !Do not modify the following lines by hand. 4433 #undef ABI_FUNC 4434 #define ABI_FUNC 'onestep' 4435 !End of the abilint section 4436 4437 implicit none 4438 4439 !Arguments ------------------------------------ 4440 !scalars 4441 integer,intent(in) :: npmax 4442 integer,intent(out) :: np 4443 real(dp),intent(in) :: hh 4444 real(dp),intent(out) :: chg,deltar 4445 !arrays 4446 real(dp),intent(inout) :: vv(3) 4447 real(dp),intent(out) :: grho(3) 4448 4449 !Local variables ------------------------------ 4450 !scalars 4451 integer :: iat,ii,ipos,jj 4452 real(dp) :: dt,rr 4453 !arrays 4454 real(dp) :: hrho(3,3),pom(3),vinter(3,200),vk(3),vkold(3) 4455 4456 !************************************************************************ 4457 dt=hh 4458 np=1 4459 deltar=1._dp 4460 vk(1:3)=vv(1:3) 4461 4462 4463 do while((np<3).or.((np<=npmax).and.(deltar>aim_deltarmin))) 4464 np=np*2 4465 dt=dt*0.5_dp 4466 vkold(1:3)=vk(1:3) 4467 call vgh_rho(vk,chg,grho,hrho,rr,iat,ipos,0) 4468 vinter(1:3,1)=vv(1:3)+dt*grho(1:3) 4469 do jj=2,np 4470 call vgh_rho(vinter(1,jj-1),chg,grho,hrho,rr,iat,ipos,0) 4471 if(jj.eq.2) then 4472 vinter(1:3,2)=vv(1:3)+2.0*dt*grho(1:3) 4473 else 4474 vinter(1:3,jj)=vinter(1:3,jj-2)+2.0*dt*grho(1:3) 4475 end if 4476 end do 4477 4478 call vgh_rho(vinter(1,np),chg,grho,hrho,rr,iat,ipos,0) 4479 vinter(1:3,np+1)=vinter(1:3,np-1)+dt*grho(1:3) 4480 4481 deltar=0._dp 4482 do ii=1,3 4483 vk(ii)=(vinter(ii,np)+vinter(ii,np+1))*0.5_dp 4484 deltar=deltar+(vkold(ii)-vk(ii))*(vkold(ii)-vk(ii)) 4485 end do 4486 end do 4487 4488 pom(:)=vk(:)-vv(:) 4489 deltar=vnorm(pom,0) 4490 vv(1:3)=vk(1:3) 4491 4492 call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0) 4493 if(deb) write(std_out,*) ':VKf ',np,vk 4494 4495 end subroutine onestep
m_bader/ordr [ Functions ]
[ Top ] [ m_bader ] [ Functions ]
NAME
ordr
FUNCTION
INPUTS
(to be filled)
OUTPUT
(to be filled)
PARENTS
critic
CHILDREN
SOURCE
2479 subroutine ordr(aa,dd,nn,cff) 2480 2481 2482 !This section has been created automatically by the script Abilint (TD). 2483 !Do not modify the following lines by hand. 2484 #undef ABI_FUNC 2485 #define ABI_FUNC 'ordr' 2486 !End of the abilint section 2487 2488 implicit none 2489 2490 !Arguments ---------------------------- 2491 !scalars 2492 integer,intent(in) :: cff,nn 2493 !arrays 2494 real(dp),intent(inout) :: aa(nn),dd(nn,nn) 2495 2496 !Local variables ---------------------- 2497 !scalars 2498 integer :: ii,jj,kk 2499 real(dp) :: uu 2500 2501 ! ********************************************************************* 2502 2503 do ii=1,nn-1 2504 kk=ii 2505 uu=aa(ii) 2506 do jj=ii+1,nn 2507 if (cff==1) then 2508 if (aa(jj) >= uu+tol12) then 2509 kk=jj 2510 uu=aa(jj) 2511 end if 2512 else 2513 if (aa(jj) <= uu-tol12) then 2514 kk=jj 2515 uu=aa(jj) 2516 end if 2517 end if 2518 end do 2519 if (kk /= ii) then 2520 aa(kk)=aa(ii) 2521 aa(ii)=uu 2522 do jj=1,nn 2523 uu=dd(jj,ii) 2524 dd(jj,ii)=dd(jj,kk) 2525 dd(jj,kk)=uu 2526 end do 2527 end if 2528 end do 2529 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)
PARENTS
drvaim
CHILDREN
vgh_rho
SOURCE
4526 subroutine plint() 4527 4528 4529 !This section has been created automatically by the script Abilint (TD). 4530 !Do not modify the following lines by hand. 4531 #undef ABI_FUNC 4532 #define ABI_FUNC 'plint' 4533 !End of the abilint section 4534 4535 implicit none 4536 4537 !Arguments ------------------------------------ 4538 4539 !Local variables ------------------------------ 4540 !scalars 4541 integer,parameter :: nd=150,ng=300 4542 integer :: cod,iat,ii,ipos,jj,kk,nn 4543 real(dp) :: dd,ee,ff,gg,hh,igr,rho,ss 4544 logical :: prep 4545 !arrays 4546 real(dp) :: grho(3),hrho(3,3),vv(3),xl(nd+1),xs(nd) 4547 real(dp),allocatable :: uu(:) 4548 4549 ! ********************************************************************* 4550 4551 ff=rprimd(1,1)/nd 4552 ss=2._dp/sqrt(3._dp)*rprimd(2,2)/rprimd(1,1)*nd 4553 nn=int(ss) 4554 gg=sqrt(3._dp)/2.*ff 4555 hh=rprimd(2,2)-nn/nd*sqrt(3._dp)/2.*rprimd(1,1) 4556 ee=hh/sqrt(3._dp) 4557 hh=hh/2. 4558 ss=sqrt(3._dp)*ff*ff/4. 4559 dd=ee*ff/2. 4560 4561 do ii=1,nd 4562 xl(ii)=ii*ff 4563 xs(ii)=ff/2.+ii*ff 4564 end do 4565 xl(nd+1)=rprimd(1,1) 4566 4567 ABI_ALLOCATE(uu,(nn+3)) 4568 4569 uu(1)=0._dp 4570 uu(nn+3)=rprimd(2,2) 4571 do ii=2,nn+2 4572 uu(ii)=hh+(ii-1)*gg 4573 end do 4574 igr=0._dp 4575 prep=.true. 4576 do kk=1,ng 4577 igr=0._dp 4578 vv(3)=(kk-1)*rprimd(3,3)/ng 4579 do ii=1,nn+3 4580 vv(2)=uu(ii) 4581 do jj=1,nd 4582 if (prep) then 4583 vv(1)=xl(jj) 4584 prep=.false. 4585 else 4586 vv(1)=xs(jj) 4587 prep=.true. 4588 end if 4589 call vgh_rho(vv,rho,grho,hrho,dd,iat,ipos,cod) 4590 if ((ii==1).or.(ii==nn+3)) then 4591 igr=igr+dd*rho 4592 elseif ((ii==2).or.(ii==nn+2)) then 4593 igr=igr+(dd+ss)*rho 4594 else 4595 igr=igr+ss*2*rho 4596 end if 4597 end do 4598 end do 4599 write(untp,'(2E16.8)') vv(3), igr 4600 end do 4601 ABI_DEALLOCATE(uu) 4602 4603 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
PARENTS
drvaim,surf
CHILDREN
aim_follow,timein,vgh_rho
SOURCE
4641 subroutine rsurf(aim_dtset,rr,grho,theta,phi,rr0,iatinit,npmax,srch) 4642 4643 4644 !This section has been created automatically by the script Abilint (TD). 4645 !Do not modify the following lines by hand. 4646 #undef ABI_FUNC 4647 #define ABI_FUNC 'rsurf' 4648 !End of the abilint section 4649 4650 implicit none 4651 4652 !Arguments ------------------------------------ 4653 !scalars 4654 integer,intent(in) :: iatinit,npmax 4655 real(dp),intent(in) :: phi,rr0,theta 4656 real(dp),intent(out) :: rr 4657 logical,intent(in) :: srch 4658 !arrays 4659 real(dp),intent(out) :: grho(3) 4660 !no_abirules 4661 type(aim_dataset_type),intent(in) :: aim_dtset 4662 4663 !Local variables ------------------------------ 4664 !scalars 4665 integer :: iat,ii,ipos,iposinit,jj,nstep 4666 real(dp),parameter :: mfkt=1.d1 4667 real(dp) :: aa,dmax,dr,drr,rho,rr1,rr2,t1,t2,wall 4668 logical :: cross,deb_tmp,in,in1,in2,low,srch_tmp 4669 !arrays 4670 real(dp) :: hrho(3,3),unvec(3),vv(3) 4671 4672 ! ********************************************************************* 4673 4674 srch_tmp=srch 4675 deb_tmp=deb 4676 4677 !unity vecteur in the direction (theta,phi) 4678 4679 unvec(1)=sin(theta)*cos(phi) 4680 unvec(2)=sin(theta)*sin(phi) 4681 unvec(3)=cos(theta) 4682 4683 4684 rr=rr0 4685 rr1=rr 4686 rr2=rr 4687 drr=1._dp 4688 if (abs(rr0-r0)<1.0d-12) then 4689 dr=aim_dtset%dr0*mfkt 4690 else 4691 dr=aim_dtset%dr0 4692 end if 4693 4694 vv(1)=xatm(1,aim_dtset%batom) 4695 vv(2)=xatm(2,aim_dtset%batom) 4696 vv(3)=xatm(3,aim_dtset%batom) 4697 4698 4699 iposinit=batcell 4700 write(std_out,'("ATOM iat=",i4," ipos=",i4)') aim_dtset%batom,batcell 4701 jj=0 4702 4703 cross=.false. 4704 4705 in=.true. 4706 low=.false. 4707 4708 dmax=h0 4709 4710 in1=.true. 4711 in2=in1 4712 4713 do while((drr>aim_drmin).or.(jj<2)) 4714 call timein(t1,wall) 4715 jj=jj+1 4716 do ii=1,3 4717 vv(ii)=xatm(ii,aim_dtset%batom)+rr*unvec(ii) 4718 end do 4719 4720 ! VACUUM CONDITION 4721 4722 call vgh_rho(vv,rho,grho,hrho,aa,iat,ipos,0) 4723 if (rho < aim_rhomin) exit 4724 4725 ldeb=.false. 4726 4727 call aim_follow(aim_dtset,vv,npmax,srch_tmp,iatinit,iposinit,iat,ipos,nstep) 4728 4729 call timein(t2,wall) 4730 t2=t2-t1 4731 4732 write(std_out,'(a,i4,a,f12.8,a,i4,a,i4,a,f10.5,a,i4)') & 4733 & ' :STEP ',jj,' r=',rr,' iat=',iat,' ipos=',ipos,' time(sec)=',t2,' nstep=',nstep 4734 4735 if ((iat.eq.iatinit).and.(ipos.eq.iposinit)) then 4736 in=.true. 4737 else 4738 in=.false. 4739 end if 4740 4741 ! 4742 ! NEW RADIUS 4743 ! 4744 4745 if ((jj.eq.1).or.((in1.eqv.in).and.(.not.cross))) then 4746 if (in) then 4747 rr2=rr1 4748 rr1=rr 4749 rr=rr+dr 4750 else 4751 rr2=rr1 4752 rr1=rr 4753 rr=rr-dr 4754 end if 4755 if ((jj>2).and.(dr<(0.6))) then 4756 ! modification of the step 4757 dr=dr*aim_fac 4758 if (deb_tmp) write(std_out,*) ':DR ',dr 4759 end if 4760 else 4761 if (.not.cross) then 4762 cross=.true. 4763 rr2=rr1 4764 else 4765 if (in2) then 4766 if (in) then 4767 rr2=rr1 4768 else 4769 in1=in2 4770 end if 4771 else 4772 if (in) then 4773 in1=in2 4774 else 4775 rr2=rr1 4776 end if 4777 end if 4778 end if 4779 rr1=rr 4780 rr=(rr2+rr1)/2.0 4781 end if 4782 4783 in2=in1 4784 in1=in 4785 drr=abs(rr2-rr1)/rr 4786 if (deb_tmp) write(std_out,*) ':DRR ',jj,rr2,rr1,drr 4787 end do 4788 4789 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)
PARENTS
drvaim
CHILDREN
coeffs_gausslegint,rsurf,timein,xmpi_sum
SOURCE
4826 subroutine surf(aim_dtset) 4827 4828 4829 !This section has been created automatically by the script Abilint (TD). 4830 !Do not modify the following lines by hand. 4831 #undef ABI_FUNC 4832 #define ABI_FUNC 'surf' 4833 !End of the abilint section 4834 4835 implicit none 4836 4837 !Arguments ------------------------------------ 4838 !scalars 4839 type(aim_dataset_type) :: aim_dtset 4840 4841 !Local variables ------------------------------ 4842 !scalars 4843 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 4844 real(dp) :: ct1,ct2,phi,rr,rsmax,rsmin,rthe,rthe0,t1,t2,theta,tt0,vcth,vph,vth 4845 real(dp) :: wall,xy,xyz 4846 logical :: srch,stemp 4847 !arrays 4848 real(dp) :: grho(3),vr(3),vv(3) 4849 real(dp),allocatable :: rs_computed(:,:) 4850 4851 !************************************************************************ 4852 4853 comm = xmpi_world 4854 me=xmpi_comm_rank(comm) 4855 nproc=xmpi_comm_size(comm) 4856 4857 ttsrf=zero 4858 4859 rewind(unts) 4860 4861 nth=aim_dtset%nth 4862 nph=aim_dtset%nph 4863 4864 !Coefficients for spherical Gauss quadrature 4865 4866 ct1=cos(aim_dtset%themin) 4867 ct2=cos(aim_dtset%themax) 4868 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 4869 call coeffs_gausslegint(aim_dtset%phimin,aim_dtset%phimax,ph,wph,nph) 4870 4871 !DEBUG 4872 !write(std_out,*)' surf : wcth=',wcth(1:nth) 4873 !write(std_out,*)' surf : wph=',wph(1:nth) 4874 !ENDDEBUG 4875 4876 do ijj=1,nth 4877 th(ijj)=acos(cth(ijj)) 4878 if (aim_dtset%isurf/=-1) then 4879 do jj=1,nph 4880 rs(ijj,jj)=zero 4881 end do 4882 end if 4883 end do 4884 4885 npmax=aim_npmaxin 4886 rsmax=0.0 4887 rsmin=100.0 4888 rthe0=r0 4889 srch=.false. 4890 4891 do ijj=1,3 4892 vv(ijj)=xatm(ijj,aim_dtset%batom) 4893 end do 4894 4895 4896 write(std_out,*) 4897 write(std_out,*) "BADER SURFACE DETERMINATION" 4898 write(std_out,*) "===========================" 4899 write(std_out,*) 4900 4901 write(untout,*) 4902 write(untout,*) "BADER SURFACE DETERMINATION" 4903 write(untout,*) "===========================" 4904 write(untout,*) 4905 4906 write(std_out,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 4907 write(std_out,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4908 write(std_out,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4909 4910 write(untout,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 4911 write(untout,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4912 write(untout,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4913 4914 write(unts,'(i3,3F15.10)') aim_dtset%batom,vv 4915 write(unts,'(i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 4916 write(unts,'(i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 4917 4918 !write(std_out,*) 'npmax in surf= ',npmax 4919 4920 ith=0 4921 iph=0 4922 tt0=0._dp 4923 call timein(tt0,wall) 4924 4925 write(untout,*) 4926 write(untout,*) "DEVELOPMENT OF THE RADII DETERMINATIONS" 4927 write(untout,*) "========================================" 4928 write(untout,*) 4929 write(untout,*) "Determination near the CPs:" 4930 4931 !Determination of the CP neighbouring radii 4932 4933 if (aim_dtset%isurf/=-1) then 4934 4935 ! Precomputation of the value of the radii (for parallelisation) 4936 ! To make the output independent of the number of processors, but still 4937 ! cut down the CPU time, use a multigrid technique 4938 srch=.true. 4939 ABI_ALLOCATE(rs_computed,(nth,nph)) 4940 rs(:,:)=zero 4941 rs_computed(:,:)=zero 4942 kk=0 ; init=0 4943 do level=3,0,-1 4944 incr=2**level 4945 if(incr<nth .and. incr<nph)then 4946 rs_computed(:,:)=rs(1:nth,1:nph) 4947 rs(1:nth,1:nph)=zero 4948 do ijj=1,nth,incr 4949 do jj=1,nph,incr 4950 if(rs_computed(ijj,jj)<1.0d-12) then 4951 kk=kk+1 4952 if(mod(kk,nproc)==me)then 4953 ! Find an approximate starting radius, from the already computed ones 4954 if(init==0)then 4955 rthe=r0 4956 else 4957 ijj_exist=ijj ; if(mod(ijj-1,2*incr)>=incr)ijj_exist=ijj-incr 4958 jj_exist=jj ; if(mod(jj-1,2*incr)>=incr)jj_exist=jj-incr 4959 rthe=rs_computed(ijj_exist,jj_exist) 4960 if(rthe<1.0d-12)then 4961 write(std_out,*)' surf : there is a bug ! rthe=',rthe 4962 MSG_ERROR("Aborting now") 4963 end if 4964 end if 4965 call timein(t1,wall) ; t2=zero 4966 call rsurf(aim_dtset,rr,grho,th(ijj),ph(jj),rthe,aim_dtset%batom,npmax,srch) 4967 rs(ijj,jj)=rr 4968 if (deb) then 4969 call timein(t2,wall) ; t2=t2-t1 4970 write(std_out,*) ':CALCULATED NP',ijj,jj,th(ijj),ph(jj),rthe,npmax,rs(ijj,jj),t2 4971 end if 4972 end if 4973 end if 4974 end do ! jj 4975 end do ! ijj 4976 call xmpi_sum(rs,comm,ierr) 4977 ! Combine the set of already computed radii and the set of the newly computed, to obtain all computed. 4978 rs(1:nth,1:nph)=rs(1:nth,1:nph)+rs_computed(:,:) 4979 init=1 4980 end if 4981 end do 4982 ABI_DEALLOCATE(rs_computed) 4983 4984 srch=.true. 4985 4986 do ijj=1,nbcp 4987 ! if ((icpc(ijj) == -1)) then 4988 rthe0=vnorm(pc(:,ijj),0) 4989 do jj=1,3 4990 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 4991 end do 4992 xy=vr(1)*vr(1)+vr(2)*vr(2) 4993 xyz=xy+vr(3)*vr(3) 4994 xyz=sqrt(xyz) 4995 4996 if (xy < aim_xymin) then 4997 vcth=1._dp 4998 if (vr(3) < 0._dp) vcth=-vcth 4999 vph=0._dp 5000 else 5001 vcth=vr(3)/xyz 5002 vph=atan2(vr(2),vr(1)) 5003 end if 5004 5005 vth=acos(vcth) 5006 write(untout,'(/," BCP: (index,theta,phi)",I4,2E16.8)') ijj,vth,vph 5007 5008 if (vth < th(1)) then 5009 ith=0 5010 else 5011 if (vth > th(nth)) then 5012 ith=nth 5013 else 5014 do ii=2,nth 5015 if (vth < th(ii)) then 5016 ith=ii-1 5017 exit 5018 end if 5019 end do 5020 end if 5021 end if 5022 5023 if (vph < ph(1)) then 5024 iph=0 5025 else 5026 if (vph > ph(nph)) then 5027 iph=nph 5028 else 5029 do ii=2,nph 5030 if (vph < ph(ii)) then 5031 iph=ii-1 5032 exit 5033 end if 5034 end do 5035 end if 5036 end if 5037 5038 write(untout,*) "ATOMIC RADII (ith,iph,theta,phi,radius)" 5039 do jj=-1,2 5040 do kk=-1,2 5041 ith2=ith+jj 5042 iph2=iph+kk 5043 stemp=(iph2 > 0).and.(iph2 < nph+1) 5044 stemp=(stemp.and.((ith2 > 0).and.(ith2 < nth+1))) 5045 if (stemp) then 5046 theta=th(ith2) 5047 phi=ph(iph2) 5048 if (abs(rs(ith2,iph2))<1.0d-12) then 5049 rthe=rthe0 5050 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 5051 call timein(t1,wall) 5052 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5053 call timein(t2,wall) 5054 t2=t2-t1 5055 rs(ith2,iph2)=rr 5056 end if 5057 rr=rs(ith2,iph2) 5058 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 5059 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 5060 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 5061 rthe0=rr 5062 end if 5063 5064 end do ! kk 5065 end do ! jj 5066 5067 ! end if 5068 5069 end do ! ijj (loop on BCP) 5070 5071 ! DEBUG 5072 ! write(std_out,*)' surf : near BCP ' 5073 ! do ijj=1,nth 5074 ! do jj=1,nph 5075 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5076 ! end do 5077 ! end do 5078 ! ENDDEBUG 5079 5080 5081 srch=.true. 5082 do ijj=nbcp+1,nbcp+nrcp ! Loop on RCP 5083 ! if ((icpc(ijj) == 1)) then 5084 rthe0=max(rminl(aim_dtset%batom),r0) 5085 do jj=1,3 5086 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 5087 end do 5088 xy=vr(1)*vr(1)+vr(2)*vr(2) 5089 xyz=xy+vr(3)*vr(3) 5090 xyz=sqrt(xyz) 5091 5092 if (xy < aim_xymin) then 5093 vcth=1._dp 5094 if (vr(3) < 0._dp) vcth=-vcth 5095 vph=0._dp 5096 else 5097 vcth=vr(3)/xyz 5098 vph=atan2(vr(2),vr(1)) 5099 end if 5100 vth=acos(vcth) 5101 write(untout,'(/,";RCP: (index,theta,phi)",I4,2E16.8)') ijj-nbcp,vth,vph 5102 5103 if (vth < th(1)) then 5104 ith=0 5105 else 5106 if (vth > th(nth)) then 5107 ith=nth 5108 else 5109 do ii=2,nth 5110 if (vth < th(ii)) then 5111 ith=ii-1 5112 exit 5113 end if 5114 end do 5115 end if 5116 end if 5117 5118 if (vph < ph(1)) then 5119 iph=0 5120 else 5121 if (vph > ph(nph)) then 5122 iph=nph 5123 else 5124 do ii=2,nph 5125 if (vph < ph(ii)) then 5126 iph=ii-1 5127 exit 5128 end if 5129 end do 5130 end if 5131 end if 5132 5133 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 5134 do jj=-1,2 5135 do kk=-1,2 5136 ith2=ith+jj 5137 iph2=iph+kk 5138 stemp=(iph2 > 0).and.(iph2 < nph+1) 5139 stemp=stemp.and.(ith2 > 0).and.(ith2 < nth+1) 5140 5141 if (stemp) then 5142 theta=th(ith2) 5143 phi=ph(iph2) 5144 if ((abs(rs(ith2,iph2))<1.0d-12)) then 5145 rthe=rthe0 5146 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 5147 call timein(t1,wall) 5148 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5149 call timein(t2,wall) 5150 t2=t2-t1 5151 rs(ith2,iph2)=rr 5152 end if 5153 rr=rs(ith2,iph2) 5154 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 5155 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 5156 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 5157 rthe0=rr 5158 end if 5159 5160 end do ! kk 5161 end do ! jj 5162 ! end if 5163 5164 end do ! ijj (Loop on RCP) 5165 5166 ! DEBUG 5167 ! write(std_out,*)' surf : near RCP ' 5168 ! do ijj=1,nth 5169 ! do jj=1,nph 5170 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5171 ! end do 5172 ! end do 5173 ! ENDDEBUG 5174 5175 ! Boundary angles 5176 rthe0=r0 5177 srch=.true. 5178 write(untout,*) 5179 write(untout,*) "The boundary angles:" 5180 write(untout,*) "====================" 5181 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 5182 5183 ! Must have sufficient angular sampling 5184 if ((nth > 8).and.(nph > 8)) then 5185 rthe=r0 5186 do ijj=1,2 5187 theta=th(ijj) 5188 if (ijj==2) rthe=rs(1,1) 5189 do jj=1,nph 5190 phi=ph(jj) 5191 call timein(t1,wall) 5192 if (abs(rs(ijj,jj))<1.0d-12) then 5193 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5194 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5195 rs(ijj,jj)=rr 5196 end if 5197 rr=rs(ijj,jj) 5198 call timein(t2,wall) 5199 t2=t2-t1 5200 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5201 write(untout,'(a,2i3,3E16.8)') '- ',ijj,jj,theta,phi,rr 5202 rthe=rs(ijj,jj) 5203 end do ! jj 5204 end do ! ijj 5205 5206 write(untout,*) 5207 5208 rthe=rs(2,1) 5209 do jj=1,2 5210 phi=ph(jj) 5211 if (jj==2) rthe=rs(2,2) 5212 do ijj=3,nth 5213 theta=th(ijj) 5214 t2=0.0 5215 call timein(t1,wall) 5216 if (abs(rs(ijj,jj))<1.0d-12) then 5217 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5218 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5219 rs(ijj,jj)=rr 5220 end if 5221 rr=rs(ijj,jj) 5222 call timein(t2,wall) 5223 t2=t2-t1 5224 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5225 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 5226 rthe=rs(ijj,jj) 5227 end do ! ijj 5228 end do ! jj 5229 5230 write(untout,*) 5231 5232 rthe=rs(nth-1,2) 5233 do ijj=nth-1,nth 5234 theta=th(ijj) 5235 if (ijj==nth) rthe=rs(nth,2) 5236 do jj=3,nph 5237 phi=ph(jj) 5238 call timein(t1,wall) 5239 if (abs(rs(ijj,jj))<1.0d-12) then 5240 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5241 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5242 rs(ijj,jj)=rr 5243 end if 5244 rr=rs(ijj,jj) 5245 call timein(t2,wall) 5246 t2=t2-t1 5247 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5248 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 5249 rthe=rs(ijj,jj) 5250 end do ! jj 5251 end do ! ijj 5252 5253 rthe=rs(2,nph-1) 5254 do jj=nph-1,nph 5255 phi=ph(jj) 5256 if (jj==nph) rthe=rs(2,nph) 5257 do ijj=3,nth-2 5258 theta=th(ijj) 5259 t2=0.0 5260 call timein(t1,wall) 5261 if (abs(rs(ijj,jj))<1.0d-12) then 5262 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5263 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5264 rs(ijj,jj)=rr 5265 end if 5266 rr=rs(ijj,jj) 5267 call timein(t2,wall) 5268 t2=t2-t1 5269 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5270 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 5271 rthe=rs(ijj,jj) 5272 end do ! ijj 5273 end do ! jj 5274 write(untout,*) 5275 5276 ! Complementary bands for boundary angles 5277 nn=int(real(nth)/1.4d1) 5278 if (nn > 1) then 5279 do ii=1,nn-1 5280 mm=int(nth/nn)*ii 5281 do kk=0,1 5282 mm=mm+kk 5283 theta=th(mm) 5284 rthe=rs(mm,2) 5285 do jj=3,nph-2 5286 phi=ph(jj) 5287 call timein(t1,wall) 5288 if (abs(rs(mm,jj))<1.0d-12) then 5289 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5290 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5291 rs(mm,jj)=rr 5292 end if 5293 rr=rs(mm,jj) 5294 call timein(t2,wall) 5295 t2=t2-t1 5296 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(mm)*wph(jj),t2 5297 write(untout,'(2i3,3E16.8)') mm,jj,theta,phi,rr 5298 rthe=rs(mm,jj) 5299 end do ! jj 5300 end do ! kk 5301 end do ! ii 5302 end if ! nn>1 5303 5304 write(untout,*) 5305 5306 nn=nint(real(nph)/1.2d1) 5307 if (nn > 1) then 5308 do ii=1,nn-1 5309 mm=int(nph/nn)*ii 5310 do kk=0,1 5311 mm=mm+kk 5312 phi=ph(mm) 5313 rthe=rs(2,mm) 5314 5315 do jj=3,nth-2 5316 theta=th(jj) 5317 call timein(t1,wall) 5318 if (abs(rs(jj,mm))<1.0d-12) then 5319 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 5320 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 5321 rs(jj,mm)=rr 5322 end if 5323 rr=rs(mm,jj) 5324 call timein(t2,wall) 5325 t2=t2-t1 5326 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(jj)*wph(mm),t2 5327 write(untout,'(2i3,3E16.8)') jj,mm,theta,phi,rr 5328 rthe=rs(jj,mm) 5329 end do ! jj 5330 5331 end do ! kk 5332 end do ! ii 5333 end if ! nn>1 5334 5335 end if ! sufficient sampling to determine boundary angles 5336 5337 write(untout,*) 5338 5339 ! DEBUG 5340 ! write(std_out,*)' surf : after boundary angles ' 5341 ! do ijj=1,nth 5342 ! do jj=1,nph 5343 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5344 ! end do 5345 ! end do 5346 ! ENDDEBUG 5347 5348 ! Output the complete Bader surface 5349 5350 write(untout,*) "The complete Bader surface:" 5351 write(untout,*) "===========================" 5352 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 5353 rthe0=r0 5354 srch=.true. 5355 5356 ! Write all the values 5357 5358 do ijj=1,nth 5359 theta=th(ijj) 5360 do jj=1,nph 5361 phi=ph(jj) 5362 rr=rs(ijj,jj) 5363 write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 5364 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 5365 write(untout,'(a,2i3,3E16.8)') ' ',ijj,jj,theta,phi,rr 5366 if (rr < rsmin) rsmin=rr 5367 if (rr> rsmax) rsmax=rr 5368 end do ! jj 5369 end do ! ijj 5370 write(unts,'(2F15.10)') rsmin,rsmax 5371 write(untout,'(/," The minimal and maximal radii:",/,/," ",2F15.10)') rsmin,rsmax 5372 5373 ! DEBUG 5374 ! write(std_out,*)' surf : final output ' 5375 ! do ijj=1,nth 5376 ! do jj=1,nph 5377 ! write(std_out,*)ijj,jj,rs(ijj,jj) 5378 ! end do 5379 ! end do 5380 ! ENDDEBUG 5381 5382 end if ! determination of the critical surface 5383 5384 call timein(ttsrf,wall) 5385 ttsrf=ttsrf-tt0 5386 5387 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
6032 function vec_prod(uu,vv) 6033 6034 6035 !This section has been created automatically by the script Abilint (TD). 6036 !Do not modify the following lines by hand. 6037 #undef ABI_FUNC 6038 #define ABI_FUNC 'vec_prod' 6039 !End of the abilint section 6040 6041 implicit none 6042 6043 !Arguments ------------------------------------ 6044 !arrays 6045 real(dp) :: vec_prod(3) 6046 real(dp),intent(in) :: uu(3),vv(3) 6047 6048 !Local variables------------------------------- 6049 6050 ! ************************************************************************* 6051 6052 vec_prod(1)=uu(2)*vv(3)-vv(2)*uu(3) 6053 vec_prod(2)=uu(3)*vv(1)-vv(3)*uu(1) 6054 vec_prod(3)=uu(1)*vv(2)-vv(1)*uu(2) 6055 6056 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
PARENTS
addout,aim_follow,cpdrv,critic,critics,integrho,onestep,plint,rsurf
CHILDREN
bschg1,bschg2
SOURCE
5427 subroutine vgh_rho(vv,rho,grho,hrho,rdmin,iat,ipos,chs) 5428 5429 5430 !This section has been created automatically by the script Abilint (TD). 5431 !Do not modify the following lines by hand. 5432 #undef ABI_FUNC 5433 #define ABI_FUNC 'vgh_rho' 5434 !End of the abilint section 5435 5436 implicit none 5437 5438 !Arguments ------------------------------------ 5439 !scalars 5440 integer,intent(in) :: chs 5441 integer,intent(inout) :: iat,ipos 5442 real(dp),intent(out) :: rdmin,rho 5443 !arrays 5444 real(dp),intent(in) :: vv(3) 5445 real(dp),intent(out) :: grho(3),hrho(3,3) 5446 5447 !Local variables ------------------------------ 5448 !scalars 5449 integer :: ii,inmax,inmin,inx,jj,kk,ll,nn,oii,omm,onn 5450 integer :: selct 5451 ! real(dp),save :: cumul_cpu=0.0_dp,cumul_cpu_old=0.0_dp 5452 real(dp),save :: tcpui,tcpuo,twalli 5453 real(dp),save :: twallo 5454 real(dp) :: aa,bb,cc,cgrad1_rr_inv,coeff,dd,rr,rr2,rr_inv 5455 real(dp) :: rrad2_nn,rrad_nn,ss,uu,uu_inv,val,vt1,vt2,vt3,vw1,vw2 5456 ! real(dp) :: ss_inv 5457 real(dp) :: vw3 5458 !arrays 5459 integer :: indx(3),inii(4,3) 5460 real(dp) :: cgrad(3),ches(3,3),cof(2,3),ddstar(6),ddu(2),grd(4) 5461 real(dp) :: hh(4,2),hrh(2),lder(4),pom2sq(2,3),pomsq(2,3) 5462 real(dp) :: rhstar(6),sqder(6,4),sqvlr(6,4),trsf(3,3),xx(3) 5463 real(dp),pointer :: ptddx(:,:,:),ptddy(:,:,:),ptddz(:,:,:),ptrho(:,:,:) 5464 5465 !************************************************************************ 5466 tcpui=0.0_dp 5467 tcpuo=0.0_dp 5468 twalli=0.0_dp 5469 twallo=0.0_dp 5470 5471 nullify(ptddx,ptddy,ptddz,ptrho) 5472 5473 selct=chs 5474 5475 if (selct/=2) then 5476 5477 ! call timein(tcpui,twalli) 5478 5479 ! TRANSFORMATION TO THE REDUCED COORD. 5480 5481 xx(:)=vv(:) 5482 call bschg1(xx,-1) 5483 5484 ! call timein(tcpuo,twallo) 5485 ! cumul_cpu=cumul_cpu+(tcpuo-tcpui) 5486 5487 ! REDUCTION TO THE PRIMITIVE CELL 5488 5489 do ii=1,3 5490 if (xx(ii) >= one-tol12 ) then 5491 xx(ii)=xx(ii)-aint(xx(ii)) 5492 elseif (xx(ii) < -tol12 ) then 5493 xx(ii)=xx(ii)-floor(xx(ii)) 5494 end if 5495 end do 5496 5497 5498 ! DETERMINATION OF THE INDEX IN THE GRID 5499 5500 do ii=1,3 5501 indx(ii)=aint(xx(ii)*ngfft(ii)) 5502 bb=(xx(ii)-indx(ii)*dix(ii))*ngfft(ii) 5503 if (indx(ii)==ngfft(ii)) then 5504 indx(ii)=1 5505 xx(ii)=0._dp 5506 else 5507 indx(ii)=indx(ii)+1 5508 end if 5509 5510 ! Explicit handling to avoid numeric problems 5511 5512 if (bb > 1._dp+tol12 ) then 5513 cof(1,ii)=0._dp 5514 cof(2,ii)=1._dp 5515 elseif (bb < -tol12 ) then 5516 cof(1,ii)=1._dp 5517 cof(2,ii)=0._dp 5518 else 5519 cof(1,ii)=1._dp-bb 5520 cof(2,ii)=bb 5521 end if 5522 end do 5523 5524 ! 3D INTERPOLATION OF THE VALENCE DENSITY 5525 5526 ! determination of the values of density and of its second derivative 5527 ! at the "star" = constructed at vv with primitive directions 5528 ! To interpolation the values at the faces of the grid cell are needed 5529 5530 rhstar(:)=0._dp 5531 sqder(:,:)=0._dp 5532 sqvlr(:,:)=0._dp 5533 ddstar(:)=0._dp 5534 pomsq(:,:)=0._dp 5535 pom2sq(:,:)=0._dp 5536 5537 oii=1; onn=1; omm=1 5538 if (indx(1)==ngfft(1)) oii=1-ngfft(1) 5539 if (indx(2)==ngfft(2)) onn=1-ngfft(2) 5540 if (indx(3)==ngfft(3)) omm=1-ngfft(3) 5541 5542 ! the values in the corners of the grid cell 5543 5544 ptddx=>ddx(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5545 ptddy=>ddy(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5546 ptddz=>ddz(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5547 ptrho=>dvl(indx(1):indx(1)+oii:oii,indx(2):indx(2)+onn:onn,indx(3):indx(3)+omm:omm) 5548 5549 ! the coefficients for spline interpolation of density and its derivation 5550 do ii=1,3 5551 do jj=1,2 5552 pomsq(jj,ii)=(cof(jj,ii)*cof(jj,ii)*cof(jj,ii)-cof(jj,ii))/6._dp*dix(ii)*dix(ii) 5553 pom2sq(jj,ii)=(3._dp*cof(jj,ii)*cof(jj,ii)-1._dp)/6._dp*dix(ii) 5554 if (jj==1) pom2sq(jj,ii)=-pom2sq(jj,ii) 5555 end do 5556 end do 5557 5558 5559 do ii=1,2 5560 do jj=1,2 5561 do kk=1,2 5562 ddstar(ii)=ddstar(ii)+cof(jj,2)*cof(kk,3)*ptddx(ii,jj,kk) 5563 ddstar(ii+2)=ddstar(ii+2)+cof(jj,3)*cof(kk,1)*ptddy(kk,ii,jj) 5564 ddstar(ii+4)=ddstar(ii+4)+cof(jj,1)*cof(kk,2)*ptddz(jj,kk,ii) 5565 sqder(ii,jj)=sqder(ii,jj)+cof(kk,2)*ptddz(ii,kk,jj) 5566 sqder(ii,jj+2)=sqder(ii,jj+2)+cof(kk,3)*ptddy(ii,jj,kk) 5567 sqder(ii+2,jj)=sqder(ii+2,jj)+cof(kk,3)*ptddx(jj,ii,kk) 5568 sqder(ii+2,jj+2)=sqder(ii+2,jj+2)+cof(kk,1)*ptddz(kk,ii,jj) 5569 sqder(ii+4,jj)=sqder(ii+4,jj)+cof(kk,1)*ptddy(kk,jj,ii) 5570 sqder(ii+4,jj+2)=sqder(ii+4,jj+2)+cof(kk,2)*ptddx(jj,kk,ii) 5571 sqvlr(ii,jj)=sqvlr(ii,jj)+cof(kk,2)*ptrho(ii,kk,jj)+pomsq(kk,2)*ptddy(ii,kk,jj) 5572 sqvlr(ii,jj+2)=sqvlr(ii,jj+2)+cof(kk,3)*ptrho(ii,jj,kk)+pomsq(kk,3)*ptddz(ii,jj,kk) 5573 sqvlr(ii+2,jj+2)=sqvlr(ii+2,jj+2)+cof(kk,1)*ptrho(kk,ii,jj)+pomsq(kk,1)*ptddx(kk,ii,jj) 5574 end do 5575 end do 5576 end do 5577 5578 do ii=1,2 5579 do jj=1,2 5580 sqvlr(ii+2,jj)=sqvlr(jj,ii+2) 5581 sqvlr(ii+4,jj)=sqvlr(jj+2,ii+2) 5582 sqvlr(ii+4,jj+2)=sqvlr(jj,ii) 5583 end do 5584 end do 5585 5586 do ii=1,2 5587 do jj=1,2 5588 rhstar(ii)=rhstar(ii)+cof(jj,3)*sqvlr(ii,jj)+pomsq(jj,3)*sqder(ii,jj)+& 5589 & cof(jj,2)*sqvlr(ii,jj+2)+pomsq(jj,2)*sqder(ii,jj+2) 5590 rhstar(ii+2)=rhstar(ii+2)+cof(jj,1)*sqvlr(ii+2,jj)+pomsq(jj,1)*sqder(ii+2,jj)+& 5591 & cof(jj,3)*sqvlr(ii+2,jj+2)+pomsq(jj,3)*sqder(ii+2,jj+2) 5592 rhstar(ii+4)=rhstar(ii+4)+cof(jj,2)*sqvlr(ii+4,jj)+pomsq(jj,2)*sqder(ii+4,jj)+& 5593 & cof(jj,1)*sqvlr(ii+4,jj+2)+pomsq(jj,1)*sqder(ii+4,jj+2) 5594 end do 5595 end do 5596 rhstar(:)=rhstar(:)/2._dp 5597 5598 rho=0._dp 5599 grho(:)=0._dp 5600 hrho(:,:)=0._dp 5601 kk=1; nn=1 5602 do ii=1,5,2 5603 do jj=1,2 5604 nn=-nn 5605 rho=rho+cof(jj,kk)*rhstar(ii+jj-1)+pomsq(jj,kk)*ddstar(ii+jj-1) 5606 grho(kk)=grho(kk)+pom2sq(jj,kk)*ddstar(ii+jj-1) 5607 hrho(kk,kk)=hrho(kk,kk)+cof(jj,kk)*ddstar(ii+jj-1) 5608 grho(kk)=grho(kk)+nn*rhstar(ii+jj-1)/dix(kk) 5609 end do 5610 kk=kk+1 5611 end do 5612 rho=rho/3._dp 5613 5614 ! Off-diagonal elements of the hessian 5615 5616 ! for the speed reasons the polynomial interpolation 5617 ! for second derivation fields is used in this case 5618 ! but the last step is always done by spline interpolation. 5619 5620 5621 do ii=1,3 5622 do jj=-1,2 5623 inii(jj+2,ii)=indx(ii)+jj 5624 if (inii(jj+2,ii) < 1) inii(jj+2,ii)=inii(jj+2,ii)+ngfft(ii) 5625 if (inii(jj+2,ii) > ngfft(ii)) inii(jj+2,ii)=inii(jj+2,ii)-ngfft(ii) 5626 end do 5627 end do 5628 5629 ! Not very nice 5630 5631 do ii=1,3 5632 select case (ii) 5633 case (1) 5634 do jj=1,4 5635 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)) 5636 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)) 5637 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))+& 5638 & 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)) 5639 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))+& 5640 & 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)) 5641 hh(jj,2)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 5642 5643 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)) 5644 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)) 5645 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))+& 5646 & 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)) 5647 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))+& 5648 & 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)) 5649 hh(jj,1)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 5650 end do 5651 case (2) 5652 do jj=1,4 5653 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)) 5654 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)) 5655 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))+& 5656 & 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)) 5657 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))+& 5658 & 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)) 5659 hh(jj,2)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 5660 5661 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)) 5662 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)) 5663 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))+& 5664 & 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)) 5665 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))+& 5666 & 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)) 5667 hh(jj,1)=(hrh(2)-hrh(1))/dix(3)+pom2sq(1,3)*ddu(1)+pom2sq(2,3)*ddu(2) 5668 end do 5669 case (3) 5670 do jj=1,4 5671 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)) 5672 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)) 5673 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))+& 5674 & 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)) 5675 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))+& 5676 & 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)) 5677 hh(jj,2)=(hrh(2)-hrh(1))/dix(2)+pom2sq(1,2)*ddu(1)+pom2sq(2,2)*ddu(2) 5678 5679 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)) 5680 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)) 5681 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))+& 5682 & 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)) 5683 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))+& 5684 & 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)) 5685 hh(jj,1)=(hrh(2)-hrh(1))/dix(1)+pom2sq(1,1)*ddu(1)+pom2sq(2,1)*ddu(2) 5686 end do 5687 end select 5688 do jj=-2,1 5689 grd(jj+3)=(indx(ii)+jj)*dix(ii) 5690 end do 5691 5692 ! write(std_out,'("hh: ",/,4F16.8,/,4F16.8)') ((hh(kk,jj),kk=1,4),jj=1,2) 5693 ! write(std_out,'("grad: ",3F16.8)') (grho(kk),kk=1,3) 5694 ! write(std_out,'("dix: ",3F16.8)') (dix(kk),kk=1,3) 5695 ! write(std_out,'("grd: ",4F16.8)') (grd(kk),kk=1,4) 5696 ! write(std_out,'("inii: ",4I4)') (inii(kk,ii),kk=1,4) 5697 5698 do jj=1,2 5699 5700 ! polynomial interpolation 5701 5702 do kk=1,3 5703 do ll=4,kk+1,-1 5704 hh(ll,jj)=(hh(ll,jj)-hh(ll-1,jj))/(grd(ll)-grd(ll-1)) 5705 end do 5706 end do 5707 lder(4)=hh(4,jj) 5708 do kk=3,1,-1 5709 lder(kk)=hh(kk,jj)+(xx(ii)-grd(kk))*lder(kk+1) 5710 end do 5711 do kk=1,2 5712 do ll=3,kk+1,-1 5713 lder(ll)=lder(ll)+(xx(ii)-grd(ll-kk))*lder(ll+1) 5714 end do 5715 end do 5716 nn=ii+jj 5717 if (nn > 3) nn=nn-3 5718 hrho(ii,nn)=hrho(ii,nn)+lder(2) 5719 hrho(nn,ii)=hrho(nn,ii)+lder(2) 5720 end do 5721 end do 5722 5723 ! averaging of the mixed derivations obtained in different order 5724 5725 do ii=1,3 5726 do jj=1,3 5727 if (ii /= jj) hrho(ii,jj)=hrho(ii,jj)/2._dp 5728 end do 5729 end do 5730 5731 5732 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 5733 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 5734 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 5735 ! stop 5736 ! write(std_out,'("xx:",3F16.8)') (xx(ii),ii=1,3) 5737 ! write(std_out,'(":GRAD pred tr ",3F16.8)') grho 5738 ! write(std_out,'(":HESSIAN pred tr",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5739 5740 5741 ! Transformation back to Cart. coordonnes 5742 5743 call bschg1(grho,2) 5744 call bschg2(hrho,2) 5745 5746 ! write(std_out,'("hrho: ",/,3F16.8,/,3F16.8,/,3F16.8)') & 5747 ! & ((hrho(ii,jj),ii=1,3),jj=1,3) 5748 ! stop 5749 5750 nullify(ptddx,ptddy,ptddz,ptrho) 5751 5752 if (selct==1) return 5753 5754 end if 5755 5756 !write(51,'(":GRADv ",3F16.8)') grho 5757 !write(52,'(":LAPv ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 5758 !write(52,'(":HESNv ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5759 5760 !INTERPOLATION OF THE CORE DENSITY 5761 5762 if (selct/=1) then 5763 5764 if (selct==2) then 5765 grho(:)=0._dp 5766 hrho(:,:)=0._dp 5767 rho=0._dp 5768 end if 5769 5770 ! SEARCH OF THE NEIGHBOUR ATOMS 5771 5772 if (selct /= -2) then 5773 iat=0 5774 ipos=0 5775 end if 5776 rdmin=20._dp 5777 5778 do jj=1,natom 5779 nn=typat(jj) 5780 rrad_nn=rrad(corlim(nn),nn) 5781 rrad2_nn=rrad_nn*rrad_nn 5782 vw1=vv(1)-xatm(1,jj) 5783 vw2=vv(2)-xatm(2,jj) 5784 vw3=vv(3)-xatm(3,jj) 5785 5786 do kk=1,nnpos 5787 5788 vt1=vw1-atp(1,kk) 5789 vt2=vw2-atp(2,kk) 5790 vt3=vw3-atp(3,kk) 5791 rr2=vt1*vt1+vt2*vt2+vt3*vt3 5792 5793 ! rr=vnorm(vt,0) 5794 5795 ! Only contribution > rhocormin (adhoc.f90) are considered 5796 5797 if (rr2 < rrad2_nn .and.(.not.((selct==-2).and.(iat==jj).and.(ipos==kk)))) then 5798 ! if (rr /= 0.0_dp) then ! XG020629 : never test a real number against zero (not portable) 5799 if (rr2 > 1.0d-28) then ! SEARCH INDEX 5800 5801 rr=sqrt(rr2) 5802 rr_inv=1.0_dp/rr 5803 5804 if (rr < rrad(1,nn)) then 5805 inx=-1 5806 elseif (rr > rrad(ndat(nn),nn)) then 5807 inx=ndat(nn) 5808 else 5809 ! Find the index of the radius by bissection 5810 inmin=1 5811 inmax=ndat(nn) 5812 inx=1 5813 do 5814 if(inmax-inmin==1)exit 5815 inx=(inmin+inmax)/2 5816 if(rr>=rrad(inx,nn))then 5817 inmin=inx 5818 else 5819 inmax=inx 5820 end if 5821 end do 5822 inx=inmin 5823 5824 ! XG020629 : old coding, slower 5825 ! inx=0 5826 ! do while (rr >= rrad(inx+1,nn)) 5827 ! inx=inx+1 5828 ! end do 5829 5830 end if 5831 5832 ! Transformation matrix radial -> cart. coord 5833 ss=sqrt(vt1*vt1+vt2*vt2) 5834 ! if (ss /=0._dp) then ! XG020629 : never test a real number against zero (not portable) 5835 if (ss*ss > 1.0d-28) then ! ss non-zero 5836 ! XG020629 : very strange : only trsf(:,1) is needed in what follows ? ! 5837 ! ss_inv=1.0_dp/ss 5838 trsf(1,1)=vt1*rr_inv 5839 ! trsf(1,2)=-vt2*ss_inv 5840 ! trsf(1,3)=vt3*vt1*rr_inv*ss_inv 5841 trsf(2,1)=vt2*rr_inv 5842 ! trsf(2,2)=vt1*ss_inv 5843 ! trsf(2,3)=vt3*vt2*rr_inv*ss_inv 5844 trsf(3,1)=vt3*rr_inv 5845 ! trsf(3,2)=0._dp 5846 ! trsf(3,3)=-ss*rr_inv 5847 ! XG020629 Not needed 5848 ! do ii=1,3 5849 ! do ll=1,3 5850 ! ches(ii,ll)=0._dp 5851 ! end do 5852 ! cgrad(ii)=0._dp 5853 ! end do 5854 else ! ss zero 5855 do ii=1,3 5856 do ll=1,3 5857 trsf(ii,ll)=0._dp 5858 end do 5859 trsf(ii,4-ii)=1._dp 5860 end do 5861 end if ! ss zero or non-zero 5862 5863 if (inx == -1) then ! LEFT EXTRAPOLATION y=a*x^2+b (a<0) 5864 val=sp2(1,nn)*0.5_dp*rr*rr/rrad(1,nn)+crho(1,nn)-sp2(1,nn)*rrad(1,nn)*0.5_dp 5865 cgrad(1)=sp2(1,nn)*rr/rrad(1,nn) 5866 ches(1,1)=sp2(1,nn)/rrad(1,nn) 5867 elseif (inx == ndat(nn) ) then ! RIGHT EXTRAPOLATION y=a*exp(b*x) 5868 val=rrad(ndat(nn),nn)*exp(sp2(ndat(nn),nn)*(rr-rrad(ndat(nn),nn))/crho(ndat(nn),nn)) 5869 cgrad(1)=val*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 5870 ches(1,1)=cgrad(1)*sp2(ndat(nn),nn)/crho(ndat(nn),nn) 5871 else ! INTERPOLATION 5872 uu=rrad(inx+1,nn)-rrad(inx,nn) 5873 uu_inv=1.0_dp/uu 5874 aa=(rrad(inx+1,nn)-rr)*uu_inv 5875 bb=(rr-rrad(inx,nn))*uu_inv 5876 cc=(aa*aa*aa-aa)*uu*uu*0.16666666666666666_dp 5877 dd=(bb*bb*bb-bb)*uu*uu*0.16666666666666666_dp 5878 val=aa*crho(inx,nn)+bb*crho(inx+1,nn)+cc*sp3(inx,nn)+dd*sp3(inx+1,nn) 5879 cgrad(1)=(crho(inx+1,nn)-crho(inx,nn))*uu_inv& 5880 & -(3._dp*aa*aa-1._dp)*uu*0.16666666666666666_dp*sp3(inx,nn)+& 5881 & (3._dp*bb*bb-1._dp)*uu*0.16666666666666666_dp*sp3(inx+1,nn) 5882 ches(1,1)=aa*sp3(inx,nn)+bb*sp3(inx+1,nn) 5883 5884 end if ! TRANSFORMATION TO CARTEZ. COORD. 5885 5886 cgrad1_rr_inv=cgrad(1)*rr_inv 5887 coeff=(ches(1,1)-cgrad1_rr_inv)*rr_inv*rr_inv 5888 cgrad(3)=trsf(3,1)*cgrad(1) 5889 cgrad(2)=trsf(2,1)*cgrad(1) 5890 cgrad(1)=trsf(1,1)*cgrad(1) 5891 ches(1,1)=coeff*vt1*vt1+cgrad1_rr_inv 5892 ches(2,2)=coeff*vt2*vt2+cgrad1_rr_inv 5893 ches(3,3)=coeff*vt3*vt3+cgrad1_rr_inv 5894 ches(1,2)=coeff*vt1*vt2 ; ches(2,1)=coeff*vt1*vt2 5895 ches(1,3)=coeff*vt1*vt3 ; ches(3,1)=coeff*vt1*vt3 5896 ches(2,3)=coeff*vt2*vt3 ; ches(3,2)=coeff*vt2*vt3 5897 5898 else ! case rr==0 5899 5900 val=crho(1,nn)-sp2(1,nn)*rrad(1,nn)/2._dp 5901 do ii=1,3 5902 do ll=1,3 5903 ches(ii,ll)=0._dp 5904 end do 5905 cgrad(ii)=0._dp 5906 ches(ii,ii)=sp2(1,nn)/rrad(1,nn) 5907 end do 5908 5909 end if ! rr>0 or rr==0 5910 5911 do ii=1,3 5912 do ll=1,3 5913 hrho(ii,ll)=hrho(ii,ll)+ches(ii,ll) 5914 end do 5915 grho(ii)=grho(ii)+cgrad(ii) 5916 end do 5917 rho=rho+val 5918 5919 end if ! rr2< rrad_nn*rrad_nn 5920 5921 if (selct==-1) then 5922 if (rr2 < rminl(jj)*rminl(jj) ) then 5923 iat=jj 5924 ipos=kk 5925 rdmin=sqrt(rr2) 5926 end if 5927 elseif (selct==-2) then 5928 cycle 5929 else 5930 if (rr2 < rdmin*rdmin) then 5931 iat=jj 5932 ipos=kk 5933 rdmin=sqrt(rr2) 5934 end if 5935 end if 5936 5937 end do 5938 end do 5939 5940 end if 5941 5942 !write(51,'(":GRADt ",3F16.8)') grho 5943 !write(52,'(":LAPt ",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3) 5944 !write(52,'(":HESNt ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3) 5945 5946 !if(abs(cumul_cpu-cumul_cpu_old)>0.499)then 5947 !write(std_out,'(a,f7.1)' )' vgh_rho : cumul_cpu=',cumul_cpu 5948 !cumul_cpu_old=cumul_cpu 5949 !end if 5950 5951 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
5974 function vnorm(vv,dir) 5975 5976 5977 !This section has been created automatically by the script Abilint (TD). 5978 !Do not modify the following lines by hand. 5979 #undef ABI_FUNC 5980 #define ABI_FUNC 'vnorm' 5981 !End of the abilint section 5982 5983 implicit none 5984 5985 !Arguments ------------------------------------ 5986 !scalars 5987 integer,intent(in) :: dir 5988 real(dp) :: vnorm 5989 !arrays 5990 real(dp),intent(in) :: vv(3) 5991 5992 !Local variables------------------------------- 5993 !scalars 5994 integer :: ii 5995 !arrays 5996 real(dp) :: vt(3) 5997 5998 ! ************************************************************************* 5999 6000 vnorm=zero 6001 if (dir==1) then 6002 do ii=1,3 6003 vt(ii)=rprimd(ii,1)*vv(1)+rprimd(ii,2)*vv(2)+rprimd(ii,3)*vv(3) 6004 vnorm=vnorm+vt(ii)*vt(ii) 6005 end do 6006 elseif (dir==0) then 6007 do ii=1,3 6008 vnorm=vnorm+vv(ii)*vv(ii) 6009 end do 6010 else 6011 MSG_ERROR('vnorm calcul') 6012 end if 6013 vnorm=sqrt(vnorm) 6014 end function vnorm