TABLE OF CONTENTS
ABINIT/surf [ 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.
COPYRIGHT
Copyright (C) 2002-2017 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
43 #if defined HAVE_CONFIG_H 44 #include "config.h" 45 #endif 46 47 #include "abi_common.h" 48 49 50 subroutine surf(aim_dtset) 51 52 use defs_basis 53 use defs_aimprom 54 use defs_parameters 55 use defs_abitypes 56 use m_profiling_abi 57 use m_errors 58 use m_xmpi 59 60 use m_numeric_tools, only : coeffs_gausslegint 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'surf' 66 use interfaces_18_timing 67 use interfaces_63_bader, except_this_one => surf 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 type(aim_dataset_type) :: aim_dtset 75 76 !Local variables ------------------------------ 77 !scalars 78 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 79 real(dp) :: ct1,ct2,phi,rr,rsmax,rsmin,rthe,rthe0,t1,t2,theta,tt0,vcth,vph,vth 80 real(dp) :: wall,xy,xyz 81 logical :: srch,stemp 82 !arrays 83 real(dp) :: grho(3),vr(3),vv(3) 84 real(dp),allocatable :: rs_computed(:,:) 85 86 !************************************************************************ 87 88 comm = xmpi_world 89 me=xmpi_comm_rank(comm) 90 nproc=xmpi_comm_size(comm) 91 92 ttsrf=zero 93 94 rewind(unts) 95 96 nth=aim_dtset%nth 97 nph=aim_dtset%nph 98 99 !Coefficients for spherical Gauss quadrature 100 101 ct1=cos(aim_dtset%themin) 102 ct2=cos(aim_dtset%themax) 103 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 104 call coeffs_gausslegint(aim_dtset%phimin,aim_dtset%phimax,ph,wph,nph) 105 106 !DEBUG 107 !write(std_out,*)' surf : wcth=',wcth(1:nth) 108 !write(std_out,*)' surf : wph=',wph(1:nth) 109 !ENDDEBUG 110 111 do ijj=1,nth 112 th(ijj)=acos(cth(ijj)) 113 if (aim_dtset%isurf/=-1) then 114 do jj=1,nph 115 rs(ijj,jj)=zero 116 end do 117 end if 118 end do 119 120 npmax=aim_npmaxin 121 rsmax=0.0 122 rsmin=100.0 123 rthe0=r0 124 srch=.false. 125 126 do ijj=1,3 127 vv(ijj)=xatm(ijj,aim_dtset%batom) 128 end do 129 130 131 write(std_out,*) 132 write(std_out,*) "BADER SURFACE DETERMINATION" 133 write(std_out,*) "===========================" 134 write(std_out,*) 135 136 write(untout,*) 137 write(untout,*) "BADER SURFACE DETERMINATION" 138 write(untout,*) "===========================" 139 write(untout,*) 140 141 write(std_out,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 142 write(std_out,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 143 write(std_out,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 144 145 write(untout,'(" Atom: ",i3,3F15.10)') aim_dtset%batom,vv 146 write(untout,'(" Theta: ",i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 147 write(untout,'(" Phi: ",i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 148 149 write(unts,'(i3,3F15.10)') aim_dtset%batom,vv 150 write(unts,'(i3,2F15.10)') nth,aim_dtset%themin,aim_dtset%themax 151 write(unts,'(i3,2F15.10)') nph,aim_dtset%phimin,aim_dtset%phimax 152 153 !write(std_out,*) 'npmax in surf= ',npmax 154 155 ith=0 156 iph=0 157 tt0=0._dp 158 call timein(tt0,wall) 159 160 write(untout,*) 161 write(untout,*) "DEVELOPMENT OF THE RADII DETERMINATIONS" 162 write(untout,*) "========================================" 163 write(untout,*) 164 write(untout,*) "Determination near the CPs:" 165 166 !Determination of the CP neighbouring radii 167 168 if (aim_dtset%isurf/=-1) then 169 170 ! Precomputation of the value of the radii (for parallelisation) 171 ! To make the output independent of the number of processors, but still 172 ! cut down the CPU time, use a multigrid technique 173 srch=.true. 174 ABI_ALLOCATE(rs_computed,(nth,nph)) 175 rs(:,:)=zero 176 rs_computed(:,:)=zero 177 kk=0 ; init=0 178 do level=3,0,-1 179 incr=2**level 180 if(incr<nth .and. incr<nph)then 181 rs_computed(:,:)=rs(1:nth,1:nph) 182 rs(1:nth,1:nph)=zero 183 do ijj=1,nth,incr 184 do jj=1,nph,incr 185 if(rs_computed(ijj,jj)<1.0d-12) then 186 kk=kk+1 187 if(mod(kk,nproc)==me)then 188 ! Find an approximate starting radius, from the already computed ones 189 if(init==0)then 190 rthe=r0 191 else 192 ijj_exist=ijj ; if(mod(ijj-1,2*incr)>=incr)ijj_exist=ijj-incr 193 jj_exist=jj ; if(mod(jj-1,2*incr)>=incr)jj_exist=jj-incr 194 rthe=rs_computed(ijj_exist,jj_exist) 195 if(rthe<1.0d-12)then 196 write(std_out,*)' surf : there is a bug ! rthe=',rthe 197 MSG_ERROR("Aborting now") 198 end if 199 end if 200 call timein(t1,wall) ; t2=zero 201 call rsurf(aim_dtset,rr,grho,th(ijj),ph(jj),rthe,aim_dtset%batom,npmax,srch) 202 rs(ijj,jj)=rr 203 if (deb) then 204 call timein(t2,wall) ; t2=t2-t1 205 write(std_out,*) ':CALCULATED NP',ijj,jj,th(ijj),ph(jj),rthe,npmax,rs(ijj,jj),t2 206 end if 207 end if 208 end if 209 end do ! jj 210 end do ! ijj 211 call xmpi_sum(rs,comm,ierr) 212 ! Combine the set of already computed radii and the set of the newly computed, to obtain all computed. 213 rs(1:nth,1:nph)=rs(1:nth,1:nph)+rs_computed(:,:) 214 init=1 215 end if 216 end do 217 ABI_DEALLOCATE(rs_computed) 218 219 srch=.true. 220 221 do ijj=1,nbcp 222 ! if ((icpc(ijj) == -1)) then 223 rthe0=vnorm(pc(:,ijj),0) 224 do jj=1,3 225 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 226 end do 227 xy=vr(1)*vr(1)+vr(2)*vr(2) 228 xyz=xy+vr(3)*vr(3) 229 xyz=sqrt(xyz) 230 231 if (xy < aim_xymin) then 232 vcth=1._dp 233 if (vr(3) < 0._dp) vcth=-vcth 234 vph=0._dp 235 else 236 vcth=vr(3)/xyz 237 vph=atan2(vr(2),vr(1)) 238 end if 239 240 vth=acos(vcth) 241 write(untout,'(/," BCP: (index,theta,phi)",I4,2E16.8)') ijj,vth,vph 242 243 if (vth < th(1)) then 244 ith=0 245 else 246 if (vth > th(nth)) then 247 ith=nth 248 else 249 do ii=2,nth 250 if (vth < th(ii)) then 251 ith=ii-1 252 exit 253 end if 254 end do 255 end if 256 end if 257 258 if (vph < ph(1)) then 259 iph=0 260 else 261 if (vph > ph(nph)) then 262 iph=nph 263 else 264 do ii=2,nph 265 if (vph < ph(ii)) then 266 iph=ii-1 267 exit 268 end if 269 end do 270 end if 271 end if 272 273 write(untout,*) "ATOMIC RADII (ith,iph,theta,phi,radius)" 274 do jj=-1,2 275 do kk=-1,2 276 ith2=ith+jj 277 iph2=iph+kk 278 stemp=(iph2 > 0).and.(iph2 < nph+1) 279 stemp=(stemp.and.((ith2 > 0).and.(ith2 < nth+1))) 280 if (stemp) then 281 theta=th(ith2) 282 phi=ph(iph2) 283 if (abs(rs(ith2,iph2))<1.0d-12) then 284 rthe=rthe0 285 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 286 call timein(t1,wall) 287 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 288 call timein(t2,wall) 289 t2=t2-t1 290 rs(ith2,iph2)=rr 291 end if 292 rr=rs(ith2,iph2) 293 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 294 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 295 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 296 rthe0=rr 297 end if 298 299 end do ! kk 300 end do ! jj 301 302 ! end if 303 304 end do ! ijj (loop on BCP) 305 306 ! DEBUG 307 ! write(std_out,*)' surf : near BCP ' 308 ! do ijj=1,nth 309 ! do jj=1,nph 310 ! write(std_out,*)ijj,jj,rs(ijj,jj) 311 ! end do 312 ! end do 313 ! ENDDEBUG 314 315 316 srch=.true. 317 do ijj=nbcp+1,nbcp+nrcp ! Loop on RCP 318 ! if ((icpc(ijj) == 1)) then 319 rthe0=max(rminl(aim_dtset%batom),r0) 320 do jj=1,3 321 vr(jj)=pc(jj,ijj)-vv(jj)+xatm(jj,aim_dtset%batom) 322 end do 323 xy=vr(1)*vr(1)+vr(2)*vr(2) 324 xyz=xy+vr(3)*vr(3) 325 xyz=sqrt(xyz) 326 327 if (xy < aim_xymin) then 328 vcth=1._dp 329 if (vr(3) < 0._dp) vcth=-vcth 330 vph=0._dp 331 else 332 vcth=vr(3)/xyz 333 vph=atan2(vr(2),vr(1)) 334 end if 335 vth=acos(vcth) 336 write(untout,'(/,";RCP: (index,theta,phi)",I4,2E16.8)') ijj-nbcp,vth,vph 337 338 if (vth < th(1)) then 339 ith=0 340 else 341 if (vth > th(nth)) then 342 ith=nth 343 else 344 do ii=2,nth 345 if (vth < th(ii)) then 346 ith=ii-1 347 exit 348 end if 349 end do 350 end if 351 end if 352 353 if (vph < ph(1)) then 354 iph=0 355 else 356 if (vph > ph(nph)) then 357 iph=nph 358 else 359 do ii=2,nph 360 if (vph < ph(ii)) then 361 iph=ii-1 362 exit 363 end if 364 end do 365 end if 366 end if 367 368 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 369 do jj=-1,2 370 do kk=-1,2 371 ith2=ith+jj 372 iph2=iph+kk 373 stemp=(iph2 > 0).and.(iph2 < nph+1) 374 stemp=stemp.and.(ith2 > 0).and.(ith2 < nth+1) 375 376 if (stemp) then 377 theta=th(ith2) 378 phi=ph(iph2) 379 if ((abs(rs(ith2,iph2))<1.0d-12)) then 380 rthe=rthe0 381 if (deb) write(std_out,*) ':CALCULATING NP',theta,phi,rthe,npmax 382 call timein(t1,wall) 383 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 384 call timein(t2,wall) 385 t2=t2-t1 386 rs(ith2,iph2)=rr 387 end if 388 rr=rs(ith2,iph2) 389 ! write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 390 write(std_out,'(":RSUR PC ",3i3,4E16.8,F10.4)') ijj,jj,kk,theta,phi,rr,wcth(ith2)*wph(iph2),t2 391 write(untout,'(a,2i3,3E16.8)') '- ',jj,kk,theta,phi,rr 392 rthe0=rr 393 end if 394 395 end do ! kk 396 end do ! jj 397 ! end if 398 399 end do ! ijj (Loop on RCP) 400 401 ! DEBUG 402 ! write(std_out,*)' surf : near RCP ' 403 ! do ijj=1,nth 404 ! do jj=1,nph 405 ! write(std_out,*)ijj,jj,rs(ijj,jj) 406 ! end do 407 ! end do 408 ! ENDDEBUG 409 410 ! Boundary angles 411 rthe0=r0 412 srch=.true. 413 write(untout,*) 414 write(untout,*) "The boundary angles:" 415 write(untout,*) "====================" 416 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 417 418 ! Must have sufficient angular sampling 419 if ((nth > 8).and.(nph > 8)) then 420 rthe=r0 421 do ijj=1,2 422 theta=th(ijj) 423 if (ijj==2) rthe=rs(1,1) 424 do jj=1,nph 425 phi=ph(jj) 426 call timein(t1,wall) 427 if (abs(rs(ijj,jj))<1.0d-12) then 428 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 429 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 430 rs(ijj,jj)=rr 431 end if 432 rr=rs(ijj,jj) 433 call timein(t2,wall) 434 t2=t2-t1 435 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 436 write(untout,'(a,2i3,3E16.8)') '- ',ijj,jj,theta,phi,rr 437 rthe=rs(ijj,jj) 438 end do ! jj 439 end do ! ijj 440 441 write(untout,*) 442 443 rthe=rs(2,1) 444 do jj=1,2 445 phi=ph(jj) 446 if (jj==2) rthe=rs(2,2) 447 do ijj=3,nth 448 theta=th(ijj) 449 t2=0.0 450 call timein(t1,wall) 451 if (abs(rs(ijj,jj))<1.0d-12) then 452 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 453 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 454 rs(ijj,jj)=rr 455 end if 456 rr=rs(ijj,jj) 457 call timein(t2,wall) 458 t2=t2-t1 459 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 460 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 461 rthe=rs(ijj,jj) 462 end do ! ijj 463 end do ! jj 464 465 write(untout,*) 466 467 rthe=rs(nth-1,2) 468 do ijj=nth-1,nth 469 theta=th(ijj) 470 if (ijj==nth) rthe=rs(nth,2) 471 do jj=3,nph 472 phi=ph(jj) 473 call timein(t1,wall) 474 if (abs(rs(ijj,jj))<1.0d-12) then 475 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 476 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 477 rs(ijj,jj)=rr 478 end if 479 rr=rs(ijj,jj) 480 call timein(t2,wall) 481 t2=t2-t1 482 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 483 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 484 rthe=rs(ijj,jj) 485 end do ! jj 486 end do ! ijj 487 488 rthe=rs(2,nph-1) 489 do jj=nph-1,nph 490 phi=ph(jj) 491 if (jj==nph) rthe=rs(2,nph) 492 do ijj=3,nth-2 493 theta=th(ijj) 494 t2=0.0 495 call timein(t1,wall) 496 if (abs(rs(ijj,jj))<1.0d-12) then 497 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 498 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 499 rs(ijj,jj)=rr 500 end if 501 rr=rs(ijj,jj) 502 call timein(t2,wall) 503 t2=t2-t1 504 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 505 write(untout,'(2i3,3E16.8)') ijj,jj,theta,phi,rr 506 rthe=rs(ijj,jj) 507 end do ! ijj 508 end do ! jj 509 write(untout,*) 510 511 ! Complementary bands for boundary angles 512 nn=int(real(nth)/1.4d1) 513 if (nn > 1) then 514 do ii=1,nn-1 515 mm=int(nth/nn)*ii 516 do kk=0,1 517 mm=mm+kk 518 theta=th(mm) 519 rthe=rs(mm,2) 520 do jj=3,nph-2 521 phi=ph(jj) 522 call timein(t1,wall) 523 if (abs(rs(mm,jj))<1.0d-12) then 524 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 525 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 526 rs(mm,jj)=rr 527 end if 528 rr=rs(mm,jj) 529 call timein(t2,wall) 530 t2=t2-t1 531 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(mm)*wph(jj),t2 532 write(untout,'(2i3,3E16.8)') mm,jj,theta,phi,rr 533 rthe=rs(mm,jj) 534 end do ! jj 535 end do ! kk 536 end do ! ii 537 end if ! nn>1 538 539 write(untout,*) 540 541 nn=nint(real(nph)/1.2d1) 542 if (nn > 1) then 543 do ii=1,nn-1 544 mm=int(nph/nn)*ii 545 do kk=0,1 546 mm=mm+kk 547 phi=ph(mm) 548 rthe=rs(2,mm) 549 550 do jj=3,nth-2 551 theta=th(jj) 552 call timein(t1,wall) 553 if (abs(rs(jj,mm))<1.0d-12) then 554 if (deb) write(std_out,*) ':CALC NP',theta,phi,rthe,npmax 555 call rsurf(aim_dtset,rr,grho,theta,phi,rthe,aim_dtset%batom,npmax,srch) 556 rs(jj,mm)=rr 557 end if 558 rr=rs(mm,jj) 559 call timein(t2,wall) 560 t2=t2-t1 561 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(jj)*wph(mm),t2 562 write(untout,'(2i3,3E16.8)') jj,mm,theta,phi,rr 563 rthe=rs(jj,mm) 564 end do ! jj 565 566 end do ! kk 567 end do ! ii 568 end if ! nn>1 569 570 end if ! sufficient sampling to determine boundary angles 571 572 write(untout,*) 573 574 ! DEBUG 575 ! write(std_out,*)' surf : after boundary angles ' 576 ! do ijj=1,nth 577 ! do jj=1,nph 578 ! write(std_out,*)ijj,jj,rs(ijj,jj) 579 ! end do 580 ! end do 581 ! ENDDEBUG 582 583 ! Output the complete Bader surface 584 585 write(untout,*) "The complete Bader surface:" 586 write(untout,*) "===========================" 587 write(untout,*) "ATOMIC RADIUS (ith,iph,theta,phi,radius)" 588 rthe0=r0 589 srch=.true. 590 591 ! Write all the values 592 593 do ijj=1,nth 594 theta=th(ijj) 595 do jj=1,nph 596 phi=ph(jj) 597 rr=rs(ijj,jj) 598 write(unts,'(2F12.8,2E16.8)') theta,phi,rr,wcth(ijj)*wph(jj) 599 write(std_out,'(":RSUR ",2F12.8,2E16.8,F10.4)') theta,phi,rr,wcth(ijj)*wph(jj),t2 600 write(untout,'(a,2i3,3E16.8)') ' ',ijj,jj,theta,phi,rr 601 if (rr < rsmin) rsmin=rr 602 if (rr> rsmax) rsmax=rr 603 end do ! jj 604 end do ! ijj 605 write(unts,'(2F15.10)') rsmin,rsmax 606 write(untout,'(/," The minimal and maximal radii:",/,/," ",2F15.10)') rsmin,rsmax 607 608 ! DEBUG 609 ! write(std_out,*)' surf : final output ' 610 ! do ijj=1,nth 611 ! do jj=1,nph 612 ! write(std_out,*)ijj,jj,rs(ijj,jj) 613 ! end do 614 ! end do 615 ! ENDDEBUG 616 617 end if ! determination of the critical surface 618 619 call timein(ttsrf,wall) 620 ttsrf=ttsrf-tt0 621 622 end subroutine surf