TABLE OF CONTENTS
ABINIT/cpdrv [ 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.
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
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
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 subroutine cpdrv(aim_dtset) 49 50 use defs_basis 51 use defs_aimprom 52 use defs_parameters 53 use defs_datatypes 54 use defs_abitypes 55 use m_xmpi 56 use m_errors 57 use m_sort 58 use m_profiling_abi 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'cpdrv' 64 use interfaces_18_timing 65 use interfaces_63_bader, except_this_one => cpdrv 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 type(aim_dataset_type),intent(in) :: aim_dtset 73 74 !Local variables ------------------------------ 75 !scalars 76 integer :: iat,iatinit,ii,inxat,inxcell,ipair,ipos,iposinit,ires,jj,kk,nb,nb_now 77 integer :: nn,nstep,nvs,me,nproc,ierr 78 real(dp) :: candidate,chg,diff1,diff2,diff3,dist,prj,rtdiff,ss,tt0,wall 79 logical :: srch=.false. 80 !arrays 81 integer :: ibat(nnpos*natom),inatm(nnpos*natom),incell(nnpos*natom) 82 integer :: ipibat(nnpos*natom) 83 integer,allocatable :: indexcp(:),nr(:) 84 real(dp) :: bmin(natom),dif(3),dists(nnpos*natom),ev(3),evec(3,3),grho(3) 85 real(dp) :: hrho(3,3),pom(3),rr(3),uu(3),vv(3),xorig(3) 86 real(dp),allocatable :: buffer(:,:),sortguide(:) 87 !no_abirules 88 !Warning : bcp_type should be transformed to cp_type 89 type(bcp_type),allocatable :: bcp(:),ccp(:),cp_tmp(:),rcp(:) 90 91 !************************************************************************ 92 93 me=xmpi_comm_rank(xmpi_world) 94 nproc=xmpi_comm_size(xmpi_world) 95 96 !Consider the critical points starting from atom #batom 97 inxat=aim_dtset%batom 98 slc=-1 99 rminl(:)=aim_dtset%rmin 100 bmin(:)=0._dp 101 ttcp=0._dp 102 103 write(std_out,*) 104 write(std_out,*) "CRITICAL POINTS ANALYSIS" 105 write(std_out,*) "========================" 106 write(std_out,*) 107 108 write(untout,*) 109 write(untout,*) "CRITICAL POINTS ANALYSIS" 110 write(untout,*) "========================" 111 write(untout,*) 112 113 114 xorig(:)=xatm(:,inxat) 115 116 call timein(tt0,wall) 117 118 !Searching the neighbouring atoms 119 120 if (aim_dtset%crit > 0) then 121 nvs=0 122 do ii=1,nnpos 123 do jj=1,natom 124 dist=0._dp 125 dif(:)=xatm(:,inxat)-xatm(:,jj)-atp(:,ii) 126 dif(:)=dif(:)/aim_dtset%scal(:) 127 dist=vnorm(dif,0) 128 if (dist < tol6 ) then 129 inxcell=ii 130 elseif (dist < maxatdst) then 131 nvs=nvs+1 132 dists(nvs)=dist 133 inatm(nvs)=jj 134 incell(nvs)=ii 135 end if 136 end do 137 end do 138 139 write(std_out,*) "ATOM:" 140 write(std_out,*) 'inxat :', inxat, 'inxcell :', inxcell 141 write(std_out, '(3es16.6)' ) (xorig(ii),ii=1,3) 142 write(std_out,*) 143 144 write(untout,*) "ATOM:" 145 write(untout,*) 'inxat :', inxat, 'inxcell :', inxcell 146 write(untout, '(3es16.6)') (xorig(ii),ii=1,3) 147 write(untout,*) 148 149 ABI_ALLOCATE(nr,(nvs)) 150 do ii=1,nvs 151 nr(ii)=ii 152 end do 153 154 ! Ordering of the nearest neighbouring atoms 155 call sort_dp(nvs,dists,nr,tol14) 156 157 nb=0 158 write(std_out,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 159 write(untout,*) "NEIGHBORING ATOMS (atindex,cellindex,distance(in bohr)):" 160 do ii=1,nvs 161 nn=nr(ii) 162 if (dists(ii) < maxatdst) then 163 nb=nb+1 164 ibat(nb)=inatm(nn) 165 ipibat(nb)=incell(nn) 166 write(std_out,*) ':NEIG ',inatm(nn),incell(nn),dists(ii) 167 write(untout,'(" ",2I6,F16.8)')inatm(nn),incell(nn),dists(ii) 168 else 169 exit 170 end if 171 end do 172 173 ! SEARCHING BCP 174 ABI_DATATYPE_ALLOCATE(bcp,(nb)) 175 nbcp=0 176 iatinit=inxat 177 iposinit=inxcell 178 bcp(:)%iat=0 179 bcp(:)%ipos=0 180 181 write(std_out,*) 182 write(std_out,*) "BONDING CRITICAL POINTS (BCP)" 183 write(std_out,*) "=============================" 184 write(std_out,*) 185 186 write(untout,*) 187 write(untout,*) "BONDING CRITICAL POINTS (BCP)" 188 write(untout,*) "=============================" 189 write(untout,*) 190 191 srbcp: do ii=1,nb 192 193 ! Start the search for BCP from the midistance between the atom 194 ! and his neighbor. 195 vv(:)=(xatm(:,inxat)+xatm(:,ibat(ii))+atp(:,ipibat(ii)))/2._dp 196 197 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,-1) 198 199 if (ires==0) then 200 ! Testing if CP is already known 201 if (nbcp > 0) then 202 do jj=1,nbcp 203 pom(:)=vv(:)-bcp(jj)%rr(:)-xorig(:) 204 dist=vnorm(pom,0) 205 if (dist < aim_dtset%dpclim) then 206 write(std_out,*) 'BCP already known !' 207 cycle srbcp 208 end if 209 end do 210 end if 211 rr(:)=vv(:)-xorig(:) 212 ss=vnorm(rr,0) 213 if (ss > maxcpdst) then 214 write(std_out, '(a,es16.6,a,es16.6)' ) 'BCP distance from atom,',ss,', exceed maxcpdst =',maxcpdst 215 cycle srbcp 216 end if 217 nn=0 218 do jj=1,3 219 nn=nn+ev(jj)/abs(ev(jj)) 220 end do 221 write(std_out, '(a,3es16.6,i4)') ' vv(1:3), nn',(vv(jj), jj=1,3), nn 222 write(std_out, '(a,3es16.6)') 'ev: ', (ev(jj), jj=1,3) 223 if (nn /= -1) then 224 write(std_out,*) ' The trial critical point is not a BCP !' 225 cycle srbcp 226 end if 227 write(std_out, '(a,3es16.6)' ) 'evec(:,1): ',(evec(jj,1), jj=1,3) 228 pom(:)=evec(:,1) 229 dist=vnorm(pom,0) 230 prj=dot_product(evec(:,1),rr) 231 write(std_out,*) 'prj:', prj, vnorm(evec(:,1),0) 232 dist=vnorm(evec(:,1),0) 233 uu(:)=vv(:)-sign(aim_epstep,prj)*evec(:,1)/dist 234 235 ! Testing whether this BCP "is bonded" to the considered atom 236 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 237 ! write(std_out,*) 'do', iat, ipos 238 ! if ((iat==0).or.(ipos==0)) cycle 239 ! write(std_out,*) 'APOS: ',(xatm(jj,iat)+atp(jj,ipos), jj=1,3) 240 if ((iat/=inxat).or.(inxcell/=ipos)) then 241 write(std_out,*) ' The trial BCP is not bonded to the Bader atom' 242 cycle srbcp 243 end if 244 245 ! A new BCP has been found ! 246 nbcp=nbcp+1 247 248 ! Searching for the second bonded atom 249 ss=vnorm(rr,0) 250 diff1=ss 251 diff3=dists(ii) 252 uu(:)=vv(:)+sign(aim_epstep,prj)*evec(:,1)/dist 253 if ((abs(bmin(iat))<1.0d-12).or.( ss<bmin(iat))) then 254 bmin(iat)=ss 255 end if 256 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 257 if ((iat==0).or.(ipos==0)) then 258 write(std_out,*) ' The trial BCP is not bonded to a bonding atom !' 259 ! cycle srbcp 260 end if 261 pom(:)=vv(:)-xatm(:,iat)-atp(:,ipos) 262 ss=vnorm(pom,0) 263 diff2=ss 264 pom(:)=xorig(:)-xatm(:,iat)-atp(:,ipos) 265 diff3=vnorm(pom,0) 266 rtdiff=diff1/diff3 267 if ((abs(bmin(iat))<1.0d-12).or.(ss<bmin(iat))) then 268 bmin(iat)=ss 269 end if 270 pom(:)=xatm(:,iat)+atp(:,ipos) 271 272 ! Store more results, for coherent, and portable output 273 bcp(nbcp)%iat=iat 274 bcp(nbcp)%ipos=ipos 275 bcp(nbcp)%chg=chg 276 bcp(nbcp)%diff(1)=diff1 277 bcp(nbcp)%diff(2)=diff2 278 bcp(nbcp)%diff(3)=diff3 279 bcp(nbcp)%ev(:)=ev(:) 280 bcp(nbcp)%pom(:)=pom(:) 281 bcp(nbcp)%rr(:)=rr(:) 282 bcp(nbcp)%vec(:,:)=evec(:,:) 283 bcp(nbcp)%vv(:)=vv(:) 284 ! Warning : iat, ipos might be modified by this call 285 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 286 bcp(nbcp)%chg=chg 287 288 end if ! ires==0 289 end do srbcp 290 291 if(nbcp>0)then 292 293 ! Order the BCP. CPs should appear by increasing values of x,y,z , the latter 294 ! varying the fastest 295 ABI_ALLOCATE(sortguide,(nbcp)) 296 ABI_ALLOCATE(indexcp,(nbcp)) 297 ABI_DATATYPE_ALLOCATE(cp_tmp,(nbcp)) 298 do ii=3,1,-1 299 ! DEBUG 300 ! write(std_out,*)' cpdrv : sort on index ii=',ii 301 ! ENDDEBUG 302 303 do jj=1,nbcp 304 ! DEBUG 305 ! write(std_out,*)bcp(jj)%vv(:) 306 ! ENDDEBUG 307 sortguide(jj)=bcp(jj)%vv(ii) 308 indexcp(jj)=jj 309 end do 310 311 ! Try to be platform-independent. Might need a larger tolerance. 312 call sort_dp(nbcp,sortguide,indexcp,tol3) 313 do jj=1,nbcp 314 cp_tmp(jj)=bcp(indexcp(jj)) 315 end do 316 do jj=1,nbcp 317 bcp(jj)=cp_tmp(jj) 318 end do 319 end do 320 ! DEBUG 321 ! write(std_out,*)' cpdrv : after the sort ' 322 ! do jj=1,nbcp 323 ! write(std_out,*)bcp(jj)%vv(:) 324 ! end do 325 ! ENDDEBUG 326 327 328 ! Output the info about the BCP 329 do jj=1,nbcp 330 write(untout,'(" Bonded atom (BAT) (indxatm,indxcell,position): ",/,2I6,3F16.8)')& 331 & bcp(jj)%iat,bcp(jj)%ipos,bcp(jj)%pom(:) 332 write(untout,'("%Bonding CP: ",3F16.8)') bcp(jj)%vv(:) 333 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') bcp(jj)%ev(:) 334 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 335 & ' Eigenvec. of Hessian:',char(10),& 336 & '-',bcp(jj)%vec(1,:),char(10),& 337 & '-',bcp(jj)%vec(2,:),char(10),& 338 & '-',bcp(jj)%vec(3,:),char(10) 339 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 340 & bcp(jj)%chg, bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 341 write(untout,'("%Relative position of BCP (AT-CP,BAT-CP,AT-BAT,relative(AT): ",/,4F16.8)') & 342 & bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 343 write(untout,*) "********************************************************************" 344 write(std_out,'(/," BCP: ",3F10.6,3E12.4,E12.4,/)') & 345 & bcp(jj)%rr(:),bcp(jj)%ev(:),bcp(jj)%ev(1)+bcp(jj)%ev(2)+bcp(jj)%ev(3) 346 write(std_out,'(":DISPC ",4F12.6)') bcp(jj)%diff(:),bcp(jj)%diff(1)/bcp(jj)%diff(3) 347 end do 348 349 ABI_DATATYPE_DEALLOCATE(cp_tmp) 350 ABI_DEALLOCATE(indexcp) 351 ABI_DEALLOCATE(sortguide) 352 353 end if ! nbcp>0 354 355 if (abs(bmin(inxat))>1.0d-12) then 356 rminl(inxat)=aim_dtset%coff1*bmin(inxat) 357 r0=bmin(inxat) 358 else 359 r0=0._dp 360 end if 361 362 ! !AD-HOC PARAMETER 363 364 do ii=1,natom 365 if ((abs(bmin(ii))>1.0d-12).and.(ii /= inxat)) rminl(ii)=aim_dtset%coff2*bmin(ii) 366 end do 367 368 ! END WARNING 369 370 write(std_out,*) ' number of BCP:', nbcp 371 write(untout,'(" Number of BCP found: ",I4)') nbcp 372 nn=nbcp*(nbcp-1)*(nbcp-2)/6 373 if (bit_size(ii) <= nbcp+1) then 374 MSG_ERROR("b-test!") 375 end if 376 377 ! SEARCHING RCP 378 379 write(std_out,*) 380 write(std_out,*) "RING CRITICAL POINTS (RCP)" 381 write(std_out,*) "=============================" 382 write(std_out,*) 383 384 write(untout,*) 385 write(untout,*) "RING CRITICAL POINTS (RCP)" 386 write(untout,*) "=============================" 387 write(untout,*) 388 389 nrcp=0 390 if(aim_dtset%crit==1)nb_now=nbcp 391 if(aim_dtset%crit==2)nb_now=nb 392 ! DEBUG 393 ! nb_now=nbcp 394 ! ENDDEBUG 395 nn=nb_now*(nb_now-1)/2 396 ABI_DATATYPE_ALLOCATE(rcp,(nn)) 397 398 ! Loop on pairs of BCP or atoms 399 ipair=0 400 ABI_ALLOCATE(buffer,(16,nn)) 401 buffer=zero 402 403 ! DEBUG 404 ! write(std_out,*)ch10,ch10,' drvcpr : enter loop to search for RCPs,nb_now,nn=',nb_now,nn 405 ! ENDDEBUG 406 407 do ii=1,nb_now-1 408 srcp1: do jj=ii+1,nb_now 409 ipair=ipair+1 410 if(mod(ipair,nproc)==me)then 411 if (aim_dtset%crit==1) then 412 vv(:)=xorig(:)+(bcp(ii)%rr(:)+bcp(jj)%rr(:))/2._dp 413 else if (aim_dtset%crit==2) then 414 vv(:)=xorig(:)*half+(xatm(:,ibat(ii))+atp(:,ipibat(ii))+xatm(:,ibat(jj))+atp(:,ipibat(jj)))*quarter 415 end if 416 417 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,1) 418 419 if(ires==1)then 420 cycle srcp1 421 end if 422 423 ! Check that it is within the maximum allowed distance for a CP 424 rr(:)=vv(:)-xorig(:) 425 ss=vnorm(rr,0) 426 if (ss > maxcpdst) then 427 write(std_out,*) 'RCP distance from atom exceed maxcpdst !' 428 cycle srcp1 429 end if 430 ! Check that it is a RCP 431 nn=0 432 do kk=1,3 433 nn=nn+ev(kk)/abs(ev(kk)) 434 end do 435 if (nn /= 1) then 436 write(std_out,*) ' the critical point that is found is not a RCP ' 437 cycle srcp1 438 end if 439 ! Might be the same RCP than one already found on the same processor 440 if (nrcp > 0) then 441 do kk=1,nrcp 442 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 443 dist=vnorm(pom,0) 444 if (dist < aim_dtset%dpclim) then 445 write(std_out,*) ':RCP already known' 446 cycle srcp1 447 end if 448 end do 449 end if 450 ! If crit==2, check that it is on the Bader surface 451 if (aim_dtset%crit==2) then 452 uu(:)=vv(:)-aim_epstep*rr(:)/ss 453 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 454 if ((iat/=inxat).or.(inxcell/=ipos))then 455 write(std_out,*) ' RCP is not on the Bader surface (outside of it)' 456 cycle srcp1 457 end if 458 end if 459 nrcp=nrcp+1 460 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 461 462 buffer(1:3,ipair)=vv 463 buffer(4:6,ipair)=ev 464 buffer(7:9,ipair)=evec(:,1) 465 buffer(10:12,ipair)=evec(:,2) 466 buffer(13:15,ipair)=evec(:,3) 467 buffer(16,ipair)=one 468 469 ! DEBUG 470 ! write(std_out,*)ch10,ch10,' drvcpr : ipair,candidate=',ipair,candidate 471 ! ENDDEBUG 472 end if 473 end do srcp1 474 end do 475 call xmpi_sum(buffer,xmpi_world,ierr) 476 477 nrcp=0 478 ipair=0 479 do ii=1,nb_now-1 480 srcp: do jj=ii+1,nb_now 481 ipair=ipair+1 482 candidate=buffer(16,ipair) 483 484 ! One CP has been found, must make tests to see whether it is a new RCP 485 if (nint(candidate)==1) then 486 487 vv=buffer(1:3,ipair) 488 ev=buffer(4:6,ipair) 489 evec(:,1)=buffer(7:9,ipair) 490 evec(:,2)=buffer(10:12,ipair) 491 evec(:,3)=buffer(13:15,ipair) 492 493 ! Check that it is not the same as a previous one 494 if (nrcp > 0) then 495 do kk=1,nrcp 496 pom(:)=vv(:)-rcp(kk)%rr(:)-xorig(:) 497 dist=vnorm(pom,0) 498 if (dist < aim_dtset%dpclim) then 499 write(std_out,*) ':RCP already known' 500 cycle srcp 501 end if 502 end do 503 end if 504 505 ! A new RCP has been found ! 506 nrcp=nrcp+1 507 508 ! DEBUG 509 ! write(std_out,*)' drvcpr : A new RCP has been found, for kk=',kk 510 ! ENDDEBUG 511 512 513 rcp(nrcp)%iat=iat 514 rcp(nrcp)%ipos=ipos 515 rcp(nrcp)%rr(:)=vv(:)-xorig(:) 516 rcp(nrcp)%vec(:,:)=evec(:,:) 517 rcp(nrcp)%ev(:)=ev(:) 518 rcp(nrcp)%vv(:)=vv(:) 519 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 520 rcp(nrcp)%chg=chg 521 522 end if ! ires==0 523 end do srcp ! jj=ii+2,nb_now 524 end do ! ii=1,nb_now-1 525 526 ABI_DEALLOCATE(buffer) 527 528 if(nrcp>0)then 529 530 ! Order the RCP. CPs should appear by increasing values of x,y,z , the latter 531 ! varying the fastest 532 ABI_ALLOCATE(sortguide,(nrcp)) 533 ABI_ALLOCATE(indexcp,(nrcp)) 534 ABI_DATATYPE_ALLOCATE(cp_tmp,(nrcp)) 535 do ii=3,1,-1 536 ! DEBUG 537 ! write(std_out,*)' cpdrv : sort on index ii=',ii 538 ! ENDDEBUG 539 do jj=1,nrcp 540 541 ! DEBUG 542 ! write(std_out,*)rcp(jj)%vv(:) 543 ! ENDDEBUG 544 545 ! Try to be platform-independent. Might need a larger tolerance. 546 sortguide(jj)=rcp(jj)%vv(ii) 547 indexcp(jj)=jj 548 end do 549 call sort_dp(nrcp,sortguide,indexcp,tol3) 550 do jj=1,nrcp 551 cp_tmp(jj)=rcp(indexcp(jj)) 552 end do 553 do jj=1,nrcp 554 rcp(jj)=cp_tmp(jj) 555 end do 556 end do 557 558 ! DEBUG 559 ! write(std_out,*)' cpdrv : after the sort ' 560 ! do jj=1,nrcp 561 ! write(std_out,*)rcp(jj)%vv(:) 562 ! end do 563 ! ENDDEBUG 564 565 566 ! Write the Ring Critical Point information 567 do jj=1,nrcp 568 write(untout,'(";Ring CP: ",3F16.8)') rcp(jj)%vv(:) 569 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') rcp(jj)%ev(:) 570 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 571 & ' Eigenvec. of Hessian:',char(10),& 572 & '-',rcp(jj)%vec(1,:),char(10),& 573 & '-',rcp(jj)%vec(2,:),char(10),& 574 & '-',rcp(jj)%vec(3,:),char(10) 575 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 576 & rcp(jj)%chg, rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 577 write(untout,*) "********************************************************************" 578 write(std_out,'(/," RCP: ",3F10.6,3E12.4,E12.4,/)') & 579 & rcp(jj)%rr(:),rcp(jj)%ev(:),rcp(jj)%ev(1)+rcp(jj)%ev(2)+rcp(jj)%ev(3) 580 end do 581 582 ABI_DATATYPE_DEALLOCATE(cp_tmp) 583 ABI_DEALLOCATE(indexcp) 584 ABI_DEALLOCATE(sortguide) 585 586 end if ! nrcp>0 587 588 write(untout,'(" Number of RCP found: ",I4)') nrcp 589 write(std_out,*) ' Number of RCP:', nrcp 590 591 ! SEARCHING CCP 592 593 write(std_out,*) 594 write(std_out,*) "CAGE CRITICAL POINTS (CCP)" 595 write(std_out,*) "=============================" 596 write(std_out,*) 597 598 write(untout,*) 599 write(untout,*) "CAGE CRITICAL POINTS (CCP)" 600 write(untout,*) "=============================" 601 write(untout,*) 602 603 604 nn=nrcp*(nrcp-1)/2 605 ABI_DATATYPE_ALLOCATE(ccp,(nn)) 606 607 nccp=0 608 do ii=1,nrcp-1 609 srccp: do jj=ii+1,nrcp 610 vv(:)=xorig(:)+(rcp(ii)%rr(:)+rcp(jj)%rr(:))/2._dp 611 call critic(aim_dtset,vv,ev,evec,aim_dmaxcs,ires,3) 612 if (ires==0) then 613 rr(:)=vv(:)-xorig(:) 614 ss=vnorm(rr,0) 615 if (ss > maxcpdst) then 616 write(std_out,*) 'CCP distance from atom exceed maxcpdst !' 617 cycle srccp 618 end if 619 nn=0 620 do kk=1,3 621 nn=nn+ev(kk)/abs(ev(kk)) 622 end do 623 if (nn /= 3) then 624 write(std_out,*) ' the critical point that is found is not a CCP ' 625 cycle srccp 626 end if 627 628 if (nccp > 0) then 629 do kk=1,nccp 630 pom(:)=vv(:)-ccp(kk)%rr(:)-xorig(:) 631 dist=vnorm(pom,0) 632 if (dist < aim_dtset%dpclim) then 633 write(std_out,*) ':CCP already known' 634 cycle srccp 635 end if 636 end do 637 end if 638 if (aim_dtset%crit==2) then 639 uu(:)=vv(:)-aim_epstep*rr(:)/ss 640 call aim_follow(aim_dtset,uu,aim_npmaxin,srch,iatinit,iposinit,iat,ipos,nstep) 641 if ((iat/=inxat).or.(inxcell/=ipos)) then 642 write(std_out,*) ' This CCP is not on the Bader surface (outside of it)' 643 cycle srccp 644 end if 645 end if 646 647 nccp=nccp+1 648 649 ccp(nccp)%iat=iat 650 ccp(nccp)%ipos=ipos 651 ccp(nccp)%rr(:)=vv(:)-xorig(:) 652 ccp(nccp)%vec(:,:)=evec(:,:) 653 ccp(nccp)%ev(:)=ev(:) 654 ccp(nccp)%vv(:)=vv(:) 655 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 656 ccp(nccp)%chg=chg 657 658 end if 659 end do srccp 660 end do 661 662 if(nccp>0)then 663 664 ! Order the CCP. CPs should appear by increasing values of x,y,z , the latter 665 ! varying the fastest 666 ABI_ALLOCATE(sortguide,(nccp)) 667 ABI_ALLOCATE(indexcp,(nccp)) 668 ABI_DATATYPE_ALLOCATE(cp_tmp,(nccp)) 669 do ii=3,1,-1 670 do jj=1,nccp 671 ! Try to be platform-independent. Might need a larger tolerance. 672 sortguide(jj)=ccp(jj)%vv(ii) 673 indexcp(jj)=jj 674 end do 675 call sort_dp(nccp,sortguide,indexcp,tol3) 676 do jj=1,nccp 677 cp_tmp(jj)=ccp(indexcp(jj)) 678 end do 679 do jj=1,nccp 680 ccp(jj)=cp_tmp(jj) 681 end do 682 end do 683 684 ! Write the Cage Critical Point information 685 do jj=1,nccp 686 write(untout,'("%Cage CP: ",3F16.8)') ccp(jj)%vv(:) 687 write(untout,'("%Eigenval. of Hessian: ",3F16.8)') ccp(jj)%ev(:) 688 write(untout,'(a,a,a,3f16.8,a,a,3f16.8,a,a,3f16.8,a)') & 689 & ' Eigenvec. of Hessian:',char(10),& 690 & '-',ccp(jj)%vec(1,:),char(10),& 691 & '-',ccp(jj)%vec(2,:),char(10),& 692 & '-',ccp(jj)%vec(3,:),char(10) 693 write(untout,'("%Density and laplacian in CP: ",2F16.8)') & 694 & ccp(jj)%chg, ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 695 write(untout,*) "********************************************************************" 696 write(std_out,'(/," CCP: ",3F10.6,3E12.4,E12.4,/)') & 697 & ccp(jj)%rr(:),ccp(jj)%ev(:),ccp(jj)%ev(1)+ccp(jj)%ev(2)+ccp(jj)%ev(3) 698 end do 699 700 ABI_DEALLOCATE(sortguide) 701 ABI_DEALLOCATE(indexcp) 702 ABI_DATATYPE_DEALLOCATE(cp_tmp) 703 704 end if ! nccp>0 705 706 write(untout,'(" Number of CCP found: ",I4)') nccp 707 write(std_out,*) 'Number of CCP:', nccp 708 write(std_out,*) 709 write(untout,*) 710 write(std_out, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 711 write(untout, '(a,3i8)' ) 'BCP-RCP-CCP', nbcp,nrcp,nccp 712 713 write(std_out,*) 714 write(std_out,*) "===============================" 715 write(std_out,*) "END OF CRITICAL POINTS ANALYSIS" 716 write(std_out,*) 717 718 write(untout,*) 719 write(untout,*) "===============================" 720 write(untout,*) "END OF CRITICAL POINTS ANALYSIS" 721 write(untout,*) 722 723 724 ! Output of the CPs 725 726 write(untc,'(I4, " :BCP''s, coordinates, laplacian eigs, type of bonding at., sum of lap.eigs., density")') nbcp 727 do ii=1,nbcp 728 write(untc,'(3F10.6,3E12.4,I4,2E12.4)') & 729 & bcp(ii)%rr(:),bcp(ii)%ev(:),bcp(ii)%iat,bcp(ii)%ev(1)+bcp(ii)%ev(2)+bcp(ii)%ev(3),bcp(ii)%chg 730 end do 731 732 write(untc,'(I4, " :RCP''s, coordinates, laplacian eigenvalues, sum of these, density")') nrcp 733 do ii=1,nrcp 734 vv(:)=rcp(ii)%rr(:)+xorig(:) 735 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 736 write(untc,'(3F10.6,3E12.4,2E12.4)') & 737 & rcp(ii)%rr(:),rcp(ii)%ev(:),rcp(ii)%ev(1)+rcp(ii)%ev(2)+rcp(ii)%ev(3),rcp(ii)%chg 738 end do 739 740 write(untc,'(I4, " :CCP''s coordinates, laplacian eigenvalues, sum of these, density")') nccp 741 do ii=1,nccp 742 vv(:)=ccp(ii)%rr(:)+xorig(:) 743 call vgh_rho(vv,chg,grho,hrho,dist,iat,ipos,0) 744 write(untc,'(3F10.6,3E12.4,2E12.4)') & 745 & ccp(ii)%rr(:),ccp(ii)%ev(:),ccp(ii)%ev(1)+ccp(ii)%ev(2)+ccp(ii)%ev(3),ccp(ii)%chg 746 end do 747 748 end if ! End the condition on aim_dtset%crit > 0 749 750 !Reading of the CPs from the file 751 752 if (aim_dtset%crit==-1) then 753 read(untc,*) nbcp 754 ABI_DATATYPE_ALLOCATE(bcp,(nbcp)) 755 do ii=1,nbcp 756 read(untc,*) bcp(ii)%rr(:) 757 end do 758 read(untc,*) nrcp 759 ABI_DATATYPE_ALLOCATE(rcp,(nrcp)) 760 do ii=1,nrcp 761 read(untc,*) rcp(ii)%rr(:) 762 end do 763 read(untc,*) nccp 764 ABI_DATATYPE_ALLOCATE(ccp,(nccp)) 765 do ii=1,nccp 766 read(untc,*) ccp(ii)%rr(:) 767 end do 768 end if 769 770 do ii=1,nbcp 771 pc(:,ii)=bcp(ii)%rr(:) 772 icpc(ii)=-1 773 end do 774 do ii=1,nrcp 775 pc(:,nbcp+ii)=rcp(ii)%rr(:) 776 icpc(nbcp+ii)=1 777 end do 778 do ii=1,nccp 779 pc(:,nbcp+nrcp+ii)=ccp(ii)%rr(:) 780 icpc(nbcp+nrcp+ii)=3 781 end do 782 npc=nbcp+nrcp+nccp 783 784 !Checking 785 786 if (allocated(bcp)) then 787 do ii=1,nbcp 788 do jj=1,npc 789 iat=bcp(ii)%iat 790 ipos=bcp(ii)%ipos 791 if ((iat/=0).and.(ipos/=0)) then 792 pom(:)=pc(:,jj)+xorig(:)-xatm(:,iat)-atp(:,ipos) 793 ss=aim_dtset%coff2*vnorm(pom,0) 794 if (rminl(iat) >= ss) rminl(iat)=ss 795 end if 796 end do 797 end do 798 ABI_DATATYPE_DEALLOCATE(bcp) 799 end if 800 do ii=1,natom 801 write(std_out,*) 'atom: ', ii, rminl(ii) 802 end do 803 804 if(allocated(rcp)) then 805 ABI_DATATYPE_DEALLOCATE(rcp) 806 end if 807 if(allocated(ccp)) then 808 ABI_DATATYPE_DEALLOCATE(ccp) 809 end if 810 811 !END CP ANALYSIS 812 813 call timein(ttcp,wall) 814 ttcp=ttcp-tt0 815 816 end subroutine cpdrv